Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Next Article in Journal
Decoding Quantum Gravity Information with Black Hole Accretion Disk
Previous Article in Journal
The Mother’s Day Solar Storm of 11 May 2024 and Its Effect on Earth’s Radiation Belts
Previous Article in Special Issue
Optical Quasi-Periodic Oscillation of Blazar PKS 1440-389 in the TESS Light Curve
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Searching for Hadronic Signatures in the Time Domain of Blazar Emission: The Case of Mrk 501

by
Margaritis Chatzis
1,2,*,
Stamatios I. Stathopoulos
1,2,*,
Maria Petropoulou
1,2 and
Georgios Vasilopoulos
1,2
1
Department of Physics, National and Kapodistrian University of Athens, University Campus Zografos, 15784 Athens, Greece
2
Institute of Accelerating Systems & Applications, University Campus Zografos, 15784 Athens, Greece
*
Authors to whom correspondence should be addressed.
Universe 2024, 10(10), 392; https://doi.org/10.3390/universe10100392
Submission received: 16 August 2024 / Revised: 9 September 2024 / Accepted: 26 September 2024 / Published: 10 October 2024
(This article belongs to the Special Issue Blazar Bursts: Theory and Observation)

Abstract

:
Blazars—a subclass of active galaxies—are intrinsically time-variable broadband sources of electromagnetic radiation. In this contribution, we explored relativistic proton (hadronic) signatures in the time domain blazar emission and searched for those parameter combinations that unveil their presence during flaring epochs. We generated time series for key model parameters, like magnetic field strength and the power-law index of radiating particles, which were motivated from a simulated time series with statistical properties describing the observed GeV gamma-ray flux. We chose the TeV blazar Mrk 501 as our test case, as it had been the study ground for extensive investigations during individual flaring events. Using the code LeHaMoC, we computed the electromagnetic and neutrino emissions for a period of several years that contained several flares of interest. We show that for both of those particle distributions the power-law index variations that were tied to moderate changes in the magnetic field strength of the emitting region might naturally lead to hard X-ray flares with very-high-energy γ -ray counterparts. We found spectral differences measurable by the Cherenkov Telescope Array Observatory at sub-TeV energies, and we computed the neutrino fluence over 14.5 years. The latter predicted ∼0.2 muon and anti-muon neutrinos, consistent with the non-detection of high-energy neutrinos from Mrk 501.

1. Introduction

Blazars, a prominent subclass of active galaxies with jets aligned with our line of sight [1], are sources of broadband electromagnetic radiation that exhibit intrinsic flux variability. The flux variability is observed across the electromagnetic spectrum, from radio frequencies to very-high-energy gamma rays (for a review, see [2]). The flux may vary on both short (minutes to hours) and long (months to years) timescales [3,4].
The Spectral Energy Distribution (SED) of blazars typically features two broad components: a low-energy peak attributed to synchrotron radiation from relativistic pairs, and a high-energy peak often explained by the inverse Compton scattering of synchrotron photons [5,6] or by external photon fields [7,8] in leptonic models. In hadronic models, on the other hand, high-energy emission results from the interaction of relativistic protons with ambient photon fields, producing secondary particles, such as pions and muons, that decay into secondary pairs, gamma-rays, and neutrinos. While time-averaged SEDs provide valuable insights into the average physical conditions within the blazar’s emission region, they do not capture the dynamic aspects of blazar behavior. During variability epochs, we can obtain more precise constraints on source parameters, such as the radius of the emission region or the magnetic field strength, rather than average values. This approach allows for a more detailed understanding of the physical conditions inside the source [9,10].
The correlation between flux variations in different energy bands is a tool for distinguishing between emission models. In leptonic models, the Synchrotron Self-Compton (SSC) mechanism predicts a quadratic relationship between the X-ray and TeV gamma-ray fluxes, as both emissions originate from the same electrons. For example, ref. [11] was able to reproduce the SED of Mrk 501 during its April–May 1997 outburst with a time-dependent SSC model by varying the electron luminosity, the maximal Lorentz factor of the electrons, and the Doppler factor. In leptohadronic models, the correlations can be more complex, due to the additional processes involved, such as photohadronic interactions, synchrotron radiation from secondary particles, and photon–photon annihilation. These models may predict different time lags and flux–flux correlations compared to leptonic models [12], providing a means to distinguish between the two scenarios.
Mrk 501 is a well-known TeV blazar that has been extensively studied, due to its remarkable flaring activity and the wealth of observational data available. Previous studies, such as the detailed analysis of the 1997 flaring period [13], have revealed significant correlations between X-ray and TeV gamma-ray fluxes, consistent with SSC predictions [11]. During X-ray flares, narrow spectral features in the very high energy (VHE, E > 100 GeV) can be explained by the emission of neutral pion decay, as shown in ref. [14]. However, the variability when combined with the multi-wavelength campaigns suggests that the observed emission from Mrk 501 could be due to a complex superposition of two different emission zones [15].
In this study, we aimed to explore the time-dependent behavior of blazar emission, focusing on Mrk 501 as our test case. Using the LeHaMoC code [16], we simulated time series for key model parameters, such as electron–proton distribution and magnetic field strength, motivated by the observed gamma-ray flux variability. Our primary objective was to identify parameter combinations that revealed hadronic signatures during flaring epochs and to search for energy windows where hadronic emissions might dominate. We analyzed electromagnetic and neutrino emissions over several years, encompassing two simulated flaring events, to investigate flux–flux correlations between different energy bands.
The structure of this paper is as follows: In Section 2, we describe briefly the physical problem and the numerical code. Then, we focus on the case study of Mrk 501, detailing the observational data and the modeling of the time-average SED that served as the baseline of our time-dependent model. In Section 2.4, we describe how we created synthetic light curves. In Section 3, we present our results on the variability of the broadband emission, discussing the identified hadronic signatures and flux correlations. Section 4 provides predictions for future observations with the Cherenkov Telescope Array (CTA), and Section 5 discusses our findings in the context of neutrino observations by IceCube. In Section 6, we discuss our findings, we present some caveats of our analysis, and we outline potential improvements. We conclude in Section 7.

2. Methods

2.1. Numerical Code

We modeled our emitting region as a spherical blob of size R within the jet of the AGN, and we assumed a homogeneous magnetic field of strength B. In blazars, the angle θ between the line of sight to the observer and the bulk plasma velocity is a few degrees [17,18]. Emissions produced in the blob’s rest frame will therefore appear to be Doppler-boosted to an external observer. For a bulk plasma motion with Lorentz factor Γ and velocity β = v / c , we defined the Doppler factor as δ = Γ 1 ( 1 β cos θ ) 1 . Throughout this work, all parameters were measured in the rest frame of the blob, unless explicitly stated otherwise. We assumed the presence of an acceleration region from which particles (primary electrons and protons) escaped upon their acceleration into the blob where they radiated. The injected distribution of the accelerated particles was described by a power law with slope p j , d N j / d γ γ p j , between Lorentz factors γ min , j and γ max , j , where j = e , p . Non-thermal radiation was produced by the particle distributions inside the blob, due to various processes, such as synchrotron radiation and inverse Compton scattering. Inelastic collisions between relativistic protons and photons also resulted in the creation of secondary particles, such as relativistic electron–positron pairs, neutrinos, anti-neutrinos, and  π 0 which decay into γ γ pairs. To quantify and track the time evolution of a population j in the blob, with  differential distribution N j ( γ ) , we solved a system of integrodifferential equations, also known as kinetic equations. The equation of each particle species can be parameterized as
N j ( γ , t ) t + N j ( γ , t ) τ e s c + i L j i ( γ , t ) = i Q j i ( γ , t ) ,
assuming a physical escape from the source on a timescale τ e s c = R / c . In the above, symbols Q j i and L j i denote the mathematical operators for the sources/injections and sinks/losses of particles species j due to the physical process i, respectively. In our calculations, we accounted for electron and proton synchrotron emission and absorption, inverse Compton scattering (for electrons only), γ γ pair production, photopair (Bethe–Heitler pair production), and photopion production processes. For primary electrons and protons, there was an additional source term Q e , p e x t that described the rate of particle injection into the blob. The associated energy injection rate in relativistic electrons or protons was then obtained as L e , p ( t ) = m e , p c 2 d γ γ Q e , p e x t ( γ , t ) . To solve the kinetic equations of the stable particles in the blob—namely, electrons (and positrons), protons, photons, and neutrinos (and anti-neutrinos)—we utilized the code LeHaMoC1 [16], a versatile time-dependent leptohadronic modeling code.

2.2. Mrk 501 as a Case Study

Mrk 501 is a High-Synchrotron Peaked (HSP) blazar at a luminosity distance of d L = 149.4 Mpc that corresponds to a redshift of z = 0.034 [19,20]. Due to its proximity and time variability, it has been the focus of numerous observational campaigns ([21,22,23,24,25,26]). Its well-defined flares in X-rays and VHE (Very-High Energy; >100 GeV) γ -rays allow for precise dissections in the time domain for modeling its low and high activity. Furthermore, Mrk 501 features different flux variability across the multi-wavelength spectrum and exhibits spectral changes, especially in X-rays, between low-activity (quiescence) and high-activity (flaring) epochs. It is also characterized as a transient extreme HSP blazar. The 1997 flare is an excellent example of such behavior. During this outburst, at X-ray energies, the peak synchrotron energy exceeded 100 keV, whereas at VHEs, Mrk 501 exhibited drastic flux variations confirming a spectrum above 10 TeV (e.g., [27,28,29,30]). Characteristics of an extreme blazar were also reported during non-flaring epochs in Mrk 501 [31], suggesting that this was not a permanent state of the source.
For our analysis, we used the archival observations available to us (accessed through the SED builder website2), supplemented by additional X-ray spectral data. Our final data set for Mrk 501 is presented in Figure 1, where we use color to differentiate between archival observations (shown in gray) and specific datasets, which will be discussed in the next paragraph.
The VHE spectral data of the 1997 outburst observed by HEGRA [30] alongside observations from the ARGO-YBJ experiment during the flare of 2011 [32] define the VHE flaring state. These are highlighted in Figure 1 as transparent gray points. Non-flaring-(quiescent)-state VHE data were provided for Mrk 501 during 2006 by the MAGIC [33] and TACTIC [34] observatories, during 2009 by VERITAS [35], between 2008 and 2012 by the ARGO-YBJ experiment [32], and between 2008 and 2015 by the Fermi Large-Area Telescope (LAT) [36]. For the HE range (High Energies; >0.1 GeV), we used data provided by the various releases of the Fermi Gamma-ray Space Telescope. In particular, we display data from the 1-, 2-, 4-, and 14-year Fermi LAT Sources Catalogs [37,38,39,40], as well as the monthly spectral data of the 2-year release. We highlight the 14-year observations of 4FGL with green color, as these provide the most reliable long-term representation of Mrk 501 in the HE range.
In the X-ray range, the archival data consists mainly of BeppoSAX observations between 1997 and 2001 [41,42], including the extreme outburst of 1997 (highlighted in light gray in Figure 1). The rest of the X-ray archival observations are shown with a faded black color. The archival data in the X-rays are only used for display purposes. To obtain the time-average broadband X-ray spectrum that we used for modeling, we supplemented our analysis with long-term observations from the X-ray Telescope (XRT, 0.2–10 keV) and the Burst Alert Telescope (BAT, 14–195 keV) on board the Neil Gehrels Swift Observatory. In particular, we used the Swift-BAT 157-Month Hard X-ray Survey Catalogue3 and XRT data from the Neil Gehrels Swift Observatory analysis tool4. We divided the XRT Light Curve (LC) (see Appendix A) into 3 segments, T < 55,500 MJD, 55,500  MJD < T < 58,000  MJD, and  T > 58,000 MJD. The first region exhibited low activity while the latter 2 capture the rise and decay of a high-activity period. We also used only the Windowed Timing (WT) mode data and not the Photon Counting (PC) mode, as the latter was saturated. For each region, we built the spectrum using online tools5, and we fitted an unabsorbed power law to the data, resulting in 15 energy flux data points for each epoch. The mean and standard deviation of the epochs was included in our data set as the representative XRT observations of the average spectrum. We also fitted an unabsorbed power law to the BAT spectrum provided by the BAT database. The Swift XRT and BAT data are highlighted with pink and orange colors, respectively, in Figure 1.
Lastly, we remove the optical bulge of the host galaxy from the archival data, which dominates the total emission at eV energies. One way of removing the bulge is to use SED templates of elliptical galaxies, which are fitted to the data and then removed. We opted for a model-independent way of removing the bulge. We isolated the optical, infrared, and radio parts of the spectrum, and we fitted a third-degree polynomial to the data outside of the visually selected bounds of the bulge assuming a non-thermal broadband component (e.g., attributed to the synchrotron emission from the jet). To define flux points belonging to the non-thermal radiation in the region of the optical bulge, we calculated the mean value of the squared residuals of the archival data and compared it to the individual residual value. We only included the data point if its flux was smaller than the mean, and we highlight the excluded observations with a faded black color in Figure 1.

2.3. Computation of Time-Average SED

To derive the baseline blazar parameters for our time variability analysis with LeHaMoC we followed a two-step approach. Recent scientific findings support that the average state of blazar emission is mostly of leptonic origin without distinct hadronic signatures [14,43,44]. Therefore, we first derived the parameters of a purely leptonic model, where the emission from infrared to X-rays is attributed to synchrotron radiation of primary electrons, and the HE to VHE γ rays are produced via SSC. The emission from relativistic protons, if present in the jet at all times, remains hidden below the leptonic spectral components. In the second step, we motivated the parameters for such a proton population.
For the leptonic average state, we did not include any flare data in our fitting procedure [21,30,32]. We also excluded low-frequency data (e.g., ≲300 GHz) for the fit, assuming that this emission originated from a more extended jet region6. This treatment was consistent with previous studies where excluded observations have been used as upper limits [44]. Additionally, observational studies have found that the radio region and the blob region are distinct [45,46], while other studies have modeled the radio region as a separate emission region [47].
We ran the Markov Chain Monte Carlo (MCMC) sampler emcee [48] for 32 parallel walkers on the prepared data. All our parameter priors used uniform distributions in logarithmic space except for the power-law index, which used a uniform distribution in linear space. We varied the step number of our fitting process throughout our testing, to allow for a thorough investigation of the parameter space. We concluded that 10,000 steps and a burn-in phase of 2000 steps were sufficient for the model to converge and have a meaningful posterior sample. In this fitting process, we chose to restrict the upper bound of the uniform radius prior to log( R m a x [cm]) = 16.5. This approach was conservative, in the sense that it predicted a variability timescale of t v a r = 1 d utilizing δ 20 , the median value of the Doppler factor posterior distribution (see Appendix B for a discussion of an unbound radius prior). We present in Table 1 our parameter results and in Figure 2 a representative sample of 100 leptonic SEDs computed from the posterior distribution of parameters. A corner plot showing the posterior distributions of the model parameters can be found in Appendix B.
Describing the relativistic proton population required four additional parameters. To reduce the number of free parameters, we assumed that both populations were accelerated before injection by the same mechanism into a power law distribution with the same slope and the same extent in Lorentz factor, namely, p p = p e and  γ p m i n / m a x = γ e m i n / m a x . In our treatment, we fixed these hadronic parameters to the median values of the posterior distributions listed in Table 1 (in Section 6, we discuss how another choice for relating the power law indices or the electron and proton energies would have affected our results). From our leptonic posterior analysis, we calculated the upper 1% limit, meaning 99% of our solutions at each frequency were below this SED (depicted in Figure 2 with a dashed line). To set an upper bound on the injection power of the relativistic proton component we expressed the latter as L p = 10 α L e and we increased α 0 until the computed SED exceeded by 10% the maximal leptonic SED (at any frequency). Our choice was conservative in the sense that it did not result in a leptohadronic SED with a distinct hadronic component. Our parameters for the underlying proton population are summarized in Table 2, while the SED components of the resulting average-state leptohadronic model are presented in Figure 3. As expected, the observed SED was dominated by synchrotron and IC emission of primary electrons. Synchrotron emission from secondary pairs from γ γ pair production and photopion production contributed to about 3% of the HE γ -ray flux. Similarly, the synchrotron emission of secondary pairs from Bethe–Heitler pair production, which peaked at about 100 keV, was a small fraction (∼6%) of the BAT flux. We note that the synchrotron photons from primary electrons provided the main seed photons for IC scattering by the secondary pairs.

2.4. Generation of Time Series

To model the multi-wavelength variability of Mrk 501 we first derived its statistical properties in the Fermi-LAT energy range and then created synthetic Light Curves (LCs) that were mapped into model parameter time series. The software celerite2 [51,52] was used to extract the statistical properties from the observational data, namely, the logarithmic energy fluxes (without upper limits) taken from the Fermi Light Curve Repository (LCR, [53]) with a mean cadence of 7 days, spanning a period of 15.4 years. Key features of celerite2 are speed and scalability, while its employed method, Gaussian process fitting, has been used to study the behavior of a variety of blazar sources [54,55,56,57,58]. Following the available literature on general celerite usage [51], as well as on applications on blazar time series specifically [54], we modeled the observed LC with a stochastically driven damped Simple Harmonic Oscillator (SHO), which is effective in detecting periodic signatures within a stochastic process. The differential equation of an SHO is
d 2 d t 2 + ω 0 2 Q d d t + ω 0 2 y ( t ) = ϵ ( t ) ,
where ω 0 is the undamped oscillator (or natural) frequency and Q is a dampening factor (also named quality factor). Assuming ϵ ( t ) to be a white noise term, the Power Spectral Density (PSD) is given by
S ( ω ) = 2 π S 0 ω 0 4 ( ω 2 ω 0 2 ) 2 + ω 2 ω 2 / Q 2 ,
where S 0 is proportional to the power at ω 0 and is given by S ( ω 0 ) = 2 / π S 0 Q 2 .
Although a mixture of SHOs can be used to describe and capture stochastic processes in an astrophysical setting we limited ourselves to the simplest case of one SHO. To derive its parameters, ω 0 , S 0 , and Q, we ran the Markov Chain Monte Carlo (MCMC) sampler emcee [48] for 32 parallel walkers with an initial 2000-step burn-in phase followed by a 5000-step phase for the posterior sampling used in our analysis. We found the walker and step size sufficient for the solutions to converge. As such, for the starting guesses for the MCMC sampler we utilized an L-BFGS-B optimization routine on a random starting guess.
The resulting parameters of the SHO alongside the mean μ of the analyzed time series can be seen in Table 3, while a visual representation of the fit and PSD is showcased in Figure 4. Furthermore, a detailed exploration of the goodness of our fit is presented in Appendix C. Thus, it is possible now to create a synthetic LC that accurately reflects the behavior of Mrk 501 in the Fermi energy range.

3. Results on Variability

To investigate time-variable SED models of Mrk 501 we used the celerite parameters listed in Table 3 and created a synthetic LC with statistical properties consistent with the observed LC in the Fermi energy band of 0.1–100 GeV. We selected a number of points equal to that of the real LC and multiplied them by the mean cadence of 7 days. However, since we did not include the upper limits of the Fermi LC in our analysis, the total number of points and, consequently, the total duration of the derived LC were shorter than those observed for about 230 days. As such, our method introduces uncertainties by neglecting the information provided by upper-limit detections. Additionally, while the generated LC accurately reflects the overall behavior of Mrk 501 in the Fermi energy band, it may display features that surpass the observed maximum or minimum flux within this range. Therefore, the LC produced by celerite2 should be interpreted as a potential scenario for the blazar’s behavior rather than a precise representation of its current state. Mitigation of this uncertainty will be discussed in Section 6.
With black stars, we showcase in Figure 5 the entirety of the zero-mean synthetic LC. As we were interested in the relative changes of the logarithmic flux (indicated as y), the flux units were disregarded. The horizontal line in all three panels indicates the average state of the source. Additionally, we highlight two flares, using gray bands in the upper panel, and we show zoomed-in versions in the lower panels. The first flaring event, labeled “extreme”, had a peak flux that exceeded by about a factor of two the maximum observed flux of Mrk 501. In contrast, the second flaring event, labeled “typical”, reached a peak flux comparable to the maximum observed flux amplitude in the variability of Mrk 501. With blue markers, we indicate the interpolation of this LC at time steps of one light crossing time, t c r = 0.55  days (in the observer’s frame), creating our Time Series (TS). Thus, we accurately captured changes in the parameter space of our source7.
In our parametric study of variability, therefore, we examined long- and short-term changes with the complete TS and the two flaring events, respectively. For each case, the respective interpolated TS started close to zero with an initial 10-step “burn-in” phase where no changes occurred in any parameter value. This allowed the system to accurately converge to its average state before starting deviations from it.

3.1. Particle Energy Injection Rate and Magnetic Field Strength

To motivate variations of the injected luminosity L e , p of the primary electron and proton populations or the magnetic field strength B of the emission region, we connected them to the observed flux variability in HE γ -rays probed by Fermi-LAT. In particular, we followed the principles outlined in the previous section and assumed a predominantly leptonic description of Mrk 501 based on the SSC model. In this description, we correlated the peak flux of each component to our desired parameters. In the δ -function approximation for the synchrotron and Inverse Compton (in the Thomson limit) emissivities, the synchrotron component at peak energy ϵ and the SSC component at peak energy ϵ s scaled with the model parameters as [59]
f ϵ s y n B 2 L e , f F e r m i f ϵ s S S C ( f e s y n ) 2 B 2 f F e r m i ( B L e ) 2 .
where the symbol f ϵ denoted the flux in ϵ F ( ϵ ) units. We then adopted the same scaling relation between the γ -ray flux and L p , assuming both populations were injected into the source via the same mechanism. Therefore, we could connect the previously derived average-state parameters of Mrk 501 listed in Table 1 and Table 2 to the logarithmic flux variations of Figure 5, as follows:
l o g ( L e , p ( t i ) ) = l o g ( L e , p ) + 1 2 y ( t i ) , or l o g ( B ( t i ) ) = l o g ( B ) + 1 2 y ( t i ) ,
with t i = i · t c r being the i-th time step of the time series. Both variations (when applied independently) failed to create flux changes that approximated the observed features of Mrk 501. We illustrate this via two arguments.
First, we present in Figure 6 the SED components for the extreme flare at the peak time of the injected TS. The average-state SED (Figure 3) from which the variations were initiated is highlighted with a black dotted line. The complete time evolution of the SEDs during this event is discussed in Appendix D, while the SED density profiles of the complete numerical runs for both parameter variations are provided in Appendix E. Furthermore, we provide histograms of the parameter variations in Appendix F. With the applied parameter variations, we could not recover the spectral hardening required to explain the most luminous X-ray flare ever observed [41]. Moreover, as demonstrated in Figure 6 for the peak of the flare (and in Appendix D for the entirety of the event), we found no distinct hadronic signatures in the SED of Mrk 501 throughout the 14.6 year-long period that we studied. All hadronic-related spectral components were hidden below the synchrotron and SSC emission of primary electrons, as illustrated in Figure 6. Moreover, the  π 0 bump was suppressed when the attenuation by the extragalactic background light (EBL) was taken into account.
Secondly, we computed the flux histograms of the complete numerical runs, comparing our time-variation models to observations in different energy bands. We used the LC from the Fermi LCR for the time series analysis and, thus, chose the Fermi band of 0.1–100 GeV for the HE γ -ray band. Having utilized Swift XRT and BAT data in our average-state analysis, we used their respective energy ranges of 0.2–10 keV and 14–195 keV to describe the soft and hard X-ray emissions. For the VHE range, we used 0.3–3 TeV to describe the region of high flux sensitivity of the upcoming Cherenkov Telescope Array Observatory [60,61]. Lastly, we chose the R-band (138–658 nm) to represent the low-energy emission in the optical. In particular, we utilized observations from the GASP program of the Whole Earth Blazar Telescope (WEBT) and data from the Tuorla observatory using the KVA telescope. The resulting LC had a total length of 4 years, a binning time of 1 day, and it has been adopted by ref. [45]. The shorter duration of this observational LC is a caveat that needs to be taken into consideration when interpreting our results. Specifically, the R-band observations are from a historically low-activity period in X-rays and γ -rays for Mrk 501. As a result, with fewer fluctuations between low- and high-states their variability remained confined within narrower flux ranges. For this reason, they might not have been representative of the LC’s behavior during high-activity periods or while studying long-term behavior.
We present in Figure 7 and Figure 8 the flux histograms of the full (14.6 yr) leptohadronic runs for L e , p and B variations, respectively. In both figures, the observational histograms (whenever available) for each energy band are overplotted for comparison. We also included the histogram for the synthetic LC used for the time variability. There is agreement in the Fermi band between the model and observations, as expected (the time series were created in a way that reproduced the observed GeV variability). However, we did not recover the expected variability in the other energy bands. The observational X-ray histograms were broader, while our models predicted narrower flux distributions. Changes in the injection power of radiating particles and/or magnetic field strength alone did not result in large spectral changes in X-rays, which are needed to produce a larger spread in the derived X-ray fluxes. On the contrary, the predicted flux distribution in the R-band was broader than the observed one. This was possibly due to the shorter duration of the R-band observations used for this plot. Including data from longer time periods might have resulted in a wider range of R-band fluxes, similar to what is observable in the other energy bands. If, however, the R-band histogram had remained narrow, this would suggest a variability mechanism operating in Mrk 501 that produces less variable optical fluxes than in the other ranges in the long term.

3.2. Power-Law Slope of Particle Distribution

Changes in the power-law slope of the radiating particle distribution are typically invoked to explain spectral changes (especially during flares) that cannot be attributed to cooling effects [26].
While we did not model the acceleration zone, we note that the power-law index of the injected particles into the radiation zone is determined by the acceleration and escape processes from the former zone. To motivate variations of the power law index in our study, we noted its relation to the acceleration and escape timescales [62]:
p e , p = 1 + t a c c t e s c .
Expressions for t a c c depended on the underlying acceleration mechanism. For example, in a Fermi acceleration scenario, we expected dependence on the magnetic field, either implicitly or explicitly through the gyroradius of a relativistic particle, r g m e , p γ c 2 / ( e B ) . In particular, for first- and second-order Fermi acceleration the dependencies of the acceleration timescales are t F I r g B 1 and t F I I r g / B 2 B 3 , respectively [63]. Assuming the magnetic field strength in the acceleration and radiation zones to be similar, we were able to parameterize the power-law index of the injected particle distributions as
p e , p ( t ) = 1 + p e , p B ( t ) B m ,
where denoted the parameter values of the time-average model. Therefore, the power-law index variations were tied to small amplitude changes in the magnetic field strength. Given that we did not model the acceleration physics, we explored the impact of different values of m. Therefore, utilizing Equation (5) and choosing a value of m we transformed the TS of Figure 5 into a time series of the power law index p e , p .
To demonstrate the impact of the parameter m we present in Figure 9 two scenarios. We compare mild ( m = 0.5 ) to extreme ( m = 4 ) power-law index variability for the extreme flare introduced in Figure 5. We note that in the former scenario, the values of p e , p vary between 2 and 2.5, while the range of the latter is 1.05–2.5. We contrast the leptohadronic SED (left column) to the leptonic SED (right column) to exemplify the clear hadronic impact in this scenario. For the scenario with m = 0.5 the spectral hardening is not enough to describe the hard spectrum of the extreme X-ray flare of 1997 (see lower panels). On the contrary, for the extreme slope variations with m = 4 the electron synchrotron spectrum becomes as hard as the observed one. Due to this spectral change, the number of low-energy seed photons used for IC scattering decreases compared to the average state (dotted line), leading to suppressed IC emission (from the primary electrons) at the HE and VHE bands. However, secondary electrons produced from photomeson and Bethe–Heitler interactions, as well as γ γ pair creation have a dominant contribution to the HE and VHE emissions, which exceeds by 5.6 and 6.7 times the maximum luminosity observed by Fermi in that range (0.1–100 GeV) for  m = 4 and m = 0.5 , respectively. Furthermore, when m = 4 we notice a prominent pion bump emerging at ϵ γ = 4 · 10 16  eV, due to the low intrinsic γ γ opacity at the peak time of the X-ray flare8. The hardening of the electron synchrotron spectrum during the flare reduces the number of low-energy photons ( ϵ t 10 5  eV) that are targets for γ γ pair creation. However, this hadronic spectral component would not be detectable due to the EBL attenuation (not included here). Still, an accompanying high-energy neutrino component is expected. We discuss our model predictions in the context of the IceCube observations in Section 5.
Motivated by the results above, we are sought slope variations that could lead to significant spectral hardening, but without producing extreme hadronic contributions at (V)HE. Therefore, we set m = 2 and reduced the injected proton luminosity L p by two orders of magnitude compared to the upper limit listed in Table 2. The choice of m = 2 was consistent with second-order Fermi acceleration ( t F I I r g / B 2 ) being at work in the acceleration region Mrk 501 [25] and with an escape timescale (from the acceleration region) being proportional to the particle gyroradius. Lastly, in the long-term numerical run, we imposed an upper limit of 2.6 on the power-law index of the injected distributions whenever Equation (7) predicted very soft power laws. This ensured that our analysis was not affected by pathological cases, e.g., the distributions were effectively mono-energetic. From this point on, we will refer to these simulations as the “modified power-law index variations”. Our results at the peak of the injected TS are shown for both flares in Figure 10 and as a density plot for the complete numerical run in Appendix E. Moreover, we present the histogram of the power-law index variability in Appendix F. By reducing the injected proton luminosity the emission of secondary electrons from Bethe–Heitler, photo-meson, and  γ γ pair creation interactions was reduced to being comparable to or below the primary electron IC emission. The notable exception is the photo-meson contribution in the GeV range during the extreme flaring event during the peak of the injected TS (compare the left columns of Figure 9 and Figure 10). A detailed discussion on the behavior of these three SED components will be presented in Section 6.
In Figure 11, we present the result for the flux histograms for the complete numerical run. We report discrepancies in the flux ranges between model predictions and observations at all energies except the Fermi band. To highlight this result, we integrated the 1997 flare observations in the XRT and BAT energy range, and we present the respective integrated flux values with a black dashed line in the appropriate histogram panels. We thus conclude that our model over-predicted the observed flux in the X-ray range even when accounting for the major X-ray flare of 1997. However, our analysis did not attempt to fit the observational flux histograms. Rather, our time variations explored the behavior of Mrk 501 in case of a future extreme flare. As such, small deviations from the 1997 flare, as seen in Figure 11, may not be indicative of disagreement between our model and observations, as discussed in Section 6. Furthermore, the distribution shapes deviated at HE. Compared to the observations, we noticed a peak at higher flux values when the power-law index was fixed at p = 2.6 , whereas the distribution peaked at lower flux values when the IC contribution of primary electrons was reduced compared to the average-state emission. This is exemplified in Figure 12, where we highlight the LCs in different energy bands as obtained from our model for the two individual flares of interest. We report a convoluted LC picture for the flaring states. This was expected as during periods of high activity (as defined by our injected TS, shown in black) the targets for IC scattering were reduced, leading to a decrease in the flux of the SSC component. Simultaneously, we see that the expected extreme flux increased in the X-ray band, as a result of the spectral hardening, and the resulting increased synchrotron emission. These two results can also be contrasted by comparing the upper to the lower panel of Figure 12. We see how the VHE flux initially “mimics” the behavior of the X-ray flux corresponding to the initial increase of high-energy electrons (due to the hardening of the injected power law). The sudden decrease around the peak of the injected TS corresponded to the spectral hardening in the low-energy synchrotron region and the resulting decrease of available SSC targets.
In Figure 13, we present the flux–flux diagrams, comparing (a) the Fermi flux to the flux in the VHE range, (b) the harder X-ray range probed by the BAT, and (c) the softer X-ray range probed by XRT; furthermore, in (d) we compare the BAT flux to the R-band flux, comparing hard X-ray to low-energy photons. We present the variability for the complete leptohadronic run (blue circles) and the results for the extreme (yellow stars) and typical flare (purple squares). Each of these diagrams consists of a linear segment (at low fluxes) and a parabolic-like shape (at higher fluxes). In each panel, the parabola section to the right of its vertex corresponds to flaring states, as a result of reduced seed photons for SSC emission. Increasing the available L p is a method of counteracting this decrease. Furthermore, we suggest that the parabola occurs during joint magnetic field and power-law index variations while the linear segment corresponds to the magnetic field variations when p e , p = 2.6 (upper bound on slope).

4. Predictions for CTAO

The Cherenkov Telescope Array Observatory (CTAO) will be the new-generation ground-based observatory for γ -ray astronomy [61]. It will cover a wide range of energies, from 20 GeV to 300 TeV, and it will offer better sensitivity and spectral resolution than the VHE instruments currently in operation [61]. With an expected start of operations in 2025, model predictions in the spectral and time domain for CTAO are timely and informative.
For this section, we were interested in creating a simulated LC for the CTAO North array and determining the impact of the chosen exposure time. Furthermore, we have created simulated spectra to test for statistically significant differences between a leptohadronic and leptonic scenario. As our test case, we chose the modified power-law index variations and focused on the extreme flare, as it exhibited strong anti-correlation behavior during the peak of the injected TS (Figure 12). Additionally, at peak time, the SED exhibited a dominant hadronic component in the HE range (synchrotron from photomeson secondaries; Figure 10).
For this analysis, we used the standard procedures of Gammapy9 [64], a community-developed open-source Python library and core package for CTA. Thus, first, we assigned the coordinates of Mrk 501 as the source position in the sky and set an 11-degree extraction region centered on this location. To this region, we assigned our leptohadronic model (with modified power-law index variations) multiplied by e τ E B L ( z , E ) . Here, τ E B L ( z , E ) was the optical depth to γ γ pair production on EBL photons and was adopted from [65]). We chose as the energy axis the range 0.05–10 TeV, in order to study the source’s behavior in both the HE and VHE bands. Then, we simulated the γ -ray spectra for each time step of the extreme flare of our variation model, using the available instrument response functions of the CTAO northern array10. Finally, we extracted the predicted LC of this event, and we present our results in Figure 14 for exposure times of 30 min, 1 h, and 5 h. For visual clarity purposes, we shifted the 1 h and 5 h points by 0.25 and 0.5 days to the right, respectively. Highlighted with gray dashed lines are two indicative time values (t = 11, 43 days) that are used in our subsequent discussion of the spectra. We conclude that all three exposure times yielded comparable results in the 0.05–10 TeV range, while increasing it only resulted in marginally lower error values. Thus, within our model framework, low exposure times are sufficient to observe Mrk 501 during a flaring event of brightness similar to our extreme flare.
To investigate differences between leptonic and leptohadronic scenarios within our model description, we fixed the exposure time to 5 h and extracted for both scenarios the LCs in the “lower” energy range of 0.05–1 TeV and the “higher” range of 1–10 TeV. We compare our findings in Figure 15. We found no significant discrepancy between the two models in the higher energy range, as expected; the 1–10 TeV range was dominated by the SSC emission of primary electrons, as shown in Figure 10. On the contrary, the LCs of the leptonic and leptohadronic models in the 0.05–1 TeV range exhibited deviations in the decay part of the flare (see, e.g., dip at 43 days).
To further exemplify and study this hadronic impact, we generate spectra for the complete energy range of 0.05–10 TeV, for a 5 h exposure, and we fitted the simulated data with an empirical model of a Log Parabola11 (LP):
ϕ L P ( E ) = ϕ 0 E E 0 α β log E E 0 .
This model had a reference energy frozen at E 0 = 1 TeV during the fitting process and an amplitude ϕ 0 in units of (cm−2 s−1 TeV−1). The fitting of our simulated spectral data to the empirical models utilized cash statistics [66] integrated within the standardized procedures of gammapy. This approach involved performing a log-likelihood minimization, returning a fit statistic (−2LogL), and calculating the residuals (residuals = (data-model)/model). An example of simulated spectra computed for both leptonic and leptohadronic models at t = 11 and 43 days is presented in Figure 16 (left and right panels, respectively), while their respective best-fit values are summarized in Table 4. No difference between the two models is observable during the onset of the plateau at t = 11 days. In contrast, at the dip of the simulated LC at t = 43 days, a clear hadronic contribution is present at energies below 1 TeV. This was expected from our previous LC analysis (Figure 15), and it is confirmed by visual inspection of the left panel of Figure 10. This SED highlights the dominant hadronic component in the HE region at lower fluxes, and the negligible contribution to the VHE region (which is dominated by primary electron IC emission at all times).
To quantify the statistical significance of this deviation between the two models, we created 1000 leptonic and leptohadronic synthetic CTA spectra for t = 43  days, and t = 11  days recording their best-fit log-parabola parameters. The resulting parameter distributions are presented as histograms and contour plots in Figure 17. Based on the contour plots, for  t = 43 days the model parameters were different at least at a 3 σ level, indicating a clear discrepancy between the leptonic and hadronic models. On the contrary, the distributions for t = 11  days overlapped significantly and, thus, the two models must be assumed indistinguishable.
The pure leptonic model resulted in a larger range of observed flux (between day 11 and 43) and had a more evident change in spectral shape (see the range of α and β parameters), whereas there was less evident spectral change in the leptohadronic model. We note, however, that from γ -ray spectral modeling alone it was challenging to distinguish between models and to secure a hadronic contribution. The reason is that the hadronic emission in our model does not have a unique spectral feature, as in the SSC+ π 0 model introduced by [14]. To that end, simultaneous X-ray and HE/VHE γ -ray observations were needed.

5. Predictions for IceCube

IceCube is a cubic-km neutrino detector located at the South Pole. In 2013, IceCube made a breakthrough with the discovery of a very high energy (TeV-PeV) flux of cosmic neutrinos [67], which led to the detection of the first likely source of high-energy neutrinos, the blazar TXS 0506+056 [68]. So far, no other blazar has been associated with high-energy neutrinos at a confidence level greater than 3 σ . The non-detection of neutrinos from the nearby and electromagnetically bright blazar, Mrk 501, can place strong constraints on leptohadronic models. For this section, therefore, we computed the number of events expected from Mrk 501 during the IceCube lifetime.
Using the neutrino energy flux at time t, F ν + ν ¯ ( ε ν , t ) , where ε ν was the observed neutrino energy, computed during the modified power-law index variations (see Section 3.2), we could make predictions regarding the number of neutrino events detected by IceCube and assess the potential presence of hadronic signatures expected from this model. During these modifications in the slope of the particles p e , p (see Equation (7) for m = 2 ) and the magnetic field B ( t ) the power injected into the electron and proton distributions remained constant. Therefore, the applied changes altered the particle number per energy bin.
The expected number of muon plus antimuon neutrinos from Mrk 501 could be calculated as
N ν μ + ν ¯ μ = 1 3 t ini t end d t E ν , min E ν , max d ε ν A eff ( ε ν , D E C ) F ν + ν ¯ ( ε ν , t ) ε ν ,
where we assumed vacuum neutrino mixing and used 1/3 to convert the all-flavour flux to muon neutrino flux. Moreover, t ini and t end defined the duration of the simulated LC, and A eff ( ε ν , D E C ) was the energy-dependent and declination-dependent point-source effective area of IceCube [69,70].
To account for the varying IceCube configurations over different operation seasons (see Table 5), we initialized our calculation by assuming that the last day of the Light Curve shown in Figure 5 corresponded to 1 January 2024. In the bottom panel of Figure 18, we present the expected muon neutrino number as a function of time (see blue solid line). When considering only the most recent IceCube configuration (IC86 II), the final neutrino number showed negligible variation ∼3% (see the magenta dashed line in the bottom panel of Figure 18). Our analysis indicates that the neutrino signal over this approximately 14.5-year period was less than one event ( N ν μ + ν ¯ μ 0.2 ). In this synthetic LC, we set the proton luminosity L p lower by two orders of magnitude compared to the upper limit listed in Table 2. An increase in the proton luminosity to this value can result in about ∼20 high-energy muon neutrinos in 14.5 years (or 1.4 neutrinos per year, on average). Mrk 501 belongs to a list of predefined γ -ray sources used in IceCube point-source searches (e.g., [69,71]). By analyzing 10 years of IceCube data [69], the Collaboration obtained 10.3 best-fit signal events from Mrk 501 (see Table III in ref. [69]), which, however, did not secure a statistically significant detection. Instead, they corresponded to an upper limit of one neutrino per year. The model predictions shown in Figure 18 are consistent with the IceCube upper limit (see also [45]).
The top panel of Figure 18 shows the ratio between the flux at 100 MeV (Fermi-LAT) and the flux at 14–195 keV (Swift-BAT) as a function of time. The power-law index of the two distributions p e , p is plotted against time in the middle panel. We found that during the hardening of the distributions, the ratio L 100 MeV / L 14 195 keV decreased while the number of neutrinos was increasing (see, for example, the time window from 10 to 11 years). The emission at 14–195 keV primarily arose from primary electrons, and the hardening of their distribution resulted in an increased flux in this energy band. Additionally, the hardening reduced the number of low-energy target photons available for upscattering to the sub-GeV band, as described in Section 3.2. During these hardening phases, this band was predominantly influenced by the synchrotron emission of secondaries produced through photomeson interactions (see, in the left panel of Figure 3, the sub-GeV band). The increase in high-energy protons during these periods also enhanced the source’s neutrino production, a trend evident when comparing the middle and bottom panels of Figure 18.

6. Summary and Discussion

In this work, we derived a leptonic-average-state description of the HSP blazar Mrk 501, using publicly available multi-wavelength archival observations, and utilizing data from the BAT 157-Month Hard X-ray Survey, the 18.5 yr Swift-XRT LC, and the 14-year Fermi-LAT 4FGL catalog. To derive the parameters of the baseline leptonic model, we used the MCMC sampler emcee and the numerical code LeHaMoC. In our subsequent arguments for the proton population parameters, we assumed p p = p e , γ p m i n / m a x = γ e m i n / m a x , and we searched for a conservative upper limit on the injected proton luminosity L p . Using celerite2, we analyzed the 0.1–100 GeV LC of Mrk 501 and created TS to describe key parameter variations ( L e , p , B, and  p e , p ). Based on these, we modeled the long-term evolution of the multi-wavelength jet emission of Mrk 501 in the context of a leptohadronic model, contrasting our findings, in the form of flux histograms, to the available observations. Lastly, we searched for detectable hadronic signatures on simulated VHE γ -ray spectra for CTAO in the investigated power-law index variations and discussed our model in light of the neutrino non-detection of Mrk 501 by IceCube.
To test the validity of our initial average-state parameters, we calculated the average jet power of Mrk 501 in the context of a leptohadronic model as P j π R 2 δ 2 β c ( u e + u p + u B + u r a d ) π R 2 δ 2 β c u p [72], where Γ δ , β = 1 1 / Γ 2 1 , u e , p were the energy densities of the relativistic particle populations, u B was the magnetic energy density, and  u r a d was the radiation energy density. In the above, all energy densities were measured in the jet co-moving frame, while the relativistic proton energy density, defined as u p = 3 L p / ( 4 π R 2 c ) , was the dominant component. In particular, for the derived parameters listed in Table 1 and Table 2 we calculated P j 3 · 10 49  erg/s, or  P j 3 · 10 47  erg/s for L p = 10 3.61 L e used in the modified power-law index variations of Section 3.2. The latter was comparable to the Eddington luminosity, L Edd 1.3 · 10 47 erg/s, of a black hole with a mass of M B H 10 9 M [73] and, as such, was favored over the preceding jet power estimate. We also noted that the inferred jet powers were comparable to those found in previous studies of blazar leptohadronic emission models [74,75].
To limit the number of free parameters in the leptohadronic model we assumed identical values for the power-law slopes and extent of the electron and proton distributions. An alternative choice for γ p m i n , m a x was to relate the the minimum and maximum energies of the particle populations as E p m i n / m a x = E e m i n / m a x or γ p m i n / m a x = ( m e / m p ) γ e m i n / m a x . This choice resulted in γ p m i n = 1 and γ p m a x = 10 4.12 , where the lower value was limited by the minimum physically acceptable Lorentz factor. Repeating the method outlined in Section 2.3, we found the limit on the injected proton luminosity to be L p 10 9 L e . Consequently, the upper bound on the jet power was P j = 7.4 · 10 52  erg/s, which was non-physical, based on energetic arguments for the central engine. We note, however, that combining a reduced proton luminosity L p = 10 3.61 L e (as in Section 3.2) with this alternative choice for the proton Lorentz factors resulted in a model with jet power comparable to the Eddington luminosity that can create a pion bump at VHE during the modified power-law index variations discussed in this work. This model bore many similarities to the SSC+ π 0 model presented recently by ref. [14] to explain a narrow spectral feature in the flaring VHE spectrum of Mrk 501 [19]. The implications of relaxing the assumption p e = p p would strongly depend on the specifics of our choice. For example, if  p e > 2 (as found by the leptonic fit, see Table 1) and p p < 2 (hard power law), then the limit obtained for L p would be lower than the one listed in Table 2. For a fixed target photon field (which is produced by the primary electrons), the energy injected into secondaries produced via proton–photon interactions increases with decreasing p p , because more power is carried by protons of the distributions with higher energies that are relevant for Bethe–Heitler and photomeson interactions (see, e.g., Figure 12 in ref. [76] and Appendix B in ref. [77]).
To initiate parameter variations from the average-state description we used a synthetic LC of Mrk 501, motivated by its observational behavior in the Fermi energy range of 0.1–100 GeV (as this is the best observationally sampled LC over year-long timescales). Specifically, we conducted a time series analysis with celerite2 on data from the Fermi LCR with a mean cadence of 7 days. A similar methodology has been recently employed by ref. [78] to study the flux variability in single parameter variations of a one-zone SSC model of blazar sources. In this study, the authors used the Emmanoulopoulos method [79], a combination of the Timmer and Koenig procedure [80] and the Schreiber and Schmitz Fourier transform algorithm [81], which produces a TS with the same PSD and PDF as the input TS. In this approach, the input data must be an equally spaced TS. However, astrophysical observations are rarely without gaps and, thus, the input TS must be either artificially connected at the gaps or interpolated between the gap points. The first can result in stiff jumps in the synthetic LCs [82], while the latter has been shown to contaminate the power spectrum of the TS [83]. In contrast, celerite is a Gaussian process method that does not have the above restrictions. Thus, we used the provided Fermi LC without further modifications that could insert uncertainties into our problem description. Furthermore, this difference would allow for investigations similar to ours using, instead of γ -ray LCs, X-ray LCs that are known to have gaps [84,85]. Observations, such as from the Swift-XRT or Swift-BAT instruments (see Appendix A), can be modeled by celerite2 and used to create synthetic LCs motivating key parameter variations. Additionally, unlike in Fourier transform methods, it is possible to increase the model complexity of a Gaussian process, e.g., by increasing the number of SHOs in celerite2, and to sample the resulting parameter space with nested sampling methods, such as Ultranest12 [86]. This could allow for the search of secondary features in the PSD and degeneracies in the parameter space. Cross-correlations between such studies and ours could provide insights and constraints into the underlying variability mechanisms in Mrk 501. However, we note that celerite2 is restrictive in the models it can use. For example, the SHO model of Equation (3) cannot be evaluated with celerite for the choice of Q = 1 / 2 , corresponding to a critically damped oscillator, as the elements of the covariance matrix of the Gaussian log-likelihood become infinite (Appendix C for log-likelihood; see [51] for more details). As such, this can impose restrictions on the celerite method when exploring the parameter space, which must be considered.
Having created a synthetic LC for Mrk 501, we instigated parameter variations from their previously derived average-state values. In particular, we studied the effect of variations in electron and proton injection luminosity, the magnetic field strength, and the power-law index of both populations motivated by the magnetic field strength variations. We found that none of those scenarios could explain the complete observational behavior of Mrk 501. In particular, the first two scenarios failed to produce the broad observational X-ray flux distributions of the source and the hardness of the X-ray spectrum during bright X-ray flares, while they over-predicted the variability in the optical R-band. In the power-law index variations, we investigated the impact of the parameter m, which modeled the hardness of the power law in terms of the magnetic field strength, and we showed the hadronic impact, using the upper limit of the injected proton luminosity L p up . lim . . Subsequently, we modeled the predicted long-term behavior of our source, using a physically motivated choice of m = 2 , a conservative value for the injected proton luminosity of L p = 10 2 L p up . lim . , and an upper limit on p e , p (≤2.6) to ensure we were not affected by pathological cases. We showed that this model can produce bright hard X-ray flares with a VHE counterpart, with only moderate changes in the magnetic field strength (see Figure 8). A hadronic component was needed, to avoid a large depletion of the GeV flux at the peak of the X-ray/VHE flare due to the reduction of the primary SSC emission.
When considering the predictions of the power-law index variation scenario on year-long timescales, we found that the obtained flux distributions in the XRT and BAT bands extended to higher values than the observed ones (see Figure 11). This was related to the presence of the extreme flare, and similar ones to it, in the synthetic LC we used to construct our TS (see Figure 5). The extreme flare exceeded the maximum-ever-observed flux in the Fermi-LAT energy range by a factor of 2. While this could be used as a predictive tool for future flaring activity of Mrk 501, caution must be exercised when making conclusive statements comparing the upper X-ray flux ranges between observations and the model in Figure 11. A more meaningful comparison against the observational flux histograms would require the creation of a sample of synthetic LCs and repetition of our analysis for each entry of this sample. This approach would mitigate the effect of having a single synthetic LC that might be characterized by extreme flux occurrences. The impact of the extreme flare on the R-band flux histogram of our model was its extension to lower values; the extreme flare translated to a hardening of the injected electron distribution and, hence, a hardening of the synchrotron spectrum (see, e.g., Figure 10 and Figure A7). Moreover, the observational flux distribution was clustered to values at the upper end of our distribution (see Figure 11), even though these observations corresponded to a historically low activity period of Mrk 501 [45] and represented the lowest recorded flux values of Mrk 501 in this range. As discussed in the next paragraph, this discrepancy may suggest that the optical emission originated from a different region than the one producing the X-ray/VHE flares.
One-zone models are often used to describe the broadband blazar emission, especially when there is evidence of correlated flux variability between energy bands in the low- and high-energy SED humps. Given the small number of free parameters, one-zone models can be constrained by the data (in the spectral and time domain). In several cases, emission from a single zone (or a single particle population) is not sufficient to explain the whole SED [15,31,45,47]. In particular, one-zone models for Mrk 501 fail to capture the observations below the optical band when modeling simultaneous X-ray and γ -ray data. Another leading argument for a second emission zone is the inability to explain the entirety of the GeV spectrum (in the LAT band) with a one-zone SSC description. Specifically, in ref. [47] the authors utilized hard power laws for the electron distribution, to model the simultaneous X-ray and (V)HE emissions of the extreme flare of Mrk 501 in 1997. Simultaneously, they argued in favor of a continuous jet model (as well as the contribution of the host galaxy) to account for the emission from radio frequencies up to ultraviolet. In that work, an SSC model was sufficient to explain the entirety of the γ -ray spectrum, a fact attributed, in retrospect, to the poor quality of the available instruments of the time (i.e., comparison between EGRET and Fermi-LAT). In more recent studies, such as ref. [15], the authors have proposed two distinct regions to explain the flux states of Mrk 501 during 2011. They attributed the X-ray and VHE emission to a compact “inner” blob, while the larger “outer” blob explained the radio and GeV observations. Similarly, in the work of ref. [31] the authors modeled Mrk 501 while it exhibited extreme HSP behavior in 2012. In addition to the region responsible for the X-ray and VHE data, a second slower emission zone with variability timescales of weeks was suggested, to explain the optical and GeV emission. Lastly, ref. [45] investigated the behavior of our source during a historically low period of X-ray and γ -ray activity. They reported the larger second region in the two-zone description to correspond to the baseline emission of Mrk 501. They also found the size of this region to be in agreement with millimeter-VLBI observations. We report that our work arrived at similar conclusions for the radio emission and the GeV component. In our average-state modeling, we did not include measurements below 300 GHz. Moreover, the model SED of the most probable parameters could not capture the entirety of the Fermi observations (see the high-energy flux points of the 4FGL spectrum in Figure 3). In our framework, precise modeling of the acceleration region could provide a physically motivated, second emission zone (similar to the outer/inner blob scenario) while simultaneously constraining the variations imposed on the compact zone. Furthermore, to probe whether the emission from the historically low activity period is attributable to one zone or two components of highly variable and baseline emission one could repeat our methodology and initiate variations from the parameters of a leptohadronic model describing this region. Lastly, we suggest that an alternative to the above could be simultaneous uncoupled variations of the primary electron and proton parameters in the one-zone leptohadronic model. A scenario of constant electron and variable proton parameters could lead to a description similar to the two-zone SSC models.
In the leptohadronic models we explored, hadronic contributions to the emitted electromagnetic spectrum were sub-dominant by construction, except during hard X-ray flares (see Section 3.2). Thanks to the LeHaMoC code we were able to display all hadronic-related spectral components of the SED separately, despite them being most of the times hidden below the primary leptonic emission (see, e.g., Figure 3 and Figure 6). This feature of the code allowed us to explore the temporal evolution of the various spectral components13 and to compare them against the primary leptonic emission components. For the Bethe–Heitler component, the secondary pair distribution extends from γ e γ p m i n up to γ e 10 11 , with the shape of the distribution at injection depending on both the shape of the proton power-law index p p and the photon target spectrum, as shown in ref. [87] (see, for example, Figure 6 therein). As a result, the synchrotron spectra of primary and secondary pairs from the Bethe–Heitler process are generally different, as exemplified in Figure 6, Figure 9, and Figure 10. Additionally, we note how the secondary pair emission from photo-meson and γ γ absorption interactions are of similar shape and peak at similar energies in the γ -ray range. This can be understood by examining the energies at which the secondary electrons are injected into the source. During neutral pion production and subsequent decay, the energy transfer from the parent proton of energy ϵ p to each high-energy photon is given by ϵ γ h i g h = 1 2 K p γ ϵ p [88]. These high-energy photons interact with the low-energy emission ( ϵ γ l o w ) close to the threshold of the γ γ absorption process, creating electron–positron pairs with total energy approximately equal to the energy of the photon from the π 0 decay, namely, 2 ϵ e = ϵ γ h i g h [88]. Thus, secondary pairs produced through the γ γ process are injected with typical energies of ϵ e = 1 4 K p γ ϵ p , which corresponds to the energy of secondary electrons from charged pion decays (in photo-meson interactions) [88]. Given that both channels inject secondary electrons with similar energies, their synchrotron spectra will also peak at similar energies. Lastly, additional features in the γ γ pair synchrotron spectrum below its peak energy (e.g., see bump at X-rays) are attributable to γ γ attenuation of SSC photons from primary electrons. This is also highlighted by their presence in the leptonic SSC models of Figure 9.
To further investigate the effect of hadronic interactions on the SEDs, we simulated results for the upcoming CTAO creating LCs and spectra for our power-law index variation model. Subsequently, we investigated the impact of exposure time in the former and the significance of the hadronic component in the latter. We found that within our model description, short exposure times (i.e., 30 min) are sufficient to observe Mrk 501 while increasing the livetime (e.g., 5 h) only marginally reduces the resulting errors. Moreover, we investigated the hadronic impact on the simulated spectra and found a statistically significant contribution in the GeV range when a decrease of the primary IC emission is present (i.e., a reduction of seed photons). To test the statistical significance of our findings, we generated 1000 leptonic and leptohadronic spectra. For each data set, we fitted the empirical model of a log-parabola and recorded the resulting best-fit parameters. Contrasting these results for the two indicative periods of t = 11 and t = 43 days revealed a clear hadronic impact during the latter, where the contour plots between the two scenarios did not overlap at a 3- σ level. Simultaneously, in the former case of t = 11 days, both histograms and contour plots featured near-identical values and, thus, the two datasets must be treated as indistinguishable.
Lastly, we computed the time-dependent neutrino energy flux F ν + ν ¯ ( ε ν , t ) during variations into the power-law index, to predict the number of muon neutrino events detected by IceCube and assess the potential hadronic signatures from Mrk 501. While we maintained constant luminosity for electrons and protons we altered the spectral shape of the particle distributions. During particle distribution hardening, we observed an increase in neutrino production and a decrease in the flux ratio L 100 MeV / L 14 195 keV shown in the top panel of Figure 18, primarily due to the increased flux at 14–195 keV from primary electrons and reduced low-energy targets for upscattering to the sub-GeV band, which was also the band at which we could detect a hadronic component in the SED during flaring episodes. Our analysis, accounting for varying IceCube configurations, indicated that the neutrino signal over 14.5 years was less than one event ( N ν μ + ν ¯ μ 0.2 ), aligning with the non-detection of high-energy neutrinos from Mrk 501. This result was derived assuming a proton luminosity L p two orders of magnitude lower than the upper limit in Table 2, which was also discussed in Section 3.1. A higher proton luminosity would result in a neutrino expectation in tension with IceCube observations. This was also exemplified in ref. [85], where the neutrino expectation from blazars, including Mrk 501, was computed under the assumption of hadronic X-ray flares [89]. In this scenario, the X-ray flares were powered by proton synchrotron radiation, and neutrinos were produced through photomeson interactions between protons with their synchrotron X-ray photons. The conditions in the flaring region were such as to maximize the neutrino production rate. The 0.5–10 keV fluence of each flare in the X-ray LC of Mrk 501 was used as a proxy for the all-flavour neutrino fluence, and the rate of neutrinos was found to be ( 1990 ± 80 ) × 10 4 yr 1 (the prediction for 14.5 years was 2.9 ± 0.1 muon neutrino events). These results illustrate how neutrino observations can be used to constrain the relativistic hadronic component of blazar jets.

7. Conclusions

We explored the predictions of a leptohadronic model for blazar emission, searching for hadronic signatures in the spectral and time domains, and using Mrk 501 as a test case. By utilizing a synthetic LC motivated by Fermi-LAT observations in the 0.1–100 GeV range over 14.5 years, we created time series for the particle injection luminosity, magnetic field strength, and power-law index of the injected particle distributions. None of the individual scenarios could fully capture the observed variability across the electromagnetic spectrum. However, power-law index variations that are tied to moderate changes in the magnetic field strength of the emitting region may naturally lead to hard X-ray flares with VHE counterparts, as with those observed in Mrk 501. Our analysis revealed that secondary synchrotron emission from hadronic interactions can dominate the sub-TeV energy band during hard X-ray flares. This would have a measurable difference when compared to the primary SSC emission below 1 TeV. Coordinated observations of CTAO and X-ray satellites during X-ray/VHE flares can provide valuable insights into the origin of the HE and VHE γ -ray emissions of HSP blazars.   

Author Contributions

Conceptualization, M.P.; data curation, M.C. and G.V.; funding acquisition, M.P.; investigation, M.C.; methodology, M.C. and M.P.; resources, M.C.; software, M.C. and S.I.S.; supervision, M.P.; visualization, M.C.; writing—original draft, M.C.; writing—review and editing, S.I.S., M.P. and G.V. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the “2nd call for H.F.R.I. Research Projects to support Faculty members and Researchers” through the project UNTRAPHOB (Project ID 3013). G.V. also acknowledges support by the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the “3rd Call for H.F.R.I. Research Projects to support Postdoctoral Researchers” through the project ASTRAPE (Project ID 7802).

Data Availability Statement

SED data are available through the SED website (https://tools.ssdc.asi.it/SED/, accessed on 20 September 2024) maintained by the Space Science Data Center (SSDC, https://www.ssdc.asi.it/, accessed on 20 September 2024). X-ray and gamma-ray data from Swift and Fermi are publicly available through The High-Energy Astrophysics Science Archive Research Center (HEASARC, https://heasarc.gsfc.nasa.gov/, accessed on 20 September 2024). The code LeHaMoC is available via a GitHub repository at https://github.com/mariapetro/LeHaMoC, accessed on 20 September 2024. Scripts for reproducing results and simulated datasets are available upon reasonable request to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:
ACFAutocorrelation Function
BATBurst Alert Telescope
CTACherenkov Telescope Array
EBLExtragalactic Background Light
HEHigh Energies
HSPHigh-Synchrotron Peaked
ICInverse Compton
KSKolmogorov–Smirnov
LATLarge Area Telescope
LCLight Curve
LCRLight Curve Repository
LPLog Parabola
MCMCMarkov Chain Monte Carlo
PCPhoton Counting
PSDPower Spectral Density
SEDSpectral Energy Distribution
SHOStochastically driven damped simple Harmonic Oscillator
SSCSynchrotron-Self Compton
TSTime Series
VHEVery High Energies
WEBTWhole-Earth Blazar Telescope
WTWindowed Timing
XRTX-ray Telescope

Appendix A

Figure A1. Multi-wavelength Light Curves of Mrk 501 (from top to bottom): (first row) 0.1–100 GeV from the Fermi LCR; (second row) 157-month Swift BAT rate in the 14–195 keV range; (third row) Swift XRT counts in the 0.2–10 keV range for the Windowed Timing (WT) and Photon Counting (PC) mode; (fourth row) spectral flux for the R-Band from observations from the GASP program of the Whole-Earth Blazar Telescope (WEBT) and data from the Tuorla Observatory using the KVA telescope. Time units of HJD can be compared to MJD for large timescales.
Figure A1. Multi-wavelength Light Curves of Mrk 501 (from top to bottom): (first row) 0.1–100 GeV from the Fermi LCR; (second row) 157-month Swift BAT rate in the 14–195 keV range; (third row) Swift XRT counts in the 0.2–10 keV range for the Windowed Timing (WT) and Photon Counting (PC) mode; (fourth row) spectral flux for the R-Band from observations from the GASP program of the Whole-Earth Blazar Telescope (WEBT) and data from the Tuorla Observatory using the KVA telescope. Time units of HJD can be compared to MJD for large timescales.
Universe 10 00392 g0a1

Appendix B

In this appendix, we show additional results from the modeling of the average SED of Mrk 501. A corner plot showing the posterior distributions of the leptonic model parameters we obtained is presented in Figure A2 alongside an error value “log(f)” that we have added to the standard deviation of the Gaussian likelihood as [82,90]
σ 2 = σ d a t a 2 + f 2
Thus, we account for additional sources of uncertainty, such as use of data that are not averaged over the same time period, that cannot be captured by the statistical uncertainties of the observations/measurements. Most parameters are well-constrained, except for the blob radius. Using even larger upper bounds in the prior distribution of this parameter (e.g., l o g ( R m a x ) = 19 ) yields solutions that describe slower and larger emitting regions with weaker magnetic field strengths, with median values (of the parameter posteriors) of R = 10 18.3  cm, B = 0.003  G, and  δ = 7.4 . These parameters lead to variability on timescales of months (in the observer’s frame), t v a r ( 1 + z ) R / ( δ c ) 100  days, and, as such, are not a physically meaningful choice when modeling flaring events on day-long timescales. The parameters of Figure A2 are self-consistent in this regard, since a fiducial variability timescale of t v a r = 1  day imposes a maximum radius of R m a x = 10 16.7  cm, above our imposed limit on the posterior.
Figure A2. Posterior distributions obtained from a leptonic fit to the average SED of Mrk 501 with emcee and LeHaMoC.
Figure A2. Posterior distributions obtained from a leptonic fit to the average SED of Mrk 501 with emcee and LeHaMoC.
Universe 10 00392 g0a2

Appendix C

To assess the goodness of our celerite fit in Section 2.4 we employed several statistical tests, which are summarized in Figure A3 and Table A1.
  • In the context of Gaussian fitting methods, the covariance matrix K α is used in the generalized N-dimensional Gaussian log-likelihood for datasets y and X, parametrized by θ ,   α . This likelihood is given by
    ln L ( θ , α ) = ln p ( y | X , θ , α ) = 1 2 r θ T K α 1 r θ 1 2 ln ( det K α ) N 2 l n ( 2 π ) ,
    with
    r θ = x 1 μ θ x n μ θ .
    Its elements, [ K α ] n m = k α ( x n , x m ) , are called kernel functions and are calculated by the celerite method.
    We compute the standardized residuals of the Gaussian fit, which must correspond to white noise. By connecting the celerite process to moving averages we express the residuals in terms of the upper triangular Cholesky factorization of the covariance matrix, as provided by the celerite process, and the zero-meaned observational data as
    y * = c h o l ( K α ) T w .
    Fitting a normal distribution to the resulting histogram via the chi-square method, the desired result should be a Gaussian with mean μ = 0 , a standard deviation of σ = 1 , and a χ 2 value of 1.
  • We calculated the deviation of our results from a normal distribution via the Kolmogorov–Smirnov (KS) test. The null-hypothesis in this test is that the results do not originate from the same distribution. This is quantified by p-values, which describe the likelihood of an event. Values smaller than 0.05 imply that the null hypothesis is true only 5% of the time. Therefore, p values greater than p > 0.05 confirm the hypothesis that our residuals stem from a normal distribution.
  • We showed the residuals’ Autocorrelation Function (ACF) and the 95% confidence limits of white noise, proving that almost the entirety of the ACF resided within that limit. Thus, our model selection appropriately captured the correlation behavior.
  • We plotted the PSDs of the MCMC process and calculated the slope above the breakpoint. Values of p = 2 and p = 0 corresponded to red noise or Brownian motion and white noise, respectively. Such combinations have been previously observed in the study of the time variability of blazars with an empirical break frequency of around 150 days for OJ 287 [91], 25 days for 3C 66A [92], and 43 days for PKS2155-304 [92]. We report a value of t b r = 620 d for Mrk 501.
Figure A3. (Left): Autocorrelation function of the standardized residuals. Indicated with a gray band is the 95% limit of white noise. (Right): Histogram of the standardized residuals (black) and Gaussian fit (blue).
Figure A3. (Left): Autocorrelation function of the standardized residuals. Indicated with a gray band is the 95% limit of white noise. (Right): Histogram of the standardized residuals (black) and Gaussian fit (blue).
Universe 10 00392 g0a3
Table A1. Statistical tests of goodness of fit. The mean, standard deviation, and χ 2 value of the Gaussian fit on the standardized residuals. The p-value of the KS test. The break frequency and the slope of the PSD above it.
Table A1. Statistical tests of goodness of fit. The mean, standard deviation, and χ 2 value of the Gaussian fit on the standardized residuals. The p-value of the KS test. The break frequency and the slope of the PSD above it.
Source Name μ ( 10 2 ) σ χ 2 KS p-ValuePSD Slope t br [Days]
Mrk 5012.950.980.960.2 1 . 90 0.03 + 0.02 620 9 + 12

Appendix D

In this Appendix, we derive the SED time evolution for the “extreme flare” in the variation of L e , p ,   B ,   p e , p and for the “typical flare” in the variations of p e , p . These events were identified in the long-term TS and introduced in Section 3. For the power-law index variability, specifically, we focused on the choices of m = 2 and L p = 10 3.61 L e (see Section 3.2 for details). In Figure A4, we present the components of the SEDs in the leptohadronic scenario at each time step and for a period of about 15 years. For visual clarity, we present only the emission components of Bethe–Heitler, photo-meson, and γ γ pair-creation secondary electrons. The initial average-state SED is highlighted with a black dotted line. The resulting figure resembles a density plot differentiating between regions of frequent and rare occurrences. Links to movies with all components of each scenario are provided in the figure caption.
Figure A4. SED time evolution of Mrk 501 for the following scenarios (from top to bottom): (upper left) p e , p variations in the extreme flaring event (YouTube); (upper right) p e , p variations in the typical flaring event (YouTube); (lower left) B variations in the extreme flaring event (YouTube); (lower right) L e , p variations in the extreme flaring event (YouTube).
Figure A4. SED time evolution of Mrk 501 for the following scenarios (from top to bottom): (upper left) p e , p variations in the extreme flaring event (YouTube); (upper right) p e , p variations in the typical flaring event (YouTube); (lower left) B variations in the extreme flaring event (YouTube); (lower right) L e , p variations in the extreme flaring event (YouTube).
Universe 10 00392 g0a4

Appendix E

To represent the density profiles of the complete numerical runs, we interpolated each SED curve for 10,000 points and calculated the 2D histograms. To describe the flux density, we chose a linear grid of ( 1000 × 1000 ) bins and normalized the color map with a power law of index p c m a p = 0.3 . Our results are summarized below:
Figure A5. L e , p time variation SED density plot for the complete numerical run of Section 3.1.
Figure A5. L e , p time variation SED density plot for the complete numerical run of Section 3.1.
Universe 10 00392 g0a5
Figure A6. B time variation SED density plot for the complete numerical run of Section 3.1.
Figure A6. B time variation SED density plot for the complete numerical run of Section 3.1.
Universe 10 00392 g0a6
Figure A7. p e , p time variation SED density plot for the complete numerical run of Section 3.2.
Figure A7. p e , p time variation SED density plot for the complete numerical run of Section 3.2.
Universe 10 00392 g0a7

Appendix F

Figure A8. Parameter histograms of the variations in the particle energy injection rate (left, center left), the magnetic field strength (center right), and the modified power-law index variations (right). The variations in the power-law index are driven by changes in the magnetic field strength, ensuring the histogram accurately reflects both scenarios.
Figure A8. Parameter histograms of the variations in the particle energy injection rate (left, center left), the magnetic field strength (center right), and the modified power-law index variations (right). The variations in the power-law index are driven by changes in the magnetic field strength, ensuring the histogram accurately reflects both scenarios.
Universe 10 00392 g0a8

Notes

1
https://github.com/mariapetro/LeHaMoC/ (accessed on 15 January 2024).
2
https://tools.ssdc.asi.it/SED/ (accessed on 15 January 2024).
3
4
https://www.swift.ac.uk/user_objects/ (accessed on 15 May 2024).
5
See note 4 above.
6
Attempts fitting data from 10 GHz to VHE γ -rays with a single-zone model do not yield meaningful solutions.
7
Both the mean cadence and time step are in the observer’s frame, as they are based on direct observational results. LeHaMoC performs calculations in the comoving frame of the blob in time-steps of 1 intrinsic light-crossing time. However, transformations between the comoving and the observer’s frame, do not change the relative values of our generated TS. Therefore, we can use them despite the different reference frames.
8
The optical depth τ γ γ transitions from ≈10 at the initial average state to τ γ γ 10 2 at the peak of the TS.
9
https://docs.gammapy.org/1.2/ (accessed on 1 June 2024).
10
We assume the source always to be centered at the field of view of the array.
11
Other empirical models such as a power law or an exponential cut-off power law could not adequately describe the data (resulting in larger residuals at the lower and higher energy cutoffs).
12
13
Links to videos showing the temporal evolution of the decomposed SED: video1, video2, video3, video4 (accessed on 10 August 2024).

References

  1. Urry, C.M.; Padovani, P. Unified Schemes for Radio-Loud Active Galactic Nuclei. Publ. Astron. Soc. Pac. 1995, 107, 803. [Google Scholar] [CrossRef]
  2. Padovani, P.; Alexander, D.M.; Assef, R.J.; De Marco, B.; Giommi, P.; Hickox, R.C.; Richards, G.T.; Smolčić, V.; Hatziminaoglou, E.; Mainieri, V.; et al. Active galactic nuclei: What’s in a name? Astron. Astrophys. Rev. 2017, 25, 2. [Google Scholar] [CrossRef]
  3. Aharonian, F.; Akhperjanian, A.G.; Bazer-Bachi, A.R.; Behera, B.; Beilicke, M.; Benbow, W.; Berge, D.; Bernlöhr, K.; Boisson, C.; Bolz, O.; et al. An Exceptional Very High Energy Gamma-ray Flare of PKS 2155-304. Astrophys. J. 2007, 664, L71–L74. [Google Scholar] [CrossRef]
  4. Abdo, A.A.; Ackermann, M.; Ajello, M.; Antolini, E.; Baldini, L.; Ballet, J.; Barbiellini, G.; Bastieri, D.; Bechtol, K.; Bellazzini, R.; et al. Gamma-ray Light Curves and Variability of Bright Fermi-detected Blazars. Astrophys. J. 2010, 722, 520–542. [Google Scholar] [CrossRef]
  5. Bloom, S.D.; Marscher, A.P. An Analysis of the Synchrotron Self-Compton Model for the Multi–Wave Band Spectra of Blazars. Astrophys. J. 1996, 461, 657. [Google Scholar] [CrossRef]
  6. Mastichiadis, A.; Kirk, J.G. Variability in the synchrotron self-Compton model of blazar emission. Astron. Astrophys. 1997, 320, 19–25. [Google Scholar] [CrossRef]
  7. Dermer, C.D.; Schlickeiser, R.; Mastichiadis, A. High-energy gamma radiation from extragalactic radio sources. Astron. Astrophys. 1992, 256, L27–L30. [Google Scholar]
  8. Ghisellini, G.; Madau, P. On the origin of the gamma-ray emission in blazars. Mon. Not. R. Astron. Soc. 1996, 280, 67–76. [Google Scholar] [CrossRef]
  9. Mastichiadis, A.; Moraitis, K. On the rapid TeV flaring activity of Markarian 501. Astron. Astrophys. 2008, 491, L37–L40. [Google Scholar] [CrossRef]
  10. Aharonian, F.A.; Barkov, M.V.; Khangulyan, D. Scenarios for Ultrafast Gamma-ray Variability in AGN. Astrophys. J. 2017, 841, 61. [Google Scholar] [CrossRef]
  11. Krawczynski, H.; Coppi, P.S.; Aharonian, F. Time-dependent modelling of the Markarian 501 X-ray and TeV gamma-ray data taken during 1997 March and April. Mon. Not. R. Astron. Soc. 2002, 336, 721–735. [Google Scholar] [CrossRef]
  12. Mastichiadis, A.; Petropoulou, M.; Dimitrakoudis, S. Mrk 421 as a case study for TeV and X-ray variability in leptohadronic models. Mon. Not. R. Astron. Soc. 2013, 434, 2684–2695. [Google Scholar] [CrossRef]
  13. Djannati-Atai, A.; Piron, F.; Barrau, A.; Iacoucci, L.; Punch, M.; Tavernet, J.P.; Bazer-Bachi, R.; Cabot, H.; Chounet, L.M.; Debiais, G.; et al. Very High Energy Gamma-ray spectral properties of MKN 501 from CAT Čerenkov telescope observations in 1997. Astron. Astrophys. 1999, 350, 17–24. [Google Scholar] [CrossRef]
  14. Petropoulou, M.; Mastichiadis, A.; Vasilopoulos, G.; Paneque, D.; Becerra González, J.; Zanias, F. TeV pion bumps in the gamma-ray spectra of flaring blazars. Astron. Astrophys. 2024, 685, A110. [Google Scholar] [CrossRef]
  15. Shukla, A.; Chitnis, V.R.; Singh, B.B.; Acharya, B.S.; Anupama, G.C.; Bhattacharjee, P.; Britto, R.J.; Mannheim, K.; Prabhu, T.P.; Saha, L.; et al. Multi-frequency, Multi-epoch Study of Mrk 501: Hints for a Two-component Nature of the Emission. Astrophys. J. 2015, 798, 2. [Google Scholar] [CrossRef]
  16. Stathopoulos, S.I.; Petropoulou, M.; Vasilopoulos, G.; Mastichiadis, A. LeHaMoC: A versatile time-dependent lepto-hadronic modeling code for high-energy astrophysical sources. Astron. Astrophys. 2024, 683, A225. [Google Scholar] [CrossRef]
  17. Blandford, R.; Meier, D.; Readhead, A. Relativistic Jets from Active Galactic Nuclei. Annu. Rev. Astron. Astrophys. 2019, 57, 467–509. [Google Scholar] [CrossRef]
  18. Hovatta, T.; Valtaoja, E.; Tornikoski, M.; Lähteenmäki, A. Doppler factors, Lorentz factors and viewing angles for quasars, BL Lacertae objects and radio galaxies. Astron. Astrophys. 2009, 494, 527–537. [Google Scholar] [CrossRef]
  19. Acciari, V.A. et al. [MAGIC Collaboration] Study of the variable broadband emission of Markarian 501 during the most extreme Swift X-ray activity. Astron. Astrophys. 2020, 637, A86. [Google Scholar] [CrossRef]
  20. Albert, A.; Alfaro, R.; Alvarez, C.; Angeles Camacho, J.R.; Arteaga-Velázquez, J.C.; Arunbabu, K.P.; Avila Rojas, D.; Ayala Solares, H.A.; Baghmanyan, V.; Belmont-Moreno, E.; et al. Long-term Spectra of the Blazars Mrk 421 and Mrk 501 at TeV Energies Seen by HAWC. Astrophys. J. 2022, 929, 125. [Google Scholar] [CrossRef]
  21. Biteau, J.; Williams, D.A. The Extragalactic Background Light, the Hubble Constant, and Anomalies: Conclusions from 20 Years of TeV Gamma-ray Observations. Astrophys. J. 2015, 812, 60. [Google Scholar] [CrossRef]
  22. Aleksić, J.; Ansoldi, S.; Antonelli, L.A.; Antoranz, P.; Babic, A.; Bangale, P.; Barres de Almeida, U.; Barrio, J.A.; Becerra González, J.; Bednarek, W.; et al. Multiwavelength observations of Mrk 501 in 2008. Astron. Astrophys. 2015, 573, A50. [Google Scholar] [CrossRef]
  23. Ahnen, M.L.; Ansoldi, S.; Antonelli, L.A.; Antoranz, P.; Babic, A.; Banerjee, B.; Bangale, P.; Barres de Almeida, U.; Barrio, J.A.; Becerra González, J.; et al. Multiband variability studies and novel broadband SED modeling of Mrk 501 in 2009. Astron. Astrophys. 2017, 603, A31. [Google Scholar] [CrossRef]
  24. Arbet-Engels, A.; Baack, D.; Balbo, M.; Biland, A.; Bretz, T.; Buss, J.; Dorner, D.; Eisenberger, L.; Elsaesser, D.; Hildebrand, D.; et al. Long-term multi-band photometric monitoring of Mrk 501. Astron. Astrophys. 2021, 655, A93. [Google Scholar] [CrossRef]
  25. Liodakis, I.; Marscher, A.P.; Agudo, I.; Berdyugin, A.V.; Bernardos, M.I.; Bonnoli, G.; Borman, G.A.; Casadio, C.; Casanova, V.; Cavazzuti, E.; et al. Polarized blazar X-rays imply particle acceleration in shocks. Nature 2022, 611, 677–681. [Google Scholar] [CrossRef]
  26. Abe, S. et al. [MAGIC Collaboration] Insights into the broadband emission of the TeV blazar Mrk 501 during the first X-ray polarization measurements. Astron. Astrophys. 2024, 685, A117. [Google Scholar] [CrossRef]
  27. Pian, E.; Vacanti, G.; Tagliaferri, G.; Ghisellini, G.; Maraschi, L.; Treves, A.; Urry, C.M.; Fiore, F.; Giommi, P.; Palazzi, E.; et al. BeppoSAX Observations of Unprecedented Synchrotron Activity in the BL Lacertae Object Markarian 501. Astrophys. J. 1998, 492, L17–L20. [Google Scholar] [CrossRef]
  28. Tavecchio, F.; Maraschi, L.; Pian, E.; Chiappetti, L.; Celotti, A.; Fossati, G.; Ghisellini, G.; Palazzi, E.; Raiteri, C.M.; Sambruna, R.M.; et al. Theoretical Implications from the Spectral Evolution of Markarian 501 Observed with BeppoSAX. Astrophys. J. 2001, 554, 725–733. [Google Scholar] [CrossRef]
  29. Costamante, L.; Ghisellini, G.; Giommi, P.; Tagliaferri, G.; Celotti, A.; Chiaberge, M.; Fossati, G.; Maraschi, L.; Tavecchio, F.; Treves, A.; et al. Extreme synchrotron BL Lac objects. Stretching the blazar sequence. Astron. Astrophys. 2001, 371, 512–526. [Google Scholar] [CrossRef]
  30. Aharonian, F.A.; Akhperjanian, A.G.; Barrio, J.A.; Bernlöhr, K.; Bolz, O.; Börst, H.; Bojahr, H.; Contreras, J.L.; Cortina, J.; Denninghoff, S.; et al. Reanalysis of the high energy cutoff of the 1997 Mkn 501 TeV energy spectrum. Astron. Astrophys. 2001, 366, 62–67. [Google Scholar] [CrossRef]
  31. Ahnen, M.L.; Ansoldi, S.; Antonelli, L.A.; Arcaro, C.; Babić, A.; Banerjee, B.; Bangale, P.; Barres de Almeida, U.; Barrio, J.A.; Becerra González, J.; et al. Extreme HBL behavior of Markarian 501 during 2012. Astron. Astrophys. 2018, 620, A181. [Google Scholar] [CrossRef]
  32. Bartoli, B.; Bernardini, P.; Bi, X.J.; Bleve, C.; Bolognino, I.; Branchini, P.; Budano, A.; Calabrese Melcarne, A.K.; Camarri, P.; Cao, Z.; et al. Long-term Monitoring of Mrk 501 for its Very High Energy γ Emission and a Flare in 2011 October. Astrophys. J. 2012, 758, 2. [Google Scholar] [CrossRef]
  33. Anderhub, H.; Antonelli, L.A.; Antoranz, P.; Backes, M.; Baixeras, C.; Balestra, S.; Barrio, J.A.; Bastieri, D.; Becerra González, J.; Becker, J.K.; et al. Simultaneous Multiwavelength Observation of Mkn 501 in a Low State in 2006. Astrophys. J. 2009, 705, 1624–1631. [Google Scholar] [CrossRef]
  34. Godambe, S.V.; Rannot, R.C.; Chandra, P.; Yadav, K.K.; Tickoo, A.K.; Venugopal, K.; Bhatt, N.; Bhattacharyya, S.; Chanchalani, K.; Dhar, V.K.; et al. Very high energy γ-ray observations of Mrk 501 using the TACTIC imaging γ-ray telescope during 2005 06. J. Phys. G Nucl. Phys. 2008, 35, 065202. [Google Scholar] [CrossRef]
  35. Acharyya, A.; Adams, C.B.; Archer, A.; Bangale, P.; Bartkoske, J.T.; Batista, P.; Benbow, W.; Brill, A.; Brose, R.; Buckley, J.H.; et al. VTSCat: The VERITAS Catalog of Gamma-ray Observations. Res. Notes Am. Astron. Soc. 2023, 7, 6. [Google Scholar] [CrossRef]
  36. Ackermann, M.; Ajello, M.; Atwood, W.B.; Baldini, L.; Ballet, J.; Barbiellini, G.; Bastieri, D.; Becerra Gonzalez, J.; Bellazzini, R.; Bissaldi, E.; et al. 2FHL: The Second Catalog of Hard Fermi-LAT Sources. Astrophys. J. Suppl. Ser. 2016, 222, 5. [Google Scholar] [CrossRef]
  37. Abdo, A.A.; Ackermann, M.; Ajello, M.; Allafort, A.; Antolini, E.; Atwood, W.B.; Axelsson, M.; Baldini, L.; Ballet, J.; Barbiellini, G.; et al. Fermi Large Area Telescope First Source Catalog. Astrophys. J. Suppl. Ser. 2010, 188, 405–436. [Google Scholar] [CrossRef]
  38. Nolan, P.L.; Abdo, A.A.; Ackermann, M.; Ajello, M.; Allafort, A.; Antolini, E.; Atwood, W.B.; Axelsson, M.; Baldini, L.; Ballet, J.; et al. Fermi Large Area Telescope Second Source Catalog. Astrophys. J. Suppl. Ser. 2012, 199, 31. [Google Scholar] [CrossRef]
  39. Acero, F.; Ackermann, M.; Ajello, M.; Albert, A.; Atwood, W.B.; Axelsson, M.; Baldini, L.; Ballet, J.; Barbiellini, G.; Bastieri, D.; et al. Fermi Large Area Telescope Third Source Catalog. Astrophys. J. Suppl. Ser. 2015, 218, 23. [Google Scholar] [CrossRef]
  40. Abdollahi, S.; Acero, F.; Ackermann, M.; Ajello, M.; Atwood, W.B.; Axelsson, M.; Baldini, L.; Ballet, J.; Barbiellini, G.; Bastieri, D.; et al. Fermi Large Area Telescope Fourth Source Catalog. Astrophys. J. Suppl. Ser. 2020, 247, 33. [Google Scholar] [CrossRef]
  41. Giommi, P.; Capalbi, M.; Fiocchi, M.; Memola, E.; Perri, M.; Piranomonte, S.; Rebecchi, S.; Massaro, E. A Catalog of 157 X-ray Spectra and 84 Spectral Energy Distributions of Blazars Observed with BeppoSAX. In Proceedings of the Blazar Astrophysics with BeppoSAX and Other Observatories; Giommi, P., Massaro, E., Palumbo, G., Eds.; ESA-ESRIN: Frascati, Italy, 2002; p. 63. [Google Scholar] [CrossRef]
  42. Verrecchia, F.; Giommi, P.; Santolamazza, P.; in’t Zand, J.J.M.; Granata, S.; Schuurmans, J.J. The BeppoSAX WFC source catalogue. In Proceedings of the Multicolored Landscape of Compact Objects and Their Explosive Origins; di Salvo, T., Israel, G.L., Piersant, L., Burderi, L., Matt, G., Tornambe, A., Menna, M.T., Eds.; American Institute of Physics Conference Series; AIP: Long Island, NY, USA, 2007; Volume 924, pp. 923–927. [Google Scholar] [CrossRef]
  43. Thiersen, H.; Zacharias, M.; Böttcher, M. Characterising the Long-Term Variability of Blazars in Leptonic Models. Galaxies 2019, 7, 35. [Google Scholar] [CrossRef]
  44. Rodrigues, X.; Paliya, V.S.; Garrappa, S.; Omeliukh, A.; Franckowiak, A.; Winter, W. Leptohadronic multi-messenger modeling of 324 gamma-ray blazars. Astron. Astrophys. 2024, 681, A119. [Google Scholar] [CrossRef]
  45. Abe, H.; Abe, S.; Acciari, V.A.; Agudo, I.; Aniello, T.; Ansoldi, S.; Antonelli, L.A.; Arbet-Engels, A.; Arcaro, C.; Artero, M.; et al. Multimessenger Characterization of Markarian 501 during Historically Low X-ray and γ-ray Activity. Astrophys. J. Suppl. Ser. 2023, 266, 37. [Google Scholar] [CrossRef]
  46. Liu, L.; Jiang, Y.; Deng, J.; Chen, Z.; Ma, C. Unveiling the Emission and Variation Mechanism of Mrk 501: Using the Multi-Wavelength Data at Different Time Scale. Universe 2024, 10, 114. [Google Scholar] [CrossRef]
  47. Katarzyński, K.; Sol, H.; Kus, A. The multifrequency emission of Mrk 501. From radio to TeV gamma-rays. Astron. Astrophys. 2001, 367, 809–825. [Google Scholar] [CrossRef]
  48. Foreman-Mackey, D.; Hogg, D.W.; Lang, D.; Goodman, J. emcee: The MCMC Hammer. Publ. Astron. Soc. Pac. 2013, 125, 306. [Google Scholar] [CrossRef]
  49. Guilbert, P.W.; Fabian, A.C.; Rees, M.J. Spectral and variability constraints on compact sources. Mon. Not. R. Astron. Soc. 1983, 205, 593–603. [Google Scholar] [CrossRef]
  50. Hillas, A.M. The Origin of Ultra-High-Energy Cosmic Rays. Annu. Rev. Astron. Astrophys. 1984, 22, 425–444. [Google Scholar] [CrossRef]
  51. Foreman-Mackey, D.; Agol, E.; Ambikasaran, S.; Angus, R. Fast and Scalable Gaussian Process Modeling with Applications to Astronomical Time Series. Astron. J. 2017, 154, 220. [Google Scholar] [CrossRef]
  52. Foreman-Mackey, D. Scalable Backpropagation for Gaussian Processes using Celerite. Res. Notes Am. Astron. Soc. 2018, 2, 31. [Google Scholar] [CrossRef]
  53. Abdollahi, S.; Ajello, M.; Baldini, L.; Ballet, J.; Bastieri, D.; Becerra Gonzalez, J.; Bellazzini, R.; Berretta, A.; Bissaldi, E.; Bonino, R.; et al. The Fermi-LAT Lightcurve Repository. Astrophys. J. Suppl. Ser. 2023, 265, 31. [Google Scholar] [CrossRef]
  54. Yang, S.; Yan, D.; Zhang, P.; Dai, B.; Zhang, L. Gaussian Process Modeling Fermi-LAT γ-ray Blazar Variability: A Sample of Blazars with γ-ray Quasi-periodicities. Astrophys. J. 2021, 907, 105. [Google Scholar] [CrossRef]
  55. Covino, S.; Landoni, M.; Sandrinelli, A.; Treves, A. Looking at Blazar Light-curve Periodicities with Gaussian Processes. Astrophys. J. 2020, 895, 122. [Google Scholar] [CrossRef]
  56. Burke, C.J.; Shen, Y.; Chen, Y.C.; Scaringi, S.; Faucher-Giguere, C.A.; Liu, X.; Yang, Q. Optical Variability of the Dwarf AGN NGC 4395 from the Transiting Exoplanet Survey Satellite. Astrophys. J. 2020, 899, 136. [Google Scholar] [CrossRef]
  57. Zhang, H.; Yan, D.; Zhang, L. Characterizing the γ-ray Variability of Active Galactic Nuclei with the Stochastic Process Method. Astrophys. J. 2022, 930, 157. [Google Scholar] [CrossRef]
  58. Zhang, H.; Yan, D.; Zhang, L. Gaussian Process Modeling Blazar Multiwavelength Variability: Indirectly Resolving Jet Structure. Astrophys. J. 2023, 944, 103. [Google Scholar] [CrossRef]
  59. Finke, J.D.; Dermer, C.D.; Böttcher, M. Synchrotron Self-Compton Analysis of TeV X-ray-Selected BL Lacertae Objects. Astrophys. J. 2008, 686, 181–194. [Google Scholar] [CrossRef]
  60. Abdalla, H.; Abe, H.; Acero, F.; Acharyya, A.; Adam, R.; Agudo, I.; Aguirre-Santaella, A.; Alfaro, R.; Alfaro, J.; Alispach, C.; et al. Sensitivity of the Cherenkov Telescope Array for probing cosmology and fundamental physics with gamma-ray propagation. J. Cosmology Astropart. Phys. 2021, 2021, 48. [Google Scholar] [CrossRef]
  61. Cherenkov Telescope Array Consortium; Acharya, B.S.; Agudo, I.; Al Samarai, I.; Alfaro, R.; Alfaro, J.; Alispach, C.; Alves Batista, R.; Amans, J.P.; Amato, E.; et al. Science with the Cherenkov Telescope Array; World Scientific: Singapore, 2019. [Google Scholar] [CrossRef]
  62. Kirk, J.G.; Rieger, F.M.; Mastichiadis, A. Particle acceleration and synchrotron emission in blazar jets. Astron. Astrophys. 1998, 333, 452–458. [Google Scholar] [CrossRef]
  63. Tammi, J.; Duffy, P. Particle-acceleration time-scales in TeV blazar flares. Mon. Not. R. Astron. Soc. 2009, 393, 1063–1069. [Google Scholar] [CrossRef]
  64. Donath, A.; Terrier, R.; Remy, Q.; Sinha, A.; Nigro, C.; Pintore, F.; Khélifi, B.; Olivera-Nieto, L.; Ruiz, J.E.; Brügge, K.; et al. Gammapy: A Python package for gamma-ray astronomy. Astron. Astrophys. 2023, 678, A157. [Google Scholar] [CrossRef]
  65. Finke, J.D.; Razzaque, S.; Dermer, C.D. Modeling the Extragalactic Background Light from Stars and Dust. Astrophys. J. 2010, 712, 238–249. [Google Scholar] [CrossRef]
  66. Cash, W. Parameter estimation in astronomy through application of the likelihood ratio. Astrophys. J. 1979, 228, 939–947. [Google Scholar] [CrossRef]
  67. Aartsen, M.G. et al. [IceCube Collaboration] Evidence for High-Energy Extraterrestrial Neutrinos at the IceCube Detector. Science 2013, 342, 1242856. [Google Scholar] [CrossRef]
  68. Aartsen, M.G. et al. [IceCube Collaboration] Neutrino emission from the direction of the blazar TXS 0506+056 prior to the IceCube-170922A alert. Science 2018, 361, 147–151. [Google Scholar] [CrossRef]
  69. Aartsen, M.G.; Ackermann, M.; Adams, J.; Aguilar, J.A.; Ahlers, M.; Ahrens, M.; Alispach, C.; Andeen, K.; Anderson, T.; Ansseau, I.; et al. Time-Integrated Neutrino Source Searches with 10 Years of IceCube Data. Phys. Rev. Lett. 2020, 124, 051103. [Google Scholar] [CrossRef]
  70. Abbasi, R. et al. [IceCube Collaboration] IceCube Data for Neutrino Point-Source Searches Years 2008–2018. arXiv 2021, arXiv:2101.0983633. [Google Scholar]
  71. Aartsen, M.G.; Abraham, K.; Ackermann, M.; Adams, J.; Aguilar, J.A.; Ahlers, M.; Ahrens, M.; Altmann, D.; Andeen, K.; Anderson, T.; et al. All-sky Search for Time-integrated Neutrino Emission from Astrophysical Sources with 7 yr of IceCube Data. Astrophys. J. 2017, 835, 151. [Google Scholar] [CrossRef]
  72. Celotti, A.; Ghisellini, G. The power of blazar jets. Mon. Not. R. Astron. Soc. 2008, 385, 283–300. [Google Scholar] [CrossRef]
  73. Barth, A.J.; Ho, L.C.; Sargent, W.L.W. Stellar Velocity Dispersion and Black Hole Mass in the Blazar Markarian 501. Astrophys. J. 2002, 566, L13–L16. [Google Scholar] [CrossRef]
  74. Petropoulou, M.; Dimitrakoudis, S.; Padovani, P.; Mastichiadis, A.; Resconi, E. Photohadronic origin of γ-ray BL Lac emission: Implications for IceCube neutrinos. Mon. Not. R. Astron. Soc. 2015, 448, 2412–2429. [Google Scholar] [CrossRef]
  75. Zdziarski, A.A.; Bottcher, M. Hadronic models of blazars require a change of the accretion paradigm. Mon. Not. R. Astron. Soc. 2015, 450, L21–L25. [Google Scholar] [CrossRef]
  76. Dimitrakoudis, S.; Mastichiadis, A.; Protheroe, R.J.; Reimer, A. The time-dependent one-zone hadronic model. First principles. Astron. Astrophys. 2012, 546, A120. [Google Scholar] [CrossRef]
  77. Petropoulou, M.; Psarras, F.; Giannios, D. Hadronic signatures from magnetically dominated baryon-loaded AGN jets. Mon. Not. R. Astron. Soc. 2023, 518, 2719–2734. [Google Scholar] [CrossRef]
  78. Polkas, M.; Petropoulou, M.; Vasilopoulos, G.; Mastichiadis, A.; Urry, C.M.; Coppi, P.; Bailyn, C. A numerical study of long-term multiwavelength blazar variability. Mon. Not. R. Astron. Soc. 2021, 505, 6103–6120. [Google Scholar] [CrossRef]
  79. Emmanoulopoulos, D.; McHardy, I.M.; Papadakis, I.E. Generating artificial light curves: Revisited and updated. Mon. Not. R. Astron. Soc. 2013, 433, 907–927. [Google Scholar] [CrossRef]
  80. Timmer, J.; König, M. On generating power law noise. Astron. Astrophys. 1995, 300, 707. [Google Scholar]
  81. Schreiber, T.; Schmitz, A. Improved Surrogate Data for Nonlinearity Tests. Phys. Rev. Lett. 1996, 77, 635–638. [Google Scholar] [CrossRef]
  82. Karaferias, A.S.; Vasilopoulos, G.; Petropoulou, M.; Jenke, P.A.; Wilson-Hodge, C.A.; Malacaria, C. A Bayesian approach for torque modelling of BeXRB pulsars with application to super-Eddington accretors. Mon. Not. R. Astron. Soc. 2023, 520, 281–299. [Google Scholar] [CrossRef]
  83. Wise, D.K.; Bristow-Johnson, R. Performance of low-order polynomial interpolators in the presence of oversampled input. In Proceedings of the Audio Engineering Society Convention 107, New York, NY, USA, 24–27 September 1999; Audio Engineering Society: New York, NY, USA, 1999. [Google Scholar]
  84. Giommi, P.; Perri, M.; Capalbi, M.; D’Elia, V.; Barres de Almeida, U.; Brandt, C.H.; Pollock, A.M.T.; Arneodo, F.; Di Giovanni, A.; Chang, Y.L.; et al. X-ray spectra, light curves and SEDs of blazars frequently observed by Swift. Mon. Not. R. Astron. Soc. 2021, 507, 5690–5702. [Google Scholar] [CrossRef]
  85. Stathopoulos, S.I.; Petropoulou, M.; Giommi, P.; Vasilopoulos, G.; Padovani, P.; Mastichiadis, A. High-energy neutrinos from X-rays flares of blazars frequently observed by the Swift X-ray Telescope. Mon. Not. R. Astron. Soc. 2022, 510, 4063–4079. [Google Scholar] [CrossRef]
  86. Buchner, J. UltraNest—A robust, general purpose Bayesian inference engine. J. Open Source Softw. 2021, 6, 3001. [Google Scholar] [CrossRef]
  87. Karavola, D.; Petropoulou, M. A closer look at the electromagnetic signatures of Bethe-Heitler pair production process in blazars. J. Cosmology Astropart. Phys. 2024, 2024, 006. [Google Scholar] [CrossRef]
  88. Mastichiadis, A.; Kirk, J.G. Self-consistent particle acceleration in active galactic nuclei. Astron. Astrophys. 1995, 295, 613. [Google Scholar]
  89. Mastichiadis, A.; Petropoulou, M. Hadronic X-ray Flares from Blazars. Astrophys. J. 2021, 906, 131. [Google Scholar] [CrossRef]
  90. Chatzis, M.; Petropoulou, M.; Vasilopoulos, G. Radio emission from colliding outflows in high-mass X-ray binaries with strongly magnetized neutron stars. Mon. Not. R. Astron. Soc. 2022, 509, 2532–2550. [Google Scholar] [CrossRef]
  91. Goyal, A.; Stawarz, Ł.; Zola, S.; Marchenko, V.; Soida, M.; Nilsson, K.; Ciprini, S.; Baran, A.; Ostrowski, M.; Wiita, P.J.; et al. Stochastic Modeling of Multiwavelength Variability of the Classical BL Lac Object OJ 287 on Timescales Ranging from Decades to Hours. Astrophys. J. 2018, 863, 175. [Google Scholar] [CrossRef]
  92. Sobolewska, M.A.; Siemiginowska, A.; Kelly, B.C.; Nalewajko, K. Stochastic Modeling of the Fermi/LAT γ-ray Blazar Variability. Astrophys. J. 2014, 786, 143. [Google Scholar] [CrossRef]
Figure 1. Spectral energy distribution of Mrk 501. Gray points indicate archival observations taken from SED builder. We highlight VHE flares and the major X-ray outburst of 1997 with a faded gray color; both were excluded from our average-state analysis. Black points represent the archival X-ray data and the optical bulge of the host galaxy that were also not used in our modeling process. Colored markers indicate the long-term average X-ray and HE γ -ray spectra included in the fitting process (see inset legend).
Figure 1. Spectral energy distribution of Mrk 501. Gray points indicate archival observations taken from SED builder. We highlight VHE flares and the major X-ray outburst of 1997 with a faded gray color; both were excluded from our average-state analysis. Black points represent the archival X-ray data and the optical bulge of the host galaxy that were also not used in our modeling process. Colored markers indicate the long-term average X-ray and HE γ -ray spectra included in the fitting process (see inset legend).
Universe 10 00392 g001
Figure 2. One hundred leptonic SEDs computed from the posterior parameter space shown in Figure A2. The shaded band has been excluded from the fit assuming emission originating from a more extended jet region. The green dashed line represents the 1% limit used to calculate the upper limit on L p . With (faded) gray points, the archival (flare) observations are shown. Black points represent the archival X-ray data not used in our modeling process. Colored markers indicate the long-term average X-ray and HE γ -ray spectra included in the fitting process (see inset legend).
Figure 2. One hundred leptonic SEDs computed from the posterior parameter space shown in Figure A2. The shaded band has been excluded from the fit assuming emission originating from a more extended jet region. The green dashed line represents the 1% limit used to calculate the upper limit on L p . With (faded) gray points, the archival (flare) observations are shown. Black points represent the archival X-ray data not used in our modeling process. Colored markers indicate the long-term average X-ray and HE γ -ray spectra included in the fitting process (see inset legend).
Universe 10 00392 g002
Figure 3. Spectral components of the long-term average SED of Mrk 501 according to a leptohadronic model for the parameters listed in Table 1 and Table 2. Highlighted are the leptonic components (synchrotron and Inverse Compton (IC) of primary electrons) and the hadronic components (synchrotron and IC of secondary electrons, proton synchrotron, and emission from π 0 -decay) alongside the synchrotron and IC emission created from γ γ -injected pairs.
Figure 3. Spectral components of the long-term average SED of Mrk 501 according to a leptohadronic model for the parameters listed in Table 1 and Table 2. Highlighted are the leptonic components (synchrotron and Inverse Compton (IC) of primary electrons) and the hadronic components (synchrotron and IC of secondary electrons, proton synchrotron, and emission from π 0 -decay) alongside the synchrotron and IC emission created from γ γ -injected pairs.
Universe 10 00392 g003
Figure 4. (Left): Celerite fit (blue line) of Mrk 501 overplotted on the logarithmic energy flux LC from the Fermi LCR (black markers). (Right): Posterior PSD of the celerite fit with emcee. Indicated is the break frequency that marks the transition from a white-noise to a red-noise-like PSD (see Appendix C for slope discussion).
Figure 4. (Left): Celerite fit (blue line) of Mrk 501 overplotted on the logarithmic energy flux LC from the Fermi LCR (black markers). (Right): Posterior PSD of the celerite fit with emcee. Indicated is the break frequency that marks the transition from a white-noise to a red-noise-like PSD (see Appendix C for slope discussion).
Universe 10 00392 g004
Figure 5. (Upper panel): Zero-mean synthetic LC for Mrk 501 (black stars) and the interpolated LC on time steps of 1 t c r = 0.55 days (blue circles). Gray bands designate analyzed flares. Lower panel: Zoom-in of the two gray bands of the upper panel, showing the extreme flare (left) and the typical flare (right).
Figure 5. (Upper panel): Zero-mean synthetic LC for Mrk 501 (black stars) and the interpolated LC on time steps of 1 t c r = 0.55 days (blue circles). Gray bands designate analyzed flares. Lower panel: Zoom-in of the two gray bands of the upper panel, showing the extreme flare (left) and the typical flare (right).
Universe 10 00392 g005
Figure 6. Components of the SED of Mrk 501 at the peak of the TS describing the extreme flare (see Figure 5). (Upper panel): Variations of the particle energy luminosity. (Lower panel): Variations of the magnetic field strength. Highlighted with a black dotted line is the average-state emission of Figure 3.
Figure 6. Components of the SED of Mrk 501 at the peak of the TS describing the extreme flare (see Figure 5). (Upper panel): Variations of the particle energy luminosity. (Lower panel): Variations of the magnetic field strength. Highlighted with a black dotted line is the average-state emission of Figure 3.
Universe 10 00392 g006
Figure 7. Flux histograms in different energy bands for the full (14.6 yr) numerical run with time variations in the particle injection luminosity L e , p . Upper panel: Energy ranges of 0.3–3 TeV and 0.1–100 GeV. Middle panel: Energy ranges of 14–195 keV and 0.2–10 keV. Lower panel: Energy range of the R-band and the histogram for the simulated LCs. The observational flux histograms are overplotted in green whenever available (for more details, see text).
Figure 7. Flux histograms in different energy bands for the full (14.6 yr) numerical run with time variations in the particle injection luminosity L e , p . Upper panel: Energy ranges of 0.3–3 TeV and 0.1–100 GeV. Middle panel: Energy ranges of 14–195 keV and 0.2–10 keV. Lower panel: Energy range of the R-band and the histogram for the simulated LCs. The observational flux histograms are overplotted in green whenever available (for more details, see text).
Universe 10 00392 g007
Figure 8. Same as Figure 7 for time variations of the magnetic field strength B.
Figure 8. Same as Figure 7 for time variations of the magnetic field strength B.
Universe 10 00392 g008
Figure 9. Components of the SED of Mrk 501 at the peak of the TS describing the extreme flare (see Figure 5). Contrasted are leptohadronic (left column) to leptonic models (right column) for mild (lower row) to extreme (upper row) power-law index variations, with m = 0.5 and m = 4 , respectively. Highlighted with a black dotted line is the average-state emission of Figure 3.
Figure 9. Components of the SED of Mrk 501 at the peak of the TS describing the extreme flare (see Figure 5). Contrasted are leptohadronic (left column) to leptonic models (right column) for mild (lower row) to extreme (upper row) power-law index variations, with m = 0.5 and m = 4 , respectively. Highlighted with a black dotted line is the average-state emission of Figure 3.
Universe 10 00392 g009
Figure 10. Components of the SED of Mrk 501 at the peak of the injected TS describing the extreme (left panel) and typical (right panel) flare during the modified power-law index variations (for details, see text). Highlighted with a black dotted line is the average-state emission of Figure 3.
Figure 10. Components of the SED of Mrk 501 at the peak of the injected TS describing the extreme (left panel) and typical (right panel) flare during the modified power-law index variations (for details, see text). Highlighted with a black dotted line is the average-state emission of Figure 3.
Universe 10 00392 g010
Figure 11. Same as Figure 7 for time variations of the power-law index p e , p .
Figure 11. Same as Figure 7 for time variations of the power-law index p e , p .
Universe 10 00392 g011
Figure 12. LCs in three energy ranges of the extreme (solid line) and typical (dashed line) flares produced in a leptohadronic model with the modified power-law index variations (for more details, see text). The TS of the flares are rescaled to fit the display range and are overplotted for comparison (black lines).
Figure 12. LCs in three energy ranges of the extreme (solid line) and typical (dashed line) flares produced in a leptohadronic model with the modified power-law index variations (for more details, see text). The TS of the flares are rescaled to fit the display range and are overplotted for comparison (black lines).
Universe 10 00392 g012
Figure 13. Flux–flux diagrams for the modified power-law index variations. We plotted the Fermi flux against (a) the VHE flux, (b) the harder X-ray flux probed by BAT, and (c) the softer X-ray flux probed by XRT. In panel (d), we plot the hard X-ray flux against the R-band flux. We present the variability for the complete 14.6-year-long leptohadronic run (blue circles) and highlight the results for the extreme (yellow stars) and typical (purple squares) flares.
Figure 13. Flux–flux diagrams for the modified power-law index variations. We plotted the Fermi flux against (a) the VHE flux, (b) the harder X-ray flux probed by BAT, and (c) the softer X-ray flux probed by XRT. In panel (d), we plot the hard X-ray flux against the R-band flux. We present the variability for the complete 14.6-year-long leptohadronic run (blue circles) and highlight the results for the extreme (yellow stars) and typical (purple squares) flares.
Universe 10 00392 g013
Figure 14. Simulated LCs of the extreme flare of Mrk 501 for CTAO North, using exposure times of 30 min, 1 h, and 5 h. Here, we used the EBL attenuated modified power-law index variation model. The 1 h and 5 h points were shifted by 0.25 and 0.5 days to the right, respectively, for clarity.
Figure 14. Simulated LCs of the extreme flare of Mrk 501 for CTAO North, using exposure times of 30 min, 1 h, and 5 h. Here, we used the EBL attenuated modified power-law index variation model. The 1 h and 5 h points were shifted by 0.25 and 0.5 days to the right, respectively, for clarity.
Universe 10 00392 g014
Figure 15. Comparison of CTAO simulated LCs based on the leptohadronic (red) and leptonic (black) models for the energy ranges of 1–10 TeV (left panel) and 0.05–1 TeV (right panel). The leptonic points are shifted by 0.25 days to the right, for clarity.
Figure 15. Comparison of CTAO simulated LCs based on the leptohadronic (red) and leptonic (black) models for the energy ranges of 1–10 TeV (left panel) and 0.05–1 TeV (right panel). The leptonic points are shifted by 0.25 days to the right, for clarity.
Universe 10 00392 g015
Figure 16. Log-parabola fits and 1 σ contours to the leptohadronic (red) and leptonic (black) simulated CTAO spectra, assuming 5 h exposure time, extracted at 11 days (left panel) and 43 days (right panel) from the simulated extreme flare. The light curves are displayed in Figure 15.
Figure 16. Log-parabola fits and 1 σ contours to the leptohadronic (red) and leptonic (black) simulated CTAO spectra, assuming 5 h exposure time, extracted at 11 days (left panel) and 43 days (right panel) from the simulated extreme flare. The light curves are displayed in Figure 15.
Universe 10 00392 g016
Figure 17. Contour plots of best-fit log-parabola parameters of 1000 leptonic and leptohadronic simulated CTA spectra, showing ϕ 0 (left panel), α (middle panel), and β (right panel) for t = 11 (left figure) and t = 43 (right figure) days.
Figure 17. Contour plots of best-fit log-parabola parameters of 1000 leptonic and leptohadronic simulated CTA spectra, showing ϕ 0 (left panel), α (middle panel), and β (right panel) for t = 11 (left figure) and t = 43 (right figure) days.
Universe 10 00392 g017
Figure 18. Upper panel: Luminosity ratio between the MeV L 100 MeV and 14–195 keV L 14 195 keV bands using the produced SEDs when varying the slope of the particles p e , p (see middle panel) according to Equation (7). Bottom panel: Expected muon neutrino events observed by IceCube using Equation (9). The solid blue line indicates the result obtained using different IceCube configurations over different operation seasons (Table 5), while the dashed magenta line shows the result when considering only the most recent IceCube configuration (IC86 II) for all times. Gray areas designate analyzed flares, as described in Section 3, while t 0 = 55,029.0 MJD (17 July 2009) is the assumed starting time of the simulated LC.
Figure 18. Upper panel: Luminosity ratio between the MeV L 100 MeV and 14–195 keV L 14 195 keV bands using the produced SEDs when varying the slope of the particles p e , p (see middle panel) according to Equation (7). Bottom panel: Expected muon neutrino events observed by IceCube using Equation (9). The solid blue line indicates the result obtained using different IceCube configurations over different operation seasons (Table 5), while the dashed magenta line shows the result when considering only the most recent IceCube configuration (IC86 II) for all times. Gray areas designate analyzed flares, as described in Section 3, while t 0 = 55,029.0 MJD (17 July 2009) is the assumed starting time of the simulated LC.
Universe 10 00392 g018
Table 1. MCMC fitting results for the leptonic model describing the average SED of Mrk 501. The errors indicate the 68% range of values from the posterior distributions. The priors are uniform distributions in log 10 space (except for the p e , where the distribution is in linear space).
Table 1. MCMC fitting results for the leptonic model describing the average SED of Mrk 501. The errors indicate the 68% range of values from the posterior distributions. The priors are uniform distributions in log 10 space (except for the p e , where the distribution is in linear space).
Parameter DescriptionParameter SymbolValueUniform Priors
Blob radius log ( R [ cm ] ) 16.46 ± 0.03 [ 13 , 16.5 ]
Magnetic field strength log ( B [ G ] ) 1.42 ± 0.05 [ 8 , 8 ]
Doppler factor log ( δ ) 1.31 ± 0.02 [ 0 , 5 ]
Min. electron Lorentz factor log ( γ e m i n ) 2.75 ± 0.05 [ 0 , 5 ]
Max. electron Lorentz factor log ( γ e m a x ) 7.4 ± 0.4 [ 6 , 9 ]
Electron compactness 1 log ( l e ) 4.75 ± 0.03 [ 6 , 1 ]
Electron power law index p e 2.50 ± 0.02 [ 1 , 5 ]
1 The injected isotropically distributed luminosity of the electron population can be quantified via the compactness [49], a measure comparing the energy contents of a region to its spatial dimensions, l e = σ T L e / ( 4 π R m e c 3 ) .
Table 2. Parameters describing the relativistic proton component in the leptohadronic model describing the average SED of Mrk 501.
Table 2. Parameters describing the relativistic proton component in the leptohadronic model describing the average SED of Mrk 501.
Parameter DescriptionParameter SymbolValue
Min. proton Lorentz factor log ( γ p m i n ) 2.75
Max. proton Lorentz factor 1 log ( γ p m a x ) 7.4
Proton power law index p p 2.50
Injected proton luminosity L p [erg/s] 10 5.61 L e
1 The Hillas criterion [50] imposes a geometrical limit on the maximum particle energy in a region of a given size R and magnetic field strength B. In particular, it requires the Larmor radius of the particle to be smaller than the confinement region, r L R . Thus, γ H = q B R / ( m c 2 ) gives the maximum allowed Lorentz factor. For Mrk 501, using the lower 68% ranges of the radius and magnetic field strength listed in Table 1 yields the result of l o g ( γ p H ) = 8.5 > l o g ( γ p m a x ) .
Table 3. celerite parameters for the Fermi-LAT LC of Mrk 501.
Table 3. celerite parameters for the Fermi-LAT LC of Mrk 501.
Source Name μ   1 ln ( S 0 ) ln ( Q ) ln ( ω 0 )   2
Mrk 501 4 . 12 0.03 + 0.03 0 . 6 0.4 + 0.5 3 . 2 1.1 + 0.8 1 . 4 0.8 + 1.1
1 Where μ is the mean of the log 10 (Flux) in the 0.1–100 GeV range in units of MeV cm−2 s−1. 2 Where ω 0 is in units of rad/days.
Table 4. Best-fit parameters of the log-parabolic model shown in Figure 16.
Table 4. Best-fit parameters of the log-parabolic model shown in Figure 16.
Log-Parabola ParameterLeptohadronic (43 Days)Leptonic (43 Days)Leptohadronic (11 Days)Leptonic (11 Days)
ϕ 0 1 5.2 ± 0.1 4.8 ± 0.1 39.9 ± 0.3 39.8 ± 0.3
α 2.19 ± 0.02 2.11 ± 0.03 2.353 ± 0.001 2.351 ± 0.01
β 0.10 ± 0.01 0.14 ± 0.02 0.121 ± 0.0001 0.123 ± 0.0001
1 in units of 10 12 cm−2 s−1 TeV−1.
Table 5. Point-source effective areas for different configurations used in our analysis. Data are taken from [69,70].
Table 5. Point-source effective areas for different configurations used in our analysis. Data are taken from [69,70].
IceCube ConfigurationSeason (MJD)
IC4054,562–54,971
IC5954,971–55,347
IC7955,347–55,694
IC86-I55,694–56,062
IC86-II>56,062
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Chatzis, M.; Stathopoulos, S.I.; Petropoulou, M.; Vasilopoulos, G. Searching for Hadronic Signatures in the Time Domain of Blazar Emission: The Case of Mrk 501. Universe 2024, 10, 392. https://doi.org/10.3390/universe10100392

AMA Style

Chatzis M, Stathopoulos SI, Petropoulou M, Vasilopoulos G. Searching for Hadronic Signatures in the Time Domain of Blazar Emission: The Case of Mrk 501. Universe. 2024; 10(10):392. https://doi.org/10.3390/universe10100392

Chicago/Turabian Style

Chatzis, Margaritis, Stamatios I. Stathopoulos, Maria Petropoulou, and Georgios Vasilopoulos. 2024. "Searching for Hadronic Signatures in the Time Domain of Blazar Emission: The Case of Mrk 501" Universe 10, no. 10: 392. https://doi.org/10.3390/universe10100392

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop