© The Authors 2024

1. Introduction

Blazars are among the brightest objects in the γ-ray sky and are known to emit radiation over a broad range of wavebands, from the radio up to the very-high-energy (VHE; > 0.1 TeV) regime. They are jetted active galactic nuclei (AGNs) with a small angle between our line of sight and the jet axis.

Even though blazars have been observed for decades across the entire electromagnetic spectrum, their main acceleration and emission mechanisms are still debated. Polarization measurements in the different wavebands could help us distinguish between different emission mechanisms (Zhang & Böttcher 2013) or allow properties such as the underlying acceleration mechanisms or the magnetic field configurations to be probed (Tavecchio et al. 2018; Tavecchio 2021). Until recently, polarization measurements of blazars had only been possible in the radio to optical bands. However, these bands are known for originating from multiple regions or regions in the jet that are more extended compared to the high-energy emission and can be contaminated by the host galaxy (see, e.g., Acciari et al. 2020). In December 2021, the Imaging X-ray Polarimetry Explorer (IXPE; Weisskopf 2022) was launched with the goal of measuring linear polarization in the 2–8 keV band for various sources.

For the first time, X-ray polarization of the archetypal blazar Markarian 501 (Mrk 501; z = 0.034; Ulrich et al. 1975) was detected (Liodakis et al. 2022). Mrk 501 is one of our closest and brightest blazars and therefore a prime object to explore using new observational techniques. As is usual for blazars, its spectral energy distribution (SED) has two peaks. The low-energy peak can be attributed to synchrotron radiation of relativistic electrons inside the magnetized plasma of the jet. It is used to classify blazars, and Mrk 501 is usually attributed to the class of high-synchrotron-peaked (HSP) blazars with a synchrotron peak frequency ν_s >  10¹⁵ Hz (Abdo et al. 2010). However, it has also shown extreme high-synchrotron-peaked (EHSP) behavior, with ν_s ≥ 2.4 × 10¹⁷ Hz (∼1 keV; Costamante et al. 2001; Abdo et al. 2010), during extended periods of time that encompassed both flaring and non-flaring activity (Ahnen et al. 2018; Furniss et al. 2015). The origin of the high-energy peak is less clear. It could either originate from electron inverse-Compton scattering off low-energy photons (leptonic scenarios; see, e.g., Maraschi et al. 1992; Ghisellini & Maraschi 1996; Tavecchio et al. 1998) or be induced by relativistic protons inside the jet that spark different cascades or hadronic synchrotron radiation (hadronic scenarios; see, e.g., Mannheim 1993; Aharonian 2000; Mücke & Protheroe 2001).

For Mrk 501, the IXPE energy range covers either the falling part of the synchrotron bump, during its usual HSP states, or the peak, during EHSP behavior. IXPE thus probes the emission from the most energetic electrons freshly accelerated in the jet.

The first X-ray polarization measurements of Mrk 501 were conducted during two pointings, carried out March 8-10 and March 27-29, 2022 (MJD 59646 to 59648 and MJD 59665 to 59667). In what follows, these two periods are referred to as IXPE-1 and IXPE-2. They reveal a polarization degree in the 2–8 keV band of ∼10%, which is a factor of 2 higher than in the optical regime (Liodakis et al. 2022). Additionally, the polarization angle was aligned with the radio jet. All these properties point to shocks being the prime acceleration mechanism; they also suggest an energy-stratified jet where the most energetic particles emit from a region that is magnetically more ordered and closer to the acceleration site (Liodakis et al. 2022; Angelakis et al. 2016). These energetic particles subsequently cool and stream away from the shock to populate broader regions. An additional IXPE pointing was conducted July 9-12, 2022 (MJD 59769 to 59772), revealing a drop in polarization degree to ∼7%, but with a stable polarization angle compared to the March measurement (Lisalda, in prep.). This observing epoch is dubbed IXPE-3 in the rest of this paper.

In this work, we present the multiwavelength (MWL) picture around these three X-ray polarization measurements, including, for the first time, data up to the VHE regime taken by the Florian Goebel Major Atmospheric Gamma Imaging Cherenkov (MAGIC) telescopes. To characterize the high-energy emission in detail, the campaign involved instruments such as the Neil Gehrels Swift Observatory (Swift) and the Nuclear Spectroscopic Telescope Array (NuSTAR) in the X-rays and the Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope (Fermi) in the high-energy γ-rays. Both the IXPE pointings and the whole time epoch from March 1 to July 19, 2022 (MJD 59639 to MJD 59779) were accompanied by extensive MWL coverage from the radio to the VHE band, which is summarized in the following sections.

This paper is structured as follows: in Sect. 2 we describe all the instruments involved together with their dedicated data analysis methods. Section 3 provides a detailed characterization of the MWL emission during the three IXPE observations, including the different flux states, spectral behaviors, and polarization. The MWL characteristics of the full campaign are then summarized in Sect. 4. Section 5 presents a theoretical modeling of the three IXPE time epochs using a two-zone leptonic model. Finally, Sect. 6 summarizes and discusses the MWL results together with the obtained theoretical models.

2. Observations and analysis

A multitude of observations from VHE to radio were performed to investigate the full MWL behavior of Mrk 501 around the first IXPE observations. These were conducted in the framework of the extensive multi-instrument campaigns organized since 2008 for Mrk 501 (Aleksić et al. 2015) and are stretching from March 1 to July 19, 2022 (MJD 59639 to MJD 59779) with the instruments and analyses summarized in this section.

2.1. MAGIC

MAGIC is a stereoscopic system of two imaging air Cherenkov telescopes located at the Roque de los Muchachos Observatory, on the Canary island of La Palma at 2243 m above sea level. MAGIC covers an energy range between tens of GeVs to tens of TeVs and has a sensitivity above 100 GeV (300 GeV) of about 2% (about 1%) of the Crab Nebula flux at low zenith angles (< 30°) after 25 h of observations (see Fig. 19 of Aleksić et al. 2016), which makes it well suited to perform blazar observations in the VHE range.

Covering the time period around the first three IXPE pointings, MAGIC observed Mrk 501 for around 38 h from March 1 to July 19, 2022 (MJD 59639 to MJD 59779), which results in 26.5 h after data selection cuts based on atmospheric transmission. The data are analyzed using the MAGIC Analysis and Reconstruction Software (MARS) package (Zanin et al. 2013; Aleksić et al. 2016).

The VHE flux light curve (LC) was computed for the full energy range above 0.2 TeV, as well as for a low-energy range of 0.2–1 TeV and a high-energy range above 1 TeV as shown in the top panel of Fig. 1 using nightly bins. For the bins showing a significance of less than 2σ, upper limits are computed according to the Rolke method (Rolke et al. 2005) with 95% confidence. The spectral parameters around the three IXPE pointings were determined using a forward folding method assuming a power-law model, which provides a good description of the data (see Table 1). For the IXPE-1 and IXPE-2 spectra, the Tikhonov unfolding method was applied to compute the spectral data points (Albert et al. 2007), while for the IXPE-3 spectrum again the forward folding method was used due to limited statistics. We checked the preference of a log parabola over a simple power-law model for each spectrum using a likelihood ratio test, but no significant (> 3σ) preference is found. Extragalactic background light absorption was corrected for using the Domínguez et al. (2011) model.

Fig. 1. MWL LC for Mrk 501 between March 1, 2022, and July 19, 2022 (MJD 59639 to MJD 59779). The gray areas mark the IXPE observations from March 8-10, 2022 (MJD 59646 to MJD 59648), March 27-29, 2022 (MJD 59665 to MJD 59667), and July 9-12, 2022 (MJD 59769 to MJD 59772). Top to bottom: (i) MAGIC fluxes in daily bins, with the arrows displaying the upper limits (95% confidence level) for the non-significant bins (< 2 σ); (ii) Fermi-LAT fluxes in 7-day bins; (iii) X-ray fluxes in daily bins, including Swift-XRT, NuSTAR, and IXPE; (iv) hardness ratio between the high- and low-energy fluxes of Swift-XRT; (v) Swift-UVOT; (vi) optical R-band data from RoboPol, Calar Alto, T090, MAPCAT, T150, and AZT-8 in daily bins; (vii) radio data including OVRO, Metsähovi, IRAM, and SMA; and polarization degree (viii) and polarization angle (ix) observations in the X-rays from IXPE, the optical R band from NOT, RoboPol, Calar Alto, MAPCAT, T150, and AZT-8; and the radio band from IRAM and SMA. The IXPE and AZT-8 data are taken from Liodakis et al. (2022) and Lisalda (in prep.). In the polarization angle panel, we plot the average direction angle determined at 43 GHz by Weaver et al. (2022) with a horizontal gray band. See Sect. 4 for further details.

2.2. Fermi-LAT

The LAT is a pair conversion instrument on board Fermi, which operates in the high-energy regime and is sensitive to an energy range from 20 MeV to beyond 300 GeV. It covers the entire sky every 3 h (Atwood et al. 2009; Ackermann et al. 2012).

To analyze the LAT data, we used the unbinned-likelihood tools from the FERMITOOLS software package¹ (v2.0.8). We employed P8R3_SOURCE_V3_v1 as the instrument response function and gll_iem_v07 and iso_P8R3_SOURCE_V3_v1 as the diffuse background² and used the event class of 128, considering events that interacted in the front and back sections of the tracker (i.e., evtype = 3). The radius of the region of interest (ROI) was set to 10° and the maximum zenith to 100° because of the chosen energy range. We chose a higher minimum energy than is conventionally chosen (0.3 GeV instead of 0.1 GeV) to ensure a better angular resolution (2° for the 68% containment for photons above 0.3 GeV compared to 5° for photons above 0.1 GeV), which leads to an improvement in the signal-to-noise ratio (S/N) since the analysis is less affected by diffuse backgrounds or potential (transient) unaccounted for neighboring sources. This is possible owing to the hard-spectrum of Mrk 501 because the reduction of detected photons due to the higher threshold is relatively small. For the maximum energy, we chose a value of 500 GeV that allows an overlap with the MAGIC energy range when reconstructing spectra.

To build a model, we used the fourth Fermi-LAT source catalog (4FGL-DL3; Abdollahi et al. 2022) using all sources within the ROI plus an annulus of 5°. The obtained model was then fit to the dataset in the time range March 1 to July 19, 2022 (MJD 59639 to MJD 59779). Subsequently, very weak components with a test statistic (Mattox et al. 1996) below 3 were removed from the model as well as those with an expected source count below 1. We then divided the dataset into 20 bins of 7-day durations and redid the fit in each bin, varying the normalization of the diffuse backgrounds and that of sources within < 3° and a detection with a test statistic > 10. Furthermore, we allowed the spectral parameters of Mrk 501 to vary, assuming a power-law shape. To evaluate the impact of very variable sources in our ROI, we conducted the same analysis freeing the normalization and spectral indices of all very-variable sources (variability index > 100). Both analyses agree within the statistical uncertainties. In addition, spectral analyses are performed for 2-week intervals centered around each of the three IXPE observations. The same ROI model and approach as before was applied. For each spectrum, we computed the likelihood ratio between a power-law model and log parabolic power-law model. No significant (> 3σ) preference is found for the log parabolic fit, and therefore a power law was chosen. The number of spectral bins for the SEDs is chosen according to the time intervals and flux levels to assure at least four SED points for each of the three states described in Sect. 3.3. The obtained parameters are summarized in Table 1.

2.3. NuSTAR

This work includes four multi-hour exposures from NuSTAR (Harrison et al. 2013). The NuSTAR instrument consists of two co-aligned X-ray telescopes focusing on two independent focal plane modules, FPMA and FPMB, and provides unprecedented sensitivity in the 3–79 keV band. In this work, we analyze all NuSTAR observations from 2022, which all took place on March 9, 22, 24, and 27, 2022 (MJD 59647, MJD 59660, MJD 59662, MJD 59665; observation ID 60701032002, 60502009002, 60502009004, and 60702062004, respectively).

The raw data were processed using the NuSTAR Data Analysis Software (NuSTARDAS) package v.2.1.1 and CALDB version 20220912. The events were screened in the nupipeline process with the flags tentacle=yes and saamode=optimized to remove potential background increase caused by the South Atlantic Anomaly passages. The source counts are obtained from a circular region centered around Mrk 501 with a radius of ≈140″. The background events are extracted from a nearby source-free circular region that has the same radius. The spectra were then grouped with the grppha task to obtain at least 20 counts in each energy bin.

For all exposures, the source spectra dominate over the background up to roughly ≈30 keV. Hence, we decided to quote fluxes only up to 30 keV, and in two separate energy bands: 3–7 keV and 7–30 keV. The best-fit spectral parameters were obtained in the full NuSTAR bandpass, 3–79 keV, and averaged over the respective observations. We fitted the spectra using XSPEC (Arnaud 1996) and assuming a log-parabolic function with a normalization energy fixed to 1 keV. A log-parabola model is significantly preferred over a simple power law. Here, and for the rest of the X-ray analysis performed in this work, a photoelectric absorption component was added to the model assuming an equivalent hydrogen column density fixed to N_H = 1.7  ×  10²⁰ cm⁻² (HI4PI Collaboration 2016). The fluxes and spectral parameters were computed by simultaneously fitting FPMA and FPMB. The cross-calibration factor between the two focal module planes for all bins is below 5%, thus well within the expected systematics (Madsen et al. 2015).

For all NuSTAR pointings, the variability is small and the fluxes vary at most by ≈10% within each observation. Thus, the fluxes and spectral parameters were averaged over the entire exposure time from each observation.

2.4. IXPE

The IXPE telescope (Weisskopf 2022) was launched in December 2021 and measures polarization in the 2–8 keV regime. In this work, we used the first three pointing conducted on Mrk 501. The first two took place from March 8-10, 2022 (MJD 59646 to MJD 59648) and from March 27-29, 2022 (MJD 59665 to MJD 59667) and the results are published in Liodakis et al. (2022), from which we extracted the flux and polarization values from Tables 1 and 2. A third pointing took place from July 9-12, 2022 (MJD 59769 to MJD 59772), with the results published in Lisalda (in prep.). We extracted the polarization values from Table 2 of this paper and integrated the spectral log parabola model stated to obtain the corresponding flux values. During all pointings, no variability is observed in the polarization degree or angle, and the polarization and flux values are obtained averaged over the full exposures. More details on the observations, data reduction, and so on are available in Liodakis et al. (2022) and Lisalda (in prep.).

Fig. 2. Broadband SED for the three IXPE windows as described in Sect. 3.3. For IXPE-2, the X-ray spectrum corresponding to the highest flux state in the time window is shown with open markers. The typical state (Abdo et al. 2011) and low state (Abe et al. 2023) of Mrk 501 are shown in gray for comparison.

2.5. Swift-XRT

Simultaneous to the MAGIC observations, several pointing of the X-ray Telescope (XRT; Burrows et al. 2005) on board Swift (Gehrels et al. 2004) were organized. Data were taken both in the windowed timing (WT) and photon counting (PC) readout modes and processed with the XRTDAS software package (v3.7.0), which was developed by the ASI Space Science Data Center³ (SSDC), and released in the HEASoft package (v.6.30.1) by NASA High Energy Astrophysics Archive Research Center (HEASARC). For calibrating and cleaning the events, calibration files from Swift-XRT CALDB (version 20210915) are processed within the xrtpipeline.

For each observation, we extracted the X-ray spectrum from the summed cleaned event file. In WT readout mode, the spectral analysis was performed after selecting events from a circle centered at the source position with a radius of 20 pixels (∼46″; ∼90% of the point spread function). In the case of data taken in the PC mode, they are affected by pile-up in the inner part of the point spread function (PSF). Therefore, events from the most inner region around the source position within a radius of 4-6 pixels are excluded and an outer radius of 30 pixels is adopted to select signal events. We then used the shape of the Swift-XRT PSF to calculate the incident X-ray flux⁴. To estimate the background, we used nearby circular regions with radii of 20 and 40 pixels for WT and PC data. The task xrtmkarf was used to produce ancillary response files that included corrections for PSF losses and CCD defects using the cumulative exposure map. The X-ray spectra in the 0.3–10 keV range were binned with grppha with a minimum of 20 counts per bin. Afterward, we modeled them in XSPEC using power-law and log-parabola models. Energy fluxes were computed in the 0.3–2 keV and 2–10 keV bands.

2.6. Swift-UVOT

Simultaneous to the XRT pointings, the Ultraviolet/Optical Telescope (UVOT; Roming et al. 2005) also observed in the optical/UV wavebands. In this work, we considered images acquired with the three UV filters, W1 (2600 Å), M2 (2246 Å), and W2 (1928 Å), which are less affected by the host galaxy emission than the optical ones. The standard UVOT software within the HEASoft package (v.6.23) was used to perform aperture photometry for all filters and the calibration was performed using CALDB (20201026). We reviewed images to avoid unstable attitudes and then, following Poole et al. (2008), we used a circular aperture of 5″ to extract source counts and an annular aperture with 26″ (34″) for the inner (outer) radii to estimate background counts. Using the standard zero points reported by Breeveld et al. (2011), the count rates were converted to fluxes and then dereddened considering the E(B − V) value 0.017 (Schlegel et al. 1998; Schlafly & Finkbeiner 2011) for the UVOT filters effective wavelengths and the mean galactic interstellar extinction curve from Fitzpatrick (1999).

2.7. Optical

Optical data, including flux and polarization measurements, were obtained from the 2.5 m Nordic Optical Telescope (NOT), the 2.2 m telescope of the Calar Alto Observatory as part of the Monitoring AGN with Polarimetry at the Calar Alto Telescopes (MAPCAT)⁵, the 1.5 m (T150) and the 0.9 m (T090) telescopes at the Sierra Nevada Observatory, the 1.3 m RoboPol telescope (Skinakas observatory, Greece; King et al. 2014; Ramaprakash et al. 2019), the 1.5 m KANATA (Higashi-Hiroshima observatory, Japan) telescopes, and the network of robotic telescopes, the Burst Observer and Optical Transient Exploring System (BOOTES; Castro-Tirado et al. 1999). Except for BOOTES, which uses Sloan Digital Sky Survey filters, Johnson-Cousins filters are used. To convert magnitudes between the different filters, Lupton et al. (2005) was used.

Additionally, Tuorla blazar monitoring data of Mrk 501 for the period between March 1, 2022, and July 19, 2022 (MJD 59639 to MJD 59779) were obtained using the 80 cm Joan Oró Telescope (TJO) at Montsec Observatory, Spain. The observations were made in the Cousins R-filter. The flux density values were obtained following standard photometry analysis (Nilsson et al. 2018) using a signal extraction aperture of 7.5″. Additionally, to correct for the flux density contamination coming from the host galaxy and nearby sources a 12.3 mJy subtraction was implemented following Nilsson et al. (2007).

Flux and polarization degree values were corrected for the host galaxy contribution following Nilsson et al. (2007), Weaver et al. (2020), or Hovatta et al. (2016), taking the apertures of the telescopes and the seeing conditions into account. Additionally, the galaxy extinction correction of 0.041 mag was implemented following Schlafly & Finkbeiner (2011) for the flux values.

To also cover the second IXPE pointing, we added the data taken by the 70 cm AZT-8 telescope of the Crimean Astrophysical Observatory⁶, which took both photometric and polarimetric data with the Cousins R filter March 25-28, 2022 (MJD 59663 to MJD 59666; Liodakis et al. 2022) and applied the same host galaxy and the Galactic extinction correction to the optical flux. For the polarization degree, values that had already been corrected were used.

2.8. Radio

In the radio band, observations were taken on Mrk 501 by the single-dish telescopes at the Owens Valley Radio Observatory (OVRO) operating at 15 GHz, the Metsähovi Radio Observatory at 37 GHz, the Institut de Radioastronomie Millimétrique (IRAM) at 86 GHz and 230 GHz and the Submillimeter Array (SMA) at 226 GHz.

The data from OVRO were taken with the 40 m telescope blazar monitoring program using a central frequency of 15 GHz and 3 GHz bandwidth as detailed in Richards et al. (2011). 37 GHz data were taken by the Metsähovi Observatory and analyzed according to Teraesranta et al. (1998) using NGC 7027, 3C 274, and 3C 84 as secondary calibrators. Estimates on the error of the flux density contain the contribution from the measurement root mean square and the uncertainty of the absolute calibration.

Additional flux, but also polarization data, in the radio band were collected by the IRAM 30 m telescope under the Polarimetric Monitoring of AGN at Millimeter Wavelengths (POLAMI) program⁷ (Agudo et al. 2018a,b) in the 3.5 mm (86.24 GHz) and 1.3 mm (230 GHz) wavebands. The XPOL procedure (Thum et al. 2008) was used to simultaneously measure the four Stokes parameters (I, Q, U and V) at the two observing bands. The POLAMI pipeline described in Agudo et al. (2018a) was adopted to reduce and calibrate the data.

The SMA (Ho et al. 2004) was used to obtain polarimetric millimeter radio measurements at ∼1.3 mm (∼226 GHz) within the framework of the SMA Monitoring of AGNs with POLarization (SMAPOL) program. SMAPOL follows the polarization evolution of forty γ-ray loud blazars, including Mrk 501, on a biweekly cadence as well as other sources in a target-of-opportunity mode. The Mrk 501 observations reported here were conducted in June and July 2022. The SMA observations use two orthogonally polarized receivers, tuned to the same frequency range in the full polarization mode, and use the SWARM correlator (Primiani et al. 2016). These receivers are inherently linearly polarized but are converted to circular using the quarter-wave plates of the SMA polarimeter (Marrone & Rao 2008). The lower sideband (LSB) and upper sideband (USB) covered 209–221 and 229–241 GHz, respectively. Each sideband was divided into six chunks, with a bandwidth of 2 GHz, and a fixed channel width of 140 kHz. The SMA data were calibrated with the mid-infrared software package⁸. Instrumental polarization leakage was calibrated independently for the USB and LSB using the MIRIAD task gpcal (Sault et al. 1995) and removed from the data. The polarized intensity, position angle, and polarization percentage were derived from the Stokes I, Q, and U visibilities. MWC 349 A, Callisto, and Titan were used for the total flux calibration and the calibrator 3C 286, which has high linear polarization degree and stable polarization angle, was observed as a cross-check of the polarization calibration.

3. Detailed VHE and X-ray analysis during IXPE observations

The MWL LCs from March 1, 2022, to July 19, 2022 (MJD 59639 to MJD 59779) are shown in Fig. 1. The three IXPE observations conducted during this time range are highlighted with a gray band. Two observations took place March 8-10, 2022 (MJD 59646 to MJD 59648) and March 27-29, 2022 (MJD 59665 to MJD 59667). The third one took place July 9-12, 2022 (MJD 59769 to MJD 59772). These three periods are referred to as IXPE-1, IXPE-2, and IXPE-3, respectively. This section discusses the MWL states during these observations with a focus on the VHE and X-ray regimes.

3.1. IXPE observations in March 2022

The first two IXPE observations, IXPE-1 and IXPE-2 (MJD 59646 to MJD 59648 and MJD 59665 to MJD 59667), were reported in Liodakis et al. (2022). The measured polarization degree in the 2–8 keV band for those observations are 10%±2% and 11%±2%, respectively, which is around 2 times higher than the one measured at optical wavelengths. The corresponding polarization angles, 134 ± 5° and 115 ± 4°, match the optical and radio measurements and are in line with the radio jet direction of 119.7 ± 11.8° (Weaver et al. 2022, depicted as a horizontal gray band in Fig. 1). The 2–8 keV flux of the second pointing (18.0 ± 0.4 × 10⁻¹¹ erg cm⁻² s⁻¹) is a factor of 2 higher than for the first one (8.8 ± 0.1 × 10⁻¹¹ erg cm⁻² s⁻¹).

Close to the IXPE-1 epoch, four MAGIC observations are available (see Fig. 1). Two observations took place strictly simultaneous to the IXPE pointing, while the two others happened one day before and after. The nightly LC is consistent with a constant flux hypothesis, with a corresponding average flux of 6.07 ± 0.29 × 10⁻¹¹ ph cm⁻² s⁻¹ (∼30% C.U.) above 200 GeV.

Due to bad weather, no MAGIC simultaneous data are available during IXPE-2. However, one observation took place one day before the IXPE window and another one two days after (Fig. 1). Between those two MAGIC observations encompassing the IXPE-2 epoch, the MWL fluxes show little variability, in particular in the X-ray band (see Fig. 1), which is typically well correlated with the VHE emission (Acciari et al. 2020; Furniss et al. 2015; Ahnen et al. 2018; Abe et al. 2023). We thus averaged the two MAGIC exposures to obtain a reasonable approximation of the VHE state simultaneous to the IXPE-2 exposure. The resulting average flux is 9.53  ±  0.44 × 10⁻¹¹ ph cm⁻² s⁻¹ above 200 GeV (∼50% C.U.), which is ∼40% higher than during IXPE-1. Overall, the IXPE-1 and IXPE-2 VHE levels are within the typical dynamical flux range for Mrk 501 as they remain close to the average state (Abdo et al. 2011) and are persistently above its low state flux observed from 2017 to 2019 (Abe et al. 2023).

Similar to the VHE band, the Swift-XRT and IXPE measurements also reveal X-ray flux levels that are around a factor of 2 higher during the IXPE-2 epoch compared to IXPE-1. On the other hand, while the VHE emission is around the average source activity, the fluxes in the Swift-XRT bands are systematically above average. The 2–10 keV fluxes lie between ≈10⁻¹⁰ erg cm⁻² s⁻¹ and ≈3 × 10⁻¹⁰ erg cm⁻² s⁻¹, which is significantly higher than the levels reported in previous works that studied the source during an average VHE activity such as a flux of 7.8 ± 0.2 × 10⁻¹¹ erg cm⁻² s⁻¹ in 2009 (Abdo et al. 2011) or a weighted mean of 4.24 ± 0.01 × 10⁻¹¹ erg cm⁻² s⁻¹ for 12-year flux dataset from 2008 to 2020 (Abe et al. 2023). Regarding the variability, we find that the Swift-XRT data display only slight flux changes (at the level of ∼10–30%) within each IXPE window. No intra-night variability is detected, in either the VHE or the X-ray regimes.

Additionally, four NuSTAR pointings were conducted on March 9, 22, 24, and 27, 2022 (MJD 59647, MJD 59660, MJD 59662, MJD 59665). One simultaneous to each IXPE-1 and IXPE-2 window, while two other observations happened shortly before the IXPE-2 pointing. The flux ranges from 7.8 to 11.4 × 10⁻¹¹ erg cm⁻² s⁻¹ in the 3–7 keV band, which is at least a factor of 5 higher than during the previous low-activity measurements in 2017–2018 (Abe et al. 2023), but similar to the observations reported in 2013 (Furniss et al. 2015) for Mrk 501. The NuSTAR exposures, which are multi-hour-long, offer the possibility to search for short timescale variability. No strong X-ray variability is found on hourly timescale, and throughout each exposure the flux varies by at most 10%.

In the radio, optical, Swift-UVOT, and Fermi-LAT wavebands, the IXPE-1 and IXPE-2 epochs are in a rather averaged and stable state (Abdo et al. 2011). No significant flux variations on daily timescales are detected in these bands.

3.2. IXPE observation in July 2022

The third IXPE pointing, IXPE-3 (MJD 59769 to MJD 59772), is characterized by a slightly lower flux level (7.4 ± 0.4 × 10⁻¹¹ erg cm⁻² s⁻¹ in the 2–8 keV band) and polarization degree (7 ± 2%) in the X-rays than for the IXPE-1 or IXPE-2 pointings (Lisalda, in prep.). For MAGIC, no strictly simultaneous data are available, but a single observation took place less than one day before IXPE-3. It shows a lower flux state compared to the ones in March with a flux of 4.36  ±  1.45 × 10⁻¹¹ ph cm⁻² s⁻¹ above 200 GeV (∼20% C.U.). The X-ray fluxes of Swift-XRT are lower than during IXPE-2, although comparable to the ones measured during IXPE-1.

3.3. Spectral evolution during IXPE campaigns

Table 1 reports the fluxes and spectral parameters during each IXPE epoch for the MAGIC, Fermi-LAT, NuSTAR and Swift-XRT observations. For all instruments except NuSTAR, we fitted a simple power-law model given the absence of significant evidence for spectral curvature. The NuSTAR spectra show a significant curvature and are thus fitted with a log-parabola model, $\frac{\mathrm{d}N}{\mathrm{d}E}\propto {(E/E_{\mathrm{norm}})}^{-\alpha -\beta log(E/E_{\mathrm{norm}})}$ (with a normalization energy of E_norm = 1 keV). Furthermore, we fixed the curvature parameter to β = 0.3 (average value of the pointings) to remove any correlation between α and β and obtain a better assessment of the spectral hardness evolution. For each instrument, we state in the corresponding panels the time ranges over which the spectra are computed. For Swift-XRT, the best-fit parameters are shown for each observation separately given the indication of daily spectral evolution.

Comparing the three IXPE epochs, one can see the higher and harder emission for IXPE-2, in particular in the X-ray band. In fact, IXPE-2 is characterized by one of the hardest X-ray states from the campaign discussed in this work (α <  1.9). Such a hard state is not evident in the VHE band, which in fact shows a constant power-law index between the three epochs within uncertainties (α ∼ 2.4 − 2.5; see Table 1). The Swift-XRT spectra show a “harder-when-brighter” trend that is also observed when comparing the NuSTAR pointings. In the Fermi-LAT band, the spectral index softens along the three IXPE time windows, leading to a decrease in energy flux despite increasing photon flux levels.

We constructed MWL SEDs for all pointings. Except for the γ-ray regime, we used the weighted average of all data points in the corresponding IXPE time windows as stated in the first row of Table 1. The only exception was made for the IXPE-2 spectrum, for which the dedicated NuSTAR observation started shortly before the IXPE time window. Therefore, the NuSTAR observational window on March 27, 2022 (MJD 59665) was used to extract the weighted averaged SED points for the radio, UV, and X-ray bands. For the optical regime, the data simultaneous to IXPE-2 from Liodakis et al. (2022) were used because they are the only data available around this time window. This ensures an easy comparison between the two X-ray instruments and therefore allows us to characterize the synchrotron peak. For comparison the Swift-XRT spectrum corresponding to the highest flux level during the IXPE-2 window is shown with open markers in Fig. 2. For Fermi-LAT a time window of two weeks centered on the IXPE epochs was used to construct the SEDs owing to the low variability in the band and the lower instrument sensitivity, which requires longer observation times. For MAGIC, simultaneous data to the IXPE observations are only available for IXPE-1. We used all simultaneous nights and also the two surrounding ones to construct the IXPE-1 MAGIC SED because no flux or spectral variations are seen across the different days. Due to the availability of VHE observations as discussed in Sects. 3.1 and 3.2 and the absence of strong spectral and flux variability at other wavebands, the observations surrounding the IXPE windows are used to construct the IXPE-2 and IXPE-3 VHE SEDs.

The resulting SEDs are shown in Fig. 2. For comparison purposes, we plot in light gray the average source state from Abdo et al. (2011), and in dark gray the low activity discussed in Abe et al. (2023). As anticipated in Sect. 3.1, Fig. 2 highlights a drop in the Compton dominance (CD; the ratio between the νF_ν values at the high-energy and low-energy peaks) with respect to the average activity (light gray points). The entire γ-ray SED component remains close to the average activity, while the X-ray emission is significantly higher.

Using the phenomenological model of Ghisellini et al. (2017), we quantified the peak frequencies and CD for each epoch. The same phenomenological model was fitted to the average and low states from Abdo et al. (2011) and Abe et al. (2023) for comparison purposes. The results are shown in Table 2. The obtained best-fit value place both IXPE-1 and IXPE-2 in the EHSP (Costamante et al. 2001; Biteau et al. 2020) family with ν_s >  2.4 × 10¹⁷ Hz (≈1 keV).

For the third pointing, one can clearly see that the synchrotron peak moved to lower frequency, into the HSP state as verified by the phenomenological fits (see Table 2) yielding a ν_s of ∼1 × 10¹⁷ Hz. While the MAGIC spectrum shows only a change in amplitude (without spectral evolution), the Fermi-LAT spectrum hints at softer emission, leading to a shift in the high-energy peak toward lower frequencies by an order of magnitude.

For all three spectra, the fit of the phenomenological model confirms a general lower CD than using the average state of Mrk 501, which has a CD around ∼0.5. For IXPE-1 and IXPE-2, the CD was determined to be ∼0.3 for both the phenomenological model and the synchrotron-self-compton (SSC) model described in Sect. 5. For IXPE-3, it was determined to be 0.2 and 0.3 using the phenomenological model and the SSC model, respectively.

3.4. Evolution of the polarization degree between the IXPE epochs

Figure 3 shows the polarization degree of all three IXPE pointings in different wavebands, from the radio to the X-ray band. Additionally, we show in the middle panel the ratio of the polarization degree to the one observed in the X-rays. For all three epochs, the polarization degree increases with frequency showing a ratio that is slightly more than two between the X-ray and optical regimes. This behavior was first reported for the IXPE-1 and IXPE-2 epochs by Liodakis et al. (2022). IXPE-3 follows the same trend, but it is interesting to see that while the absolute polarization degree drops with respect to IXPE-1 and IXPE-2 in optical and X-rays, the ratio is roughly conserved and stable over time (the derived ratios remain consistent within ∼1σ).

Fig. 3. MWL polarization and flux properties of Mrk 501 during the three IXPE observations. Top: MWL polarization degree as a function of frequency. For the radio and optical data, weighted averages over the IXPE time windows are used. Middle: ratio of the frequency-dependent polarization degree to the corresponding X-ray polarization degree. Bottom: flux values associated with the polarization values presented above.

As described in Liodakis et al. (2022) and Di Gesu et al. (2022), the increase in the polarization degree with energy, the absence of fast variability in the polarization behavior as well as the alignment of the radio jet with the polarization angle in optical and X-rays strongly suggest that electrons are accelerated in shocks. The most energetic electrons (those emitting X-ray photons) are confined near the shock front, where the magnetic field is more ordered, leading to a stronger polarization degree. The electrons subsequently cool, emit lower energy photons, and advect (or diffuse) toward larger regions in which the degree of magnetic field ordering drops significantly. Such an energy-stratified scenario qualitatively explains the broadband polarization degree evolution (Tavecchio et al. 2020), and will be considered for the MWL modeling of the SEDs presented in Sect. 5.

The lower panel of Fig. 3 shows the average flux states for the three epochs in each energy bands. No obvious coherent trend is apparent between the evolution of the polarization degree and the flux level. In fact, while in the optical the polarization degree tends to be anticorrelated with the flux, the X-ray hints at an opposite trend.

4. Long-term multiwavelength evolution

This section discusses the MWL evolution throughout the entire period between March 1, 2022, and July 19, 2022 (MJD 59639 to MJD 59779) to provide a broader overview of the source evolution between the various IXPE epochs.

From Fig. 1, the VHE flux shows no significant outbursts and varies around ∼20–50% that of the Crab Nebula, which is close to the typical state of Mrk 501 (Abdo et al. 2011). The low (0.2–1 TeV) and high (> 1 TeV) energy ranges follow similar variability patterns. No harder-when-brighter trend is found; this, however, might be due to the lower instrumental sensitivity of MAGIC to slight spectral changes compared to, for example, X-ray instruments.

The Swift-XRT bands reveal clear hardening with increasing flux. In Fig. 4, we show the hardness ratio (defined as the ratio of the 2–10 keV flux to the one in the 0.3–2 keV band) as a function of the 2–10 keV flux. The gray data points are from individual measurements, while the black points are the hardness ratio binned in flux. The error bars represents the standard deviation of the hardness ratio in each bin. A positive slope of 0.032  ±  0.002 is found when fitting the individual data points with a linear function. For the hardness ratio binned in flux, the best fit slope is comparable, 0.024  ±  0.003. These slopes are significantly lower than what was obtained in Abe et al. (2023) during a much lower activity of Mrk 501, possibly indicating a saturation of the hardness ratio above a certain flux level similarly to what was noted in Mrk 421 in Acciari et al. (2021).

Fig. 4. Swift-XRT hardness ratio (flux ratio of the 2–10 keV range over the 0.3–2 keV band) compared to flux at 2–10 keV for the LCs shown in Fig. 1. The gray data points are from individual measurements, and the black points are the hardness ratio binned in flux.

For the UV band, we see a change in flux happening in April that slightly increases the average UV flux level; this can be explained by the shift in the synchrotron peak to lower frequencies. This explanation is also suggested by the comparison of the three IXPE SEDs in Fig. 2, which show a mild increase in the optical/UV during the IXPE-3 period and thus a shift in the synchrotron component toward lower energies.

We quantified the variability strength throughout the spectrum using the fractional variability (F_var) defined following Eq. (10) in Vaughan et al. (2003):

$$\begin{array}{c}\hfill F_{\mathrm{var}}=\sqrt{\frac{S^2-⟨{\sigma }_{\mathrm{err}}^2⟩}{{⟨F_{\gamma }⟩}^2}},\end{array}$$(1)

where S² is the variance of the signal, $⟨{\sigma }_{\mathrm{e}\mathrm{r}\mathrm{r}}^2⟩$ the mean square error of the measurement uncertainties, and ⟨F_γ⟩² the arithmetic mean of the measured flux. Its uncertainty is computed following Poutanen et al. (2008) as

$$\begin{array}{c}\hfill \mathrm{\Delta }{F}_{\mathrm{var}}=\sqrt{F_{\mathrm{var}}^2+\mathrm{err}({\sigma }_{\mathrm{NXS}}^2)}-F_{\mathrm{var}},\end{array}$$(2)

with

$$\begin{array}{c}\hfill \mathrm{err}({\sigma }_{\mathrm{NXS}}^2)=\sqrt{{\left(\sqrt{\frac{2}{N}}\frac{⟨{\sigma }_{\mathrm{err}}^2⟩}{{⟨F_{\gamma }⟩}^2}\right)}^2+{\left(\sqrt{\frac{⟨{\sigma }_{\mathrm{err}}^2⟩}{N}}\frac{2F_{\mathrm{var}}}{⟨F_{\gamma }⟩}\right)}^2}.\end{array}$$(3)

Considering the whole time epoch from March 1 to July 19, 2022 (MJD 59639 to MJD 59779), all wavebands show small variability scales. As shown in Fig. 5 and reported in Table A.1, F_var is always below 0.3 for the full time epoch (closed markers). The high-energy X-rays and VHE γ-rays display the highest degrees of variability. As reported with data from other multi-instrument campaigns, Metsähovi systematically shows a higher F_var than OVRO (Aleksić et al. 2015; Furniss et al. 2015; Ahnen et al. 2018; Madsen et al. 2022); however, we also find a higher degree of variability of the 230 GHz IRAM data compared with the 86 GHz data. The high variability in the 37 GHz is not caused by a flaring event, but rather by a flickering in the radio fluxes, which may perhaps be produced by the somewhat limited sensitivity of Metsähovi, and the relatively low brightness of Mrk 501 at radio. To check that the higher variability degree in Metsähovi is not caused by a threshold effect, we computed it using different signal quality cuts, namely S/Ns >  5, 7, and 10. We found that, for all the abovementioned signal qualities, the computed F_var values are consistent with the one computed with the full dataset (S/N >  4), which is higher than the F_var obtained for low-frequency radio instruments. Additionally, we investigate the variability for only the IXPE time windows shown by the turquoise markers in Fig. 5. Due to a lack of simultaneous data, it is not computed for the γ-ray data. For the X-rays, a higher degree of variability is found than in the full time period, while for the lower wavebands (radio to UV), a very similar behavior across all time epochs is found.

Fig. 5. Fractional variability for the LCs displayed in Fig. 1. Nightly bins are used for MAGIC and 7-day bins for Fermi-LAT; for all other wavebands, the single observations without further binning are used. For the radio data, frequencies are stated to distinguish between the different datasets (from left to right: the orange circles depict OVRO at 15 GHz and Metsähovi at 37 GHz, the lighter orange square the IRAM data at 86 GHz, the orange diamond SMA at 226 GHz, and the darker orange square and open turquoise square IRAM at 230 GHz). The open purple markers refer to the NuSTAR (square) and IXPE (round) observations, which only consist of two or three measurements and therefore far fewer data points than considered for the full time epoch. Due to the different sampling, not all instruments are directly comparable with each other.

We searched for intra-band correlation and found a hint of positive correlation only between the VHE and X-ray LCs, which also correspond to the energy bands with the strongest variability and the locations of the two SED peaks. The correlation is quantified using the discrete correlation function (DCF; Edelson & Krolik 1988), which we computed between the MAGIC fluxes above 200 GeV and the 0.3–2 keV and 2–10 keV fluxes from Swift-XRT. The DCF was calculated for a series of time lags from −20 days to +20 days between the LCs using time lag steps of 5 days, which is in agreement with the timescale taken from the auto-correlation behavior of Mrk 501 during low activity in Abe et al. (2023). The significance of the correlation was estimated based on Monte Carlo simulations in a similar fashion as the one described in MAGIC Collaboration (2021). We simulated 10⁴ pairs of uncorrelated LCs that respect the sampling and binning of the real data. The DCF significance band in each time lag was then obtained from the distribution of DCF values of the simulations. To simulate the LCs, we assumed that for each energy band the power spectral density function used as an input for the simulation follows a power-law model with an index of 1.4, which is the same index derived in Aleksić et al. (2015). As shown in Figs. 6 and 7, the significance is slightly above 2σ at zero time lag.

Fig. 6. Correlation between the MAGIC flux above 200 GeV and the Swift-XRT flux, 0.3–2 keV. See the main text for further details.

Fig. 7. Correlation between the MAGIC flux above 200 GeV and the Swift-XRT flux, 2–10 keV. See the main text for further details.

In summary, we find that during the full campaign Mrk 501 has showed very typical behaviors considering its MWL flux values. It does not just stay around its average state (Abdo et al. 2011), but also depicts the highest F_var in the X-rays and VHE as previously found in Furniss et al. (2015) and Abe et al. (2023) and shows evidence for the commonly observed correlation between the X-ray and VHE band (Acciari et al. 2020; Furniss et al. 2015; Ahnen et al. 2018; Abe et al. 2023).

As for the long-term evolution in the polarization, we note an increase in the polarization degree on a few day timescale in the R band around 2022-05-21 (MJD 59720). The R-band polarization increases by a factor of ∼2 and reaches a level of 12%, which is above the degree in the X-ray during the IXPE epochs. No significant flare is observed simultaneous to this elevated optical polarization state. The lack of simultaneous MWL polarization data prevents us from making more detailed investigations. As for the polarization angle, moderate variability by a few tens of degrees can be seen in the radio and optical bands. A significant offset between the R-band and radio (SMA; 226 GHz) angle temporarily occurs from 2022-05-31 (MJD 59730) to 2022-06-30 (MJD 59760). The angle deviates by ∼60° around 2022-05-31 (MJD 59730). This behavior contrasts with the one during the IXPE epochs where a rough alignment (within ∼30°) of the optical/radio/X-ray polarization angle can be noticed.

5. Modeling of the SEDs during the IXPE observing windows

We exploited the broadband SEDs shown in Fig. 2 to perform a modeling of the emission assuming a leptonic scenario. We considered an energy-stratified emission region, as suggested by the optical-to-X-ray polarization properties (see Sect. 3.4). The region consists of two overlapping blobs, each of them with different radius R′. The two blobs contribute dominantly to distinct energy regions of the SED. The X-ray and > 100 GeV γ-ray fluxes are dominated by a “compact zone”, located close to the shock. The second zone, dubbed as “extended zone”, contributes to the radio, optical/UV and MeV-GeV fluxes. The extended zone is significantly larger than the compact zone, and expands over a larger volume downstream from the shock front. The theoretical framework that we used here is in line with the energy-stratified model of Liodakis et al. (2022), Angelakis et al. (2016), Itoh et al. (2016), where the X-ray emission is produced closer to the shock front than the optical emission. In Fig. 8 we show a sketch depicting the jet morphology of our modeling framework. While considering two distinct zones is surely a simplification of the system (one would rather expect a continuum of several regions contributing to the different parts of the SED), a more realistic description of the system would require a much more detailed treatment of particle cooling and diffusion/advection that is beyond the scope of this work. This simplified model already allows one to constrain important properties of the emission states of Mrk 501 during the three IXPE pointings.

Fig. 8. Simplified sketch representing the morphology of the two-zone leptonic model discussed in Sect. 5.

In order to limit the degrees of freedom, we made the following assumptions. Parameters referred to the blob frame are marked as primed, while unmarked ones relate to the observer frame.

• The Doppler factor, δ, was fixed to 11 consistent with the one used in Abe et al. (2023). For simplicity, the same value was chosen for both regions.

• In order to estimate the relative size between the two zones, we proceeded as follows: we first assumed that the compact region can be modeled by N turbulent plasma cells of identical magnetic field strength, but with random orientation. In such a configuration, the expected average polarization degree from the zone can be approximated as $P_{\text{deg}}\approx 75\%/\sqrt{N}$ (Marscher 2014; Tavecchio 2021). Given that IXPE-1 and IXPE-2 are characterized by P_{deg, X-ray} ≈ 10% in the X-ray band, this implies N ≈ 60 during these two epochs. Regarding IXPE-3, P_{deg, X-ray} ≈ 7%, corresponds to N ≈ 110. The size of the extended zone, which dominates the optical emission, was then assumed to be a factor of l greater than the compact region such that the number of turbulent cells in the extended zone matches the optical polarization degree according to $P_{\text{deg}}\approx 75\%/\sqrt{lN}$. For IXPE-1 and IXPE-2, P_{deg, opt.} ≈ 5%, yielding l ≈ 5. For IXPE-3, P_{deg, opt.} ≈ 2%, which implies l ≈ 10. Assuming that turbulent cells roughly span equal volumes, the radius of the extended zone should be a factor of l^1/3 larger than the one of the compact zone. For all epochs, we set the radius of the compact zone to R′=2.9 × 10¹⁶ cm. This is in agreement with a light crossing time R′/(cδ)≈1 day, which is the variability timescale observed in the X-ray band.

• The electron distribution inside the compact zone was modeled with a broken power-law function normalized such that the electron energy density is given by $U_{\mathrm{e}}^{\mathrm{\prime }}$. The distribution is defined between a minimum and maximum Lorentz factor of ${\gamma }_{\text{min}}^{\mathrm{\prime }}$ and ${\gamma }_{\text{max}}^{\mathrm{\prime }}$. The break is located at ${\gamma }_{\mathrm{b}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{k}}^{\mathrm{\prime }}$. The slopes below and above the break are given by n₁ and n₂, respectively. ${\gamma }_{\text{min}}^{\mathrm{\prime }}$ was set to a value such that the compact zone remains subdominant in the optical/UV. With the choice of δ and R′ mentioned above, we derive ${\gamma }_{\text{min}}^{\mathrm{\prime }}\sim {10}^4$. Within an ion-electron plasma, mildly relativistic shocks, in which the shock front is moving at a Lorentz factor γ_sh ∼ 1 − 3 relative to the un-shocked plasma, are expected to generate somewhat lower ${\gamma }_{\text{min}}^{\mathrm{\prime }}$, on the order of a few times 10³ (Zech & Lemoine 2021). Values of ${\gamma }_{\text{min}}^{\mathrm{\prime }}\sim {10}^4$ may be reached for fully relativistic shocks with γ_sh of a few tens. In order to constrain ${\gamma }_{\text{max}}^{\mathrm{\prime }}$, we equated the acceleration and cooling timescales. The acceleration timescale in shock acceleration was estimated as

  $$\begin{array}{c}\hfill t_{\mathrm{acc}}^{\prime }=\frac{20\lambda ({\gamma }^{\prime })c}{3u_{\mathrm{s}}^2},\end{array}$$(4)

  where λ(γ′) = (ξγ′m_ec²)/(eB′) is the mean free path of electrons, parameterized as a fraction, ξ, of the Larmor radius. The ξ is an acceleration parameter, which we fixed to 10⁴. This value is in agreement with previous estimates from Zhang (2002) based on spectral hysteresis observations of a similar HSP, Mrk 421. B′ is the magnetic field strength, e the electron charge, m_e the electron mass, and u_s the speed of the shock, which we assumed to be relativistic, u_s ∼ c. The synchrotron cooling timescale is

  $$\begin{array}{c}\hfill t_{\mathrm{cool},\mathrm{synch}}^{\prime }=\frac{6\pi m_{\mathrm{e}}c}{{\sigma }_T{B}^{\prime 2}{\gamma }^{\prime }},\end{array}$$(5)

  where σ_T is the Thomson cross section.

• The electron distribution inside the extended zone was modeled with a simple power-law function because the available data prevent a precise determination of breaks or other distributions with more degrees of freedom. The slope of the distribution was assumed to be p = 2.2. The minimum and maximum Lorentz factor of the distribution, ${\gamma }_{\text{min}}^{\mathrm{\prime }}$ and ${\gamma }_{\text{max}}^{\mathrm{\prime }}$ are constrained by the radio-to-UV data. The distribution is normalized to an energy density of $U_{\mathrm{e}}^{\mathrm{\prime }}$. We kept ${\gamma }_{\text{max}}^{\mathrm{\prime }}$ within a factor of 2 of the expected cooling break that is given by equating the cooling timescale due to synchrotron and inverse-Compton processes with the escape timescale R′/c.

We started by first fitting the X-ray and γ-ray data to determine the remaining free parameters of the compact zone (i.e., magnetic field B′, $U_{\mathrm{e}}^{\mathrm{\prime }}$, n₁, n₂). In a second step, the extended zone was added in order to get a good description of the radio-to-UV and Fermi-LAT data. Between the IXPE-1 and IXPE-2 epochs, the data are satisfactorily described with the same parameters for the extended zone, which is not surprising given the lower variability in the radio, optical/UV and Fermi-LAT data (see Figs. 1 and 5), and the relatively short time period between IXPE-1 and IXPE-2 of a few weeks, while IXPE-3 took place several month later.

The adopted model is shown in Fig. 9 for each epoch. The dashed blue lines and dashed-dotted violet lines are the contributions from the compact zone and extended zone, respectively. The solid light blue line is the sum of the two regions. We also computed the interaction between the two zones (dotted pink line), which results from electron inverse-Compton scattering off the photon field that the two regions feed to each others. The derived parameter values are shown in Table 3. Overall, the model describes all three broadband SEDs well. A more detailed interpretation of the results is given in Sect. 6.

Fig. 9. Modeling of the three IXPE epochs: IXPE-1 (top), IXPE-2 (middle), and IXPE-3 (bottom). The dashed blue curve is the emission originating from the compact zone, located near the shock front. The dashed-dotted violet curve is the contribution from the extended zone, which spans a larger volume downstream from the shock. The emission resulting from the interaction between the two zones is plotted as a dotted pink curve. The sum of all components is given by the solid blue line. The light gray and dark gray markers show the average state (from Abdo et al. 2011) and low state (from Abe et al. 2023), respectively. The parameters of the model are listed in Table 3.

6. Summary and discussion

In this work, we present the MWL picture around the first X-ray polarization measurements of Mrk 501. For the first time, simultaneous information in the VHE regime is published and compared to the X-ray measurements, allowing the first broadband SEDs up to the VHE band to be constructed, covering the three IXPE time windows (March 8-10, 2022, March 27-29, 2022, and July 9-12, 2022, i.e., MJD 59646-59648, MJD 59665-59667, and MJD 59769-59772).

Both the VHE and X-ray regimes showed the highest flux level during IXPE-2 (a factor of ∼2 higher than during the other IXPE windows). Comparing the spectral properties, a spectral hardening with increasing flux levels is observed in the X-rays between the different pointings, which is supported by the harder-when-brighter trend observed over the full campaign (see Fig. 4). The opposite behavior can be seen in the high-energy γ-ray data, where the spectral index softens with higher flux levels. No spectral change in the VHE is seen between the three pointings. Interestingly, the power-law indices in the X-rays are rather hard during IXPE-1 and IXPE-2 (≲2), as is typically observed for elevated flux states, while the VHE spectral hardness remains close to the average of Mrk 501 activity (power-law index of ∼2.5; Abdo et al. 2011).

Noteworthy is not just the change in spectral index, but also the broadband characteristics of the different SEDs (see Table 2), which revealed an EHSP behavior for Mrk 501 in the IXPE-1 and IXPE-2 states, with a synchrotron peak frequency well beyond 1 keV. During IXPE-3, the blazar went back to its more typical HSP behavior, shifting the synchrotron peak by an order of magnitude toward lower energies. The shift in synchrotron peak frequencies can explain the positive change in flux observed in the UV band of the MWL LC in April (see Fig. 1). The accompanying shift in the high-energy peak can explain the different flux behavior of the high-energy γ-rays, which increases steadily between the three states, compared to the rise and then fall of the fluxes in the VHE and X-rays. During all three states, the CD is lower by a factor of 2 (see Table 2) compared to the typical state of Mrk 501 (Abdo et al. 2011). We show the same parameters obtained from the leptonic modeling in Table 4 for comparison. While the parameters for IXPE-1 and IXPE-2 match very well, both peak frequencies shift toward lower energies for IXPE-3. However, when considering the increased uncertainties of the parameters for the IXPE-3 state, a clear shift can only be claimed for the synchrotron peak. This shows the poorer coverage of this emission state in our data, and highlights that the perceived drop in CD for IXPE-3, as recorded in Table 2, should be considered with particular caution due to the dependence on the peak frequency.

Mrk 501 previously showed EHSP behavior during periods of high flaring activity in 1997 (Tavecchio et al. 2001) and 2013 (Furniss et al. 2015). This has also been seen in other HSP blazars (MAGIC Collaboration 2020b,a). Additionally, similarly to the HSP blazar 1ES 2344+514 (Abe et al. 2024), Mrk 501 has shown EHSP behavior during non-flaring activity as well, as reported in Ahnen et al. (2018) based on the extensive multi-instrument campaign carried out in 2012. Compared to the data from 2022 reported in this manuscript, the abovementioned studies from past campaigns showed a far higher variability in both VHE and X-rays. While similar average flux levels are found in this work and the X-rays in 2012, the average VHE flux levels were close to the typical flux levels of Mrk 501 during the EHSP behavior. However, both bands depict larger flux changes in 2012 than in 2022, together with higher degrees of variability, reaching up to F_var ∼ 1 in the VHE. The main difference in behavior in 2022 compared to that reported in 2012–2013 is that only the X-rays show a hard spectrum, while in 2012–2013 both X-rays and VHE revealed atypically hard SEDs.

In addition to the shift in peak frequencies, the polarization degree drops for the IXPE-3 pointing, both in the X-rays and optical regime. However, the ratio between the optical and X-ray band polarization stays rather constant during all three pointings (see Fig. 3). In fact, the X-ray polarization is systematically a factor of ∼2 − 3 higher than in the optical band for each IXPE epoch. No obvious trend between the changes in flux in the corresponding regimes and the change in polarization is seen. Therefore, the spectral changes and broadband SEDs more probably explain the differences in polarization degree between the different states, as discussed below using our theoretical modeling.

The polarization angle during the IXPE epochs shows a rough alignment between the radio, optical, and X-ray bands (within ∼30°). The values are close to the average jet angle of 119 ± 12° determined by Weaver et al. (2022) using Very Long Baseline Array imaging at 43 GHz (see Fig. 1, bottom panel). Nonetheless, in the middle of the campaign, both the radio (SMA; 226 GHz) and optical (R-band) polarization angles show some variability and deviate from the jet angle in an incoherent manner. Indeed, around May 31, 2022 (MJD 59730), the offset between the radio and optical bands reaches ∼60°. This behavior is very suggestive of a separation (at least partially) between the optical and radio regions. The deviation of the polarization angle from the jet axis may arise if the emitting regions are located in a portion of the jet where bending on a short timescale occurs (Lyutikov et al. 2005; Myserlis et al. 2018; Tavecchio 2021). Furthermore, it points to multiple emission regions being responsible for the long-term behavior of Mrk 501, with a stable part of the magnetic field regulated by shocks but at times dominated by more variable and/or turbulent contributions, as proposed in Abe et al. (2023).

Over the full time period, the blazar showed a rather stable behavior with typical flux levels (Abdo et al. 2011) in the VHE regime (∼20–50% C.U.). In contrast, the XRT flux values are constantly above 10⁻¹⁰ erg cm⁻² s⁻¹. Based on the long-term LC presented in Abe et al. (2023), this indicates an elevated state with respect to the average X-ray activity (Abdo et al. 2011). It confirms that the atypically low CD during the IXPE windows is a feature present throughout the entire campaign presented in this work. For all bands, the fractional variability stays below F_var = 0.3 when considering the full time epoch, with slightly higher variability in the X-rays for only the IXPE time windows (see Fig. 5). Even though different behaviors are reported in the X-rays and VHE γ-rays, we identify evidence for a correlation between the two bands without a time lag (see Figs. 6 and 7). The correlations are found with significances of between 2σ and 3σ; this is limited by the short time epoch considered, the low flux levels, the very low variability, and the somewhat limited sensitivity to flux changes in the VHE γ-ray band in comparison with the X-ray band.

To further investigate the emission states during the first IXPE pointings, we modeled the broadband SEDs with a two-zone leptonic model (see Sect. 5). We assumed that a compact region, which is near the shock front and populated by freshly accelerated electrons, is responsible for the X-ray and VHE γ-ray band. The compact region is embedded into an extended region that expands downstream from the shock. The extended region dominates the emission in the radio, optical, UV, and Fermi-LAT bands. Such a morphology of the emitting region is supported by the observed energy dependence of the polarization degree between the radio and the X-rays (Liodakis et al. 2022), which points to an energy-stratified jet. In addition, previous long-term correlation studies (e.g., Abe et al. 2023) point to the same emitting region dominating the X-ray and γ-ray bands. While the extended zone is assumed to be constant across the three pointings, the compact one varies over time, explaining the observed spectral changes.

A good description of the observations from the radio to VHEs is achieved for each epoch. The X-ray data from IXPE-1 and IXPE-2 constrain both the rising and falling edges of the synchrotron SED close to the peak energy. This in turn provides strong constraints on the spectral shape of the electron distribution within the region near the shock front. Below the break, which results from either cooling or acceleration effects, the derived index of the electron distribution is n₁ = 2.25 for the compact zone. This is in close agreement with expectations from diffusive acceleration in relativistic shocks (Kirk et al. 2000). Regarding the IXPE-3 epoch, similar n₁ is obtained, although the n₁ is less constrained given the shift in the synchrotron peak toward lower energies.

The emission resulting from the interaction between two zones contributes significantly to the Fermi-LAT band. Consequently, the X-rays should not only correlate with the VHE fluxes (as reported several times in the past for Mrk 501), but also with the one from Fermi-LAT. Using 12 years of data, Abe et al. (2023) detected a positive correlation between the Swift-XRT and Fermi-LAT emission, which is thus in good agreement with our model.

Across the different epochs, most of the parameters of the compact zone are quite similar. A main difference is the drop in ${\gamma }_{\mathrm{b}\mathrm{r}}^{\mathrm{\prime }}$ by a factor of ≈4 during IXPE-3 with respect to IXPE-1 and IXPE-2, which explains the shift in the synchrotron component toward lower energies. The obtained ${\gamma }_{\mathrm{b}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{k}}^{\mathrm{\prime }}$ for IXPE-1 and IXPE-2 are a factor of 4 higher than that expected for a self-consistent equilibrium between the cooling due to synchrotron and inverse-Compton processes (injection and escape; Tavecchio et al. 1998). For IXPE-3, the difference relaxes to a factor of 2. Additionally, the magnetic field, B′, increases by ≈40% during IXPE-3, from 0.050 G to 0.068 G.

The drop in the optical/X-ray polarization degree during the IXPE-3 epoch compared to the earlier IXPE epochs suggests a general evolution of the shock properties. The increase in B′ in the compact zone for IXPE-3 points to an evolution of the magnetization. Within a multi-cell scenario, the shift in peak frequencies can be attributed to an increase in synchrotron cooling due to a stronger magnetic field. At the same time, the drop in polarization implies that the magnetic field becomes less ordered and that the emitting region is composed of a larger number of turbulent plasma cells. It has been shown before that turbulences created by relativistic shock propagation could lead to an amplification of the magnetic field (Mizuno et al. 2014). Furthermore, as discussed in Kirk et al. (2000) and Komissarov & Lyutikov (2011), a higher magnetization can lead to a decrease in the shock compression ratio, κ (defined as the ratio between the upstream and downstream plasma velocity in the shock reference frame). Since the polarization degree is positively correlated with κ, an increase in the magnetization caused by a more turbulent plasma can be expected to result in an overall decrease in the polarization degree. Such a scenario may thus explain why IXPE-3 is characterized by a larger number of turbulent cells compared to IXPE-1 and IXPE-2, as pointed out in Sect. 5, where we derive ∼60 cells for IXPE-1 and IXPE-2 and ∼110 cells for IXPE-3 in the compact zone.

For simplicity, the radius of the compact zone is the same for all epochs in our model, and is set to the highest value allowed by the daily timescale variability observed in the X-rays. However, one may argue that a higher number of turbulent cells during IXPE-3 actually suggests a larger emitting zone with respect to the other epochs. If one assumes a stable magnetic field of 5 × 10⁻² G in the compact zone throughout the observing campaign (which is the value derived for IXPE-1 and IXPE-2), we find that the radius of the compact zone would have to increase by ∼20% with respect to IXPE-1 and IXPE-2 to properly describe the SED of IXPE-3. Such a relative radius increase leads to 1.7 times more turbulent cells in the emitting region, under the assumption that turbulent cells have a roughly constant volume over time. This is in good agreement with the fact that the IXPE-3 compact zone should be composed of approximately twice the number of turbulent cells compared to IXPE-1 and IXPE-2 (see Sect. 5). In summary, based on our modeling, the decrease in the polarization degree throughout the IXPE epochs may be qualitatively explained by a change in the magnetization and/or a change in the emitting region size.

In conclusion, with this 2022 dataset we could, for the first time, combine broadband MWL data up to the VHE band together with polarization measurements up to the X-rays, which allowed us to better constrain the emission and acceleration responsible for the observed radiation. This shows the crucial role of further MWL monitoring in combination with coverage of the X-ray polarization of different emission states of Mrk 501 to extend our knowledge on the mechanisms behind the blazar’s emission. Further insights could be gained from polarization measurements at even higher energies that cover the high-energy peak of the SED. In this context, e-ASTROGAM (De Angelis et al. 2017), COSI (Tomsick et al. 2021), or AMEGO (McEnery et al. 2019), which are being constructed or considered for construction in the next few years, could provide even more vital insights.

Acknowledgments

The MAGIC Collaboration would like to thank the Instituto de Astrofísica de Canarias for the excellent working conditions at the Observatorio del Roque de los Muchachos in La Palma. The financial support of the German BMBF, MPG and HGF; the Italian INFN and INAF; the Swiss National Fund SNF; the grants PID2019-104114RB-C31, PID2019-104114RB-C32, PID2019-104114RB-C33, PID2019-105510GB-C31, PID2019-107847RB-C41, PID2019-107847RB-C42, PID2019-107847RB-C44, PID2019-107988GB-C22 funded by the Spanish MCIN/AEI/ 10.13039/501100011033; the Indian Department of Atomic Energy; the Japanese ICRR, the University of Tokyo, JSPS, and MEXT; the Bulgarian Ministry of Education and Science, National RI Roadmap Project DO1-400/18.12.2020 and the Academy of Finland grant no. 320045 is gratefully acknowledged. This work was also been supported by Centros de Excelencia “Severo Ochoa” y Unidades “María de Maeztu” program of the Spanish MCIN/AEI/ 10.13039/501100011033 (SEV-2016-0588, CEX2019-000920-S, CEX2019-000918-M, CEX2021-001131-S, MDM-2015-0509-18-2) and by the CERCA institution of the Generalitat de Catalunya; by the Croatian Science Foundation (HrZZ) Project IP-2016-06-9782 and the University of Rijeka Project uniri-prirod-18-48; by the Deutsche Forschungsgemeinschaft (SFB1491 and SFB876); the Polish Ministry Of Education and Science grant No. 2021/WK/08; and by the Brazilian MCTIC, CNPq and FAPERJ. The Imaging X-ray Polarimetry Explorer (IXPE) is a joint US and Italian mission. The US contribution is supported by the National Aeronautics and Space Administration (NASA) and led and managed by its Marshall Space Flight Center (MSFC), with industry partner Ball Aerospace (contract NNM15AA18C). The Italian contribution is supported by the Italian Space Agency (Agenzia Spaziale Italiana, ASI) through contract ASI-OHBI-2017-12-I.0, agreements ASI-INAF-2017-12-H0 and ASI-INFN-2017.13-H0, and its Space Science Data Center (SSDC), and by the Istituto Nazionale di Astrofisica (INAF) and the Istituto Nazionale di Fisica Nucleare (INFN) in Italy. This research used data products provided by the IXPE Team (MSFC, SSDC, INAF, and INFN) and distributed with additional software tools by the High-Energy Astrophysics Science Archive Research Center (HEASARC), at NASA Goddard Space Flight Center (GSFC). The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l’Énergie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d’Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515. The corresponding authors of this manuscript, namely Axel Arbet-Engels, Lea Heckmann and David Paneque, acknowledge support from the Deutsche Forschungs gemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. The IAA-CSIC group acknowledges financial support from the grant CEX2021-001131-S funded by MCIN/AEI/10.13039/501100011033 to the Instituto de Astrofísica de Andalucía-CSIC and through grant PID2019-107847RB-C44. The POLAMI observations were carried out at the IRAM 30 m Telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). Some of the data are based on observations collected at the Observatorio de Sierra Nevada, owned and operated by the Instituto de Astrofísica de Andalucía (IAA-CSIC). Further data are based on observations collected at the Centro Astronómico Hispano en Andalucía (CAHA), operated jointly by Junta de Andalucía and Consejo Superior de Investigaciones Científicas (IAA-CSIC). Some of the data reported here are based on observations made with the Nordic Optical Telescope, owned in collaboration with the University of Turku and Aarhus University, and operated jointly by Aarhus University, the University of Turku, and the University of Oslo, representing Denmark, Finland, and Norway, the University of Iceland and Stockholm University at the Observatorio del Roque de los Muchachos, La Palma, Spain, of the Instituto de Astrofisica de Canarias. E. L. was supported by Academy of Finland projects 317636 and 320045. The data presented here were obtained [in part] with ALFOSC, which is provided by the Instituto de Astrofisica de Andalucia (IAA) under a joint agreement with the University of Copenhagen and NOT. We acknowledge funding to support our NOT observations from the Finnish Centre for Astronomy with ESO (FINCA), University of Turku, Finland (Academy of Finland grant nr 306531). The research at Boston University was supported in part by National Science Foundation grant AST-2108622, NASA Fermi Guest Investigator grants 80NSSC23K1507 and 80NSSC22K1571, and NASA Swift Guest Investigator grant 80NSSC22K0537. T. H. was supported by the Academy of Finland projects 317383, 320085, 322535, and 345899. This research has made use of data from the RoboPol program, a collaboration between Caltech, the University of Crete, IA-FORTH, IUCAA, the MPIfR, and the Nicolaus Copernicus University, which was conducted at Skinakas Observatory in Crete, Greece. D.B., S.K., R.S., N.M., acknowledge support from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program under grant agreement No. 771282. C.C. acknowledges support from the European Research Council (ERC) under the HORIZON ERC Grants 2021 program under grant agreement No. 101040021. We acknowledge the use of public data from the Swift data archive. Based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. This research has made use of data from the OVRO 40-m monitoring program (Richards et al. 2011), supported by private funding from the California Institute of Technology and the Max Planck Institute for Radio Astronomy, and by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G and NSF grants AST-0808050 and AST- 1109911. This publication makes use of data obtained at Metsähovi Radio Observatory, operated by Aalto University in Finland. AJCT EFG HYD acknowledge support from Spanish MICINN project PID2020-118491GB-I00. This work was supported by JST, the establishment of university fellowships toward the creation of science and technology innovation, Grant Number JPMJFS2129. This work was supported by Japan Society for the Promotion of Science (JSPS) KAKENHI Grant Numbers JP21H01137. This work was also partially supported by the Optical and Near-Infrared Astronomy Inter-University Cooperation Program from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan. We are grateful to the observation and operating members of the Kanata Telescope. The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica. Maunakea, the location of the SMA, is a culturally important site for the indigenous Hawaiian people; we are privileged to study the cosmos from its summit. This work was supported by NSF grant AST-2109127. The Joan Oró Telescope (TJO) of the Montsec Observatory (OdM) is owned by the Catalan Government and operated by the Institute for Space Studies of Catalonia (IEEC). H.K. acknowledges NASA support through the grants NNX16AC42G, 80NSSC20K0329, 80NSSC20K0540, NAS8- 03060, 80NSSC21K1817, 80NSSC22K1291, and 80NSSC22K1883. Author contribution: A. Arbet Engels: project management, coordination of MWL data analysis, NuSTAR analysis, correlation analysis, theoretical modeling and interpretation, paper drafting; L. Linhoff: MAGIC analysis cross-check; I. Liodakis: organization of IXPE and MWL observations, coordination of IXPE and MWL data; L. Heckmann: project management, coordination of MWL data analysis, MAGIC and Fermi data analysis, variability analysis, theoretical interpretation, paper drafting; D. Paneque: organization of the MWL observations, theoretical interpretation, paper drafting; The rest of the authors have contributed in one or several of the following ways: design, construction, maintenance and operation of the instrument(s) used to acquire the data; preparation and/or evaluation of the observation proposals; data acquisition, processing, calibration and/or reduction; production of analysis tools and/or related Monte Carlo simulations; overall discussions about the contents of the draft, as well as related refinements in the descriptions.

References

Appendix A: Additional information: Long-term multiwavelength evolution

Table A.1 summarizes the values depicted in Fig. 5.

All Tables

All Figures

Fig. 1. MWL LC for Mrk 501 between March 1, 2022, and July 19, 2022 (MJD 59639 to MJD 59779). The gray areas mark the IXPE observations from March 8-10, 2022 (MJD 59646 to MJD 59648), March 27-29, 2022 (MJD 59665 to MJD 59667), and July 9-12, 2022 (MJD 59769 to MJD 59772). Top to bottom: (i) MAGIC fluxes in daily bins, with the arrows displaying the upper limits (95% confidence level) for the non-significant bins (< 2 σ); (ii) Fermi-LAT fluxes in 7-day bins; (iii) X-ray fluxes in daily bins, including Swift-XRT, NuSTAR, and IXPE; (iv) hardness ratio between the high- and low-energy fluxes of Swift-XRT; (v) Swift-UVOT; (vi) optical R-band data from RoboPol, Calar Alto, T090, MAPCAT, T150, and AZT-8 in daily bins; (vii) radio data including OVRO, Metsähovi, IRAM, and SMA; and polarization degree (viii) and polarization angle (ix) observations in the X-rays from IXPE, the optical R band from NOT, RoboPol, Calar Alto, MAPCAT, T150, and AZT-8; and the radio band from IRAM and SMA. The IXPE and AZT-8 data are taken from Liodakis et al. (2022) and Lisalda (in prep.). In the polarization angle panel, we plot the average direction angle determined at 43 GHz by Weaver et al. (2022) with a horizontal gray band. See Sect. 4 for further details.

In the text

	Fig. 2. Broadband SED for the three IXPE windows as described in Sect. 3.3. For IXPE-2, the X-ray spectrum corresponding to the highest flux state in the time window is shown with open markers. The typical state (Abdo et al. 2011) and low state (Abe et al. 2023) of Mrk 501 are shown in gray for comparison.
In the text

Fig. 3. MWL polarization and flux properties of Mrk 501 during the three IXPE observations. Top: MWL polarization degree as a function of frequency. For the radio and optical data, weighted averages over the IXPE time windows are used. Middle: ratio of the frequency-dependent polarization degree to the corresponding X-ray polarization degree. Bottom: flux values associated with the polarization values presented above.

In the text

	Fig. 4. Swift-XRT hardness ratio (flux ratio of the 2–10 keV range over the 0.3–2 keV band) compared to flux at 2–10 keV for the LCs shown in Fig. 1. The gray data points are from individual measurements, and the black points are the hardness ratio binned in flux.
In the text

Fig. 5. Fractional variability for the LCs displayed in Fig. 1. Nightly bins are used for MAGIC and 7-day bins for Fermi-LAT; for all other wavebands, the single observations without further binning are used. For the radio data, frequencies are stated to distinguish between the different datasets (from left to right: the orange circles depict OVRO at 15 GHz and Metsähovi at 37 GHz, the lighter orange square the IRAM data at 86 GHz, the orange diamond SMA at 226 GHz, and the darker orange square and open turquoise square IRAM at 230 GHz). The open purple markers refer to the NuSTAR (square) and IXPE (round) observations, which only consist of two or three measurements and therefore far fewer data points than considered for the full time epoch. Due to the different sampling, not all instruments are directly comparable with each other.

In the text

	Fig. 6. Correlation between the MAGIC flux above 200 GeV and the Swift-XRT flux, 0.3–2 keV. See the main text for further details.
In the text

	Fig. 7. Correlation between the MAGIC flux above 200 GeV and the Swift-XRT flux, 2–10 keV. See the main text for further details.
In the text

	Fig. 8. Simplified sketch representing the morphology of the two-zone leptonic model discussed in Sect. 5.
In the text

Fig. 9. Modeling of the three IXPE epochs: IXPE-1 (top), IXPE-2 (middle), and IXPE-3 (bottom). The dashed blue curve is the emission originating from the compact zone, located near the shock front. The dashed-dotted violet curve is the contribution from the extended zone, which spans a larger volume downstream from the shock. The emission resulting from the interaction between the two zones is plotted as a dotted pink curve. The sum of all components is given by the solid blue line. The light gray and dark gray markers show the average state (from Abdo et al. 2011) and low state (from Abe et al. 2023), respectively. The parameters of the model are listed in Table 3.

In the text