Multimessenger Characterization of Markarian 501 during Historically Low X-Ray and γ-Ray Activity

MAGIC consists of two IACTs separated by a distance of 85 m. It is located at the Roque de los Muchachos Observatory, on the Canary Island of La Palma at an altitude of 2243 m above sea level. The telescopes work in an energy range between 50 GeV and tens of TeVs, with a sensitivity above 100 GeV (300 GeV) of about 2% (about 1%) of the Crab Nebula flux at low zenith angles (< 30°) after 25 hr of observations (see Figure 19 of Aleksić et al. 2016). With this performance, the MAGIC telescopes are very well suited to perform blazar observations in the VHE range.

During the 4 yr period from 2017 to 2020, Mrk 501 was observed by MAGIC for around 160 hr, yielding around 120 hr after the data selection based on the atmospheric transmission and night-sky background (NSB) levels. The data are analyzed using the MAGIC Analysis and Reconstruction Software package (Zanin et al. 2013; Aleksić et al. 2016). Owing to Mrk 501 being one of the brightest sources in the MAGIC source catalog, it is possible to observe it even during moon conditions. The analysis is adjusted accordingly to the high NSB levels as prescribed in Ahnen et al. (2017b).

The VHE flux light curve is computed with a minimum energy of 0.2 TeV in order to minimize the impact of the various observing conditions considered in this study, which can increase by a factor ∼2 the analysis energy threshold at zenith. Our selected energy threshold is compatible with the applied data selection (zenith<50° and exclusion of bright moon levels). We bin the observations once night-wise and once on a weekly basis. For bins with a significance of less than 2σ, the upper limits (ULs) with 95% confidence according to the Rolke method (Rolke et al. 2005) are computed. For the spectral reconstruction, we use a forward folding method assuming a simple power law as the spectral model to obtain the relevant parameters, which are summarized in Table 9, and the Tikhonov unfolding method to obtain the spectral data points (Albert et al. 2007). For each spectrum described in Section 4.2, we check if a log parabolic power-law model is preferred over a simple power-law model using a likelihood ratio test performed over the aggregated data in the corresponding time intervals. The log parabola describes slightly better the spectra, but the preference for this model is not statistically significant (< 3σ).

The LAT on board the Fermi Gamma-Ray Space Telescope is constantly monitoring the high-energy sky since its launch in 2008. As a pair conversion instrument, it is sensitive to an energy range from 20 MeV to beyond 300 GeV, and covers the whole sky every ∼3 hr (Atwood et al. 2009; Ackermann et al. 2012).

Owing to the low activity of Mrk 501, as well as the need to characterize its γ-ray emission on relatively short-time intervals, we decided to use the unbinned-likelihood tools provided by the FERMITOOLS software ¹⁰⁷ (v1.0.10), which is more suitable than the binned analysis for situations with low-event statistics. Table 1 summarizes the basic analysis settings that were used. The usage of 0.3 GeV as the minimum energy for the analysis (instead of the conventional 0.1 GeV) reduces the detected number of photons from the source. However, this reduction is small for hard-spectrum sources (photon index <2) like Mrk 501. On the other hand, the angular resolution (68% containment) of photons above 0.1 GeV is about 5°, while it is about 2° for photons above 0.3 GeV. This means that an LAT analysis above 0.3 GeV is less affected by the diffuse backgrounds (which are always softer than a photon index 2), and hence will lead to an increase in the signal-to-noise ratio for hard sources. Additionally, the LAT analysis above 0.3 GeV is less sensitive to possible contamination from nonaccounted (transient) neighboring sources, in comparison to an LAT analysis above 0.1 GeV. The maximum energy range is chosen to overlap with the MAGIC energy range when reconstructing spectra. The fourth Fermi-LAT source catalog (4FGL; Abdollahi et al. 2020) is used to build a first model consisting of all sources within the region of interest (ROI) plus 5°, which is a standard ROI to investigate. We fit the obtained model to our data set covering the time range from MJD 57,754 (2017 January 1, 00:00:00) to MJD 59,126 (2020 October 4, 00:00:00). The preliminary fit result is used to remove very weak components from the model (counts <1 or TS <3). Afterward, the data set is divided in 98 bins each of 14 days duration, and each bin is fit with the model. In the fitting procedure, only the normalization of bright sources (TS >10), sources close to the ROI center (<3°), and the diffuse background are allowed to vary. Additionally, the spectral parameters of Mrk 501 are allowed to vary. From this fit we obtain the flux values per time bin and produce the light curve for Mrk 501. To check the impact of very variable sources in our ROI, we recalculated the light curve freeing the normalization and spectral indices of all sources with a variability index >100. The resulting flux values agree with the previous light curve within the statistical uncertainties. Furthermore, we perform spectral analyses for the 2 yr period with very low activity, and for two-week bins centered around each of the three NuSTAR observations conducted during the campaign (see Section 2.3). We use the same ROI model and approach as for the light curve. We checked the likelihood ratio between a power-law and log parabolic power-law model applied to the time intervals. As for none of the spectra a preference of more than 3σ is seen for the log parabolic fit, a power law is chosen for Mrk 501 for the spectral reconstruction. The parameters obtained are summarized in Table 10. For the spectral data points the number of spectral bins is chosen according to the time intervals and flux level.

For part of the correlation analysis, a long-term Fermi-LAT light curve is used including all data since its launch. We therefore apply the same procedure as described above to the data from MJD 54,688 (2008 August 10, 00:00:00) to MJD 59,126 (2020 October 4, 00:00:00), also using 14 day bins. The starting MJD is the earliest possible in 2008 that is in line with the time bins used for the 4 yr epoch featured in this paper, and results in 317 time bins for the 2008–2020 period.

NuSTAR is one of NASA's Small Explorer satellites, sensitive in the hard-X-ray band. It has two multilayer-coated telescopes, focusing the reflected X-rays on the pixilated CdZnTe focal plane modules, FPMA and FPMB. The observatory provides a bandpass of 3–79 keV with spectral resolution of ∼1 keV. The field of view of each telescope is $\sim 13^{\prime}$, and the half-power diameter of an image of a point source is $\sim 1^{\prime}$. This allows for a reliable estimate and subtraction of instrumental and cosmic backgrounds, resulting in an unprecedented sensitivity for measuring hard-X-ray fluxes and spectra of celestial sources. For more details, see Harrison et al. (2013).

The data reported here (Table 12) cover three observations in 2017 and 2018 that were planned as part of two dedicated NuSTAR proposals (from PI Balokovic and PI Paneque). After screening for the South Atlantic Anomaly passages and Earth occultation, the observations yielded roughly 5 hr of on-target data per pointing. The raw data products are processed separately for each pointing with the NuSTAR Data Analysis Software (NuSTARDAS) package v.1.3.1 (via the script nupipeline), producing calibrated and cleaned event files. Source data are extracted from a region of 45'' radius, centered on the centroid of X-ray emission, while the background is extracted from a $1\buildrel{\,\prime}\over{.} \,5$ radius region roughly $5^{\prime}$ north of the source location. The spectra are binned in order to have at least 20 counts per rebinned channel. We consider the spectral channels corresponding nominally to the 3–30 keV energy range, where the source was detected. The mean net (background-subtracted) count rates in the modules FPMA and FPMB are consistent with each other.

We perform spectral fits to the data, using the standard NuSTAR response matrices and effective area files, via the NuSTAR data analysis package nuproducts. We adopt a simple power-law model, absorbed by intervening material with solar abundances and the Galactic column of 1.55 × 10²⁰ cm⁻² (Kalberla et al. 2005). In all cases, the spectra are adequately described by such a simple model, but in all three cases, we see a modest improvement to the spectral fit when we attempt a more complex model, such as log parabola. While this is not significant in any of the three pointings reported here, this spectral behavior is fully consistent with that reported in Furniss et al. (2015).

The study reported in this paper makes use of two instruments on board the Neil Gehrels Swift Gamma-Ray Burst Observatory (Gehrels et al. 2004); namely, XRT (Burrows et al. 2005) and the Ultraviolet/Optical Telescope (UVOT; Roming et al. 2005). The related observations were organized and performed within the framework of planned extensive multi-instrument campaigns on Mrk 501, which occur yearly since 2008 (Aleksić et al. 2015a). In this study we consider all observations within the years 2017 and 2020 for both the UVOT and XRT instruments with an extension using the all data from 2005 to 2020 for the long-term studies on the XRT light curve.

The Swift-XRT observations are performed in the Windowed Timing (WT) and Photon Counting (PC) readout modes, and the data are processed using the XRTDAS software package (v.3.5.0) developed by the ASI Space Science Data Center (SSDC), and released by the NASA High Energy Astrophysics Archive Research Center (HEASARC) in the HEASoft package (v.6.26.1). The calibration files from Swift-XRT CALDB (version 20190910) are used within the xrtpipeline to calibrate and clean the events. The X-ray spectrum from each observation is extracted from the summed cleaned event file. For WT readout mode data, events for the spectral analysis are selected within a circle of 20 pixel (∼46'') radius, which contains about 90% of the point-spread function (PSF), centered at the source position. For PC readout mode data, the source count rate is above ∼0.5 counts s⁻¹ and data are significantly affected by pileup in the inner part of the PSF. We remove pileup effects by excluding events within a 4–6 pixel radius circles centered on the source position and use an outer radius of 30 pixels. The background is estimated from a nearby circular region with a radius of 20 and 40 pixels for WT and PC data, respectively. The ancillary response files (ARFs) are generated with the xrtmkarf task applying corrections for PSF losses and CCD defects using the cumulative exposure map. The 0.3–10 keV source spectra are binned using the grppha task to ensure a minimum of 20 counts per bin, and then are modeled in XSPEC using power-law and log-parabola models (with a pivot energy fixed at 1 keV) that include photoelectric absorption due to a neutral-hydrogen column density fixed to the Galactic 21 cm value in the direction of Mrk 501, namely, 1.55 × 10²⁰ cm⁻² (Kalberla et al. 2005).

The Swift-UVOT data analysis reported here relates only to all the observations with the UV filters (namely, W1, M2, and W2) performed during the Swift pointings to Mrk 501, 259 exposures. Differently to the optical bands, the emission in the UV is not affected by the emission from the host galaxy, which is very low at these frequencies. We perform aperture photometry for all filters using the standard UVOT software within the HEAsoft package (v6.23) and the calibration included in the latest release of the CALDB (20201026). The source photometry is evaluated following the recipe in Poole et al. (2008), extracting source counts from a circular aperture of 5'' radius, and the background ones from an annular aperture of 26'' and 34'' for the inner and outer radii in all filters. The count rates are converted to fluxes using the standard zero-points (Breeveld et al. 2011) and finally dereddened considering an E(B − V) value of 0.017 (Schlegel et al. 1998; Schlafly & Finkbeiner 2011) for the UVOT filters effective wavelengths and the mean galactic interstellar extinction curve from Fitzpatrick (1999).

We focus on the R band for the optical wave band, as it is often done in previous studies of Mrk 501, and HSPs, in general. The optical data are collected within the GLAST and AGILE program (GASP) of the Whole Earth Blazar Telescope (WEBT; Villata et al. 2008, 2009; Gazeas 2016; Carnerero et al. 2017; Raiteri et al. 2017) including the instruments: West Mountain (91 cm), Vidojevica (140 cm), Vidojevica (60 cm), University of Athens Observatory (UOAO), Tijarafe (40 cm), Teide (STELLA-I), Teide (IAC80), St. Petersburg, Skinakas, San Pedro Martir (84 cm), Rozhen (200 cm), Rozhen (50/70 cm), Perkins (1.8 m), New Mexico Skies (T21), New Mexico Skies (T11), Lulin (SLT), Hans Haffner, Crimean (70 cm; ST-7; pol), Crimean (70 cm; ST-7), Crimean (70 cm; AP7), Connecticut (51 cm), Burke-Gaffney, Belogradchik, AstroCamp (T7), and Abastumani (70 cm). Additional data were provided by AAVSO and by the Tuorla observatory using the Kungliga Vetenskapsakademien (KVA) telescope.

The data analysis is performed using standard prescriptions. The host galaxy contribution is subtracted according to the Nilsson et al. (2007) recipe for an aperture of 75, which was adopted by the participating instruments. The R-band flux is then corrected for Galactic extinction assuming the values reported by Schlafly & Finkbeiner (2011). In order to account for instrumental (systematic) differences among the analyses related to the various telescopes (i.e., due to different filter spectral responses and analysis procedures, combined with the strong host galaxy contribution), offsets of a few mJy have to be applied. To calculate the corresponding offsets, KVA is used as a reference due to its good time coverage taking simultaneous data within two days into account. For data sets containing majorly data collected in 2020, when KVA was not operational anymore, Hans Haffner is used as the reference. The corresponding offsets can be found in Table 11. To further account for instrumental (systematic) uncertainties, a relative error of 2% is added in quadrature to the statistical uncertainties of all the flux values, as done in previous works (Ahnen et al. 2018). Afterward, the data sets from all the instruments are combined into a single R-band light curve and binned in 1 day time intervals.

We report here radio observations from the single-dish telescopes at the Owens Valley Radio Observatory (OVRO) operating at 15 GHz; the Medicina observatory operating at 8 GHz and 24 GHz; RATAN-600 at 4.7 GHz, 11.2 GHz, and 22 GHz; the Metsähovi Radio Observatory at 37 GHz; IRAM at 100 GHz and 230 GHz; the interferometry observations from the Very Long Baseline Array (VLBA) at 43 GHz; and the Submillimeter Array (SMA) at 230 GHz and 345 GHz.

For Metsähovi the detection limit of the telescope at 37 GHz is on the order of 0.2 Jy under optimal conditions. Data points with a signal-to-noise ratio below 4 are handled as nondetections. The flux density scale is set by observations of DR 21. Sources NGC 7027, 3C 274 and 3C 84 are used as secondary calibrators. A detailed description of the data reduction and analysis is given in Teraesranta et al. (1998). The error estimate in the flux density includes the contribution from the measurement root mean square and the uncertainty of the absolute calibration. The data from OVRO and Medicina were analyzed following the prescription from Richards et al. (2011) and Giroletti & Righini (2020), and provided by the instrument teams specifically for this study. The flux density measurements with the RATAN-600 radio telescope were obtained at three frequencies, 22 GHz, 11.2 GHz, and 4.7 GHz, over several minutes per object in transit mode (Parijskij 1993; Sotnikova 2020). The data reduction procedures and main parameters of the antenna and radiometers are described in, e.g., Udovitskiy et al. (2016) and Mingaliev et al. (2017). Flux density data of Mrk 501 at 230 GHz and 345 GHz were obtained with SMA as part of a monitoring program of mm band gain calibrators (Gurwell et al. 2007). Sources are periodically observed for 3–5 minutes, and calibrated against known standards (primarily Titan, Uranus, Neptune, or Callisto). The light-curve data are updated regularly at the SMA website. ¹⁰⁸

Additional radio data have been taken by VLBA. Since 2014 June, Mrk 501 has been observed monthly with the VLBA at 43 GHz within the VLBA-BU-BLAZAR program, which includes 38 γ-ray active galactic nuclei (AGNs). Fully calibrated data of Mrk 501 are posted at the program website. ¹⁰⁹ The data reduction is described in Jorstad et al. (2017). The total VLBA intensity individual measurements during this period are compatible with a constant emission with an average flux value of 0.34 ± 0.03 Jy, and the typical flux uncertainty is about 0.2 Jy. For the sake of clarity, these values are not reported in Figure 2.

For the single-dish instruments, Mrk 501 is a point source, and therefore the measurements represent an integration of the full source extension. The size of the radio emitting region is expected to be larger than that of the region of the jet that dominates the X-ray and γ-ray emission, known to vary on much shorter timescales than the radio emission. However, as reported by Ackermann et al. (2011), there is a correlation between the 8 and 15 GHz radio and the GeV emission of blazars. In the case of Mrk 501, Ahnen et al. (2017a) showed that the radio core emission increased during a period of high γ-ray activity. Therefore, a significant fraction of the radio emission seems to be related to the γ-ray component, at least during some periods of time, and hence should be considered when studying and interpreting the overall broadband emission of Mrk 501.

Linear polarization measurements have been taken in both the optical and radio band. Through the GASP-WEBT program, optical polarization data using the R band was obtained with the 70 cm telescope at the Crimean Observatory and with the 1.8 m telescope at the Perkins Observatory with an aperture of 75. Moreover, R-band polarimetry data were collected by the Nordic Optical Telescope (NOT) using the ALFOSC instrument ¹¹⁰ with an aperture of 15. Publicly available data from the Steward Observatory ¹¹¹ collected in the 5000–7000 Å band using an aperture of 3'' completes our optical polarization data set.

The GASP-WEBT optical polarization data are corrected for interstellar polarization (ISP) using field stars (numbers 6, 1, and 4 in Villata et al. 1998). The host galaxy contribution is taken into account assuming a seeing of 2'' and using Table B.1. in Nilsson et al. (2007). The degree of polarization is corrected for the host galaxy contribution according to the prescription given in Weaver et al. (2020). The analysis of the NOT data is done as described in Hovatta et al. (2016) and MAGIC Collaboration et al. (2018). Due to the small aperture (15) no ISP correction is applied, but for the host galaxy the correction is done in the same way as mentioned above. The same is the case for the Steward observatory data. As the Stokes parameters are ambiguous with respect to 180° shifts in the electric vector polarization angle (EVPA), the data are corrected using 180° shifts if differences between neighboring measurements are more than +90 (or less than −90) for observations within 15 days.

For polarization information in the radio regime, data obtained by the VLBA are used. A brief description of the polarization analysis can be found in Marscher & Jorstad (2021). The absolute EVPA calibration is performed using either quasi-simultaneous observations with the Very Long Array of the several objects from the sample twice per year, or the D-term method provided by Gomez et al. (2002). The degree of polarization and EVPA at 43 GHz are obtained from Stokes I, Q, and U parameters integrated over Stokes corresponding images at each epoch.

The fractional variability F_var is used, as given by Equation (10) in Vaughan et al. (2003), to estimate the degree of variability in each wave band. Its uncertainty is defined using the estimates and descriptions in Poutanen et al. (2008) and Aleksić et al. (2015a), resulting in Equation (2) in Aleksić et al. (2015a).

The fractional variability is computed for the light curves shown in Figure 2 for both the whole 4 yr period, as well as the 2 yr low-activity period. The results are displayed in Figure 3 for both the whole data set as well as only quasi-simultaneous data to the weekly binned MAGIC observations. Table 2 summarizes the results for the whole data set. Most light-curve measurements used for this study are retrieved from daily bins, except Fermi-LAT, with its 14 day binning, and MAGIC, with its 7 day binning. The 7 day binning is chosen for this analysis to mitigate the impact of the low flux levels, and hence the low number of significant single-night measurements (especially during the abovementioned 2 yr low-state period). However, when using a 1 day binning for the MAGIC fluxes, and discarding the nonsignificant fraction of flux measurements, the obtained fractional variability values are ${F}_{\mathrm{var}}^{2017\ \mathrm{to}\ 2020}\,=0.52\pm 0.03$ and ${F}_{\mathrm{var}}^{\mathrm{low}-\mathrm{state}}=0.25\pm 0.05$. These are very similar to values reported in Table 2 obtained with the 7 day binning. We note that the differences in the temporal bins used to characterize the variability may affect the comparability of F_var between different wave bands.

Zoom In Zoom Out Reset image size

For the two time intervals considered, the radio, optical, and UV frequencies show only mild flux variability. As reported in previous works (Furniss et al. 2015; Ahnen et al. 2018), Metsähovi shows a larger variability compared to the other radio data. In the case of the 2 yr low state, the fractional variability increases with the energy after the eV regime, but then it reaches a plateau for the GeV and TeV regimes. On the other hand, the 4 yr epoch shows a two peak structure with the 2–10 keV data showing a similar variability level as the VHE data.

Conventional methods to compute correlations between data sets are often hard to apply to astrophysical light curves due to the differences in sampling, varying uncertainties and the often long gaps between observations. The discrete correlation function (DCF) provides an assumption-free representation of the correlation without interpolating in time (Edelson & Krolik 1988).

To compute the DCF between two light curves, an appropriate time binning is chosen and applied to both light curves in the same manner, yielding pairs of flux data points. We choose the appropriate binning for each light curve by considering the typical decay times in the autocorrelation distributions of each wave band. Additionally, the compatibility of the chosen bins with the original binning of the different light curves is taken into account. The time bins used for the correlation analysis for all instruments can be found in Table 3. We then compute the DCF using Equation (4) in Edelson & Krolik (1988) if the number of flux pairs is bigger than 10. In order to identify delayed correlations, we shift the light curves with respect to each other for certain time intervals, and then repeat the correlation computations. For the UVOT data, we choose to select only the W1 filter for the correlation analyses because the light curves of all three UV filters show very similar behavior. An easy way of estimating the significance of the DCF values may be exploited by using the 1σ errors as defined in Equation (5) in Edelson & Krolik (1988). However, when the light curves include correlated red-noise data, which is the case for blazars, these 1σ may not be suitable to determine the significance of the correlation correctly. In these situations, the calculated significance could be overestimated, as discussed in Uttley et al. (2003). Therefore, in order to better assess the reliability of the significance in our correlations, we use dedicated Monte Carlo (MC) simulations that take the actual sampling and flux measurement uncertainties into account. This strategy, described further below, allows us to determine the 1σ and 3σ confidence levels that apply for the specific light curves probed, and hence provide a robust evaluation of whether a correlation is significant or not. For the purpose of comprehensibility and completeness, both the 1σ errors from Equation (5) in Edelson & Krolik (1988) and the dedicated MC-derived 1σ and 3σ contours are displayed in our results.

In order to estimate the statistical significance of the DCF measurements, we simulate 10,000 uncorrelated light curves for each wave band using the DELCgen package (Emmanoulopoulos et al. 2013; Connolly 2016). The underlying power density spectrum (PDS) of the original light curve is used to reproduce the level of variability in the simulations and the same time sampling as for the original data sets is applied. We choose a power law as the function to estimate the PDS after no significant preferences for more complicated functions are seen in our evaluations. As our flux distributions are mostly compatible with Gaussian distributions, we select the Timmer & Koenig method (Timmer & Koenig 1995) to simulate the light curves. Uncertainties on the simulated flux points are estimated by using the minimum error of the original light curve as a minimum error, and drawing random additional errors from a distribution constructed from the original relative error distributions. This for example results in a minimum error of 1.1 × 10⁻¹² ph cm⁻² s⁻¹ being added to all simulated MAGIC data points complemented by additional errors drawn from a distribution with a distribution peaked around 2.0 × 10⁻¹² ph cm⁻² s⁻¹. For the optical data, a minimum error of 0.06 mJy is added, with additional errors being drawn from a distribution peaking around 0.1 mJy with a tail toward higher values. For the Fermi-LAT light curve, the original relative error distribution shows a different distribution at higher than at lower flux level. Therefore the original relative error distribution is divided into two samples using the flux level of 1.3 × 10⁻⁸ ph cm⁻² s⁻¹ as a limit where the differences in error distributions appear. For each simulated flux point, the relative error is then drawn using the sample corresponding to its flux level as the weight. For all other wave bands no dependency on the flux level is found. Afterward, the same binning, time shifting, and correlation calculation methods as for the real light curves are applied to the simulated ones. The resulting DCF distributions allow us to derive confidence levels at which random correlation from originally uncorrelated light curves can be excluded. The computed significance contours derived in this manner consider the specific sampling and flux measurement errors of the two light curves being used, and hence are a reliable evaluation of the significance of our correlations. However, when the same computation is performed a large number, N, of times, as we do when computing the DCF for different time lags, one should also consider the "look-elsewhere effect." The effect could artificially increase the chance probability of obtaining, for a random time lag ΔT, a DCF value above the 3σ contour. One could potentially correct for this effect by considering the trial factors, as described in Gross & Vitells (2010) and Algeri et al. (2016). However, the "look-elsewhere effect" does not affect the DCF obtained for physically meaningful time lags, such as ΔT = 0, or general trends as, e.g., the broader peaks spanning over many consecutive time lags.

A 3σ correlation (DCF value >3σ confidence level) at zero time lag is found between the high-energy (HE) X-ray (2–10 keV) range, measured with Swift-XRT, and the VHE γ-ray (>0.2 TeV) range, measured with MAGIC, as shown in Figure 4(a). Additionally, a 3σ correlation in the form of a structure can be seen with XRT lagging behind MAGIC peaking at a time lag of 28 days. For the low-energy (LE) X-ray range (0.3–2 keV), the same two features can be seen in Figure 4(b). We note here that, while significant VHE-X-ray correlations have always been observed in the emission of Mrk 501 during flaring activities, such a correlation has always been elusive during periods of very low activity (e.g., Aleksić et al. 2015c; Ahnen et al. 2017a). The improved sensitivity to characterize the VHE emission, in comparison to past measurements, and the extensive data set collected during these 4 yr, made possible the measurement of a correlated behavior with a significance above 3σ during a period of historically low VHE and X-ray activity.

Zoom In Zoom Out Reset image size

Similarly, when comparing the two energy ranges in X-rays with each other, as shown in Figure 5, 3σ correlations are seen with a clear peak at time lag zero and a broader structure around time lag −35 days, suggesting that at least a fraction of the low-energy photons are lagging behind the more energetic ones.

Zoom In Zoom Out Reset image size

Between OVRO and MAGIC, a 3σ correlation is found for most positive time lags, as shown in Figure 19(c). As these correlations vanish when detrending the light curves as described further below in more detail, we can characterize them as long-term features. Other 3σ correlations can be seen for a few negative time lags. However, they occur simultaneous with an increase in the confidence levels (due to the smaller number of data pairs used in the computation), and hence we cannot say if they are physically meaningful, or due to potentially nonaccounted uncertainties because of the limitations in the time range used. We also observe a marginally significant (< ∼3σ) correlation when comparing the HE X-rays with OVRO where the DCF values rise for positive time lags with two values reaching the 3σ level at 42 days and 84 days, as shown in Figure 19(a) (see also Table 4). This behavior is not seen when using the LE band of the X-ray data set, as shown in Figure 19(b). These correlations disappear when detrending the light curves, and hence they can be ascribed to the long-term behavior in the light curves at HE X-rays and radio.

Other 3σ correlations are found between Fermi-LAT and OVRO with a time delay of 112 days by OVRO, as displayed in Figure 20(a). As for both OVRO and Fermi-LAT data have been taken continuously since 2008, we expanded the analysis to their long-term light curves displayed in Figure 1 to refine our picture of the correlation. A similar DCF distribution can be seen in Figure 6(a), but this time with smaller statistical uncertainties, due to the much larger data set spanning over a time window that is 3 times larger. The highest degrees of correlation occur at time lags of 126 and 238 days, but there is a 3σ correlation that persistently occurs for negative time lags larger than 3 months.

Zoom In Zoom Out Reset image size

We additionally took advantage of the multiyear long observations of Mrk 501 with Swift-XRT, and quantified the correlation between X-rays and γ-rays using the XRT and LAT data from the period 2008–2020. Interestingly, while the 2017–2020 data set does not show any correlation between X-ray and γ-ray fluxes, the 12 yr data set reveals a very clear correlation between these two bands, as shown in Figure 7. We note that the correlation exists for both XRT energy bands, 0.3–2 keV and 2–10 keV. The main difference between the 4 and 12 yr data sets is easily explained. In the longer data set, Mrk 501 did show several periods of extended high activity, and hence it becomes easier to measure a correlated behavior between the fluxes. We note that the highest DCF values occur at a time lag zero, but with a very broad and somewhat asymmetric bump with a 3σ correlation, that extends over a time lag of ±1 yr. This indicates that the measured correlation is dominated by the flux variations with timescales of about 1 yr, such as the substantial decrease in the activity during the period 2017–2020, in comparison with the activity during the period 2008–2016. A very similar and expected behavior can be seen for the long-term correlation between the LE X-rays and the HE X-rays (Figure 8(a)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As the correlations in our 12 yr data set are dominated by long-term (month to years) timescales, we additionally apply a detrending as described in Section 5.1. in Lindfors et al. (2016) to the light curves. We slightly adapt step 2 of the procedure by not using the variance of the original high-energy light curve to scale the variance of the polynomial fit to the low-energy data, but rather use the variance of a polynomial fit to the high-energy data for the scaling. In this way we are matching the variability of the long-term behaviors of both light curves. For our data set this is especially relevant, as the large flares also include very variable behavior on short timescales in the X-ray range (Figure 1), we would otherwise overcorrect the light curves leading to negative dips at these variable time epochs. The procedure allows us to subtract the long-term behavior from our light curves and provides us with a view on the shorter timescales and correlations.

When the long-term behavior is removed from the light curves, one sees that the correlation between the γ-ray and radio light curves disappears, as shown in Figure 6(b). However, this is not the case for the correlation between the γ-ray and X-ray fluxes, as shown in Figures 7(c) and (d), where prominent bumps centered at time lag zero appear. When considering the correlation between the two X-ray energy bands, namely, 0.3–2 keV and 2–10 keV (see Figure 8(b)), the detrended light curves yield a DCF plot with a clear correlation peak centered at time lag zero, but also some indications of peaks that repeat with a period of about 30 days. The potential periodicity in the X-ray emission of Mrk 501 is discussed further in Section 3.3.

For the rest of our 4 yr data set, we see a clear correlation between the UV and the R band (see Figure 20(b)). The highest DCF values occur at a time lag zero, but the peak extends over ±1 week, due to the relatively slow flux variations in these two bands. The correlation analysis between the R-band and the HE γ-ray emission measured with Fermi-LAT does show a small DCF peak at about time lag of zero, but the DCF values do not reach the 3σ significance level (see Figure 20(c)).

We did not find a 3σ correlation for any other combination of untreated (not detrended) light curves, except for a marginally significant correlation with a delay of −28 days between the two long-term radio data sets from OVRO and Metsähovi (see Figure 20(d)). However, if we apply the detrending described above to our 4 yr data set, we can identify a 3σ correlation between the VHE and HE γ-rays as depicted in Figure 21.

Table 4 reports the DCF values and the 3σ contour limits for the various energy band combinations and the relevant time lags that were discussed above in the text.

As for the polarization data, we do not find any significant correlation between the polarization degree and the flux levels, neither for optical nor for the 43 GHz radio (see Figure 22). Concerning the EVPA, in the optical R band and the radio frequency probed, 43 GHz, the preferred values fluctuate around 130°, and they are independent (not correlated) of the polarization degree (see Figures 23(a) and (b)). Due to 180° ambiguity and the corrections applied, some of the EVPA values are above 180° for the optical data. Here the Von Mises distribution can be used to check the accurate description of the data, which results in a mean value of 1337 ± 171. This matches the EVPA angle of 134° ± 5 ° reported from the first polarization measurements in the X-rays of Mrk 501 (Liodakis et al. 2022) as well as the measured jet direction of 1197 ± 118 from Weaver et al. (2022). When comparing the polarization properties in the optical and radio regime with each other, no correlation can be seen between the polarization degree and angle (see Figure 24).

Both the autocorrelation of the LE and HE X-ray ranges as well as their correlation with each other (Figure 8) indicate the possibility of a periodic behavior of ∼30 days. Previously for the large 1997 flare a periodicity both in the VHE and the X-ray data was claimed at a timescale of 23 days (Osone 2006; Kranich 1999) and explained with a binary black hole system as the center of Mrk 501 (Rieger & Mannheim 2000). Additionally, a periodicity of 332 days was reported in Bhatta (2019) for the Fermi-LAT data set with a significance of 99.4%, slightly below the 3σ level, as well as at 195 days with a significance slightly above 90%.

Exploiting the availability of our 12 yr data set, we apply the Lomb–Scargle periodogram (LSP; Lomb 1976; Scargle 1982). Details on the implementation and computation can be found in Appendix D as well as on the significance estimation which is based on the same simulations as described in Section 3.2.

No significant (> 3σ) periodicity can be detected for the two X-ray energy ranges (see Figure 25). All peaks in the signal LSP including the one at f ∼ 1/30 days⁻¹ are accompanied by corresponding rises in the confidences limits and the signal does not go above the 3σ confidence limits. The peaks can be attributed to the time sampling of the X-ray light curve, which is coordinated together with Earth-based telescopes (see details in Appendix D). Similarly, for the VHE light curve no periodicity can be identified.

Furthermore, no significant periodicity is found for the Fermi-LAT (Figure 27(a)) or OVRO (Figure 27(b)) data. However, we can reproduce the reported weak hint for a "330 day" periodicity in the Fermi-LAT light curve at a significance of 2.3σ as well as at ∼200 days at 2.2σ.

O'Neill et al. (2022) laid out a standard in which they proposed a 3σ global significance cutoff, which clearly Mrk 501 does not meet. For this reason, we interpret all claimed periodicities in Mrk 501 as red noise for now, although we will continue to monitor the source for periodicities that reach our cutoff threshold.

The 2 yr long period without substantial flux variations allows one to average a large amount of data to compute a broadband SED with small statistical uncertainties, despite the historically low activity of the source. This is particularly important at γ-ray energies, where the sensitivity of the instruments is more limited, in comparison with optical or X-ray instruments, and both Fermi-LAT and MAGIC would have limitations to deliver accurate spectra for timescales of days or even weeks (for such low source activity). Therefore, the entire 2 yr data set from 2017 June 17 to 2019 July 23 (MJD 57,921 to MJD 58,687) was used to derive the Fermi-LAT and MAGIC spectra shown in Figure 9(a). For the other wave bands, we determine the average of the spectral points weighted with their uncertainties from the individual observations, which is shown together with the minimum and maximum spectral flux values within the 2 yr long epoch considered. Owing to the low flux variability in the radio, optical, and X-ray bands, the weighted average spectral points are used as a good proxy of the emission of Mrk 501 during this 2 yr long low state. The hard-X-ray spectrum derived with NuSTAR data from 2018 April 20 (MJD 58,228) is also depicted, showing a good agreement (within 30%) with the weighted average Swift-XRT spectrum in the overlapping energy range. This agreement is consistent with the relatively low flux variability at X-rays mentioned above and suggests that the NuSTAR spectrum from 2018 April 20 is a good approximation of the hard-X-ray emission of Mrk 501 during this 2 yr long epoch. As a consistency check, Figure 28 shows the broadband SED during the abovementioned 2 yr long period of low activity, together with the SED around the NuSTAR observation from 2018 April 20, using MAGIC and Fermi-LAT spectra derived with data within ±1 week of the NuSTAR observation. In radio, optical and soft X-ray data the sensitivity is good enough to use the individual observations (typically less than 1 hr long). The good agreement among the spectra further supports the use of the NuSTAR observation from 2018 April 20 as a good proxy of the hard-X-ray emission of the 2 yr long integrated (weighted average) spectrum of Mrk 501.

Zoom In Zoom Out Reset image size

For comparison purposes, Figure 9(a) also depicts the typical (nonflaring) state of Mrk 501 from Abdo et al. (2011a). This helps visualize the historically low-activity state of Mrk 501 during the 2 yr time interval that goes from 2017 June 17 to 2019 July 23 (MJD 57,921 to MJD 58,687).

For the first NuSTAR observation (NuSTAR-1) on 2017 April 28 (MJD 57,871) no simultaneous VHE observations are available. Therefore we choose to use all nights inside ±1 week for the computation of both the MAGIC as well as the Fermi-LAT spectra. For the radio, optical and X-ray frequencies, we use simultaneous data (within ±5 hr of the NuSTAR observation) for all instruments. The resulting SED is shown in Figure 9(b). As the Fermi-LAT spectral analysis yields only one flux point, the figure additionally depicts the spectral shape (with statistical uncertainties) derived with our analysis.

The second NuSTAR observation (NuSTAR-2) on 2017 May 25 (MJD 57,898) was performed simultaneous to the observations carried out with the MAGIC telescopes, which allows us to derive an SED with simultaneous (within ±4 hr) data from radio to VHE γ-rays. The only exception is the spectrum from Fermi-LAT, which in order to reduce the statistical uncertainties (and owing to the lack of significant flux variability) is derived with data integrated within ±1 week of NuSTAR-2. The corresponding multi-instrument SED is shown in Figure 9(c).

Additionally, we added infrared (IR) points taken by the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE) mission ¹¹² (Mainzer et al. 2014). For the low state all simultaneous data were averaged as done for the other low-energy wave bands. No simultaneous data were available for the two NuSTAR time intervals. We chose the closest NEOWISE observation, which took place on MJD 57,807. However, NEOWISE only makes use of two (W1, W2) of the four filters of the Wide-field Infrared Survey Explorer (WISE) spacecraft. In order to include also information about the other two filters, we added archival WISE data from 2010 (Wright et al. 2010; Cutri et al. 2021) to all three SEDs from Figure 9 and used them for our spectral studies based on the observation of very low variability in the IR band.

Thereafter, all three SEDs (low state, NuSTAR-1, and NuSTAR-2) are characterized within different theoretical scenarios, as described in the following sections. Owing to the very low variability at radio and optical frequencies and the consideration that this low-energy emission is dominated by the contributions from different, more extended and outer regions of the jet (see, e.g., Acciari et al. 2020a), the radio and optical flux points are treated as upper limits in the adjustment of the theoretical models to the data. For all theoretical models, we consider a γ-ray absorption according to the extragalactic background light (EBL) model of Franceschini et al. (2008), which for the distance of Mrk 501 and the energies considered here is perfectly compatible with that of many other EBL models (see, e.g., Domínguez et al. 2011).

It is worth noting that, even for accurately determined SEDs such as the ones presented here, there is an ample degeneracy in the model parameters from the various theoretical scenarios. Therefore, the following results should not be interpreted as unique solutions, but rather as plausible theoretical scenarios that are able to explain the broadband data in line with our physical understanding of the underlying mechanisms.

To exploit the historically low-activity broadband SED of Mrk 501 and explore different blazar scenarios to describe this sort of "baseline broadband emission," we employ various widely used theoretical frameworks that consider leptonic, hadronic, and lepto-hadronic scenarios. For all models the low-energy component is mainly produced by synchrotron radiation of relativistic electrons. As leptonic scenarios we consider models where the high-energy component is also produced purely by the emission of relativistic electrons. For the scenarios where relativistic protons are contributing to the emission, we are following the conventions defined in Cerruti (2020). We define hadronic models as the ones where purely hadronic-initiated emission processes, i.e., proton synchrotron, dominate the high-energy component. For the lepto-hadronic models, synchrotron self-Compton (SSC) processes are responsible for a significant part of the high-energy component but proton-initiated processes are present as well leading to an expected emission of neutrinos while staying in a similar parameter space as for the leptonic models.

4.2.1. Leptonic Model

To evaluate a leptonic origin of the low-state SED, we choose a stationary one-zone SSC model (see, e.g., Ghisellini & Maraschi 1996; Tavecchio et al. 1998). For the model fitting, we compare two independent frameworks.

Within the first framework, in order to constrain the model parameters more efficiently, the naima package (Zabalza 2015), is modified to derive the best-fit and uncertainty distributions of spectral model parameters through Markov Chain Monte Carlo (MCMC) sampling of their likelihood distributions. Our prior constraints of the model parameter space obtained via "fit by eye" strategy, and the data likelihood function is passed on to emcee (Foreman-Mackey et al. 2013), which is a python implementation of the Goodman & Weares Affine Invariant MCMC Ensemble sampler (Goodman & Weare 2010).

As a second framework, we use the public open source C/Python framework jetset version 1.2.2 ¹¹³ (Tramacere et al. 2009, 2011; Tramacere 2020). We first approximate the spectral shape and use it together with basic information about the electron distribution and jet properties as input for a pre-fit to constrain the parameter space. Afterward, a full spectral fit is carried out using the minuit minimizer. The resulting best fit is then used as a prior for an MCMC chain using as recommended by the instructions on the algorithm, 128 walkers, 10 steps of burn-in, and 50 run steps, varying all free parameters from the minuit fit. This enables us to both improve the best fit as well as obtain a confidence interval on the resulting model.

We use the same reference parameters for both frameworks. Our emission zone is assumed to be a spherical blob with radius $R^{\prime}$ filled with relativistic electrons. The radiating electron distribution is described by a broken power law where the high-energy power-law index α₂ is connected to the low-energy index α₁ via the relation α₂ = α₁ + 1. For all models the parameters referred to in the blob frame are marked as primed while unmarked ones are in the observer frame.

As no short-timescale variability is observed in our data, a value of $R^{\prime} =1.14\times {10}^{17}$ cm is chosen for the emission region size. This is consistent with the values used in a previous work reporting the typical (nonflaring) broadband emission of Mrk 501 (Abdo et al. 2011a) based on the observed low variability taking our choice for the Doppler factor into account. We assume a viewing angle of ≈1/Γ_b , with Γ_b being the bulk Lorentz factor. Therefore, Γ_b and the Doppler factor δ can be considered equal. For the Doppler factor, we fixed δ = 11, which showed a good agreement between the models and our data in our preliminary fits where we tried a grid of values around δ = 12 from Abdo et al. (2011a). For this, we first used a step size of 2 and then optimized further using a step size of 1 for the region between 10 and 12. The literature shows that, because of short-time flux variability and VHE γ-ray opacity arguments, values of δ > 20 are often necessary to explain adequately the data from BL Lac objects (Tavecchio et al. 1998). On the other hand, high Doppler factors imply a small angle between our line of sight and the blazar, which cause tensions with the idea of a parent population of inefficiently accreting AGN including both BL Lac objects and Fanaroff–Riley Type I (FR-I) radio galaxies (Chiaberge et al. 2000; Tavecchio 2006). The Doppler factor δ = 11 used in this study is in good agreement with blazar unification schemes, as well as being consistent with our low flux variability levels and the required transparency for the measured VHE γ-rays. However, it is still considerably higher than the Doppler factors derived in the radio regime (Finke 2019) and in strong tension with the slow jet component speeds observed in the radio (Giroletti et al. 2004; Piner et al. 2010). This supports the scenario where part of radio emission is assumed to be produced in a different, more extended region of the jet than the X-ray to TeV emission.

For the magnetic field $B^{\prime}$, we obtain $B^{\prime} \sim$ 0.025 G in our first fits, and fix it at this value. For both frameworks we fix the break energy ${\gamma }_{\mathrm{br}}^{{\prime} }$ of the electron distribution according to the cooling break defined by Equation (30) in Tavecchio et al. (1998) with the cooling time being equal to the escape time of the radiating particles ${t}_{\mathrm{cool}}^{\prime} ={t}_{\mathrm{esc}}^{\prime} =\nu \tfrac{R^{\prime} }{c}$ with our choice of ν = 1. Longer escape times with ν = 2, 3, ... are possible as well, but prior investigations do not show any major changes for our physical conclusions. A ν of two would, for example, still allow the ambient parameters to stay at the same level as for the results presented below and just the electron parameters would change to slightly lower values for the energies and the spectral indices. Additionally, we fixed the minimum energy of the electrons to ${\gamma }_{\min }^{{\prime} }=1000$ to further constrain the fitting procedure and following the discussion in Aleksić et al. (2015a). We note that our ability to constrain the model parameter ${\gamma }_{\min }^{{\prime} }$ is strongly limited by the lack of data at MeV energies, as well as the contributions from other (more extended) regions to the emission at radio and optical frequencies.

The resulting low-state SED data–model adjustment with the theoretical frameworks mentioned above is shown in Figure 10, and the corresponding model and energy parameters are reported in Table 5. The two theoretical frameworks used yield very similar results.

Zoom In Zoom Out Reset image size

Table 5. Parameter Values from the Leptonic One-zone SSC Models used to Describe the Low-state SED of Mrk 501 Derived with Data from 2017 June 17 to 2019 July 23 (MJD 57,921 to MJD 58,687), as Described in Section 4.2.1 and Shown in Figure 10

	Le	α1	${\gamma }_{\mathrm{br}}^{{\prime} }$	${\gamma }_{\max }^{{\prime} }$	${U}_{{\rm{e}}}^{{\prime} }$	${U}_{{\rm{B}}}^{{\prime} }$	${U}_{{\rm{e}}}^{{\prime} }$/${U}_{{\rm{B}}}^{{\prime} }$	Ljet
	(erg s−1)				(erg cm−3)	(erg cm−3)		(erg s−1)
Modified Naima	7.7 × 1043	2.6	2.0 × 105 a	1.2 × 106	5.2 × 10−4	2.5 × 10−5	21	8.8 × 1043
Jetset	8.4 × 1043	2.6	2.0 × 105 a	1.2 × 106	5.7 × 10−4	2.5 × 10−5	23	9.6 × 1043

Note. In the two model realizations, the radius of the emission region $R^{\prime}$ is fixed to 1.14 × 10¹⁷ cm, the Doppler factor δ to 11, and the magnetic field $B^{\prime}$ to 0.025 G. For the radiating electron distribution, a broken power law is used with the spectral indices fulfilling α₂ = α₁ + 1. The minimum energy of the electrons ${\gamma }_{\min }^{{\prime} }$ is set to 1000. The table reports, for the broken power law describing the shape of the electron distribution, the first spectral index α₁, the break energy ${\gamma }_{\mathrm{br}}^{\prime}$ and the maximum energy ${\gamma }_{\max }^{{\prime} }$ as well as the electron luminosity L_e, the energy densities held by the electron population ${U}_{{\rm{e}}}^{{\prime} }$ and the magnetic field ${U}_{{\rm{B}}}^{{\prime} }$, their ratio, and the total jet luminosity L_jet. For the EBL γ-ray absorption at a redshift z = 0.034, the model from Franceschini et al. (2008) is used.

^a Fixed to the cooling break.

Download table as:  ASCIITypeset image

4.2.2. Hadronic Model

From a purely electromagnetic perspective, leptonic radiative models and hadronic radiative models can fit blazar SEDs equally well (see, e.g., Cerruti 2020). While for luminous blazars hadronic models face the difficulty of requiring often super-Eddington powers, and can thus be disfavored from an energetic point of view, this is not true for HSP BL Lac type objects, and in particular when low flux states are studied. In this case, leptonic and hadronic emission scenarios are truly degenerate in their photon emission and can only be distinguished via the detection of neutrinos, naturally produced in hadronic interactions while absent in leptonic ones. With this in mind, we investigate here a hadronic modeling of the SED of Mrk 501 in its historically low-activity state described above (see Figure 9(a)).

In order to apply this modeling, we again employ two different numerical frameworks. The first numerical code, the LeHa code, used for the hadronic modeling is described in Cerruti et al. (2015). The code simulates photon and neutrino emission from a spherical plasmoid (with radius $R^{\prime}$) in the jet, moving with Doppler factor δ and filled with a homogeneous and entangled magnetic field $B^{\prime}$. The plasmoid contains a population of primary electrons and protons. Proton–photon interactions via photomeson and Bethe–Heitler channels (of protons on synchrotron photons produced by primary electrons) inject secondary particles in the emitting region; these secondary particles trigger synchrotron-supported pair cascades for which the synchrotron emission at equilibrium is computed.

The second code is SOPRANO ¹¹⁴ (see, e.g., Gasparyan et al. 2022), a time-dependent code including leptonic and hadronic processes designed to study the particle interaction mechanisms in different astrophysical objects. The code follows the temporal evolution of the isotropic distribution functions of primary injected particles and the secondaries produced in photopair and photopion interactions, alongside the evolution of photon and electron/positron distribution functions. In this case, the final spectrum is computed by evolving the kinetic equations for several dynamical timescales to guarantee that the steady-state condition is achieved.

It is important to note that when we are comparing the model parameter values between the different frameworks, the first gives the parameters of the steady-state solution of the particle distributions while the second gives the particle properties at injection. Nonetheless, the results are compatible between the different codes.

We first make the same assumptions as for the leptonic case and choose the radius $R^{\prime}$ to be 1.14 × 10¹⁷ cm and δ = 11. Simple power-law distributions are chosen for the radiating electron and proton distributions. The maximum proton energy is determined by equating the acceleration timescale (see, e.g., Rieger et al. 2007) and the escape timescale (${t}_{\mathrm{esc}}^{\prime} \sim R^{\prime} /c$). The minimum proton energy is not constrained by the data, we therefore fix it to ${\gamma }_{\min ,{\rm{p}}}^{{\prime} }=1$ following the discussion in Cerruti et al. (2015). To be in line with a slow variability, the magnetic field is decreased as much as possible while looking for a fitting scenario (a higher magnetic field being associated with a faster synchrotron cooling). The resulting models from the two frameworks are shown in Figure 11. For simplicity, all individual components of the emission are only shown for the LeHa framework while for the SOPRANO result only the total emission in photons and neutrinos is indicated. A more detailed comparison including the distinguishing components of our hadronic and lepto-hadronic scenarios is shown in Figure 29. The parameter values are reported in Table 6 for the LeHa framework stating the steady-state particle properties as is done for the leptonic case, while Table 14 states the parameters of the SOPRANO framework displaying the injected particle properties. For the injected spectrum a simple power law is used to describe the distribution while it evolves into a broken power law for the steady-state spectrum. For some cases, the break energy ${\gamma }_{\mathrm{br}}^{{\prime} }$ obtained in the optimization of the model is higher than the maximum particle energy ${\gamma }_{\max }^{{\prime} }$ constrained by the acceleration timescale. For these cases the results are shown for a simple power law. In these hadronic scenarios the γ-ray emission is ascribed to proton-synchrotron radiation. Emission from pair cascades (triggered by pion decay and Bethe–Heitler pair production) and from muons (synchrotron radiation) is subdominant and largely absorbed by the EBL. The total jet powers required to produce the observed photon flux are 3.1 × 10⁴⁶ erg s⁻¹ and 3.0 × 10⁴⁶ erg s⁻¹ for the two frameworks, which is sub-Eddington for typical supermassive black hole masses.

Zoom In Zoom Out Reset image size

Table 6. Parameter Values Obtained by the LeHa Code for the Hadronic and Lepto-hadronic One-zone Models Used to Describe the Low-state SED of Mrk 501 Derived with Data from 2017 June 17 to 2019 July 23 (MJD 57,921 to MJD 58,687), as Described in Section 4.2.2 and Shown in Figure 11 (hadronic), and in Section 4.2.3 and Figure 12 (lepto-hadronic)

	Hadronic	Lepto-hadronic
$B^{\prime}$ [G]	3	0.025
$R^{\prime}$ [cm]	1.14 × 1017	1.14 × 1017
${{N}}_{\mathrm{total}}\,^{\prime}$ [1/cm3]	3	2.3 × 104
${n}_{{\rm{e}}}^{{\prime} }$/${n}_{{\rm{p}}}\,^{\prime}$	2.2	0.01
αe,1	2.5	2.6
αe,2	⋯	3.6
${\gamma }_{\min ,{\rm{e}}}\,^{\prime}$	400	1 × 103
${\gamma }_{\mathrm{br},{\rm{e}}}\,^{\prime}$	⋯	2.0 × 105
${\gamma }_{\max ,{\rm{e}}}^{{\prime} }$	3.5 × 104	1.2 × 106
αp	2.2	2.0
${\gamma }_{\min ,\ {\rm{p}}}\,^{\prime}$	1	1
${\gamma }_{\max ,\ {\rm{p}}}\,^{\prime}$	1.1 × 1010	2 × 107
${U}_{{\rm{e}}}^{{\prime} }$ [erg cm−3]	2.2 × 10−7	5.3 × 10−4
${U}_{{\rm{B}}}^{{\prime} }$ [erg cm−3]	0.36	2.5 × 10−5
${U}_{{\rm{p}}}\,^{\prime}$ [erg cm−3]	0.05	5.6
${U}_{{\rm{e}}}^{{\prime} }$/${U}_{{\rm{B}}}^{{\prime} }$	6.1 × 10−7	21.2
${U}_{{\rm{p}}}\,^{\prime}$/${U}_{{\rm{B}}}^{{\prime} }$	0.14	2.3 × 105
Ljet [erg s−1]	3.1 × 1046	1.7 × 1048

Note. In these model realizations, the radius of the emission region $R^{\prime}$ is fixed to 1.14 × 10¹⁷ cm and the Doppler factor to δ = 11. For the radiating particle distributions simple and broken power laws are used. The table reports the magnetic field $B^{\prime}$ and the radius $R^{\prime}$ of the emitting region and, for the assumed power-law distribution for the electrons and protons, the density ratio ${n}_{{\rm{e}}}^{{\prime} }/{n}_{{\rm{p}}}\,^{\prime}$; the total density of radiating particle N_total'; the slopes α_e,1, α_e,2, α_p; and the minimum, maximum, and break energies ${\gamma }_{\min ,{\rm{e}}}\,^{\prime}$, ${\gamma }_{\min ,{\rm{p}}}^{{\prime} }$, ${\gamma }_{\max ,{\rm{e}}}^{{\prime} }$, ${\gamma }_{\max ,{\rm{p}}}^{{\prime} }$, ${\gamma }_{\mathrm{br},{\rm{e}}}\,^{\prime}$. We also show the energy densities held by the electron population ${U}_{{\rm{e}}}^{{\prime} }$, the proton population ${U}_{{\rm{p}}}\,^{\prime}$, and the magnetic field ${U}_{{\rm{B}}}^{{\prime} }$, their ratios and the total jet luminosity L_jet. For the EBL γ-ray absorption at a redshift z = 0.034, the model from Franceschini et al. (2008) is used.

Download table as:  ASCIITypeset image

The hadronic scenario enables the prediction of an estimated neutrino spectrum, which can then be converted into an IceCube detection rate per year using the instrument effective area by IceCube Collaboration (2021). For our model, the expected neutrino rate is 1.1 × 10⁻⁵ events per year for the LeHa model and 1.5 × 10⁻⁵ events per year for the SOPRANO model. As this result does not come from a fit, and is not a unique solution, hadronic models exploiting other parts of the parameter space could lead to different explanations for the SED, and different neutrino spectra can be produced. A brighter neutrino emission can be achieved if the emitting region is more compact (smaller size and higher particle density), resulting in a higher proton–photon interaction rate, but this is limited by our choice of $R^{\prime}$ based on the observed low variability.

4.2.3. Lepto-hadronic Model

In recent years special attention has been given to mixed lepto-hadronic models, in which the high-energy SED component is associated to a combination of both leptonic (inverse Compton) and hadronic (emission by cascades triggered by hadronic interactions) processes. It is of particular interest due to the fact that the first evidence for the detection of joint photon and neutrino emission from the direction of a blazar, the 2017 flare of TXS 0506 + 056 (IceCube Collaboration et al. 2018), supported this kind of emission scenario, disfavoring a proton-synchrotron one (see, e.g., Keivani et al. 2018; Gao et al. 2019; Cerruti et al. 2019). In the lepto-hadronic solutions considered for this work, the bulk of the high-energy SED component is due to SSC, while the hadronic components are subdominant and can emerge (and dominate the SED) in hard X-rays, filling in the SED dip, and in the VHE band. In this scenario, the proton-synchrotron emission is very suppressed, mainly due to the lower magnetization of the emitting region with respect to proton-synchrotron solutions. We explore this lepto-hadronic model starting from the SSC solution described in Section 4.2.1, and adding a proton distribution with index α_p = 2.0 and ${\gamma }_{\max ,{\rm{p}}}^{\prime} =2\times {10}^{7}$. Again both the LeHa code as well as the SOPRANO code as described in Section 4.2.2 are utilized. The models are shown in Figure 12 with the more detailed comparison given in Figure 30. The parameter values are reported in Table 6 for the LeHa framework stating the steady-state particle properties as is done for the leptonic case, while Table 14 states the parameters from the SOPRANO framework displaying the injected particle properties. Compared to the proton-synchrotron solution, this scenario is much more demanding in terms of energetics, with a total jet power of 1.7 × 10⁴⁸ erg s⁻¹ and 4.4 × 10⁴⁷ erg s⁻¹ for the two frameworks, which are already super-Eddington for a black hole mass of 10⁹ M_⊙ (L_Edd ∼ 10⁴⁷ erg s⁻¹). However, this is not a strong constraint on the model: the solution shown here is not the result of the fit and not achieved by minimizing the jet power; the proton distribution is also conservative with respect to the total power, and harder distributions, or the introduction of a low-energy cutoff (${\gamma }_{\min ,{\rm{p}}}^{\prime} \gt 1$) would significantly lower the jet power. As an example, a similar electromagnetic emission can be achieved by hardening the proton distribution to α_p = 1.8 and lowering the proton power to 6.6 × 10⁴⁷ erg s⁻¹. The expected IceCube neutrino rates for the solutions shown in Figure 12 are 5.8 × 10⁻³ events per year for the LeHa model and 6.3 × 10⁻³ events per year for the SOPRANO model. This value is larger than the proton-synchrotron solution, closer to the IceCube detection energy range. It is a rate that remains however low (less than one event in two decades) and consistent with the nondetection of Mrk 501 as a point-like source in the IceCube data.

Zoom In Zoom Out Reset image size

The two NuSTAR observations performed in 2017 April and May, right before the 2 yr long period of historically low activity, provide us with the opportunity to characterize the overall decrease in the broadband emission of Mrk 501 down to this sort of "baseline emission." The variable behavior of Mrk 501 is often ascribed to variations in a radiatively efficient electron population (see, e.g., Ahnen et al. 2017a, 2018; Acciari et al. 2020a), and hence this temporal evolution study concentrates on SSC scenarios. In these theoretical frameworks, hard X-rays contain information from the dynamics of the highest-energy electrons. Because of the very low activity of Mrk 501 during the time interval considered here, hard-X-ray instruments such as Swift-BAT or the International Gamma-Ray Astrophysics Laboratory do not have sufficient sensitivity to detect the X-ray emission, and thus only NuSTAR has the capability to provide an accurate characterization of the hard-X-ray emission of Mrk 501.

First, we study the temporal evolution of the broadband SED using the same one-zone SSC scenario employed to describe the historically low-activity state reported in Section 4.2.1. We take the model parameters reported in Table 5 as a starting point and evaluate what parameters need to vary to explain the data from the previous months. As a second approach, we consider the 2 yr long low-state SED reported in Section 4.2.1 as sort of steady (or very slowly variable) "baseline broadband emission," and evaluate the presence of an additional region (located somewhere else along the jet of Mrk 501) whose emission is variable on timescales of weeks and months, and is responsible for the blazar activity during the NuSTAR-1 and NuSTAR-2 observations.

4.3.1. One-zone Model

As a starting point, we adopt the SSC scenario shown in Figure 10, with the model parameters reported in Table 5. Subsequently, we change the model parameter values to describe the SEDs from the months before the 2 yr long low state. The environmental parameters, including the magnetic field, are kept as close as possible to the values from the low-state model, while most of the model parameters describing the electron distribution are allowed to vary. The model parameters that describe well the SEDs during the NuSTAR-1 and NuSTAR-2 observations are reported in Table 7, and the model curves (together with the SED data) are displayed in Figure 13.

Zoom In Zoom Out Reset image size

Table 7. Parameter Values from the Leptonic One-zone SSC Models Used to Describe the Broadband SED of Mrk 501 During the NuSTAR-1 and NuSTAR-2 Observations on 2017 April 28 (MJD 57,898) and 2017 May 25 (MJD 57,871), as Shown in Figure 13

a) Model for NuSTAR-1 for a Magnetic Field of $B^{\prime} =$ 0.01 G
	Le	α1	${\gamma }_{\mathrm{br}}^{{\prime} }$	${\gamma }_{\max }^{{\prime} }$	${U}_{{\rm{e}}}^{{\prime} }$	${U}_{{\rm{B}}}^{{\prime} }$	${U}_{{\rm{e}}}^{{\prime} }$/${U}_{{\rm{B}}}^{{\prime} }$	Ljet
	(erg s−1)				(erg cm−3)	(erg cm−3)		(erg s−1)
Modified Naima	1.1 × 1044	2.3	6.6 × 105 a	7.2 × 106	7.1 × 10−4	4.0 × 10−6	178	1.1 × 1044
Jetset	1.1 × 1044	2.3	6.6 × 105 a	7.3 × 106	7.1 × 10−4	4.0 × 10−6	178	1.1 × 1044
b) Model for NuSTAR-2 for a Magnetic Field of $B^{\prime} =$ 0.025 G
	Le	α1	${\gamma }_{\mathrm{br}}^{{\prime} }$	${\gamma }_{\max }^{{\prime} }$	${U}_{{\rm{e}}}^{{\prime} }$	${U}_{{\rm{B}}}^{{\prime} }$	${U}_{{\rm{e}}}^{{\prime} }$/${U}_{{\rm{B}}}^{{\prime} }$	Ljet
	(erg s−1)				(erg cm−3)	(erg cm−3)		(erg s−1)
Modified Naima	7.8 × 1043	2.5	1.9 × 105 a	1.5 × 106	5.3 × 10−4	2.5 × 10−5	21	8.9 × 1043
Jetset	8.2 × 1043	2.6	1.9 × 105 a	1.6 × 106	5.5 × 10−4	2.5 × 10−5	22	9.3 × 1043

Note. See Section 4.3.1 for further details. In the various model realizations, the radius of the emission region $R^{\prime}$ is fixed to 1.14 × 10¹⁷ cm, and the Doppler factor δ to 11. For the radiating electron distribution, a broken power law is used with the spectral indices fulfilling α₂ = α₁ + 1. The minimum energy of the electrons ${\gamma }_{\min }^{{\prime} }$ is set to 1000. The table reports, for the broken power law describing the shape of the electron distribution, the first spectral index α₁, the break energy ${\gamma }_{\mathrm{br}}^{\prime}$, and the maximum energy ${\gamma }_{\max }^{\prime}$ as well as the electron luminosity L_e, the energy densities held by the electron population ${U}_{{\rm{e}}}^{{\prime} }$, and the magnetic field ${U}_{{\rm{B}}}^{{\prime} }$, their ratio, and the total jet luminosity L_jet. For the EBL γ-ray absorption at a redshift z = 0.034, the model from Franceschini et al. (2008) is used.

^a Fixed to the cooling break.

Download table as:  ASCIITypeset image

The SED related to the NuSTAR-2 observation differs from that of the baseline emission only in the X-ray domain (both soft and hard X-rays). Hence one can describe it with a set of model parameters that are very similar to those from the baseline emission. The magnetic field strength can be kept constant at 0.025 G between the two states, and the parameters describing the electron distribution are almost identical, with a slightly higher ${\gamma }_{\max }^{{\prime} }$ for the NuSTAR-2 state. We note that the two leptonic frameworks employed agree well with each other, as shown in Figure 13(b).

On the other hand, the substantially higher emission during the NuSTAR-1 observation requires a decrease in the magnetic field strength to 0.01 G connected to an increase in ${\gamma }_{\mathrm{br}}^{{\prime} }$, so that the HE peak is not underestimated. Additionally, we need a slightly harder spectral index and higher ${\gamma }_{\max }^{{\prime} }$ with the higher flux state. The values of the corresponding model parameters are reported in Table 7, and the model curves are depicted in Figure 13(a).

The simple one-zone leptonic approach does explain the evolution of the SEDs reasonably well. However, we still want to test our hypothesis of the low state of Mrk 501 being its baseline emission and therefore employ a two-zone scenario building on the hypothesis.

4.3.2. Two-zone Model

In our second approach, we investigate the temporal evolution of the SED assuming the existence of two independent emission zones. The first is a stable and always present part of the SED as represented by the model describing the 2 yr long low-activity state described in Section 4.2.1. This baseline emission would be often outshone by other emission regions occurring at different positions in the jet, which show a larger degree of brightness and variability. Therefore, we fit the various SEDs assuming two independent zones: one with the properties fixed to the ones given in Table 5 and another zone that is variable.

For the additional region that contributes to the two broadband SEDs during the NuSTAR-1 and NuSTAR-2 observations, we again assume a minimum electron energy of ${\gamma }_{\min }^{\prime} ={10}^{3}$. For the NuSTAR-1 SED we choose a radius of $R^{\prime} =5\times {10}^{15}$ cm after trying a grid of $R^{\prime}$ expanding up to 5 × 10¹⁶ cm. The radius is limited to be smaller than 5 × 10¹⁶ cm to not interfere with the baseline region. This is based upon the assumption that the emission region could extend over a large fraction of the cross section of the jet, and if that is the case, the smaller and more active region should be closer to the central engine than the baseline region. As over time the active region travels along the jet, an upper limit has to be set for the radius, which relates to a minimum distance between the active and baseline regions to ensure our independent treatment of the two zones during the time interval considered in this study.

Assuming a conical jet model, we follow Zdziarski et al. (2022) and assume a constant bulk Lorentz factor and therefore δ = 11, as obtained for the low-state model, along the jet.

For the NuSTAR-2 SED we allow $R^{\prime}$ to expand with an upper limit related to the expected maximum change for a jet opening angle of 1/Γ_b ≈ 1/δ and the 27 days between the two spectra. The electron density ${n}_{{\rm{e}}}^{{\prime} }$ is made dependent on the change in radius, assuming spherical symmetry for the blob and a constant number of electrons for the different models. Again, a broken power-law distribution is used for the radiating electrons, the break energy is fixed to the cooling break, and the two spectral indices, α₁ and α₂, are assumed to have a difference of 1. As the first investigations showed a strong preference of the magnetic field to be around the same value as for the NuSTAR-1 state, we fixed it to the same value during the evolution.

Owing to the very good agreement in the results obtained with the two leptonic software packages employed above, for the sake of simplicity, this time we decided to use only the jetset package to conduct the modeling of the data. The model parameters that describe well the temporal evolution of the data are reported in Table 8, and the model curves (for the two zones together and separated) are depicted in Figures 14(a) and (b).

Zoom In Zoom Out Reset image size

Table 8. Parameter Values from the Active Zone in the Two-zone Leptonic Scenario Used to Describe the Broadband SED of Mrk 501 During the NuSTAR-1 Observation on 2017 May 28 (MJD 57,898) Shown in Figure 14(a), the NuSTAR-2 Observation on 2017 May 25 (MJD 57,871) Shown in Figure 14(b), the Low-activity State of Mrk 501 Derived with Data from 2017 June 17 to 2019 July 23 (MJD 57,921 to MJD 58,687) Shown in Figure 14(c), and the Typical (Nonflaring) Activity State Derived with Data from 2009 March 15 to 2009 August 1 (MJD 54,905 to MJD 55,044) as Reported in Abdo et al. (2011a) and shown in Figure 14(d)

	$R^{\prime}$ (cm)	${n}_{{\rm{e}}}^{{\prime} }$ (cm−3)	$B^{\prime}$ (G)	α1	${\gamma }_{\min }^{{\prime} }$	${\gamma }_{\mathrm{br}}^{{\prime} }$	${\gamma }_{\max }^{{\prime} }$	${U}_{{\rm{e}}}^{{\prime} }$ (erg cm−3)	${U}_{{\rm{B}}}^{{\prime} }$ (erg cm−3)	${U}_{{\rm{e}}}^{{\prime} }$/${U}_{{\rm{B}}}^{{\prime} }$	Ljet (erg s−1)
NuSTAR-1	5.0 × 1015	6.5	7.4 × 10−2	2.0	1 × 103	4.8 × 105 a	2.8 × 106	3.6 × 10−2	2.6 × 10−4	139	1.2 × 1043
NuSTAR-2	2.9 × 1016 b	3.4 × 10−2 b	7.4 × 10−2	2.1	1 × 103	1.5 × 105 a	1.8 × 106	1.5 × 10−4	2.2 × 10−4	0.7	3.7 × 1042
Low State	5.0 × 1016 c	6.5 × 10−3 b	7.4 × 10−2	2.5	1 × 103	8.6 × 104 a	2.6 × 106	1.5 × 10−5	2.2 × 10−4	0.1	6.6 × 1042
Typical State	1.8 × 1016	33.5	2.7 × 10−2	2.1	100	5.1 × 105 a	5.1 × 106	1.8 × 10−2	2.9 × 10−5	621	8.7 × 1043

Note. The baseline broadband emission (first zone) is described with the one-zone SSC model of the low state reported in Table 5. For the active (variable) emission zone, the radiating electron distribution is assumed to be a broken power law with the spectral indices fulfilling α₂ = α₁ + 1. The Doppler factor δ is set to 11. The table reports the magnetic field $B^{\prime}$ and radius $R^{\prime}$ of the emitting region and, for the assumed broken power law describing the shape of the electron distribution, the electrons density ${n}_{{\rm{e}}}^{{\prime} }$, the first slope α₁, and the minimum, maximum, and break energies ${\gamma }_{\min }^{{\prime} }$, ${\gamma }_{\max }^{{\prime} }$, ${\gamma }_{\mathrm{br}}^{{\prime} }$. Additionally, the table also reports the energy densities held by the electron population ${U}_{{\rm{e}}}^{{\prime} }$ and the magnetic field ${U}_{{\rm{B}}}^{{\prime} }$, their ratio, and the total jet luminosity L_jet. For the EBL γ-ray absorption at a redshift z = 0.034, the model from Franceschini et al. (2008) is used.

^a Fixed to the cooling break. ^b Fixed to expanding $R^{\prime}$ assuming a spherical blob. ^c Fixed to same expanding velocity as obtained from the fits above.

Download table as:  ASCIITypeset image

We use this scenario to further expand the region for 24 more days until the historically low activity (baseline) starts, and fit the low-state SED with a radius fixed to the same expansion velocity as determined with the previous fits. The same assumptions for the electron density and the Doppler factor as before are used. The result is shown in Figure 14(c) and agrees with the data. Besides the change in the value of $R^{\prime}$, the injection mechanism becomes less energetic between the three states: ${\gamma }_{\max }^{{\prime} }$ stays rather constant, but the spectral index softens with decreasing flux.

Additionally, we use the abovementioned two-zone theoretical scenario to describe the typical broadband SED of Mrk 501 reported in Abdo et al. (2011a). For this purpose we use the broadband SED data points depicted in Figure 8 of Abdo et al. (2011a), but this time excluding the time interval MJD 54,952–54,982 when deriving the Fermi-LAT spectrum to avoid a spectral hardening above 10 GeV caused by a flaring episode in 2009 May (see Figure 9 of Abdo et al. 2011a, and discussion in Section 5.3). We assume the same connection between δ and $R^{\prime}$ as above but leave all other parameters free to vary. The radius is again limited to be smaller than 5 × 10¹⁶ cm. For the electron distribution, we used again a broken power law. The result is shown in Figure 14(d). The theoretical model describes the SED data quite well, except for the peaky structure at the lowest end of the X-ray data, whose exact description (within less than 15% accuracy) would require a complex shape for the electron distribution. We consider that, given the relatively small magnitude of the model–data disagreement, together with systematic uncertainties for the X-ray data (which is already at the level of 10%–15%, as reported in Abdo et al. 2011a), this extra complexity in the theoretical model is not justified.

Section 3 discusses the variability and correlations connected to the low activity of Mrk 501. Deduced from the obtained fractional variability in Figure 3, the main changes between the 4 yr period 2017–2020 and the 2 yr period with historically low activity take place mainly in the HE X-rays and the VHE γ-rays. Usually, an increase in F_var with energy is seen for Mrk 501. However, this behavior is more prominent when flaring episodes are included in the evaluated data set (Ahnen et al. 2017a, 2018; Acciari et al. 2020a) than when looking at nonflaring activities (Aleksić et al. 2015a; Ahnen et al. 2018). The F_var from the 4 yr data set shows a double-peak structure with the HE X-rays and VHE γ-rays showing the highest variability at a similar level. The same was reported in Furniss et al. (2015) for Mrk 501. F_var computed with the 2 yr low-state period also increases with the energy, but it reaches a plateau at X-rays and γ-rays. A similar behavior, but limited by the sensitivity at that time, is reported in Aleksić et al. (2015a).

What should be noted is that, for all mentioned previous results, far shorter time periods are taken into account than our 2 yr and 4 yr data sets. Recently, the FACT collaboration (Arbet-Engels et al. 2021) has reported a study performed with a Mrk 501 data set that spans over 5.5 yr, from 2012 December to 2018 April, and therefore overlapping with the 2017–2020 data set featured in this paper. The study performed by the FACT collaboration also reports a double-bump structure in the fractional variability versus energy, with F_var values that are very similar to the ones reported here, except for higher values in the VHE regime. The latter result is not surprising because of the very large VHE activity shown by Mrk 501 in the year 2014 (Cologna et al. 2017; Acciari et al. 2020a).

Hence, we can conclude that the variability pattern of Mrk 501 is dominated by variations of the VHE γ-rays and X-rays, with certain periods showing large and different variability levels, while other time epochs show comparable variability levels. This observation would be consistent with the broadband emission described by a multiple zone scenario where, for the higher flux states, different variability patterns are influenced by different active regions. This variability pattern would be naturally expected from the two-zone model described in Section 4.3.2, where the active (smaller) region would dominate the emission at X-rays and VHE, and hence the most variable parts of the SED of Mrk 501.

Our study of possible correlations between the different wave bands, identifies a clear correlation without time lag between the LE/HE X-rays and VHE γ-rays (see Section 3.2). This correlation has been reported multiple times for Mrk 501 during flaring activities (see, e.g., Furniss et al. 2015; Ahnen et al. 2018; Acciari et al. 2020a), but is reported here for the first time with a statistical significance above 3σ for an extensive period of low activity. It further supports that those scenarios where X-ray and VHE photons are produced by the same particle population, such as in the case for an SSC model, and dominate the broadband emission of Mrk 501 during all kinds of activity states.

In addition to the correlation with zero time lag, we measure a 3σ (DCF value >3σ confidence level) X-ray versus VHE correlation for a time delay of roughly 30 days. This correlation exists also when comparing the 0.3–2 keV with the 2–10 keV X-ray fluxes measured with Swift-XRT, with the HE band preceding the LE band. This feature is visible (although with marginal significance) in both the 4 yr and the 12 yr data sets, once the light curves are detrended to remove the long-term behavior. It is the first time that such a observation is made indicating that a fraction of the X-rays are produced with some delay, in particular those at the lowest energies. This suggests that the acceleration of the high-energy particles may be produced by more than one process, with a potential connection between them. However, a word of caution is in order. The timescale of about 30 days coincides with the timescales where an increase in the confidence levels in our periodicity study (see Section D) is apparent. Such increase in the confidence levels is probably related to the observing sampling of the source, which is affected by the moon periods that constrain the observations from the MAGIC telescopes (and hence all MWL observations that were coordinated with MAGIC observations). While we tried to include all these effects in our MC simulations used to determine the confidence levels, we cannot exclude the existence of some nonaccounted for effects that may introduce small artifacts at this timescale.

For the UV versus R band, the DCF analysis yields a clear peak centered at zero time lag, but relatively broad, extending over ±1 week. This implies that these two neighboring energy ranges have the same behavior, with a variability timescale on the order of about one week, probably due to this emission being produced by relatively low-energy electrons, which have longer radiation timescales.

The extensive data set spanning from 2008 to 2020 at HE γ-rays with Fermi-LAT, X-rays with Swift-XRT, and radio with OVRO, allowed us to study potential correlations among these bands with unprecedented precision. While the 4 yr data set does not yield conclusive results, the 12 yr data set shows 3σ correlations between the HE γ-rays and the X-rays, as well as between HE γ-rays and radio (See Figures 6 and 7). Such correlated behavior had not been reported previously at a statistical significance above 3σ.

In the case of HE γ-rays versus X-rays, the DCF analysis shows a correlation that is highest at about zero time lag (for both X-ray energy bands, 0.3–2 keV and 2–10 keV), but with a very broad peak that extends ±1 yr. This positive correlation, extending over a large range of time lags, is ascribed to the long-term flux increase for multiple years, together with a decrease around the year 2017, which is observed in the LCs from all these instruments (see Figure 1). A time shift of several months in these LCs would still keep the overall long-term behavior, and hence the positive correlation obtained by our study. When the long-term flux variations in our light curves are removed, as described in Section 3.2, the correlation plots yield a single bump centered at a time lag zero, and with a width of about ±10 days (see Figures 7(c) and (d)). The DCF value for a time lag zero is clearly above the 3σ confidence level and hence statistically significant. Overall, the correlation plots from Figure 7 show that the keV and the GeV emissions are clearly correlated on both long (months, years) and short (weeks) timescales. This indicates that the radiation at these two energy bands is produced, at least partially, by the same population of particles, which further supports the SSC scenarios for the variable emission of this source. Further support for the SSC scenario is given by the correlation between VHE and HE γ-rays revealed after removing the long-term trend from the 4 yr data set.

On the other hand, in the correlation between OVRO and Fermi-LAT, the radio lags behind the γ-rays by more than 100 days, with slightly higher DCF values at ∼126 days and 238 days. Using a 5.5 yr data set from Mrk 501, from 2012 to 2018, the FACT collaboration (Arbet-Engels et al. 2021) has recently shown a positive correlation between the Fermi-LAT fluxes and those from OVRO, with a relatively flat behavior and with time lags extending from −300 days to +300 days; the significance of such correlated behavior was not computed (see Figure 9 of Arbet-Engels et al. 2021). In our study, that employs a 12 yr data set, we show that the correlation is statistically significant (> 3σ) only for time lags larger than −100 days, implying the radio emission is connected to the γ-rays with a delay larger than 3 months. This correlation, however, disappears when the light curves are detrended (see Figure 6(b)), hence indicating that the radio and γ-ray emission are related only when considering the flux variations with timescales of a few months. A delay between radio and γ-ray fluxes with time delays of several tens or even hundred days has been reported with a statistical significant larger than 3σ for other blazars (Max-Moerbeck et al. 2014; Acciari et al. 2021). Such delays are usually explained by moving disturbances that travel along the jet. The time delay of the radio emission can then be converted to the distance between the locations dominating the radio and γ-ray emission using Equation (1) in Max-Moerbeck et al. (2014). With a maximum jet speed of β = 0.9 and a maximum Doppler factor in the radio regime of δ = 2 (Lister et al. 2021) we obtain the bulk Lorentz factor Γ = 1.45 using Equation (4) in Hovatta et al. (2009). Hence, for the time delay of 126 days (238 days) the emission regions are at maximum 0.27 pc (0.51 pc) apart. Following the recipe for the calculations of the distance from the central engine of the radio core emission d_core for Mrk 421 in Max-Moerbeck et al. (2014) replacing the radio core size with 0.13 mas for Mrk 501 (Weaver et al. 2022) we obtain ${d}_{\mathrm{core}}=2.05$ pc. Hence, the γ-ray emission region is at least 1.78 pc (1.54 pc) away from the central engine, and located closer to the radio emission than the turbulent inner regions of the blazar, as already found for Mrk421 (Max-Moerbeck et al. 2014; Acciari et al. 2021). However, the Doppler factor of 11 obtained for the γ-ray emission site in Section 4.2 is in contradiction to assuming a constant δ = 2 for the whole region between the two emission sites. Assuming δ = 11 for the whole region would give a distance between the sites of 5.66 pc (10.69 pc), which is bigger than d_core and can therefore be excluded. If we consider a linear decrease from δ = 11 to δ = 2 between the two sites, we obtain a separation of 2.59 pc (4.88 pc). Already this simplified assumption relaxes the overshooting of d_core. More developed models such as decelerating jets have been proposed before (Meyer et al. 2011) to explain the tension between Doppler factors measured in the radio regime and the ones needed to explain MWL SEDs.

The DCF analysis derived with the HE γ-ray data from Fermi-LAT and the optical R band yields a hint of correlation, although not significant. In the SSC models, one would expect a direct correlation between the eV and the GeV emission, as these two energy bands are produced by the same particle population, and hence are good tracers of the dynamics of these particles. The lack of a clear correlation could be explained by the existence of additional contributions to these energy bands, perhaps coming from different regions which are not physically connected. This would yield an uncorrelated behavior that would worsen the correlated behavior expected from the most simple one-zone SSC scenarios. Additionally, it could point to a hadronic nature of the emission, which would not show a correlation between the two wave bands.

The absence of correlation observed between the polarization degree and the multiband fluxes could be a sign of different mechanisms than shock acceleration being at work in the jet (Jorstad et al. 2007). Additionally, the missing correlation between optical and radio polarization degrees and angles is hinting toward different emission zones producing the radiation. Nonetheless, as both optical and radio emission are known to have additional components, apart from the main blazar emission, the correlation could be partially washed out. Together with the rather sparse radio coverage, this does not allow us to draw significant conclusions. For the polarization angle, the preferred angles in the optical R-band and the 43 GHz radio measurements coincide with the jet direction previously determined for Mrk 501 (Weaver et al. 2022). The agreement between the measured polarization angles and the obtained jet direction could point toward a magnetic field perpendicular to the jet direction contributing to the collimation of the jet. However, the situation might be more complex taking into account relativistic effects (Lyutikov et al.2005).

The 2017–2020 data set includes a 2 yr time interval with the historical low activity of Mrk 501 from mid-2017 to mid-2019. This provides us with a remarkable opportunity to study the baseline emission without disturbances from other contributions that are more variable and normally dominate the broadband emission of an active Mrk 501. From previous correlation studies, as well as the ones included in this work, the preferred scenario for explaining the variable part of the blazar emission is of leptonic origin. However, for the baseline and thus stable part of the emission, this does not necessarily apply. Therefore, we consider both electrons and protons to be possible emitters for the observed low-state emission.

From the modeling results in Section 4 we can see that both relativistic electrons as well as protons can explain reasonably well the observed low-activity SED. Indeed, even the hadronic model can be constructed with a low magnetic field of 3 G preserving the required low variability. The consideration of including hadronic components in the blazar emission origin is of particular importance as it allows us to estimate the possible contribution from blazars to the flux of neutrinos and ultra-high-energy cosmic rays (UHECRs; E>10¹⁸ eV). Taking the multimessenger picture that is nowadays provided by neutrino telescopes like IceCube into account, we can use the upper limits in Aartsen et al. (2020) to check our predicted neutrino emission. The obtained neutrino rates of 10⁻⁵ per year in our hadronic models are well below the upper limit of 1 neutrino per year (10.3 best-fit neutrino in the first 10 yr of IceCube data). For the lepto-hadronic models, higher neutrino rates are expected. However, our lepto-hadronic models prediction of 10⁻³ neutrinos per year is perfectly consistent with the nondetection of Mrk 501 by IceCube. The higher neutrino flux is accompanied with higher estimated jet powers compared to the hadronic case. This can be accounted for by fine-tuning the low-energy part of the proton distribution. Further, the hadronic scenario provides very high proton energies up to ${\gamma }_{\max }^{\prime} \sim {10}^{10}$ potentially indicating BL Lac objects as UHECR accelerators. For the available data set, leptonic, hadronic, and lepto-hadronic models are valid scenarios to explain the observed low-activity state of Mrk 501. As one of the currently missing pieces of the multimessenger blazar picture, data in the MeV energy range would be extremely efficient in reducing the degeneracy among the various emission models. Furthermore, future X-ray/γ-ray polarization measurements would help in distinguishing between theoretical models due to the different predictions in polarization properties between the models.

To explain such a long-lasting stable emission, a standing shock scenario (see, e.g., Marscher et al. 2008; Marscher 2014) provides a reasonable scenario as previously discussed for the quiescent behavior of Mrk421 in Abdo et al. (2011b). Shock acceleration has been shown before to be a plausible scenario explaining the behavior of Mrk 501 (Baring et al. 2016) and would be in line with the obtained spectral indices in all our applied models. It has been further strengthened for Mrk 501 by the first X-ray polarization measurements reported in Liodakis et al. (2022), whose higher polarization degree compared to the optical measurements supports a shock acceleration scenario. In this scenario, the particles are accelerated when the jet flow crosses the standing shock, and subsequently radiation is emitted. As long as the particle flow and shock properties remain stable, constant acceleration and emission take place.

In what concerns the time evolution of the SED in the months before the historically low activity that starts in mid-2017, only leptonic scenarios (one-zone and two-zone models) are considered to explain the variations in the broadband emission, because of the tight correlations between X-rays and γ-rays (both HE and VHE).

In earlier studies, transitions between different states were commonly attributed to changes in the break energy ${\gamma }_{\mathrm{br}}^{\prime}$ for Mrk 501, and therefore to the injection of electrons (Anderhub et al. 2009; Ahnen et al. 2018; Acciari et al. 2020a). Within the one-zone leptonic scenario described in Section 4.3.1, the most important parameter change occurs in the magnetic field $B^{\prime}$, which increases from 0.01 to 0.025 G as Mrk 501 transitions toward the historical low activity in mid-2017. The observed changes in $B^{\prime}$ go hand in hand with changes in ${\gamma }_{\mathrm{br}}^{{\prime} }$ as they are linked through the cooling break in our scenario. Furthermore, the electron distribution requires small adjustments with the power-law indices α₁ (which is linked to α₂) becoming softer, and the maximum energy ${\gamma }_{\max }^{{\prime} }$ decreasing with time, as Mrk 501 reaches the low activity. We propose that a small increase in the ambient parameter $B^{\prime}$ explains the observed broadband SED time evolution when the injected electron distribution flowing through the shock adjusts to the surrounding.

For all three emission states in the 4 yr data set, we determined the synchrotron peak frequency ν_s . First we used the phenomenological description in Ghisellini et al. (2017) applied to all data points to determine ν_s . Additionally, we exploited the one-zone leptonic and hadronic modeling results, which are described further below, for further estimates of ν_s . Table 13 summarized the peak frequencies, which for the low state are on the order of 5 × 10¹⁵ to 5 × 10¹⁶ Hz. The differences between the estimates have their root in the treatment of the low-energy points as upper limits for the theoretical models while for the phenomenological fits this assumption is not made so as to be able to compare it to previously published results more easily. However, the synchrotron peak is not very well covered by our data set and can therefore not be identified more specifically. While for the NuSTAR-2 the peak frequencies hardly shift, the NuSTAR-1 state shows a clear shift to higher frequencies between 2 × 10¹⁶ and 3 × 10¹⁷ Hz. The phenomenological evaluation clearly places Mrk 501 in the HSP regime.

In order to test our assumption of the low state being the baseline emission of Mrk 501, our second scenario assumes two emission zones as described in Section 4.3.2. We consider that there is one region that is responsible for the historical low-activity broadband SED, and a distinct and independent region that is responsible for the main variations in the SED. These variations can occur on timescales of days or weeks, and are particularly dramatic in the X-ray and γ-ray energy range. The low-activity broadband emission could be produced by high-energy electrons or high-energy protons, as demonstrated in Section 4.2. As it is assumed to be steady (in reality there could be small long-term flux variations, which are not considered here), it is irrelevant whether we use a leptonic or hadronic framework in our two-zone scenario. For simplicity, we used the one-zone SSC scenario described in Section 4.2.1. Because of the abovementioned flux variability and correlations, as we do with the one-zone scenario, we consider that the broadband emission produced in the second region is dominated by high-energy electrons. We assume a region with a size $R^{\prime}$ smaller than that of the baseline one, and with a higher $B^{\prime}$ field, potentially indicating a position closer to the central engine. Over time, the active region expands, the density of radiating particles becomes smaller, and the spectral index softens. If we assume a dependence of the jet magnetic field on the size of the blob using the relation in Tramacere et al. (2022), we can conclude that the variable region stays at a stable position in the jet when expanding. Additionally, it is only valid if the radius of the variable region makes up only a fraction of the low-state radius, which suggests that at least the variable region does not cover the whole jet radius. However, the magnetic field could also be dominated by local contribution traveling together with the shock region inside the jet.

The radio emission is not reproduced by the models as it is assumed to originate from a more extended region in the outer part of the jet than considered here. Therefore, at least a third if not further multiple or more complex zones would be required to reconstruct the full SEDs as, for example, demonstrated by Lucchini et al. (2019). Other alternative models could potentially also reproduce the low Doppler factor observed in the radio regime as shown by Ghisellini et al. (2005) using a structured jet or with a decelerating jet by Georganopoulos & Kazanas (2003). It is, however, remarkable that the flux density predicted by the baseline model matches reasonably well with what is observed with millimeter-wavelength VLBI. Using global millimeter-VLBI array observations, Giroletti et al. (2008) measured a flux of S_{86 GHz} ∼ 45 mJy for the central component seen at the jet base, corresponding to ν F_ν ∼ 4.0 × 10⁻¹⁴ erg cm⁻² s⁻¹. The deconvolved size of such a emission region is smaller than ∼5 × 10¹⁶ cm (in the observer's frame). It is therefore possible that at least the low-state baseline VHE emission has a counterpart that is directly accessible with millimeter-VLBI (while the more compact components responsible for the variable emission remain self-absorbed).

This two-zone scenario can also be used to describe the typical (nonflaring) SED of Mrk 501 derived with data from 2009 and reported in Abdo et al. (2011a). To account for the higher radio flux (at 15 GHz) in 2009 with respect to the low-activity state after mid-2017 (see Figure 1), while considering the correlations between the radio and the γ-rays reported in Section 3.2, we had to decrease the minimum energy of the high-energy electrons, down to ${\gamma }_{\min }^{{\prime} }=100$. The active region produces a radio emission that is comparable to that of the baseline region (see Figure 14(d)), and hence account for variations in radio, X-rays, and γ-rays. In this realization of the model, the active region is quite far from equipartition (${U}_{{\rm{e}}}^{{\prime} }/{U}_{{\rm{B}}}^{{\prime} }\gt 5\times {10}^{2}$), but well within the values that are considered possible for HSPs such as Mrk 501 (Ahnen et al. 2017a, 2018; Acciari et al. 2020b). Such high deviation from equilibrium has been, for example, consistently reproduced using relativistic, oblique, magnetohydrodynamic shocks by Baring et al. (2016).