© ESO, 2017

1. Introduction

The BL Lac type object Markarian (Mrk) 501 is among the most prominent members of the class of blazars. Owing to its brightness, almost the entire broadband spectral energy distribution (SED) of Mrk 501 can be measured accurately with the current instrumentation. It is also known as one of the most active blazars, showing very strong and fast variability on timescales as short as a few minutes (Albert et al. 2007a). Moreover, because of its low redshift of z = 0.034, even the multi-TeV γ rays are influenced only weakly by the absorption on the extragalactic background light (EBL). Altogether, this makes Mrk 501 an excellent candidate source to study flux and spectral variability in the broadband emission of blazars.

Mrk 501 was the second extragalactic object to be detected in very high energy (>100 GeV, hereafter VHE) γ rays (Quinn et al. 1996; Bradbury et al. 1997), and it has been the subject of extensive studies in the different accessible energy bands over the last two decades. Based on its SED, it has been classified as a high frequency peaked BL Lac-type source (HBL) according to Padovani & Giommi (1995), or high-synchrotron peaked BL Lac (HSP) if following the classification given in Abdo et al. (2010b).

In 1997 Mrk 501 was found to be in an exceptionally high state, with the emission at VHE energies being up to 10 times the flux of the Crab Nebula (Protheroe et al. 1997; Djannati-Ataiet al. 1999). During this large flare, the synchrotron bump appeared to peak at or above 100 keV, indicating a shift of the peak position f quiescent state by at least two orders of magnitude (Catanese et al. 1997; Pian et al. 1998; Villata & Raiteri 1999; Tavecchio et al. 2001a). Over the following years, the source was intensively monitored at X-rays and VHE γ rays (e.g. Kataoka et al. 1999; Quinn et al. 1999; Sambruna et al. 2000; Aharonian et al. 2001; Massaro et al. 2004), and additional studies were done with the collected data a posteriori (e.g. Gliozzi et al. 2006). The observations could be well reproduced in the scope of one-zone synchrotron self-Compton (SSC) models. In 2005, the source showed another strong flaring event, for which flux-doubling times down to two minutes were measured at VHE (Albert et al. 2007a). This fast variability is a strong argument for a comparatively small emission region (with R ≈ 10¹⁵ cm), while the typical activity of the source could still be accommodated in models assuming a radius of the emission region which is larger by one to two orders of magnitude (e.g. Abdo et al. 2011a). Throughout the observations, the SED at the highest energies appeared to be harder in higher flux states (e.g. Albert et al. 2007a). Together with the observed shift of the synchrotron peak during the 1997 event, this suggests a change in the electron energy distribution as the cause for flaring events (Pian et al. 1998), but long-term changes in the Doppler factor or the size of the emission region are also being discussed as a possibility (Mankuzhiyil et al. 2012).

High-resolution radio images revealed a comparatively slow moving jet that features a limb brightening structure (Piner et al. 2008, 2009; Giroletti et al. 2008). The radio core position of Mrk 501 has been found to be stationary within 2 parsec (pc), using observations from the observing campaign in 2011 with the VLBI Exploration of Radio Astrometry (VERA, Koyama et al. 2015a), although variations in its location on year timescales cannot be excluded. High-resolution Global mm-VLBI Array (GMVA) observations at 86 GHz during the observing campaign in 2012 detected a new feature in the jet of Mrk 501, located 0.75 milliarcseconds (mas) southeast of the radio core (which corresponds to ~0.5 pc de-projected distance), and one order of magnitude dimmer than the core (Koyama et al. 2015b). This radio feature is consistent with the one reported in Giroletti et al. (2008) using GMVA data from 2005. This confirms that there are several distinct regions in the jet of Mrk 501, possibly stationary on year timescales, with the presence of high-energy electrons which could potentially produce optical, X-ray, and γ-ray emission, in addition to the emission detected with these high-resolution radio instruments.

Even though Mrk 501 has been studied over a comparatively long time, clear constraints on the properties of the highest activity regions, and on the particle populations involved, are still to be set. In this paper we present an extensive multi-instrument campaign on Mrk 501 that was conducted in 2009 in order to shed light on some of these open questions. This paper is a sequel to Abdo et al. (2011a), where, among other things, the averaged broadband SED from the campaign was studied in detail. A study focused on the flaring activity of May 1 (MJD 54 952), which includes very fast variability detected with the Whipple 10 m telescope, VERITAS light curves and spectra, and some measurements of the optical polarization performed by the Steward Observatory are reported in Pichel & Paneque (2011) and Aliu et al. (2016). In the work presented here, we address the variability seen during the full campaign, possible interband correlation of flux changes, and the characterization of the measured SED during two states of increased activity. While Aliu et al. (2016) looks at the average X-ray spectrum for a low state covering three weeks and a high state covering three days of the first VHE enhancement, we perform a detailed investigation characterizing the X-ray spectra for each pointing available for the campaign, hence providing a better quantification of the X-ray spectral variability. Furthermore, we consider an expanded data set, which also includes radio observations performed with the Very Long Baseline Array (VLBA), measuring the radio flux coming from the entire source and the radio flux from the compact core region only, and additional measurements of the optical polarization performed by the Steward and St. Petersburg observatories before and after the flaring activity of May 1.

This paper is structured as follows. In Sect. 2 an overview of the multi-instrument campaign is given, and updates with respect to the information provided in Abdo et al. (2011a) are discussed. In Sect. 3 the collected light curves and spectra are assessed for variability and interband correlation. The discussion of the broadband spectral energy distributions and a quantification of these measured spectra within synchrotron self-Compton scenarios by means of a novel technique based on a scan over the full parameter range is reported in Sect. 4. Finally, the results are discussed in Sect. 5, and a short summary and concluding remarks are given in Sect. 6.

2. Multi-instrument observing campaign performed in 2009

The presented multiwavelength (MWL) campaign was conducted over 4.5 months in 2009. The aim of this campaign was to sample the SED over all wavelengths every ~5 days. This way, the intrinsic flux variability of the source could be probed during non-flaring activity, hence reducing the observational bias towards states of high activity, which are the main focus of target of opportunity (ToO) campaigns. The covered frequency range spans from radio to VHE γ rays, including data from ~30 different instruments. The campaign took place from 2009 March 15 (MJD 54 905) to 2009 August 1 (MJD 55 044). Good coverage was achieved, while the sampling density varies among the different wavelengths because of different duty cycles and observational constraints of the participating instruments. The individual data sets and the data reduction are described in detail in Table 1 and Sect. 5 of Abdo et al. (2011a), and are not reported again in this paper. In this section we only briefly mention the various observations performed, and report on the updates of some data analyses and on extended data sets.

In the radio band, several single-dish instruments were used, namely the Effelsberg 100 m radio telescope, the 32 m Medicina radio telescope, the 14 m Metsähovi radio telescope, the 32 m Noto radio telescope, the Owens Valley Radio Observatory (OVRO) 40 m telescope, the 26 m University of Michigan Radio Astronomy Observatory (UMRAO), and the 600 m ring radio telescope RATAN-600. The mm-interferometer Submillimeter Array (SMA) and the Very Long Baseline Array (VLBA) were also used during the campaign. These single-dishes and the SMA monitored the total flux of Mrk 501 as a point-like unresolved source at frequencies between 2.6 GHz and 225 GHz. The VLBA took data ranging from 5 GHz to 43 GHz through various programs (BP143, BK150, and MOJAVE). Owing to the better angular resolution of MOJAVE, in addition to the total flux of the source, measurements of the flux from the compact (~10^-3 pc) core region of the jet could be obtained through two-dimensional Gaussian fits to the observed data.

Observations in optical frequency ranges have been performed by numerous instruments distributed all over the globe. In the R band, the Abastumani, Lulin, Roque de los Muchachos (Kungliga Vetenskaplika Academy, KVA), St. Petersburg, Talmassons, and Valle d’Aosta observatories performed observations as part of GASP-WEBT, the GLAST-AGILE Support Program of the Whole Earth Blazar Telescope (e.g. Villata et al. 2008, 2009). Additional data with several optical filters were provided by the Goddard Robotic Telescope (GRT), the Remote Observatory for Variable Object Research (ROVOR), and the Multicolor Imaging Telescopes for Survey and Monstrous Explosions (MITSuME). At near-infrared wavelengths, measurements performed by the Guillermo Haro Observatory (OAGH) have been included in the data set. Also within the GASP-WEBT program, the Campo Imperatore took measurements in near-infrared frequencies (JHK bands). The data obtained in the optical and near-infrared regime used the calibration stars reported in Villata et al. (1998), and have been corrected for Galactic extinction following Schlegel et al. (1998).

Through various observing proposals related to this extensive MWL campaign, 29 pointing observations were performed with the Rossi-X-ray Timing Explorer (RXTE), and 44 pointing observations performed with the Swift satellite¹. These observations provided coverage in the ultraviolet frequencies with the Swift Ultraviolet/Optical Telescope (UVOT), and in the X-ray regime with the RXTE Proportional Counter Array (PCA) and the Swift X-ray Telescope (XRT). Swift/XRT performed 41 snapshot observations in Windowed Timing (WT) mode throughout the whole campaign, and three observations in Photon Counting (PC) mode around MJD 54 952. The PC observations had not been used in Abdo et al. (2011a). For PC mode data, events for the spectral analysis were selected within a circle of 20 pixel (~46 arcsec) radius, which encloses about 80% of the point spread function (PSF), centered on the source position. The source count rate was above ~5 counts s^-1 and data were significantly affected by pile-up in the inner part of the PSF. After comparing the observed PSF profile with the analytical model derived by Moretti et al. (2005), pile-up effects were removed by excluding events within a 4 pixel radius circle centered on the source position, and an outer radius of 30 pixels was used. Occasionally, during the first ~100 s of a WT mode observation, Swift/XRT data can display a deviation in the light curve that is not due to the source variability, but is instead due to the settling of the spacecraft pointing causing a hot column to come in and out of either the source or background region. We inspected these data for any such deviations that could significantly impact our analysis, and none were found.

While Mrk 501 can be significantly detected with XRT and PCA for each single observation (~0.3 h), integration times of ~30 days are required in order to obtain significant detections with the RXTE All-Sky Monitor (ASM) and the Swift Burst Alert Telescope (BAT). The advantage of “all-sky instruments” like RXTE/ASM and Swift/BAT is that they can observe Mrk 501 without specifically pointing to the source, and hence provide a more uniform and continuous coverage than pointed instruments like Swift/XRT and RXTE/PCA. Details on the analysis of the RXTE/ASM and Swift/BAT data were given in Abdo et al. (2011a).

The range of high-energy γ rays was covered with the Fermi Large Area Telescope (Fermi-LAT). As is the case with RXTE/ASM and Swift/BAT, the sensitivity of Fermi-LAT is quite moderate, but it provides a more uniform temporal coverage than the pointing instruments; to detect Mrk 501 typically it is necessary to integrate over ~15–30 days in order to have significant detections. In addition to the observations from the coordinated MWL campaign, here we also report on the X-ray/γ-ray activity of Mrk 501 measured with RXTE/ASM, Swift/BAT, and Fermi-LAT for a time interval spanning from MJD 54 800 to MJD 55 100, which exceeds the time span of the campaign.

The Fermi-LAT data were reanalyzed using the Pass 8 SOURCE class events, and the ScienceTools software² package version v10r1p1. We used all events (from MJD 54 800 to MJD 55 100) with energies from 200 MeV to 300 GeV and within a 10° region of interest (RoI) centered at the position of Mrk 501. In order to avoid contamination from the Earth limb γ rays, only events with zenith angles below 100° were used. We used the P8R2_SOURCE_V6 instrument response functions, and the gll_iem_v06 and iso_P8R2_SOURCE_V6_v06 models to parameterize the Galactic and extragalactic diffuse emission (Acero et al. 2016)³. Given that Mrk 501 is a relatively hard source, we only used events above 300 MeV for the spectral analysis, as was done in Abdo et al. (2011a). All point sources in the third Fermi-LAT source catalog (3FGL, Acero et al. 2015) located in the 10° RoI and an additional surrounding 5° wide annulus (called “source region”) were modeled in the fits, with the spectral parameters set to the values from the 3FGL, and the normalization parameters kept free only for the nine sources identified as variable (in the 3FGL) and located within 10° of Mrk 501. The normalization parameters for the two diffuse components were also kept free. The spectral analysis performed on 15- and 30-day time intervals from MJD 54 800 to MJD 55 100 led to spectra successfully described by a power-law (PL) function with an index compatible⁴ with Γ = 1.75. For the determination of the light curves in the two energy bands 0.2–2 GeV and >2 GeV that are reported in Sect. 3.1, we decided to fix the value of the PL index to Γ = 1.75.

MAGIC observations were carried out with a single telescope, as the second telescope was under construction during the campaign period. Owing to a scheduled upgrade, no data were taken with MAGIC between MJD 54 948 and MJD 54 960. All observations were carried out in “wobble” mode (Fomin et al. 1994). For the work presented here, the data underwent a revised quality check and were reanalyzed with an improved analysis pipeline with respect to the one presented in Abdo et al. (2011a). Compared to the analysis presented in the first publication, the data set has been expanded by several nights (MJD 54 937, 54 941, 54 944, 54 945, 54 973, 54 975, 55 035, 55 038). Three nights were rejected because of revised quality criteria (MJD 54 919, 54 977, 55 026). After all data selection and analysis cuts, the effective observation time covered by the data comprises 17.4 h, while the first analysis yielded 16.2 h of selected data.

VERITAS observed Mrk 501 with different telescope configurations over the duration of this campaign. The data presented here amounts to 9.7 h of effective time, and are identical to those presented in Abdo et al. (2011a). However, the work in this paper presents the VERITAS light curve for the first time.

The Whipple 10 m telescope observed Mrk 501 for 120 h throughout the campaign, separately from the VERITAS array. The data taken with the Whipple 10 m were not used in the first publication which focused on the average state of the source throughout the campaign (Abdo et al. 2011a). However, the Whipple 10 m data over a flaring period around May 1 have been recently reported in a separate paper (Aliu et al. 2016). For better comparison to the other VHE instruments, Whipple 10 m fluxes, originally computed as flux in Crab Units (C.U.) above 400 GeV, were converted into fluxes above 300 GeV using the Crab flux above 300 GeV of F_{> 300 GeV} = 1.2 × 10^-10 cm^-2 s^-1 (Aleksić et al. 2012).

For more details on the observation strategy, list of instruments, and analysis procedures performed for the different instruments, the reader is referred to Abdo et al. (2011a) and references therein.

In addition to the MWL observations conducted as part of the campaign, the data set was expanded with measurements of the optical polarization performed by the Steward (Bok telescope) and St. Petersburg (LX-200) observatories from February to September 2009. The LX-200 polarization measurements were obtained from R-band imaging polarimetry, while the measurements from the Steward Observatory were derived from spectropolarimetry between 4000 and 75 00 Å with a resolution of ~15 Å, and the reported values are constructed from the median Q/I and U/I in the 5000−7000 Å band. The effective wavelength of this bandpass is similar to the Kron-Cousins R band, and the wavelength dependence in the polarization of Mrk 501 seen in the spectro-polarimetry is small and does not significantly affect the results. The details related to the observations and analysis of the polarization data is reported in Larionov et al. (2008) and Smith et al. (2009). The Steward observations are part of the public Steward Observatory program to monitor γ-ray bright blazars during the Fermi-LAT mission⁵, and a fraction of these polarization observations have been recently reported in Aliu et al. (2016).

3. Multi-instrument flux and spectral variability

During the 4.5-month MWL campaign, Mrk 501 was observed with numerous instruments covering the entire broadband SED. In the following section, we report the measured multiband flux and spectral variability, as well as multiband correlations.

3.1. Multi-instrument light curves

The light curves which were derived from pointed observations in the different energy bands, spanning from radio to VHE γ rays, are shown in Fig. 1. Figure 2 presents the X-ray and γ-ray activity as measured with the all-sky instruments RXTE/ASM, Swift/BAT and Fermi-LAT.

The light curves obtained during pointing observations in the radio regime exhibit a nearly constant flux at a level of ~1.2 Jy. The well-sampled light curve taken with the OVRO telescope shows constant emission of 1.158 ± 0.003 Jy.

The measurements performed with the VLBA at a frequency of 43 GHz are presented in Fig. 3. A constant fit delivers a reduced χ² of 8.4/3 for the total flux and 15.6/3 for the core flux, yielding a probability for the data points to be well described by a constant fit of 3.8% and 0.14%, respectively. Although marginally significant, this suggests an increase in the radio flux in May 2009 (dominated by the core emission) in comparison to that measured during the other months.

For the near-infrared observations in Fig. 1, flux levels of ~40−50 mJy (J and K bands) and ~50−60 mJy (H band) were measured. Only small variations can be seen, even though the sampling is less dense and the uncertainties of the measurements are comparatively large. For the extensive data sample in the optical regime, a nearly constant flux was measured at flux levels of ~ 6 mJy (B band), 11 mJy (V), ~16 mJy (R), and 24−29 mJy (I/Ic). No correction for emission by the host galaxy was applied. At ultraviolet frequencies, a flux level of ~2 mJy with flux variations of about 25% over timescales of about 25 to 40 days can be seen.

The average Swift/XRT measured fluxes during the entire campaign are F_{0.3−2 keV} = (9.2 ± 0.3) × 10^-11 erg cm^-2 s^-1 in the energy range between 0.3 and 2 keV and F_{2−10 keV} = (7.2 ± 0.3) × 10^-11 erg cm^-2 s^-1 in the range 2−10 keV, while RXTE/PCA, due to a slightly different temporal coverage, measured an average 2−10 keV flux of F_{2−10 keV} = (7.8 ± 0.2) × e10^-11 erg cm^-2 s^-1.

The Fermi-LAT measured a variable flux in the two probed γ-ray bands, with an average flux of F_{0.2−2 GeV} = (2.75 ± 0.14) × 10^-8 ph cm^-2 s^-1 between 200 MeV and 2 GeV and F_{> 2 GeV} = (5.3 ± 0.4) × 10^-9 ph cm^-2 s^-1 at energies above 2 GeV (shown in Fig. 2). The highest emission is seen in the 15-day time interval between MJD 54 967 and MJD 54 982.

The VHE γ-ray light curves are shown in the upper panel of Fig. 1. The average flux above 300 GeV of Mrk 501 during the campaign, including the flaring activities, is about 5 × 10^-11 ph cm^-2 s^-1 (~0.4 C.U.)⁶. Flux variability is evident throughout the VHE light curve, in addition to flaring episodes of a few days occurring in MJD 54 952 (2009 May 1) and MJD 54 973 (2009 May 22).

In the following paragraphs we review the first VHE flare in a MWL context, and include additional details specifically on the X-ray data. We then provide details on the second VHE flare.

3.1.1. VHE γ-ray flaring event starting at MJD 54 952

On 2009 May 1, the Whipple 10 m telescope observed Mrk 501 for 2.3 h and detected a VHE flux (>300 GeV) increase from ~1.0−1.5 C.U. to ~4.5 C.U. in the first 0.5 h (from MJD 54 952.35 to MJD 54 952.37), which implies a flux increase of about one order of magnitude with respect to the average VHE flux level recorded during the full campaign. Following the alert by the Whipple 10 m, VERITAS started to observe Mrk 501 after 1.4 h (at MJD 54 952.41) and detected the source at a VHE flux of 1.5 C.U. without statistically significant flux variations during the full observation (from MJD 54 952.41 to MJD 54 952.48). This VHE flux level was also measured by the Whipple 10 m telescope in approximately the same time window (from MJD 54 952.41 to MJD 54 952.47), and corresponds to a VHE flux ~4 times larger than the typical flux level of 0.4 C.U. measured during the full campaign. The peak of the flare (which occurred at MJD 54 952.37) was caught only by the Whipple 10 m. Still, the Mrk 501 VHE γ-ray flux remained high for the rest of the night and the following two nights (until MJD 54 955), which was measured by VERITAS and the Whipple 10 m with very good agreement. Further details about the VERITAS and Whipple 10 m intra-night variability measured on 2009 May 1, and about the enhanced activity during the first days of May, can be found in Pichel & Paneque (2011) and Aliu et al. (2016).

During the period of the considered VHE γ-ray flare, no substantial increase in the X-ray regime can be claimed based on the Swift/XRT observations: the 0.3−2 keV and the 2−10 keV fluxes during this flaring episode are about ~ 8 × 10^-11 erg cm^-2 s^-1 and ~ 1 × 10^-10 erg cm^-2s^-1, which are about ~10% lower and ~30% higher than the average X-ray flux values reported above. However, the Swift/XRT observations started seven hours after the Whipple 10 m and VERITAS observations of this very high VHE state on MJD 54 952. The reason is that the XRT observations were taken within a ToO activated by the enhanced VHE activity measured by the Whipple 10 m and VERITAS unlike most of the planned X-ray observations from the MWL campaign which were coordinated with the VHE observations.

In the two lowest panels in Fig. 1, the evolution of the optical polarization degree and orientation are shown. The degree of polarization during the few days around the first VHE flaring activity is measured at 5% compared to a 1% measurement during several observations before and after this flaring activity. There is also a rotation of the EVPA by 15 degrees, which comes to a halt at the time of the VHE outburst when the degree of polarization drops from 5.4% to 4.5% (see further details in Pichel & Paneque 2011; Aliu et al. 2016).

3.1.2. VHE γ-ray flaring event starting at MJD 54 973

The MAGIC telescope observed Mrk 501 for 1.7 h on 2009 May 22 (MJD 54 973) and measured a flux of 1.2 C.U., which corresponds to ~3 times the low flux level. At the next observation on May 24 (MJD 54 975.00 to MJD 54 975.12), the flux had already decreased to a level of ~0.5 C.U. The Whipple 10 m observed Mrk 501 later on the same date (from MJD 54 975.25) and measured a flux of ~0.7 C.U., while the following day (from MJD 54 976.23) it measured a flux increase to 1.1 C.U. No VERITAS observations of Mrk 501 took place at this time, due to scheduled telescope maintenance.

The MAGIC data of the flaring night were probed for variations on timescales down to minutes, but no significant intra-night variability was found. Moreover, tests for spectral variability during the night in terms of hardness ratios vs. time in different energy bands showed no significant variations either.

Unfortunately, there are no X-ray observations that are strictly simultaneous with the MAGIC observations on MJD 54 973. The closest RXTE/PCA observations took place on MJD 56 969 and MJD 54 974, and the closest Swift/XRT observations are from MJD 54 970 and MJD 54 976, all of which show a flux increase (up to a factor of ~2) with respect to the average X-ray flux measured during the campaign.

Under the assumption that no unobserved intra-day variability occurred in the X-ray band, it can be inferred that Mrk 501 was in a state of increased X-ray and VHE activity over a period of up to 5 days. During this period there were no flux changes observed at optical or radio frequencies.

3.2. Spectral variability in individual energy bands

In this section we report on the spectral variability observed during the two VHE flaring episodes around the peaks of the two SED bumps, namely at X-ray and γ rays, where most of the energy is being emitted and where the flux variability is highest.

3.2.1. VHE γ rays

The VHE spectra measured with MAGIC and VERITAS, averaged over the entire campaign between 2009 March 15 (MJD 54 905) and 2009 August 1 (MJD 55 044), were reported in Abdo et al. (2011a). Only the time span MJD 54 952−54 955, where VERITAS recorded VHE flaring activity, was excluded for the average spectrum and was presented as a separate high-state spectrum (see Fig. 8 of Abdo et al. 2011a). The resulting average spectra relate to a VHE flux of about 0.3 C.U., which is the typical non-flaring VHE flux level of Mrk 501. Additionally, two spectra were obtained with the Whipple 10 m for that night: a very-high-state spectrum spanning MJD 54 952.35−54 952.41, which seems to cover the peak of the flare, and a high-state spectrum derived from the time interval MJD 54 952.41−54 955.00, which is simultaneous with the observations performed with VERITAS. These spectra were reported in Pichel & Paneque (2011) following the general Whipple analysis technique described in Horan et al. (2007), and further details from these spectra are reported in Aliu et al. (2016).

The reanalysis of the MAGIC data (see Sect. 2), which contains some additional data compared to the analysis presented in Abdo et al. (2011a), revealed a flaring state on MJD 54 973, for which a dedicated spectrum was computed. An averaged spectrum was derived based on the remaining data set. The energy distribution of the differential photon flux can be well described by a PL function of the form (1)yielding F₀ = ⁽9.3 ± 0.4⁾ × 10^-8 ph m^-2s^-1TeV^-1 and Γ = 2.40 ± 0.05. This new MAGIC averaged spectrum was found to be in agreement with the previously presented value where a power-law fit gave F₀ = (9.0 ± 0.5) × 10^-8 phm^-2s^-1TeV^-1 and Γ = 2.51 ± 0.05 (Abdo et al. 2011a). Here we only quote statistical uncertainties of the measurements. The systematic errors affecting data taken by the MAGIC telescope at the time of the presented campaign are discussed in Albert et al. (2008) and are valid for both analyses. They are estimated as an energy scale error of 16%, a systematic error on the flux normalization of 11%, and an error on the obtained spectral slope of ± 0.2. In the following, the more recent analysis result are used.

All the VHE γ-ray spectra described above are presented in Fig. 4. The spectra shown in the figure were corrected for absorption by the EBL using the model from Franceschini et al. (2008). Given the proximity of Mrk 501, the impact of the EBL on the spectrum is relatively weak: the attenuation of the flux reaches 50% at an energy of 5 TeV. Many other EBL models (e.g. Finke et al. 2010; Domínguez et al. 2011) provide compatible results at energies below 5 TeV, hence the results do not depend significantly on the EBL model used. The power-law fit parameters (see Eq. (1)) of the measured spectra (i.e. the spectra not corrected for EBL) can be found in Table 1. For spectra measured with MAGIC, the presented fits also take into account the correlation between the individual spectral points which is introduced by the unfolding of the spectrum, while no explicit unfolding has been applied for the other instruments. The average-state spectra measured by the three instruments (after subtracting the time intervals with strong flaring activity in the VHE) agree very well, despite the somewhat different observing periods. This suggests that these VHE spectra are a good representation of the typical VHE spectrum of Mrk 501 during this MWL campaign. The high-state spectra show a spectral slope that is harder than that from the non-flaring state, hence indicating a “harder when brighter” behavior, as has been reported previously (e.g. Albert et al. 2007a).

3.2.2. GeV γ rays

The two short VHE flaring episodes discussed in this paper occurred within the time interval MJD 54 952−54 982, which is the 30-day time interval with the highest flux and hardest GeV γ-ray spectrum reported in Abdo et al. (2011a). The flux above 300 MeV F_{> 300 MeV} and photon index Γ for this 30-day time interval computed using the ScienceTools software package version v9r15p6 and the P6_V3_DIFFUSE instrument response functions are F_{> 300 MeV} = (3.6 ± 0.5) × 10^-8 ph cm^-2 s^-1 and Γ = 1.64 ± 0.09, while values for the Fermi-LAT spectrum averaged for the entire MWL campaign are F_{> 300 MeV} = (2.8 ± 0.2) × 10^-8 ph cm^-2 s^-1 and Γ = 1.74 ± 0.05 (for further details, see Abdo et al. 2011a). Performing the analysis with the ScienceTools software package version v10r1p1 and the Pass 8 data (which implies somewhat different photon candidate events), as described in Sect. 2, led to a photon flux (above 300 MeV) of F_{> 300 MeV} = (4.2 ± 0.5) × 10^-8 ph cm^-2 s^-1 and a PL index of Γ = 1.68 ± 0.07 for the time interval MJD 54 952−54 982, and a flux (above 300 MeV) of F_{> 300 MeV} = (3.0 ± 0.2) × 10^-8 ph cm^-2 s^-1 and a PL index of Γ = 1.75 ± 0.04 for the entire campaign. The spectral results derived with Pass 6 and Pass 8 are compatible, and show a marginal increase in the flux and the hardness of the spectra during the time interval MJD 54 952−54 982 with respect to the full campaign period.

The Pass 8 Fermi-LAT data analysis is more sensitive than the Pass 6 data analysis, and allows us to detect Mrk 501 significantly (TS > 25)⁷ and to determine the spectra around these two flares in time intervals as short as 2 days centered at MJD 54 952 and 54 973. Additionally, for comparison purposes, we also computed the spectra for 7-day time intervals centered at MJD 54 952 and 54 973⁸. The Fermi-LAT spectral results for the various time intervals in May 2009 are reported in Table 2. For the first flare, for both the 2-day and 7-day time intervals, the LAT analysis yields a signal with TS ~ 40. This shows that increasing the time interval from 2 days to 7 days did not increase the γ-ray signal, and hence indicates that the 2-day time interval centered at MJD 54 952 dominates the γ-ray signal from the 7-day time interval. The spectrum is marginally harder than the average spectrum from the time interval MJD 54 952−54 982. For the second flare, the 7-day time interval yields a signal significance () 2.6 times larger than that of the 2-day time interval, showing that, contrary to the first flare, increasing the time interval from 2 days to 7 days enhanced the γ-ray signal considerably. The Fermi-LAT spectrum around the second flare is very similar to the average spectrum obtained for the 30-day time interval MJD 54 952−54 982.

For the MWL SEDs presented in Fig. 9, we show the Fermi-LAT spectral results for these two flares performed on three and five differential energy bins (starting from 300 MeV). Here, the shape of the spectrum was fixed to that obtained for the full range for each temporal bin. Upper limits at 95% confidence level were computed whenever the TS value (for the γ-ray signal of the bin) was below six and/or the uncertainty was equal to or larger than the energy flux value.

3.2.3. X-rays

In the X-ray band, individual spectra could be derived for each pointing of the two instruments Swift/XRT and RXTE/PCA. Both indicated significant variability in flux and spectral index during the course of the campaign. Figure 5 shows the XRT and PCA spectra around the times of the first and second flux increase in the VHE range. For the first flare, the variability in flux and spectral shape is greater for XRT than for PCA, but mostly because many of the XRT observations were performed within a ToO program, and so they provide a better characterization of the enhanced activity (see Sect. 2).

Around the first VHE flare, the XRT spectra tend to be much harder and appear to show an upward curvature towards higher energies. The hardening of the spectrum is confirmed by a spectral analysis performed using a power-law spectral model with the hydrogen density N_H fixed to the Galactic value. Figure 6 shows the spectral index light curve (see also Table 3) and the reduced χ² of the individual fits. Based on the reduced χ² values, the representation by a simple power-law function is sufficient for most spectra. Around MJD 54 952−54 953, which roughly corresponds to the time of the first VHE flare, a peak in the hardness of the spectrum can be seen.

Around the second flux increase in the VHE γ-ray band, variability was seen by both Swift/XRT and RXTE/PCA, with flux changes of up to a factor of 2 with respect to the flux average of ~(7−8) × 10^-11 erg cm^-2 s^-1 in the 2−10 keV band (see Fig. 1). However, no particular hardening of the spectrum was found (see Fig. 6), as was observed for the first flare.

3.3. Quantification of the multi-instrument variability

As a quantitative study of the underlying variability seen at different wavelengths, the fractional variability F_var was determined for each instrument according to Eq. (10) in Vaughan et al. (2003), (2)where S² represents the variance, specifies the mean square error stemming from measurement uncertainties, and ⟨ F_γ ⟩ is the arithmetic mean of the measured flux. The term under the square root is also known as the normalized excess variance .

The uncertainty of F_var is calculated following the prescription in Poutanen et al. (2008), as described in Aleksić et al. (2015a), so that they are also valid in the case when ΔF_var ~ F_var, (3)with the error of the normalized excess variance as defined in Eq. (11) in Vaughan et al. (2003).

This method for quantifying the variability comes with the caveat that the resulting F_var and related uncertainty depend very much on instrument sensitivity and the observing sampling, which is different for the different energy bands. In other words, a densely sampled light curve with small uncertainties in the flux measurements may allow us to see flux variations that are hidden otherwise, and hence may yield a larger F_var and/or smaller uncertainties in the calculated values of F_var. Some practical issues in the application of this methodology in the context of multiwavelength campaigns are elaborated in Aleksić et al. (2014, 2015b,a).

For Swift/XRT and RXTE/PCA in the X-ray band, and MAGIC, VERITAS, and the Whipple 10 m in the VHE regime, the fractional variability has been calculated for the full data set and also after removal of the temporal intervals related to the two flaring episodes (MJD 54 952−54 955, MJD 54 973−54 978). The fractional variability specifically computed for the period around the first flaring episode has been recently reported in Aliu et al. (2016). For measurements in the optical R band, F_var has additionally been calculated for optical fluxes corrected for the host galaxy emission as derived in Nilsson et al. (2007). For data sets containing fewer than five data points, no F_var was calculated. The results are presented in Fig. 7.

A negative excess variance was obtained for data sets from the following instruments: UMRAO (at 5 GHz and 8 GHz), Noto (at 8 GHz and 43 GHz), Medicina (at 8 GHz), Effelsberg (all bands), and the near-IR measurements within the GASP-WEBT program (all bands). Such a negative excess variance is interpreted as an absence of flux variability within the sensitivity range of the instrument. These data sets have not been included in Fig. 7.

At low frequencies, from radio to optical, no substantial variability was detected, with F_var ranging from ≈ 0.02−0.06 in radio to 0.01−0.1 in optical. In the X-ray band, we find F_var ≈ 0.3, indicating substantial variation in the flux during the probed time interval. After removal of the flaring times, variabilities of F_var ≈ 0.2−0.25 are still seen. The fractional variability in the γ-ray band covered by Fermi-LAT is on the order of F_var ≈ 0.3−0.4; yet the Fermi-LAT F_var values are not directly comparable to the other instruments, as GeV variability on day timescales, which could be higher than that computed (separately) for the 15-day and the 30-day timescales, cannot be probed. Strong variability can be noted at VHE with F_var ≥ 0.4 for the data sets without the flares and F_var ≥ 0.6 (0.9 for Whipple 10 m) for observations including the flaring episodes.

All in all, Mrk 501 showed a large increase in variability with increasing energy, ranging from an almost steady behavior at the lowest frequencies to the highest variability observed in the VHE band.

3.4. Multi-instrument correlations

To study possible cross-correlations of flux changes between the different wavelengths, we determined the discrete correlation functions (DCF), following Edelson & Krolik (1988), based on the light curves obtained by the various instruments. The DCF allows a search for correlations with possible time lags, which could result for example from a spatial separation of different emission regions. We probed time lags in steps of 5 days up to a maximum shift of 65 days. The step size corresponds to the overall sampling of the light curve and thus to the objective of the MWL campaign itself, which was to probe the source activity and spectral distribution every ≈ 5 days. The maximum time span is governed by the duration of the campaign, as a good fraction of the light curve should be available for the calculation of cross-correlations. We chose a maximum of 65 days, which corresponds to roughly half the time span of the entire campaign. Because of the uneven sampling and varying exposure times, the significance of the correlations derived from the prescription given in Edelson & Krolik (1988) might be overestimated (Uttley et al. 2003). We derived an independent assessment of the significance of the correlation by means of dedicated Monte Carlo simulations as described in Arévalo et al. (2009) and Aleksić et al. (2015b,a).

In this study, possible cross-correlations between instruments of different wavelengths were examined. As already suggested by the low level of variability in the radio and optical band throughout the campaign, no correlations with any other wavelengths were found for these instruments. A correlation with flux changes in the MeV−GeV range could not be probed on timescales of days due to the integration time of 15−30 days required by Fermi-LAT for a significant detection. A similar situation occurs in the X-ray bands from Swift/BAT and RXTE/ASM, which also need integration times of the same order, and are thus also neglected for day-scale correlation studies.

Therefore, the study focuses on the highly sensitive X-ray and VHE γ-ray observations, namely the ones performed with Swift/XRT, RXTE/PCA, MAGIC, VERITAS, and Whipple 10 m, which are also the ones that report the highest variability (see Fig. 7). In the VHE γ-ray band, the number of observations is relatively small (in comparison to the number of X-ray observations performed with Swift and RXTE), and hence we compile a single light curve with a dense temporal sampling of Mrk 501, including the measured flux points from all three participating VHE γ-ray telescopes. This procedure is straightforward as VERITAS and MAGIC both measured the flux above 300 GeV and the Whipple 10 m measurements have been scaled to report a flux in the same energy range (see Sect. 2). We also combined measurements by Swift/XRT in the 2−10 keV band and data points from RXTE/PCA to a single light curve, as the same energy range is covered by the two instruments. The light curve in the 0.3−2 keV band only consists of measurements performed by Swift/XRT. The DCF vs. time-shift distributions for the two X-ray bands and the VHE γ-ray measurements are shown on the left-hand side of Fig. 8.

At a time lag on the order of −20 to −25 days, a hint of correlation at the level of 2σ between fluxes in the soft X-ray band and the VHE γ-ray band is seen in the top left panel of Fig. 8. This feature is dominated by the two flaring events , as the dominant flare in VHE γ rays occurred around MJD 54 952, while the largest flux increase in soft X-rays was seen around MJD 54 977 , with a separation of 24−25 days. The right-hand side of Fig. 8 reports the evaluation of the correlations after the flaring episodes have been excluded from the X-ray and VHE γ-ray light curves. The above-mentioned feature at 20−25 days is no longer present.

The large growth of the confidence intervals apparent at time shifts of ΔT ≈ 40 days are caused by sparsely populated regions in the VHE γ-ray light curve, mainly towards the end of the campaign. When the light curves are shifted by ≈ 40 days with respect to each other, these regions overlap with densely populated regions in the X-ray light curves, which results in a larger uncertainty of the determined DCF.

Overall, no significant correlation between X-ray and VHE γ-ray fluxes is found for any of the combinations probed.

4. Evolution of the spectral energy distribution

The time-averaged broadband SED measured during this MWL campaign (from MJD 54 905 to MJD 55 044) was reported and modeled satisfactorily in the context of a one-zone synchrotron self-Compton (SSC) scenario (Abdo et al. 2011a). In this model, several properties of the emission region are defined, such as the size of the region R, the local magnetic field B, and the Doppler factor δ, which describes the relativistic beaming of the emission towards the observer. Furthermore, the radiating electron population is described by a local particle density n_e and the spectral shape. For the averaged data set of this campaign, the underlying spectrum of the electron population was parameterized with a power-law distribution from a minimum energy γ_min to a maximum energy γ_max, with two spectral breaks γ_break,1 and γ_break,2. The two breaks in the electron energy distribution (EED) were required in order to properly model the entire broadband SED. Because of the relatively small multiband variability during the 4.5-month observing campaign (once the first VHE flare was removed) and the large number of observations performed with all the instruments, the average SED could be regarded as a high-quality representation of the typical broadband emission of Mrk 501 during the time interval covered by the campaign, and hence the one-zone SSC model was constrained to describe all the data points (including 230 GHz SMA and interferometric 43 GHz VLBA observations).

In this work, we focus on the characterization of the broadband SED during the two flaring episodes occurring in May 2009. As reported in Sect. 3.1, these two flaring episodes start on MJD 54 952 and MJD 54 973, and last for approximately three and five days, respectively. There is some flux and spectral variability throughout these two flaring episodes, but for the sake of simplicity in this section we attempt to model only the SEDs related to the VHE flares on MJD 54 952 and 54 973, which are the first days of these two flaring activities. We try to model these two SEDs with the simplest leptonic scenarios, namely a one-zone SSC and a two-independent-zone SSC model. In the latter we assume that the quiescent or slowly changing emission is dominated by one region that is described by the SSC model parameters used for the average/typical broadband emission from the campaign (see Abdo et al. 2011a), while the flaring emission (essentially only visible in the X-ray and γ-ray bands) is dominated by a second independent and spatially separated region.

The assumption of a theoretical scenario consisting of one (or two) steady-state homogenous emission zone(s) could be an oversimplification of the real situation. The blazar emission may be produced in inhomogeneous regions, involving stratification of the emitting plasma both along and across a relativistic outflow, and the broadband SED may be the superposition of the emission from all these different regions, characterized by different parameters and emission properties, as reported by various authors (e.g. Ghisellini et al. 2005; Graff et al. 2008; Giannios et al. 2009; Chen et al. 2011; Zhang et al. 2014; Chen et al. 2015). In this paper we decided to continue using the same theoretical scenario used in Abdo et al. (2011a), which we adopted as the reference paper for this data set. We also kept the discussion of the model parameters at a basic level, and did not attempt to perform a profound study of the implications of these parameters.

In this work we used the SSC model code described in Takami (2011), which is qualitatively the same as the one used in Abdo et al. (2011a), with the difference that the EED is parameterized as (4)where n_e is the electron number density. For reasons of comparability, only this definition is applied in all the SED modeling results in this section, including that of the quiescent, averaged SED obtained over the full MWL campaign. The corresponding one-zone SSC model parameter values defining the averaged SED from the full 2009 multi-instrument campaign are listed in Table 4. The parameter values are identical to those from the “Main SSC fit” reported in Table 2 of Abdo et al. (2011a), with the only difference being the usage of the electron number density n_e, instead of the equipartition parameter. The contribution of starlight from the host galaxy can be approximately described with the template from Silva et al. (1998), as was done in Abdo et al. (2011a).

For the characterization of the SEDs collected during the two flaring states, we allow for an EED with two spectral breaks in the case of one-zone SSC models. For the second zone in the two-zone SSC scenario, we keep the somewhat simpler description of the electron energy distribution as a power law with only one spectral break (i.e. α₂ ≡ α₁ in Eq. (4)).

4.1. Grid-scan strategy for modeling the SED

In contrast to the commonly used method of adjusting the model curve to the measured SED data points (e.g. Tavecchio et al. 1998, 2001b; Albert et al. 2007a, in this study we applied a novel variation on the grid-scan approach in the space of model parameters. Given a particular theoretical scenario (e.g. the one-zone or two-zone SSC model), we make a multi-dimensional grid with the N model parameters that we want to sample. For each parameter, we define a range of allowed values and a step size for the variation within this range. Theoretical (SSC) model curves are calculated for each point on the grid, i.e. for each combination of the N parameter values. Subsequently, the goodness of the resulting model curves in reconstructing the data points is quantified by means of the χ² between data and model, which takes into account the statistical uncertainties of the individual measurements. At the moment, systematic uncertainties are not considered for the evaluation of the agreement. This would require performing the entire procedure for various shifts in the flux and energy scale for each instrument, as well as for possible distortions in the individual spectra. The net impact of including systematic uncertainties in the single-instrument spectra would be a larger tolerance for the agreement between the experimental data and the theoretical model curves, which would yield a larger degeneracy in the parameter values that can model the data. While this will be investigated in the future, it is beyond the scope of this paper. Therefore, the data-model agreements reported in this manuscript, which are based on the χ² analysis using only the statistical uncertainties, provide a lower limit to the actual agreement between the presented experimental data and the theoretical model curves being tested, and we mostly use them to judge the relative agreement of the various theoretical model curves.

Depending on the complexity of the model itself, the model calculations for an entire grid can be very intensive in computing power. For instance, one of the simplest SSC scenarios, involving only one emission zone with an electron energy distribution with one spectral break, already leads to a grid spanning a nine-dimensional parameter space. With the ranges and grid spacings we are using in this work, the number of model curves to calculate and evaluate amounts to tens of million. For this reason, the access to cluster computing becomes essential for this grid-scan modeling approach. The model calculations in this work have been performed using the computing farms at SLAC⁹ and TU Dortmund¹⁰.

After the evaluation of all models regarding their level of agreement with the data, individual models can be chosen for the final set, according to the achieved probability of agreement (derived from the χ² and the number of degrees of freedom). These sets of models can then be visualized both in the SED representation and in the space of parameter values defining the models, which could populate non-continuous regions in the parameter space.

One aim of the grid-scan strategy is to keep the range of model parameters as wide as possible. By sampling a large parameter phase space we can reduce the bias which is usually introduced into the model by adopting a set of assumptions or educated choices. In addition to the obvious aim of finding parameter values which describe the data in the best way, another advantage is that the “grid-scan” approach also offers the possibility of investigating the degeneracy of the model-to-data agreement regarding each individual model parameter. In order to do this, sets of models within bands of achieved fit probabilities are compiled and their distributions in each of the model parameters are visualized. Based on such plots, interesting regions in the parameter space can be selected for a deeper search, which leads to models with an even better agreement with the data and to a more thorough study of the degeneracy of individual model parameters. Finally, the grid-scan method can find multiple clusters or regions in the model parameter phase space that could be related to different physical scenarios, which can be equally applicable to the data set at hand, but might be missed by statistical methods aiming at only “one best” solution.

The concept of grid-scan SED modeling has already been presented in Cerruti et al. (2013), where model curves are computed for each point on the parameter space grid, but the assessment of the agreement between model and data is performed in a different way: the authors evaluate the agreement based on seven observables (i.e. the frequency and luminosity of the synchrotron peak, the measured X-ray spectral slope and the GeV and TeV spectral slopes and flux normalizations), which are derived from the model curves and are compared to the data. They also provide a family of solutions involving any uncertainties in the observables. In the work presented here, the model-to-data agreement criterion, which is used to select a set of models, is derived directly from the χ²-distances between each data point and each model curve without computing any secondary characteristics of the SED which may introduce additional uncertainties. Cerruti et al. (2013) also determine this distance for the models picked by their algorithm, but apply it only as an a posteriori check of their result. Furthermore, the authors have reduced the dimensionality of the parameter space from nine to six, and used only five steps for each parameter, which implies the creation of a grid with 5⁶ = 15 625 SED realizations. In the work presented here, the smallest grid-scan implies the creation of more than 40 × 10⁶ SED realizations. Additionally, after selecting interesting regions in the various model parameters with the grid-scan, we went one step further and performed a second (dense) grid-scan focused only on those regions, and using a smaller step size.

The objective of finding uncertainty ranges of model parameters has also been addressed by Mankuzhiyil et al. (2011) and Zabalza (2015). Here, a Markov chain Monte Carlo procedure is used to fit emission model curves (for a number of different emission models) to the observational results. While this approach delivers uncertainties or probability distribution functions for the particle distribution parameters, this is done only for one particular solution. Disjointed regions of equally good model configurations, i.e. “holes” in the probability distribution for the individual parameters, are not found following this method.

A three-dimensional parameter grid with 9504 (48 × 22 × 9) steps was used by Petry et al. (2000) to find the most suitable model parameter set to describe weekly averaged SEDs of Mrk 501, where the “best” model was selected as the one with the smallest data-model difference, quantified with a χ² approach. Although a parameter grid was used, the goal and merits of that work differ from those of the methodology presented here. While Petry et al. (2000) used the three-dimensional parameter grid to find the best model (as in Mankuzhiyil et al. 2011, with a χ² minimization procedure), in this work the nine-dimensional grid is used to find the family (or families) of parameter values that give a good representation of the broadband SED, and to show the large degeneracy of the model parameters to describe the SED.

For the theoretical SED modeling of the two flaring states of Mrk 501, following the grid-scan strategy outlined above, the parameter ranges given in Tables 5 and 6 have been investigated for the one-zone and two-zone scenarios described at the beginning of this section. Given that we aim to sample a wide range of parameter values with a relatively coarse step (for each parameter), we denote these scans as “coarse grid-scans”. The general orientation for the choice of parameter ranges is based on previous works on modeling of the SED of Mrk 501, e.g. Albert et al. (2007a), Anderhub et al. (2009), Abdo et al. (2011a), Mankuzhiyil et al. (2012). Based on these values¹¹, one-zone SSC models have been built and second zones for the two-zone scenario as well. In the latter, the first zone is described by the model reproducing the average emission seen over the entire campaign (see Abdo et al. 2011a), while only the second zone is varied as described by the model parameters from the grid presented here. The phase space of the grid-scan could have been reduced by imposing a relation between the locations of the breaks (γ_break) and the size R and magnetic field B values, and by forcing the change of index before and after the breaks to be one (i.e. Δα = 1). However, cooling breaks with a spectral change two times larger than the canonical value of one were necessary to describe the broadband SED of Mrk 421 within a SSC homogeneous model scenario (see Sect. 7.1 of Abdo et al. 2011b), and the breaks needed by the SSC models are not always related to the cooling of the electrons, but instead could be related to the acceleration mechanism, as reported for Mrk 501 in Abdo et al. (2011a). Internal breaks (related to the electron acceleration) have been reported for various blazars (e.g. Abdo et al. 2009; Abdo et al. 2010a). The origin of these internal breaks, as well as large spectral changes at the EED breaks, may be related to variations in the global field orientation, turbulence levels sampled by particles of different energy, or gradients in the physical quantities describing the system. These characteristics are not taken into account in the relatively simple homogenous SSC models, and argue for more sophisticated theoretical scenarios like the ones mentioned above. In order to keep the range of allowed model parameter values as broad as possible, in this exercise we did not impose constraints on the location of the EED breaks or in the index values before or after the breaks. The hardest index we use in this study is 1.7, which is harder than the canonical index values > 2 derived from shock acceleration mechanisms and used very often to parameterize the broadband SEDs of blazars. But this is actually not a problem as various authors have shown that indices as hard as 1.5 can be produced through stochastic acceleration (e.g. Virtanen & Vainio 2005) or through diffusive acceleration in relativistic magnetohydrodynamic shocks, as reported in Stecker et al. (2007), Summerlin & Baring (2012) and Baring et al. (2016). We also use γ_min values extending up to 10⁶, substantially higher than those used in conventional SSC models (which typically go up to ~ 10³), but such high γ_min values have already been used by various authors (e.g. Katarzyński et al. 2006; Tavecchio et al. 2009; Lefa et al. 2011a,b).

In the evaluation of the models, we used two other constraints in addition to the requirement of presenting a good agreement with the SED data points. Equipartition arguments impose the condition that the energy densities held by the electron population (u_e) and the magnetic field (u_B) should be of comparable order. Typically, the parameterization of the broadband SED of Mrk501 (and all TeV blazars in general) within SSC theoretical scenarios require u_e ~ 10²⁻³u_B, which implies higher energy in the particles than in the magnetic field, at least locally, where the broadband blazar emission is produced (see e.g. Tavecchio et al. 2001b; Abdo et al. 2011a; Aleksić et al. 2015b). There is no physical reason for any specific (somewhat arbitrary) cut value in the quantity u_e/u_B; however, driven by previous works in the literature, in this study we only consider models fulfilling the requirement of u_e/u_B< 10³. Secondly, the observed variability timescales have to be taken into account. Following causality arguments, the observed variability should not happen on timescales that are shorter than the time needed to distribute information throughout the emitting region. Based on the given Doppler factor δ and the size of the emitting region R, the implied minimum variability timescale quantity for each model is derived according to (5)For the first flare we observed large flux changes (up to factors of ~4) within 0.5 h (Pichel & Paneque 2011; Aliu et al. 2016); instead, the second flare shows substantial flux changes (~2) on timescales of several days. Consequently, we consider only models that yield a minimum variability timescale of t_varmin ≤ 0.25 h and t_varmin ≤ 1 day for the first and second VHE flare, respectively.

4.2. First VHE flare

All spectral points that were obtained at the time or close to the time of the VHE flare measured by VERITAS and the Whipple 10 m telescope at MJD 54 952 are shown in the top panel of Fig. 9 (see Sects. 3.1.1 and 3.1.1 for details on the individual observation times).

The attempt to apply the grid-scan to the broadband SED from this flaring episode is affected by the flux variability on sub-hour timescales and the lack of strictly-simultaneous multiwavelength observations, as discussed in the previous sections. Therefore, the SED reported in this section is not necessarily a good representation of the true SED for this flaring episode, and hence any modeling results have to be regarded as inconclusive.

In this SED modeling exercise, the data used are the measurements from Swift (UV and X-rays), Fermi-LAT (two-day spectrum), and Whipple 10 m very high state. The optical and infrared, as well as the radio points, are not taken into account for the evaluation of the agreement of the SSC model curves with the data. The first two are strongly dominated by emission from the host galaxy, and the last only serve as upper limits for the SSC flux as the radio emission shows substantially lower variability timescales and is widely assumed to stem from a larger region than the emitting blob responsible for the few-day long flaring activity.

Exploiting the entire parameter grid space, neither the one-zone SSC model nor the two-zone SSC model can reconstruct the measured broadband SED, with the data-model agreement quantified with χ²/ d.o.f.> 300/20, which would imply a probability of agreement P (or p-value)¹² between the SSC model curves and the data points of P< 10^-50. When removing the tight constraint given by the cut in t_varmin, the best agreement obtained with the one-zone SSC scenario from the grid-scan defined by Table 5 is χ²/ d.o.f. = 180/20 (P ~ 10^-27). The two-zone scenario with the quiescent emission characterized by the model parameters from the average SED reported in Table 4 and the (spatially independent) region responsible for the flaring activity modeled based on the coarse grid parameter values shown in Table 6 provides at best an agreement given by χ²/ d.o.f. = 225/20 (P ~ 10^-36). Since the X-ray spectrum at low energies is already accounted for with the “quiescent” zone, the contribution from the “flaring” zone (which is needed to explain the increase in the flux at VHE) exceeds the measured flux at X-ray energies, and hence yields a bad agreement with the data points.

In addition to trying with the grid-scan defined in Table 5, we also evaluated the model-to-data agreement when using a one-zone scenario with the grid-scan defined in Table 6, which provides a simpler theoretical scenario (only one break in the EED instead of two), but with somewhat extended ranges probed for the parameters γ_min, γ_max, γ_break, n_e, and . We found a few models with data-model agreement given by χ²/ d.o.f. = 95/20 (P ~ 10^-11). But as soon as the requirement for fast variability is applied, all these models (mostly featuring large emission regions with R ≥ 10^16.5) are no longer applicable, and the agreement between the SSC model curves and the data points become χ²/ d.o.f.> 300/20.

One of the difficulties in modeling these data with a one-zone scenario is that it is difficult to describe the emission in the UV and the X-ray range with a synchrotron component. These UV flux points cannot be modeled only with the host galaxy template, and the one-zone models that could potentially describe well the shape of the X-ray spectrum would produce a flux that is many times below the measured UV flux, and hence would give a very bad data-model agreement. Contrary to the mentioned caveat of a time offset between the X-ray and VHE γ-ray observations, the UV and the X-ray observations were performed simultaneously and thus should be reconstructed consistently. The difficulty in modeling the UV and X-ray measurements in a consistent way suggests that a more complex scenario is needed to explain this emission. In Aliu et al. (2016), the host galaxy was modeled using a different template with respect to the one in Abdo et al. (2011a) that is used in this paper. The host galaxy template used in Aliu et al. (2016) describes approximately the measured UV flux level from the three-day broadband SED considered in Aliu et al. (2016), but it would not be consistent with the variability in the data set presented here. Figure 1 shows the relative variations in the UV flux of more than 50% (peak to peak), which cannot occur if this UV emission is dominated by the steady emission from the host galaxy.

4.3. Second VHE flare

The SED of Mrk 501 built from spectra around the time of the second flux increase seen by MAGIC on May 22 (MJD 54 973) is shown in the bottom panel of Fig. 9 (see Sects. 3.1.2 and 3.1.1 for details on the individual observation times). The data related to the second flare were not taken strictly simultaneously. However, here the resulting caveat is not as strong as for the first flare. On the one hand, the observed variability occurs on timescales of days, rather than tens of minutes, and the RXTE/PCA measurements were performed within a day of the VHE observations. While this is not true for the Swift/XRT measurements, the overall flux changes are relatively small, and the derived Swift/XRT spectrum is in very good agreement with the one derived from RXTE/PCA, as can be seen in the bottom panel of Fig. 9.

The results obtained for the one-zone scenario following the grid-scan from Table 5 gave a best probability of agreement with the data points of P_best ≈ 4 × 10^-10 (χ²/ d.o.f. ≈ 123/41). We found that there are 14 additional SSC model curves with a model-to-data probability higher than 0.1% of the best-matching model (i.e. P> 10^-3 × P_best), which we set as a generous probability threshold to consider the model-to-data agreement comparable. Given that P_best ≈ 4 × 10^-10, even those models with a best agreement of P> 10^-3 × P_best do not adequately describe the measured broadband SED, yet this relatively bad model-to-data agreement is not worse than some of the agreements between (simple) models and SED data shown in some studies (e.g. Abdo et al. 2010b; Giommi et al. 2012; Domínguez et al. 2013; Ahnen et al. 2016). This occurs because, in most studies involving broadband SEDs, the models are adjusted “by eye” to the data without any rigorous mathematical procedure that quantifies the model-to-data agreement. Differences on the order of 20−30% in a log-log plot spanning many orders of magnitude do not “appear to be problematic”, although these differences could be (statistically) significant owing to the small errors from some of the data points (e.g. optical/UV and X-ray). If the differences between the data and model are not substantial (regardless of the statistical agreement), the models are considered to approximately describe the data and can be used to extract some physical properties of the source and its environment.

Figure 10 depicts the best SSC model curves from the one-zone scenario, with the model featuring the best agreement to the data shown with a red curve, and the other 14 SSC models with comparable (down to 0.1%) model-to-data agreement shown with dark gray curves. Given the very low number of SSC model curves in this group, we decided to also depict those SSC models with model-to-data probability of agreement higher than 10^-6 × P_best and 10^-9 × P_best with lighter gray shades (see legend), which increased the number of SSC model curves depicted to 34. The thresholds used of 10^-6 × P_best and 10^-9 × P_best are somewhat arbitrary, and could be changed without any major qualitative impact on the reported results. The inclusion of these additional 20 models in the figure helps illustrate the behavior of the SSC model curves that start being worse than the best-matching model. To guide the eye, the SSC model describing the average state is also shown (from Abdo et al. 2011a, dash-dotted black line): the most significant deviations of the model curves from the data points stem from the Swift region. Therefore, while the hard X-ray and γ-ray bands can be satisfactorily modeled with a one-zone SSC scenario, this model realization fails at reconstructing both the soft X-ray data points and the UV emission at the same time. Figure 11 displays how many model curves produced for each point on the parameter grid yield a model-to-data agreement probability P better than 10^-3 × P_best, which are the models that are considered to be comparable. This is shown for each of the parameters separately. Some parameters are more constrained than others: e.g. γ_break,1, γ_break,2, and α₂ show a narrower distribution than, for instance, γ_max or α₃, which lead to equally good models over essentially the entire range of values probed. Additionally, as done for Fig. 10, with lighter gray shades we also report the parameter values for P> 10^-6 × P_best and P> 10^-9 × P_best. The SSC models that are not comparable to the best-matching model (i.e. those with P< 10^-3 × P_best) have a similar distribution for those parameters that are not constrained, like γ_max or α₃. On the other hand, of the parameters that can be constrained, like γ_break,1 and α₂, these additional models extend the range of parameter values with respect to the distributions for the models with P> 10^-3 × P_best. The parameter γ_break,2 seems to be quite well constrained, and even the models with P< 10^-3 × P_best converge to the same value of 3.2× 10⁵. The implications of these distributions on the possibility to constrain the different model parameters is discussed further in Sect. 5.

We also evaluated the model-to-data agreement for the one-zone scenario that uses the more simple grid-scan defined by Table 6, which is related to a grid of 9 parameters (instead of 11), but with a somewhat extended region for some of these parameters. We found that this grid-scan did not provide any additional SSC model with P> 10^-3 × P_best, and only five additional SSC models with P> 10^-9 × P_best. Hence this grid-scan did not bring any practical improvement with respect to that from Table 5, which led to 14 SSC models with P> 10^-3 × P_best and 34 SSC models with P> 10^-9 × P_best.

When using the above-mentioned two-zone SSC scenario, with the quiescent emission characterized by the model parameters from the average SED reported in Table 4 and the spatially independent region responsible for the flaring activity modeled based on the coarse grid parameter values reported in Table 6, we find a substantial improvement with respect to the one-zone models in describing the measured broadband SED (including the UV emission), with a best model-to-data probability of P_best ≈ 2.5 × 10^-3 (χ²/ d.o.f. ≈ 71/41). The two-zone model provides a better description of the SED than the one zone model, but it still does not reproduce the data perfectly.

Figure 12 displays the 69 SSC model curves with model-to-data agreement probability P better than 10^-3 × P_best, which is the generous probability threshold that we adopted to consider the probability of agreement comparable. Because of the relatively large number of SSC model curves (in comparison to those surviving the same selection criteria in the one-zone scenario), we decided to split those models into three groups according to their model-to-data probability P being better than 10^-1 × P_best, 10^-2 × P_best, and 10^-3 × P_best. Since P_best ≈0.25%, these models start providing an acceptable representation of the data, with the different bands reporting slightly different levels of success in the model-to-data agreement. The parameter values for these models are depicted in Fig. 13. As occurred for the one-zone scenario, some parameters are better constrained than others; for example, γ_break shows a narrow distribution, while γ_min and γ_max show a rather flat distribution. Although the parameter γ_min was probed up to 10⁶, the highest γ_min values used in the SSC models that can adequately describe the broadband SED only go up to 10⁴, which for the highest B field values reported in Fig. 13 (~0.15 G) relate to a cooling time of 3.5 × 10⁶ s. This is one order of magnitude longer than the dynamical timescale set by the highest R values reported in Fig. 13 (10¹⁶ cm), hence ensuring the existence of a low-energy cutoff. See Sect. 5 for further discussions of this topic.

In order to refine the adjustment of the different model parameters even further, a second iteration of the grid-scan modeling, referred to as a dense grid-scan, is performed. The dense grid-scan focuses on the parameter ranges that provide the best model-to-data agreement in the coarse grid-scan, which are depicted in Fig. 13. Following this strategy, the chosen parameter ranges can be narrowed in favor of a smaller step size in the individual parameters, while keeping the computing time at a reasonable amount. The new dense grid ranges and number of steps for each of the parameters are given in Table 7.

The model with the highest probability of model-to-data agreement in the dense grid-scan yields P_best ≈ 6.6 × 10^-2 (χ²/ d.o.f. ≈ 55.4/41), which implies an order of magnitude improvement with respect to the best-matching model obtained with the coarse grid-scan. If this model curve had been obtained through a regular mathematical fit, and conservatively considering that the nine dimensions of the grid relate to nine independent and free parameters in the fitting procedure, we would have obtained a p-value ≈ 6.3 × 10^-3 (χ²/ d.o.f. ≈ 55.4/32). The dense grid-scan focused on relatively good regions of the parameter space, which yielded a large number of SSC curves with a good model-to-data agreement although the parameter values of this dense grid-scan still vary widely (implying very different physical conditions in the source). Because of the large number of model curves, we can be more demanding with the probability threshold for considering the probability of agreement comparable to that of the best-matching model: a probability threshold of 0.1 × P_best still keeps 1684 SSC models, which is a large increase in statistics, in comparison to the results obtained with the coarse scan. Given that P_best ≈ 6.6%, all the models above this probability threshold provide a decent representation of the data. We split these models into three groups according to their model-to-data agreement being P> 0.9 × P_best, P> 0.5 × P_best, and P> 0.1 × P_best, hence reporting somewhat different levels of success in describing the measured broadband SED.

Here too we investigate the spread – or degeneracy – of the different models within the dense grid space of model parameters. Figure 14 shows again the distribution of the best model (red) and the models with P> 0.9 × P_best, P> 0.5 × P_best, and P> 0.1 × P_best over the entire dense grid parameter space. In comparison to Fig. 13, an apparent larger degree of degeneracy can be seen; the distributions have entries in most of the probed parameter ranges depicted in the figure. The wider spread in the parameter values shown in Fig. 14 is caused by the selected parameter range, which is narrower and intentionally only covers regions with an already reasonable agreement between model and data, as derived from the coarse grid-scan. Despite the large spread, one can see that there are regions with slightly better model-to-data agreement, like the region around γ_break ~ 5 × 10⁵ or the region around α₁ ~ 1.9. The results are discussed further in Sect. 5.

The SED models that were picked as a result of the dense grid-scan for two-zone SSC models are presented in Fig. 15. The figure highlights three SSC model curves: the model that gave the best agreement with the SED data points, a model featuring a prominent high-energy component in the EED, and a model that features a low Doppler factor of δ = 5. The parameter values for these three specific SSC model curves are given in Table 8, showing once more that three very distinct sets of SSC model parameters can provide comparable agreement with the experimental data.

5. Discussion

5.1. Variability and correlations

For Mrk 501, an increase in the fractional variability with energy has been reported in the past within the X-ray and VHE band (Gliozzi et al. 2006; Albert et al. 2007a; Aliu et al. 2016). In the work presented here, we extend this trend throughout all wavelengths from radio to VHE γ rays, showing that the source is relatively steady at radio/optical frequencies, but variable (F_var ≥ 0.2) and very variable (F_var ≥ 0.4) in the X-ray and VHE γ-ray bands, respectively, with a clear increase in the fractional variability with energy (observed in all the bands we can measure). A similar variability pattern was reported in Aleksić et al. (2015b) and, during the preparation of this study, also in Ansoldi et al. (in prep.) in relation to the extensive campaigns on Mrk 501 performed in 2008 and 2012, respectively. This suggests that this variability vs. energy behavior is an intrinsic characteristic in Mrk 501. On the other hand, Furniss et al. (2015) has recently reported a different fractional variability vs. energy pattern based on observations taken in 2013 where the observed variability at X-rays is similar to that at VHE.

The multiband variability pattern that has been observed in Mrk 501 is quite different from that observed in Mrk 421 during the multi-instrument campaigns from 2009, 2010, and 2013, as reported in Aleksić et al. (2015a, (2015c) Baloković et al. (2016). In these works a double-bump structure in the fractional variability plot was found (instead of a continuous increase with energy) which relates to the two bumps in the broadband SED, and where the highest variability occurs at X-rays and VHE at comparable levels.

A clear correlation of the X-ray and VHE γ-ray emission was observed during the large and long γ-ray activity from 1997 (e.g. Pian et al. 1998; Gliozzi et al. 2006), but this correlation was only marginally detected during the γ-ray flare observed in 2005 (Albert et al. 2007a). The low significance in the X-ray-to-VHE correlation during the flares in 2005 was ascribed to the lack of sensitive X-ray measurements during this observing campaign; only RXTE/ASM data, which has limited sensitivity to detect Mrk 501, was available for this study. A positive correlation between X-ray and VHE γ rays was also reported – for the first time – during very low X-ray and VHE activity, but only at a 99% confidence level (Aleksić et al. 2015b). The marginally significant correlation observed during this low activity, using data from the multi-instrument campaign in 2008, was ascribed to the very low variability during that campaign, where the measured F_var values were about 0.1 for X-rays and 0.2 for VHE. As reported above, during the multi-instrument campaign in 2009, Mrk 501 was mostly in its low/typical state, but we also measured two flaring activities in May 2009. The measured F_var is about 0.3 for X-rays and 0.8 for VHE, while if we exclude the two flaring episodes we obtain F_var values of about 0.2 for X-rays and 0.5 for VHE. However, despite the larger variability observed in 2009 (with respect to 2008), we did not observe any significant correlation between the X-ray and the VHE emission (including and excluding the flaring episodes). This may appear to be a controversial result, but we would like to stress that a very significant correlation with past data was only observed during the very large and long flare in 1997. Recently, Furniss et al. (2015) reported a significant X-ray-to-VHE correlation using data from the multi-instrument campaign in 2013. This correlation is dominated by the large X-ray and VHE activity observed over four consecutive days in July 2013, although it still remains at 2 σ (for the 0.3−3 keV energy band) and 5σ (for the 3−7 keV band) when removing the flaring activity. In conclusion, some multi-instrument campaigns on Mrk 501 do not show a clear X-ray-to-VHE correlation when the source is not flaring strongly or persistently high. However, for the other archetypal TeV blazar, Mrk 421, the X-ray-to-VHE correlation is significantly detected during both low- (e.g. Aleksić et al. 2015a; Baloković et al. 2016) and high-activity states (e.g. Fossati et al. 2008; Acciari et al. 2011a; Aleksić et al. 2015c).

The X-ray-to-VHE correlation and the fractional variability vs. energy pattern observed in Mrk 421 suggests that the X-ray and VHE emissions are produced by the same electrons within the framework of SSC scenarios, and that the highest variability is produced by the highest energy and most-variable electrons that dominate the emission at the keV and the TeV bands, respectively. Instead, in Mrk 501 we observe a continuous increase in the variability with energy and absence of persistent correlation between the keV and TeV emissions. This suggests that the highest energy electrons, in the framework of SSC scenarios, are not responsible for the keV emission, while they are responsible (at least partially) for the TeV emission. Alternatively, there could be an additional (and very variable) component contributing to the γ-ray emission in addition to that coming from the SSC scenario, like inverse-Compton of the high-energy electrons off some external low-energy photon fields (Dermer et al. 1992; Sikora et al. 1994; Finke 2016).

5.2. VHE flaring state SEDs

The first flaring event (MJD 54 952) is characterized by a fast and large outburst in the VHE band, which was apparently not accompanied by a substantial increase in the X-ray flux, and hence appeared to be like an “orphan flare” (see e.g. Krawczynski et al. 2004). In fact, based on these observations, this event was tentatively categorized as an orphan flare event (Pichel & Paneque 2011; Neronov et al. 2012), which would substantially challenge the currently favored SSC emission models (for HBLs). However, a detailed look at the SED of the flaring episode reveals a hardening of the X-ray spectrum measured by Swift/XRT (see Sect. 3), which more likely corresponds to a shift of the synchrotron bump towards higher energies. During the outstanding activity in 1997, the synchrotron peak was shifted to beyond 100 keV, as accurately measured by BeppoSax (Pian et al. 1998). Such a large increase in the location of the synchrotron peak position could have occurred in the MJD 54 952 flare discussed here. Additionally, the peak of the high-energy γ-ray bump at the time of this flare also appears to shift towards higher energies, as occurred in 1997. This suggests a more general appearance of such phenomena, and that – even though the measured keV and TeV flux are not correlated during this flaring activity – the overall broadband X-ray and VHE emission may still be correlated, which could have been measured if X-ray observations at several tens of keV had been available during this flaring episode. Such a shift of the entire SED has been interpreted as a shift in the energy distribution of the radiating electron population (e.g. Pian et al. 1998; Albert et al. 2007a). In this context, the small change in the inverse Compton peak position compared to that of the synchrotron peak location could be ascribed to Klein-Nishina effects. High-energy electrons can efficiently produce high-energy synchrotron photons; however, they are not able to upscatter photons as effectively as the lower-energy electrons because the Klein-Nishina cross-section is smaller than the Thomson cross-section (Tavecchio et al. 1998; Acciari et al. 2011b).

We tried to parameterize the broadband SED during this first flare (MJD 54 952) using a wide range of SSC emission scenarios following the grid-scan strategy defined in Sect. 4.1, allowing for models with one or two (independent) emission zones and covering a wide range in the space of model parameters. We found that none of the tested models could satisfactorily reproduce the changes observed in the spectral distribution. This broadband SED can be explained with more sophisticated theoretical models, like the inhomogeneous time-dependent models reported in Ghisellini et al. (2005), Graff et al. (2008), Chen et al. (20112015, 2016), which provide a more elaborate physical scenario, at the expense of an increase in the number of degrees of freedom of the model. However, a caveat has to be taken into account when interpreting these results, which is the lack of strict simultaneity of the different data sets, in particular the Swift/XRT and the VHE γ-ray data. The individual exposure times are separated by seven hours, while we see flux changes of a factor of ~4 on sub-hour timescales, as reported in Pichel & Paneque (2011) and Aliu et al. (2016). Therefore, it is somewhat uncertain whether the measurements of the synchrotron peak and the high-energy peak probe the same source state. In the recent study reported in Aliu et al. (2016), the broadband SED derived for the three-day time interval MJD 54 952–55 could be satisfactorily parameterized with a one-zone SSC scenario that differs from the one used here and in Abdo et al. (2011a) in various aspects, including the template used to describe the host galaxy contribution.

Compared to the first flaring event, the second flare (≈MJD 54 973) occurs during VHE flux changes of factors of ~2 on timescales of a few days, and hence the lack of strict simultaneity in the X-ray/VHE observations is a much smaller caveat than for the first flaring event. In this case, again following the grid-scan modeling approach, one-zone SSC models were unable to describe the measured SED (reaching best probabilities of agreement ~ 10^-10). The two-zone SSC models were able to reproduce the experimental observations better (reaching best probabilities of agreement ~ 10^-3). Therefore, the two-zone scenario appears to be favored compared to the one-zone scenario considered here. Building on the range of two-zone model parameters providing decent data-model agreement, a fine grid-scan was performed, yielding hundreds of two-zone SSC models with probabilities of agreement ~ 10^-2. The obtained set of two-zone SSC models providing the best agreement comprises several setups with quite different implications for the parameters defining the EED and the surrounding region of the second emission region (see Sect. 4.3).

Comparing the configurations obtained for the emission region responsible for the second flare with the parameter values describing the emission region assumed to create the quiescent emission, some general trends can be stated: while the parameters describing the EED and the Doppler factor are found to populate roughly the same ranges of values, the electron density n_e is increased by 1−2 orders of magnitude, the magnetic field is larger by ≈ 1 order of magnitude, and the size of the emission region R is found to be smaller by 1−2 orders of magnitude. The last result is affected by the requirement of a minimum variability timescale of a day in order to account for the variability seen in the data.

In addition to the general observations made above, some interesting model configurations stand out from the set of adequate scenarios: models that feature a prominent high-energy component in the EED and models with Doppler factors as low as δ = 5 can be used to adequately model the flaring SED. In the paragraphs below we discuss the benefits of these two families of models.

Synchrotron self-Compton models with a strong high-energy component are interesting not only to explain the SED collected during the presented campaign, but also in the context of other observations of Mrk 501. During the extreme flare seen in 1997, a strong increase in the regime of hard X-rays, around 100 keV, was observed (Bradbury et al. 1997). This increase can be interpreted as the emergence of a strong high-energy component adding to the overall SED, which only sometimes becomes visible during extreme flaring states. Moreover, Cherenkov telescope observations often give hints of an additional hard component in the EED during flaring times: in Albert et al. (2007a) a significant spectral hardening during flaring states was observed and reported for the first time, and in the course of several more observational campaigns this “harder when brighter” behavior has been established as typical for Mrk 501. Ultimately, a tendency for this behavior was also seen during the campaign presented in this paper. In this light, SSC models with such a high-energy contribution to the EED could be favored as they can also explain such mentioned observations. Naturally, with the data set at hand and the lack of hard X-ray observations above ≈ 20 keV, they can neither be confirmed nor discarded.

The finding that models with δ = 5 can also adequately reconstruct the data is particularly interesting. Quite high values (above δ = 10, up to δ = 50 or more) are usually required to model the SEDs of blazars (Tavecchio et al. 1998; Krawczynski et al. 2001; Saugé & Henri 2004). These high Doppler factors are in tension with regard to expectations from the small (typically less than 2c) apparent velocities observed in the 43 GHz radio emission of various high-peaked BL Lac objects, and particularly with that measured for Mrk 501 (Edwards & Piner 2002; Piner & Edwards 2004; Piner et al. 2010). This has posed a common problem for TeV sources, which has been dubbed the “bulk Lorentz factor crisis” (Henri & Saugé 2006), and requires the radio and TeV emission to be produced in regions with different bulk Lorentz factors. Debates on this problem (see e.g. Georganopoulos & Kazanas 2003; Levinson 2007; Stern & Poutanen 2008) have led to a series of sophisticated models, for example the “spine-sheath” model from Ghisellini et al. (2005), in which the jet is structured transverse to its axis into a fast “spine” that produces the high-energy emission, and a slower “layer” which dominates the radio emission. The modeling results presented in this paper show that it is actually possible to model the SED of a flaring activity of Mrk 501 using a relatively small Doppler factor to describe the flaring emission, hence alleviating the tension with the radio interferometric observations. Naturally, such low Doppler factors cannot be used for broadband SEDs related to periods with fast (sub-hour) variability, such as the one from MJD 54 952, but they can be be used for broadband SEDs related to flaring activities with day timescales, such as the one from MJD 54 973, which are more commonly observed in high-peaked BL Lac objects.

In addition to individual models which give a good reproduction of the data, the degree of degeneracy of well-fitting models in the individual parameters has also been studied, revealing a wide range of equally good models in the SSC parameter space, and showing that some model parameters can be constrained much better than others. While this can be seen already for the one-zone SSC models, where γ_break,1 and γ_break,2 show a narrower distribution than γ_max or α₃, it is particularly interesting to study this for the more applicable two-zone scenario, which is the one that suitably describes the data. We find that for both the coarse and the dense grid-scan, the distribution of parameter values giving a good agreement with the data is quite well constrained for some parameters, such as the Lorentz factor at the break energy of the electrons γ_break, while other parameters show a rather broad distribution, like γ_max or the index of the EED after the break α₂, which points to a real degeneracy in these parameters. To some extent this degeneracy can be explained by the unequal sampling of the SED: the density and accuracy of measurements at or around the positions of the synchrotron and the IC peak is rather dense, which leads to a good definition of the spectral break in the EED. However, moving from the peak positions up to higher energies, the uncertainties of the measurements increase (especially for the synchrotron peak) and parameters such as the spectral index after the break α₂ or the Lorentz factor where the EED is cut off γ_max cannot be constrained equally well.

This result has several implications for the modeling of SEDs in general. On the one hand, it shows that an actual fitting procedure, which moves along the direction of the steepest gradient in the parameter space towards a minimum in the χ² of the model-to-data agreement, does not necessarily reveal the entire picture of possible descriptions of the data in the context of the applied model. Usually one best solution is quoted as the result, while most of the time a wide range of models explain the data equally well. We also see that in order to be able to put stronger constraints on the parameters defining SED models, we need data sets that are characterized by a better coverage in energy and by smaller uncertainties in flux. We see in Fig. 15 that especially the hard X-ray regime, but also the HE and VHE γ-ray regime, allow for a wide range of possible model curves.

Unfortunately, the exercise presented here indicates that the SED modeling results that are performed for less constrained data sets, e.g. for “weak sources” that are sampled with a smaller energy coverage and with spectral measurements with larger statistical uncertainties, should be taken with caution because they are likely to have substantial degeneracies in the model parameters. Such modeling exercises can demonstrate that a particular scenario (e.g. one-zone or two-zone SSC) is capable of reproducing the measured data, but they certainly cannot claim the exclusiveness or even the prominence of the particular set of parameter values that has been chosen or found to be “best”.

5.3. Change in optical polarization during VHE γ-ray flare

The first VHE flare (MJD 54 952) was found to coincide with an observed change in the optical polarization. While simulations of turbulent processes in blazar jets show that a rotation of this dimension can be ascribed to random behavior (Marscher 2014), the coinciding occurrence of the change in rotation and a flare of the VHE γ-ray flux suggests a common origin of these events. Such combined events have already been seen in low-frequency peaked BL Lac-type sources (LBL) and flat spectrum radio quasars (FSRQ), but it was observed for the first time for an HBL in the course of the 2009 campaign, and already reported in Pichel & Paneque (2011) and Aliu et al. (2016).

These observations show similarities to double or multiple flaring events seen in the LBL BL Lacertae in 2005 and in the FSRQ PKS 1510-089 in 2009, which were discussed by Marscher et al. (2008) and Marscher et al. (2010), respectively.

Exhibiting different peak frequencies for the synchrotron and the IC bump, the optical variability seen in BL Lac could be seen as corresponding to the X-ray variability in Mrk 501. While a strong flare in the VHE band was observed during the first flare of Mrk 501, BL Lac gives hints for activity in that band during the first optical outburst (Albert et al. 2007b). A coincidence of a flaring event and a change in the optical polarization is seen in all three data sets. The observed degree in optical polarization in Mrk 501 of ≈ 5% appears smaller than that in BL Lacertae (up to 18%). Still, the optical flux in Mrk 501 is strongly dominated by the host galaxy, so that the jet contribution amounts to only ~1/3. Therefore, the measured degree of polarized light in Mrk 501 corresponds to a fraction of ≈ 15−20% of polarized emission from the jet, which is comparable to BL Lacertae. The second episode of high activity in both Mrk 501 and BL Lac was characterized by an increased flux at the synchrotron bump over a longer time span.

In the case of BL Lac, Marscher et al. (2008) suggested that the first flare and the change in polarization may have occurred when a blob of highly energetic particles traveled along the last spiral arm of a helical path within the acceleration and collimation zone of the jet, and finally left this zone to enter a more turbulent region. The second flare seen in BL Lac has been identified with the passage of the feature through the shocked region of the radio core. The observed behavior of Mrk 501 suggests that the discussed scenario could be applicable here. Despite the lack of simultaneous interferometric radio observations during the two flares, an enhancement of the activity in the VLBA 43 GHz core emission in May 2009 (with respect to the previous months) was observed, supporting the interpretation of the blob traversing a standing shock region during the second flaring episode.

The polarization data collected during this campaign could also hint at a different physical scenario. After the VHE flare on MJD 54 952, not only did the EVPA rotation stop, but it also started rotating in the reverse direction during the following two days, and the polarization degree did not drop to the “typical low-state values” of about 1−2%, but only decreased from 5.4% to 4.5% (see the bottom two panels of Fig. 1). The characteristics of the polarization data after the VHE flare may be better explained as resulting from light-travel-time effects in a straight shock-in-jet model with helical magnetic fields, as proposed by Zhang et al. (2014, 2015). This shock-in-jet model, which uses a full three-dimensional radiation transfer code and takes into account all light-travel-time and other geometric effects (for some assumed geometries), may be more successful in explaining the broadband SED (and variability patterns) observed during the VHE flare from MJD 54 952, which could not be explained with the relatively simple one-zone and two-zone SSC scenarios described in Sect. 4. However, the lack of strictly simultaneous X-ray/VHE data during the MJD 54 952 VHE flare, and the relatively scarce polarization observations after the VHE flare would be an important limitation in the full application of this theoretical scenario to the multi-instrument data set presented here.

6. Summary and concluding remarks

We presented a detailed study of the MWL variability of the HBL Mrk 501, based on a multi-instrument campaign that was conducted over 4.5 months in 2009, with the participation of MAGIC, VERITAS, the Whipple 10 m, Fermi-LAT, RXTE, Swift, GASP-WEBT, and several optical and radio telescopes. Mrk 501 shows an increase in the fractional variability with energy, from a steady flux at radio and optical frequencies to fast and prominent flux changes in the VHE γ-ray band. Overall, no significant correlation was found between any of the measured energy bands, particularly no correlation was seen between X-rays and VHE γ rays despite the relatively large variability measured in these two energy bands. This suggests that the highest energy (and most variable) electrons that are responsible for the VHE γ rays measured by MAGIC, VERITAS, and Whipple 10 m do not have a dominant contribution to the ~1 keV emission measured by Swift/XRT. These high-energy electrons may have a dominant contribution to the hard X-ray emission above 10–50 keV where the instrumentation used in this campaign did not provide sensitive data. Alternatively, there could be a component contributing to the VHE γ-ray emission in addition to the component coming from the SSC scenario (e.g. external Compton), which is highly variable and further increases the variability of Mrk 501 at VHE γ rays with respect to that expected from the pure SSC scenario.

This paper discusses two prominent flaring events at VHE γ rays with different characteristics that were seen during the campaign. The first flare is dominated by a fast outburst in the VHE range, which does not appear to be accompanied by a large flux increase in the X-ray band, but shows a hardening in the X-ray spectrum that can be associated with a shift of the synchrotron bump to higher energies. On the other hand, the second flare is characterized by a flux increase in both the VHE and the X-ray band. For the parameterization of the broadband SEDs from these two VHE γ-ray flares, we applied a novel variation of the grid-scan approach in the space of model parameters. For the two theoretical scenarios investigated, the one-zone and two-zone SSC models, we probed multi-dimensional grids with the various model parameters, evaluating the model-to-data agreement for tens of millions of SSC models. This strategy allowed us to identify disjointed regions of equally good model configurations, and provided a quantification of the degeneracy in the model parameters that describe the measured broadband SEDs. The presented methodology provides a less biased interpretation than the commonly used “single-curve model adjustment procedure” typically reported in the literature.

A lack of strict simultaneity in the X-ray/VHE observations of the first flare, which is characterized by large VHE flux changes in sub-hour timescales, does not permit us to draw final conclusions on the underlying mechanism; but the SED modeling with the grid-scan suggests that a simple one-zone or two-independent-zone SSC model is not sufficient to explain the measured broadband emission. The broadband SED derived for the second flare also lacks strictly simultaneous observations, but the flux changes here are smaller and on longer timescales, and hence substantially less problematic than for the first flare. The overall SED from the second flare cannot be properly described by a one-zone SSC model (with an EED with two spectral breaks), while it can be reproduced satisfactorily within a two-independent-zone SSC scenario. In the two-zone models applied here, one zone is responsible for the quiescent emission from the averaged 4.5-month observing period, while the other one, which is spatially separated from the first, dominates the flaring emission occurring at X-rays and VHE γ rays. The grid-scan shows that there is a large number of SSC model realizations that describe the broadband SED data equally well, and hence that there is substantial degeneracy in the model parameters despite the relatively well-measured broadband SEDs. For instance, regarding the features of the EED, the position of the break(s) appear to be well constrained, while the highest Lorentz factor and the high-energy spectral index vary more strongly within the best-fitting model realizations. While the few models with the best relative agreement to the data feature Doppler factors δ in the range 10−20, the data can also be reproduced using substantially lower Doppler factors of δ = 5 while still reaching fit probabilities higher than 10% P_best. This shows that it is possible to reproduce the observed SED from Mrk 501 assuming boost factors well below the usually required values of δ ≈ 10−50, which may loosen a bit the tension posed between large values of δ required for modeling and low values imposed from radio velocity measurements, which has been dubbed the “bulk Lorentz factor crisis”.

A change in the rotation of the EVPA was measured in temporal coincidence with the first VHE flare, at MJD 54 952, as reported in Pichel & Paneque (2011) and Aliu et al. (2016). Here we also show that during the first VHE flare, the degree of polarization increased by a factor of ~3 with respect to the polarization measured before and after this flaring activity. This is the first time that this behavior was observed in Mrk 501, or in any other HBL object, and suggests a common origin of the VHE flare and the optical polarized emission. With the coincidence seen of a VHE flare and changes in the trend of the optical polarization, this two-flare event resembles prior events observed in the LBL BL Lacertae and the FSRQ PKS 1510-089, which were discussed in Marscher et al. (2008) and Marscher et al. (2010), respectively. The common features suggest a similar interpretation of the flaring event of Mrk 501 as an emission region which is traveling upstream of the radio core along a spiral path in a helical magnetic field, entering a region of turbulent plasma and crossing the standing shock region of the radio core. After the VHE flare at MJD 54 952, the polarization degree decreased from 5.4% to only 4.5% (instead of the typical low-state value of 1–2%), and there is also a small EVPA rotation in the reverse direction during the following two days, which may be difficult to explain with the typical helical pattern motions mentioned above; they may be better explained as light-travel-time effects in a shock-in-jet model in a straight, axisymmetric jet embedded in a helical magnetic field, as reported in Zhang et al. (2014, 2015). Beyond the interpretation of the flaring event itself, the observational results obtained in the course of this MWL campaign reveal phenomena that have not been seen for any HBL before, but have already been studied for LBLs and FSRQs. This gives a strong indication of the intrinsic similarity of these blazar subclasses, even though they show different jet characteristics in general, such as apparent jet speed and the overall power output. Observations of rapid variability in the LBL BL Lacertae (Arlen et al. 2013) support this further, as such fast flux changes had only been seen in HBL observations.

Additional multi-instrument observations of Mrk 501 will be crucial to confirm and extend several of the observations presented here. First of all, the large degeneracy in the model parameter values providing an acceptable description of the broadband SED is largely dominated by the poor coverage at hard X-rays, as well as the somewhat limited resolution at VHE γ rays. Since mid-2012, NuSTAR¹³ provides excellent sensitivity above 10 keV, which narrows the large gap in the SED between soft X-rays and the Fermi-LAT regime for campaigns from 2012 onwards. In the coming years, observations with ASTROSAT¹⁴ will also be possible, hence facilitating a much more accurate characterization of the evolution of the synchrotron bump of Mrk 501, including the determination of hard components in the EED whose synchrotron emission may peak at 50 or 100 keV. Moreover, in the regime of VHE γ rays, data sets of much higher quality are already being collected, yielding a better resolution and an extended energy coverage. This is achieved on the one hand by the operation of the MAGIC telescopes as a stereo system, which gave a remarkable improvement in the overall performance compared to the mono mode which was still operational during the campaign reported in this paper (Aleksić et al. 2012). On the other hand, both MAGIC and VERITAS underwent major upgrades in the years 2011 and 2012, which gave a further substantial push to the performance (Aleksić et al. 2016a,b; Zitzer et al. 2013). In the future, the Cherenkov Telescope Array (CTA) promises to deliver further substantial improvement both in terms of the energy coverage and the resolution of the flux measurement (Actis et al. 2011). Additionally, the temporal coverage extending over many years will permit variability studies (including many flares) over large portions of the electromagnetic spectrum and with good sensitivity, which will permit an evaluation of whether the association of EVPA rotations and polarization degree changes with VHE γ-ray flares are rare or regular events, whether these events occur together with a second flaring activity with contemporaneous enhancement of the VLBA radio core emission, and whether the measured multiband variability and lack of 1 keV–1 TeV correlation is a typical characteristic in Mrk 501 that is repeated over time.

---

Acknowledgments

The authors thank the anonymous referee for providing a very detailed and constructive list of remarks that helped us to clarify and improve some of the results reported in the manuscript. 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. VERITAS is supported by grants from the US Department of Energy Office of Science, the US National Science Foundation and the Smithsonian Institution, by NSERC in Canada, by Science Foundation Ireland (SFI 10/RFP/AST2748), and by STFC in the UK The VERITAS collaboration acknowledges the excellent work of the technical support staff at the Fred Lawrence Whipple Observatory and at the collaborating institutions in the construction and operation of the instrument. The Fermi LAT Collaboration acknowledges support from a number of agencies and institutes for both development and the operation of the LAT as well as scientific data analysis. These include NASA and DOE in the United States; CEA/Irfu and IN2P3/CNRS in France; ASI and INFN in Italy; MEXT, KEK, and JAXA in Japan; and the K. A. Wallenberg Foundation, the Swedish Research Council, and the National Space Board in Sweden. Additional support from INAF in Italy and CNES in France for science analysis during the operations phase is also gratefully acknowledged. We acknowledge the use of public data from the Swift and RXTE data archive. The Metsähovi team acknowledges the support from the Academy of Finland to our observing projects (numbers 212656, 210338, 121148, and others). This research has made use of data obtained from the National Radio Astronomy Observatory’s Very Long Baseline Array (VLBA), projects BK150, BP143, and MOJAVE. The St. Petersburg University team acknowledges support from Russian RFBR grant 15-02-00949 and St. Petersburg University research grant 6.38.335.2015. AZT-24 observations are made within an agreement between the Pulkovo, Rome, and Teramo observatories. The Abastumani Observatory team acknowledges financial support by the Shota Rustaveli National Science Foundation under contract FR/577/6-320/13. This research is partly based on observations with the 100 m telescope of the MPIfR (Max-Planck-Institut für Radioastronomie) at Effelsberg, as well as with the Medicina and Noto telescopes operated by INAF–Istituto di Radioastronomia. The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica. The OVRO 40 m program was funded in part by NASA Fermi Guest Investigator grant NNX08AW31G, and the NSF grant AST-0808050, and the Steward observatory by the NASA Fermi Guest Investigator grant NNX08AW56G.

References

All Tables

All Figures