© The Authors 2025

1. Introduction

Blazars are a class of jetted active galactic nuclei (AGN). The relativistic plasma jet is oriented with a small angle with respect to the line of sight. Blazars emit across the full electromagnetic spectrum, ranging from radio to high-energy (HE; E > 100 MeV) and very-high-energy (VHE; E > 100 GeV) gamma rays. Blazars that show no or very faint emission lines in their optical emission are referred to as BL Lac-type objects (Urry & Padovani 1995).

The broadband emission of BL Lac-type objects is dominated by non-thermal radiation from the jet. The spectral energy distribution (SED) exhibits two large features in the form of two large bumps (see e.g. Abdo et al. 2011). Measurements of spectral and polarization characteristics strongly indicate that the first bump originates from synchrotron radiation produced by relativistic electrons and positrons moving in the magnetic field within the jet. The origin of the second bump is still under debate and is more difficult to determine. Possible leptonic scenarios include electron inverse Compton (IC) scattering on synchrotron photons originating from the first bump, so called synchrotron self-Compton (SSC) (Maraschi et al. 1992; Ghisellini et al. 1998; Madejski et al. 1999), or in certain cases additionally on external target photons (Dermer et al. 1992; Sikora et al. 1994). Hadronic scenarios can also provide explanations for the gamma-ray emission (Mannheim 1993; Mücke & Protheroe 2001; Cerruti et al. 2015). BL Lac-type objects can be classified by their peak frequency of the synchrotron bump (Urry & Padovani 1995; Padovani et al. 2017; Abdo et al. 2010a). Blazars with a peak frequency of ν_s < 10¹⁴ Hz are referred to as low synchrotron peaked blazars (LSPs), with a peak frequency 10¹⁴ Hz < ν_s < 10¹⁵ Hz as intermediate synchrotron peaked blazars (ISPs) and with ν_s > 10¹⁵ Hz as high synchrotron peaked blazars (HSPs).

Markarian 421 (Mrk 421; RA = 11^h4′27.31″, Dec = 38°12′31.8″, J2000, z = 0.031) is an archetypal and very close HSP. It is among the most studied sources in the VHE sky. Mrk 421 was found in states of extreme activity in the X-ray and VHE gamma-ray bands during previous campaigns (e.g. Gaidos et al. 1996; Fossati et al. 2008; Aleksić et al. 2015a; Acciari et al. 2020; Abeysekara et al. 2020). Among all campaigns, 2010 stands out as the most active and violent year seen in the VHE emission from Mrk 421. In February 2010, the Very Energetic Radiation Imaging Telescope Array System (VERITAS) detected the highest VHE flux recorded to date from Mrk 421 (Abeysekara et al. 2020). The highest flux measurement reached a level of ∼27 Crab Units (C.U.) above 1 TeV, making it the brightest VHE gamma-ray flux ever recorded for an AGN. The high photon statistic during such flaring activity allowed for a binning of the flux in 2-minute intervals. A cross-correlation study with the optical band revealed a significant positive correlation with a time lag of ∼25−55 minutes.

Additional correlation studies between the VHE and X-ray bands showed a complex and fast-varying correlation from linear to quadratic behavior down to no correlation at times. Shortly after, in March 2010, a decaying flare was detected by the Florian Goebel Major Atmospheric Gamma Imaging Cherenkov (MAGIC) and VERITAS telescopes (Aleksić et al. 2015a). The detailed and broad multiwavelength (MWL) coverage enabled the construction of 13 consecutive daily SEDs. The emission could be successfully modeled by a single-SSC zone. However, as described in Aleksić et al. (2015a), a better model-data agreement and a more natural explanation of the SED evolution could be achieved with a two-zone scenario. The variation of only a few model parameters could successfully describe the SED variation and made the underlying particle population a very plausible cause of the variable emission.

For the first time, this work exploits the full campaign covering a time period spanning from November 2009 until June 2010 with a broad MWL coverage. High-resolution radio data taken with the Very Long Baseline Interferometry (VLBI) technique add observations of the jet structure to the campaign from mid January 2010 to early July 2012. This wealth of data allows for a detailed study of the variability and correlation behavior of Mrk 421 during its most violent recorded year.

2. Multiwavelength light curves

The study involves data from several instruments covering the emission from radio frequencies to VHE in the time period from November 5, 2009, (MJD 55140) until July 3, 2010, (MJD 55380). The key instruments used are MAGIC, VERITAS, the Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope (Fermi-LAT), the X-Ray telescope (XRT) and the Ultra-violet Optical Telescope (UVOT) on board the Neil Gehrels Swift Observatory (Swift) and the GLAST-AGILE Support Program – Whole Earth Blazar Telescope (GASP-WEBT) network of optical telescopes. Additional telescopes in the X-ray, optical and radio support these. The details of the data analysis of the individual instruments are given in Appendix A.

The top row of Fig. 1 shows the light curves provided by the MAGIC telescopes in the 0.2−1 TeV and > 1 TeV band. For comparison, multiples of the flux of the Crab Nebula¹ are given with dashed-dotted lines. At the beginning of the campaign in November 2009, Mrk 421 was already in an elevated state of activity with an emission of around 1 C.U. in both bands. The average emission state of the source estimated by Whipple over a time span of 14 years is around 0.45 C.U. (Acciari et al. 2014). After a small gap of observations, Mrk 421 showed a strong flaring activity reaching over 3 C.U. above 1 TeV in January 2010. The emission showed strong variability by up to a factor of almost three on a daily timescale. We found significant intranight variability on January 15 (see Appendix C). After a short phase at around 1 C.U., another observational gap due to adverse weather conditions is present in the MAGIC telescope data. Following the gap, a decaying flare starting at around 2 C.U. could be observed, which is discussed in great detail in Aleksić et al. (2015a). Mrk 421 showed average emission for the following two months until it entered an especially low emission state around June, going below 0.25 C.U.

Fig. 1. MWL light curves covering the time period from November 5, 2009, (MJD 55140) to July 3, 2010, (MJD 55380). Top to bottom: MAGIC fluxes in daily bins for two energy bands (note the two different y-axes); VHE fluxes obtained from MAGIC and VERITAS above 0.2 TeV (in log scale); Fermi-LAT fluxes in 3-day bins in two energy bands; X-ray fluxes in 1-day bins from the all-sky monitors Swift-BAT and RXTE-ASM; X-ray fluxes from the pointing instruments Swift-XRT and RXTE-PCA; hardness ratio between the high and low-energy fluxes of Swift-XRT and between the two VHE bands of MAGIC (note the two different y-axes); optical R-band data from GASP-WEBT; radio data from Metsähovi, UMRAO, SMA, OVRO and VLBA; polarization degree and polarization angle observations in the optical from the Steward and Perkins observatories and in radio from VLBA.

The second row shows the combined VHE light curve provided by MAGIC and VERITAS on a logarithmic scale. Due to the additional VERITAS data, the gap in February 2010 is covered. These observations yielded the largest VHE flux observed to date in an AGN, reaching almost 10 C.U. for the daily average. Using a shorter binning, the peak flux values reached ∼15 C.U. above 200 GeV (Abeysekara et al. 2020).

The third row displays the emission in high-energy gamma rays (HE; E > 100 MeV) measured by Fermi-LAT in a three-day binning in the 0.3 GeV to 3 GeV and 3 GeV to 300 GeV band. During the large VHE flare in February, both bands show increased emission and reached their highest flux coincident with the highest VHE flux point. The other two VHE flares during January and March do not show increased activity in HE gamma rays. We additionally checked if daily binning reveals any activity on a shorter scale during these two flares that might be hidden with a 3-day binning but did not find significant variability.

In the fourth and fifth rows, the results of the X-ray are shown. A short outburst is detected in November 2009, around MJD 55145. The flare faints over a few days, and the source remains at a comparably lower activity level without any noticeable outbursts. Coinciding with the VHE flare in February, all wavebands show a sharp increase in flux, which slowly decreases in the following month. At the end of the campaign, the X-ray emission also entered a remarkably low state.

Using the two separate VHE and X-ray light curves, hardness ratios (i.e. the ratio of the flux in the higher-energy band to the flux in the lower-energy band) can be obtained and are shown in row six. Overall, the source showed a harder-when-brighter trend in VHE and X-ray, which has been observed frequently in the past (e.g. Acciari et al. 2021; MAGIC Collaboration 2021). This behavior is investigated in more detail in Appendix D.

Row seven and eight contain the UV and optical fluxes, respectively. Due to the closeness in frequency, the trends are quite similar. Both light curves reach their highest emission coincident to the January VHE flare. When the subsequent February flare starts, there is no clear visible flare in the UV or optical.

The ninth row shows the fluxes obtained from radio observations. The instruments cover a variety of frequencies, which makes flux comparisons challenging. Overall, the source showed a low activity and a low variability in all radio bands.

The second to last row reports the polarization degree, which shows some strong fluctuations throughout the campaign, going from values as low as 0.3% up to ∼11% around a mean of ∼5%. The corresponding electric vector polarization angles (EVPA) in the optical are shown in the last row. The data were adjusted for the intrinsic 180° ambiguity as reported in Carnerero et al. (2017). The EVPA shows a constant behavior at the beginning at around ∼150° followed by a wide rotation happening around MJD 55230 (early February). For the remaining campaign, it fluctuates around ∼350°. Since the EVPA rotates by around 200° during a time interval with very limited coverage (only one measurement), one cannot completely exclude that this rotation arises from an improper correction of the 180° ambiguity. If the rotation is real, there might be a possible connection with the large flare in February at higher energies.

3. VLBA observations of the jet evolution

Fig. 2 shows the total intensity VLBA images of Mrk 421 from May 2010 to August 2011. The images are shown from May onwards when the ejection of the first feature is observed. The source structure was modeled using a number of emission components (knots) with circular Gaussian brightness distributions. Based on their parameters, we identify knots across the epochs. The brightest knot located at the southern end of the jet is the VLBI core at 43 GHz, which we place at zero distance in Fig. 2. The core, designated as A0, is assumed to be a stationary physical structure of the jet. Another (quasi-)stationary feature is located at 0.38 ± 0.06 mas from the core, labeled A1. At the distance of Mrk 421, one can use the conversion scale 0.638 pc mas⁻¹² to convert the angular distance of 0.38 mas into a physical distance of 0.24 pc. During the observations, two new components emerged from the core, K1 and K2. The components were already reported in Jorstad et al. (2017) (K1 and K2 are referred to as B1 and B2, respectively). As reported there, the existence of B1 (K1 in this manuscript) as a real knot is not absolutely certain. All average parameters of the components (i.e. names of the knot, number of epochs at which the knot is detected, flux density, distance from the core, position angle with respect to the core, and the size of the knot) are listed in Table 1. The corresponding kinematic properties of K1 and K2 are listed in Table 2 (Phi refers to the speed direction). The newly emerging component K1 moves at a speed of 0.78 mas/yr, which corresponds to an apparent speed of 1.56 ± 0.45 c, making it just about a superluminal knot (note the large uncertainty). K2 moves much slower with a speed of 0.16 mas/yr corresponding to 0.32 ± 0.07 c.

Fig. 2. Total intensity VLBA images from May 2010 to August 2011. The black and yellow lines show the average positions of the stationary features A0 and A1. The yellow, red and blue circles show the fitted positions of A1, K1, and K2, respectively.

Fig. 3 shows the angular distance from the core of the components A1, K1, and K2, as a function of time. Since A1 is identified as a stationary feature, the distance versus time of A1 with respect to A0 is fit with a horizontal line, while separations of two new components from the core are approximated by linear fits with the best χ² value. This allows us to determine times of ejection (passage through the center of the core) of K1 and K2 by extrapolation of the motion of knots back to the core position (see Table 2).

Fig. 3. Angular distance of A1 (yellow), K1 (red), and K2 (blue) from the jet core (black). At the location of Mrk 421, the angular distance of one mas is equivalent to a physical distance of 0.638 pc. The dotted yellow line indicates a constant fit of the quasi-stationary component. The red and blue lines show linear fits to determine the time of ejection of the components. The uncertainties on the ejection times are given by the red and blue bands. The campaign of this work is given by the light gray band with the three flares marked with gray lines.

K1 is ejected at MJD 55112±88 and K2 at MJD 55400±157. Both ejection times fall well within the main period of this work, shown in light gray, from MJD 55140 to 55380 within their one-sigma band. Because of the large uncertainties, however, it is impossible to establish a connection with a particular event, such as one of the three VHE flares, marked with solid gray lines.

Nevertheless, it is remarkable that in the period Mrk 421 exhibited the most violent behavior with three large VHE flares, including the brightest flare ever, the ejection of two new components was seen right after. Despite VLBA components having been identified for a few TeV-emitting blazars (mostly FSRQs), a potential association of gamma-ray flares with the ejection and propagation of knots is rarely observed in HSPs (Weaver et al. 2022). To our knowledge, this is the only time the ejection of new components in the jet of Mrk 421 was observed. Previous works only found stationary or already present subluminal knots (e.g. Piner & Edwards 2005; Lico et al. 2012; Richards et al. 2013; Blasi et al. 2013).

4. Variability

4.1. Fractional variability

In order to estimate the degree of variability in each energy band, we use the fractional variability F_var as it is described in Vaughan et al. (2003). The estimation of the uncertainties follows the description in Aleksić et al. (2015b), which uses the approach from Poutanen et al. (2008). Fig. 4 shows the resulting fractional variability using the full light curves in open markers. The full markers show the results for quasi-simultaneous data. We define quasi-simultaneity as the temporal agreement with VHE data within 6 h for X-ray and UV data. The choice of 6 h is a compromise between having sufficient data and not too much loss of temporal coincidence. We use a window of 1 day for optical data and 3 days for radio and Fermi-LAT data.

Fig. 4. Fractional variability Fvar as a function of the frequency for the light curves shown in Fig. 1. Open markers show the results using the whole campaign for each light curve. Full markers only include quasi-simultaneous data to the VHE data (quasi-simultaneity is defined as temporal agreement with VHE data within 6 h for X-ray and UV data, within 1 day for optical data, and within 3 days for radio and Fermi-LAT data).

In both cases, a pronounced two-peak structure is visible, indicating the highest variability in the Xray and VHE. Similar behavior has been observed multiple times for Mrk 421 (e.g. Aleksić et al. 2015c; Baloković et al. 2016; MAGIC Collaboration 2021). As mentioned in Sect. 1, the typical SED from a blazar such as Mrk 421 shows a double bump structure (Abdo et al. 2011), where the first bump is electron synchrotron emission. The falling edge of that bump, covered by the X-rays, is emitted by the most freshly accelerated and energetic electrons, which also have the shortest cooling times, meaning high variability. At energies below the X-rays, the electrons responsible for the synchrotron emission have lower energies and hence longer cooling times, resulting in a strong decrease in variability. In the context of a SSC model, the same electrons produce the second bump of the SED in the gamma rays via inverse Compton emission. A similar variability pattern is therefore also expected in the variability of the gamma-ray peak. Fig. 4 shows the same degree of variability at TeV energies as in keV energies. This result already suggests the existence of a correlation between the two wavebands, which is investigated in more detail in Sect. 5.

4.2. Power spectral density

Additionally, we investigate the variability of Mrk 421 with the use of the power spectral density (PSD). The PSD quantifies the amplitude of the variability as a function of the timescale of the variations. It is based on the discrete Fourier transform of a light curve. The typical shape of the PSD for blazars follows a simple power law P_ν ∝ ν^−a with the spectral index a ranging between 1 and 2 (e.g. Uttley et al. 2002; Chatterjee et al. 2008; Abdo et al. 2010b). A falling power law indicates a large variability at long timescales (small frequencies), and the power of this variability decreases as one considers shorter and shorter timescales (high frequencies).

Our method of estimating the PSD indices is described in Appendix B, and the results obtained are reported in Table 3. We find values for the PSD indices that are well compatible within their uncertainties with the ones derived in Aleksić et al. (2015c). Due to the larger uncertainties of the light curve, the PSD indices for the Fermi-LAT bands are not well constrained and show noticeably smaller values. We did not find evidence of a spectral break of the PSD index in any energy band. In addition, we used the Lomb-Scargle method (Lomb 1976; Scargle 1982) to search for signs of periodicity but found no evidence in our data set.

5. Correlation studies

5.1. VHE gamma rays versus X-rays

The VHE gamma rays and the X-rays are the two energy bands where Mrk 421 shows the highest flux variations (as displayed in Fig. 4), which can occur on timescales as short as hours (see Appendix C). Because of that, the observing campaign was organized with a special focus on maximizing the simultaneity between the VHE gamma rays and X-ray measurements. Hence, the collected dataset allows us to investigate the correlations with a large number of VHE gamma rays and X-ray measurements taken within a time window of only a few hours. Previous works investigated the source behavior either during a low state of activity (Baloković et al. 2016), a typical state of activity (MAGIC Collaboration 2021) or during a flaring activity (Acciari et al. 2020). The large variability shown by Mrk 421 during the MWL campaign in 2010 allows us to probe the correlation behavior across all states of emission within a single campaign lasting about half a year. Since we did not find a time delay between the two bands, we evaluate the data at zero time lag and with a flux-flux plot within a time window of only 6 hours. Fig. 5 shows the fluxes (in decimal logarithmic scale) of all three VHE bands versus the fluxes of the two Swift-XRT bands and the one provided by Swift-BAT.

Fig. 5. VHE flux versus X-ray flux obtained by MAGIC/VERITAS and Swift-XRT/BAT. Only pairs of observations within 6 hours are considered. If more than one Swift observation falls within that window, the weighted mean is computed. VHE fluxes are in the > 1 TeV band (top panels), in the 0.2−1 TeV band (middle panels), and in the > 0.2 TeV band (bottom panels). Swift-XRT fluxes are computed in the 0.3−2 keV (left panels) and 2−10 keV bands (middle panels). Swift-BAT provides the flux in the 15−50 keV band (right panels). The top left corner of each panel shows the Pearson coefficient of the flux pairs, with the significance of the correlations given in parentheses. The gray dashed and dotted lines depict a fit with slope fixed to 0.5 and 2, respectively, and the black line is the best-fit line to the data, with the slope quoted in the legend at the bottom right of each panel.

In order to evaluate the degree of correlation, we use the Pearson coefficient using the pairs of flux data shown in Fig. 5, but without applying the decimal logarithm. We quantify the significance of the correlation using dedicated Monte-Carlo simulations. The methods are described in Appendix E. We find strong correlations between all VHE and Swift-XRT energy bands. In all six panels, we find Pearson coefficients equal or above 0.75 with significances ranging from slightly above 3 up to 4σ. The correlation between the VHE and the Swift-BAT band is much weaker, reaching coefficients around 0.4−0.5 and significances around 2σ. The significance of the correlation may be reduced by the large flux uncertainties of the BAT measurements (in comparison to the small uncertainties in the flux measurements from XRT).

We investigate the correlation slopes by fitting lines to the decimal logarithm of the data shown in the scatter plots. Since we did not find a significant change in the correlation slope over time by fitting individual ∼1 month periods, we show a single line fit for the whole campaign. It is noteworthy that in the subplots between the > 0.2 TeV flux with both X-ray bands, the observations of the highest and lowest X-ray flux lie well on the linear fit. This indicates a linear trend reaching over a full order of magnitude of emission and on the timescale of multiple months. The steepest slope (1.5 ± 0.2) is obtained for the highest VHE gamma-ray band versus the lowest X-ray band, while the flattest slope (0.3 ± 0.1) is obtained for the lowest VHE gamma-ray band versus highest X-ray band. Overall, the slope is increasing with rising VHE energy, showing a greater scaling of the > 1 TeV flux with rising X-ray fluxes. Additionally, the slope decreases with rising X-ray energy. The same trends have been found in MAGIC Collaboration (2021) as well, but the absolute slope values were found to be greater, reaching a cubic relation in some cases.

5.2. VHE gamma rays versus UV

The Swift-UVOT instrument has the same data coverage in the UV as Swift-XRT in X-rays. This allows us to follow the same approach as in the previous section to investigate a possible correlation between the VHE and UV band pairs. Fig. 6 shows the decimal logarithm of the fluxes of all three VHE bands versus the flux of the W1 filter (all three UV filters provide almost identical results, and we have selected W1 as a representative for all correlations). The values for the Pearson coefficient range from around 0.45 to 0.65 for all three panels, indicating a much weaker correlation. However, the corresponding significances of these correlations are all around 2σ, and hence, this result should be considered as a hint of correlation. If the correlation were real, this would be the first time, to our knowledge, that a positive correlation between VHE gamma rays and the UV band is found.

Fig. 6. VHE flux versus UV flux obtained by MAGIC/VERITAS and Swift-UVOT. Only pairs of observations within 6 hours are considered. If more than one Swift observation falls within that window, the weighted mean is computed. VHE fluxes are in the > 1 TeV band (top panels), in the 0.2−1 TeV band (middle panels) and in the > 0.2 TeV band (bottom panels). The Swift-UVOT fluxes are taken with the W1 filter. The top left corner of each panel shows the Pearson coefficient of the flux pairs, with the significance of the correlations given in parentheses. The gray dashed and dotted lines depict a fit with slope fixed to 0.5 and 2, respectively, and the black line is the best-fit line to the data, with the slope quoted in the legend at the bottom right of each panel.

For the linear fits, we obtain compatible values of around m = 2 for all three panels. Due to the higher dispersion of the data, the uncertainties of the fit are rather large. The higher slope indicates a steeper correlation of the VHE with UV compared to X-rays. Since the degree of variability is much higher in the VHE than the UV (see Fig. 4), the flux varies substantially more, and a steeper slope is expected in the case of a correlation.

5.3. VHE versus HE gamma rays

Next, we investigate the correlation between the two energy bands in gamma rays using the combined light curve above 0.2 TeV from MAGIC and VERITAS. We correlate the bands using the DCF without rebinning the LCs. We use a 3-day binned time lag in the range of −40 to +40 days. The significance of the resulting DCF is estimated with simulations similar to the approach in MAGIC Collaboration (2021): We first simulate a set of 10 000 uncorrelated light curves for a pair of energy bands as before. We then compute the DCF of the 10 000 simulated light curve pairs. The 1σ, 2σ and 3σ confidence bands are derived from the distribution of the simulated DCF values in each time lag bin.

The result for the DCFs between the > 0.2 TeV band with the 3−300 GeV band is shown in Fig. F.1a. The figure shows the DCF obtained from the data in dark blue. The DCF peaks at a time lag of 2 ± 2 days, where it crosses the 3σ line. To estimate the uncertainty of the time lag, we follow the method outlined in Peterson et al. (1998).

We note that the marginally significant (a little over 3σ) peak at zero time lag is primarily driven by the big flare in February 2010 because the highest flux points in the VHE and HE light curves occur during this flare. If the highest flux in the HE light curve (that relates to a 3-day time interval from MJD 55242.0 to MJD 55245.0 or 2010-02-15 00:00 to 2010-02-18 00:00), together with the corresponding VHE fluxes in this 3-day time interval, are removed from the DCF study, the correlation at zero time lag vanishes, as shown in Fig. F.1b. A similar behavior is seen in the correlation between the > 0.2 TeV band and the 0.3−3 GeV band, as it is also shown in Figs. F.1c and F.1d. But this time, the significance of the peak at zero time lag is just over the 2σ line, even including the 2010 February flare.

5.4. HE gamma rays versus X-rays

For the correlation between HE gamma rays from Fermi-LAT and X-rays from Swift-XRT, we again use the DCF. We correlate each of the two HE gamma-ray bands with each of the X-ray bands from Swift-XRT. Here, we only display the strongest correlation as an example, which is found between the 3−300 GeV and the 0.3−2 keV band, shown in Fig. F.2a. A correlation peak lying above 3σ is found at a time lag of −1 ± 2 days. Similar to the VHE versus HE gamma-ray case, the sharp peak at zero time lag is likely dominated by the outstanding flaring activity on February 17. Fig. F.2b shows the same correlation but with the flare removed from both the light curves taken by Fermi-LAT, as well as the simultaneous flux in the X-ray light curve. The high peak at exactly zero time lag vanishes, but a broad peak remains.

The other three combinations of energy bands (0.3−3 GeV versus 0.2−3 keV, 3−300 GeV versus 2−10 keV, 0.3−3 GeV versus 2−10 keV) are also shown in Fig. F.2. In all combinations, one can see a marginally significant peak at a time lag of about zero that disappears when removing the 3-day LAT flux bin (and corresponding X-ray fluxes) from the outstanding flaring activity around 2010 February 17. The period in early 2010, including the flare, marks the first time such a correlation was observed, as previously reported in Abeysekara et al. (2020).

5.5. HE gamma rays versus UV

Fig. F.3a shows the DCF between the 3−300 GeV and the UV bands. It exhibits a broad peak crossing the 3σ line, and it reaches its highest value at a time lag of −7 ± 6 days. Between the 0.3−3 GeV band and the UV in Fig. F.3c, the highest value at −1 ± 6 days also reaches a significance level above 3σ. The positive correlation between both HE gamma-ray regions and the UV band is a novel behavior for Mrk 421 with no previous reports to our knowledge. Acciari et al. (2021) found a positive correlation compatible with a time lag of zero between the 0.3−300 GeV band, and the R-band, which frequency is close to the UV. The correlation was found by investigating the time period 2007–2016 with a 15-day binning resulting in a lower temporal resolution than the one used here.

Contrary to the previous cases, the correlation between HE gamma-rays and the UV increases if the flare in February is excluded. Figs. F.3b and F.3d show the results of correlating the HE gamma-ray light curves with the UV light curve without the flare. In the case of 3−300 GeV, the peak is clearly enlarged, going well above 3σ, indicating a stronger correlation than before. For 0.3−3 GeV, the improvement is only marginal. The peaks remain at the same time lag as before.

5.6. HE gamma rays versus R-band

We additionally test whether the correlation between the HE gamma rays and the optical as previously reported in Acciari et al. (2021) is also present in our data set. Fig. F.4 shows the DCF for the R-band and the two HE bands, which does not lead to any significantly correlated behavior, apart from a marginally significant (∼3σ) feature at about −35 days.

Removing the 3-day interval that contains the large VHE and X-ray flare in February 2010 does not change substantially the shape and the significance of the DCF results for these two bands. Only a marginal increase to 2.5σ is seen for the 3−300 GeV versus optical (see Fig. F.4b).

5.7. Other energy bands

We also investigated the correlation behavior between all other wavebands not mentioned in the previous sections. We found no correlation (positive or negative) between any bands in the X-ray with the UV or the optical. Previous works have found anti-correlated behavior at times (e.g. MAGIC Collaboration 2021; Abe et al. 2024), which was not present in our data set. Naturally, there is a strong correlation between the UV and optical, as the energy is close to each other. Since we established a hint of a correlation between VHE and UV, we also checked if a correlation is present between the VHE and the optical. Applying the same method, we find no correlations at all between VHE and optical. We also found no correlation of any sort between radio and other energy bands. Due to the limited duration of the campaign (about seven months), the overall number of radio observations is also relatively small, in comparison to studies that used multi-year datasets (see e.g. Carnerero et al. 2017; Acciari et al. 2021). In addition, the synchrotron self-absorption of radio radiation might reduce any detectable correlation.

6. Discussion

6.1. First detection of ejection of radio knots contemporaneous to a flare of Mrk 421

VLBA observations reveal the ejection of two features in the jet that are close in time to the observed VHE gamma-ray flares (see Sect. 3). While the association of the ejection and propagation of bright knots in the jet with gamma-ray flares has been observed repeatedly in blazars, the vast majority of these detections were reported in FSRQs (like PKS 1510–089) or LSP and ISP BL Lacertae (see Jorstad et al. 2001; Marscher et al. 2008; MAGIC Collaboration 2018). The association has rarely been observed in HSPs, and this is the first time for Mrk 421 (see e.g. the sample in Weaver et al. 2022), hence making the dataset presented in this paper unique on its own. However, the large uncertainties of the ejection time prevent us from significantly correlating the appearance of the radio features with the specific VHE gamma-ray flares during this campaign. We also note that previous works on Mrk 421 reported a possible connection between GeV gamma-ray and radio flares detected with single-dish telescopes (Hovatta et al. 2015). However, this large GeV flaring activity, which was measured with Fermi-LAT in the year 2012, and the exceptionally large radio flare, which was (mostly) measured with OVRO also in 2012, was not associated with the ejection of features in the VLBA images (Richards et al. 2013). Hence, they can have a different nature, and are not comparable to the ejection of features in the VLBA images during the 2010 campaign that is reported in this manuscript.

Assuming the flare in February is indeed correlated with the ejection of the radio features, the VHE intranight variability implies a small, compact zone responsible for the emission based on causality arguments. As discussed in Abeysekara et al. (2020), following Dondi & Ghisellini (1995), the intranight variability down to 22 min at VHE in February 2010 allows us to derive a lower limit on the Doppler factor of δ_min ≳ 33 based on opacity arguments.

Jorstad et al. (2017) estimate the average Doppler factor derived from VLBA observations to be around 24, hence below the lower limit mentioned above using the VHE observations. This regularly observed discrepancy between the observed jet kinematics in the radio and the Doppler factors derived from the fast variability at high energies is referred to as the Doppler crisis. This could potentially be resolved by the presence of structured jets, in which different regions show different Lorentz factors, such as a fast spine responsible for the gamma-ray emission and a slower sheath layer responsible for the radio emission (Ghisellini et al. 2005). Recent works also proposed a scenario in which the discrepancy is reconciled by assuming, on the one hand, a large viewing angle relative to the jet explaining the slow-moving knot in radio, while on the other hand, the gamma-ray flare originates from a magnetic reconnection event within a misaligned layer that effectively generates large Doppler factors (Jormanainen et al. 2023). Finally, jet deceleration occurring between the time of the gamma-ray flare and the detection of the radio feature may partly solve this crisis (Georganopoulos & Kazanas 2003).

The FWHM of the radio knots reported in Sect. 3 is in the order of 0.15 mas (both for K1 and K2). At the distance of Mrk 421, this corresponds to a linear size of roughly 0.1 pc ≈ 3 × 10¹⁷ cm (0.638 pc/mas). Assuming a Doppler factor of 33 and a variability timescale of 22 minutes, the upper limit to the size of the emitting region during the VHE flare is in the order of 10¹⁵ cm (Abeysekara et al. 2020). If one again postulates that the radio knot corresponds to the same region that produced the VHE flare in February 2010, this implies an expansion velocity in the comoving frame of the emitting region in the order of β_exp = 10⁻²c (using a time difference between the gamma-ray flare and the radio detection of roughly 300 days, as one can infer from Fig. 3). It is interesting to note that such expansion velocity is well in line with the estimates of Tramacere et al. (2022), which are obtained when trying to explain the delayed gamma-ray and radio response in Mrk 421 by an adiabatic expansion of the blob. The adiabatic expansion of the blob is expected to induce a decrease in the electron density on timescales of $T_{\mathrm{ad}}(t)=\frac{R(t)}{3{\beta }_{\exp }c}$ (Gould 1975). At the time of the VHE flares R ∼ 10¹⁵ cm, hence T_ad ≈ 3 × 10² h, corresponding to ≈10 h in the observers frame, being much longer than the minimum variability timescale noted at VHE (22 min). Hence, the VHE flux variability at the shortest timescales is likely dominantly driven by acceleration and cooling mechanisms and negligibly affected by the expansion of the blob.

6.2. Multiband correlations as a probe of SSC in the jet

We find a close correlation between the emission in the VHE gamma rays at zero time delay and the X-rays over an order of magnitude in flux. We also find a comparable level of variability in both energy bands, similar to the results in MAGIC Collaboration (2021) and Acciari et al. (2021). The tight correlation and similar variability behavior suggest a cospatial origin, which is in good agreement with the SSC scenario (e.g. Maraschi et al. 1992). The strongest correlation is found between the 0.3−2 keV and > 200 GeV bands reaching a level of 4 σ. Assuming standard values for the jet parameters (B = 0.1 G and δ = 35 see e.g. Aleksić et al. 2015a), one derives that electrons with Lorentz factor of around 10⁵ (source reference frame) are required to emit ∼keV photons (in observer’s frame). Using Eq. (14) of Tavecchio et al. (1998), which takes into account the Klein-Nishina effects, one expects that electrons with such energy would emit ∼0.5 TeV photons via IC scattering of ∼keV photons. This is well in agreement with the correlation trends.

We also report marginally significant correlations between the UV and HE gamma rays for the first time in Mrk 421, with both bands showing a similar degree of variability. The pattern is again consistent with the SSC model, but in this case caused by electrons of lower energies populating a larger emission zone. Following an analogous approach as in the previous case, we derive estimates of the necessary electron Lorentz factors. Synchrotron photons in the UV W1 band are produced by electrons with a Lorentz factor of ∼10⁴. These electrons produce IC emission at a few tens of GeV, falling well within the 3−300 GeV band, for which the highest correlation is observed.

The correlation strength decreases to a non-significant level when going from the UV to the optical band. Even though the UV and R-band are close together in frequency, these results imply that the underlying particle population responsible for the HE gamma rays are dominantly radiating in the UV rather than in the R-band (assuming a leptonic model). Carnerero et al. (2017) and Acciari et al. (2021) reported significant positive correlations between the HE gamma rays and the R-band. These two studies considered much longer periods (∼8 years) where, in relative terms, Mrk 421 showed overall lower activity. The slightly different correlation behavior might be explained by the different overall magnitude and timescales of the flux variability considered, as well as the underlying mechanisms driving the emission variability (such as the evolution of the emitting region environment, change in the acceleration and cooling efficiencies etc...)

6.3. Evidence of multiple emission zones

We find no correlation between the X-ray and the UV. The most simple explanation of this result is with the existence of two separate particle populations, of which a compact zone with higher particle energies is responsible for the keV (and TeV) emission, and a larger and more extended zone dominates the emission in the eV (and possibly GeV) band. If the compact zone is embedded into the larger zone, this two-zone scenario would also be consistent with the marginal significant (∼2.0 − 2.3σ) correlation between VHE gamma rays and the UV. Should the UV emission increase due to a change in the underlying particle population, synchrotron photons from the extended zone would enter the compact zone and provide additional seed photons, which would then be IC scattered and hence produce an increase in the observed VHE gamma-ray flux.

The preference for a multiple-zone scenario was also noted in earlier works that included X-ray data from this campaign (see e.g. Kapanadze et al. 2018). It is also further motivated by new findings of the Imaging X-ray Polarimetry Explorer (IXPE), which suggest an energy dependency of the polarization degree from radio to the X-ray. The energy dependency points towards an energy-stratified jet creating emission zones of different extent (Di Gesu et al. 2022, 2023; Abe et al. 2024). The proposed scenario of two separate emission zones is a simplified implementation of an energy-stratified jet. Similar polarization results have also been found for the HSP Mrk 501 (Liodakis et al. 2022). The outlined two-zone leptonic model with interaction between the zones has been successfully used to fit the SEDs during the IXPE observations (MAGIC Collaboration 2024).

6.4. Peculiarities of the exceptional flare in February 2010

We find a positive correlation between the VHE and HE gamma rays over the whole campaign. A similar correlation was first reported for a dataset taken in 2015–2016, during which Mrk 421 was in a historically low-activity state (Acciari et al. 2021). However, the correlation in our data set is non-significant when removing the 3-day interval related to the large February 2010 flare, and hence these two bands were directly connected only during this exceptional flaring event. In other words, the VHE versus HE correlation is not representative of the behavior of Mrk 421 during the 8-month long 2010 campaign (which covers different activity levels).

The same phenomenon is observed when considering the HE gamma rays and the X-rays, for which the correlation also becomes non-significant when the 3-day time interval of the February 2010 flare is removed. Including the flare, we find the strongest correlation for the highest-energy band of Fermi-LAT (3−300 GeV) and the lowest-energy band of Swift-XRT (0.3−2 keV). While our study shows that the X-ray and VHE radiation share a common emission region and the UV and HE bands may originate from a different region, the flare could be driven by a single compact zone whose activity is so high that it not only dominates in the VHE radiation, but also contributes significantly to the emission of HE gamma rays on these days. The X-ray and HE gamma-ray bands might then be close to the peak of the SED bumps caused by synchrotron and IC emission of the same particles. Matching this scenario, the 0.3−3 GeV band, sitting at the rising edge of the second bump, and the 2−10 keV band, sitting at the falling edge of the synchrotron bump, show the weakest correlation.

Data availability

Data shown in Fig. 1 are available available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr ( 130.79.128.5 ) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/694/A195

Acknowledgments

Author contribution: J. Abhir: MAGIC analysis cross-check; A. Arbet-Engels: variability and correlation analysis, discussion and interpretation, paper drafting; D. Paneque: coordination of MWL observations and coordination of the MWL data reduction, RXTE-PCA analysis, discussion and interpretation, paper drafting; F. Schmuckermaier: project management, MAGIC and Fermi-LAT data analysis, variability and correlation analysis, discussion and 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. The MAGIC collaboration would like to thank the Instituto de Astrofísica de Canarias for the excellent working conditions at the Observatorio del Roque de los Muchachos in La Palma. The financial support of the German BMBF, MPG and HGF; the Italian INFN and INAF; the Swiss National Fund SNF; the grants PID2019-104114RB-C31, PID2019-104114RB-C32, PID2019-104114RB-C33, PID2019-105510GB-C31, PID2019-107847RB-C41, PID2019-107847RB-C42, PID2019-107847RB-C44, PID2019-107988GB-C22, 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. A.A.E. and D.P. acknowledge support from the Deutsche Forschungsgemeinschaft (DFG; German Research Foundation) under Germany’s Excellence Strategy EXC-2094-390783311. The research at Boston University was supported in part by several NASA Fermi Guest Investigator grants, the latest of which are 80NSSC23K1507 and 80NSSC23K1508. This study was based in part on observations conducted using the 1.8m Perkins Telescope Observatory (PTO) in Arizona, which is owned and operated by Boston University. The VLBA is an instrument of the National Radio Astronomy Observatory. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated by Associated Universities, Inc. Research at UMRAO was supported by a series of grants from the NSF (most recently AST-0607523) and from NASA (including Fermi G.I. awards NNX09AU16G, NNX10AP16G, NNX11AO13G, and NNX13AP18G). Funds for the operation of UMRAO were provided by the University of Michigan. 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.

References

Appendix A: Observations and data processing

This section contains a detailed description of all used instruments and their corresponding data analysis. Since the same campaign is exploited in parts in Aleksić et al. (2015a), some instrument sections closely follow the descriptions given in there.

A.1. VHE gamma rays

The MAGIC telescopes are two Imaging Atmospheric Cherenkov Telescopes (IACTs), MAGIC I and MAGIC II, each with a diameter of 17 m, located at Observatorio del Roque de los Muchachos (ORM, 28.762°N 17.890°W, 2200 m above sea level) on the Canary Island of La Palma. With the start of stereoscopic observations in 2009 and substantial hardware upgrades completed in 2012, MAGIC is capable of detecting gamma rays with energies from about 30 GeV up to ≳100 TeV (Aleksić et al. 2016; MAGIC Collaboration 2020).

This work covers the time period from November 11, 2009, (MJD 55149) until June 16, 2010, (MJD 55363). In total, the MAGIC telescopes observed Mrk 421 for 62.4 h in the zenith angle range between 5° and 50°. The data are analyzed using the MAGIC Analysis and Reconstruction Software, MARS (Zanin et al. 2013; Aleksić et al. 2016). The final data are selected based on quality criteria to exclude periods with unfavorable weather conditions or too bright moon. The data were taken under low moonlight conditions in order to limit contamination from night sky background light (Ahnen et al. 2017). After applying quality cuts, 50.1 h of data remained.

This work focuses on the MWL variability and correlations, and hence, we use light curves that report the VHE gamma-ray emission of Mrk 421. Due to the high gamma-ray brightness of Mrk 421, we can construct two separate light curves in the VHE band covering two energy ranges: 0.2-1 TeV and > 1 TeV. The analysis was performed for the whole campaign including the already published periods in Aleksić et al. (2015a) and Abeysekara et al. (2020). This ensures that the low-level analysis remains consistent throughout and the results are well compatible. In addition to the data from MAGIC, we also used VERITAS light curves with VHE gamma-ray fluxes above 0.2 TeV published in Aleksić et al. (2015a) and Abeysekara et al. (2020). In order to create a joint light curve, a third light curve was constructed from the MAGIC data for this energy range.

A.2. HE gamma rays

The Fermi Gamma-ray Space Telescope satellite carries the LAT detector on board. It is a pair-conversion telescope surveying the gamma-ray sky in the 20 MeV to > 300 GeV energy range (Atwood et al. 2009; Ackermann et al. 2012). For the construction of light curves, we perform a binned-likelihood analysis using tools from the FERMITOOLS software³ v2.2.0. We use the instrument response function P8R3_SOURCE_V3 and the diffuse background models⁴ gll_iem_v07 and iso_P8R3_SOURCE_V3_v1.

We create two light curves in the range from 0.3 GeV to 3 GeV and 3 GeV to 300 GeV by selecting Source class events in a circular region of interest (ROI) with a radius of 20° around Mrk 421 in the respective energy band. All events with a zenith angle above 90° are discarded to reduce the contribution from limb gamma rays. The analysis threshold energy of 0.3 GeV was chosen over the more usual 0.1 GeV in order to make use of the improved angular resolution of Fermi-LAT at higher energies. A higher energy threshold additionally reduces background contamination, which leads to an overall improvement of the signal-to-noise ratio for hard sources (photon index < 2) such as Mrk 421. For the source model, we include all sources from the fourth Fermi-LAT source catalog Data Release 3 (4FGL-DR3; Abdollahi et al. 2020; Ajello et al. 2022) that are found within the ROI plus an additional 5°. We fit the obtained model to our data covering the time period from November 2, 2009, (MJD 55141) to July 2, 2010, (MJD 55379). The initial fit results are used to remove weak components from the model (counts < 1 or TS < 3⁵). After the first optimization, each time bin is fit with the model. During the fitting procedure, the normalization of bright sources (TS > 10), sources close to the ROI center (< 3°), the diffuse background, and the spectral parameters of Mrk 421 itself are allowed to vary. We produce the light curves with a 3-day binning. In all time bins, the source is detected with TS > 25 (i.e., > 5σ).

A.3. X-rays

We scheduled observations with the X-Ray Telescope (XRT; Burrows et al. 2005) on board the Neil Gehrels Swift Observatory (Swift) throughout the full campaign to achieve the best possible coincidence with VHE observations. All Swift-XRT observations were carried out both in Windowed Timing (WT) and Photon Counting (PC) readout modes. The data were then processed using 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.30.1). The data were reprocessed with the xrtpipeline script and using calibration files from Swift-XRT CALDB (version 20210915). The X-ray spectrum was constructed from the calibrated and cleaned event file for each pointing. The events were selected within a radius of 20 pixels (∼47 arcsec) in both WT and PC modes. The background was extracted from a nearby circular region with a radius of 40 pixels. The ancillary response files were generated with the xrtmkarf task applying the corrections for PSF losses and CCD defects using the cumulative exposure map. The 0.3 − 10 keV source spectra were binned using the grppha task by requiring at least 20 counts per energy bin. We used XSPEC with both a power-law and log-parabola model (with a pivot energy fixed at 1 keV and an added photoelectric absorption component assuming an equivalent hydrogen column density fixed to the Galactic value along the line of sight). In the overall majority of the observations, the preference for a log-parabola model is statistically significant (i.e., > 5σ). The intrinsic fluxes were extracted in the 0.3-2 keV, and 2-10 keV energy bands.

Data from the Swift Burst Alert Telescope (BAT) are publicly accessible online⁷. In this work, the daily light curve is used from 15-50 keV.

The Rossi X-ray Timing Explorer (RXTE; Bradt et al. 1993) satellite performed almost daily pointing observations of Mrk 421 during the time period from December 12, 2009, (MJD 55177) to April 29, 2010, (MJD 55315). The data were analyzed using FTOOLS v6.9 following the settings and procedures recommended by the NASA RXTE Guest Observer Facility⁸. We produce a light curve from 2-10 keV. For more details on the analysis settings see the instrument description given in Aleksić et al. (2015a).

Data from the RXTE All-Sky Monitor (ASM) are publicly accessible online⁹. In this work, the one-day average light curve in the energy range from 2-12 keV is used.

A.4. Ultraviolet

Swift also provides coverage in the ultraviolet (UV) band from the UV/Optical Telescope (UVOT, Roming et al. 2005). We selected observations in the W1 (251 nm), M2 (217 nm), and W2 (188 nm) filters by applying standard quality checks to all observations in the chosen time interval, excluding those with unstable attitudes or affected by contamination from a nearby starlight (51 UMa). For each individual observation, we performed photometry over the total exposures in each filter to extract flux values. The same apertures for source counts (the standard with 5 arcsec radius) and background estimation (mostly three-four circles of ∼16 arcsec radii off the source) were applied to all. We performed the photometry extraction with the official software included in the HEAsoft 6.23 package, from HEASARC, and then applied the official calibrations (Breeveld et al. 2011) included in the CALDB release (20201026). As a last step, the source fluxes were dereddened according to a mean interstellar extinction curve (Fitzpatrick 1999) and the mean Galactic E(B − V) value of 0.0123 mag (Schlegel et al. 1998; Schlafly & Finkbeiner 2011).

A.5. Optical

Observations in the R-band were performed within the GLAST-AGILE Support Program (GASP) of the Whole Earth Blazar Telescope (WEBT; e.g. Villata et al. 2008, 2009), including multiple optical telescopes around the globe. We use the data published in Carnerero et al. (2017), for which the contribution of the host galaxy has been subtracted.

Optical polarization measurements are included from the Lowell (Perkins), Crimean, Calar Alto, and Steward observatories, also taken from Carnerero et al. (2017). The contribution of the host galaxy is also taken into account for the polarization data.

A.6. Radio

The radio data were taken with the 14 m Metsähovi Radio Observatory at 37 GHz, the 40 m Owens Valley Radio Observatory (OVRO) telescope at 15 GHz, and the 26 m University of Michigan Radio Astronomy Observatory (UMRAO) at 14.5 GHz. Observations at 225 GHz (1.3 mm) were performed at the Submillimeter Array (SMA) near the summit of Mauna Kea (Hawaii). For more details on the instruments and observations see the description given in Aleksić et al. (2015a).

In addition, the Very Long Baseline Array (VLBA) performed observations at 43 GHz. The data were provided by the Boston University blazar group¹⁰ and are part of their monitoring program of gamma-ray blazars. We obtained the images of the parsec-scale jet using the Astronomical Image Processing System (AIPS) and Difmap software packages (Jorstad et al. 2017).

Appendix B: Fitting the power spectral density

The PSD cannot be directly estimated from real data sets. Real observations suffer from unevenly sampled data with large gaps in the coverage. The limited time coverage as well as the temporal binning of the data can cause a transfer of variability power from lower to higher frequencies and vice versa. These distortions need to be accounted for obtaining an accurate estimate of the true PSD.

We estimate the true PSD index using a simulation-based forward-folding procedure, using the code described in Arbet Engels (2021), MAGIC Collaboration (2021). We cover the parameter space from 0.3 to 2.5 using steps of 0.1. For each assumed spectral index, we simulate a set of 3000 light curves using the assumption of a power-law-shaped PSD. We use the underlying probability distribution described in Emmanoulopoulos et al. (2013) to match the observed flux distribution as well as possible. We simulate light curves 10 times longer than the real light curves and apply the same temporal binning to account for the previously mentioned leakage effects. In order to compare the simulated with the real light curves, we use the multiple fractions variance function (MFVF) as a proxy for the PSD (Kastendieck et al. 2011). The MFVF is computed by splitting a given light curve in the temporal middle and taking the flux variance of both halves. The two halves are then split again into two, and the variances are estimated. This is repeated until a minimum time scale of one day is reached (which is the minimum temporal resolution of the light curves). The resulting variances as a function of time scale give a robust alternative representation of the PSD. The MFVF is computed for all simulated light curves, which results in a probability density function p(a, f_i) for each out of N frequency bins. It gives the probability of measuring a MFVF in the frequency bin f_i, assuming the underlying PSD with the index a. The best-fit value for the spectral index of the real light curve can then be estimated with a maximum likelihood estimation by finding the index a maximizing

$$\begin{array}{c}\hfill \mathcal{L}(a)={\displaystyle \sum _{i=0}^N}\ln (p(a,f_i)).\end{array}$$(B.1)

Fig. B.1 shows an example likelihood profile for the 2-10 keV band. The resulting values for a selected set of light curves are given in Tab. 3.

Fig. B.1. Likelihood profile ℒ(a) (see Eq. B.1) for the 2-10 keV band with its best fit value of 1.4 denoted by a vertical blue dotted line.

Following Arbet Engels (2021), MAGIC Collaboration (2021) we estimated the uncertainties of the indices by simulating 100 light curves using the best-fit index for each light curve and then refit all 100 created light curves with the same method. This creates a histogram of the resulting indices, where the uncertainties are given by the 68% (1 σ) containment region. As an example, Fig. B.2 shows the resulting histogram for the 2-10 keV band. Since the histograms can show a skewed distribution, we separately provide upper and lower uncertainties.

Fig. B.2. Histogram of the best-fit indices derived from simulations using as input a simulated light curve that has a known PSD index a = 1.4, in agreement with the real data for the 2-10 keV band.

Appendix C: Intranight VHE variability

Fig. C.1 shows the intranight light curves in the 0.2-1 TeV band of four selected nights in 10 min bins. The upper panel shows the January 14, 2010, (MJD 55210) light curve, where the highest fluxes are recorded by MAGIC during the flare in January. We fit the data with a constant model in order to test the hypothesis of a non-variable emission. The source shows no significant variability but a stable and high emission throughout the full exposure of around 3 hours.

In the second panel, the following night of January 15 (MJD 55211) is shown. The hypothesis of a constant emission is rejected at 3.7 σ, indicating a significant variability on a time-scale of well below 1 hour in VHE gamma rays. Missing data points in between were caused by technical interruptions.

During the flare in January, the VHE flux increases again around January 20 (MJD 55216). We see no significant variability but again a strong but constant emission throughout the exposure.

Lastly, we show April 18 (MJD 55304) in the lowest panel. We observe a quick decay and small rise of the flux, which corresponds to significant variability at 4.5 σ.

Fig. C.1. Intranight VHE variability in the two energy bands 0.2-1 TeV and > 1 TeV. The data are binned in 10 min. The dotted line shows the constant fits for the 0.2-1 TeV band.

Appendix D: Hardness ratios of VHE gamma rays and X-rays

Fig. D.1 shows the hardness ratios (HR) of all MAGIC observations as a function of the flux between 0.2-1 TeV in the left plot and above 1 TeV in the right one. Both plots show a clear harder-when-brighter trend. The rise of the HR with the flux is steep for low fluxes but flattens at higher flux levels. In the right plot, the HR seems to remain constant with a rising flux above ∼3 × 10⁻¹¹ cm⁻² s⁻¹. The plot on the left shows similar behavior, but the overall scatter of the data is slightly higher. This flattening harder-when brighter-trend at VHE gamma rays is consistent with previous reports (e.g. Acciari et al. 2021).

Fig. D.1. Hardness ratios as a function of the flux 0.2-1 TeV (left) and above 1 TeV (right) obtained by MAGIC. The color indicates the time of the observation in MJD.

Fig. D.2 shows the HRs of the Swift-XRT observations as a function of the fluxes between 0.3-2 keV in the left and 2-10 keV in the right plot. At the end of the campaign, when the source is in a low state of activity, the HR goes as low as ∼0.26. Similar to the VHE data, a clear harder-when-brighter trend is visible. On the contrary, the HR does not show a clear flattening, but a rather linear behavior. This is especially pronounced in the right plot. In the left plot for the lower-energy fluxes, the data is more scattered. The highest HR of ∼1.55, directly at the start of the campaign, still seems to follow an almost linear trend. The HR ratio corresponding to the highest flux values, however, is more compatible with a flattening of the trend. Overall, the flattening of the HR is much less pronounced than in the results from Acciari et al. (2021).

Fig. D.2. Hardness ratios as a function of the flux between 0.3-2 keV (left) and 2-10 keV (right) obtained by Swift-XRT. The color indicates the time of the observation in MJD.

Appendix E: Estimating the significances of correlations

We simulate a set of 100.000 uncorrelated pairs of light curves for each combination of energy bands, following the previous prescription, in Appendix B. To get the same degree of variability at different frequencies as the real data, we use the PSDs obtained in Sec. 4.2. The simulated light curves are generated with a temporal precision matching the typical observation time of the observations and then binned with the same temporal sampling as the real data. We then compute the coefficient for each pair. The resulting distributions are fitted with a Gaussian Kernel model to approximate the probability density function (PDF). Integration of the PDF above the coefficient of the real dataset provides a p-value. It indicates the probability of finding uncorrelated datasets that have a correlation at least as extreme as the one computed from the real dataset. We then translate the p-value into a significance expressed in levels of 1 σ. It must be emphasized that the standard approach for assessing the significance of the Pearson coefficient under the assumption of two Gaussian-distributed data sets is not applicable. The given datasets show flux distributions with strong tails towards higher fluxes. Additionally, measurement uncertainties are not taken into account in the standard method. We, therefore, rely on simulations to estimate the significance, which gives a more robust and conservative estimate.

Besides the Pearson correlation, many studies also make use of the discrete correlation function (DCF; Edelson & Krolik 1988) to quantify the correlation between fluxes in two energy bands. However, in this particular case of correlating the X-ray and the VHE gamma-ray light curves, we found that the DCF approach is less sensitive and yields somewhat underestimated values for the correlation strength. The resulting underestimation of the correlation occurs because the VHE and X-ray light curves have very different temporal coverage, with the X-rays having a much denser sampling with substantially fewer gaps. This is visible, for instance, between December 2009 and mid-January 2010, or between mid-February 2010 and early March 2010 (see Fig. 1), where no VHE gamma-ray observations are available, while the densely Swift-XRT light curve unveils strong variability. Since the normalization of the DCF depends on the standard deviation from the entire light curve and not only from the VHE/X-ray simultaneous measurements (as is the case for the Pearson coefficient), the DCF strategy, for this specific dataset is biased towards lower values, given the strong X-ray variations in periods without a VHE coverage. The mismatch in the coverage between the X-ray and VHE light curves also generates a larger spread in the DCF of the simulated light curves (used to determine its significance), further reducing the significance of the measured correlations. These effects were already reported in Max-Moerbeck et al. (2014). Therefore, we concluded that, for this specific case of the X-ray and VHE gamma-ray light curves, the DCF strategy is not adequate to properly quantify the magnitude of the correlation and its related significance. On the other hand, the local cross-correlation function (LCCF; Welsh 1999) provides an alternative approach to that of the DCF, since it also takes into account the measurement uncertainties, but differently to the DCF, the values are normalized with the standard deviation from only the coincident measurements (analogous to the Pearson coefficient). As a cross-check, we evaluated all correlations estimated from the flux-flux plots in this work with the LCCF, and we found that the values and significances are almost identical to those obtained for the Pearson coefficient.

Appendix F: Discrete correlation functions

This section presents the results of the DCF analysis referenced in the main text.

Fig. F.1. Discrete correlation function computed between the VHE gamma-ray fluxes above 0.2 TeV, as measured by MAGIC and VERITAS, and the HE gamma-ray fluxes measured by Fermi-LAT. The DCF is computed using a time bin of 3-days for a range of time lags between -40 to +40 days. The 1 σ, 2 σ, and 3 σ confidence levels obtained by simulations are shown by the yellow, red, and green lines, respectively. (a) > 0.2 TeV versus 3-300 GeV; (b) > 0.2 TeV versus 3-300 GeV without the flare in February 2010; (c) > 0.2 TeV versus 0.3-3 GeV; (d) > 0.2 TeV versus 0.3-3 GeV without the flare in February 2010.

Fig. F.2. Discrete correlation function computed between two energy ranges provided by Fermi-LAT and Swift-XRT with and without the big flare in February using a binning of 3 days. It is computed for a range of time lags between -40 to +40 days. The 1 σ, 2 σ, and 3 σ confidence levels obtained by simulations are shown by the yellow, red, and green lines, respectively. (a) 3-300 GeV versus 0.3-2 keV; (b) 3-300 GeV versus 0.3-2 keV without the flare in February 2010; (c) 0.3-3 GeV versus 0.3-2 keV; (d) 0.3-3 GeV versus 0.3-2 keV without the flare in February 2010; (e) 3-300 GeV versus 2-10 keV; (f) 3-300 GeV versus 2-10 keV without the flare in February 2010; (g) 0.3-3 GeV versus 2-10 keV; (h) 0.3-3 GeV versus 2-10 keV without the flare in February 2010.

Fig. F.3. Discrete correlation function computed between the two energy ranges provided by Fermi-LAT and the W1 filter by Swift-UVOT without the big flare in February using a binning of 3 days. It is computed for a range of time lags between -40 to +40 days. The 1 σ, 2 σ, and 3 σ confidence levels obtained by simulations are shown by the yellow, red, and green lines, respectively. (a) 3-300 GeV versus W1; (b) 3-300 GeV versus W1 without the flare in February 2010; (c) 0.3-3 GeV versus W1; (d) 0.3-3 GeV versus W1 without the flare in February 2010.

Fig. F.4. Discrete correlation function computed between the two energy ranges provided by Fermi-LAT and the R-band by GASP-WEBT without the big flare in February using a binning of 3 days. It is computed for a range of time lags between -40 to +40 days. The 1 σ, 2 σ, and 3 σ confidence levels obtained by simulations are shown by the yellow, red, and green lines, respectively. (a) 3-300 GeV versus R-band; (b) 3-300 GeV versus R-band without the flare in February 2010; (c) 0.3-3 GeV versus R-band; (d) 0.3-3 GeV versus R-band without the flare in February 2010.

All Tables

All Figures

Fig. 1. MWL light curves covering the time period from November 5, 2009, (MJD 55140) to July 3, 2010, (MJD 55380). Top to bottom: MAGIC fluxes in daily bins for two energy bands (note the two different y-axes); VHE fluxes obtained from MAGIC and VERITAS above 0.2 TeV (in log scale); Fermi-LAT fluxes in 3-day bins in two energy bands; X-ray fluxes in 1-day bins from the all-sky monitors Swift-BAT and RXTE-ASM; X-ray fluxes from the pointing instruments Swift-XRT and RXTE-PCA; hardness ratio between the high and low-energy fluxes of Swift-XRT and between the two VHE bands of MAGIC (note the two different y-axes); optical R-band data from GASP-WEBT; radio data from Metsähovi, UMRAO, SMA, OVRO and VLBA; polarization degree and polarization angle observations in the optical from the Steward and Perkins observatories and in radio from VLBA.

In the text

	Fig. 2. Total intensity VLBA images from May 2010 to August 2011. The black and yellow lines show the average positions of the stationary features A0 and A1. The yellow, red and blue circles show the fitted positions of A1, K1, and K2, respectively.
In the text

Fig. 3. Angular distance of A1 (yellow), K1 (red), and K2 (blue) from the jet core (black). At the location of Mrk 421, the angular distance of one mas is equivalent to a physical distance of 0.638 pc. The dotted yellow line indicates a constant fit of the quasi-stationary component. The red and blue lines show linear fits to determine the time of ejection of the components. The uncertainties on the ejection times are given by the red and blue bands. The campaign of this work is given by the light gray band with the three flares marked with gray lines.

In the text

Fig. 4. Fractional variability Fvar as a function of the frequency for the light curves shown in Fig. 1. Open markers show the results using the whole campaign for each light curve. Full markers only include quasi-simultaneous data to the VHE data (quasi-simultaneity is defined as temporal agreement with VHE data within 6 h for X-ray and UV data, within 1 day for optical data, and within 3 days for radio and Fermi-LAT data).

In the text

Fig. 5. VHE flux versus X-ray flux obtained by MAGIC/VERITAS and Swift-XRT/BAT. Only pairs of observations within 6 hours are considered. If more than one Swift observation falls within that window, the weighted mean is computed. VHE fluxes are in the > 1 TeV band (top panels), in the 0.2−1 TeV band (middle panels), and in the > 0.2 TeV band (bottom panels). Swift-XRT fluxes are computed in the 0.3−2 keV (left panels) and 2−10 keV bands (middle panels). Swift-BAT provides the flux in the 15−50 keV band (right panels). The top left corner of each panel shows the Pearson coefficient of the flux pairs, with the significance of the correlations given in parentheses. The gray dashed and dotted lines depict a fit with slope fixed to 0.5 and 2, respectively, and the black line is the best-fit line to the data, with the slope quoted in the legend at the bottom right of each panel.

In the text

Fig. 6. VHE flux versus UV flux obtained by MAGIC/VERITAS and Swift-UVOT. Only pairs of observations within 6 hours are considered. If more than one Swift observation falls within that window, the weighted mean is computed. VHE fluxes are in the > 1 TeV band (top panels), in the 0.2−1 TeV band (middle panels) and in the > 0.2 TeV band (bottom panels). The Swift-UVOT fluxes are taken with the W1 filter. The top left corner of each panel shows the Pearson coefficient of the flux pairs, with the significance of the correlations given in parentheses. The gray dashed and dotted lines depict a fit with slope fixed to 0.5 and 2, respectively, and the black line is the best-fit line to the data, with the slope quoted in the legend at the bottom right of each panel.

In the text

	Fig. B.1. Likelihood profile ℒ(a) (see Eq. B.1) for the 2-10 keV band with its best fit value of 1.4 denoted by a vertical blue dotted line.
In the text

	Fig. B.2. Histogram of the best-fit indices derived from simulations using as input a simulated light curve that has a known PSD index a = 1.4, in agreement with the real data for the 2-10 keV band.
In the text

	Fig. C.1. Intranight VHE variability in the two energy bands 0.2-1 TeV and > 1 TeV. The data are binned in 10 min. The dotted line shows the constant fits for the 0.2-1 TeV band.
In the text

	Fig. D.1. Hardness ratios as a function of the flux 0.2-1 TeV (left) and above 1 TeV (right) obtained by MAGIC. The color indicates the time of the observation in MJD.
In the text

	Fig. D.2. Hardness ratios as a function of the flux between 0.3-2 keV (left) and 2-10 keV (right) obtained by Swift-XRT. The color indicates the time of the observation in MJD.
In the text

Fig. F.1. Discrete correlation function computed between the VHE gamma-ray fluxes above 0.2 TeV, as measured by MAGIC and VERITAS, and the HE gamma-ray fluxes measured by Fermi-LAT. The DCF is computed using a time bin of 3-days for a range of time lags between -40 to +40 days. The 1 σ, 2 σ, and 3 σ confidence levels obtained by simulations are shown by the yellow, red, and green lines, respectively. (a) > 0.2 TeV versus 3-300 GeV; (b) > 0.2 TeV versus 3-300 GeV without the flare in February 2010; (c) > 0.2 TeV versus 0.3-3 GeV; (d) > 0.2 TeV versus 0.3-3 GeV without the flare in February 2010.

In the text

Fig. F.2. Discrete correlation function computed between two energy ranges provided by Fermi-LAT and Swift-XRT with and without the big flare in February using a binning of 3 days. It is computed for a range of time lags between -40 to +40 days. The 1 σ, 2 σ, and 3 σ confidence levels obtained by simulations are shown by the yellow, red, and green lines, respectively. (a) 3-300 GeV versus 0.3-2 keV; (b) 3-300 GeV versus 0.3-2 keV without the flare in February 2010; (c) 0.3-3 GeV versus 0.3-2 keV; (d) 0.3-3 GeV versus 0.3-2 keV without the flare in February 2010; (e) 3-300 GeV versus 2-10 keV; (f) 3-300 GeV versus 2-10 keV without the flare in February 2010; (g) 0.3-3 GeV versus 2-10 keV; (h) 0.3-3 GeV versus 2-10 keV without the flare in February 2010.

In the text

Fig. F.3. Discrete correlation function computed between the two energy ranges provided by Fermi-LAT and the W1 filter by Swift-UVOT without the big flare in February using a binning of 3 days. It is computed for a range of time lags between -40 to +40 days. The 1 σ, 2 σ, and 3 σ confidence levels obtained by simulations are shown by the yellow, red, and green lines, respectively. (a) 3-300 GeV versus W1; (b) 3-300 GeV versus W1 without the flare in February 2010; (c) 0.3-3 GeV versus W1; (d) 0.3-3 GeV versus W1 without the flare in February 2010.

In the text

Fig. F.4. Discrete correlation function computed between the two energy ranges provided by Fermi-LAT and the R-band by GASP-WEBT without the big flare in February using a binning of 3 days. It is computed for a range of time lags between -40 to +40 days. The 1 σ, 2 σ, and 3 σ confidence levels obtained by simulations are shown by the yellow, red, and green lines, respectively. (a) 3-300 GeV versus R-band; (b) 3-300 GeV versus R-band without the flare in February 2010; (c) 0.3-3 GeV versus R-band; (d) 0.3-3 GeV versus R-band without the flare in February 2010.

In the text