FIRST NuSTAR OBSERVATIONS OF MRK 501 WITHIN A RADIO TO TeV MULTI-INSTRUMENT CAMPAIGN

Markarian 501 (Mrk 501) is a nearby, bright X-ray-emitting blazar at $z=0.034,$ also known to emit very-high-energy (VHE; $E\geqslant 100$ GeV) gamma-ray photons (Quinn et al. 1996). Blazars are among the most extreme astrophysical sources, displaying highly variable emission at nearly every wavelength and timescale probed thus far. These objects are understood to be active galactic nuclei that are powered by accretion onto supermassive black holes and have relativistic jets pointed along the Earth's line of sight (Urry & Padovani 1995). Relativistic charged particles within blazar jets are responsible for the non-thermal spectral energy distribution (SED), which is characterized by two broad peaks in the $\nu {F}_{\nu }$ spectral representation. The origin of the lower-energy peak is relatively well understood, resulting from the synchrotron radiation of relativistic leptons in the presence of a tangled magnetic field (Marscher 2008). Within the leptonic paradigm, the higher-energy SED peak is attributed to inverse-Compton up-scattering by the relativistic leptons within the jet of either the synchrotron photons themselves, namely synchrotron self-Compton (SSC) emission (Maraschi et al. 1992), or a photon field external to the jet, namely external Compton emission (e.g., Dermer et al. 1992; Sikora et al. 1994). Alternatively, hadronic models attribute the higher-energy peak of blazar emission to proton synchrotron emission and/or synchrotron emission by secondary leptons produced in p–γ interactions (Bednarek 1993; Aharonian et al. 2002).

Along with the other nearby VHE blazar Mrk 421, Mrk 501 represents one of the most comprehensively studied VHE blazars. The blazar has been the subject of multiple broadband observation campaigns (e.g., Catanese et al. 1997; Kataoka et al. 1999; Petry et al. 2000; Abdo et al. 2011a). Mrk 501 is one of the brightest X-ray sources in the sky, and has been observed by RXTE to display significant X-ray variability up to 20 keV (Gliozzi et al. 2006). During a phase of high activity at VHE energies in 1997, this source was also observed by BeppoSAX to display unusually hard, correlated X-ray emission up to $\gt 100\;{\rm{keV}},$ with a photon index of ${\rm{\Gamma }}\lt 2$ (Pian et al. 1998).

Observations of Mrk 501 have so far lacked sufficient sensitivity at the hard X-ray energies (10–100 keV). Observations at hard X-ray energies provide direct insight into the highest energy particles through detection of synchrotron emission. There is also the possibility for insight into the lower energy particles through the detection of inverse-Compton emission from photon up-scattering by the lower-energy electrons. As a relativistic synchrotron emitter, the falling edge of the synchrotron peak mimics the energy distribution of the emitting particles, allowing the highest energy particles to be directly probed through hard X-ray observations. The energy-dependent cooling timescale can lead to more rapid variability at hard X-ray energies than at soft X-ray energies. Gliozzi et al. (2006) reported independent soft (2–10 keV) and hard (10–20 keV) X-ray variability of Mrk 501 using RXTE.

Other hard X-ray observations have previously been performed with BeppoSAX (Massaro et al. 2004a) and Suzaku HXD (Anderhub et al. 2009). Due to the rapid X-ray variability displayed by blazars such as Mrk 501, the long integration time required for significant detection and spectral reconstruction by the aforementioned X-ray instruments was not ideal for extracting information about hard X-ray variability. Much more sensitive hard X-ray observations of blazars, however, are now possible with Nuclear Spectroscopic Telescope Array NuSTAR.

NuSTAR is a hard X-ray (3–79 keV) observatory launched into a low Earth orbit in 2012 June (Harrison et al. 2013). It features the first focusing hard X-ray telescope (XRT) in orbit that allows high sensitivity beyond the 10 keV cutoff shared by all other currently active focusing soft X-ray telescopes. The inherently low background associated with concentrating the X-ray light enables NuSTAR to achieve approximately a one-hundred-fold improvement in sensitivity over the collimated and coded-mask instruments that operate in the same spectral range.

NuSTAR observed Mrk 501 four times in 2013 as part of a simultaneous multiwavelength (MWL) campaign, including VHE observations by MAGIC and VERITAS, high-energy (HE; 100 MeV−100 GeV) gamma-ray observations by the Fermi Large Area Telescope (LAT), soft X-ray and UV observations with Swift X-ray Telescope (XRT) and Ultraviolet Optical Telescope (UVOT), optical observations from a number of ground-based instruments including the GASP-WEBT program, as well as radio observations by the Owens Valley Radio Observatory (OVRO; 15 GHz), Metsähovi (37 GHz), and the F-Gamma monitoring program, providing measurements between 2.64 and 228.39 GHz. The NuSTAR observations took place on 2013 April 13, May 8, and July 12 and 13 (MJD 56395, 56420, 56485, and 56486, respectively), with the latter two observations resulting from target of opportunity (ToO) exposures triggered by an elevated state observed by the Swift XRT and the MAGIC telescopes.

We use these observations to study the hard X-ray spectral behavior of Mrk 501 in detail over multiple flux states. The NuSTAR observations, analysis, and results are detailed in Section 2, with the contemporaneous MWL observations, analysis and results shared in Section 3. After comparing the simultaneous Swift XRT and NuSTAR observations in Section 4, we investigate variability of the source in Section 5. The MWL SEDs are constructed over the multiple observed states and investigated in terms of a single-zone equilibrium SSC model in Section 6, with discussion and conclusions provided in Section 7.

In order to maximize the strictly simultaneous overlap of observations by NuSTAR and ground-based VHE observatories during this broadband campaign of Mrk 501, the observations were arranged according to visibility of the blazar at the MAGIC and VERITAS sites. The NuSTAR coordinated observations involving both VERITAS and MAGIC were performed on 2013 April 13 and May 8, with the NuSTAR ToO observations (initiated by Swift and MAGIC) performed on 2013 July 12 and 13. The NuSTAR observations typically spanned 10 hr, resulting in 10–30 ks of source exposure after removing periods of orbital non-visibility. The observation details are summarized in Table 1. The data were reduced using the standard NuSTARDAS software package⁹⁵ v1.3.1.

The spectral analysis was performed with XSPEC ⁹⁶ Version 12.7.1. The data were binned to require 20 counts per bin, and fit with three spectral models via ${\chi }^{2}$ minimization. The first model applied to the data is a power law

$$\begin{eqnarray}&&A{(E)}_{{\rm{PL}}}=K{(E/{E}_{0})}^{-{\rm{\Gamma }}},\end{eqnarray} \tag{ 1 }$$

referred to as the PL model for the remainder of this work, where F(E) is the flux at energy E, ${\rm{\Gamma }}$ is the index, K is the normalization parameter (in units of photons keV⁻¹ cm⁻² s⁻¹) and E₀ is fixed at 10 keV.

The second spectral model applied to the data is a broken power law, referred to as BKNPL model for the remainder of this work. The model is made up of two power-law photon indices, meeting at a break energy E_break

$$\begin{eqnarray}&&A{(E)}_{{\rm{BKNPL}}}=K{(E/{E}_{\mathrm{break}})}^{-{{\rm{\Gamma }}}_{\mathrm{1,2}}}\end{eqnarray} \tag{ 2 }$$

where ${{\rm{\Gamma }}}_{1}$ and ${{\rm{\Gamma }}}_{2}$ represent the photon indices below and above the break energy E_break, respectively.

The third spectral model applied to the data is a log parabola, referred to as the LP model for the remainder of this work. This model has been suggested to better represent the X-ray spectra of TeV-detected blazars between 0.2 and 100 keV (e.g., Massaro et al. 2004b; Tramacere et al. 2007). This model allows the spectral index to vary as a function of energy according to the expression

$$\begin{eqnarray}&&A{(E)}_{{\rm{LP}}}=K{(E/{E}_{0})}^{-({\rm{\Gamma }}+\beta {\rm{log}}(E/{E}_{0}))},\end{eqnarray} \tag{ 3 }$$

with a curvature parameter β. The spectral data, model fits, and data-to-model ratios for each NuSTAR observation are shown in Figure 1. The spectral fitting results for each model as applied to the NuSTAR observations are summarized in Table 2. The errors for each parameter are found using a value of ${\rm{\Delta }}{\chi }^{2}$ = 2.706, corresponding to a 90% confidence level for one parameter.

Zoom In Zoom Out Reset image size

Table 2.  NuSTAR Spectral Fit Summary, with Integral Flux Values (in Units of $\times {10}^{-11}$ erg cm⁻² s⁻¹) Derived From the Log-parabolic Fits

	Power law		Broken Power law				Log Parabola
Obs.	Index	PL	Index	Index	${E}_{{\rm{break}}}$	BKNPL	Index	Curvature	LP	3–7 keV	7–30 keV
ID	Γ	${\chi }^{2}/$dof	${{\rm{\Gamma }}}_{1}$	${{\rm{\Gamma }}}_{2}$	(keV)	${\chi }^{2}/$dof	Γ	β	${\chi }^{2}/$dof	Flux	Flux
002	2.216 ± 0.009	831/700	2.04 ± 0.03	2.34 ± 0.02	6.3 ± 0.4	747/698	2.290 ± 0.010	0.26 ± 0.03	729/699	3.72 ± 0.02	4.81 ± 0.03
004	2.191 ± 0.006	1204/889	1.25 ± 0.20	2.21 ± 0.01	3.1 ± 0.1	1211/887	2.250 ± 0.008	0.21 ± 0.02	1051/888	5.19 ± 0.02	6.98 ± 0.05
006	2.060 ± 0.006	1246/924	1.92 ± 0.02	2.22 ± 0.02	7.9 ± 0.4	1057/922	2.115 ± 0.008	0.24 ± 0.02	1024/923	12.08 ± 0.09	18.6 ± 0.1
008	2.081 ± 0.007	1152/863	1.90 ± 0.02	2.25 ± 0.02	7.4 ± 0.3	914/861	2.149 ± 0.008	0.32 ± 0.02	892/862	10.75 ± 0.05	16.4 ± 0.1

Notes. Data, models, and ratios, are shown in Figure 1. The indices of the LP fits are derived at 10 keV. The errors for each parameter are found using a value of ${\rm{\Delta }}{\chi }^{2}$ = 2.706, corresponding to a 90% confidence level for a parameter. Observation IDs are shortened by removing the first 60002024 identifier in column one.

Download table as:  ASCIITypeset image

For all four NuSTAR observations, the X-ray emission of Mrk 501 is best represented with a log parabola. A statistical F-test (Snedecor & Cochran 1989) using the ${\chi }^{2}$ and degrees of freedom (dof) of the PL versus LP fit results in F-statistics of 97.8, 129.3, 200.1 and 251.3 for the observations 002, 004, 006, and 008, respectively, corresponding to probabilities of $1.1\times {10}^{-21},$ $4.6\times {10}^{-28},$ $2.9\times {10}^{-41}$, and $7.9\times {10}^{-50}$ for being consistent with the null PL hypothesis. The broken power-law fit to the second NuSTAR observation, ID 004, produces a break energy at the lower limit of the NuSTAR sensitivity window, and is interpreted as a failed fit. The other three observations fit the break energy near ${E}_{{\rm{break}}}$ = 7 keV, motivating the decision to present the NuSTAR flux values in the 3–7 and 7–30 keV bands throughout this work. The upper bound of 30 keV is the typical orbit-timescale detection limit for the Mrk 501 observations.

The NuSTAR observations show the blazar to be in a relatively low state for the first two observations, and a relatively high state during the last two observations, with the 3–7 keV integral fluxes derived from the log-parabolic fits 2–4 times higher than found for the first two observations. More specifically, the average 3–7 keV integral flux values (in units of 10⁻¹¹ erg cm⁻² s⁻¹) were 3.72 ± 0.02 and 5.19 ± 0.02, respectively, for the observations occurring on MJD 56395 and 56420, and 12.08 ± 0.09 and 10.75 ± 0.05, respectively, for the observations starting on MJD 56485 and 56486. In the same flux units, the 7–30 keV integral flux values for the first two observations are similarly 3–4 times lower than the flux states observed in the last two observations (4.81 ± 0.03 and 6.98 ± 0.05 on MJD 56395 and 56420 as compared to 18.6 ± 0.1 and 16.4 ± 0.1 on MJD 56485 and 56486). These integral flux values are summarized in Table 2.

The NuSTAR observations extend across multiple occultations by the Earth, and the integral flux and index (Γ) light curves for the orbits of each extended observation are shown in Figure 2. The periods with simultaneous observations with the ground-based TeV instruments of MAGIC and VERITAS are highlighted by gray and brown bands in the upper portion of each light curve. The observations and results from MAGIC and VERITAS for these time periods are summarized in Section 3.1.

Zoom In Zoom Out Reset image size

The 3–7 and 7–30 keV integral flux values of the first exposure (Observation ID 002) show low variability (${\chi }^{2}=7.0$ and 13.4 for 5 dof), while the trend of increasing flux in both the 3–7 and 7–30 keV bands is clear during the second observation (Observation ID 004). The 7–30 keV flux increases from (5.1 ± 0.1) × 10⁻¹¹ erg cm⁻² s⁻¹ to (8.8 ± 0.1) × 10⁻¹¹ erg cm⁻² s⁻¹ in fewer than 16 hr. The 7–30 keV increases from (1.7 ± 0.1) × 10⁻¹⁰ erg cm⁻² s⁻¹ to (2.0 ± 0.1) × 10⁻¹⁰ erg cm⁻² s⁻¹ in fewer than 7 hr on MJD 56485 (Observation ID 006) and significantly decreases from (1.9 ± 0.1) × 10⁻¹⁰ erg cm⁻² s⁻¹ to (1.4 ± 0.1) × 10⁻¹⁰ erg cm⁻² s⁻¹, again in fewer than 7 hr on MJD 56486 (Observation ID 008).

The relation between the log-parabolic photon indices and 7–30 keV flux values resulting from the fits to the NuSTAR observations of Mrk 501 are shown for each observation separately in Figure 3. The curvature β was not seen to change significantly from orbit to orbit and therefore was fixed at the average value found for each observation (see Table 2 for values). The count rate light curves show no indications of variability on a timescale of less than an orbit period (∼90 minutes). As observed previously in the X-ray band for Mrk 501 (Kataoka et al. 1999), the source was displaying a harder-when-brighter trend during this campaign. This has also been observed in the past for Mrk 421 (Takahashi et al. 1996).

Zoom In Zoom Out Reset image size

3.1. VHE Gamma-rays

3.1.1. MAGIC

MAGIC is a VHE instrument composed of two imaging atmospheric Cherenkov telescopes (IACTs) with mirror diameters of 17 m, located at 2200 m above sea level at the Roque de Los Muchachos Observatory on La Palma, Canary Islands, Spain. The energy threshold of the system is 50 GeV and it reaches an integral sensitivity of 0.66% of the Crab Nebula flux above 220 GeV with a 50-hr observation (Aleksić et al. 2015a).

MAGIC observed Mrk 501 in 2013 from April 9 (MJD 56391) to August 10 (MJD 56514). On July 11 (MJD 56484), ToO observations were triggered by the high count rate of ∼15 counts s⁻¹ observed by Swift XRT (see Section 3.3). The flaring state was observed intensively for five consecutive nights until July 15 (MJD 56488). After that the observations continued with a lower cadence until August 10.

The source was observed during 17 nights, collecting a total of 22 hr of data with zenith angles between 10° and 60°. Only five hours survived the standard quality cuts for regular MAGIC data analysis because many observations were taken during the presence of a Saharan sand–dust layer in the atmosphere known as "Calima." As we explain below, using the Light Detection and Ranging (LIDAR) information we could recover 10 of the 17 hr which would have been rejected otherwise. The telescopes were operated in the so-called wobble mode (Fomin et al. 1994), where the pointing direction is changed every 20 (or 15) minutes among 2 (or 4) positions with an offset of 04 from the source position.

All the data were analyzed following the standard procedure (Aleksić et al. 2012) using the MAGIC Analysis and Reconstruction Software (MARS; Zanin et al. 2013). An image cleaning was applied based on information of signal amplitude and timing of each pixel, and the shower images were parametrized using the Hillas parameters (Hillas 1985). For the reconstruction of the gamma-ray direction and the gamma-hadron separation, the random forest method is applied using the image parameters and the stereoscopic parameters. (Albert et al. 2008; Aleksić et al. 2010). The energy reconstruction utilizes look-up tables. The analysis steps were confirmed independently with data from the Crab Nebula and dedicated Monte Carlo simulations of gamma-ray showers.

A fraction of the data set (10.4 of 15.1 hr, specifically the observations between MJD 56485 and MJD 56514) was affected by "Calima," a Saharan sand–dust layer in the atmosphere. A correction within the framework of the MARS software is applied to account for the absorption due to Calima using LIDAR measurements taken simultaneously with the MAGIC observations (Fruck et al. 2013). The correction was carried out in two steps. Due to the dust attenuation during Calima, the estimated energy is shifted toward low energies, and thus is corrected event by event, as the first step. Then, to account for the shift of the energy estimation, a correction to the collection area is applied as a second step, due to the energy dependence in the collection area. The atmospheric transmission values for this method were obtained from the temporally closest LIDAR measurement. During the observations affected by Calima the atmospheric transmission ranged from 85% down to 60%, being relatively stable within a timescale of one day, which is a typical feature of a Calima layer (unlike a cloudy sky). The precision on the energy correction is estimated to be around 5% of the attenuation (40%–15%), which corresponds to $\lt 2\%$ of the estimated energy, at most. After the Calima correction, the energy threshold increases inversely proportional to the transmission value. This correction method was tested independently on a Crab Nebula data set observed under similarly hazy weather conditions (Fruck & Gaug 2015). Details of the method can be found in Fruck (2015). This is the first time an event-by-event atmospheric correction is applied to MAGIC data.

The analysis results of the MAGIC data taken during good weather conditions have a systematic uncertainty in the flux normalization and in the energy scale. For both of them, the component changing run-by-run is estimated to be ∼11% using Crab Nebula observations (Aleksić et al. 2015a). It is attributed mainly to the atmospheric transmission of the Cherenkov light, which can change on a daily basis (even during so-called good weather conditions) and the mirror reflectivity, which can change also on a daily basis due to the deposition of dust. The atmospheric correction applied in the analysis of the data taken during Calima increases this run-by-run systematic error from 11% to 15% due to the uncertainty in the correction. Since the systematic uncertainty can be different according to the atmospheric correction,we have added 15% or 11% (with or without the atmospheric correction) to the statistical errors of the flux in quadrature for the evaluation of flux variability.

The summary of the MAGIC analysis results for observations occurring simultaneously with NuSTAR is provided in Table 3. The derived spectra are shown in Figure 4, where the spectral points are drawn with statistical errors only. The resultant flux values above 200 GeV range from $(2.39\pm 0.51)\times {10}^{-11}$ ph cm⁻² s⁻¹ (0.11 Crab Nebula flux) on MJD 56395 to $(5.52\pm 0.87)\times {10}^{-10}$ ph cm⁻² s⁻¹ (2.5 times the Crab Nebula flux) on MJD 56484. As seen in the overall light curve (top panel of Figure 5, shown again only with statistical errors), MAGIC observations indicate a significant variability around MJD 56484. A hint of intra-night variability was observed on MJD 56486 and 56487 simultaneously with the NuSTAR observations, as shown in the zoomed-in light curve (top panel of Figure 6). During these two nights the VHE emission is consistent with a constant flux, resulting in a ${\chi }^{2}$/dof of 7.3/4 (12% probability) with the inclusion of the systematic error. Without accounting for the additional systematic error, the constant fit to the flux results in a ${\chi }^{2}$/dof of 57/4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  MAGIC and VERITAS Observations, Analysis, and Spectral Fit Summary for NuSTAR-simultaneous Observations

Exposure	Exposure	Exposure	Instrument	Zenith	Detection	Power-law	Integral Flux	${\chi }^{2}$	dof
Start MJD	Stop MJD	Length		Angle	Significance	Index	> 200 GeV
		(hr)		(deg)	(σ)		${\rm{(}}\times {10}^{-11}$ ph cm−2 s−1)
56395.179	56395.223	1.0	MAGIC	10–14	7.8	2.50 ± 0.24	2.39 ± 0.44	0.58	6
56395.336	56395.493	2.5	VERITAS	15–35	8.3	3.1 ± 0.4	1.85 ± 0.38	0.76	5
56421.142	56421.209	1.1	MAGIC	12–28	12.5	2.24 ± 0.08	5.08 ± 0.54	15.5	13
56421.340	56421.462	1.0	VERITAS	20–32	14.7	2.25 ± 0.15	4.45 ± 0.61	6.9	9
56485.972	56486.014	1.0	MAGIC	12–24	20.4	2.19 ± 0.07	20.8 ± 1.2	10.0	12
56486.039	56486.083	1.0	MAGIC	28–43	20.7	2.39 ± 0.08	25.2 ± 1.3	26.5	10
56486.106	56486.148	1.0	MAGIC	48–60	14.3	2.71 ± 0.12	32.4 ± 2.0	11.9	11
56485.972	56486.148	2.9	MAGIC	12–60	32.3	2.28 ± 0.04	24.3 ± 0.8	24.1	15
56486.966	56487.022	1.3	MAGIC	12–27	25.2	2.37 ± 0.06	24.9 ± 1.1	20.3	12
56487.050	56487.091	0.9	MAGIC	33–46	18.5	2.23 ± 0.09	17.8 ± 1.0	14.5	11
56486.966	56487.091	2.2	MAGIC	12–46	31.8	2.31 ± 0.05	20.9 ± 0.7	30.4	12

Notes. Observations occurring on the same day are grouped with horizontal lines. Daily average values of MAGIC observations are shown in bold, below the results for each observation occurring on that day. Statistical (1σ) error bars are provided for the power-law indices and the integral fluxes. The flux value between MJD 56486.106 and 56486.148 (shown in italics) is estimated with fitting parameters due to an energy threshold above 200 GeV. The significance of the observed gamma-ray signals is computed according to Equation (17) in Li & Ma (1983).

Download table as:  ASCIITypeset image

3.1.2. VERITAS

VERITAS is a VHE instrument comprised of four 12-m IACTs and is sensitive to gamma-rays between ∼100 GeV and ∼30 TeV (Holder et al. 2006; Kieda 2013). This instrument can detect 1% Crab Nebula flux in under 25 hr. VERITAS observed Mrk 501 fourteen times between 2013 April 7 (MJD 56389) and 2013 June 18 (MJD 56461), with 2.5 and 1.0 hr quality-selected exposures occurring simultaneously with NuSTAR on MJD 56395 and MJD 56421, respectively. On days without simultaneous NuSTAR observations, the exposure times ranged between 0.5 and 1.5 hr. The observations occurring simultaneously with NuSTAR are summarized in Table 3. Due to an annual, ∼2 month long monsoon season in southern Arizona where VERITAS is located, no VERITAS observations were possible for this campaign after 2013 June 18.

The VERITAS observations were taken with 05 offset in each of the four cardinal directions to enable simultaneous background estimation (Fomin et al. 1994). Events were reconstructed following the procedure outlined in Acciari et al. (2008a). The recorded shower images were parameterized by their principal moments, giving an efficient suppression of the far more abundant cosmic-ray background. Cuts were applied to the mean scaled width, mean scaled length, apparent altitude of the maximum Cherenkov emission (shower maximum), and θ, the angular distance between the position of Mrk 501 and the reconstructed origin of the event. The results were independently reproduced with two analysis packages (Cogan 2008; Prokoph 2013). The uncertainty on the energy calibration of VERITAS is estimated at 20%. Additionally, the systematic uncertainty on the spectral index is estimated at 0.2, appearing to be relatively independent of the source slope (Madhavan 2013).

A differential power law is fit to the data (${dN}/{dE}\propto {E}^{-{\rm{\Gamma }}}$) to characterize the VHE spectrum of the source. VERITAS observed Mrk 501 to vary by no more than a factor of three in flux throughout the observations, with the integral flux ranging from (1.85 ± 0.38) × 10⁻¹¹ ph cm⁻² s⁻¹ above 200 GeV (8% Crab Nebula flux above the same threshold) on MJD 56395 to (4.45 ± 0.61) × 10⁻¹¹ ph cm⁻² s⁻¹ (20% Crab Nebula flux) on MJD 56421. The source displayed low spectral variability, ranging between ${\rm{\Gamma }}=3.1\pm 0.4$ in the low flux state to ${\rm{\Gamma }}=2.19\pm 0.07$ in the higher flux state. The observation and analysis results are summarized in Table 3 (for NuSTAR simultaneous observations only), with the VHE spectra of the NuSTAR simultaneous observations shown in Figure 4. Day-to-day uncertainties in flux calculations that might be introduced by different atmospheric conditions (even under strictly good weather conditions) are not included in Table 3 and are estimated at less than 10%.

3.1.3. VHE Results

The full light curve of VHE observations from MAGIC and VERITAS is shown in Figure 5, with a zoom into the period of elevated flux in Figure 6. The flux values are shown with statistical errors only. The MAGIC and VERITAS observations of Mrk 501 in 2013 show the source in states which are consistent with the range of states observed in the past. The observations of VERITAS, occurring primarily in the beginning of the campaign, detected the source in a 5%–10% Crab state, in agreement with the early MAGIC observations. Later on in the campaign, MAGIC observed a flux elevated state of order ∼2.5 times the Crab flux.

3.2. HE Gamma-rays

3.3. Swift X-Ray and UV Telescope Observations

The XRT onboard Swift (Gehrels et al. 2004) is a focusing X-ray telescope sensitive to photons with energies between 0.3 and 10 keV. The Swift satellite observed Mrk 501 59 times between 2013 January 1 and September 5 (MJD 56293–56540). All XRT observations were carried out using the Windowed Timing readout mode. The data set was first processed with the XRTDAS software package (v.2.9.0) developed at the ASI Science Data Center and distributed by HEASARC within the HEASoft package (v. 6.13). Event files were calibrated and cleaned with standard filtering criteria with the xrtpipeline task using the calibration files as available in the Swift CALDB version 20140120.

The spectrum from each observation was extracted from the summed and cleaned event file. Events for the spectral analysis were selected within a circle of 20 pixel ($\sim 46^{\prime\prime}$) radius, which encloses about 80% of the Swift XRT point-spread function (PSF), centered on the source position. The background was extracted from a nearby circular region of 40 pixel radius. The ancillary response files were generated with the xrtmkarf task, applying corrections for PSF losses and CCD defects using the cumulative exposure map. The latest response matrices (v.014) available in the Swift CALDB were used. Before the spectral fitting, the 0.3–10 keV source energy spectra were binned to ensure a minimum of 20 counts per bin.

The data were fit with an absorbed power-law model, with index Γ, as well as an absorbed log-parabolic model, where in both cases the neutral hydrogen column density was set at 1.55 × 10²⁰ cm⁻², taken from Kalberla et al. (2005). The summary of the XRT observations and spectral analysis results are provided in Table 4. The light curve of the observations, including 0.3–3 and 3–7 keV integral flux bands, is shown in Figure 5, with a zoom into the period of elevated flux in Figure 6. The 3–7 keV band is not traditionally quoted for Swift XRT data, but is motivated by direct comparison to the 3–7 keV band computed for the NuSTAR observations.

Table 4.  Swift XRT Observations and Analysis Results for NuSTAR-simultaneous Periods

Observation	Date	Exp	Flux	Flux	Flux	Flux	Index	${\chi }^{2}/$dof	Γ	β	${\chi }^{2}/$dof
ID	(MJD)	(s)	2–10 keV	0.5–2 keV	3–7 keV	0.3–3 keV	Γ		LP	LP
00080176001	56395.06	9636.0	6.9 ± 0.1	6.41 ± 0.06	3.6 ± 0.1	11.0 ± 0.1	2.05 ± 0.01	403.5/416	2.06 ± 0.02	−0.02 ± 0.04	402.6/415
00091745001	56485.84	250.7	21.1 ± 1.7	12.7 ± 0.4	10.9 ± 0.9	22.3 ± 0.7	1.77 ± 0.05	108.1/94	1.74 ± 0.08	0.10 ± 0.16	107.0/93
00030793235	56485.98	709.1	24.3 ± 1.1	14.6 ± 0.2	13.1 ± 0.9	24.1 ± 0.4	1.77 ± 0.03	228.7/222	1.75 ± 0.05	0.03 ± 0.09	227.6/221
00030793236	56486.31	1002.0	24.0 ± 0.7	14.1 ± 0.3	13.4 ± 0.6	23.4 ± 0.4	1.73 ± 0.03	291.6/270	1.68 ± 0.04	0.13 ± 0.08	285.1/269
00030793237	56487.04	949.5	19.1 ± 0.9	12.0 ± 0.2	10.4 ± 0.4	18.9 ± 0.3	1.76 ± 0.03	229.9/237	1.73 ± 0.05	0.07 ± 0.08	228.9/236

Notes. Integral flux values are calculated according to the PL model, and are provided in ×10⁻¹¹ erg cm⁻² s⁻¹ units. The errors for each parameter are found using a value of ${\rm{\Delta }}{\chi }^{2}$ = 2.706, corresponding to a 90% confidence level for a parameter.

Download table as:  ASCIITypeset image

Mrk 501 displays a relatively steady flux state until after MJD 56480, when the flux increases to (38.3 ± 1.5) × 10⁻¹¹ erg cm⁻² s⁻¹ on MJD 56483 (corresponding to the day with the XRT count rate of 15 counts s⁻¹ which triggered MAGIC and NuSTAR observations). This high X-ray state was followed by a general drop in flux, continuing through the last XRT observation included in this work (2013 September 1; MJD 56540).

The power-law fitted indices and 3–7 keV flux derived from the power-law fits are plotted in Figure 7 for all 59 observations. The source clearly displays the harder-when-brighter trend found previously in other TeV blazars, such as Mrk 421 (Takahashi et al. 1996). This behavior is similar to that displayed in the hard X-ray band 7–30 keV observed by NuSTAR and shown in Figure 3. Notably, the photon indices in the soft X-ray band are systematically harder than those observed by NuSTAR in the 7–30 keV band. The spectral index observed by Swift XRT (Γ, determined at 1 keV) ranges between 1.4 and 2.2 (Figure 7) while the NuSTAR index, determined at 10 keV, ranges from 2.1 to 2.4 (Figure 3).

Zoom In Zoom Out Reset image size

Additionally, UV/optical observations were collected with the UVOT onboard Swift. These observations were carried out using the "filter of the day," i.e., one of the six lenticular filters (V, B, U, UVW1, UVM2, and UVW2), unless otherwise specified in the ToO request, so images are not always available for all filters. There are 50 observations included in this Mrk 501 campaign, 18 of which included exposures in all filters while the remaining 32 observations contain UV imaging only.

For each filter observation, we performed aperture photometry analysis using the standard UVOT software distributed within the HEAsoft 6.10.0 package and the calibration included in the latest release of CALDB. Counts were extracted from apertures of 5'' radius for all filters and converted to fluxes using the standard zero points from Poole et al. (2008). The flux values were then de-reddened using the value of $E(B-V)$ = 0.017 (Schafly & Finkbeiner 2011) with ${A}_{\lambda }/E(B-V)$ ratios calculated for UVOT filters using the mean Galactic interstellar extinction curve from Fitzpatrick (1999). No variability was detected to occur within single exposures in any filter. The processing results were verified, checking for possible contamination from nearby objects falling within the background apertures.

3.4. Optical

3.5. Radio

Since Mrk 501 is highly variable, detailed inferences regarding the broadband SED and its temporal evolution require simultaneous observations of multiple bands. In particular, for the determination of the low-energy peak ${E}_{{\rm{syn}}},$ and the flux at ${E}_{{\rm{syn}}},$ $F({E}_{{\rm{syn}}}),$ Swift XRT and NuSTAR observations must be simultaneous. There are five periods within the campaign for Mrk 501 where the observations by NuSTAR and Swift occurred within one hour of each other. The Swift exposure IDs for these quasi-simultaneous periods are summarized in Table 4. For Mrk 501, ${E}_{{\rm{syn}}}$ is located in the X-ray band and can be determined reliably (except for the first NuSTAR observation where ${E}_{{\rm{syn}}}$ is $\leqslant 0.85$ keV) since there is no evidence of X-ray variability of Mrk 501 on a timescale shorter than a NuSTAR orbit (∼90 minutes).

As a precursor to the joint fitting of XRT and NuSTAR data, we confirm agreement between the 3–7 keV flux values derived from the Swift XRT and NuSTAR fitted models. There is a residual discrepancy (not a uniform offset) at the level of $\lt 10\%.$ Using XSPEC, we performed simultaneous fitting to the data sets using the absorbed log-parabolic model as done in Section 2 for the NuSTAR data alone. During the fitting process, we allowed the normalizations of the data sets to vary, but required the same spectral shape parameters. A representative plot of the simultaneous fit for XRT and NuSTAR data collected on MJD 56485 is provided in Figure 8. The model spectrum is shown as a solid line in Figure 8. The agreement between XRT and NuSTAR was studied and found to be within the calibration uncertainties.⁹⁹

Zoom In Zoom Out Reset image size

For the determination of the spectral parameters characterizing the synchrotron peak (namely the energy ${E}_{{\rm{syn}}}$ and $F({E}_{{\rm{syn}}}$)) with the simultaneous NuSTAR and Swift XRT observations, we apply the log-parabolic model modified by the photoelectric absorption due to our Galaxy, with a (fixed) neutral hydrogen column density of $1.55\times {10}^{20}$ cm⁻², taken from Kalberla et al. (2005). The procedure to search for ${E}_{{\rm{syn}}}$ involves the variation of the "normalization energy" parameter (in the ${\mathtt{logpar}}$ model in XSPEC) until the local index Γ returns a value of 2—then ${E}_{{\rm{syn}}}$ corresponds to the peak in the $E\times F(E)$ representation. This procedure correctly accounts for the effect of the soft X-ray absorption by Galactic column density as the absorption is included in the model fitted to the data. For the determination of the error on ${E}_{{\rm{syn}}},$ we freeze the "local index"—defined at energy ${E}_{{\rm{syn}}}$—to a value of 2, and then step the value of ${E}_{{\rm{syn}}}$ keeping all other parameters free. We then search for the value of the ${E}_{{\rm{syn}}}^{\prime }$ which corresponds to the departure of ${\chi }^{2}$ from the minimum by ${\rm{\Delta }}{\chi }^{2}=2.7.$ The error quoted is the difference between ${E}_{{\rm{syn}}}$ and ${E}_{{\rm{syn}}}^{\prime }.$ The ${E}_{{\rm{syn}}}$ and curvature parameters (β) for each of the simultaneous data sets are summarized in Table 5. We quote the value of $F({E}_{{\rm{syn}}})$ inferred from the NuSTAR module FPMA (Focal Plane Module A).

Table 5.  Fitting Results for Swift XRT and NuSTAR Simultaneous Observations

Observation	Date	Orbit	${E}_{{\rm{syn}}}$	$F({E}_{{\rm{syn}}})$	Curvature	${\chi }^{2}/$dof
ID	(MJD)	Number	(keV)	($\times {10}^{-11}$ erg cm−2 s−1)	β
60002024002	56395.1	1	<0.85	4.1	0.061	669/673
60002024006	56485.9	1	4.9 ± 0.7	13.8	0.21	596/577
60002024006	56486.0	2	5.1 ± 0.9	13.7	0.22	697/715
60002024006	56486.2	4	7.0 ± 0.8	14.6	0.2	877/848
60002024008	56487.1	4	3.3 ± 0.9	11.2	0.17	832/851

Note. The data were simultaneously fit with a log-parabolic function.

Download table as:  ASCIITypeset image

The combination of Swift XRT and NuSTAR observations provides an unprecedented view of the synchrotron peak variability. From Table 5, it is evident that the synchrotron peak moves by a factor of about ten during this campaign, with the highest synchrotron peak occurring during the elevated X-ray and gamma-ray state.

5.1. Fractional Variability

In order to quantify the broadband variations we utilize the fractional variability, ${F}_{{\rm{var}}}.$ We follow the description given in Vaughan et al. (2003), where ${F}_{{\rm{var}}}$ is calculated as

$$\begin{eqnarray}&&{F}_{{\rm{var}}}=\sqrt{{\displaystyle \frac{{S}^{2}-\langle {\sigma }^{2}\rangle }{\langle {F}_{\gamma }\rangle }}^{2}}\end{eqnarray} \tag{ 4 }$$

where $\langle {F}_{\gamma }\rangle$ is the average photon flux, S is the standard deviation of the flux measurements, and $\langle {\sigma }^{2}\rangle$ is the mean squared error of the measurement.

${F}_{{\rm{var}}}$ was determined for the temporal binning and sampling presented in Figure 5 and Table 3 (for MJD 56485 and 56486, the bold lines in Table 3 are used). The value of ${F}_{{\rm{var}}}$ is known to be dependent on sampling and should be interpreted with caution. For example, a well sampled light curve with small temporal bins will allow us to probe the variability on small timescales (e.g., NuSTAR), which could be hidden if the variability is computed with fluxes obtained with relatively coarse temporal bins (e.g., Fermi LAT).

The fractional variability for each band (from 15 GHz radio through VHE) is shown in Figure 9. For the period of observations covered in this work, the fractional variability shows a double-peaked shape with the highest variability in the X-ray and VHE bands. A similar broadband variability pattern has recently been reported for Mrk 501 (Doert 2013; Aleksić et al. 2015c), for Mrk 421 (Aleksić et al. 2015b; M. Baloković et al. 2015, in preparation), and for other high-synchrotron-peaked blazars in, for example, Aleksić et al. (2014). This double-peaked shape of ${F}_{{\rm{var}}}$ from radio through VHE can be interpreted as resulting from a correlation between the synchrotron and inverse-Compton peaks.

Zoom In Zoom Out Reset image size

${F}_{{\rm{var}}}$ is below ∼5% at 15 GHz and optical/UV frequencies, while at 37 GHz the fractional variability is ∼20%. The relatively high fractional variability at 37 GHz is not produced by any single flaring event, but rather by a consistent flickering in the radio flux. Such flickering is not typically observed in blazars, but has been reported for Mrk 501 in Aleksić et al. (2015c). At X-ray frequencies, ${F}_{{\rm{var}}}$ gradually increases with energy, reaching the largest value (∼0.6) in the 7–30 keV band measured by NuSTAR. The ${F}_{{\rm{var}}}$ computed for the Swift XRT 3–7 keV observations is higher than for the NuSTAR 3–7 keV fluxes due to the larger temporal coverage of the Swift observations, allowing for observation of Mrk 501 during high activity levels that were not observed with NuSTAR.

The Swift XRT ${F}_{{\rm{var}}}$ for Mrk 501 published in Stroh & Falcone (2013) was 0.15 or 0.18, depending on the timescale used for calculation, illustrating that the value of ${F}_{{\rm{var}}}$ is dependent on sampling. In Abdo et al. (2011a), RXTE-ASM (2–10 keV) and Swift BAT (15–50 keV) show ${F}_{{\rm{var}}}$ values between 0.2 and 0.3, although it should be noted that due to the limited sensitivity of RXTE-ASM and Swift BAT (in comparison with Swift XRT and NuSTAR), the variability was studied on timescales larger than 30 days.

5.2. Cross Correlations

5.3. X-Ray/VHE Correlation

The light curve of the broadband observations is shown in Figure 5, with a zoom of the period showing an elevated X-ray and VHE state in Figure 6. The VERITAS and MAGIC flux points within the light curve are shown with statistical errors only. Correlation studies using the VHE flux values are completed with statistical and systematic errors included, as described below. The radio, optical, and UV observations show relatively steady flux over the campaign period, while the largest amplitude of variability can be seen in the X-ray and VHE gamma-ray bands. An elevated state in both the X-ray and VHE bands can be seen to occur on MJD 56483 (Swift Observation ID 00030793232 in Table 4). Zooming in on this epoch (Figure 6), shows that the NuSTAR observations occurring on MJD 56485 and 56486 occurred after the highest state observed by MAGIC and Swift. The XRT observations show an elevated X-ray flux in both the 0.3–3 and 3–7 keV bands on MJD 56483.

A comparison between the NuSTAR-observed X-ray photon flux values (derived from XSPEC) in the 3–7 and 7–30 keV bands and the epochs of simultaneous VHE observations is shown in Figure 10. During this campaign, 10 observations occurred within one hour between either NuSTAR and MAGIC (seven observations) or NuSTAR and VERITAS (three observations). The simultaneous X-ray and VHE data, where the VHE data include both statistical and systematic errors, were fit with both a linear and a quadratic function.

Zoom In Zoom Out Reset image size

Within the one-zone SSC emission paradigm, there is a physical motivation for a quadratic relationship between the X-ray and VHE flux values (Marscher & Gear 1985). More specifically, the inverse-Compton flux depends not only on the density of photons, but also on the density of the electron population producing those photons. If, however, the particle population is energetic enough for the inverse-Compton scattering to occur in the Klein–Nishina regime, the relationship between the X-ray and VHE fluxes can be complex and will depend in detail on the energy bands considered, the particle energy loss mechanisms and the magnetic field evolution. In particular, Katarzyński et al. (2005) suggest that a roughly linear relationship may arise during the declining part of a flare when the emitting region expands adiabatically, leading to a decrease of both the particle number density and the magnetic field strength.

A quadratic relationship provides a better fit than the linear fit for the 3–7 keV flux values measured simultaneously by NuSTAR, with ${\chi }^{2}$ of 11.4 and 87.3, respectively, for 9 dof. The 3–7 keV flux and the $\gt 200\;{\rm{GeV}}$ flux are highly correlated, with a Pearson correlation coefficient (r) of 0.974. Similarly, for the 7–30 keV band, the quadratic relation fits the data better than the linear relation, with ${\chi }^{2}$ of 17.5 and 79.1, respectively, for 9 dof. The r-value for the 7–30 keV flux and the $\gt 200\;{\rm{GeV}}$ flux is 0.979.

A comparison between the Swift-observed X-ray photon flux values (derived from XSPEC) in the 0.3–3 and 3–7 keV bands and the epochs of simultaneous VHE observations is shown in Figure 11. These data are not simultaneous with the NuSTAR observations shown in Figure 10 and therefore the results cannot be directly compared. During this campaign, 12 absolutely simultaneous observations occurred between Swift and MAGIC (10) and Swift and VERITAS (2), shown in Figure 11. Similarly as done for the NuSTAR bands, the simultaneous Swift X-ray and VHE data were fit with both a linear and a quadratic function with an offset fixed to zero. For the 0.3–3 keV flux values measured simultaneously by Swift, a quadratic relationship provides a better fit than the linear fit, with ${\chi }^{2}$ of 81.8 and 162.0, respectively, for 11 dof. The 0.3–3 keV flux and the $\gt 200\;{\rm{GeV}}$ flux are highly correlated, with a Pearson correlation coefficient (r) of 0.958. For the 3–7 keV band, the quadratic function fits the data better than the linear function, with ${\chi }^{2}$ of 58.0 and 114.0, respectively, for 11 dof. The r-value for the 3–7 keV flux as measured with Swift and the $\gt 200\;{\rm{GeV}}$ flux is 0.954.

Zoom In Zoom Out Reset image size

Previous MWL campaigns on Mrk 501 have been sufficiently characterized with a one-zone SSC model (Abdo et al. 2011a; Acciari et al. 2011), although there are a few notable instances where a one-zone SSC model was found not to be appropriate for the broadband emission (Pian et al. 1998; Kataoka et al. 1999). In this study we decided to use the simplest approach, which is provided by a leptonic model with a single emitting region. The broadband spectral data were modeled with an equilibrium version of the single-zone SSC model from Böttcher & Chiang (2002) and Böttcher et al. (2013). This model has been used to describe the broadband emission from various other VHE-detected blazars (e.g., Acciari et al. 2009a; Abdo et al. 2011b; Aliu et al. 2011, 2013).

Within this equilibrium model, the emission originates from a spherical region of relativistic leptons with radius R. This emission region moves down the jet with a Lorentz factor Γ. We set the Doppler factor δ to 15 for all model representations. Notably, it has been shown that when using least-squares fitting of emission models to broadband data of Mrk 501, the Doppler factor can vary widely from state to state (Mankuzhiyil et al. 2012). We do not complete least-squares fitting in this work and instead choose to fix the Doppler factor to 15 for the representation of all states, limiting the number of free parameters of the SSC model. The Doppler factor of 15 is similar to the Doppler factor used in previous studies of Mrk 501 (Abdo et al. 2011a; Acciari et al. 2011; Mankuzhiyil et al. 2012). In order to reduce the number of free parameters, the jet axis is aligned toward the line of sight with the critical angle $\theta =3\buildrel{\circ}\over{.} 8.$ At the critical angle, the jet Lorentz factor is equal to the Doppler factor (${\rm{\Gamma }}=\delta$).

Within this emission model, relativistic leptons are injected into this emission region continuously according to a power-law distribution $Q={Q}_{0}{\gamma }^{-q}$ between ${\gamma }_{{\rm{min}}}$ and ${\gamma }_{{\rm{max}}}.$ The injected population of particles is allowed to cool. The simulation accounts for synchrotron emission due to a tangled magnetic field B₀, Compton up-scattering of synchrotron photons, $\gamma \;\gamma$ absorption and the corresponding pair production rates (via the general solution in Böttcher & Schlickeiser 1997). The cooling of the injected electrons is dominated by radiative losses, which are balanced by injection and particle escape from the system. This particle escape is characterized with an escape efficiency factor η = 100, where ${t}_{{\rm{esc}}}=\eta$R/c. The use of η = 100 is motivated by success in representing SEDs of TeV blazars in previous studies using the same model (e.g., Aliu et al. 2013). The electron cooling rates and photon emissivity and opacity are calculated using similar routines of the code for jet radiation transfer described in Böttcher et al. (1997). Together, the particle injection, cooling and escape mechanisms lead to an equilibrium particle population.

A key result of the equilibrium that occurs between continual particle injection, particle escape, and radiative cooling is a break in the electron distribution ${\gamma }_{b}$ (referred to as ${\gamma }_{c}$ within Böttcher et al. 2013), where ${t}_{\mathrm{esc}}={t}_{\mathrm{cool}}({\gamma }_{b}).$ As described in Equations (1) and (2) of Böttcher et al. (2013), if ${\gamma }_{b}$ is smaller than ${\gamma }_{\mathrm{min}},$ the system will be in a fast cooling regime. If ${\gamma }_{b}$ is greater than ${\gamma }_{\mathrm{min}},$ the system will be in a slow cooling regime. Within the fast cooling regime, the equilibrium particle distribution is a broken power law, with an index of 2 for particles with Lorentz factors less than ${\gamma }_{\mathrm{min}},$ and an index of (q + 1) for Lorentz factors above ${\gamma }_{\mathrm{min}}.$ In the slow cooling regime, the resulting broken power law of the equilibrium particle distribution is equal to the injected spectrum (q) for particles with Lorentz factor below ${\gamma }_{b},$ and $(q+1)$ above ${\gamma }_{b}.$ It is known that a hard injected electron spectrum would lead to a small amount of pile-up, followed by a smooth cut-off toward the HE end of the distribution (for details, see, e.g., Kardeshev 1962 and Stawarz et al. 2008). More specifically, the equilibrium electron spectrum slightly deviates from the ($q+1$) approximation at the HE end ($\gamma \sim {\gamma }_{{\rm{max}}}$) due to pile-up effects that increase as the injected spectrum becomes harder (i.e., $q\lt 1.5$). Notably, although scattering in the KN regime is appropriately accounted for within the SSC model, neither the pile-up at the highest energy nor the energy loss (Compton cooling) of the electrons participating in scattering within the KN regime is accounted for within the model. The two aforementioned effects, however, are expected to result in a negligible deviation of the equilibrium electron spectrum from the approximated index of $q+1$.

L_e is the kinetic power in the relativistic electrons and L_B is the power in the Poynting flux carried by the magnetic field of the equilibrium particle distribution. The L_e and L_B parameters allow the calculation of the equipartition parameter ${L}_{B}/{L}_{e}.$ A state with an equipartition near unity minimizes the total (magnetic field + particle) energy requirement to produce a given synchrotron flux. Therefore, from an energetics point of view a situation near equipartition is usually favored. If the jet is powered by a Blandford–Znajek type mechanism, it is expected to be initially Poynting-flux dominated, and this luminosity is then (through an unknown mechanism, possibly magnetic reconnection) converted partially into particle energy. This conversion is expected to stop at an approximately equipartition situation as an equilibrium is reached between the conversion of magnetic energy to particle energy, and vice versa (via turbulent charged-particle motion generating small-scale, turbulent magnetic fields). For examples of blazar modeling based on equipartition, see Cerruti et al. (2013), Dermer et al. (2014). Alternatively, a sub-equipartition magnetic field may be expected in an MHD-driven, initially particle-dominated jet, where magnetic fields could be self-generated (amplified) by, e.g., shocks. The sub-equipartition magnetic fields that are often found in blazar SED modeling might therefore favor this latter scenario. Sub-equipartition states are a common result in the application of single-zone SSC emission scenarios to VHE blazars, as in Aliu et al. (2011, 2012a, 2012b), Acciari et al. (2008b, 2009a, 2009b, 2009c) and Abdo et al. (2011c).

The broadband data and model representations for five days from the MWL observation campaign are shown in Figure 12. The flux resulting from the model simulation (solid line) is corrected for absorption by interaction with the extragalactic background light (EBL) for the redshift of $z=0.03,$ assuming the EBL model outlined in Domínguez et al. (2011). The model thus represents the observed VHE emission as opposed to the intrinsic VHE emission. When applying the model to the data, the radio flux is likely to include a significant portion of extended radio emission and is therefore taken as an upper limit, as done in Abdo et al. (2011a).

Zoom In Zoom Out Reset image size

The parameters used to represent the data with the equilibrium model are summarized in Table 6. The data in this work are represented with the emission model within the fast cooling regime, where the emitting equilibrium particle population follows $n(e)\propto {\gamma }^{-2}$ for ${\gamma }_{b}\lt \gamma \lt {\gamma }_{\mathrm{min}}$ and $n(e)\propto {\gamma }^{-(q+1)}$ for ${\gamma }_{\mathrm{min}}\lt \gamma \lt {\gamma }_{\mathrm{max}}.$ A particle population with an injection index of $q=1.8-1.9$ provides a reasonable representation of the synchrotron emission on MJD 56395 (red; top panel) and 56420 (green; second panel from top). There are no Swift data for observations on MJD 56420. Each of these epochs (MJD 56395 and 56420) can be sufficiently described with similar parameters, although the SED on MJD 56420 requires a slightly more energetic electron population and lower magnetic field to account for the marginally elevated X-ray and VHE emission as compared to what is observed on MJD 56395.

Table 6.  Single-zone SSC Model Parameter Values (See Section 6 for Overview of Model and Parameters)

Parameter	MJD 56395	MJD 56420	MJD 56485.0	MJD 56485.9	MJD 56486.9
${\gamma }_{\mathrm{min}}$ (×104)	2.0	2.1	2.0	2.0	2.0
${\gamma }_{\mathrm{max}}$ (×106)	1.0	1.4	1.4	1.7	1.4
${\gamma }_{\mathrm{break}}$ (×103)	4.1	4.6	2.8	3.3	3.4
q	1.9	1.8	1.3	1.3	1.3
η	100	100	100	100	100
δ	15	15	15	15	15
B0 (G)	0.06	0.05	0.03	0.03	0.03
Γ	15	15	15	15	15
R (×1015 cm)	7.0	7.0	5.0	7.0	7.0
θ (degrees)	3.8	3.8	3.8	3.8	3.8
tvar (hr)	4.3	4.3	3.1	4.3	4.3
Le (×1042 erg s−1)	9	12	36	28	26
= ${L}_{B}/{L}_{e}$	1.8 × 10−2	6.1 × 10−2	5.3 × 10−4	1.3 × 10−3	1.4 × 10−3

Note. Model representations are shown along with the broadband data in Figure 12.

Download table as:  ASCIITypeset image

Although the highest VHE gamma-ray ($\geqslant 200$ GeV) flux during this campaign was observed by MAGIC on MJD 56484, a reliable spectrum from that MAGIC observation could not be reconstructed due to the presence of Calima and the lack of LIDAR data to correct for it. Swift XRT also recorded the highest X-ray flux in its observation on the same day. On the other hand, there are sufficient broadband data to model the SED on MJD 56485.0 (turquoise; middle panel Figure 12), which is less than one day later than the MAGIC and Swift observation of the highest fluxes occurred.

The light curve in Figure 5 shows that Mrk 501 displayed relatively steady emission in each band between MJD 56420 and the elevated state observed by Swift and MAGIC on MJD 56484. In moving from the relatively quiescent SED on MJD 56420 to the elevated state observed on MJD 56485, a hardening of the injection spectrum is required (q = 1.3) to match the X-ray spectrum observed by Swift XRT. With the injection index responsible for the hardness of the synchrotron emission at X-ray energies, the frequency at which the synchrotron emission peaks, is related to the spectrum of the injected particle population, and the magnetic field (B₀). When moving from the state on MJD 56420 to 56485.0, the strength of the magnetic field decreases, moving the peak of the synchrotron emission to lower energies. The decrease of the synchrotron flux resulting from a lower magnetic field is counteracted with an increase of particle luminosity L_e. Finally, to match the relative magnitudes of the synchrotron and inverse-Compton peak fluxes, the electron and photon density of the emission region was increased with a decrease of the emission region size. The decrease of the emission region size to 5.0 × 10¹⁵ cm on MJD 56485.0 provides a higher inverse-Compton flux while maintaining the synchrotron flux.

Following Blumenthal & Gould (1970), the regime at which the up-scattering is occurring can be estimated (in the observer frame) according to $4\;h{\nu }_{\mathrm{syn}\;\mathrm{pk}}\gamma /\delta {m}_{e}\;{c}^{2},$ where γ represents the energy of the electrons up-scattering ${\nu }_{\mathrm{syn}\;\mathrm{pk}}$ photons. If this quantity is less than 1, the inverse-Compton emission is occurring within the Thomson regime, while if it is greater than 1, the emission is occurring in the KN limit. With ${\nu }_{\mathrm{syn}\;\mathrm{pk}}$ at approximately 5 keV, $4\;h{\nu }_{\mathrm{syn}\;\mathrm{pk}}{\gamma }_{{\rm{min}}}/\delta {m}_{e}\;{c}^{2}\sim 25,$ indicating that, according to the model applied within this work, the inverse-Compton scattering of the photons near the synchrotron peak is far into the KN regime. We note that this is not necessarily in conflict with the quadratic relationship between the simultaneous X-ray and VHE flux measurements, but it implies a reasonably steady value of magnetic field that is supported by our SSC models; see Table 6. For a more extensive discussion, see Katarzyński et al. (2005).

The SEDs on the days MJD 56485.9 and 56486.9 are similar to MJD 56485.0. All model representations explored here result in emission scenarios which are heavily matter dominated (far below equipartition), where the majority of the energy is distributed within the particle population instead of in the magnetic field. Notably, even a single-zone SSC model is difficult to constrain, and the solutions presented here are not applied with the intent of constraining parameter space, but instead to just show that a reasonable representation of the data is possible. There are additional models (e.g., multi-zone or hadronic models) which might alternatively be used to describe the broadband emission from Mrk 501 during these epochs (e.g., Tavecchio et al. 2011; Aleksić et al. 2015b). However, these models have twice as many free parameters as single-zone leptonic models and, in this particular case, there are not strong constraints from MWL flux evolution correlations that point to the necessity of such models.

The inclusion of the hard X-ray telescope NuSTAR in this observational campaign has provided unprecedented insight into the temporal evolution of the 3–30 keV X-rays emitted by Mrk 501. Before this campaign, Mrk 501 had not been observed to display hard X-ray variability on timescales of ∼7 hr. The fractional variability of Mrk 501 observed during this campaign was highly significant for the NuSTAR 7–30 keV band (${F}_{{\rm{var}}}\;=\;0.6$).

Investigation of the DCF allows insight into possible leads or lags between the low (0.3–3 keV) and high (3–7 keV) X-ray and VHE emission. The variability between these two bands shows evidence for a zero day lag. Correlation between the X-ray and VHE bands is further supported by the correlated variability inferred from the Pearson coefficients of 0.958 and 0.954 for simultaneous observations (occurring within one hour), respectively. Correlation is also found between the NuSTAR X-ray flux values and the simultaneous >200 GeV flux values (with observations occurring within one hour), with Pearson coefficients of 0.974 and 0.979 for the 3–7 keV and 7–30 keV bands, respectively.

Correlation of variability between the X-ray and VHE flux, and more notably direct correlation without any lead or lag time, is a natural signature of SSC emission. Within the single-zone SSC paradigm, the inverse-Compton flux is emerging from the same region as the synchrotron emission, and is fundamentally derived from the same particle and photon populations as the synchrotron emission. In this way, any variability in the synchrotron photon luminosity will immediately be translated into a change in the up-scattered inverse Compton luminosity.

In applying a single-zone equilibrium SSC model to the broadband data of Mrk 501, we find that the data could be reasonably represented in each of the five simultaneous epochs. Notably, the injected particle populations on MJD 56485.0, 56485.9, and 56486.9 are very hard, with an injection index of q = 1.3. Such a hard injection index is difficult to produce with standard shock acceleration scenarios alone, but is possible through a magnetic reconnection event (e.g., as explained in Romanova & Lovelace 1992; Guo et al. 2014; Sironi & Spitkovsky 2014). The increase in energy of the particle population (with an additional hardening to the injection index of q = 1.3) between the SED derived for MJD 56485.0 as compared to MJD 56420 indicates an introduction of additional energetic particles to the emission region, requiring some source of energy input. The decrease of the magnetic field, similar to what would naturally occur after a magnetic reconnection event, is capable of accelerating particles near the point of reconnection and producing the newly injected q = 1.3 particle population. Additionally, the decrease in the emission region size is consistent with a magnetic reconnection event that affects a more localized region as compared to a larger, more steady non-thermal emission region. More information on particle acceleration via magnetic reconnection can be found in Werner et al. (2014) and Guo et al. (2015).

The variability timescale for these model representations, quoted in Table 6, is determined from the light-crossing timescale of the emission region according to ${t}_{{\rm{var}}}=R/c\delta (1+z).$ For the emission region sizes and Doppler factor of δ = 15 used within the model, the predicted variability timescales of a couple of hours are compatible with the variability timescale observed during the broadband observations. The radiative cooling timescale is approximately equal to the synchrotron cooling timescale, ${t}_{{\rm{sync}}}\sim 1.4\times {10}^{4}{({B}_{0}/0.06\;G)}^{-2}{\gamma }_{6}^{-1}$ s, where ${\gamma }_{6}=\gamma /({10}^{6}).$ With a minimum light crossing time, corresponding to the minimum variability timescale of ${t}_{{\rm{var}}}\sim 1.6\times {10}^{4}\;{\rm{s}}$ (in the observer frame), all but the most energetic electrons within the emitting region cool on timescales that are longer than the crossing timescale, showing that the observed variability is likely a reflection of changes in the particle acceleration and/or injection processes directly.

Notably, faster variability timescales have been observed from Mrk 501 in the past (e.g., Albert et al. 2007) and so the model parameters shown here cannot be generalized to all Mrk 501 flux variability episodes. NuSTAR observations show the hard X-ray flux to significantly decrease by more than 10% between its 90-minute orbits. Moreover, on MJD 56420 the source hard X-ray flux was observed to change by a factor of greater than 40% in the 7–30 keV band during a 7 hr exposure.

In an attempt to describe a possible emission scenario which might result in the broadband SED variability observed for Mrk 501 in 2013, the parameter changes were made to the single zone equilibrium SSC model monotonically. With a degeneracy between several of the input parameters, the model applied here cannot be used for conclusive studies regarding which changes occur within the emitting region from one state to the next. Instead, through the study of band-to-band spectral variability, leads and/or lags and fractional variability, as well as broadband modeling of various flaring episodes, we find compelling evidence to support a single zone SSC emission scenario for Mrk 501 during the broadband observations in this campaign.

The collection of simultaneous broadband observations is a necessity for the study of the relativistic emission mechanisms at work within blazars such as Mrk 501. It is known that these sources vary continually, with characteristics that significantly change between different flaring episodes, requiring the continuation of deep broadband observations such as those presented in this work.

This work was supported under NASA Contract No. NNG08FD60C, and made use of data from the NuSTAR mission, a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory, and funded by the National Aeronautics and Space Administration. We thank the NuSTAR Operations, Software and Calibration teams for support with the execution and analysis of these observations. This research has made use of the NuSTAR Data Analysis Software (NuSTARDAS) jointly developed by the ASI Science Data Center (ASDC, Italy) and the California Institute of Technology (USA).

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, and the Japanese JSPS and MEXT is gratefully acknowledged. This work was also supported by the Centro de Excelencia Severo Ochoa SEV-2012-0234, CPAN CSD2007-00042, and MultiDark CSD2009-00064 projects of the Spanish Consolider-Ingenio 2010 programme, by grant 268740 of the Academy of Finland, by the Croatian Science Foundation (HrZZ) Project 09/176 and the University of Rijeka Project 13.12.1.3.02, by the DFG Collaborative Research Centers SFB823/C4 and SFB876/C3, and by the Polish MNiSzW grant 745/N-HESS-MAGIC/2010/0.

This research is also supported by grants from the U.S. Department of Energy Office of Science, the U.S. National Science Foundation, and the Smithsonian Institution, by NSERC in Canada, by Science Foundation Ireland (SFI 10/RFP/AST2748), and by STFC in the U.K. We acknowledge 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 VERITAS instrument. The VERITAS Collaboration is grateful to Trevor Weekes for his seminal contributions and leadership in the field of VHE gamma-ray astrophysics, which made this study possible.

The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden.

Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France.

We thank the Swift team duty scientists and science planners and we acknowledge the use of public data from the Swift data archive. This research has made use of the Swift XRT Data Analysis Software (XRTDAS) developed under the responsibility of the ASI Science Data Center (ASDC), Italy.

This research is partly based on observations with the 100-m telescope of the MPIfR (Max-Planck-Institut für Radioastronomie) at Effelsberg and with the IRAM 30-m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). V. Karamanavis and I. Myserlis are funded by the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne.

The OVRO 40-m monitoring program is supported in part by NASA grants NNX08AW31G and NNX11A043G, and NSF grants AST-0808050 and AST-1109911.

T.G. acknowledges support from Istanbul University (Project numbers 49429 and 48285), Bilim Akademisi (BAGEP program) and TUBITAK (project numbers 13AT100-431, 13AT100-466, and 13AT60-430).

St. Petersburg University team acknowledges support from Russian RFBR grant 15-02-00949 and St. Petersburg University research grant 6.38.335.2015.

The research at Boston University (BU) was funded in part by NASA Fermi Guest Investigator grant NNX14AQ58G and Swift Guest Investigator grant NNX14AI96G. The PRISM camera at Lowell Observatory was developed by K. Janes et al. at BU and Lowell Observatory, with funding from the NSF, BU, and Lowell Observatory.

This research was partially supported by Scientific Research Fund of the Bulgarian Ministry of Education and Sciences under grant DO 02-137 (BIn-13/09). The Skinakas Observatory is a collaborative project of the University of Crete, the Foundation for Research and Technology—Hellas, and the Max-Planck-Institut für Extraterrestrische Physik.

The Abastumani Observatory team acknowledges financial support by the by Shota Rustaveli National Science Foundation under contract FR/577/6-320/13.

We are grateful to the American Association of Variable Star Observers (AAVSO) for assisting us with some of the optical data acquisition.

G. Damljanovic and O. Vince gratefully acknowledge the observing grant support from the Institute of Astronomy and Rozhen National Astronomical Observatory, Bulgarian Academy of Sciences. This work is in accordance with the Projects No 176011 ("Dynamics and kinematics of celestial bodies and systems"), No 176004 ("Stellar physics") and No 176021 ("Visible and invisible matter in nearby galaxies: theory and observations") supported by the Ministry of Education, Science and Technological Development of the Republic of Serbia.

Facilities: NuSTAR - The NuSTAR (Nuclear Spectroscopic Telescope Array) mission, MAGIC - Major Atmospheric Gamma Imaging Cherenkov Telescope, VERITAS - Very Energetic Radiation Imaging Telescope Array System, Fermi LAT - , Swift XRT/UVOT - , AAVSOnet - , FLWO - , WEBT, CrAO: 1.25m, Perkins - Lowell Observatory's 72in Perkins Telescope, Rozhen - , Skinakas - , OVRO - , Metsahovi - .