© The Authors 2025

1. Introduction

Blazars are a subclass of active galactic nuclei (AGNs), each characterized by a powerful relativistic plasma jet (Blandford et al. 2019) whose axis is aligned at a small angle with the observer’s line of sight, leading to strong relativistic aberration of the observed radiation. Their broadband emission is dominated by nonthermal radiation from the jet that goes from the radio to the very high-energy (VHE; > 0.1 TeV) gamma ray band.

One of the key features of their emission is a high degree of variability observed over the full spectrum of timescales from years to hours (see, e.g., Fossati et al. 2008; Chatterjee et al. 2008). During high-flux states, commonly called flares, extreme flux variations up to an order of magnitude on the timescale of minutes have also been reported (Albert et al. 2007; Aharonian et al. 2007). Flares are not only accompanied by strong flux changes, but also by large spectral variations (Pian et al. 1998), hence implying that particle acceleration mechanisms play a central role in the origin of the variability. However, all those phenomena, and in particular the underlying processes that accelerate particles to highly relativistic energies, are still a topic of debate.

Among the acceleration mechanisms commonly considered for blazars there are shock acceleration (Marscher 1978; Crumley et al. 2019) and magnetic reconnection (Sironi & Spitkovsky 2014). For the first, the source of acceleration is the interaction with a collisionless shock wave, where turbulence in the magnetic field of the plasma are responsible for the multiple crossings at the shock front, leading to an efficient energy gain. In the case of magnetic reconnection, it is the magnetic field lines themselves that are the root of the acceleration. Through instabilities in the magnetic field, field lines of opposite polarity can reconnect with each other and transfer magnetic energy to kinetic energy efficiently. Which of these two mechanisms is dominating in blazars is still unknown.

It is also dependent on where we expect the observed emission taking place in the jet due to different dependences on the magnetization of the jet.

An important tool for probing acceleration mechanisms are polarization measurements since they directly relate to the structures of the magnetic fields at the emission sites, which in turn lead to constraints on the possible acceleration mechanisms (Tavecchio 2021). As discussed in Tavecchio et al. (2020), Di Gesu et al. (2022a), in the scenario of shock acceleration, the polarization degree is expected to show a strong chromatic behavior. The polarization increases with energy because the highest-energy particles are located closer to the shock front where they are freshly accelerated, and where the magnetic field is more ordered as it gets compressed. The polarization degree is expected to be slowly variable, and the angle parallel to the shock normal (i.e., parallel to the jet axis in the case of a shock normal aligned with the jet). In the case of a highly turbulent plasma crossing a shock front, both the polarization degree and angle can exhibit strong variability dictated by the stochastic nature of the magnetic field in the plasma cells crossing the shock. The average polarization angle is nevertheless expected to be roughly aligned with the shock normal due to the partial ordering and compression of the field by the shock. In this context, Marscher (2014), Marscher et al. (2017) also developed the Turbulent Extreme Multi-Zone Model for Blazar Variability (TEMZ), in which a turbulent plasma crosses a standing shock (e.g., recollimation shock). This model predicts a chromatic behavior of the polarization driven by the turbulent nature of the plasma, and is consistent with the radio to optical range polarization behaviors seen in several blazars (e.g., MAGIC Collaboration 2018). For magnetic reconnection, the polarization degree and angle can show fast (below the light-crossing time of the reconnection layer) and chaotic variability driven by the complex geometry of the current sheets (Zhang et al. 2022).

Until recently, polarization observations of blazar jets were only available in the optical and radio wavebands. However, since the end of 2021, the Imaging X-ray Polarimetry Explorer (IXPE; Weisskopf 2022) has been providing us with linear polarization observations in the 2–8 keV band. For high synchrotron peaked blazars (HSPs), which have a synchrotron spectral energy distribution (SED) peaking above 10¹⁵ Hz (Abdo et al. 2010), the IXPE energy regime is measuring the synchrotron emission from the most energetic, freshly accelerated electrons. Therefore, X-ray polarization measurements directly probe the conditions close to the acceleration regions for HSPs, while optical and radio observation are most probably connected to regions farther downstream the jet due to the longer cooling timescale (Tavecchio 2021). Combining radio-to-X-ray polarization observations provides a unique opportunity to disentangle the acceleration mechanisms in blazars, and thus better understand the origin of flaring events.

During its first years of operations, the HSPs observed by IXPE had been probed in nonflaring activity. In 2022, IXPE announced the first detection of X-ray polarization from a blazar (Liodakis et al. 2022), the archetypal HSP Markarian 501 (hereafter Mrk 501), revealing significantly higher polarization degrees than in the optical or radio regimes, but similar polarization angles. The same was found for another archetypal HSP, Markarian 421 (hereafter Mrk 421; Di Gesu et al. 2022b, 2023). This indicated (at least in nonflaring activity) that shock acceleration in an energy stratified jet is preferred. This scenario is further supported by the results derived with the data from the multi-instrument campaigns organized for these two sources during the year 2022 (Abe et al. 2024; MAGIC Collaboration 2024).

In December 2023 another IXPE observation of Mrk 421 was performed, which lasted from December 6 (MJD 60284) to December 21 (MJD 60300), around which we organized a dense monitoring campaign from radio to VHE. This time, the source was found in a flaring state reaching VHE fluxes of more than twice that of the Crab Nebula, and in a bright X-ray state compatible with prominent archival flares (Aleksić et al. 2015a). In this work we present the first broadband study from radio to VHE of a gamma-ray flare of a blazar in combination with simultaneous radio-to-X-ray polarization measurements. For this study we coordinated and combined extensive observations by the Florian Goebel Major Atmospheric Gamma Imaging Cherenkov (MAGIC) telescopes (Aleksić et al. 2016), the Large Area Telescope (LAT) on board Fermi Gamma-ray Space Telescope (Fermi-LAT; Atwood et al. 2009; Ackermann et al. 2012), the Neil Gehrels Swift Observatory (Swift; Gehrels et al. 2004), the X-ray Multi-Mirror Mission (XMM-Newton), and various ground-based telescopes to cover the optical and radio frequencies.

The paper is structured as follows. The observations and data processing are described in Sect. 2. The multi-wavelength behavior and polarization characteristics during the campaign are discussed in Sect. 3; in Sect. 4 we characterize the intra-band correlations. In Sect. 5 we describe the modeling of the broadband emission throughout the IXPE window using constraints from the observed multi-wavelength polarization, and the discussion and summary are presented in Sect. 6.

2. Observations and data processing

Most of the observations were carried out by the same instruments and using the same data processing as for the multi-wavelength campaigns recently published in Abe et al. (2024) and MAGIC Collaboration (2024). We refer to the latter works for more details, but we provide below a summary for completeness.

Observations in the VHE band were performed by the MAGIC telescopes at the Observatorio del Roque de los Muchachos (ORM) on the Canary Island of La Palma, Spain. The data were analyzed with the standard procedures and the MAGIC Analysis and Reconstruction Software (MARS; Zanin et al. 2013; Aleksić et al. 2016) software package as in Abe et al. (2024). We extracted daily light curves for the full campaign in the 0.2–1 TeV and > 1 TeV bands. Moreover, a deep MAGIC exposure of 4.3 hours (15.6 ks) was organized during the multi-hour observations from XMM-Newton that took place on December 13, 2023 (MJD 60291), which was used to derive an intranight light curve above 0.4 TeV in 25-minutes intervals. The latter minimum energy is higher than the threshold adopted for the daily light curve because, on that night, the MAGIC observations were performed at zenith angles from ≈10° to ≈60°. At zenith angles close to ≈60°, the threshold of MAGIC is ≈350 GeV. Hence, using an increased threshold of 0.4 TeV allows us to include all the 25-min bins, even those with data taken at the highest zenith angles up to 60°. Nightly SEDs and spectral parameters were obtained with a forward folding method. Since for the long exposure night, a log-parabola model was preferred by more than 3σ, it is employed for all nights to ease comparability.

For the high-energy (HE) gamma rays, data from Fermi-LAT were obtained and analyzed following the same procedure as in MAGIC Collaboration (2024). We produced a light curve with a binning of 3 days and SEDs centered around each MAGIC observation also integrated over 3 days. A simple power-law model was used to produce the SEDs. We note that, when using the 15 days considered in this study (December 6–21), a log parabola is not preferred (with respect to a power law) by more than 3σ.

A dense X-ray monitoring campaign with the Swift-XRT telescope (Burrows et al. 2005) was organized to support the MAGIC observations. A special effort was put to schedule the observations simultaneously with the MAGIC observations. The data were reduced as in Abe et al. (2024) using an updated version of the XRTDAS software package (v.3.7.0) developed by the ASI Space Science Data Center¹ (SSDC), released by the NASA High Energy Astrophysics Archive Research Center (HEASARC) in the HEASoft package (v.6.32.1). We extracted fluxes in the 0.3–2 keV and 2–10 keV bands by fitting a log-parabola model assuming a hydrogen column density fixed to N_H = 1.34 × 10²⁰ cm⁻² (HI4PI Collaboration 2016). The same N_H will be used throughout this work. In the vast majority of the cases (> 97% of the fits), the log-parabola model is preferred over the power-law model with a significance above 3σ.

On December 13, 2023 (MJD 60291), we performed a multi-hour observation with the XMM-Newton, which carries on board several coaligned X-ray instruments. One of them is the European Photon Imaging Camera (EPIC), consisting of the Metal Oxide Semiconductor cameras (EPIC-MOS; Turner et al. 2001) and the pn junction camera (EPIC-pn; Strüder et al. 2001) that both operate in the 0.2–10 keV band. In this work, we focus on the data from the EPIC-pn camera considering its better sensitivity that allows us to resolve spectral changes on short timescale, which is of importance for this work (see Sect. 3.2). Within instrumental systematics, both the EPIC-pn and EPIC-MOS show consistent results. The data were taken in TIMING mode with the THICK filter, for an overall exposure of 16.9 ks (4.6 h) after good time interval screening. We extracted the source and background spectra following the same analysis procedures as in de la Calle Pérez et al. (2021), Abe et al. (2024). Due to a count rate close to the threshold above which pile-up may occur, we removed the central column around the source to suppress any potential pile-up artifacts. The fluxes in the 0.3–2 keV and 2–10 keV bands were computed by fitting a log-parabola model, which is significantly preferred over a simple power law.

The IXPE telescope (Weisskopf 2022) is the first instrument capable of resolving the X-ray polarization degree and angle in blazars. In December 2023 four observations took place, from December 6 (MJD 60284) to December 21 (MJD 60300), for a total exposure of 514 ks spread over four observations. The data processing was performed using the ixpeobssim software, version 30.6.3 (Baldini et al. 2022). As in Maksym et al. (2024), the I, Q, and U spectra were determined in the 2–8 keV band and the pcube algorithm was used to obtain polarization angle and degree. We refer to Maksym et al. (2024) for more details on the source and background region selection. The data were binned into ≈12 h intervals, for which a significant detection of the X-ray polarization is obtained in all intervals, revealing a strong polarization variability. To investigate variability on shorter timescale, the data were also binned into 6 h intervals, and farther down to 3 h intervals contemporaneous to the long XMM-Newton observations.

In the UV band, we analyzed the Swift-UVOT images in the W1, M2, and W2 filters, for the observations in the time interval of interest. We applied a reduction and data analysis procedure similar to those of Abe et al. (2024). The identical HEAsoft software version was used for the aperture photometry task, as well as the same CALDB release to apply the standard calibrations (Breeveld et al. 2011). We included calibrations systematic errors to the magnitude statistical uncertainties (Poole et al. 2008, 2005; Poole & Breeveld 2005). However, various observations in this time interval are affected by attitude instabilities, due to increased noise to one of the spacecraft gyroscopes. We checked carefully image photometry, executing photometry in few cases to single exposure slices to recover some observations, and we ended discarding about one-sixth of the observations.

Optical photometry and polarimetry observations in the R-band were performed by the 2.2 m telescope of the Calar Alto Observatory as part of the Monitoring AGN with Polarimetry at the Calar Alto Telescopes (MAPCAT²; Agudo et al. 2012), the 1.5 m (T150), and the 0.9 m (T090) telescopes at the Sierra Nevada Observatory. We also make use of observations from the Nordic Optical Telescope at the ORM (NOT; Nilsson et al. 2018), the KANATA telescope (Higashi-Hiroshima observatory, Japan Kawabata et al. 1999; Akitaya et al. 2014), the Liverpool Telescope (Jermak et al. 2016; Shrestha et al. 2020), and the Boston University’s Perkins telescope (Perkins Telescope observatory, Flagstaff, AZ). Finally, we obtained R-band photometry data from the Tuorla blazar monitoring program using the 80 cm Joan Oró Telescope (TJO) at Montsec Observatory, Spain. The details of the data analysis an reduction for the different telescopes can be found in Kouch et al. (2024), Abe et al. (2024), MAGIC Collaboration (2024), Liodakis et al. (2022), Nilsson et al. (2018), Otero-Santos et al. (2024), Escudero Pedrosa et al. (2024). Additional R-band data were acquired with a Moravian G4-16000 CCD camera equipped with Bessel filters from Chroma Technology at the 0.5m PlaneWave CDK astrograph of the Hans-Haffner-Sternwarte in Hettstadt, Germany, as part of the long-term AGN observation program of the science laboratory for students at the Friedrich-Koenig-Gymnasium (FKG), the University of Würzburg, and the TU Dortmund University. To extract the fluxes, photometric settings, and comparison stars were taken from the Glast AGILE Support Program (GASP) list of the Whole Earth Blazar Telescope³. The polarization data were corrected for the contribution of the host galaxy with the same method as Hovatta et al. (2016) and using the host fluxes reported in Nilsson et al. (2007). In order to build broadband SEDs, the flux densities were also corrected for a galactic extinction of 0.033 mag according to the NASA/IPAC Extragalactic Database (NED)⁴.

In the radio band, the flux density and polarization were measured in the mm band (225.5 GHz) by the SMA (Ho et al. 2004) telescope within the framework of the SMA Monitoring of AGNs with POLarization (SMAPOL) program (Myserlis et al., in preparation) and the data reduction was performed as in MAGIC Collaboration (2024). We also obtained observations in the cm range (4.85 GHz, 10.45 GHz and 14.25 GHz) thanks to the Monitoring the Stokes Q, U, I, and V Emission of AGN jets in Radio (QUIVER) program (Myserlis et al. 2018; Kraus et al. 2003) using the Effelsberg 100 m telescope. We refer to Myserlis et al. (2018) for more details on the analysis methods.

3. Characterization of the broadband emission and polarization behavior

Fig. 1 depicts the multi-wavelength light curves from the radio up to the VHE gamma-ray range, that span from November 7, 2023 (MJD 60255), and January 16, 2024 (MJD 60325). These observations were performed around the long IXPE observation taking place from December 6, 2023 (MJD 60284), to December 22, 2023 (MJD 60300).

Fig. 1. Multi-wavelength light curve between November 7, 2023 (MJD 60255), and January 16, 2024 (MJD 60325). The red line marks the XMM-Newton observation taking place on December 13, 2023 (MJD 60291), and the vertical gray bands highlight the IXPE observing windows. From top to bottom: MAGIC fluxes in the > 1 TeV (first panel) and 0.2–1 TeV (second panel) bands with nightly bins (the typical nightly exposure is ≈40 min). The gray dashed line depicts the average flux of Mrk 421 from Acciari et al. (2014); Fermi-LAT fluxes in three-day bins; X-ray fluxes binned per observation including Swift-XRT and XMM-Newton; Swift-UVOT and optical R-band data using the telescopes listed in Sect. 2; radio data from the SMAPOL and QUIVER programs; polarization degree and polarization angle observations in the X-rays from IXPE (in black and empty gray markers the data are binned over 12 h and 6 h bins, respectively), the optical R-band and the radio band. In the polarization angle panel, we plot with a horizontal blue band the average direction angle determined at 43 GHz by Weaver et al. (2022) (see text in Section 3 for further details).

The MAGIC light curves in the top two panels show that Mrk 421 was in a high VHE state during most of the observation epoch. Notably, two flaring periods exceeding 1 Crab Nebula units⁵ (C.U.; see red dashed line in Fig. 1) can be seen, one in November 2023 and another one in December 2023. For the flare in December 2023, which was simultaneous with the IXPE observations, the 0.2–1 TeV flux reached more than 2 C.U. toward the end of the flare. We note that the average flux of Mrk 421 at those energies, plotted as a gray dashed line in Fig. 1, is around 0.45 C.U. (Abdo et al. 2011; Acciari et al. 2014), which is more than a factor of 4 smaller. The highest VHE flux levels in December 2023 are four to eight times brighter than during the first IXPE observations of Mrk 421 in 2022 where the simultaneous 0.2–1 TeV flux ranged from 0.25 to 0.5 C.U. (Abe et al. 2024).

The Swift-XRT light curves reveal a bright X-ray state, in particular during the IXPE window. The 0.3–2 keV fluxes are regularly above 10⁻⁹ erg cm⁻² s⁻¹, which is higher by at least a factor of 2 when considering previous campaigns with close-to-average X-ray states (see, e.g., MAGIC Collaboration 2021). The 0.3–2 keV fluxes are comparable to the bright flare of March 2010 (Aleksić et al. 2015a). Differently from the 0.3–2 keV band, the 2–10 keV fluxes are closer to average values (≈0.4 × 10⁻⁹ erg cm⁻² s⁻¹), if one refers to the 6-month campaign of 2017 discussed in MAGIC Collaboration (2021). This difference between the 0.3–2 keV and 2–10 keV fluxes relative to the average states from archival data points toward a softer X-ray spectrum than usually observed. This is confirmed in the spectral study presented in Sect. 3.2. In the other wavebands, radio, optical, UV, and HE gamma rays, the fluxes show an absence of strong variability. The emission is comparable to the typical state of Mrk 421 (Abdo et al. 2011), as well as the one noticed in 2022 during the first IXPE observations of Mrk 421.

The fractional variability (F_var) is reported in Fig. 2 for all the available wavebands with sufficient number of observations to be representative of the campaign. For instance, the UV bands are excluded because of the very low number of measurements, as well as the highly non-uniform distribution of the observations, for any single UV band from Swift-UVOT. F_var quantifies the variability strength following Eq. (10) in Vaughan et al. (2003), and the corresponding (statistical) uncertainty is estimated following Poutanen et al. (2008). For the full time period, the strongest variability is seen in the VHE gamma rays and the X-rays, with an increased variability (F_var > 0.5) in the higher sub-energy range in these two bands (> 1 TeV and 2–10 keV). During the IXPE window (shown with green markers), the VHE gamma rays and X-ray F_var is ≈0.2–0.3, and it is very similar between the two sub-energy bands in the X-ray and VHE regimes. This is different from the strong chromatic trend observed over the full campaign. The lower variability and essentially achromatic behavior of the fractional variability during the IXPE observations may be indicative of a different state. However, we cannot exclude that this difference is due to the shorter time interval that is being considered (15 days vs 70 days), hence reducing the probability to catch flux variability in the source.

Fig. 2. Fractional variability, Fvar, for the light curves displayed in Fig. 1. Nightly bins are used for MAGIC, three-day bins for Fermi-LAT, and for all other wavebands the single observations are without further binning. The filled orange to violet markers depict Fvar when data over the full campaign are considered, and the open turquoise markers only include measurements during the IXPE time window. For Fermi-LAT, Fvar cannot be computed during the IXPE window due to a variability level below the statistical uncertainties, leading to a negative value in the square root of the numerator of Fvar (see Vaughan et al. 2003).

3.1. Multi-wavelength polarization evolution

In addition to the multi-wavelength fluxes, polarization measurements in the radio, optical, and X-ray bands are shown in the two bottom panels of Fig. 1. In the different radio bands, the polarization degree varies by about $\pm 0.5\times {\overline{P}}_{\mathrm{radio},\mathrm{deg}.}$, with a mean value of ${\overline{P}}_{\mathrm{radio},\mathrm{deg}.}\approx 2\%$. A similar degree of variability is found in the R-band but the mean value is higher, ${\overline{P}}_{R-\mathrm{band},\mathrm{deg}.}\approx 4.5\%$, both for the full campaign as well as during the IXPE observation window. The variability in the X-ray polarization degree is much stronger. The IXPE polarization degree binned over 12 h intervals (black markers) show variations by more than a factor of 4 between the minimum (6 ± 2%) and maximum (26 ± 2%). Over the 6 h intervals (gray empty circles) the polarization degree shows a similar level of variability, implying that the polarization evolved at least down to ∼6 h timescale at some instances during the campaign. The average X-ray polarization is ≈13%, significantly larger than the optical and radio bands. The highest X-ray polarization degree coincides with simultaneous long exposures from XMM-Newton and MAGIC. The multi-wavelength evolution during this long exposure is discussed in more details in the next section.

In Fig. 3, the temporal evolution of the X-ray-to-optical (black markers) and radio-to-optical (blue markers) polarization degree ratios are shown to illustrate the chromatic behavior. In addition to the polarization being on average higher in the X-rays, a strong variability in the ratios can be seen. On MJD 60291 (December 13, 2023, coinciding with the long exposure from XMM-Newton and MAGIC), the X-ray polarization degree is above the R-band by a factor of ≈8–9. For other days (e.g., MJD 60295, December 17, 2023), the ratio is close to 1. The increase with energy of the polarization degree is in line with previous observations (Di Gesu et al. 2022b; Kim et al. 2024), but it is the first time that such variations in their respective ratios are reported (see Abe et al. 2024). In the HSP Mrk 501 (which shares very similar spectral properties with Mrk 421), the ratio between the polarization degrees stayed constant over several IXPE observations (MAGIC Collaboration 2024).

Fig. 3. Ratio of the X-ray (2–8 keV) and radio (4.85 GHz) polarization degree to that measured in the R-band during the IXPE observing window (PX − ray, deg/PR − band, deg and Pradio, deg/PR − band, deg, respectively). The vertical dashed line marks the date of the simultaneous XMM-Newton/MAGIC long exposure.

The polarization angle from the radio to optical bands shows variations by at most 40° over the full observation period (bottom panel in Fig. 1). The optical polarization angle is consistent with the results from the 2022 campaign (Di Gesu et al. 2022b; Abe et al. 2024) and aligns well with the radio jet direction reported in Weaver et al. (2022) (−14.4° ±14.2°, or 165.6° ±14.2° taking into account the 180° ambiguity). However, our radio data show systematic shifts between the measured angles at different radio frequencies. It ranges from ≈90° in the 4.85 GHz band, and gradually increases with frequency, reaching ≈150° at 225.5 GHz, which is much closer to the angle in the optical. This is indicative of a Faraday screen with a rotation measure of ∼ − 200 rad/m². Alternatively, it may indicate a bending of the jet, which is a relatively common feature in jetted AGNs (Bridle et al. 1994). The gradual shift toward lower frequencies suggests that the bending occurs in broader regions located downstream of the one dominating the emission beyond optical frequencies, because at the lowest frequencies the radio flux in compact regions is synchrotron self-absorbed (Konigl 1981). Further studies are needed to properly assess the origin of this chromatic evolution of the radio polarization angle.

Differently from the other bands, the X-ray polarization angle displays strong, but stochastic-like fluctuations of ≈100°, from ≈130° to 230°. This behavior is quite different from the systematic angle rotation of Mrk 421 reported in Di Gesu et al. (2023) which proceeded over 5 consecutive days with a roughly constant velocity. We note however that the average value of the polarization angle during December 2023 is ≈180°, aligning well with the results in the optical band and the direction of the radio jet (Weaver et al. 2022).

3.2. Spectral evolution

To investigate the spectral behavior of Mrk 421 during the period studied here and compare it with archival observations, we computed the hardness ratio in the X-ray and VHE regimes using the high- and low-energy flux in each band. Fig. 4 depicts the hardness ratios in the different bands compared to the low- and high-energy flux levels for this campaign. Points with black-filled markers highlight measurement during the IXPE window. We quantified the correlation between the hardness ratio and the fluxes using the Pearson’s R coefficient and the corresponding statistical significance in Gaussian units was derived following the prescription of Press et al. (2007). The results are shown in each panel, both for the full 2023 campaign and the IXPE window only. Additionally, for comparison purposes, we plot the observations from the 2022 campaign (Abe et al. 2024), when the first IXPE observations of Mrk 421 were taking place and with the 2017 campaign published in MAGIC Collaboration (2021), which is representative of the typical dynamical flux range of Mrk 421 in the X-rays and at VHE. In comparison to the archival campaigns, the current data lies in a softer regime than usual for the observed high flux level.

Fig. 4. Hardness ratio in the VHE (upper panels) and X-rays (lower panels) bands. The VHE hardness ratio is defined as the ratio of the > 1 TeV to the 0.2–1 TeV fluxes, while in the X-ray it is the ratio of the 2–10 keV to the 0.3–2 keV fluxes. They are plotted vs the 0.2–1 TeV and > 1 TeV bands, while in the X-ray they are plotted vs the 0.3–2 keV and 2–10 keV bands. The colored circular markers correspond to measurements from the campaign under study, and highlight the data simultaneous with the IXPE window with black-filled markers. Archival measurements from the extensive multi-instrument observing campaigns in the years 2015–2016 (Acciari et al. 2021), 2017 (MAGIC Collaboration 2021), and 2022 (Abe et al. 2024) are plotted as gray triangles, gray squares, and black diamonds, respectively.

A strong indication of a harder-when-brighter behavior is observed in the high-energy X-ray and VHE gamma-ray bands over the full 70-day dataset presented in this manuscript, from November 7, 2023, to January 16, 2024 (see the right panels in Fig. 4). The significance of the correlation is 11σ and 4.4σ for the 2–10 keV and > 1 TeV bands, respectively. For the low-energy band fluxes, 0.3–2 keV and 0.2–1 TeV (see Fig. 4 left panels), however, the significance is much lower and hints toward multiple trends as indicated by the two branches seen in the figures. One follows a steeper trend similar to the one shown in the high energies. The second one depicts an almost flat relation extending to high fluxes with low hardness ratios. This flatter branch can be associated with the December flare, indicating that no harder-when-brighter behavior is seen during the IXPE time window.

Also in the archival data the harder-when-brighter trend is more visible in the comparison with the high-energy fluxes than for the low-energy band. However, hardness ratios as low as during the IXPE window have not been observed in the 2017 and 2022 campaigns in combination with the high flux levels.

For completeness, we present in Appendix A and Appendix B all the spectral parameters for each of the Swift-XRT and MAGIC observations during the campaign. For Swift-XRT, we show the results of power-law (dN/dE ∝ (E/E₀)^−Γ) and log-parabola fits (dN/dE ∝ (E/E₀)^{−α − βlog(E/E₀)}) with a pivot energy E₀ fixed to 1 keV. In the case of MAGIC, the spectra are fitted with a log-parabola model having β fixed to 0.5, and a normalization energy set to E₀ = 300 GeV. The motivation to fix β is the following: although we find a significant preference (> 3σ) for a log-parabola model over a simple power law for several of the days, β does not exhibit any significant dependence on the flux level. By fixing β we still obtain a satisfactory description of all the days, and α can be directly used to evaluate the spectral hardness evolution given that the correlation between α and β is removed. The parameters from the MAGIC fits are the intrinsic ones, in other words after correcting for the effect of extragalactic background light (EBL) absorption using the model of Domínguez et al. (2011).

3.3. Long MAGIC/XMM-Newton exposure on MJD 60291

The organization of multi-band long exposures represents a good opportunity to probe flux and spectral variations down to sub-hour scales, being the timescale over which Mrk 421 (and blazars in general) are known to vary. Such investigations are also essential to provide constraints on the emitting region dimension using causality arguments. Figure 5 presents the intra-night multi-wavelength light curves on December 13, 2023 (MJD 60291) during the simultaneous long exposure with XMM-Newton and MAGIC. The total exposure is 16.9 ks for XMM-Newton (in the pn camera) and 15.6 ks for MAGIC. The MAGIC light curve (top panel) is computed over 20-min intervals above 0.4 TeV. The XMM-Newton fluxes (second panel from the top) are computed over ≈15 min bins, in the 0.3–2 keV and 2–10 keV bands. We also show in the third panel from the top the hardness ratio from the XMM-Newton data, which we defined as the ratio between the 2–10 keV and 0.3–2 keV fluxes. The lower panels provide the R-band fluxes, multi-wavelength polarization degree and angle, respectively. The IXPE data are binned over 3 h intervals.

Fig. 5. Multi-wavelength light curves during the multi-hour exposure from MAGIC and XMM-Newton on MJD 60291 (December 13, 2023). The MAGIC fluxes are calculated above 400 GeV and binned over 25-minute intervals. The XMM-Newton fluxes are computed in the same 25-minute bins as MAGIC, and in the 0.3–2 keV and 2–10 keV bands. In the third panel from the top we show the hardness ratio (F2 − 10 keV/F0.3 − 2 keV) from the XMM-Newton data. The bottom panels show the optical (R-band) fluxes, and the polarization in the radio and optical (R-band) with the same binning as in Fig. 1. The black markers represent the X-ray polarization from IXPE binned over 3 h intervals.

The MAGIC data show a nonsignificant (< 3σ) flux variability. By fitting the data with a constant model, we obtain χ²/d.o.f. = 17.2/11, implying that the hypothesis of a nonvariable emission is rejected at a significance of only 1.3σ. Differently, the XMM-Newton light curves reveal significant variability, although with a moderate flux amplitude at the level of 10–15% in both bands. The hardness ratio also varies significantly, implying sub-hour spectral variability in the X-rays. From a constant fit, we obtain χ²/d.o.f. = 96.5/13, and the hypothesis of a constant spectral behavior is rejected at a significance of 7.8σ. The hardness ratio is lower during the rising phase of the flux than during the decaying phase, hinting toward a delay of the high-energy flux compared to the low-energy flux. As can be seen in the light curve, the 2–10 keV band reach its maximum about 15 min after the 0.3–2 keV band.

Figure 6 presents the hardness ratio plotted versus the 2–10 keV flux. Gray arrows give the direction of time. The data show an evident loop in counter-clockwise direction. This hysteresis pattern is indicative of a delay of the higher energy flux with respect to the lower energy one. As discussed in Kirk et al. (1998) and in Sect. 6, such a lag is likely the manifestation of a regime where the radiation cooling timescale is comparable to the particle acceleration timescale.

Fig. 6. Hardness ratio (F2 − 10 keV/F0.3 − 2 keV) as function of the 2-10 keV flux during the XMM-Newton observations on December 13, 2023 (MJD 60291). The markers are color-coded with a gradient from dark blue (first time bin) to light green (last time bin). The gray arrows depict the direction of time, and show a clear counter-clockwise loop pattern.

In Fig. 5, the IXPE polarization measurements are binned over 3 h intervals to search for potential short timescale variability contemporaneous to the long MAGIC/XMM-Newton exposure. The polarization degree is not varying significantly, and is stable around a value of 20–25%. This is a factor of ≈8–9 higher than in the optical band. As for polarization angle, some slow variations in the order of 20° are visible. The average X-ray polarization angle exhibits an offset by ≈60°–70° with respect to the optical and radio polarization angles (which are within ≈10° from each other).

4. Correlated variability

This section first reports on the quantification of the intra-band flux correlations, focusing on the VHE and X-rays. Owing to the unprecedented variability seen in the IXPE data, we also searched for potential correlation of the X-ray polarization properties with the flux and the behavior in other bands.

4.1. Intra-band flux correlation

We investigated the correlation between the VHE and X-ray fluxes in light of their strong variability throughout the campaign, as well as their expected correlated variability within leptonic models (Giebels et al. 2007; Fossati et al. 2008). To do so, we correlated the MAGIC and Swift-XRT fluxes, considering all possible combinations between the respective sub-energy bands (i.e., 0.2–1 TeV and > 1 TeV for MAGIC, and 0.3–2 keV and 2–10 keV for Swift-XRT; see Fig. 1). We matched pairs of measurements that are within a time window of less than 3 h from each other, that is well below the flux doubling/halving timescales reported at those energies during this campaign. The correlation was then quantified using the discrete correlation coefficient (DCF; Edelson & Krolik 1988) assuming a time-lag t_lag = 0 between the light curves. For completeness and comparison purposes, we also show in each panel the widely used Pearson correlation coefficient R. We note the DCF is generally preferred over the Pearson correlation coefficient when dealing with two uneven time series. The DCF method also has the advantage to naturally incorporate measurement uncertainties (Edelson & Krolik 1988).

The resulting above-mentioned flux-flux plots are shown in Fig. 7. For clarity, black filled markers correspond to observations during the IXPE observing window. In each panel, the DCF and R values are provided together with their significance (in Gaussian σ units). To assess the significance of the correlation, we followed a procedure described in Abe et al. (2024) that is based on the code described in Arbet Engels (2021). In summary, we simulated 10⁵ uncorrelated light curves for each energy bands that have the same sampling and binning as the observations. The light curves were simulated assuming a power spectral density (PSD) that follows a power law with index fixed to −1.4 for all energy bands. Such a PSD index agrees with recent works on Mrk 421 (Abe et al. 2025; Aleksić et al. 2015b; Isobe et al. 2015). The simulated light curves preserve the probability density function of the data using the method from Emmanoulopoulos et al. (2013). Finally, the significance was estimated by comparing the observed DCF and R values with the distribution from the simulated light curves.

Fig. 7. VHE energy flux vs X-ray flux using MAGIC and Swift-XRT observations throughout the multi-wavelength campaign. The MAGIC data are binned nightly with a typical exposure time per night of about 40 min. The Swift-XRT observations are binned per observation (typical exposure time of ≈15 min). Only pairs of measurement within 3 hours are considered. The MAGIC fluxes are computed in the > 1 TeV band (top panels) and in the 0.2–1 TeV band (lower panels). The Swift-XRT fluxes are computed in the 0.3–2 keV (left panels) and 2–10 keV bands (right panels). The black filled squares are the measurements during the IXPE observing window. In each subpanel we display the DCF value (and uncertainty) and the corresponding significance in Gaussian σ units in parentheses. For completeness, the Pearson’s R coefficient and its significance are also given. The significance is determined based on dedicated Monte Carlo simulations (see Sect. 4 for more details).

A positive VHE/X-ray correlation is typical in Mrk 421, and has been reported multiple times both during low (see, e.g., Aleksić et al. 2015b; Baloković et al. 2016; Acciari et al. 2021) and high activity (see, e.g., Aleksić et al. 2015a; Acciari et al. 2020), suggesting a common underlying particle population. During some flaring states, a tighter correlation of the X-rays with the ≳1 TeV band than with the 0.2–1 TeV band (see, e.g., Acciari et al. 2020) has been found. However, this does not happen for all flares (see, e.g., MAGIC Collaboration 2021). Interestingly, during the 70-day campaign in 2023-2024 considered in this study, there is substantial variability in both X-rays and VHE gamma rays, but the correlation between the various bands in X-rays and VHE gamma rays is only marginally significant (see inlaid information in the panels from Fig. 7). The magnitude of the DCF and Pearson correlation is high (> 0.7) for most of the bands, but the significance of these values is not large due to the large scattering in the data points. There is an apparent overall correlation pattern in all the panels. Nonetheless, there are some groups of data points that are far away from the overall pattern. For instance, within the IXPE window (see black filled markers), one can see a few observations where the > 1 TeV (0.2–1 TeV) flux varied by more than a factor of 2 (1.5) for a very similar X-ray flux state. Such deviations from the overall trend, which indicate different emission mechanisms at work, increase the calculated error for the DCF value, and naturally also decrease the significance of the correlations.

4.2. Polarization

We searched for possible correlations between the X-ray polarization degree and angle, as well as between the X-ray polarization properties and flux. We also correlated the X-ray polarization with the VHE flux (0.2–1 TeV). The latter investigation is motivated by the fact that IXPE probes the 2–8 keV band that is originating from electrons with similar energies as the one that are also responsible for the VHE emission within leptonic models (Giebels et al. 2007; Fossati et al. 2008). The correlation was conducted using the IXPE data binned into 6 h and 12 h intervals as presented in Fig. 1. No significant (> 3σ) evidence for an underlying correlation pattern is found, neither with any of the X-ray and VHE bands (see Appendix C, Figs. C.1 and C.2), nor between the polarization degree and angle themselves (see Appendix C, Fig. C.3). The X-ray polarization (degree and angle) also does not show a significant correlation with the X-ray or VHE hardness ratios (see bottom panels in Figs. C.1 and C.2). In summary, the X-ray polarization behaves erratically, and the variations do not seem to be related to the flux level nor the spectral shape.

5. Theoretical description of the flare evolution using constrains from the X-ray polarization measurements with IXPE

We modeled the broadband evolution during the IXPE observing window assuming a leptonic scenario (Maraschi et al. 1992; Ghisellini & Maraschi 1996; Tavecchio et al. 1998). Leptonic scenarios assume that the radio to X-ray originates from electron-synchrotron radiation, while the gamma-ray photons are produced by electron inverse-Compton scattering off the synchrotron photons. Such models are supported by the correlation X-ray/VHE that has been observed in the emission of Mrk 421, even if the period considered here led to only marginally significant correlations (despite the relatively large magnitude in the correlation), as reported in Fig. 7. For this, we first computed daily SEDs from radio to VHE around each day that includes a MAGIC observation. The MAGIC SEDs were extracted after correcting for the EBL absorption using the model of Domínguez et al. (2011). We complemented the MAGIC SEDs with strictly simultaneous X-ray data from Swift-XRT and also from XMM-Newton during the long exposure on December 13, 2023 (MJD 60291). The Fermi-LAT observations were averaged over 3-day intervals centered at the MAGIC observing time because of the limited sensitivity to resolve the spectrum of Mrk 421 below daily timescales and the low variability in this waveband. We included R-band and UV observations that are the closest in time to MAGIC. Due to the sparse sampling of Swift-UVOT (Fig. 1), some UV observations have a time difference reaching a maximum of 8.5 days with respect to the MAGIC one. Due to the moderate variability in those bands (see Figs. 1 and 2), they can nevertheless be considered as a good proxy for the UV emission during the X-ray and VHE gamma-ray measurements. Finally, we added contemporaneous radio data (simultaneous on timescales of several days) for completeness, they are not included in our model given that they may receive significant contribution from multiple regions (not considered in our model for simplicity) farther downstream the jet. In total, this provides us with nine SEDs (between MJD 60287 and MJD 60299, December 9, 2023, to December 21, 2023), which we modeled within the scenario described below. We used the same code as in MAGIC Collaboration (2021, 2024), which employs routines from the naima package (Zabalza 2015) to compute the synchrotron and inverse-Compton emissivities.

Our model aims at qualitatively capturing the radio-to-X-ray polarization properties that suggest an energy stratification of the emitting region (see Sect. 3.1). In a similar approach to MAGIC Collaboration (2024), we assumed a morphology consisting of two overlapping, spherical zones: a “compact” zone near the acceleration site, and an “extended” zone that occupies a larger volume downstream the jet. The compact zone contains freshly accelerated and energetic electrons that dominate the X-ray and > 100 GeV emission. The extended zone, populated by less energetic electrons, is dominant in the UV/optical and < 100 GeV bands. In line with Liodakis et al. (2022), Angelakis et al. (2016), Marscher (2014), the observed higher polarization degree in the X-ray band compared to the radio/optical is due to the confinement of the freshly accelerated electrons in a smaller region near the shock front (mimicked by the compact zone), which compresses the magnetic field leading to a more ordered magnetic field structure. When electrons subsequently cool and advect toward larger regions (mimicked here by the extended zone) in which the degree of magnetic field turbulence increases significantly, a drop of the polarization degree at lower frequencies is expected. Given the very different polarization variability patterns observed in the optical and X-ray bands, it is required that the compact zone remains subdominant in the UV/optical compared to the extended zone. Finally, since we assumed that the two regions are overlapping, our code also includes the emission resulting from the interaction of the two zones since the (synchrotron) emission from the two regions provide an additional target photon field for each other for the inverse-Compton scattering. The up-scattering of the extended zone field by the electrons in the compact zone leads to an inverse-Compton luminosity comparable to the one resulting from the up-scattering of the compact zone field by the extended zone electrons. The only difference is that the two components are slight shifted in frequency due to the different electron and photon energies populating each zone.

Our assumption of using two distinct components is naturally a simplification of the reality because one would rather expect a continuum of several regions contributing to the different parts of the SED. A detailed treatment of the turbulence, electron advection and diffusion is beyond the scope of this work.

We exploit the polarization data to constrain the relative size between the compact and extended zones for the nine SEDs modeled here. We assume that each of the region is made of N turbulent plasma cells of identical magnetic field strength, but with random orientation. In such a configuration, the expected average polarization degree from each zone can be approximated as $P_{\mathrm{deg}}\approx 75\%/\sqrt{N}$ (Marscher 2014; Tavecchio 2021). The measured ratio between the X-ray and optical polarization degree (P_{X − ray, deg}/P_opt, deg; see Fig. 3) was then used to estimate the relative difference, l, between the number of turbulent cells in the extended zone and the compact zone. Under the assumption that cells in both zones roughly span the same spherical volume, one derives a relative difference in radius of l^1/3 between the two regions. Here, we fixed the radius of the extended zone to 1.5 × 10¹⁶ cm (in line with constraints from the light crossing time), and then determined the radius of the compact zone using the observed ratio P_{X − ray, deg}/P_opt, deg during each day. The choice of fixing the radius of the extended zone is motivated by the low (flux and polarization) variability in the UV/optical bands.

Fig. 8. Theoretical models for of the broadband SEDs during the IXPE observing window, from MJD 60287 to MJD 60299 (December 9, 2023, to December 21, 2023). The observations, spanning from the radio to VHE, are shown by the orange markers. The dashed violet curve is the emission originating from the compact zone, located near the shock front. The dot-dashed blue line is the contribution from the extended zone, which spans a larger volume downstream the shock. The emission produced by the interaction between the two zones is plotted as a pink dotted curve. The solid blue line shows the sum of all the components. For comparison purposes, we plot with gray markers the average state of Mrk 421 (taken from Abdo et al. 2011). The values for the model parameters can be found in Tables 2 and 1 (see Sect. 5 for more details on the modeling procedure).

In order to limit the degrees of freedom of the used theoretical model, we also made the following assumptions (primed quantities are expressed in the source reference frame):

• The Doppler factor δ was fixed to 60 in the compact zone. Such a Doppler factor is somewhat larger than those typically adopted in the literature for Mrk 421, which are in the range ∼20–40 (Abdo et al. 2011; Baloković et al. 2016; Banerjee et al. 2019), although δ ≈ 60 is still consistent with the modeling of the March 2008 flare presented in Aleksić et al. (2012). The rather soft X-ray spectrum gives rise to a large separation between the two SED peak frequencies, hence requiring a high Doppler factor within a leptonic scenario. Using values significantly lower than 60 (e.g., by a factor of two) leads to VHE spectra softer than the observed ones, and hence we excluded them. To illustrate this, Appendix D reports a modeling of the broadband SED from MJD 60291 (December 13, 2023), which corresponds to the day with simultaneous MAGIC/XMM-Newton observations, using δ = 30. We would like to point out that δ values below 60 for the extended zone would describe satisfactorily the optical and Fermi-LAT data in the part of the SED dominated by the extended zone. However, for the sake of simplicity, and to limit the number of free parameters, we set δ = 60 in both zones. We note that the radius of the two zones differs by less than a factor of 4, and that the usage of very different Doppler factors between the two zones would imply a strong deceleration of the flow on spatial scales comparable to the size of the emission, that is ∼10¹⁶ cm (∼10⁻² pc). This is significantly smaller than what is observed in FR I galaxies (Hardcastle et al. 2005).

• In both regions, the electron distribution follows a power-law distribution with exponential cutoff: $\mathrm{d}n/\mathrm{d}\gamma \prime \propto {\gamma }^{\prime -p}exp(-\gamma \prime /\gamma {\prime }_{\mathrm{cutoff}})$, for $\gamma {\prime }_{\min }<\gamma \prime <\gamma {\prime }_{\max }=+\infty$. In the compact zone, the slope can be well constrained by the X-ray data. For the extended zone, p was fixed to the canonical value of 2.2 for relativistic shocks (Kirk et al. 2000) given that the UV/optical and MeV/GeV coverage does not provide a strong constrain on this parameter.

• In the compact zone, we fixed $\gamma {\prime }_{\min }=2\times {10}^4$ and $\gamma {\prime }_{\mathrm{cutoff}}=1.3\times {10}^5$. Given the magnetic field strength necessary to adequately reproduce the SEDs (≈0.05–0.08 G), such a large $\gamma {\prime }_{\min }$ is needed so that the synchrotron flux from the compact zone (responsible for the X-rays) becomes subdominant at lower energies in the UV/optical. $\gamma {\prime }_{\mathrm{cutoff}}$ was selected to provide a good description of the X-ray and VHE spectral shapes at the highest energies for all days.

• The electron distribution in the extended zone is defined with $\gamma {\prime }_{\min }=2\times {10}^3$ and $\gamma {\prime }_{\mathrm{cutoff}}=2.7\times {10}^4$. $\gamma {\prime }_{\min }$ is limited from below in order to not overshoot the radio data in the 225.5 GHz band, while $\gamma {\prime }_{\mathrm{cutoff}}$ is set such that the extended zone does not dominate the overall emission beyond the UV band.

Based on those assumptions, the only remaining free parameter in our model for the extended zone is the magnetic field B′ and the electron density N′_e. In the compact zone the remaining free parameters are B′, N′_e, and p, which we varied for each day in order to properly describe the X-ray and VHE data. The results are shown in Fig. 8. The observations are plotted with orange markers. The contribution from the extended and compact zones are shown in blue and violet solid lines, respectively. The inverse-Compton emission resulting from the interaction between the two zones is plotted with a pink dotted line and the sum of all components is depicted with a solid light blue curve. For all days a good description of the observations is achieved. We summarize in Table 2 the parameters of the two emitting zones that are fixed, and in Table 1 the evolving parameters in the compact and the extended component.

The model does not evolve the particle population self-consistently by taking into account the radiation cooling. Nonetheless, according to the B′ strength required to well reproduce the SED and assuming that the particles escape the regions on timescales $t{\prime }_{\mathrm{esc}}\approx R\prime /c$, the synchrotron cooling break is above $\gamma {\prime }_{\mathrm{cutoff}}$ by more than a factor of a few (for both zones). The system is therefore in a slow-cooling regime (i.e., the particles likely escape the emitting region prior to any significant cooling). The cooling effects are thus neglected when modeling the SED evolution.

Figure 9 shows the relative evolution for each of the parameters in the compact zone (solid lines) and in the extended zone (dashed lines). The radius, R′, and the electrons density within the zone, N′_e, are the two parameters that vary the most. The radius, determined from the optical-to-X-ray polarization degree varies up to a factor of 2 along the flare, and the density evolves by a factor of 3 at least. The magnetic field, B′, changes moderately, by less than 40%. Regarding the slope of the electrons p, it varies by less than 10%, which is expected given the moderate spectral variability reported in Sect. 3.2. We provide in Sect. 6 a physical interpretation and a more detailed discussion of the proposed model.

Fig. 9. Relative evolution of the evolving parameters in the compact zone (solid lines) and extended zone (dashed line). Each parameter, plotted in a different color, is normalized to its average value.

6. Discussion and summary

In this work we studied, for the first time for a HSP blazar, a bright state of Mrk 421 from radio to VHE with simultaneous X-ray polarization measurements using IXPE. The 0.2–1 TeV flux reached a peak activity slightly above 2 C.U., which is more than a factor of 4 higher than the average state at those energies (Acciari et al. 2014). A bright X-ray flux is also detected in the 0.3–2 keV band, and is comparable to the level recorded during archival flaring states, such as the one in March 2010 (Aleksić et al. 2015a). Until now, multi-wavelength studies including X-ray polarization have only probed HSPs in low or quiescent states (see, e.g., Abe et al. 2024; MAGIC Collaboration 2024). Therefore, the study presented here, with data from this campaign, offers the most complete characterization of a HSP blazar during a period of high activity to date, hence providing crucial insights into the physical origin of blazar flares.

6.1. Interpretation of the multi-wavelength polarization variability during the VHE flare

The polarization degree shows significant chromatic behavior, indicating an energy stratification of the emitting zone. The average X-ray polarization degree is about three times higher than the R-band, and about seven times higher than in the radio. Such an energy dependence is comparable to previous results on Mrk 421, and also for the HSP Mrk 501 (Liodakis et al. 2022; Di Gesu et al. 2022b). However, the drastic variations (by a factor of ∼8) of the X-ray polarization degree happening down to ∼6 h timescales is a novel trend not yet found in earlier observations. In addition, the X-ray polarization angle varies significantly by about 100° on day timescales. These variations of the X-ray polarization do not reveal any evident structured pattern and are consistent with a stochastic trend. This differs from the steady rotation seen in 2022 (Di Gesu et al. 2023) that occurred over about 5 days with a roughly steady angular velocity and constant polarization degree. We note however that the average X-ray polarization angle during December 2023 aligns closely with that in the optical and with the radio jet axis (Weaver et al. 2022). This suggests a partial ordering of the magnetic field perpendicular to the jet, possibly caused by the presence of a shock compressing the helical field structure (Marscher et al. 2017; Tavecchio et al. 2018). The stochastic nature of the polarization angle and degree evolution is then likely due to the plasma being highly turbulent before crossing the shock front (Marscher 2014). We note that magnetic reconnection is also expected to generate rapid and strong variability of the polarization (Zhang et al. 2022), but the polarization angle is expected to reach any orientation with respect to the jet. However, the average angle measured with IXPE remains closely aligned with the jet and the optical band, suggesting that magnetic reconnection does not dominate the acceleration of the electrons responsible for the X-ray and VHE gamma-ray emission measured during these IXPE observations.

The highest X-ray polarization degree of the campaign is ≈26% on MJD 60290 (December 12, 2023). It is close to the record value measured to date in a HSP, being ≈30% in PKS 2155-304 (Kouch et al. 2024). The following day (MJD 60291, December 13, 2023), the polarization was similarly high, ∼23%, and about eight times higher than that in the optical. Simultaneously with this high polarization state, the multi-hour exposure from XMM-Newton reveals a hysteresis spectral loop in the counter-clockwise direction. Despite being a pattern observed in earlier XMM-Newton observations of Mrk 421 and other HSPs (Zhang et al. 2002; Brinkmann et al. 2003; Ravasio et al. 2004), it is the first time we can study it in combination with the recorded high X-ray polarization degree coupled with the strong chromatic behavior.

Counter-clockwise loops are indicative of a delay from the high-energy X-ray band (here 2–10 keV) with respect to the lower energy band (here 0.3–2 keV). It can be explained by the gradual acceleration of particles toward higher energy, first reaching the 0.3–2 keV radiation band and then the 2–10 keV band. As discussed in Kirk et al. (1998), such a high-energy delay becomes apparent when the radiation cooling timescale is comparable to the particle acceleration timescale, $t{\prime }_{\mathrm{cool}}\sim t{\prime }_{\mathrm{acc}}$, which occurs at the highest particle energy reachable by the system. The IXPE energy range, which largely overlaps with the 2–10 keV band, therefore probes during that day the emission from particles that are located at (or very close to) the high-energy cutoff. In a shock acceleration scenario, the most energetic particles are located close to the shock front, and probe the highest degree of magnetic field ordering due to the compression of the helical field (Marscher & Gear 1985; Tavecchio et al. 2018). A high X-ray polarization degree is thus expected simultaneously with counter-clockwise hysteresis loops, which is supported by the data from our campaign.

The polarization degree in the R-band simultaneous with the XMM-Newton observation on MJD 60291 (December 13, 2023) is drastically lower, by a factor of ≈7, and is the largest difference observed in this campaign (see Fig. 3). The lower polarization at decreasing frequencies is possibly caused by the advection and cooling of the freshly accelerated particles toward the downstream regions of the jet where a higher degree of turbulence exists (Marscher & Gear 1985; Tavecchio et al. 2018). The TEMZ model developed by Marscher (2014), which considers a number of turbulent plasma cells crossing a standing conical shock, also predicts an increase in the polarization with frequency. It is caused by the higher acceleration efficiency in plasma cells whose magnetic field is almost parallel to the shock normal. Since the cells fulfilling such conditions are less numerous than those whose magnetic field is oriented away from the shock normal, the polarization increases at higher frequencies ($P_{\mathrm{deg}}\propto 1/\sqrt{N}$, where N is the number of cells). In such a situation, a stronger variability of the polarization at higher frequencies is also expected, in agreement with the observations reported here. Based on TEMZ simulations presented in Marscher & Jorstad (2021) the X-ray polarization is expected to be on average ≈5% higher (in absolute terms) than the one in the R-band. Although this chromaticity is lower than observed in December 2023 for several of the days, and in particular during the XMM-Newton observations on MJD 60291 (December 13, 2023). Nonetheless, since the counter-clockwise loops indicate that IXPE probed the emission from particles at the high-energy cutoff, we probed an extreme case of the TEMZ model for which the chromatic trend may be more pronounced. This extreme situation would need to be investigated with a dedicated simulation to verify if TEMZ model is still in agreement with the data.

6.2. Implications of the observed variability in the X-ray and VHE bands

We find a > 3σ indication of a positive X-ray/VHE correlated variability (see Fig. 7). This is a general trend regularly reported for Mrk 421 and other HSPs (see, e.g., Abe et al. 2023; Arbet-Engels et al. 2021; MAGIC Collaboration 2020; Aharonian et al. 2009) that indicates a co-spatial origin and a common underlying population of radiating particles, supporting a leptonic origin of the emission. Nonetheless, a large scatter is measured in the correlation throughout the IXPE observing window. For a roughly equivalent X-ray flux, the VHE flux varies by almost a factor of 2, implying an evolution with time of the Compton dominance (the ratio between the inverse-Compton and synchrotron peak luminosities). As discussed in Katarzyński et al. (2005), the VHE/X-ray correlation in HSPs may exhibit very different trends (ranging from a sublinear to more-than-quadratic relationship) depending on the parameters driving the variability. The observed scatter may thus be caused by simultaneous and uncorrelated changes in multiple characteristics, such as the size of the emitting region, the magnetic field, and the particle density varying simultaneously in an uncorrelated manner.

We highlight that changes in the Compton dominance and a significant scatter in the VHE/X-ray correlation are naturally expected (and have been observed, see, e.g., MAGIC Collaboration 2021) when considering datasets covering years or months. The reason is that, on such timescales, separate regions from the jet with different environments may dominate one after the other the broadband emission (see, e.g., Hervet et al. 2019). In the present case the evolution occurs down to daily timescales, which is much shorter. It further supports the scenario in which the flare is caused by a highly turbulent plasma crossing a shock front since its turbulent nature may easily explain simultaneous rapid (and stochastic) variations of multiple source environment parameters, which in turn produce a scatter in the correlation patterns.

Regarding the VHE and X-ray spectral characteristics, rather soft spectra are observed (see Sect. 3.2). For a given X-ray and VHE flux, the hardness ratio is on average lower compared to archival data. Mrk 421 commonly shows a harder-when-brighter behavior in VHE and X-rays (Acciari et al. 2021; MAGIC Collaboration 2021), but this trend is only partly observed in the current campaign (see Fig. 4). While the hardness ratio over the entire period indeed shows a harder-when-brighter evolution, it becomes nonsignificant during the IXPE window. Focusing on the hardness ratio as function of the 0.3–2 keV flux, the observations during the IXPE window show a clearly diverging pattern with respect to the rest of the campaign, and the spectral shape stays roughly constant. The absence of harder-when-brighter trend could again be indicative of a highly turbulent plasma crossing a shock front. For instance, if there were stochastic and uncorrelated changes in the electron density and the size of the emitting region, any harder-when-brighter trend would likely be washed out.

6.3. Interpretation of the SED modeling within the context of a shock acceleration scenario

We modeled the evolution of the broadband SED across the IXPE window. The model consists of two overlapping components, one dominating the X-ray and VHE (the compact zone) and another dominating the UV/optical and MeV/GeV regimes (the extended zone). The radio band is dominated by the extended zone, but the model remains below the observed fluxes. We implicitly assumed that a significant fraction of the radio flux receives an additional contribution from broader regions farther downstream the jet. Both zones were assumed to be composed of N turbulent plasma cells that have a roughly identical magnetic field strength but a random orientation, leading to an averaged observed polarization degree of $P_{\mathrm{deg}}\approx 75\%/\sqrt{N}$. The relative size of the components were then determined based on the optical-to-X-ray polarization degree. The compact zone is located close to the acceleration site and contains freshly accelerated particles, where the magnetic field is more ordered. When particles cool and advect, they populate a larger region (mimicked by the extended zone) where a higher degree of turbulence induces a drop in the polarization.

The large X-ray polarization variability that is not observed in the optical, as well as the chromatic trend of the polarization degree, indicate that the compact zone is largely subdominant in the optical regime. Using the magnetic field that we derived from the model (≈0.05–0.08 G), this implies a high minimum Lorentz factor $\gamma {\prime }_{\min }\gtrsim {10}^4$. This is significantly higher that the values used in previous modeling of BL Lac-type object SEDs, which are usually $\gamma {\prime }_{\min }\sim {10}^2-{10}^3$ (Tavecchio et al. 2010; Baloković et al. 2016). In the case of extreme TeV BL Lac objects (that have an inverse-Compton peak located above 1 TeV) values around 10⁴ are sometimes required to model the hard VHE spectrum (Katarzyński et al. 2006; Costamante et al. 2018). In our present case, the high minimum Lorentz factor is robustly implied by the polarization data. A large minimum Lorentz factor is expected in an electron-ion plasma. While in such a scenario most of the energy is carried by the ions, a fraction of it can be transferred to the electrons, leading to a large $\gamma {\prime }_{\min }$ in the electron population. In the case of relativistic shocks, particle-in-cell (PIC) simulations indicate that $\gamma {\prime }_{\mathrm{electrons},\min }\sim 600{\gamma }_{\mathrm{sh}}$, where γ_sh is the shock velocity in the un-shocked plasma frame (Sironi & Spitkovsky 2011; Sironi et al. 2013; Zech & Lemoine 2021). Hence, the data suggest accelerations at highly relativistic shocks (γ_sh ≳ 10). The requirement of a relativistic shock is also in line with the relatively high Doppler factor (δ ∼ 60) needed to capture the X-ray/VHE spectrum.

In the model, the cutoff Lorentz factor in the compact zone is fixed to $\gamma {\prime }_{\mathrm{cutoff}}=1.3\times {10}^5$ for all days, and provides a good description of the X-ray/VHE SEDs. $\gamma {\prime }_{\mathrm{cutoff}}$ cannot be precisely constrained with the data at hand, but it is constrained from below by the X-ray spectrum since values significantly lower than 1.3 × 10⁵ cannot accommodate the highest X-ray energy points. At $\gamma {\prime }_{\mathrm{cutoff}}=1.3\times {10}^5$, the electrons radiate synchrotron photons in the ≈1–3 keV band (according to the magnetic field value used in the model). Thus, in agreement with the hysteresis pattern in the counter-clockwise direction discussed above, the X-ray data from IXPE, XMM-Newton, and Swift-XRT are probing the cutoff energy of the accelerated electrons.

The temporal evolution of the SEDs can be well described by varying in both zones the electron density N′_e and magnetic field B′, as well as as the slope of the injected electrons in the compact zone in order to capture the X-ray/VHE spectral variability. The electron density is the parameter that varies the most, up to a factor of 3 (see Fig. 9), while the magnetic field changes by less than 40%. We do not find any correlation between the fluxes and the evolution of any of the parameters, which again corroborates the stochastic nature of the physical mechanism driving the flare.

Our model shows that B′ within the compact zone tends to increase for smaller radii (see Fig. 9) when the emitting particles are on average closer to the acceleration site. During the XMM-Newton observation (MJD 60291, December 13, 2023), the elevated X-ray polarization degree observed by IXPE leads to the smallest region radius among all the days based on the constrains described earlier (R′ = 0.4 × 10¹⁵ cm), and in order to properly describe the data, the magnetic field reaches a maximum (B′ = 0.08 G). We also find that the magnetic field is always stronger in the compact zone compared to that of the extended zone. PIC simulations show that the magnetic field at a shock front is expected to be self-amplified before decaying farther downstream. The observed anti-correlation between B′ and R′ in our modeling is thus consistent with a shock-acceleration scenario (Sironi & Spitkovsky 2011; Tavecchio et al. 2018). Finally, we highlight that an anti-correlation between the radius and magnetic field is also expected to be caused by the cooling of the particles. Assuming that synchrotron radiation is the dominant cooling mechanisms, an increase in B′ will make the high-energy particles radiate their energy closer to the acceleration site, hence making the size of the emitting region smaller.

The slope of the electrons in the compact zone vary moderately throughout the days, max(p) − min(p) = 0.45. This is in line with the relatively moderate X-ray/VHE spectral variability discussed previously. We find that p ≈ 2.2 − 2.5 provides a good description of the majority of the X-ray and VHE spectra. Such slopes are close to the predictions from relativistic parallel shocks (i.e., ≈2.23), and values close to ≈2.5 can easily be produced, for instance when the shock is oblique instead of purely parallel (Kirk et al. 2000; Summerlin & Baring 2012).

To conclude, this paper reports the first measurement of X-ray polarization from a HSP blazar during an X-ray and gamma-ray flaring activity, which has provided unprecedented information on the acceleration and radiation mechanisms driving blazar flares. The fast variations of the X-ray polarization, combined with an average angle remaining closely aligned with that of the jet (and in the optical), favors a scenario in which the flare is caused by a highly turbulent electron-ion plasma crossing the shock front where the particles get accelerated. The measured increase in the polarization degree with frequency further implies that the emitting region is stratified, and that the most energetic particles are located closest to the shock. The multi-wavelength spectral properties, intra-band correlation, and the counter-clockwise loop in the hardness ratio of the X-ray emission supports this scenario, which can also be described, to first order, within a theoretical leptonic model based on two components. Further multi-wavelength observations that include IXPE during flaring episodes of blazars will be crucial to determine whether the observed trends are common characteristics of blazar flares.

Data availability

Appendix A, B, C, and D, including their Tables and Figures, are available in Zenodo at the following link: https://zenodo.org/records/14781432.

Acknowledgments

A. Arbet Engels: project management, P.I. of MAGIC observations, organization of multi-wavelength observations and data analysis, MAGIC analysis, correlation analysis, theoretical modeling and interpretation, paper drafting; I. Liodakis: organization of multi-wavelength observations and data analysis; L. Heckmann: project management, MAGIC and Fermi-LAT 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. We 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, PID2022-136828NB-C41, PID2022-137810NB-C22, PID2022-138172NB-C41, PID2022-138172NB-C42, PID2022-138172NB-C43, PID2022-139117NB-C41, PID2022-139117NB-C42, PID2022-139117NB-C43, PID2022-139117NB-C44 funded by the Spanish MCIN/AEI/ 10.13039/501100011033 and “ERDF A way of making Europe”; 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 nr. 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 (CEX2019-000920-S, CEX2019-000918-M, CEX2021-001131-S) and by the CERCA institution and grants 2021SGR00426 and 2021SGR00773 of the Generalitat de Catalunya; by the Croatian Science Foundation (HrZZ) Project IP-2022-10-4595 and the University of Rijeka Project uniri-prirod-18-48; by the Deutsche Forschungsgemeinschaft (SFB1491) and by the Lamarr-Institute for Machine Learning and Artificial Intelligence; by the Polish Ministry Of Education and Science grant No. 2021/WK/08; and by the Brazilian MCTIC, CNPq and FAPERJ. 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’Energie 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 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). I. L. and S.K. were funded by the European Union ERC-2022-STG – BOOTES – 101076343. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them Some of the data are based on observations collected at the Observatorio de Sierra Nevada; which is owned and operated by the Instituto de Astrofísica de Andalucía (IAA-CSIC); and at the Centro Astronómico Hispano en Andalucía (CAHA); which is operated jointly by Junta de Andalucía and Consejo Superior de Investigaciones Científicas (IAA-CSIC). The research at Boston University was supported in part by National Science Foundation grant AST-2108622, NASA Fermi Guest Investigator grants 80NSSC21K1917 and 80NSSC22K1571, and NASA Swift Guest Investigator grant 80NSSC22K0537. This research was conducted in part using the Mimir instrument, jointly developed at Boston University and Lowell Observatory and supported by NASA, NSF, and the W.M. Keck Foundation. We thank D. Clemens for guidance in the analysis of the Mimir data. This study used observations conducted with the 1.8m Perkins Telescope (PTO) in Arizona (USA), which is owned and operated by Boston University. This work was supported by NSF grant AST-2109127. 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 work has made use of data from the Joan Oró Telescope (TJO) of the Montsec Observatory (OdM), which is owned by the Catalan Government and operated by the Institute for Space Studies of Catalonia (IEEC). 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. We recognize that Maunakea is a culturally important site for the indigenous Hawaiian people; we are privileged to study the cosmos from its summit. The 100 m radio telescope at Effelsberg is operated by the Max-Planck-Institut für Radioastronomie (MPIfR) on behalf of the Max-Planck-Society. Observations with the 100 m telescope have received funding from the European Union’s Horizon 2020 research and innovation program under grant agreement No. 101004719 (ORP). The Liverpool Telescope is operated on the island of La Palma by Liverpool John Moores University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias with financial support from the UK Science and Technology Facilities Council.

References

All Tables

All Figures

Fig. 1. Multi-wavelength light curve between November 7, 2023 (MJD 60255), and January 16, 2024 (MJD 60325). The red line marks the XMM-Newton observation taking place on December 13, 2023 (MJD 60291), and the vertical gray bands highlight the IXPE observing windows. From top to bottom: MAGIC fluxes in the > 1 TeV (first panel) and 0.2–1 TeV (second panel) bands with nightly bins (the typical nightly exposure is ≈40 min). The gray dashed line depicts the average flux of Mrk 421 from Acciari et al. (2014); Fermi-LAT fluxes in three-day bins; X-ray fluxes binned per observation including Swift-XRT and XMM-Newton; Swift-UVOT and optical R-band data using the telescopes listed in Sect. 2; radio data from the SMAPOL and QUIVER programs; polarization degree and polarization angle observations in the X-rays from IXPE (in black and empty gray markers the data are binned over 12 h and 6 h bins, respectively), the optical R-band and the radio band. In the polarization angle panel, we plot with a horizontal blue band the average direction angle determined at 43 GHz by Weaver et al. (2022) (see text in Section 3 for further details).

In the text

Fig. 2. Fractional variability, Fvar, for the light curves displayed in Fig. 1. Nightly bins are used for MAGIC, three-day bins for Fermi-LAT, and for all other wavebands the single observations are without further binning. The filled orange to violet markers depict Fvar when data over the full campaign are considered, and the open turquoise markers only include measurements during the IXPE time window. For Fermi-LAT, Fvar cannot be computed during the IXPE window due to a variability level below the statistical uncertainties, leading to a negative value in the square root of the numerator of Fvar (see Vaughan et al. 2003).

In the text

	Fig. 3. Ratio of the X-ray (2–8 keV) and radio (4.85 GHz) polarization degree to that measured in the R-band during the IXPE observing window (PX − ray, deg/PR − band, deg and Pradio, deg/PR − band, deg, respectively). The vertical dashed line marks the date of the simultaneous XMM-Newton/MAGIC long exposure.
In the text

Fig. 4. Hardness ratio in the VHE (upper panels) and X-rays (lower panels) bands. The VHE hardness ratio is defined as the ratio of the > 1 TeV to the 0.2–1 TeV fluxes, while in the X-ray it is the ratio of the 2–10 keV to the 0.3–2 keV fluxes. They are plotted vs the 0.2–1 TeV and > 1 TeV bands, while in the X-ray they are plotted vs the 0.3–2 keV and 2–10 keV bands. The colored circular markers correspond to measurements from the campaign under study, and highlight the data simultaneous with the IXPE window with black-filled markers. Archival measurements from the extensive multi-instrument observing campaigns in the years 2015–2016 (Acciari et al. 2021), 2017 (MAGIC Collaboration 2021), and 2022 (Abe et al. 2024) are plotted as gray triangles, gray squares, and black diamonds, respectively.

In the text

Fig. 5. Multi-wavelength light curves during the multi-hour exposure from MAGIC and XMM-Newton on MJD 60291 (December 13, 2023). The MAGIC fluxes are calculated above 400 GeV and binned over 25-minute intervals. The XMM-Newton fluxes are computed in the same 25-minute bins as MAGIC, and in the 0.3–2 keV and 2–10 keV bands. In the third panel from the top we show the hardness ratio (F2 − 10 keV/F0.3 − 2 keV) from the XMM-Newton data. The bottom panels show the optical (R-band) fluxes, and the polarization in the radio and optical (R-band) with the same binning as in Fig. 1. The black markers represent the X-ray polarization from IXPE binned over 3 h intervals.

In the text

	Fig. 6. Hardness ratio (F2 − 10 keV/F0.3 − 2 keV) as function of the 2-10 keV flux during the XMM-Newton observations on December 13, 2023 (MJD 60291). The markers are color-coded with a gradient from dark blue (first time bin) to light green (last time bin). The gray arrows depict the direction of time, and show a clear counter-clockwise loop pattern.
In the text

Fig. 7. VHE energy flux vs X-ray flux using MAGIC and Swift-XRT observations throughout the multi-wavelength campaign. The MAGIC data are binned nightly with a typical exposure time per night of about 40 min. The Swift-XRT observations are binned per observation (typical exposure time of ≈15 min). Only pairs of measurement within 3 hours are considered. The MAGIC fluxes are computed in the > 1 TeV band (top panels) and in the 0.2–1 TeV band (lower panels). The Swift-XRT fluxes are computed in the 0.3–2 keV (left panels) and 2–10 keV bands (right panels). The black filled squares are the measurements during the IXPE observing window. In each subpanel we display the DCF value (and uncertainty) and the corresponding significance in Gaussian σ units in parentheses. For completeness, the Pearson’s R coefficient and its significance are also given. The significance is determined based on dedicated Monte Carlo simulations (see Sect. 4 for more details).

In the text

Fig. 8. Theoretical models for of the broadband SEDs during the IXPE observing window, from MJD 60287 to MJD 60299 (December 9, 2023, to December 21, 2023). The observations, spanning from the radio to VHE, are shown by the orange markers. The dashed violet curve is the emission originating from the compact zone, located near the shock front. The dot-dashed blue line is the contribution from the extended zone, which spans a larger volume downstream the shock. The emission produced by the interaction between the two zones is plotted as a pink dotted curve. The solid blue line shows the sum of all the components. For comparison purposes, we plot with gray markers the average state of Mrk 421 (taken from Abdo et al. 2011). The values for the model parameters can be found in Tables 2 and 1 (see Sect. 5 for more details on the modeling procedure).

In the text

	Fig. 9. Relative evolution of the evolving parameters in the compact zone (solid lines) and extended zone (dashed line). Each parameter, plotted in a different color, is normalized to its average value.
In the text