INSIGHTS INTO THE HIGH-ENERGY γ-RAY EMISSION OF MARKARIAN 501 FROM EXTENSIVE MULTIFREQUENCY OBSERVATIONS IN THE FERMI ERA

The list of all the instruments that participated in the campaign is given in Table 1, and the scheduled observations can be found online.¹⁴⁹ In some cases, the planned observations could not be performed due to bad observing conditions, while in some other occasions the observations were performed but the data could not be properly analyzed due to technical problems or rapidly changing weather conditions. In order to quantify the actual time and energy coverage during the campaign on Mrk 501, Figure 6 shows the exposure time as a function of the energy range for the instruments/observations used to produce the SED shown in Figure 8. Apart from the unprecedented energy coverage (including, for the first time, the GeV energy range from Fermi-LAT), the source was sampled quite uniformly with the various instruments participating in the campaign and, consequently, it is reasonable to consider the SED constructed below as the actual average (typical) SED of Mrk 501 during the time interval covered by this multifrequency campaign. The largest non-uniformity in the sampling of the source comes from the Cherenkov Telescopes, which are the instruments most sensitive to weather conditions. Moreover, while there are many radio/optical instruments spread all over the globe, there are only three Cherenkov Telescope observatories in the northern hemisphere we could utilize (MAGIC, VERITAS, Whipple). Hence, the impact of observing conditions was more important to the coverage at the VHE γ-ray energies.

Zoom In Zoom Out Reset image size

We note that Figure 6 shows the MAGIC, VERITAS, and Whipple coverage at VHE γ-ray energies, but only the MAGIC and VERITAS observations were used to produce the spectra shown in Figure 8. The more extensive (120 hr), but less sensitive, Whipple data (shown as gray boxes in Figure 6) were primarily taken to determine the light curve (Pichel et al. 2009) and a re-optimization was required to derive the spectrum which will be reported elsewhere.

In the following paragraphs, we briefly discuss the procedures used in the data analysis of the instruments participating in the campaign. The analysis of the Fermi-LAT data was described in Section 2, and the results obtained will be described in detail in Section 5.2.

5.1.1. Radio Instruments

5.1.2. Optical and Near-IR Instruments

5.1.3. Swift-UVOT

5.1.4. Swift-XRT

5.1.5. RXTE-PCA

5.1.6. Swift-BAT

5.1.7. MAGIC

5.1.8. VERITAS

The Mrk 501 spectrum measured by Fermi-LAT, integrated during the time interval of the multifrequency campaign, is shown in the panel (b) of Figure 7. The spectrum can be described by a power-law function with the photon index 1.74 ± 0.05. The flux data points resulting from the analysis in differential energy ranges are within 1σ–2σ of the power-law fit result; this is an independent indication that a single power-law function is a good representation of the spectrum during the multifrequency campaign. On the other hand, the shape of the spectrum depicted by the differential energy flux data points suggests the possibility of a concave spectrum. As was discussed in Sections 3 and 4 (see Figures 1 and 5), Mrk 501 showed substantial spectral variability during the time period covered by the multifrequency campaign, with some 30-day time intervals characterized by relatively soft spectra (photon index ∼2 for the 30-day intervals MJD 54892–54922 and MJD 54922–54952) and others by relatively hard spectra (photon index ∼1.6 for the 30-day intervals MJD 54952–54982, MJD 54982–55012, and MJD 55012–55042). Panel (b) of Figure 7 presents the average spectrum over those time intervals, and hence it would not be surprising to see two slopes (instead of one) in the spectrum. In order to evaluate this possibility, a broken power-law fit was applied, yielding indices of 1.86 ± 0.08 and 1.44 ± 0.14 below and above a break energy of 10 ± 3 GeV, respectively. The likelihood ratio of the broken power law and the power law is 2.2. Given that the broken power law has two additional degrees of freedom, this indicates that the broken power law is not statistically preferred over the single power-law function.

Zoom In Zoom Out Reset image size

For comparison purposes we also computed the spectra for time intervals before and after the multifrequency campaign (MJD 54683–54901 and MJD 55044–55162).¹⁵⁵ These two spectra, shown in panels (a) and (c) of Figure 7, can both be described satisfactorily by single power-law functions with photon indices 1.82 ± 0.06 and 1.80 ± 0.08. Note that the two spectra are perfectly compatible with each other, which is consistent with the relatively small flux/spectral variability shown in Figures 1 and 2 for those time periods.

The average broadband SED of Mrk 501 resulting from our 4.5 month long multifrequency campaign is shown in Figure 8. The TeV data from MAGIC and VERITAS have been corrected for the absorption in the EBL using the particular EBL model by Franceschini et al. (2008). The corrections given by the other low-EBL-level models (Kneiske et al. 2004; Gilmore et al. 2009; Finke et al. 2010) are very similar for the low redshift of Mrk 501 (z = 0.034). The attenuation factor at a photon energy of 6 TeV (the highest energy detected from Mrk 501 during this campaign) is in the range $e^{-\tau _{\gamma \gamma }} \simeq 0.4\hbox{--}0.5$, and smaller at lower energies.

Zoom In Zoom Out Reset image size

During the campaign, as already noted above, the source did not show large flux variations like those recorded by EGRET in 1996, or those measured by X-ray and TeV instruments in 1997. Nevertheless, significant flux and spectral variations at γ-ray energies occurred in the time interval MJD 54905–55044. The largest flux variation during the campaign was observed at TeV energies during the time interval MJD 54952.9–54955.9, when VERITAS measured a flux about five times higher than the average one during the campaign. Because of the remarkable difference with respect to the rest of the analyzed exposure, these observations were excluded from the data set used to compute the average VERITAS spectrum for the campaign; the three-day "flaring-state" spectrum (2.4 hr of observation) is presented separately in Figure 8. Such a remarkable flux enhancement was not observed in the other energy ranges and hence Figure 8 shows only the averaged spectra for the other instruments.¹⁵⁶

The top panel in Figure 9 shows a zoom of the high-energy bump depicted in Figure 8. The last two energy bins from Fermi (60–160 and 160–400 GeV) are systematically above (1σ–2σ) the measured/extrapolated spectrum from MAGIC and VERITAS. Even though this mismatch is not statistically significant, we believe that the spectral variability observed during the 4.5 month long campaign (see Sections 4 and 5.2) could be the origin of such a difference. Because Fermi-LAT operates in a survey mode, Mrk 501 is constantly monitored at GeV energies,¹⁵⁷ while this is not the case for the other instruments which typically sampled the source for ⩽1 hr every five days approximately. Moreover, because of bad weather or moonlight conditions, the monitoring at the TeV energies with Cherenkov telescopes was even less regular than that at lower frequencies. Therefore, Fermi-LAT may have measured high activity that was missed by the other instruments. Indeed, the 2.4 hr high-flux spectrum from VERITAS depicted in Figure 8 (which was obtained during the three-day interval MJD 54952.9–54955.9) demonstrates that, during the multifrequency campaign, there were time periods with substantially (factor of five) higher TeV activity. It is possible that the highest energy LAT observations (⩾50 GeV) include high TeV flux states which occurred while the IACTs were not observing.

Zoom In Zoom Out Reset image size

If the flaring activity occurred only at the highest photon energies, then the computed Fermi-LAT flux (>0.3 GeV) would not change very much and the effect might only be visible in the measured power-law photon index. This seems to be the case in the presented data set. As was shown in Figure 5, the 30-day intervals MJD 54922–54952 and MJD 54952–54982 have photon fluxes above 0.3 GeV of (3.9 ± 0.6) × 10⁻⁸ ph cm⁻² s⁻¹ and (3.6 ± 0.5) × 10⁻⁸ ph cm⁻² s⁻¹, while their photon indices are 2.10 ± 0.13 and 1.63 ± 0.09, respectively. Therefore, the spectral information (together with the enhanced photon flux) indicates the presence of flaring activity at the highest γ-ray energies during the second 30-day time period. Besides the factor ∼5 VHE flux enhancement recorded by VERITAS and Whipple at the beginning of the time interval MJD 54952–54982, MAGIC, and Whipple also recorded a factor ∼2 VHE flux enhancement at the end of this 30-day time interval (see preliminary fluxes reported in Paneque 2009; Pichel et al. 2009). This flux enhancement was measured for the time interval MJD 54975–54977, but there were no VHE measurements during the period MJD 54970.5–54975.0. Thus, the average Fermi-LAT spectrum could have been affected by elevated VHE activity during the 30-day time interval MJD 54952–54982, which was only partly covered by the IACTs participating in the campaign.

For illustrative purposes, in the bottom panel of Figure 9, we show separately the Fermi-LAT spectra for the 30-day time interval MJD 54952–54982 (high photon flux and hard spectrum), and for the rest of the campaign. It is interesting to note that the Fermi-LAT spectrum without the 30-day time interval MJD 54952–54982 (blue data points in the bottom panel of Figure 9) agrees perfectly with the VHE spectrum measured by IACTs. We also want to point out that the power-law fit to the Fermi-LAT spectrum without the 30-day interval MJD 54952–54982 gave a photon flux above 0.3 GeV of (2.62 ± 0.25) × 10⁻⁸ ph cm⁻² s⁻¹ with a photon index of 1.78 ± 0.07, which is statistically compatible with the results for the power-law fit to the Fermi-LAT data from the entire campaign (see panel (b) in Figure 7). As discussed above, the flaring activity occurred mostly at the highest energies, where the (relatively) low photon count has little impact on the overall power-law fit performed above 0.3 GeV.

This is the most complete quasi-simultaneous SED ever collected for Mrk 501, or for any other TeV-emitting BL Lac object (see also A. A. Abdo et al. 2011, in preparation). At the highest energies, the combination of Fermi and MAGIC/VERITAS data allows us to measure, for the first time, the high-energy bump without any spectral gap. The low-energy spectral component is also very well characterized with Swift-UVOT, Swift-XRT, and RXTE-PCA data, covering the peak of the synchrotron continuum. The only (large) region of the SED with no available data corresponds to the photon energy range 200 keV–100 MeV, where the sensitivity of current instruments is not good enough to detect Mrk 501. It is worth stressing that the excellent agreement in the overlapping energies among the various instruments (which had somewhat different time coverages) indicates that the collected data are representative of the truly average SED during the multi-instrument campaign.

Let us assume that the emitting region is a homogeneous and roughly spherically symmetric moving blob, with radius R and comoving volume V' ≃ (4π/3) R³. For this, we evaluate the comoving synchrotron and SSC emissivities, $\nu ^{\prime }\!j^{\prime }_{\nu ^{\prime }}$, assuming isotropic distributions of ultrarelativistic electrons and synchrotron photons in the rest frame of the emitting region. Thus, we use the exact synchrotron and inverse-Compton kernels (the latter one valid in both Thomson and Klein–Nishina regimes), as given in Crusius & Schlickeiser (1986) and Blumenthal & Gould (1970), respectively. The intrinsic monochromatic synchrotron and SSC luminosities are then $\nu ^{\prime }\!L^{\prime }_{\nu ^{\prime }} = 4 \pi \, V^{\prime } \, \nu ^{\prime }\!j^{\prime }_{\nu ^{\prime }}$, while the observed monochromatic flux densities (measured in erg cm⁻² s⁻¹) can be found as

$$\begin{equation} \nu F_{\nu } = {\delta ^4 \over 4 \pi \, d_L^2} \, [\nu ^{\prime }\!L^{\prime }_{\nu ^{\prime }}]_{\nu ^{\prime } = \nu \, (1+z)/\delta } \simeq {4 \pi \, \delta ^4 R^3 \over 3 \, d_L^2} \, [\nu ^{\prime }\!j^{\prime }_{\nu ^{\prime }}]_{\nu ^{\prime } = \nu \, (1+z)/\delta } \,, \end{equation} \tag{ 2 }$$

where δ is the jet Doppler factor, z = 0.034 is the source redshift, and d_L = 142 Mpc is the luminosity distance to Mrk 501. In order to evaluate the comoving emissivities $\nu ^{\prime }\!j^{\prime }_{\nu ^{\prime }}$, the electron energy distribution n'_e(γ) has to be specified. For this, we assume a general power-law form between the minimum and maximum electron energies, γ_min and γ_max, allowing for multiple spectral breaks in between, as well as for an exponential cut-off above γ_max. In fact, the broadband data set for Mrk 501 requires two different electron break energies, and hence we take the electron energy distribution in the form

$$\begin{equation} n^{\prime }_e(\gamma) \propto \left\{\begin{array}{@{}lcl}\gamma ^{-s_1} & {\rm for} & \gamma _{\rm min} \le \gamma < \gamma _{\rm {br},\,1} \\ \gamma ^{-s_2} & {\rm for} & \gamma _{\rm {br},\,1} \le \gamma < \gamma _{\rm {br},\,2} \\ \gamma ^{-s_3} \, \exp \left[-\gamma /\gamma _{\rm max}\right] & {\rm for} & \gamma _{\rm {br},\,2} \le \gamma, \end{array} \right. \end{equation} \tag{ 3 }$$

with the normalization expressed in terms of the equipartition parameter (the ratio of the comoving electron and magnetic field energy densities), namely

$$\begin{equation} \eta _e \equiv {U^{\prime }_e \over U^{\prime }_B} = {\int \, \gamma \, m_e c^2 \, n^{\prime }_e(\gamma) \, d\gamma \over B^2 / 8 \pi } \,. \end{equation} \tag{ 4 }$$

The measured SED is hardly compatible with a simpler form of the electron distribution with only one break and an exponential cutoff. However, some smoothly curved spectral shape might perhaps be an alternative representation of the electron spectrum (e.g., Stawarz & Petrosian 2008; Tramacere et al. 2009).

The model adopted is thus characterized by four main free parameters (B, R, δ, and η_e), plus seven additional ones related to the electron energy distribution (γ_min, γ_br, 1, γ_br, 2, γ_max, s₁, s₂, and s₃). These seven additional parameters are determined by the spectral shape of the non-thermal emission continuum probed by the observations, predominantly by the spectral shape of the synchrotron bump (rather than the inverse-Compton bump), and depend only slightly on the particular choice of the magnetic field B and the Doppler factor δ within the allowed range.¹⁵⁸ There is a substantial degeneracy regarding the four main free parameters: the average emission spectrum of Mrk 501 may be fitted by different combinations of B, R, δ, and η_e with little variation in the shape of the electron energy distribution. Note that, for example, [νF_ν]_syn ∝ R³ η_e, but at the same time [νF_ν]_ssc ∝ R⁴ η²_e. We can attempt to reduce this degeneracy by assuming that the observed main variability timescale is related to the size of the emission region and its Doppler factor according to the formula

$$\begin{equation} t_{\rm var} \simeq {(1+z) \, R \over c \, \delta } \,. \end{equation} \tag{ 5 }$$

The multifrequency data collected during the 4.5 month campaign (see Section 5) allow us to study the variability of Mrk 501 on timescales from months to a few days. We found that, during this time period, the multifrequency activity typically varied on a timescale of 5–10 days, with the exception of a few particular epochs when the source became very active in VHE γ-rays, and flux variations with timescales of a day or shorter were found at TeV energies. Nevertheless, it is important to stress that several authors concluded in the past that the dominant emission site of Mrk 501 is characterized by variability timescales longer than one day (see Kataoka et al. 2001, for a comprehensive study of the Mrk 501 variability in X-rays), and that the power in the intraday flickering of this source is small, in agreement with the results of our campaign. Nevertheless, one should keep in mind that this object is known for showing sporadic but extreme changes in its activity that can give flux variations on timescales as short as a few minutes (Albert et al. 2007a). In this work, we aim to model the average/typical behavior of Mrk 501 (corresponding to the 4.5 month campaign) rather than specific/short periods with outstanding activity, and hence we constrained the minimum (typical) variability timescale t_var in the model to the range 1–5 days.

Even with t_var fixed as discussed above, the reconstructed SED of Mrk 501 may be fitted by different combinations of B, R, δ, and η_e. Such a degeneracy between the main model parameters is an inevitable feature of the SSC modeling of blazars (e.g., Kataoka et al. 1999), and it is therefore necessary to impose additional constraints on the physical parameters of the dominant emission zone. Here, we argue that such constraints follow from the requirement for the electron energy distribution to be in agreement with the one resulting from the simplest prescription of the energy evolution of the radiating electrons within the emission region, as discussed below.

The idea of separating the sites for the particle acceleration and emission processes is commonly invoked in modeling different astrophysical sources of high-energy radiation, and blazar jets in particular. Such a procedure is not always justified, because interactions of ultrarelativistic particles with the magnetic field (leading to particle diffusion and convection in momentum space) are generally accompanied by particle radiative losses (and vice versa). On the other hand, if the characteristic timescale for energy gains is much shorter than the timescales for radiative cooling (t'_rad) or escape (t'_esc) from the system, the particle acceleration processes may be indeed approximated as being "instantaneous," and may be modeled by a single injection term ${\dot{Q}}(\gamma)$ in the simplified version of the kinetic equation

$$\begin{equation} {\partial n^{\prime }_e(\gamma) \over \partial t} = - {\partial \over \partial \gamma } \left[{\gamma \, n_e(\gamma) \over t^{\prime }_{\rm rad}(\gamma)}\right] - {n_e(\gamma) \over t^{\prime }_{\rm esc}} + {\dot{Q}}(\gamma), \end{equation} \tag{ 6 }$$

describing a very particular scenario for the energy evolution of the radiating ultrarelativistic electrons.

It is widely believed that the above equation is a good approximation for the energy evolution of particles undergoing diffusive (first-order Fermi) shock acceleration, and cooling radiatively in the downstream region of the shock. In such a case, the term ${\dot{Q}}(\gamma)$ specifies the energy spectrum and the injection rate of the electrons freshly accelerated at the shock front and not affected by radiative losses, while the escape term corresponds to the energy-independent dynamical timescale for the advection of the radiating particles from the downstream region of a given size R, namely t'_esc ≃ t'_dyn ≃ R/c. The steady-state electron energy distribution is then very roughly $n^{\prime }_e(\gamma) \sim t^{\prime }_{\rm dyn} \, {\dot{Q}}(\gamma)$ below the critical energy for which t'_rad(γ) = t'_dyn, and $n^{\prime }_e(\gamma) \sim t^{\prime }_{\rm rad}(\gamma) \, {\dot{Q}}(\gamma)$ above this energy. Note that in the case of a power-law injection ${\dot{Q}}(\gamma) \propto \gamma ^{-s}$ and a homogeneous emission region with dominant radiative losses of the synchrotron type, t'_rad(γ) ∝ γ⁻¹, the injected electron spectrum is expected to steepen by Δs = 1 above the critical "cooling break" energy. This provides us with the additional constraint on the free model parameters for Mrk 501: namely, we require that the position of the second break in the electron energy distribution needed to fit the reconstructed SED, γ_br 2, should correspond to the location of the cooling break for a given chosen set of the model free parameters.

Figure 10 (black curves) shows the resulting SSC model fit (summarized in Table 2) to the averaged broadband emission spectrum of Mrk 501, which was obtained for the following parameters: B = 0.015 G, R = 1.3 × 10¹⁷ cm, δ = 12, η_e = 56, γ_min = 600, γ_br, 1 = 4 × 10⁴, γ_br, 2 = 9 × 10⁵, γ_max = 1.5 × 10⁷, s₁ = 2.2, s₂ = 2.7, and s₃ = 3.65. The overall good agreement of the model with the data is further discussed in Section 6.2. Here, we note that, for these model parameters, synchrotron self-absorption effects are important only below 1 GHz, where we do not have observations.¹⁵⁹ We also emphasize that with all the aforementioned constraints and for a given spectral shape of the synchrotron continuum (including all the data points aimed to be fitted by the model, as discussed below), and thus for a fixed spectral shape of the electron energy distribution (modulo critical electron Lorentz factors scaling as $\propto 1/\sqrt{B\,\delta }$), the allowed range for the free parameters of the model is relatively narrow. Namely, for the variability timescale between one and five days, the main model parameters may change within the ranges R ≃ (0.35–1.45) × 10¹⁷ cm, δ ≃ 11–14, and B ≃ 0.01–0.04 G. The parameter η_e depends predominantly on the minimum Lorentz factor of the radiating electrons. Hence, it is determined uniquely as η_e ≃ 50 with the submillimeter flux included in the fitted data set. Only with a different prescription for the spectral shape of the electron energy distribution could the main free parameters of the model be significantly different from those given above.

Zoom In Zoom Out Reset image size

Despite the absence of any fast variability during this multifrequency campaign (apart from the already discussed isolated three-day-long flare), Mrk 501 is known for the extremely rapid flux changes at the highest observed photon energies (e.g., Albert et al. 2007a). Hence, it is interesting to check whether any shorter than few-day-long variability timescales can be accommodated in the framework of the simplest SSC model applied here for the collected data set. In order to do that, we decreased the minimum variability time scale by one order of magnitude (from four days to 0.4 days), and tried to model the data. A satisfactory fit could be obtained with those modified parameters, but only when we relaxed the requirement for the electron energy distribution to be in agreement with the one following from the steady-state solution to Equation (6), and in particular the resulting constraint for the second break in the electron spectrum to be equal to the cooling break. This "alternative" model fit is shown in Figure 10 (red curves) together with the "best" model fit discussed above. The resulting model parameters for the "alternative" fit are B = 0.03 G, R = 2 × 10¹⁶ cm, δ = 22, η_e = 130, γ_min = 300, γ_br, 1 = 3 × 10⁴, γ_br, 2 = 5 × 10⁵, γ_max = 3 × 10⁶, s₁ = 2.2, s₂ = 2.7, and s₃ = 3.5. This particular parameter set—which should be considered as an illustrative one only—would be therefore consistent with a minimum variability timescale of 0.36 days, but at the price of much larger departures from the energy equipartition (η_e > 100). The other source parameters, on the other hand, would change only slightly (see Table 2). Because of the mismatch (by a factor ∼3) between the location of the cooling break and the second break in the electron distribution, we consider this "alternative" fit less consistent with the hypothesis of the steady-state homogeneous one-zone SSC scenario, which is the framework with which we chose to model the broadband SED of Mrk 501 emerging from the campaign.

The low-frequency radio observations performed with single-dish instruments have a relatively large contamination from non-blazar emission due to the underlying extended jet component, and hence they only provide upper limits for the radio flux density of the blazar emission zone. On the other hand, the flux measurements by the interferometric instruments (such as VLBA), especially the ones corresponding to the core, provide us with the radio flux density from a region that is not much larger than the blazar emission region.

The radio flux densities from interferometric observations (from the VLBA core) are expected to be close upper limits to the radio continuum of the blazar emission component. The estimated size of the partially resolved VLBA core of Mrk 501 at 15 GHz and 43 GHz is ≃0.14–0.18 mas ≃2.9–3.7 ×10¹⁷ cm (with the appropriate conversion scale 0.67 pc mas⁻¹). The VLBA size estimation is the FWHM of a Gaussian representing the brightness distribution of the blob, which could be approximated as 0.9 times the radius of a corresponding spherical blob (Marscher 1983). That implies that the size of the VLBA core is only a factor of two to three larger than the emission region in our SSC model fit (R = 1.3 × 10¹⁷ cm). Therefore, it is reasonable to assume that the radio flux density from the VLBA core is indeed dominated by the radio flux density of the blazar emission. Forthcoming multi-band correlation studies (in particular VLBA and SMA radio with the γ-rays from Fermi-LAT) will shed light on this particular subject. Interestingly, the magnetic field claimed for the partially resolved radio core of Mrk 501 (which has a size of ≲0.2 mas) and its submilliarcsec jet, namely B ≃ (10–30) mG (Giroletti et al. 2004, 2008), is in very good agreement with the value emerging from our model fits (15 mG), assuring self-consistency of the approach adopted.

In addition to this, in the modeling we also aimed at matching the submillimeter flux of Mrk 501, given at the observed frequency of 225 GHz, assuming that it represents the low-frequency tail of the optically thin synchrotron blazar component. One should emphasize in this context that it is not clear if the blazar emission zone is in general located deep within the millimeter photosphere or not. However, the broadband variability of luminous blazars of the FSRQ type indicates that there is a significant overlap of the blazar zone with a region where the jet becomes optically thin at millimeter wavelengths (as discussed by Sikora et al. 2008, for the particular case of the blazar 3C 454.3). We have assumed that the same holds for BL Lac objects.

The IR/optical flux measurements in the range ∼(1–10) × 10¹⁴ Hz represent the starlight of the host galaxy and hence they should be excluded when fitting the non-thermal emission of Mrk 501. We modeled these data points with the template spectrum of an elliptical galaxy instead (including only the dominant stellar component due to the evolved red giants, as discussed in Silva et al. 1998), obtaining a very good match (see the dotted line in Figure 10) for the bolometric starlight luminosity L_star ≃ 3 × 10⁴⁴ erg s⁻¹. Such a luminosity is in fact expected for the elliptical host of a BL Lac object. The model spectrum of the galaxy falls off very rapidly above 5 × 10¹⁴ Hz, while the three UV data points (above 10¹⁵ Hz) indicate a prominent, flat power-law UV excess over the starlight emission. Therefore, it is reasonable to assume that the observed UV fluxes correspond to the synchrotron (blazar-type) emission of Mrk 501 and, consequently, we used them in our model fit. However, many elliptical galaxies do reveal in their spectra the so-called UV upturn, or UV excess, whose origin is not well known, but which is presumably related to the starlight continuum (most likely due to young stars from the residual star-forming activity within the central region of a galaxy) rather than to non-thermal (jet-related) emission processes (see, e.g., Code & Welch 1979; Atlee et al. 2009). Hence, it is possible that the UV data points provided here include some additional contamination from the stellar emission, and as such might be considered as upper limits for the synchrotron radiation of the Mrk 501 jet.

The observed X-ray spectrum of Mrk 501 agrees very well with the SSC model fit, except for a small but statistically significant discrepancy between the model curve and the first two data points provided by Swift-XRT, which correspond to the energy range 0.3–0.6 keV. As pointed out in Section 5.1, the Swift-XRT data had to be corrected for a residual energy offset which affects the lowest energies. The correction for this effect could introduce some systematic differences with respect to the actual fluxes detected at those energies. These low-energy X-ray data points might also be influenced by intrinsic absorption of the X-ray photons within the gaseous environment of Mrk 501 nucleus, as suggested by the earlier studies with the ASCA satellite (see Kataoka et al. 1999). As a result, the small discrepancy between the data and the model curve within the range 0.3–0.6 keV can be ignored in the modeling.

The agreement between the applied SSC model and the γ-ray data is also very good. In particular, the model returns the γ-ray photon index 1.78 in the energy range 0.3–30 GeV, which can be compared with the one resulting from the power-law fit to the Fermi-LAT data above 0.3 GeV, namely 1.74 ± 0.05. However, the last two energy bins from Fermi (60–160 and 160–400 GeV) are systematically above (2σ) the model curves, as well as above the averaged spectrum reported by MAGIC and VERITAS. A possible reason for mismatch between the average Fermi-LAT spectrum and the one from MAGIC/VERITAS was discussed in Section 5.3.

The values for the emission region size R = 1.3 × 10¹⁷ cm and the jet Doppler factor δ = 12 emerging from our SSC model fit give a minimum (typical) variability timescale of t_var ≃ (1 + z) R/c δ ∼ 4 days, which is consistent with the variability observed during the campaign and with previous studies of the X-ray activity of Mrk 501 (Kataoka et al. 2001). At this point, it is necessary to determine whether an emission region characterized by these values of R and δ is optically thin to internal two-photon pair creation γγ → e⁺e⁻ for the highest TeV energies observed during the campaign. We now affirm pair transparency due to insufficient density of soft target photons.

Since Mrk 501 is a cosmologically local object, pair conversion in the EBL is not expected to prevent its multi-TeV photons from reaching the Earth, although the impact of this process is not negligible, as mentioned in Section 5. Therefore, dealing with a nearby source allows us to focus mostly on the intrinsic absorption processes, rather than on the cosmological, EBL-related, attenuation of the γ-ray emission. Moreover, because of the absence (or weakness) of accretion-disk-related circumnuclear photon fields in BL Lac objects like Mrk 501, we only need to consider photon–photon pair production involving photon fields internal to the jet emission site. The analysis is therefore simpler than in the case of FSRQs, where the attenuation of high-energy γ-ray fluxes is dominated by interactions with photon fields external to the jet—such as those provided by the broad line regions or tori—for which the exact spatial distribution is still under debate.

Pair-creation optical depths can be estimated as follows. Using the δ-function approximation for the photon–photon annihilation cross-section (Zdziarski & Lightman 1985), σ_γγ(ε'₀, ε'_γ) ≃ 0.2 σ_T ε'₀ δ[ε'₀ − (2m²_ec⁴/ε'_γ)], the corresponding optical depth for a γ-ray photon with observed energy ε_γ = δ ε'_γ/(1 + z) interacting with a jet-originating soft photon with observed energy

$$\begin{equation} \varepsilon _0 = {\delta \, \varepsilon ^{\prime }_0 \over 1+z} \simeq {2 \, \delta ^2 m_e^2 c^4 \over \varepsilon _{\gamma } \, (1+z)^2} \simeq 50 \, \left({\delta \over 10}\right)^2 \left({ \varepsilon _{\gamma } \over {\rm TeV}}\right)^{-1}\,{\rm eV} \,, \end{equation} \tag{ 7 }$$

may be found as

$$\begin{eqnarray} \tau _{\gamma \gamma } &\simeq &\int ^R ds \, \int _{m_ec^2/\varepsilon _0} d\varepsilon _0 \, n^{\prime }_0(\varepsilon _0) \, \sigma _{\gamma \gamma }(\varepsilon _0, \varepsilon _\gamma)\nonumber \\ & \sim & 0.2 \, \sigma _{\rm T} \, R \, \varepsilon ^{\prime }_0 \, n^{\prime }_0(\varepsilon _0) \,, \end{eqnarray} \tag{ 8 }$$

where n'₀(ε'₀) is the differential comoving number density of soft photons. Noting that ε'²₀ n'₀(ε'₀) = L'₀/4π R²c, where L'₀ is the intrinsic monochromatic luminosity at photon energy ε'₀, we obtain

$$\begin{eqnarray} \tau _{\gamma \gamma } &\simeq &{\sigma _{\rm T} \, d_L^2 F_0 \, \varepsilon _{\gamma } (1+z) \over 10 \, R \, m_e^2 c^5 \, \delta ^5} \simeq 0.001 \, \left({ \varepsilon _{\gamma } \over {\rm TeV}}\right) \nonumber \\ && \times \left({F_0 \over 10^{-11}\,{\rm erg\;cm^{-2}\;s^{-1}}}\right) \, \left({R \over 10^{17}\,{\rm cm}}\right)^{-1} \left({\delta \over 10}\right)^{-5},\nonumber \\ \end{eqnarray} \tag{ 9 }$$

where F₀ = L₀/4πd²_L is the observed monochromatic flux energy density as measured at the observed photon energy ε₀. Thus, for 5 TeV γ-rays and the model parameters discussed (implying the observed ε₀ = 15 eV flux of Mrk 501 roughly F₀ ≃ 3.2 × 10⁻¹¹ erg s⁻¹ cm⁻²), one has τ_γγ(5 TeV) ≃ 0.005. Therefore, the values of R and δ from our SSC model fit do not need to be adjusted to take into account the influence of spectral modifications due to pair attenuation. Note that such opacity effects, studied extensively in the context of γ-ray bursts, generally yield a broken power law for the spectral form, with the position and magnitude of the break fixed by the pair-production kinematics (e.g., Baring 2006 and references therein). The broad-band continuum of Mrk 501, and in particular its relatively flat spectrum VHE γ-ray segment, is inconsistent with such an expected break. This deduction is in agreement with the above derived transparency of the emitting region for high-energy γ-ray photons.

Next we evaluated the "monoenergetic" comoving energy density of ultrarelativistic electrons for a given electron Lorentz factor,

$$\begin{equation} \gamma \, U^{\prime }_e(\gamma) \equiv \gamma ^2 m_e c^2 \, n^{\prime }_e(\gamma) \,, \end{equation} \tag{ 10 }$$

and this is shown in Figure 11 (solid black line). The total electron energy density is then U'_e = ∫U'_e(γ) dγ ≃ 5 × 10⁻⁴ erg cm⁻³. As shown, most of the energy is stored in the lowest energy particles (γ_min ≃ 600). For comparison, the comoving energy density of the magnetic field and that of the synchrotron photons are plotted in the figure as well (horizontal solid red line and dotted blue line, respectively). These two quantities are approximately equal, namely U'_B = B²/8π ≃ 0.9 × 10⁻⁵ erg cm⁻³ and

$$\begin{equation} U^{\prime }_{\rm syn} = {4 \pi \, R \over 3 \, c} \, \int j^{\prime }_{\nu ^{\prime },\,\rm syn} \, d\nu ^{\prime } \simeq 0.9 \times 10^{-5} \, {\rm erg\,cm^{-3}} \,. \end{equation} \tag{ 11 }$$

In Figure 11, we also plot the comoving energy density of synchrotron photons which are inverse-Compton upscattered in the Thomson regime for a given electron Lorentz factor γ,

$$\begin{equation} U^{\prime }_{\rm syn/T}(\gamma) = {4 \pi \, R \over 3 \, c} \, \int ^{\nu ^{\prime }_{KN}\!(\gamma)} j^{\prime }_{\nu ^{\prime },\,\rm syn} \, d\nu ^{\prime } \end{equation} \tag{ 12 }$$

(dashed blue line), where ν'_KN(γ) ≡ m_ec²/4γ h. Because of the well-known suppression of the inverse-Compton scattering rate in the Klein–Nishina regime, the scattering in the Thomson regime dominates the inverse-Compton energy losses.¹⁶⁰ Hence, one may conclude that even though the total energy density of the synchrotron photons in the jet rest frame is comparable to the comoving energy density of the magnetic field (U'_syn ≃ U'_B), the dominant radiative cooling for all the electrons is due to synchrotron emission, since U'_syn/T < U'_B for any γ.

Zoom In Zoom Out Reset image size

The timescale for synchrotron cooling may be evaluated as

$$\begin{equation} t^{\prime }_{\rm syn} \simeq {3 m_e c \over 4 \sigma _T \gamma \, U^{\prime }_B} \simeq 4 \, \left({\gamma \over 10^7}\right)^{-1} \, {\rm days} \,. \end{equation} \tag{ 13 }$$

Hence, t'_rad ≃ t'_syn equals the dynamical timescale of the emitting region, t'_dyn ≃ R/c, for electron Lorentz factor γ ≃ 8 × 10⁵, i.e., close to the second electron break energy γ_br, 2. Also, the difference between the spectral indices below and above the break energy γ_br, 2 determined by our modeling, namely Δs_3/2 = s₃ − s₂ = 0.95, is very close to the "classical" synchrotron cooling break Δs = 1 expected for a uniform emission region, as discussed in Section 6.1. This agreement, which justifies at some level the assumed homogeneity of the emission zone, was in fact the additional constraint imposed on the model to break the degeneracy between the main free parameters. Note that in such a case the first break in the electron energy distribution around electron energy γ_br, 1 = 4 × 10⁴ is related to the nature of the underlying particle acceleration process. We come back to this issue in Section 7.2,

Another interesting result from our model fit comes from the evaluation of the mean energy of the electrons responsible for the observed non-thermal emission of Mrk 501. In particular, the mean electron Lorentz factor is

$$\begin{equation} \langle \gamma \rangle \equiv {\int \gamma \, n^{\prime }_e(\gamma) \, d\gamma \over \int n^{\prime }_e(\gamma) \, d\gamma } \simeq 2400. \end{equation} \tag{ 14 }$$

This value, which is determined predominantly by the minimum electron energy γ_min = 600 and by the position of the first break in the electron energy distribution, is comparable to the proton-to-electron mass ratio m_p/m_e. In other words, the mean energy of ultrarelativistic electrons within the blazar emission zone of Mrk 501 is comparable to the energy of non-relativistic/mildly-relativistic (cold) protons. This topic will be discussed further in Section 7.2 as well.

The analysis presented allows us also to access the global energetics of the Mrk 501 jet. In particular, with the given energy densities U'_e and U'_B, we evaluate the total kinetic powers of the jet stored in ultrarelativistic electrons and magnetic field as

$$\begin{equation} L_e = \pi R^2 c \Gamma ^2 \, U^{\prime }_e \simeq 10^{44} \, {\rm erg\,s^{-1}} \,, \end{equation} \tag{ 15 }$$

and

$$\begin{equation} L_B = \pi R^2 c \Gamma ^2 \, U^{\prime }_B \simeq 2 \times 10^{42} \, {\rm erg\,s^{-1}} \,, \end{equation} \tag{ 16 }$$

respectively. In the above expressions, we have assumed that the emission region analyzed occupies the whole cross-sectional area of the outflow, and that the jet propagates at sufficiently small viewing angle that the bulk Lorentz factor of the jet equals the jet Doppler factor emerging from our modeling, Γ = δ. These assumptions are justified in the framework of the one-zone homogeneous SSC scenario. Moreover, assuming one electron–proton pair per electron–positron pair within the emission region (see Celotti & Ghisellini 2008), or equivalently N'_p ≃ N'_e/3, where the total comoving number density of the jet leptons is

$$\begin{equation} N^{\prime }_e = \int n^{\prime }_e(\gamma) \, d\gamma \simeq 0.26 \, {\rm cm^{-3}} \,, \end{equation} \tag{ 17 }$$

we obtain the comoving energy density of the jet protons U'_p = 〈γ_p〉 m_pc²N'_e/3, and hence the proton kinetic flux L_p = πR²cΓ² U'_p ≃ 0.3 〈γ_p〉 10⁴⁴ erg s⁻¹. This is comparable to the kinetic power carried out by the leptons for mean proton Lorentz factor 〈γ_p〉 ≃ 4 (see Equation (15)). It means that, within the dominant emission zone of Mrk 501 (at least during non-flaring activity), ultrarelativistic electrons and mildly relativistic protons, if comparable in number, are in approximate energy equipartition, and their energy dominates that of the jet magnetic field by two orders of magnitude. It is important to compare this result with the case of powerful blazars of the FSRQ type, for which the relatively low mean energy of the radiating electrons, 〈γ〉 ≪ 10³, assures dynamical domination of cold protons even for a smaller proton content N'_p/N'_e ≲ 0.1 (see the discussion in Sikora et al. 2009 and references therein).

Assuming 〈γ_p〉 ∼ 1 for simplicity, we find that the implied total jet power L_j = L_e + L_p + L_B ≃ 1.4 × 10⁴⁴ erg s^-1 constitutes only a small fraction of the Eddington luminosity L_Edd = 4π GM_BH m_pc/σ_T ≃ (1.1–4.4) × 10⁴⁷ erg s⁻¹ for the Mrk 501 black hole mass M_BH ≃ (0.9–3.5) × 10⁹M_☉ (Barth et al. 2002). In particular, our model implies L_j/L_Edd ∼ 10⁻³ in Mrk 501. In this context, detailed investigation of the emission-line radiative output from the Mrk 501 nucleus by Barth et al. (2002) allowed for an estimate of the bolometric, accretion-related luminosity as L_disk ≃ 2.4 × 10⁴³ erg s⁻¹, or L_disk/L_Edd ∼ 10⁻⁴. Such a relatively low luminosity is not surprising for BL Lac objects, which are believed to accrete at low, sub-Eddington rates (e.g., Ghisellini et al. 2009a). For a low-accretion rate AGN (i.e., those for which L_disk/L_Edd < 10⁻²) the expected radiative efficiency of the accretion disk is η_disk ≡ L_disk/L_acc < 0.1 (Narayan & Yi 1994; Blandford & Begelman 1999). Therefore, the jet power estimated here for Mrk 501 is comfortably smaller than the available power of the accreting matter L_acc.

Finally, we note that the total emitted radiative power is

$$\begin{equation} L_{\rm {em}} \simeq \Gamma ^2 \, (L^{\prime }_{\rm syn} + L^{\prime }_{\rm ssc}) = 4 \pi R^2 c \, \Gamma ^2 \, (U^{\prime }_{\rm syn} + U^{\prime }_{\rm ssc}) \sim 10^{43} \, {\rm erg\,s^{-1}} \,, \end{equation} \tag{ 18 }$$

where U'_syn is given in Equation (11) and the comoving energy density of γ-ray photons, U'_ssc, was evaluated in an analogous way as ≃1.7 × 10⁻⁶ erg cm^-3. This implies that the jet/blazar radiative efficiency was at the level of a few percent (L_em/L_j ≃ 0.07) during the period covered by the multifrequency campaign. Such a relatively low radiative efficiency is a common characteristic of blazar jets in general, typically claimed to be at the level of 1%–10% (see Celotti & Ghisellini 2008; Sikora et al. 2009). On the other hand, the isotropic synchrotron and SSC luminosities of Mrk 501 corresponding to the observed average flux levels are L_syn = δ⁴L'_syn ≃ 10⁴⁵ erg s⁻¹ and L_ssc = δ⁴L'_ssc ≃ 2 × 10⁴⁴ erg s⁻¹, respectively.

The results of the SSC modeling presented in the previous sections indicate that the energy spectrum of freshly accelerated electrons within the blazar emission zone of Mrk 501 is of the form ∝E^−2.2_e between electron energy E_e, min = γ_minm_ec² ∼ 0.3 GeV and E_e, br = γ_br, 1m_ec² ∼ 20 GeV, steepening to ∝E^−2.7_e above 20 GeV, such that the mean electron energy is 〈E_e〉 ≡ 〈γ〉 m_ec² ∼ 1 GeV. At this point, a natural question arises: is this electron distribution consistent with the particle spectrum expected for a diffusive shock acceleration process? Note in this context that the formation of a strong shock in the innermost parts of Mrk 501 might be expected around the location of the large bend (change in the position angle by 90°) observed in the outflow within the first few parsecs (projected) from the core (Edwards & Piner 2002; Piner et al. 2009). This distance scale could possibly be reconciled with the expected distance of the blazar emission zone from the center for the model parameters discussed, r ∼ R/θ_j ∼ 0.5 pc, for jet opening angle θ_j ≃ 1/Γ ≪ 1.

In order to address this question, let us first discuss the minimum electron energy implied by the modeling, E_e, min ∼ 0.3 GeV. In principle, electrons with lower energies may be present within the emission zone, although their energy distribution has to be very flat (possibly even inverted), in order not to overproduce the synchrotron radio photons and to account for the measured Fermi-LAT γ-ray continuum. Therefore, the constrained minimum electron energy marks the injection threshold for the main acceleration mechanism, meaning that only electrons with energies larger than E_e, min are picked up by this process to form the higher-energy (broken) power-law tail. Interestingly, the energy dissipation mechanisms operating at the shock fronts do introduce a particular characteristic (injection) energy scale, below which the particles are not freely able to cross the shock front. This energy scale depends on the shock microphysics, in particular on the thickness of the shock front. The shock thickness, in turn, is determined by the operating inertial length, or the diffusive mean free path of the radiating particles, or both. Such a scale depends critically on the constituents of the shocked plasma. For pure pair plasmas, only the electron thermal scale enters, and this sets E_e, min ∼ Γm_ec². In contrast, if there are approximately equal numbers of electrons and protons, the shock thickness can be relatively large. Diffusive shock acceleration can then operate on electrons only above a relatively high energy, establishing E_e, min ∼  Γm_pc². Here, represents some efficiency of the equilibration in the shock layer between shocked thermal protons and their electron counterparts, perhaps resulting from electrostatic potentials induced by charge separation of species of different masses (Baring & Summerlin 2007). Our multifrequency model fits suggest that ∼ 0.025, providing an important blazar shock diagnostic.

At even lower electron energies, other energization processes must play a dominant role (e.g., Hoshino et al. 1992), resulting in formation of very flat electron spectra (see the related discussion in Sikora et al. 2009). Which of these energy dissipation mechanisms are the most relevant, as well as how flat the particle spectra could be, are subjects of ongoing debates. Different models and numerical simulations presented in the literature indicate a wide possible range for the lowest-energy particle spectral indices (below E_e, min), from s_inj ≃ 1.0–1.5 down to s_inj ≃ −2, depending on the particular shock conditions and parameters (Amato & Arons 2006; Sironi & Spitkovsky 2009).

All in all, we argue that the relatively high minimum energy of the radiating electrons implied by the SSC modeling of the Mrk 501 broadband spectrum and the overall character of the electron energy evolution in this source are consistent with a proton-mediated shock scenario. In addition, the fact that the mean energy of the ultrarelativistic electrons is of the order of the proton rest energy, 〈E_e〉 ∼ m_pc², can be reconciled with such a model. Moreover, the constrained power-law slope of the electrons with energies E_e, min ⩽ E_e ⩽ E_e, br, namely s₁ = 2.2, seems to suggest a dominant role for diffusive shock acceleration above the injection energy E_e, min, as this value of the spectral index is often claimed in the literature for particles undergoing first-order Fermi acceleration at relativistic shock (Bednarz & Ostrowski 1998; Kirk et al. 2000; Achterberg et al. 2001). The caveat here is that this result for the "universal" particle spectral index holds only for particular conditions (namely, for ultrarelativistic shock velocities with highly turbulent conditions downstream of the shock: see the discussion in Ostrowski & Bednarz 2002), whereas in general a variety of particle spectra may result from the relativistic first-order Fermi mechanism, depending on the local magnetic field and turbulence parameters at the shock front, and the speed of the upstream flow (Kirk & Heavens 1989; Niemiec & Ostrowski 2004; Lemoine et al. 2006; Sironi & Spitkovsky 2009; Baring 2010). Nevertheless, the evidence for ultrarelativistic electrons with spectral index s₁ = 2.2 in the jet of Mrk 501 may be considered as an indication that the plasma conditions within the blazar emission zone allow for efficient diffusive shock acceleration (at least in this source), as described by the simplest asymptotic test-particle models, though only in a relatively narrow electron energy range.

If the relativistic shock acceleration plays a dominant role in the blazar Mrk 501 (as argued above), the observations and the model results impose important constraints on this mechanism, many aspects of which are still not well understood. Firstly, this process must be very efficient in the sense that all the electrons pre-accelerated/preheated to the energy E_e, min ∼ 0.3 GeV are picked up by the acceleration process so that a single electron component is formed above the injection threshold and there is no Maxwellian-like population of particles around E_e, min outnumbering the higher energy ones from the power-law tail.¹⁶¹ The second constraint is due to the presence of the spectral break Δs_2/1 = s₂ − s₁ ≃ 0.5 around electron energies E_e, br ∼ 20 GeV. As discussed in the previous section, this break cannot be simply a result of cooling or internal pair-attenuation effects, and hence must be accounted for by the acceleration mechanism. The discussion regarding the origin of this break—which may reflect variations in the global field orientation or turbulence levels sampled by particles of different energy—is beyond the scope of this paper. However, the presence of high-energy breaks in the electron energy distribution (intrinsic to the particle spectrum rather than forming due to cooling or absorption effects) seems to be a common property of relativistic jet sources observed by Fermi-LAT, such as the FSRQ objects 3C 454.3 and AO 0235+164 (Abdo et al. 2009, 2010a).

In Sections 3–4, we reported on the γ-ray flux and spectral variability of Mrk 501 as measured by the Fermi-LAT instrument during the first 16 months of operation. In this section, we discuss whether the observed spectral evolution can be accounted for by our simple one-zone SSC model.

Figure 12 (top panel) presents in more detail the GeV–TeV γ-ray spectrum of Mrk 501, together with the decomposition of the SSC model continuum. Here the contributions of different segments of the electron energy distribution are indicated by different colors. As shown, the low-energy electrons, γ_min ⩽ γ < γ_br, 1, which emit synchrotron photons up to ≃10¹⁵ Hz, dominate the production of γ-rays up to a few GeV (red line). The contribution of higher energy electrons with Lorentz factors γ_br, 1 ⩽ γ < γ_br, 2 is pronounced within the observed synchrotron range 10¹⁵ − 10¹⁸ Hz, and at γ-ray energies from a few GeV up to ∼TeV (green line). Finally, the highest energy tail of the electron energy distribution, γ ⩾ γ_br, 2, responsible for the observed hard-X-ray synchrotron continuum (>2 keV) in the fast cooling regime, generates the bulk of γ-rays with observed energies > TeV (blue line). Interestingly, even though any sharp breaks in the underlying electron energy distribution are "smeared out" into a smoothly curved spectral continuum due to the nature of the SSC emission, the average data set does support the presence of distinct low-energy and high-energy segments in the electron spectrum.

Zoom In Zoom Out Reset image size

It therefore seems reasonable to argue that the spectral variability of Mrk 501 observed by Fermi-LAT may be explained by postulating that the low-energy segment of the electron energy distribution (γ < γ_br, 1) is characterized by only small flux variations, while the high-energy electron tail (γ>γ_br, 1) varies more substantially. In such a scenario, some correlation might be expected between the fluxes in the UV-to-soft-X-ray photon energies from the synchrotron bump and the GeV–TeV fluxes from the inverse-Compton bump. This expectation is not inconsistent with the fact that we do not see any obvious relation between the ASM/BAT fluxes and the LAT (>2 GeV) fluxes (see Figure 2), because the electrons producing X-ray synchrotron photons above 2–3 keV contribute to the SSC emission mostly at the highest photon energies in the TeV range (see Figure 12). This issue will be studied in detail (on timescales of 1 month down to 5 days) in a forthcoming publication using the data from this multifrequency campaign.

The X-ray/TeV connection has been established in the past for many BL Lac objects (and for Mrk 501 in particular). However, the exact character of the correlation between the X-ray and TeV fluxes is known to vary from object to object, and from epoch to epoch in a given source, as widely discussed in, e.g., Krawczynski et al. (2004), Błazejowski et al. (2005), Gliozzi et al. (2006), and Fossati et al. (2008). Note that the data analyzed in those papers were obtained mostly during periods of high activity. Consequently, the conclusions presented were somewhat biased towards flaring activity, and hence they might not apply to the typical (average) behavior of the source, which is the main focus of this paper. Moreover, the data set from our campaign includes UV fluxes, soft-X-ray fluxes (down to 0.3 keV; Swift), and γ-ray fluxes spanning a very wide photon energy range (0.1 GeV−10 TeV; Fermi-LAT combined with MAGIC and VERITAS). This unique data set allows us to evaluate the multifrequency variability and correlations for Mrk 501 over an unprecedented range of photon energies.

Considering only the data set presented in this paper, we note that by steepening the high-energy electron continuum above the intrinsic break energy γ_br, 1 (and only slightly adjusting the other model parameters), one can effectively remove photons above 10 GeV in the SSC component, leaving a relatively steep spectrum below 10 GeV, similar to the one observed by Fermi-LAT during the time interval MJD 54862–54892 (see Section 4). Such a change should be accompanied by a decrease in the UV-to-soft-X-ray synchrotron fluxes by a factor of a few, but the data available during that time interval are not sufficient to detect this effect.¹⁶² This statement is further justified by the bottom panel in Figure 12, where the contributions of the different segments of electrons Comptonizing the different segments of the synchrotron bump to the average γ-ray emission of Mrk 501 are shown. Note that the lowest-energy electron population (γ < γ_br, 1) inverse-Compton upscattering only synchrotron photons emitted by the same population (ν < ν_br, 1 ∼ 10¹⁵ Hz; solid red line) may account for the bulk of the observed steep-spectrum γ-ray emission.

Another important conclusion from this figure is that Comptonization of the highest-energy synchrotron photons (ν>ν_br, 2 ∼ 6 × 10¹⁷ Hz) by electrons with arbitrary energies produces only a negligible contribution to the average γ-ray flux of Mrk 501 due to the Klein–Nishina suppression. Thus, the model presented here explains in a natural way the fact that the X-ray and TeV fluxes of TeV-emitting BL Lac objects are rarely correlated according to the simple scaling F_TeV ∝ F²_keV which would be expected from the class of SSC models in which the highest-energy electrons upscatter (in the Thomson regime) their own synchrotron photons to the TeV band (see, e.g., Gliozzi et al. 2006; Fossati et al. 2008). In addition, it opens a possibility for accommodating short-timescale variability (t_var < 4 days) at the highest synchrotron and inverse-Compton frequencies (hard X-rays and TeV photon energies, respectively). The reason for this is that, in the model considered here, these high-energy tails of the two spectral components are produced by the highest-energy electrons which are deep in the strong cooling regime (i.e., for which t'_rad ≪ R/c), and thus the corresponding flux changes may occur on timescales shorter than R/c δ (see in this context, e.g., Chiaberge & Ghisellini 1999; Kataoka et al. 2000).

A more in-depth analysis of the multifrequency data set (including correlation studies of the variability in different frequency ranges) will be presented in a forthcoming paper. The epoch of enhanced γ-ray activity of Mrk 501 (MJD 54952–54982; see Section 4) may be more difficult to explain in the framework of the one-zone SSC model, because a relatively flat Fermi-LAT spectrum above 10 GeV, together with an increased TeV flux as measured by the VERITAS and Whipple 10 m telescopes around this time, may not be easy to reproduce with a set of model parameters similar to that discussed in previous sections. This is mostly due to Klein–Nishina effects, which tend to steepen the high-energy tail of the SSC component, thus precluding the formation of a flat power law extending beyond the observed TeV energies. Hence, detailed modeling and data analysis will be needed to determine whether the enhanced VHE γ-ray activity period can be accommodated within a one-zone SSC model, or whether it will require a multi-zone approach.