Insights into the broad-band emission of the TeV blazar Mrk 501 during the first X-ray polarization measurements

We present the first multi-wavelength study of Mrk 501 including very-high-energy (VHE) gamma-ray observations simultaneous to X-ray polarization measurements from the Imaging X-ray Polarimetry Explorer (IXPE). We use radio-to-VHE data from a multi-wavelength campaign organized between 2022-03-01 and 2022-07-19. The observations were performed by MAGIC, Fermi-LAT, NuSTAR, Swift (XRT and UVOT), and several instruments covering the optical and radio bands. During the IXPE pointings, the VHE state is close to the average behavior with a 0.2-1 TeV flux of 20%-50% the emission of the Crab Nebula. Despite the average VHE activity, an extreme X-ray behavior is measured for the first two IXPE pointings in March 2022 with a synchrotron peak frequency>1 keV. For the third IXPE pointing in July 2022, the synchrotron peak shifts towards lower energies and the optical/X-ray polarization degrees drop. The X-ray polarization is systematically higher than at lower energies, suggesting an energy-stratification of the jet. While during the IXPE epochs the polarization angle in the X-ray, optical and radio bands align well, we find a clear discrepancy in the optical and radio polarization angles in the middle of the campaign. We model the broad-band spectra simultaneous to the IXPE pointings assuming a compact zone dominating in the X-rays and VHE, and an extended zone stretching further downstream the jet dominating the emission at lower energies. NuSTAR data allow us to precisely constrain the synchrotron peak and therefore the underlying electron distribution. The change between the different states observed in the three IXPE pointings can be explained by a change of magnetization and/or emission region size, which directly connects the shift of the synchrotron peak to lower energies with the drop in polarization degree.


Introduction
Blazars are among the brightest objects in the γ-ray sky and are known to emit radiation over a broad range of wavebands from 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 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).Up to recently, polarization measurements of blazars have only been possible in the radio to optical bands.However, these bands are known for originating from more extended or multiple regions in the jet compared to the high-energy emission or are 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, the detection of X-ray polarization of the archetypal blazar Markarian 501 (Mrk 501; z=0.034,Ulrich et al. 1975) was published by Liodakis et al. (2022).Mrk 501 is one of our closest and brightest blazars and therefore a prime object to explore new observational techniques.As usual for blazars, its spectral energy distribution (SED) depicts 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 usually is attributed to the class of high synchrotron peaked blazars (HSPs) with a synchrotron peak frequency ν s > 10 15 Hz (Abdo et al. 2010).However, it has also shown extreme HSP (EHSP) behavior with ν s ≥ 2.4 × 10 17 Hz (∼ 1 keV) (Costamante et al. 2001;Abdo et al. 2010) during extended periods of time encompassing 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 could be induced by relativistic protons inside the jet sparking 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, during its usual HSP states, or the peak, during EHSP behavior, of the synchrotron bump.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 from 2022-03-08 to 2022-03-10 and 2022-03-27 to 2022-03-29 (MJD 59646 to 59648 and MJD 59665 to 59667).In what follows, those two periods will be referred to as IXPE-1 and IXPE-2.They revealed a polarization degree in the 2-8 keV band of ∼ 10% that is a factor of two higher than in the optical regime (Liodakis et al. 2022).Additionally, the polarization angle was found to be aligned with the radio jet.All those properties point towards shocks being the prime acceleration mechanism, as well as an energy stratified jet where the most energetic particles emit from a region which 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 from 2022-07-09 to 2022-07-12 (MJD 59769 to 59772) revealing a drop in polarization degree to ∼ 7%, but with a stable polarization angle compared to the March measurement (Lisalda 2023).This observing epoch is dubbed as IXPE-3 throughout the rest of this paper.
In this work, we present the multi-wavelength (MWL) picture around these three X-ray polarization measurements, including, for the first time, also data up to the VHE regime taken by the Florian Goebel Major Atmospheric Gamma Imaging Cherenkov (MAGIC) telescopes.To characterize the highenergy emission in detail, the campaign involved instruments such as the Neil Gehrels Swift Observatory (Swift), and the Nuclear Spectroscopic Telescope Array (NuSTAR) in the Xrays, and the Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope (Fermi) in the high-energy γ-rays.Both the IXPE pointings as well as the whole time epoch from 2022-03-01 to 2022-07-19 (MJD 59639 to MJD 59779) were accompanied by extensive multi-wavelength coverage from radio to VHE, which is summarized in the following sections.
This paper is structured as follows: In Section 2 we describe all 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 Section 4. Section 5 presents a theoretical modeling of the three IXPE time epochs using a two-zone leptonic model.Finally, Section 6 summarizes and discusses the MWL results together with the obtained theoretical models.

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. 2015a) and are stretching from 2022-03-01 to 2022-07-19 (MJD 59639 to MJD 59779) with the instruments and analyses summarized in this section.

MAGIC
MAGIC is a stereoscopic system of two imaging air cherenkov telescopes (IACTs) 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 hours from 2022-03-01 to 2022-07-19 (MJD 59639 to MJD 59779), which results in 26.5 hours after data selection cuts based on atmospheric transmission.The data are analyzed using the MARS (MAGIC Analysis and Reconstruction Software) package (Zanin et al. 2013;Aleksić et al. 2016).
The VHE flux light curve (LC) is computed for the full energy range above 0.2 TeV, as well as for a low-energy range of Article number, page 3 of 19 0.2-1 TeV and a high-energy range above 1 TeV as shown in the top panel of Figure 1 using nightly bins.For the bins showing a significance of less than 2σ, upper limits (ULs) are computed according to the Rolke method (Rolke et al. 2005) with 95% confidence.The spectral parameters around the three IXPE pointings are 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 is applied to compute the spectral data points (Albert et al. 2007), while for the IXPE-3 spectrum again the forward folding method is 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 (EBL) absorption is corrected for using the model of Domínguez et al. (2011).

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 hours (Atwood et al. 2009;Ackermann et al. 2012).
To analyze the LAT data, we use the unbinned-likelihood tools from the FERMITOOLS software package1 (v2.0.8).We employ P8R3_SOURCE_V3_v1 as the instrument response function, gll_iem_v07 & iso_P8R3_SOURCE_V3_v1 as the diffuse background2 and use 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) is 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 deg for the 68% containment for photons above 0.3 GeV compared to 5 deg for photons above 0.1 GeV), which leads to an improvement in the signal-to-noise ratio since the analysis is less affected by diffuse backgrounds or potential (transient) non-accounted neighbouring 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 use 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 is then fit to the data set in the time range 2022-03-01 to 2022-07-19 (MJD 59639 to MJD 59779).Subsequently, very weak components with a test statistic (TS; Mattox et al. 1996) below 3 are removed from the model as well as those with an expected source count below 1.We then divide the data set in 20 bins of 7-day duration, and redo the fit in each bin varying only the normalization of the diffuse backgrounds, and that of sources within < 3 • and a detection with TS > 10.Furthermore, we allow 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 analysis agree within the statistical uncertainties.In addition, spectral analyses are performed for 2 week intervals cen-tered around each of the three IXPE observations.The same ROI model and approach as before is applied.For each spectrum, we computed the likelihood ratio between a power-law and log parabolic power-law model.No significant (>3σ) preference is found for the log parabolic fit, and therefore a power-law is 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 Section 3.3.The obtained parameters are summarized in Table 1.
The raw data are processed using the NuSTAR Data Analysis Software (NuSTARDAS) package v.2.1.1 and CALDB version 20220912.The events are screened in the nupipeline process with the flags tentacle=yes & 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 source-free nearby circular region having the same radius.The spectra are 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 decide 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 are obtained in the full NuSTAR band-pass, 3-79 keV, and averaged over the respective observations.We fit the spectra using XSPEC (Arnaud 1996) 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 Xray analysis performed in this work, a photoelectric absorption component is added to the model assuming an equivalent hydrogen column density fixed to N H = 1.7 × 10 20 cm −2 (HI4PI Collaboration et al. 2016).The fluxes and spectral parameters are computed by fitting simultaneously FPMA & FPMB.The cross-calibration factor between the two focal module planes is for all bins 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 are averaged over the entire exposure time from each observation.

IXPE
The IXPE telescope (Weisskopf 2022) 2023).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.For more details on the observations, data reduction, etc, see Liodakis et al. (2022) and Lisalda (2023).

Swift-XRT
Simultaneous to the MAGIC observations, several pointing of the XRT (X-ray Telescope; Burrows et al. 2005) on board of Swift (Neil Gehrels Swift Gamma-ray Burst Observatory; Gehrels et al. 2004) were organized.Data were taken both in the Windowed Timing (WT) and Photon Counting (PC) readout mode and processed with the XRTDAS software package (v3.7.0), that was developed by the ASI Space Science Data Center3 (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 extract the X-ray spectrum from the summed cleaned event file.In WT readout mode, the spectral analysis is performed after selecting events from a circle centered at the source position with a radius of 20 pixels (∼ 46 ′′ ; ∼90% of the PSF).In the case of data taken in the PC mode, they are affected by pile-up in the inner part of the 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 use the shape of the Swift-XRT PSF to calculate the incident X-ray flux4 .To estimate the background, we use nearby circular regions with radii of 20 and 40 pixels for WT and PC data.The task xrtmkarf is used to produce ancillary response files (ARFs) including corrections for PSF losses and CCD defects using the cumulative exposure map.The X-ray spectra in the 0.3-10 keV range are binned with grppha ensuring a minimum of 20 counts per bin.Afterwards, we model them in XSPEC using power-law and log-parabola models.Energy fluxes are computed in the 0.3-2 keV and 2-10 keV bands.

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) is used to perform aperture photometry for all filters and the calibration is performed using CALDB (20201026).We review images to avoid unstable attitudes, then following Poole et al. (2008), we use 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 are 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).

Optical
Optical data including flux and polarization measurements were obtained from the 2.5m Nordic Optical Telescope (NOT), the 2.2m telescope of the Calar Alto Observatory as part of the Monitoring AGN with Polarimetry at the Calar Alto Telescopes (MAPCAT)5 , the 1.5m (T150) and the 0.9 m (T090) telescopes at the Sierra Nevada Observatory, the 1.3m RoboPol telescope (Skinakas observatory, Greece; King et al. 2014;Ramaprakash et al. 2019), the 1.5m KANATA (Higashi-Hiroshima observatory, Japan) telescopes and network of robotic telescopes BOOTES (Burst Observer and Optical Transient Exploring System;Castro-Tirado et al. 1999).Except for BOOTES, which uses SDSS filters, Johnson-Cousins filters are used.To convert magnitudes between the different filters Lupton et al. ( 2005) is used.
Additionally, Tuorla blazar monitoring data of Mrk 501 for the period between 2022-03-01 and 2022-07-19 (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 are 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 into account the apertures of the telescopes and the seeing conditions.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 Observatory6 , which took both photometric and polarimetric data with the Cousins R-filter from 2022-03-25 to 2022-03-28 (MJD 59663 to MJD 59666) published in 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 are used.

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 Millimetrique (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).alyzed 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) program7 (Agudo et al. 2018b,a) in the 3.5 mm (86.24GHz) 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 PO-LAMI pipeline described in Agudo et al. (2018a) was adopted to reduce and calibrate the data.
The Submillimeter Array (SMA; Ho et al. 2004) was used to obtain polarimetric millimeter radio measurements at ∼1.3 mm (∼ 226 GHz) within the framework of the SMAPOL (SMA Monitoring of AGNs with POLarization) program.SMAPOL follows the polarization evolution of forty γ-ray loud blazars, including Mrk 501, on a bi-weekly cadence as well as other sources in a target-of-opportunity (ToO) 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 MIR software package 8 .Instrumental polarization leakage was calibrated independently for 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.

Detailed VHE & X-ray analysis during IXPE observations
The ).In the following, those three periods will be 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.
Close to the IXPE-1 epoch, four MAGIC observations are available (see Figure 1).Two observations took place strictly simultaneous to the IXPE pointing, while the two others happened one day before and after.The nightly light curve is consistent with a constant flux hypothesis, with a corresponding average flux of 6.07 ± 0.29 × 10 −11 ph cm −2 s −1 (∼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 (Figure 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 Figure 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. 2023a).We thus average 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−11 ph cm −2 s −1 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. 2023a).
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 −10 erg cm −2 s −1 and ≈ 3 × 10 −10 erg cm −2 s −1 , 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 −11 erg cm −2 s −1 in 2009 (Abdo et al. 2011) or a weighted mean of 4.24 ± 0.01 × 10 −11 erg cm −2 s −1 for 12-year flux data set from 2008 to 2020 (Abe et al. 2023a).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 during 2022-03-09, 2022-03 -22, 2022-03-24 and2022-03-27 (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 −11 erg cm −2 s −1 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. 2023a), but similar to the observations reported in 2013 (Furniss et al. 2015) for Mrk 501.The NuSTAR exposures, which are multihour 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%.2023).In the polarization angle panel, we plot with a horizontal grey band the average direction angle determined at 43 GHz by Weaver et al. (2022).See text in Section 4 for further details.
Article number, page 7 of 19 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), without significant flux variations on daily timescales.

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 fit 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, with a normalization energy of E norm = 1 keV).Furthermore, we fix 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 construct 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 is 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 2022-03-27 (MJD 59665) is 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) are used because they are the only data available around this time window.This ensures a smooth characterization between the two X-ray instruments and therefore of 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 Figure 2.For Fermi-LAT a time window of two weeks centered on the IXPE epochs is 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 between the different days.Due to the availability of VHE observations as discussed in Section 3.1 and Section 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 Figure 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. (2023a).As anticipated in Section 3.1, Figure 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 quantify the peak frequencies and CD for each epoch.The same phenomenological model is fitted to the average and low states from Abdo et al. (2011) and Abe et al. (2023a) 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 17 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 17 Hz.While the MAGIC spectrum shows only a change in amplitude (without spectral evolution), the Fermi-LAT spectrum hints towards softer emission, leading to a shift of the high-energy peak towards lower frequencies by an order of magnitude.
For all three spectra, the fit of the phenomenological model confirms a general lower CD (for IXPE-1 and IXPE-2 ∼0.3, for IXPE-3 0.2 using the phenomenological model or 0.3 using the SSC model described in Section 5) than for the average state of Mrk 501, which is determined to be around ∼0.5.

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 & IXPE-2 in optical and X-rays, the ratio is roughly conserved and stable over time (the derived ratios remain consistent within ∼ 1σ).
As described in Liodakis et al. (2022) and Di Gesu et al. ( 2022), the increase of 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 Article number, page 8 of 19 Notes.Stated are the MJD of the observations, the flux values in 10 −10 erg cm −2 s −1 (for MAGIC with E > 200 GeV; for Fermi-LAT from 0.3 − 300 GeV; for XRT from 2 − 10 keV; from 3 − 7 keV for NuSTAR), the spectral index α and the χ 2 /d.o.f. for the spectral fits.For all wavebands, a simple power law is used as a spectral model (with a normalization energy of 300 GeV for MAGIC, 1 GeV for Fermi-LAT and 1 keV for Swift-XRT) because there is no significant preference for a log-parabola function, except for the NuSTAR data.The NuSTAR fits are performed using a log-parabola model with a curvature parameter fixed to 0.3 (and normalization energy at 1 keV).[Hz] IXPE-1 pheno 5.4 ± 0.2 × 10 17 2.0 ± 0.4 × 10 25 0.30 ± 0.07 IXPE-2 pheno 7.9 ± 0.6 × 10 17 3.0 ± 0.5 × 10 25 0.28 ± 0.06 IXPE-3 pheno 1.0 ± 0.1 × 10 17 4.5 ± 5.9 × 10 24 0.20 ± 0.27 Typical pheno 2.9 ± 0.1 × 10 17 4.6 ± 0.8 × 10 25 0.49 ± 0.10 Low pheno 1.3 ± 0.1 × 10 16 1.8 ± 0.4 × 10 24 0.26 ± 0.10 photons) are confined nearby 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) towards larger regions in which the degree of magnetic field ordering drops significantly.Such energy-stratified scenario qualitatively explains the broad-band polarization degree evolution (Tavecchio et al. 2020), and will be considered for the MWL modeling of the SEDs presented in Section 5.
The lower panel of Figure 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 anti-correlated with the flux, the X-ray hints towards an opposite trend.

Long-term multi-wavelength evolution
This section discusses the MWL evolution throughout the entire period between 2022-03-01 to 2022-07-19 (MJD 59639 to MJD 59779) in order to provide a broader overview of the source evolution between the several IXPE epochs.
From Figure 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, which however might be rooted in the lower instrumental sensitivity of MAGIC to unveil slight spectral changes compared to e.g.X-ray instruments.
The Swift-XRT bands reveal clear hardening with increasing flux.In Figure 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 grey 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. (2023a) 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).
For the UV band, we see a change in flux happening around April slightly increasing the average UV flux level, which can be explained by the shift of the synchrotron peak to lower frequencies.The latter explanation is also suggested by comparing the three IXPE SEDs in Figure 2, which unveil a mild increase in the optical/UV during the IXPE-3 period that displays a shift of the synchrotron component towards lower energies.We quantify the variability strength throughout the spectrum using the fractional variability (F var ) defined following Eq. 10 in Vaughan et al. ( 2003) with the variance of the signal S 2 , the mean square error of the measurement uncertainties < σ 2 err > and the arithmetic mean of the measured flux < F γ > 2 .Its uncertainty is computed following Poutanen et al. (2008): with (3) Considering the whole time epoch from 2022-03-01 to 2022-07-19 (MJD 59639 to MJD 59779), all wavebands show small variability scales.As shown in Figure 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. 2015a;Furniss et al. 2015;Ahnen et al. 2018;Madsen et al. 2021), 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 the signal to noise ratio S/N > 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 data set (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 Figure 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 between the time epochs is found.
We searched for intra-band correlation and found a hint of positive correlation only between the VHE and X-ray light curves, 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 compute between the MAGIC fluxes above 200 GeV and the 0.3-2 keV & 2-10 keV fluxes from Swift-XRT.The DCF is calculated for a series of time lags from -20 days to +20 days between the light curves using time lag steps of 5 days, which is in agreement with the time scale taken from the auto-correlation behavior of Mrk 501 during low activity in Abe et al. (2023a).The significance of the correlation is estimated based on Monte Carlo simulations in a similar fashion as the one described in MAGIC Collaboration et al. (2021).We simulate 10 4 pairs of uncorrelated light curves that respect the sampling and binning of the real data.The DCF significance band in each time lag is then obtained from the distribution of DCF values of the simulations.Ratio to X-ray pol.deg.curves, we assume 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. (2015b).As shown in Figure 6 and Figure 7 the significance is slightly above 2σ at zero time lag.In summary, we find that during the full campaign Mrk 501 has showed very typical behaviors considering its MWL flux val-  ues.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. (2023a) 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. 2023a).
As for the long-term evolution in the polarization, we note an increase of 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 ∼ 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.

Modeling of the SEDs during the IXPE observing windows
We exploit the broad-band SEDs shown in Figure 2 to perform a modeling of the emission assuming a leptonic scenario.We consider an energy-stratified emission region, as suggested by the optical-to-X-ray polarization properties (see Section 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 the shock front.
The theoretical framework that we use 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 Figure 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.
In order to limit the degrees of freedom, we make the following assumptions.Parameters referred to the blob frame are marked as primed, while unmarked ones relate to the observer frame.
-The Doppler factor δ is fixed to 11 consistent with the one used in Abe et al. (2023a).For simplicity, the same value is chosen for both regions.-In order to estimate the relative size between the two zones, we proceed as follows: we first assume that the "compact region" can be modelled 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 Marscher 2014;Tavecchio 2021).Given that IXPE-1 and IXPE-2 are characterized by P deg, X-ray ≈ 10% in the X-ray band, it implies N ≈ 60 during the latter two epochs.Regarding IXPE-3, P deg, X-ray ≈ 7%, corresponds to N ≈ 110.The size of the "extended zone" dominating the optical emission is then assumed to be a factor 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 deg ≈ 75%/

√
l N.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, we obtain that the radius of the "extended zone" should be a factor 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 16 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" is modelled with a broken power-law function normalized such that the electron energy density is given by U ′ e .The distribution is defined between a minimum and maximum Lorentz factor of γ ′ min and γ ′ max .The break is located at γ ′ break .The slopes below and above the break are given by n 1 and n 2 , respectively.γ ′ min is set to a value such that the "compact zone" remains sub-dominant in the optical/UV.With the choice of δ and R ′ mentioned above, we derive γ ′ min ∼ 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 unshocked plasma, are expected to generate somewhat lower γ ′ min , in the order of a few 10 3 (Zech & Lemoine 2021).Values of γ ′ min ∼ 10 4 may be reached for fully relativistic shocks with γ sh of a few tens.In order to constrain γ ′ max , we equate the acceleration and cooling timescales.The acceleration timescale in shock acceleration is estimated with: where λ(γ ′ ) = (ξγ ′ m e c 2 )/(eB ′ ) is the mean free path of electrons, parameterized as a fraction ξ of the Larmor radius.ξ is an acceleration parameter, which we fix to 10 4 .The latter 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 and e the electron charge, m e the electron mass.u s is the speed of the shock, that we assume to be relativistic, u s ∼ c.The synchrotron cooling timescale is: Article number, page 12 of 19 where σ T is the Thomson cross section.-The electron distribution inside the "extended zone" is modelled 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 is assumed to be p = 2.2.The minimum and maximum Lorentz factor of the distribution, γ ′ min & γ ′ max are constrained by the radio-to-UV data.The distribution is normalized to an energy density of U ′ e .We keep γ ′ max within a factor 2 from 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 start 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 ′ e , n 1 , n 2 ).In a second step, the "extended zone" is 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 Figure 1 and Figure 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 Figure 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 compute 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 well all three broad-band SEDs.A more detailed interpretation of the results is given in Section 6.Both the VHE and X-ray regimes showed the highest flux level during IXPE-2 (∼ 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 Figure 4).The opposite behavior can be seen in the high-energy γ-ray data where the spectral index is softening 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 Mrk 501 activity (power-law index ∼2.5;Abdo et al. 2011).
Noteworthy is not just the change in spectral index, but also the broad-band characteristics of the different SEDs (see Table 2), which revealed an EHSP behavior of Mrk 501 in the IXPE-1 and IXPE-2 states with a synchrotron peak frequency well beyond 1 keV.During IXPE-3 the blazar has gone back to its more typical HSP behavior shifting the synchrotron peak by an order of magnitude towards lower energies.The shift in synchrotron peak frequencies can explain the positive change in flux observed in the UV band of the MWL LC around April (see Figure 1).The accompanying shift of the high-energy peak can explain the different flux behavior of the high-energy γ-rays, which increases steadily between the three states compared with the rise and then fall of the fluxes in the VHE and X-rays.During all three states, we find a low Compton dominance compared to the typical state of Mrk 501 (Abdo et al. 2011) by a factor of 2 (see Table 2).Additionally, 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 towards 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 especially the perceived drop in CD for IXPE-3 as observed in Table 2 should be considered with caution due to the dependence on the peak frequency.
Mrk 501 has previously shown EHSP behavior during large flaring activity in 1997 (Tavecchio et al. 2001) and 2013(Furniss et al. 2015).This EHSP behavior during large flaring ac-  (2018), that used the extensive multi-instrument campaign in 2012.Compared to the data from 2022 that is reported in this manuscript, the above-mentioned studies from past campaigns showed a far higher variability in both VHE and X-rays.While similar average flux levels are found in the X-rays in 2012 compared to this work, 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 the 2022 campaign, together with higher degrees of variability, reaching up to F var ∼ 1 in the VHE.The main difference in behavior during 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 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 Figure 3).In fact, the X-ray polarization is systematically a factor ∼ 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 broad-band SEDs are more probable to explain the differences in polarization degree between the different states as will be discussed further 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 Figure 1, bottom Article number, page 14 of 19 panel).Nonetheless, in the middle of the campaign, both the radio (SMA; 226 GHz) and optical (R-band) polarization angles depict some variability and deviate from the jet angle in an incoherent manner.Indeed, around 2022-05-31 (MJD 59730), the offset between the radio and optical bands reaches ∼ 60 • .This behaviour strongly points towards 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 short scale occurs (Lyutikov et al. 2005;Myserlis et al. 2018;Tavecchio 2021).Furthermore, it points towards 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/turbulent contributions as also proposed in Abe et al. (2023a).
Over the full time period, the blazar showed a rather stable behavior with typical (Abdo et al. 2011) flux levels in the VHE (∼20-50% C.U.).In contrast, the XRT flux values are constantly above 10 −10 erg cm −2 s −1 .Based on the long-term light curve presented in Abe et al. (2023a), this indicates an elevated state with respect to the average X-ray activity (Abdo et al. 2011).It confirms that the atypically low Compton dominance 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 ones in the X-rays for only the IXPE time windows (see Figure 5).Even though different behaviors are reported in the X-rays and VHE, we identify evidence for a correlation between the two bands without a time lag (see Figures 6 and 7).The correlations are found with significances between 2σ to 3σ, limited by the short time epoch considered, the low flux levels, the very low variability, and the somewhat limited sensitivity in the VHE γ-ray band in comparison with the higher sensitivity to measure flux changes in the X-ray band.
To further investigate the emission states behind the first IXPE pointings, we modelled the broad-band SEDs with a twozone leptonic model (see Section 5).We assume a "compact" region, which is nearby the shock front and populated by freshly accelerated electrons, responsible for the X-ray and VHE γ-ray band.The "compact region" is embedded into an "extended region" which expands downstream 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 motivated by the energy dependency of the polarization degree between the radio and the X-rays (Liodakis et al. 2022), which points towards an energy stratified jet.In addition, previous longterm correlation studies (e.g.Abe et al. 2023a) points towards the same emitting region dominating the X-ray and γ-ray bands.While the "extended" zone is assumed to be constant between the three pointing, the "compact" one varies over time to describe the observed spectral changes.
A good description of the observations from radio to VHE is achieved for each epoch.The X-ray data during IXPE-1 and IXPE-2 constrain well both the rising and falling edges of the synchrotron SED close to the peak energy.This provides in turn strong constraints on the spectral shape of the electron distribution within the region nearby the shock front.Below the break, which either results from cooling or acceleration effects, the derived index of the electron distribution is n 1 = 2.25 for the "compact zone".This is well in agreement with expectation from diffusive acceleration in relativistic shocks (Kirk et al. 2000).Regarding the IXPE-3 epoch, similar n 1 are obtained although in those cases n 1 is less constrained given the shift of the synchrotron peak towards lower energies.
The emission resulting from the interaction between two zones brings a significant contribution 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. (2023a) detected a positive correlation between the Swift-XRT and Fermi-LAT emission, which is thus in good agreement with our model.
Between the different epochs, most of the parameters of the "compact zone" are quite similar.A main difference is the drop of γ ′ br by a factor ≈ 4 during IXPE-3 with respect to IXPE-1 and IXPE-2 in order to explain the shift towards lower energies of the synchrotron component.The obtained γ ′ break for IXPE-1 and IXPE-2 are a factor of 4 higher than the one expected for a self-consistent equilibrium between 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 of B ′ in the "compact zone" for IXPE-3 points towards an evolution of the magnetization.Within a multi-cell scenario, the shift in peak frequencies can be attributed to an increase of synchrotron cooling due to stronger magnetic field strength.At the same time, the drop in polarization implies that the magnetic field becomes less ordered and that the emitting region is composed by a larger number of turbulent plasma cells.It has been shown before that turbulences created by relativistic shock propagation could lead to an amplification in 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 of 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 κ, it can be expected that an increase of the magnetization caused by a more turbulent plasma results in an overall decrease of 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 assume a stable magnetic field of 5 × 10 −2 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" should 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 by ∼twice the number of turbulent cells compared to IXPE-1 and IXPE-2 (see Sect. 5).In summary, based on our modeling, the decrease of the polarization degree throughout the IXPE epochs may be qualitatively explained by a change of magnetization and/or a change of the emitting region size.
Article number, page 15 of 19 In conclusion, with this 2022 data set we could, for the first time, combine broad-band MWL data up to the VHE together with polarization measurements up to the X-rays, which allowed us to better constrain the emission and acceleration behind 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 by polarization measurements at even higher energies covering 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 years, could add even more vital insights.

Fig. 3 :Fig. 4 :
Fig.3: Top: Multi-wavelength polarization degree as a function of frequency for the three IXPE observations.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.

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

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

Fig. 7 :
Fig. 7: Correlation between the MAGIC flux above 200 GeV and the Swift-XRT flux 2-10 keV.See text for further details.
Article number, page 6 of 19 S. Abe et al.: Radio to TeV observations of Mrk 5011 during the IXPE campaigns

Table 1 :
Parameters for the VHE and X-ray observations around the three IXPE pointings.

Table 2 :
Ghisellini et al. (2017)d ν IC , and Compton dominance (CD) for the different SEDs shown in Fig.2extracted from the maxima of the phenomenological description ofGhisellini et al. (2017).
(Abe et al. 2023a)SED for the three IXPE windows as described in Section 3.3.For the 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. 2023a)of Mrk 501 are shown in grey for comparison.
To simulate the light Article number, page 10 of 19 S. Abe et al.: Radio to observations of Mrk 5011 during the IXPE campaigns

Table 3 :
Parameters of the two-components of the leptonic model obtained for the three IXPE epochs.