Unraveling the Complex Behavior of Mrk 421 with Simultaneous X-Ray and VHE Observations during an Extreme Flaring Activity in 2013 April^*

Markarian 421 (Mrk 421), with a redshift of z = 0.0308, is one of the closest BL Lac objects (Ulrich et al. 1975), which happens to also be the first BL Lac object significantly detected at gamma-ray energies (with EGRET; Lin et al. 1992) and the first extragalactic object significantly detected at very-high-energy (VHE; >0.1 TeV) gamma-rays (with Whipple; Punch et al. 1992). Mrk 421 is also the brightest persistent X-ray/TeV blazar in the sky and among the few sources whose spectral energy distributions (SEDs) can be accurately characterized by current instruments from radio to VHE (Abdo et al. 2011). Consequently, Mrk 421 is among the few X-ray/TeV objects that can be studied with a great level of detail during both low and high activity (Fossati et al. 2008; Aleksić et al. 2015b; Baloković et al. 2016) and, hence, an object whose study maximizes our chances of understanding the blazar phenomenon in general.

Because of these reasons, every year since 2009, we organize extensive multiwavelength (MWL) observing campaigns where Mrk 421 is monitored from radio to VHE gamma-rays during the half year that it is visible with optical telescopes and Imaging Atmospheric Cherenkov Telescopes (IACTs). This multi-instrument and multi-year program provides a large time and energy coverage that, owing to the brightness and proximity of Mrk 421, yields the most detailed characterization of the broadband SED and its temporal evolution compared to any other MWL campaign on any other TeV target.

During the MWL campaign in the 2013 season, in the second week of 2013 April, we observed exceptionally high X-ray and VHE gamma-ray activity with the Neil Gehrels Swift Observatory (Swift), the Nuclear Spectroscopic Telescope Array (NuSTAR), the Large Area Telescope on board the Fermi Gamma-ray Space Telescope (Fermi-LAT), the Major Atmospheric Gamma Imaging Cherenkov telescope (MAGIC), and the Very Energetic Radiation Imaging Telescope Array System (VERITAS), as reported in various Astronomer's Telegrams (e.g., see Baloković et al. 2013; Cortina & Holder 2013; Paneque et al. 2013). Among other things, the VHE gamma-ray flux was found to be two orders of magnitude larger than that measured during the first months of the MWL campaign in 2013 January and February (Baloković et al. 2016). This enhanced activity triggered very deep observations with optical, X-ray, and gamma-ray instruments, including a modified survey mode for Fermi from April 12 (23:00 UTC) until April 15 (18:00 UTC), which increased the LAT exposure on Mrk 421 by about a factor of two.

While Mrk 421 has shown outstanding X-ray and VHE gamma-ray activity in the past (e.g., Gaidos et al. 1996; Fossati et al. 2008; Abeysekara et al. 2020), this is the most complete characterization of a flaring activity of Mrk 421 to date. An extensive multi-instrument data set was accumulated during nine consecutive days. It includes VHE observations with MAGIC, the use of public VHE data from VERITAS, and high-sensitivity X-ray observations with NuSTAR. Notably, there are 43 hr of simultaneous VHE gamma-ray (MAGIC and VERITAS) and X-ray (NuSTAR) observations. A first evaluation of the X-ray activity measured with Swift and NuSTAR was reported in Paliya et al. (2015). This manuscript reports the full multiband characterization of this outstanding event, which includes, for the first time, a report of the VHE gamma-ray data, and it focuses on an unprecedented study of the X-ray-versus-VHE correlation in 3 × 3 energy bands. This study demonstrates that there is a large degree of complexity in the variability in the X-ray and VHE gamma-ray domains, which relates to the most energetic and variable segments of Mrk 421's SED, and indicates that the broadband emission of blazars requires multi-zone theoretical models.

This paper is organized as follows. In Section 2, we briefly describe the observations that were performed, and in Section 3, we report the measured multi-instrument light curves. Section 4 provides a detailed characterization of the multiband variability, with a special focus on the X-ray and VHE gamma-ray variations observed on April 15. In Section 5, we characterize the multiband correlations observed when comparing the X-ray emission in three energy bands, with that of VHE gamma-rays in three energy bands. In Section 6, we discuss the implications of the observational results reported in this paper, and finally, in Section 7, we provide some concluding remarks.

The observations presented here are part of the multi-instrument campaign for Mrk 421 that has occurred yearly since 2009 (Abdo et al. 2011). The instruments that participate in this campaign can change somewhat from year to year, but they typically consist of more than 20 covering energies from radio to VHE gamma-rays. The 2013 campaign included observations from NuSTAR for the first time, as a part of its primary mission (Harrison et al. 2013). The instruments that participated in the 2013 campaign, as well as their performances and data analysis strategies, were reported in Baloković et al. (2016), which is our first publication with the 2013 multi-instrument data set, and it focused on the low X-ray/VHE activity observed in 2013 January–March.

During the first observations in 2013 April, Mrk 421 showed high X-ray and VHE gamma-ray activity, which triggered daily few-hour-long multi-instrument observations that lasted from April 10 (MJD 56392) to April 19 (MJD 56401). Among other instruments, this data set contains an exceptionally deep temporal coverage at VHE gamma-rays above 0.2 TeV, as the source was observed with MAGIC during nine consecutive nights, and with VERITAS during six nights (Benbow & VERITAS Collaboration 2017). The geographical longitude of VERITAS is 93° (about 6 hr) west of that of MAGIC, and hence, VERITAS observations followed those from MAGIC, sometimes providing continuous VHE gamma-ray coverage during 10 hr in a single night. The total MAGIC observation time was nearly 42 hr, while for VERITAS, it was 27 hr, yielding a total VHE observation time of 69 hr in nine days (66 hr when counting the MAGIC–VERITAS simultaneous observations just once). The time coverage of VHE data is slightly different for different VHE gamma-ray energies; they are a few hours longer above 0.8 TeV in comparison to those below 0.4 TeV. This uneven coverage is due to the increased energy threshold associated with observations taken at large zenith angles. In the case of MAGIC, the energy threshold at zenith angles of about 60° is about 0.4 TeV, and hence, the low-energy gamma-ray observations are not possible (see Aleksić et al. 2016, for the dependence of the analysis energy threshold on the zenith angle of observations). MAGIC and VERITAS observed Mrk 421 simultaneously for 2 hr 45 minutes. The simultaneous observations between these two instruments occurred when MAGIC was observing at large zenith angle (>55°), hence, yielding simultaneous flux measurements only above 0.4 TeV. The extensive VHE coverage is particularly relevant since, as it will be described in Section 3, Mrk 421 showed large variability and one of the brightest VHE flaring activities recorded to date. This unprecedented brightness allows us to match the sampling frequency of the simultaneous VHE (MAGIC and VERITAS) and X-ray (NuSTAR) observations to 15 minute time intervals in three distinct energy bands: 0.2–0.4 TeV, 0.4–0.8 TeV, and >0.8 TeV. The above-mentioned time cadence and energy bands were chosen as a good compromise between having both, a good sampling of the multiband VHE activity of Mrk 421 during the 2013 April period and reasonably accurate VHE flux measurements, with relative flux errors typically below 10%. See Benbow & VERITAS Collaboration (2017) for the VERITAS photon fluxes. In the case of the MAGIC light curves in the three energy bands, the small effect related to the event migration in energy was computed using the VHE gamma-ray spectrum from the full nine day data set, which is well represented by the following log-parabola function:

$\begin{eqnarray}&&\displaystyle \frac{d{\rm{\Phi }}}{{dE}}={\left(\displaystyle \frac{E}{0.3\mathrm{TeV}}\right)}^{-(2.14)-(0.45)\cdot \mathrm{log}10\left(\tfrac{E}{0.3\mathrm{TeV}}\right)}.\end{eqnarray} \tag{ 1 }$

However, owing to the relatively narrow energy bands, the derived photon fluxes are not significantly affected by the specific choice of the used spectral shape: the photon fluxes derived with power-law spectra with indices p = 2 and p = 3 are in agreement, within the statistical uncertainties, with those derived with the nine day log-parabolic spectral shape.

Using the simultaneous MAGIC and VERITAS observations, we noted a systematic offset of about 20% in the VHE gamma-ray flux measurements derived with these two instruments. The VHE gamma-ray fluxes from VERITAS are systematically lower than those from MAGIC by a factor that is energy dependent, about 10% in the 0.2–0.4 TeV band and about 30% above 0.8 TeV. This offset, which is perfectly consistent with the known systematic uncertainties affecting each experiment (Madhavan & VERITAS Collaboration 2013; Aleksić et al. 2016), becomes evident due to the low statistical uncertainties associated with the flux measurements reported here. Appendix A reports a characterization of this offset and describes the procedure that we followed for correcting it, scaling up the VERITAS fluxes to match those from MAGIC. The physics results reported in this manuscript do not depend on the absolute value of the VHE gamma-ray flux, and hence, one could have scaled down the MAGIC fluxes to match those of VERITAS. The correction applied is only relevant for the intra-night variability and correlation studies.

A key characteristic of this data set is the extensive and simultaneous coverage in the X-ray bands provided by Swift and, especially, by NuSTAR. Swift observed Mrk 421 for 18 hr, split into 63 observations spread out over the nine days and performed during the MAGIC and VERITAS observations. NuSTAR observed Mrk 421 for 71 hr during the above-mentioned nine days, out of which 43 hr were taken simultaneously with the VHE observations from MAGIC and VERITAS. The VHE and X-ray temporal coverage is summarized in Table 1.

The raw NuSTAR data were processed exactly as described in Baloković et al. (2016), except that in this study, the NuSTAR analysis was performed separately for each 15 minute time bin with simultaneous VHE observations, as summarized in Table 1. Using Xspec (Arnaud 1996), we calculated fluxes in the 3–7 keV, 7–30 keV, and 30–80 keV bands from a fit of a log-parabolic model to the data within each time bin. The cross-normalization between the two NuSTAR telescope modules was treated as a free parameter. The statistical uncertainties of the fluxes were calculated at 68% confidence intervals and do not include the systematic uncertainty in the absolute calibration, which is estimated to be 10%–20% (Madsen et al. 2015).

The analysis procedures used to process the Swift-XRT data are described in Baloković et al. (2016). In addition, in order to avoid additional flux uncertainties, we excluded 16 Swift-XRT observations in which Mrk 421 was positioned near the CCD bad columns (Madsen et al. 2017). Figure 1 shows a comparison of the Swift-XRT and NuSTAR X-ray fluxes in the band 3–7 keV. Overall, there is a good agreement between the two instruments, with flux differences typically smaller than 20%. Such flux differences are within the systematic uncertainties in the absolute flux calibration of NuSTAR (Madsen et al. 2015) and Swift-XRT (Madsen et al. 2017).

Zoom In Zoom Out Reset image size

Differently to the Fermi-LAT analysis reported in Baloković et al. (2016), the LAT data results shown here were produced with events above 0.3 GeV (instead of 0.1 GeV) and with Pass8 (instead of Pass7). The analysis above 0.3 GeV is less affected by systematic uncertainties, and it is also less sensitive to possible contamination from non-accounted (transient) neighboring sources. The higher minimum energy somewhat reduces the detected number of photons from the source, but, owing to its hard gamma-ray spectrum (photon index <2.0), the effect is small. Specifically, we used the standard Fermi analysis software tools version v11r07p00, and the P8R3_SOURCE_V2 response function on events with energy above 0.3 GeV coming from a 10° region of interest (ROI) around Mrk 421. We used a 100° zenith-angle cut to avoid contamination from the Earth's limb, ⁶⁷ and we modeled the diffuse Galactic and isotropic extragalactic background with the files gll_iem_v07.fits and iso_P8R3_SOURCE_V2_v1.txt, respectively. ⁶⁸ All point sources in the fourth Fermi-LAT source catalog (4FGL; Abdollahi et al. 2020) located in the 10° ROI and an additional surrounding 5°-wide annulus were included in the model. In the unbinned likelihood fit, the spectral parameters were set to the values from the 4FGL, while the normalization of the diffuse components and the normalization parameters of the 16 sources (within the ROI) identified as variable were initially left free to vary. However, owing to the short timescales considered in this analysis, only two of these sources were significantly detected in 10 days: 4FGL J1127.8+3618 and 4FGL J1139.0+4033, and hence, we fixed the normalization of the other ones to the 4FGL catalog values. The Fermi-LAT spectrum from the 10 day time period considered here (from MJD 56392 to MJD 56402) is well described by a power-law function with a photon flux above 0.3 GeV of (23.3 ± 1.6) × 10⁻⁸ cm⁻² s⁻¹ and photon index 1.79 ± 0.05. A spectral analysis over 1 day and 12 hr time intervals shows that the photon index does not vary significantly throughout the 10 day period. The data was split into 12 hr–long intervals centered at the VHE observations (e.g., simultaneous to the VHE) and their complementary time intervals (e.g., when there are no VHE observations), which are close to 12 hr–long intervals. Owing to the limited event count in the 12 hr time intervals and the lack of spectral variability throughout the 10 day period, we fixed the shape to a power-law index of 1.79 (the value from the 10 day period) to derive the photon fluxes, always keeping the normalization factor for Mrk 421 and the two above-mentioned 4FGL sources as free parameters in the log-likelihood fit of each of the 12 hr time intervals.

The characterization of the activity of Mrk 421 at optical frequencies was performed with many instruments from the GLAST-AGILE Support Program (GASP) of the Whole Earth Blazar Telescope (WEBT), hereafter GASP-WEBT (e.g., Villata et al. 2008, 2009), namely, the observatory in Roque de los Muchachos (KVA telescope), Lowell (Perkins telescope), Crimean, St. Petersburg, Abastumani, Rozhen (50/70 cm, 60 cm, and 200 cm telescopes), Vidojevica, and Lulin. Moreover, this study also uses data from the iTelescopes, the Remote Observatory for Variable Object Research (ROVOR), and the TUBITAK National Observatory (TUG). The polarization measurements were performed with four observatories: Lowell (Perkins telescope), St. Petersburg, Crimean, and Steward (Bok telescope). The data reduction was done exactly as in Baloković et al. (2016).

Other than the 15 and 37 GHz radio observations performed with the OVRO and Metsahovi telescopes, which were described in Baloković et al. (2016), here, we also present a flux measurement performed with the IRAM 30 m telescope at 86 GHz. This observation was performed under the Polarimetric Monitoring of AGNs at Millimeter Wavelengths program (POLAMI, ⁶⁹ Agudo et al. 2018b), which regularly monitors Mrk 421 in the short millimeter range. The POLAMI data was reduced and calibrated as described in Agudo et al. (2018a).

The multi-instrument light curves derived from all of the observations spanning from radio to VHE gamma-rays are shown in Figure 1. The top panel of Figure 1 shows an excellent coverage of the nine day flaring activity in the VHE regime as a result of the combined MAGIC and VERITAS observations. The peak flux at TeV energies, observed on April 13 (MJD 56395), reached up to 15 times the flux of the Crab Nebula, which is about 30 times the typical non-flaring activity of Mrk 421 and about 150 times the activity shown a few months before, on 2013 January and February, as reported in Baloković et al. (2016). Moreover, this is the highest TeV flux ever measured with MAGIC for any blazar. This is also the third highest flux ever measured from a blazar with an IACT, after the extremely large outburst from Mrk 421 detected with VERITAS in 2010 February (Abeysekara et al. 2020) and the large flare from PKS 2155-304 detected by HESS in 2006 July (Aharonian et al. 2007).

Figure 1 shows that the most extreme flux variations occur in the X-ray and VHE gamma-ray bands. At GeV energies, within the accuracy of the measurements, there is enhanced activity only for MJD 56397 (April 15), when the flux is about a factor of two larger than the flux in the previous and in the following ∼12 hr time intervals. Interestingly, on April 15, we also find the highest X-ray flux and the highest intra-night X-ray flux increase measured during this flaring activity in 2013 April.

The R-band activity is comparable to the one measured in 2013 January–March, when Mrk 421 showed very low VHE and X-ray activity (Baloković et al. 2016). The measured fluxes at optical wavelengths are large when compared to the flux levels typically seen during the period of 2007–2015 (Carnerero et al. 2017). Generally, during the observations performed in 2013 January–April, Mrk 421 was four to five times brighter in the optical than were the photometric minima that occurred in 2008–09 and at the end of 2011. Figure 1 shows that Mrk 421 faded at the R band from about 60 mJy on MJD 56393 to about 45 mJy two days later. It then varied between 45 and 50 mJy during the following week and appeared decoupled from the VHE and X-ray activity. The optical light curve is in agreement with that of the less well-sampled Swift/UVOT light curve. Other than the optical brightness of Mrk 421 in 2013, the object showed a bluer optical continuum than average. This was determined from the differential spectrophotometry obtained by the Steward Observatory monitoring program. By comparing the instrumental spectrum of Mrk 421 with that of a nearby field comparison star, it is found that for the wavelengths 475 and 725 nm that [F(475)/F(725)]_{2013 April}/[F(475)/F(725)]_average = 1.072 ±  0.002, where the average instrumental flux ratio is determined from all of the available observations from 2008 to 2018. The bluer color of Mrk 421 is consistent with a higher dominance of the nonthermal continuum over the host galaxy starlight included within the observing aperture, which has a redder spectrum. This explanation for the observed variations in the optical color of Mrk 421 is further confirmed by the trend that the continuum becomes slightly redder, as the AGN generally fades during 2013 April. The same trend in color is also seen in the long-term near-IR data (Carnerero et al. 2017).

The optical linear polarization of Mrk 421 was also monitored, and the measurements are shown in Figure 1. Again, the results are comparable to those measured during the first quarter in 2013 and reported in Baloković et al. (2016). Since 2008, the degree of polarization, P, has ranged from 0% to 15%; although, observations of P > 10% are rare, about 10 out of around 1400 observations (Carnerero et al. 2017). During 2013 April, the polarization ranged from about 1%–9%, with a large majority of measurements showing P < 5%. The largest changes in the degree of polarization on a daily timescale were an increase from P ∼ 3% to P ∼ 7% on MJD 56399/400, followed by a decrease back to about 4% on the next day. Changes of nearly as much as 5% in polarization are observed within a day, particularly for MJD 56398/99, but otherwise, variations in P are typically limited to <1% over hour timescales. The electric vector position angle (EVPA) of the optical polarization was at about −20° at the start of 2013 April. Between MJD 56394 and MJD 56395, the EVPA rotated from about −30° to about −90° while, generally, P < 2.5%. The largest daily rotation in EVPA occurs between MJD 56395 and MJD 56396, where the EVPA goes from about −90° to about −10° . Because of the daily gap in the optical monitoring, it is unfortunately not clear if the EVPA reversed its direction of rotation from MJD 56394 to MJD 56396 (i.e., 2 days) or continued in the same direction requiring a rotation >90° during one of the two observing gaps for MJD 56394/6. The variability of Mrk 421 during the densely sampled portions of the optical monitoring does not hint that such large changes in EVPA can take place on short timescales until near the the end of MJD 56398 when a counterclockwise rotation of about 50° is seen over a period of about 6 hr. Outside of this excursion, the EVPA stays near 0° from MJD 56396 onward. The single daily deviation of the EVPA to 90°–100° for MJD 56395 coincides with brightest VHE flare observed in 2013 April. However, no significant change in the EVPA is apparent during the sharp rise in VHE flux observed near the middle of MJD 56394 or during the dramatic high-energy activity at the beginning of MJD 56397. For most of the monitoring period, the optical EVPA was near the historically most likely angle for this object (EVPA = 0°; Carnerero et al. 2017); although, the one-day excursion on MJD 56395 brought the EVPA nearly orthogonal to the most likely value. For comparison, the 15 GHz VLBI maps of Mrk 421 show a jet detected out to about 5 mas at a position angle of about −40° (Lister et al. 2019). In the radio band, the activity measured during the entire nine day observing period is constant, with a flux of about 0.6 Jy. The single 86 GHz measurement with IRAM 30 m shows a polarization degree of about 3%, which is similar to that of the optical frequencies; yet, the polarization angle differs by about 70°, which suggests that the optical and radio emissions are being produced in different locations of the jet of Mrk 421. Overall, the radio and optical fluxes, as well the optical polarization variations (polarization degree and EVPA), appear completely decoupled from the large X-ray and VHE gamma-ray activity seen in 2013 April. In fact, the behavior observed at radio and optical during 2013 April is similar to that observed during the previous months, when Mrk 421 showed extremely low X-ray and VHE gamma-ray activity (see Baloković et al. 2016).

The hard X-ray and VHE gamma-ray bands covered with NuSTAR, MAGIC, and VERITAS are the most interesting ones because they exhibit the largest flux variations and because of the exquisite temporal coverage and the simultaneity in the data set. Figure 2 reports the flux measurements in these bands, each split into three distinct bands, 3–7 keV, 7–30 keV, and 30–80 keV for ⁷⁰ NuSTAR and 0.2–0.4 TeV, 0.4–0.8 TeV, and >0.8 TeV for MAGIC and VERITAS. The temporal coverage for the >0.8 TeV band is a about 2 hr longer than for the 0.2–0.4 TeV band because of the increasing analysis energy threshold with the increasing zenith angle of the observations. The exquisite characterization of the multiband flux variations in the X-ray and VHE gamma-ray bands reported in Figure 2 will be used in the sections that follow for the broadband variability and correlation studies.

Zoom In Zoom Out Reset image size

4.1. Fractional Variability

The flux variability reported in the multiband light curves can be quantified using the fractional variability parameter F_var, as prescribed in Vaughan et al. (2003):

$\begin{eqnarray}&&{F}_{\mathrm{var}}=\sqrt{\displaystyle \frac{{S}^{2}-\langle {\sigma }_{\mathrm{err}}^{2}\rangle }{\langle {F}_{\gamma }{\rangle }^{2}}}.\end{eqnarray} \tag{ 2 }$

$\langle {F}_{\gamma }\rangle$ denotes the average photon flux, S denotes the standard deviation of N flux measurements, and $\langle {\sigma }_{\mathrm{err}}^{2}\rangle$ denotes the mean squared error, all determined for a given instrument and energy band. The uncertainty on F_var is calculated using the prescription from Poutanen et al. (2008), as described in Aleksić et al. (2015a). This formalism allows one to quantify the variability amplitude, with uncertainties dominated by the flux measurement errors and the number of measurements performed. The systematic uncertainties on the absolute flux measurements ⁷¹ do not directly add to the uncertainty in F_var. The caveats in the usage of F_var to quantify the variability in the flux measurements performed with different instruments are described in Aleksić et al. (2014, 2015a, 2015b). The most important caveat is that the ability to quantify the variability depends on the temporal coverage (observing sampling) and the sensitivity of the instruments used, which are somewhat different across the electromagnetic spectrum. A big advantage of the study presented here is that the temporal coverage of the three bands in X-rays (from NuSTAR) and the three bands in VHE gamma-rays (from MAGIC and VERITAS) is exactly the same, which allows us to make a more direct comparison of the variability in these energy bands.

The fractional variability parameter F_var was computed using the flux values and uncertainties reported in the light curves from Section 3 (see Figures 1 and 2), hence, providing a quantification of variability amplitude for this nine day–long flaring activity from radio to VHE gamma-ray energies. The results are depicted in the upper panel of Figure 3, where open markers are used for the variability computed with all of the available data, and filled markers are used for simultaneous observations. Given the slightly different temporal coverage for the different VHE bands, as described in the previous section, we decided to use the 0.2–0.4 TeV band to define the time slots for simultaneous X-ray/VHE observations. This ensures that the same temporal bins are being used for the 3 × 3 X-ray and VHE bands. For comparison purposes, we added the F_var values obtained for the period from 2013 January to March, when Mrk 421 showed very low activity (see Baloković et al. 2016).

Zoom In Zoom Out Reset image size

The fractional variability plot shows the typical double-bump structure, which is analogous to the broadband SED. This plot shows that most of the flux variations occur in the X-ray and VHE bands, which correspond to the falling segments of the SED. Additionally, it also shows that, during the nine day flaring activity in 2013 April, the amplitude variability in the hard X-ray band was substantially larger than that measured during the low activity from 2013 January to March. The higher the X-ray energy, the larger the difference between the F_var values from the low and the high activity.

In addition to the study of the nine day behavior, the high photon fluxes and the deep exposures allow us to compute F_var with the single-night light curves from six consecutive nights (from April 11 to 17), ⁷² hence, allowing us to study the fractional variability on hour timescales for the three X-ray bands and three VHE gamma-ray bands. For this study, only simultaneous data (using the time bins from the 0.2–0.4 TeV band) were used, which means that the 3 × 3 X-ray/VHE bands sample exactly the same source activity. The results are depicted in the lower panels of Figure 3. In general, all F_var values computed with the single-night light curves are lower than those derived with the nine day light curve for the corresponding energy band. This is clearly visible when comparing the data points with the gray shaded regions in the upper panel of Figure 3. Despite the X-ray and VHE flux varying on sub-hour timescales, the resulting intra-night fractional variability is significantly lower than the overall fractional variability in the nine day time interval. This result is expected because, while for single days, the light curves show flux variations within a factor of about two, the nine day light curve shows flux variations larger than a factor of about 10. Unexpectedly, we find a large diversity in the variability versus energy patterns observed for the different days. On April 13, 15, 16, and 17, one finds a typical pattern of higher fractional variabilities at higher photon energies within each of the two SED bumps. On the other hand, one finds that the fractional variability is approximately constant with energy on April 11 and 14 and that the fractional variability decreases with energy on April 12. The decrease in F_var with increasing energy is only marginally significant in the VHE bands (∼2σ) but very prominent in the X-ray emission, within which the synchrotron self-Compton (SSC) scenario, provides a direct mapping to the energy of the radiating electrons. These different variability versus energy patterns suggest the existence of diverse causes (or regions) responsible for the variability in the broadband blazar emission on timescales as short as days and hours. This is the first time that the variability of Mrk 421 has been studied with this level of detail, and the implications will be discussed in Section 6.

4.2. Flux Variations on Multi-hour and Sub-hour Timescales

This section focuses on the flux variations observed in the hard X-ray and VHE gamma-ray bands, which are the ones with the largest temporal coverage and highest variability (see Section 4.1). The light curves for all nights for these 3 × 3 energy bands are reported in Appendix B. There is clear intra-night variability in all of the light curves, which can be significantly detected because of the high fluxes and the good temporal coverage, as described in the previous section (e.g., see Figure 3). The single-night light curves show a large diversity of temporal structures that relate to different timescales, from sub-hours (i.e., fast variation) to multi-hours (trends). We note that some of these fast components are present in both X-rays and VHE gamma-rays, while some others are visible only at X-rays or only at VHE gamma-rays, and, in some cases, the features are present only in specific bands (either X-rays or VHE gamma-rays) and not in the others. As it occurred with the study of the F_var versus energy, the evaluation of the single-night multiband light curves also suggests that there are different mechanisms responsible for the variability, some of them being achromatic (affecting all energies in a similar way) and others chromatic (affecting the different energy bands in a substantially different manner).

In this section, we attempt to quantify the main trends and fast features, as well as their evolution across the various energy bands. We do that by fitting with a function formed by a slow trend F_s (t) and fast feature F_f (t) components.

$\begin{eqnarray}&&F(t)={F}_{s}(t)+{F}_{f}(t)\end{eqnarray} \tag{ 3 }$

where

$\begin{eqnarray}&&{F}_{s}(t)=\mathrm{Offset}\cdot (1+\mathrm{Slope}\cdot t)\end{eqnarray} \tag{ 4 }$

and

$\begin{eqnarray}&&{F}_{f}(t)=\displaystyle \frac{2}{{2}^{-\tfrac{t-{t}_{0}}{{t}_{\mathrm{rise}}}}+{2}^{\tfrac{t-{t}_{0}}{{t}_{\mathrm{fall}}}}}\cdot A\cdot {F}_{s}({t}_{0}).\end{eqnarray} \tag{ 5 }$

Here, A is the flare amplitude, t is the time since midnight for the chosen night, t₀ is the time of the peak flux of the flare, and t_rise and t_fall are the flux-doubling timescales for the rising and falling part of the flare, respectively. This formulation, with the slope of the slow component normalized to the offset, and the flare amplitude (of the fast component) normalized to the slow component at t₀, enables a direct comparison of the parameter values among the different energy bands, for which the overall measured flux may differ by factors of a few.

In general, we find that, whenever fast flares occur, they appear to be quite symmetric and, given the relatively short duration (sub-hour timescales) and the flux measurement uncertainties, we do not have the ability to distinguish (in a statistically meaningful way) between different rise- and fall-doubling times. For the sake of simplicity, we decided to fit the light curves with a function given by Equation (5) where t_rise = t_fall = flux-doubling time. This fit function provides a fair representation of the intra-night rapid flux variations from all days except for April 16, where the flux variations have much longer (multi-hour) timescales.

This relatively simple function provides a rough description of the energy-dependent light curves and may not describe all of the data points perfectly well. For instance, in the low-energy X-ray bands, the statistical uncertainties are very small, and one can appreciate the significant and complex substructure that is not reproduced by the above-described (and relatively simple) fitting function. We do not intend to find a model that describes all of the data points accurately. Rather, we look for a model that provides a description of the main flux-variability trends and how they evolve with the X-ray and VHE energies.

The multiband flux variations during April 15 and its related quantification using Equation (3) are depicted in Figure 4, with the parameters resulting from the fits reported in Table 2. The main multiband emission varies on timescales of several hours, and hence, it is dominated by the "slow component" in Equation (3). The slope of this variation (quantified relative to the offset in each band for better comparison among all bands) has a strong energy dependence, with the parameter value for the highest energies being around a factor of two to three times larger than that for the lowest energies for both the X-ray and VHE gamma-ray bands (e.g., Slope_{>0.8 TeV} ≃ 3 · Slope_{0.2–0.4 TeV}). The second most important feature of this multiband light curve is the existence of a short flare, on the top of the slowly varying flux, in all of the energy bands for both X-rays and VHE gamma-rays. The location of the flare t₀ is the same (within uncertainties) in all three X-ray bands and VHE bands. In order to better quantify the location of the short flare at the X-ray and VHE gamma-ray energies, using the information from Table 2, we computed the weighted average separately for the three VHE bands, t_{0, VHE} = 2.44 ± 0.03 hr and the three X-ray bands, t_{0, X-ray} = 2.41 ± 0.04 hr past midnight. This indicates that, for this fast feature in the light curve, for all of the energies probed, there is no delay in between the X-ray and VHE gamma-ray emissions down to the resolution of the measurement, which, adding the errors in quadrature, corresponds to 3 minutes. The flux-doubling time is comparable among all of the energy bands, with about 0.3 hr for all of the X-ray bands and the highest VHE band (>0.8 TeV), and about 0.2 hr for the lowest and middle VHE bands. The characteristic by which the fast X-ray flare differs from the fast VHE flare is in the normalized flare amplitude A (see Table 2): it is energy independent (achromatic) for the X-ray fast flare, while it increases its value (chromatic) for the VHE fast flare (with amplitude A_{>0.8 TeV} ≃ 2 · A_{0.2–0.4 TeV}).

Zoom In Zoom Out Reset image size

Table 2. Parameters Resulting from the Fit with Equation (3) to the X-Ray and VHE Multiband Light Curves from 2013 April 15

Band	Offset a	Slope	Flare	Flare	Flare	χ2/d.o.f.
		(hr−1)	Amplitude A	Flux-doubling Time b (hr)	t0 (hr)
2013 Apr 15
3–7 keV	0.71 ± 0.01	0.153 ± 0.006	0.49 ± 0.07	0.30 ± 0.04	2.35 ± 0.06	836/24
7–30 keV	0.78 ± 0.02	0.199 ± 0.009	0.59 ± 0.11	0.30 ± 0.04	2.41 ± 0.06	889/24
30–80 keV	0.21 ± 0.01	0.241 ± 0.018	0.56 ± 0.18	0.32 ± 0.09	2.50 ± 0.10	111/24
0.2–0.4 TeV	6.60 ± 0.17	0.031 ± 0.008	0.40 ± 0.09	0.23 ± 0.07	2.41 ± 0.09	96.9/38
0.4–0.8 TeV	2.99 ± 0.07	0.042 ± 0.008	0.72 ± 0.09	0.19 ± 0.03	2.47 ± 0.04	68.1/42
>0.8 TeV	1.68 ± 0.05	0.103 ± 0.010	0.82 ± 0.08	0.27 ± 0.03	2.41 ± 0.04	90.0/45

Notes.

^a For VHE bands in 10⁻¹⁰ ph cm⁻² s⁻¹, for X-ray bands in 10⁻⁹ erg cm⁻² s⁻¹. ^b Parameters t_rise and t_fall in Equation (3) are set to be equal and correspond to the values of flare flux-doubling time.

Download table as:  ASCIITypeset image

In order to evaluate potential spectral variability throughout the ∼10 hr light curves measured on 2013 April 15, we computed the flux hardness ratios HR (=F_high-energy/F_low-energy) for several energy bands in both the X-ray and VHE gamma-ray domains. Figure 5 depicts the HR computed with the data flux measurements (in time bins of 15 minutes) and the HR expected from the fitted functions reported in Figure 4 and Table 2. For comparison purposes, we also included the HR from the fitted functions from Table 2 excluding the fast component given by Equation (5) (dotted line in Figure 5). The dashed vertical red line indicates the weighted average time of the peak of the flare t₀, calculated separately for the three X-ray bands and VHE gamma-ray bands (see above). One can see that the overall impact of the fast component in the HR temporal evolution is small and only noticeable in some panels (e.g., F_{>0.8 TeV}/F_{0.4–0.2 TeV} or F_{7–30 keV}/F_{3–7 keV}). This is due to the relatively short duration of the fast component and the relatively small magnitude of the flare amplitude, in comparison to the overall flux. Therefore, the temporal evolution of the spectral shape in both bands, X-ray and VHE, is dominated by the slow component, i.e., by the variations with timescales of several hours.

Zoom In Zoom Out Reset image size

In addition to April 15, we also performed the fit with Equation (3) to the other five consecutive nights with large X-ray and VHE gamma-ray simultaneous data sets, namely all nights from April 11 to 16 (both included). The results from these fits are reported in Appendix C (see Table C1 and Figures C1–C5). It is worth stating that, when comparing the quantification of the various light curves with the function in Equation (3), we found diversity among the fit parameter values and their energy dependencies. For April 11, we did not find any fast component, and the flux decreases monotonically through the observation with energy-independent slope for both the X-rays and VHE gamma-rays (fully achromatic flux variations). On the other hand, during April 12, the emission increased throughout the observation but with a slope that decreased with increasing energy in both the X-rays and VHE gamma-rays. This trend is also observed, from a different perspective, in the lower-right panel of Figure 3, which displays a decreasing F_var with increasing energy for both the X-ray and VHE gamma-ray emissions from April 12. This is a very interesting behavior because it is opposite to the trend reported in most data sets from Mrk 421, where the variability increases with energy. For this night, we can also see a fast X-ray flare (flux-doubling time of about 0.3 hr) whose amplitude increases with energy. Unfortunately, this fast X-ray flare occurred during a time window without VHE observations.

For April 13, we also observed a slow flux variation with an energy-independent slope, as on April 11, but this time with a flux increase instead of a decrease. Additionally, we did observe a super-fast X-ray flare (flux-doubling time of 5 ± 1 minutes) without any counterpart in the VHE light curve, i.e., an "orphan" X-ray flare (see Figure C3). As shown in Table C1, the X-ray NusTAR flare amplitude relative to the overall baseline is only about 11%, but it is significant (3–4σ depending on the energy band), and there is no correlated flux variation in the simultaneous VHE MAGIC fluxes, which have flux uncertainties of about 5%.

On April 14, we see again a monotonically decreasing flux with an energy-independent slope for both the X-rays and VHE gamma-rays, with another fast X-ray flare (flux-doubling time ∼0.5 hr) without a counterpart in the VHE light curve.

The night that differs the most is April 16, which does not show any monotonic increase or decrease and shows a largely nonsymmetric flare with flux variation timescales of hours. In order to quantify the temporal multiband evolution of the flux during April 16, we used Equation (5) (i.e., the fitting function without the slow component), with ${t}_{\mathrm{rise}}\ne {t}_{\mathrm{fall}}$. See Appendix C for further details about the quantification of the multiband flux variations during the six consecutive nights, from April 11 to 16.

In summary, during these six consecutive nights with enhanced activity and with multi-hour-long X-ray/VHE simultaneous exposures in 2013 April, we found achromatic and chromatic flux variability, with timescales spanning from multi-hours to sub-hours, and several X-ray fast flares without VHE gamma-ray counterparts. We did not see any VHE gamma-ray orphan fast flares (whenever we had simultaneous X-ray coverage). However, we did observe fast flares in some specific energy bands that are not detected in the other nearby energy bands (X-ray or VHE), which suggests the presence of flaring mechanisms affecting relatively narrow energy bands.

The temporal evolution of the X-ray and VHE emissions, and the particularity of being able to approximately describe it with a two-component function with a fast (sub-hour variability timescale) and slow (multi-hour variability), will be discussed in Section 6.

We evaluated the correlations among all of the frequencies covered during the 2013 April flare and found that the largest flux variations and the largest degree of flux correlation occurs in the X-ray and VHE gamma-ray bands. No correlation was found among the radio, optical, and gamma-ray bands, a result that was expected because of the lower activity and longer variability timescales at these energies. Apart from some variability in the GeV flux around April 15, which is the day with the highest X-ray activity, the GeV emission appears constant for the 12 hr time intervals related to flux variations by factors of a few at keV and TeV energies. If the GeV and TeV fluxes were correlated on 12 hr timescales, Fermi-LAT should have detected large flux variations, and hence, we can exclude this correlation.

The quality and extent of this data set, both in time and energy, allows for an X-ray/VHE correlation study that is unprecedented among all data sets collected from Mrk 421 and any other TeV blazar. The relation between the VHE gamma-ray and the X-ray fluxes in the 3 × 3 energy bands is shown in Figure 6 for the nine day flaring activity and in Figure 7 for April 15. The discrete correlation function (DCF) and Pearson correlation coefficients, as well as the slope of the VHE versus X-ray flux, are reported in Table 3. There is a clear pattern: the strength of the correlation increases for higher VHE bands and lower X-ray bands. The strongest correlation is observed between the 3–7 keV and >0.8 TeV bands. This combination of bands also shows a slope (from the fit in Table 3) closest to 1, among all of the 3 × 3 bands reported. Moreover, the scatter in the plots becomes smaller as we increase the VHE band and decrease the X-ray energy band. The smallest scatter, which can be quantified with the χ² of the fit (lower values of χ² relate to a smaller scatter in the data points), occurs for the combination >0.8 TeV and 3–7 keV.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3. Correlation Coefficients and Slopes of the Linear Fit to the VHE vs. X-Ray Flux (in Log Scale) Derived with the Nine Day Flaring Episode of Mrk 421 in 2013 April

VHE Band	X-Ray Band	Pearson Coeff. a	Nσ Pearson a	DCF	Linear Fit Slope	χ2/d.o.f.
0.2–0.4 TeV	3–7 keV	0.92 ± 0.01	20.2	0.93 ± 0.12	0.61 ± 0.02	1183/162
	7–30 keV	0.87 ± 0.02	17.0	0.88 ± 0.11	0.45 ± 0.03	1891/162
	30–80 keV	0.79 ± 0.03	13.6	0.81 ± 0.11	0.35 ± 0.02	2277/162
0.4–0.8 TeV	3–7 keV	${0.946}_{-0.009}^{+0.007}$	23.4	0.96 ± 0.11	0.79 ± 0.03	1038/170
	7–30 keV	0.91 ± 0.01	19.8	0.92 ± 0.11	0.58 ± 0.03	1725/170
	30–80 keV	0.84 ± 0.02	15.8	0.86 ± 0.11	0.45 ± 0.03	2160/170
>0.8 TeV	3–7 keV	${0.964}_{-0.006}^{+0.005}$	26.0	0.97 ± 0.11	1.11 ± 0.03	704/170
	7–30 keV	${0.947}_{-0.008}^{+0.007}$	23.5	0.96 ± 0.11	0.81 ± 0.03	1245/170
	30–80 keV	0.89 ± 0.02	18.6	0.91 ± 0.10	0.61 ± 0.03	1736/170

Note.

^a The Pearson correlation function 1σ errors and the significance of the correlation are calculated following Press et al. (2007).

Download table as:  ASCIITypeset image

Figure 6 reveals that the different days occupy (roughly) different regions in the VHE versus X-ray flux plots (for all of the 3 × 3 bands). This is expected because the largest flux changes occur on day-long timescales. In addition, individual days appear to show different patterns. In order to better characterize these different patterns (observed for the different days), we also computed the same quantities (DCF, Pearson, and linear fit) to the simultaneous data points from the single nights with multi-hour light curves (namely, April 11–16). The results are reported in Table D1, in Appendix D.

The main conclusions from this study performed on data from April 11 to 16 are as follows ⁷³ :

• 1.  
  During some nights, namely on April 15 and 16, Mrk 421 shows the "general trend" that is observed for the full nine day flaring activity, with the highest magnitude and significance in the correlation occurring for the >0.8 TeV versus the 3–7 keV bands.
• 2.  
  During other nights, namely on April 11, 13, and 14, the general trend from the nine day data set is less visible: there is a larger similarity in the magnitude and significance of the correlation among the various energy bands.
• 3.  
  During one night, April 12, we found no correlation between the X-ray and VHE gamma-ray bands, despite the significant variability in both bands. The lack of correlation between the X-rays and VHE gamma-rays (being both highly variable and characterized with simultaneous observations) has not been observed to date in Mrk 421; although, it has been observed in another HBL (PKS 2155-304, from 2008 August–September; Aharonian et al. 2009).

The correlation study was also performed splitting the data set into two subsets: (a) April 15, 16, and 17 (which appear somewhat away from the main trend in Figure 6), and (b) April 11, 12, 13, 14, 18, and 19. In comparison to the nine day data set, the separate analysis of the two subsets yields a reduction of the scatter of the flux points around the main trends (which show up as a noticeable reduction in the χ² values from the fit) and also a smaller dependence of the magnitude and significance of the correlations on the specific combination of VHE and X-ray energy bands. The largest linear fit slope occurs for the combination >0.8 TeV and 3–7 keV in both subsets (as with the nine day data set), but while for subset (a) we continue having the largest significance and magnitude of the correlation for >0.8 TeV and 3–7 keV, in subset (b), the change in magnitude and correlation with energy bands is much smaller, and the highest values occur for >0.8 TeV and 7–30 keV. This indicates a somewhat different physical state of the source during those three consecutive days (April 15, 16, and 17) with respect to the others.

In addition to the VHE gamma-ray and the X-ray fluxes in the 3 × 3 energy bands for April 15, Figure 7 also depicts the flux–flux values from the fitted functions reported in Section 4.2. This figure shows that multiple components in the flux evolution (e.g., the fast component on the top of the slow component) appear as "different trends" in the flux–flux plots with the flaring component having a sharper VHE flux rise (with increasing X-ray flux) than the slow component. Because of the statistical uncertainties in the flux measurements, as well as the fact that one component has a much smaller flux and shorter duration, even for very good data sets such as this one, it is not easy to recognize and separate the contribution of different components in the flux–flux plots. However, these different patterns can produce collective deviations (when considering many of these different single-trends) that are statistically significant when fitting the data points in the flux–flux plots with simple trends, such as the linear or quadratic functions in the log–log scale.

A discussion of these observational results is given in Section 6.

Although detailed SED and light-curve modeling are beyond the scope of this paper, we discuss the main results of our analysis and provide possible interpretations.

6.1. Minimum Doppler Factor

6.2. Flux–Flux Correlations

6.3. Temporal Variability

One of the main results of this work is the detection of fast-evolving flares on top of a slower evolving emission in both the X-ray and VHE gamma-ray bands. Although this temporal behavior was qualitatively discussed for some of the X-ray NuSTAR light curves in Paliya et al. (2015), here, we present a quantitative study of these characteristics in both the X-ray light curves and the VHE gamma-ray light curves (see Figure 4 in Section 4.2 and Figures C1–C5 in Appendix C). In the following paragraphs, we discuss possible interpretations for the origin of the multiband temporal variability.

6.3.1. Acceleration and Cooling Processes

The rise and decay timescale of a flare may be associated with the acceleration and cooling timescales of the radiating electrons. In this case, one can use the fact that the acceleration and cooling timescales are found to be equal in order to estimate the magnetic field of the emitting region, as follows. If electrons undergo Fermi-1 (or Fermi-2) acceleration, then the acceleration timescale, in the co-moving frame, can be written as

$\begin{eqnarray}&&{t}_{\mathrm{acc}}^{{\prime} }=\displaystyle \frac{2\pi {m}_{e}c{\gamma }^{{\prime} }{N}_{a}}{{{eB}}^{{\prime} }}\end{eqnarray} \tag{ 6 }$

(e.g., Finke et al. 2008) where B' is the tangled (co-moving) magnetic field strength, ${\gamma }^{{\prime} }$ is the electron Lorentz factor in the co-moving frame, and N_a  ≥ 1 is the number of gyrations an electron makes to double its energy. The synchrotron cooling timescale is

$\begin{eqnarray}&&{t}_{\mathrm{syn}}^{{\prime} }=\displaystyle \frac{6\pi {m}_{e}{c}^{2}}{4c{\sigma }_{{\rm{T}}}{B}^{{\prime} 2}{\gamma }^{{\prime} }}.\end{eqnarray} \tag{ 7 }$

Setting ${t}_{\mathrm{acc}}^{{\prime} }={t}_{\mathrm{syn}}^{{\prime} }$ results in

$\begin{eqnarray}&&{B}^{{\prime} }=\displaystyle \frac{3e}{4{\sigma }_{{\rm{T}}}{\gamma }^{{\prime} 2}{N}_{a}}.\end{eqnarray} \tag{ 8 }$

The Lorentz factor of electrons producing the peak of the SSC mission can be estimated as

$\begin{eqnarray}&&{\gamma }_{b}^{{\prime} }\approx \displaystyle \frac{{E}_{\mathrm{ssc}}}{\delta {m}_{e}{c}^{2}}\approx 5\times {10}^{4}\,{\left(\displaystyle \frac{\delta }{40}\right)}^{-1}\,\left(\displaystyle \frac{{E}_{\mathrm{ssc}}}{1\ \mathrm{TeV}}\right)\end{eqnarray} \tag{ 9 }$

where E_ssc is the observed energy of the peak of the SSC component. Here, we use a value of δ consistent with the lower limit from $\gamma \gamma$ pair production (Section 6.1). Here and from this point forward, we neglect factors of 1 + z, which will be quite small given the redshift of Mrk 421 (z = 0.03). To accelerate electrons to this peak requires

$\begin{eqnarray}&&{B}^{{\prime} }\approx 9\,{\left(\displaystyle \frac{\delta }{40}\right)}^{2}{\left(\displaystyle \frac{{N}_{a}}{100}\right)}^{-1}\ {\rm{G}}.\end{eqnarray} \tag{ 10 }$

Large values of the magnetic field or large values of N_a are required to accelerate electrons to ${\gamma }_{b}^{{\prime} }$ where ${t}_{\mathrm{acc}}^{{\prime} }={t}_{\mathrm{syn}}^{{\prime} }$. Similar results were found for the luminous, rapid flare from PKS 2155−304 in 2006 (Finke et al. 2008).

Acceleration and cooling may control the light-curve timescales, if ${R}_{b}^{{\prime} }/c\ll {t}_{\mathrm{syn}}^{{\prime} }={t}_{\mathrm{acc}}^{{\prime} }$ or ${R}_{b}^{{\prime} }\ll 1.4\times {10}^{12}{(\delta /40)}^{-3}{({N}_{a}/100)}^{-2}\ \mathrm{cm}$, where ${R}_{b}^{{\prime} }$ is the radius of a spherical emitting region in the co-moving frame. The required upper limit on ${R}_{b}^{{\prime} }$ is two to three orders of magnitude smaller than typical values for the size of the emitting region (e.g., Abdo et al. 2011), even for flaring episodes (e.g., Aleksić et al. 2015b). Additionally, if the flare's rise and decay times were dominated by the acceleration and cooling timescales, one would expect the flaring timescales to be energy dependent; however, they seem to have the same timescale across energy bands (see Table 2). Therefore, the timescales of the fast component of the light curve are likely controlled by the light-crossing time of a blob with a fixed size. As a result, they should appear symmetric and with timescales independent of the energy.

6.3.2. Plasmoids in Magnetic Reconnection

Magnetic reconnection is invoked as an efficient particle acceleration process in a variety of astrophysical sources of nonthermal high-energy radiation, including AGN jets (Romanova & Lovelace 1992; Giannios et al. 2009, 2010; Giannios 2013). It has been proposed that plasmoids (i.e., blobs of magnetized plasma containing energetic particles) that are formed and accelerated in the reconnection regions of jets can serve as high-energy emission sites in both blazars and radio galaxies (Giannios et al. 2009; Sironi et al. 2015). Petropoulou et al. (2016, hereafter PGS16) presented a semi-analytic model of flares powered by plasmoids in a reconnection layer, simplifying the results of detailed particle-in-cell (PIC) simulations (for a full numerical treatment, see Christie et al. 2019).

A single plasmoid produces a flare with a peak luminosity and flux-doubling timescale that depend on its size and Doppler beaming. A unique feature of this model is that the flux-doubling timescale in the rising part of the flare is mostly determined by the acceleration of the plasmoid in the layer (i.e., by its bulk motion). As a result, similar rise timescales should be observed at different energy bands of a flare powered by a single plasmoid (achromatic behavior). In contrast to the rising part of a flare, its decay is not constrained by PIC simulations. By setting the decay timescale to be approximately equal to the rise timescale of the fast flares, as observed, one can infer the declining rate at which accelerated particles are injected in the plasmoid or the decay rate of the magnetic field after the plasmoid has left the layer (PGS16).

No strong spectral evolution is expected during a flare produced by a single plasmoid (PGS16; Christie et al. 2019). At any given time though, an observer receives radiation from a large number of plasmoids in the layer, having different sizes and Doppler factors. Those plasmoids that move with mildly relativistic speeds (in the jet frame) and have intermediate sizes (in terms of the layer's size) can contribute and even dominate the overall emission. The superposition of their emission could result in a slowly varying and more luminous component of the light curve (Giannios 2013; Christie et al. 2019), which may exhibit spectral variations and drive the fractional variability on longer timescales than the fast varying component of the light curve (Christie et al. 2020).

Energetics and timescales—PGS16 provided simple formulae to estimate the flux-doubling timescale (t_1/2) and peak luminosity (L_pk) of flares produced by individual plasmoids in a reconnection layer (with half-length L') ⁷⁴ of a blazar jet (for an illustration, see Figure 3 in Christie et al. 2019). The bolometric peak flare luminosity can be written as

$\begin{eqnarray}&&{L}_{\mathrm{pk},\mathrm{bol}}=\displaystyle \frac{\pi }{2}{\beta }_{{\rm{g}}}{{cw}}_{{\rm{f}}}^{{\prime} 2}{\delta }_{{\rm{p}},{\rm{f}}}^{4}{u}_{{\rm{e}}}^{{\prime} }\end{eqnarray} \tag{ 11 }$

where β_g is the growth rate of a plasmoid, ${w}_{{\rm{f}}}^{{\prime} }$ and ${\delta }_{{\rm{p}},{\rm{f}}}$ are, respectively, the plasmoid transverse size and the Doppler factor at the end of its lifetime (i.e., when it is being advected from the layer or when it merges with another bigger plasmoid). The Doppler factor takes into account the relativistic motion of the plasmoid in the layer and the relativistic motion of the layer itself (for a definition, see Equation (8) in PGS16). In the above equation, ${u}_{{\rm{e}}}^{{\prime} }\simeq {f}_{\mathrm{rec}}{L}_{{\rm{j}}}/4\pi {\varpi }^{2}c{\beta }_{{\rm{j}}}{{\rm{\Gamma }}}_{{\rm{j}}}^{2}$, where Γ_j is the jet's bulk Lorentz factor, β_j  ∼ 1, L_j is the absolute power of a two-sided jet, $\varpi /{L}^{{\prime} }={\epsilon }_{\mathrm{rec}}/{{\rm{\Gamma }}}_{{\rm{j}}}{\theta }_{{\rm{j}}}$ is the cross-sectional radius of the jet, _rec ≃ 0.15 is the reconnection rate, and f_rec ≃ 0.5 is the fraction of dissipated magnetic energy transferred to relativistic pairs (e.g., Sironi et al. 2015). The flux-doubling timescale can be estimated as

$\begin{eqnarray}&&{t}_{1/2}\approx \displaystyle \frac{1}{c{\beta }_{{\rm{g}}}}{\int }_{{w}_{1/2}^{{\prime} }}^{{w}_{{\rm{f}}}^{{\prime} }}\displaystyle \frac{d\tilde{w}}{{\delta }_{{\rm{p}}}(\tilde{w})}\end{eqnarray} \tag{ 12 }$

where δ_p is the plasmoid's Doppler factor (see Equation (8) in PGS16), which evolves during the plasmoid lifetime as this accelerates in the layer. Here, ${w}_{1/2}^{{\prime} }$ is the size of the plasmoid at the moment the flare luminosity reaches half of its peak value. For details about the derivation of Equations (11)–(12) and assumptions therein, we refer the reader to PGS16. The free parameters of the model are: the plasma magnetization σ, the orientation of the layer with respect to the jet axis θ', the observer's angle θ_obs, L', L_j, and Γ_j.

We apply the PGS16 model to the fast flaring activity observed in Mrk 421. We use as an illustrative example the results for 2013 April 15, where fast flares with similar flux-doubling timescales (Table 2) have been detected in all energy bands (3–7 keV, 7–30 keV, 30–80 keV, 0.2–0.4 TeV, 0.4–0.8 TeV, and > 0.8 TeV). Using the peak flux of the fast flare in each energy band, we estimate the peak bolometric luminosity of the fast flare. Given the similar rise timescales in all bands, we use their mean value as the observed t_1/2 of the fast flare (see the black symbol in Figure 8).

Zoom In Zoom Out Reset image size

We consider a case with plasma magnetization σ = 10. Other than the microphysical parameters, which are benchmarked with PIC simulations of reconnection (see Table 1 in PGS16), we adopt the following values for the free parameters of the model: L' = 4 × 10¹⁵ cm, L_j = 2 × 10⁴⁶ erg s⁻¹, Γ_j = 14, and ${\theta }_{\mathrm{obs}}={\left(2{{\rm{\Gamma }}}_{j}\right)}^{-1}\simeq 2\buildrel{\circ}\over{.} 1$. Using Equations (11) and (12), we can then estimate the peak flux and flux-doubling timescale of flares produced by plasmoids with different sizes, as shown in Figure 8. We note that the curves showing different orientations of the reconnecting layer do not depend monotonically on θ' due to the kinematics of plasmoids in the jet (for details, see PGS16). For the adopted parameter values, we find that a larger range of θ' values results in fluxes and flux-doubling timescales that are compatible with the observed values. A choice of larger L' would just shift the curves horizontally to longer timescales, while an increase of Γ_j would shift the curves diagonally toward the upper-left corner of the plot (for scalings, see Equations (33)–(37) in PSG16). A higher σ value would have a similar effect as that of a higher Γ_j, as it would result in stronger relativistic motions of the plasmoids in the layer. The predicted peak flux depends strongly on the angle at shorter variability timescales, whereas it is almost independent of the orientation of the layer at longer timescales. This merely reflects the fact that flares with longer timescales are produced by the largest plasmoids in the layer that move with nonrelativistic speeds in the jet frame. For these plasmoids, the beaming of the radiation is basically determined by θ_obs and Γ_j . On the contrary, flares with short durations are produced by plasmoids that move relativistically in the jet frame and whose Doppler factor δ_p depends sensitively on both angles θ' and θ_obs.

Flux–flux correlations—Here, we examine the flux–flux correlations predicted by the model for the fast and slow components of a plasmoid-powered light curve. To do so, we adopt the results of Christie et al. (2019) and, in particular, the light curves computed for a "vanilla" model ⁷⁵ of a BL Lac–like source for σ = 10 (BL10) with θ' = 30° and θ_obs = 0°. The values for the angles correspond to the vanilla BL Lac model of Christie et al. (2019). A choice of 21 instead of 0° would have no effect on the conclusions. We focus on a 4 hr–long segment of the total light curve that displays a fast flare (with duration ∼0.3 hr) emerging on top of a less variable component (see inset plot in Figure 9). The fluxes computed in the 3–7 keV and >0.8 TeV energy bands during the selected time window are displayed in the main panel of Figure 9 (blue points). The fluxes computed during the fast flare are highlighted for clarity (orange colored points). Although our vanilla model cannot explain the details of the observed correlations shown in Figure 7, the predicted flux–flux correlations bear some similarities with the observed ones: a tight correlation with a slope close to 1 is produced by the slow component of the light curve, whereas the fast flare leads to a looser correlation with a steeper slope.

Zoom In Zoom Out Reset image size

We note that the temporal resolution of our model light curve is about 30 s. The loop-like structure produced by the fast component would likely be missed in real data due to the coarser sampling. For example, even with the 15 minute sampling used in this work, the loop would consist of just two data points. Our results are also in agreement with those shown in Figure 7, where the loop is not evident in the real data but becomes visible only when the fitted fluxes of the slow and fast components are plotted.

Caveats—We next discuss some caveats of the plasmoid-dominated reconnection model. The fast component of the optical light curve is predicted to be wider than the one seen in the X-rays and VHE gamma-rays, due to the longer cooling timescales of electrons radiating in the optical band. Nevertheless, some degree of correlation between the optical and other high-energy bands is expected in this model. Unfortunately, this cannot be tested with our current data set because we do not have proper optical coverage around the peak times of the fast (sub-hour) flares in the X-rays and VHE gamma-rays (see Figures 1 and 4, C2–C4). On longer (multi-hour) timescales, the observed optical flux appears to roughly follow the temporal trends in X-rays and VHE gamma-rays for some days (e.g., April 13 and 17), while it appears to be anticorrelated for some other days (e.g., April 11, 12, 15, and 18). Such variability patterns cannot be explained by the vanilla model we described above. If magnetic reconnection can be triggered in different jet locations, then it is possible to have the formation of multiple layers in the jet with different properties (e.g., sizes, magnetizations, and orientations; Giannios et al. 2009; Giannios & Uzdensky 2019). It is then possible that the observed optical emission is dominated by a different layer (e.g., larger and with different orientation) than the one producing the high-energy emissions (e.g., smaller and more aligned to the observer).

To fully account for all of the observed properties of the flares (i.e., absolute flux and spectra), one would have to adjust some of the parameters entering the vanilla model, such as the properties of the injected particle distribution within each plasmoid (e.g., slope, and minimum and maximum Lorentz factor). The vanilla model described above could be improved through the inclusion of continuous slow acceleration of particles within plasmoids due to plasmoid compression (Petropoulou & Sironi 2018) and/or impulsive particle injection with a harder energy spectrum during plasmoid mergers (Christie et al. 2019).

PGS16 and Christie et al. (2019) developed a radiative model based on the PIC simulations of Sironi et al. (2016). Although the radiative model takes into account the plasmoid kinematics and dynamics, it treats each plasmoid as a homogeneous source. Thus, this model cannot be used to predict the polarization signatures of a reconnection event. This issue has been recently tackled in Zhang et al. (2018), where the polarized emission was computed directly from 2D PIC simulations of reconnection, demonstrating that the latter can produce variable optical polarization. Due to the different simulation setups used in Sironi et al. (2016) and Zhang et al. (2018), it is not clear how the polarization results found in the latter study can be applied to the plasmoid-dominated model discussed here.

This manuscript presents results from an extensive multiwavelength campaign in 2013 on the blazar Mrk 421, which included, for the first time, highly sensitive hard X-ray observations with NuSTAR. This study focuses on the exceptional flaring activity observed during nine consecutive days, from April 11 to 19, and complements the results obtained with the low VHE gamma-ray activity from 2013 January to March, which was reported in Baloković et al. (2016). The 2013 April flaring activity is the brightest detected with MAGIC to date and the second brightest observed in Mrk 421 to date, after the 2010 February flaring activity reported in Abeysekara et al. (2020). More important than the bright blazar activity is the comprehensive multi-instrument data set collected during these nine days, which includes 42 hr with MAGIC, 27 hr with VERITAS, 18 hr with Swift, and 71 hr with NuSTAR, out of which 43 hr were taken simultaneously with the VHE observations from MAGIC and VERITAS.

Using simultaneous MAGIC–VERITAS observations of Mrk 421 during this high activity, we noticed that the VHE gamma-ray flux values from VERITAS are systematically lower than those from MAGIC throughout the flaring activity. The flux difference is energy dependent, about 10% in the 0.2–0.4 TeV band and about 30% in the >0.8 TeV band. These differences are within the quoted systematic uncertainties from both instruments and are probably related to the absolute energy reconstruction. In order to perform variability and correlation studies during single days, the VERITAS flux values were scaled up to match those from MAGIC (see Appendix A). The physics results reported in this manuscript relate to the variations in the VHE gamma-ray flux and its correlation to the X-ray flux and, hence, do not depend on the absolute value of the VHE gamma-ray flux. The results reported here would essentially be the same if we had scaled down the MAGIC fluxes to match the VERITAS fluxes.

Owing to the large fluxes and the unprecedented coverage provided by the simultaneous NuSTAR and MAGIC/VERITAS observations during these nine consecutive days, this data set allowed us to evaluate the variability and correlations over three X-ray bands (3–7 keV, 7–30 keV, and 30–80 keV) and three VHE gamma-ray bands (0.2–0.4 TeV, 0.4–0.8 TeV, and >0.8 TeV) on timescales of 15 minutes, producing the most detailed X-ray/VHE variability and correlation study of Mrk 421 to date. This study yielded a number of results, which we summarize as follows:

• 1.  
  The fractional variability F_var versus energy for the nine day flaring activity shows a similar pattern to that during the low X-ray and VHE gamma-ray activity from Baloković et al. (2016). The F_var versus energy shows a double-bump structure with the highest variability occurring in the X-ray and VHE gamma-ray bands, while the variability in the radio and optical bands is very low. Additionally, we find that: (a) the F_var values for the highest X-ray and VHE energies during this nine day flaring activity are much higher than those from the low activity reported in Baloković et al. (2016); (b) the F_var values obtained with data from the single nights are much smaller than those obtained with the nine day flaring activity, indicating that the processes with timescales larger than a half day dominate the variability at X-ray and VHE over the processes that have timescales of hours; and (c) F_var typically increases with energy for single nights, in the same way as observed for multi-day or multi-month timescales; however, we also find nights where F_var does not increase (or even decreases) with energy, hence, suggesting the existence of distinct mechanisms (or distinct particle populations and/or regions) responsible for the variability of Mrk 421.
• 2.  
  There is significant variability in the optical emission, as well as in the polarization degree and EVPA (see Figure 1), but they are not correlated with the X-ray and VHE gamma-ray emissions. The optical polarization variations, while statistically significant, relate to a very low polarization degree (typically less than 5%) and appear to be random and without any obvious coherent structure. This is consistent with a multi-zone scenario such as the one proposed in Marscher (2014), where the sum of the polarization vectors from many zones would result in a low level of polarization with random fluctuations in both the polarization degree and the EVPA. At 86 GHz, the polarization degree is also low (about 3%), as also happens with the optical data, but the polarization angle differs by about 70°, which indicates that the radio and optical emissions are produced, or at least dominated, by different regions along the jet of Mrk 421.
• 3.  
  The flux variability measured at GeV is dominated by the large flux variation around April 15, which is the day with the highest X-ray activity. Despite the lower sensitivity of Fermi-LAT, in comparison to that of the X-ray and VHE gamma-ray instruments, the GeV flux of Mrk 421 can be measured significantly on 12 hr timescales, and it does not show the large flux variations (by factors of a few) that are reported for the X-ray and VHE gamma-ray energies. Therefore, this observation indicates that, when considering timescales of days, there is no correlation between the GeV band and the keV and TeV bands.
• 4.  
  In contrast to the 2013 January–March results reported in Baloković et al. (2016), the light curves from the 2013 April flaring activity show substantial intra-day variability in both the X-ray and VHE gamma-ray bands. Moreover, one can clearly distinguish variability on sub-hour timescales on the top of flux variations occurring on multi-hour timescales. The intra-night flux variations were quantified with a function consisting of an exponential increase and decay (fast component) on the top of a monotonically increasing or decreasing flux (slow component). We found out that, within the X-ray and VHE gamma-ray bands, the parameters describing the fast component do not depend on energy, while those describing the slow component can depend strongly on the energy. We also found out that the fast component is symmetric (rise time ∼decay time). This suggests that the mechanisms that dominate the production of the sub-hour flux variability appear to be "achromatic", while those responsible for the production of the multi-hour flux variability can be "chromatic", at least in some cases.
• 5.  
  A lower limit to the Doppler factor of 35 was derived by requiring that the emitting region producing the fastest VHE variability is optically thin to γγ pair production (see Section 6.1). This limit implicitly assumes that the low-energy (optical/IR) and VHE gamma-ray photons are produced in the same region. However, the multiband variability and correlations derived with the nine day 2013 April flaring activity indicate that these two emissions likely originated in different regions, which would relax the derived lower limit.
• 6.  
  Using the parameters from the fast component, we find that there are no delays between the X-ray and VHE gamma-ray emissions down to the resolution of our measurement (3 minutes during April 15). This is the strongest constraint on the correlated behavior between the X-rays and gamma-rays in Mrk 421 and among the strongest constraints derived with TeV blazar data until now.
• 7.  
  The correlation between VHE gamma-rays and X-rays is positive (i.e., the VHE gamma-ray emission increases when the X-ray emission increases); but there are multiple flavors in the strength and characteristics of this correlated behavior that change both across energy (even for nearby energy bands) and over time (on day timescales). The strongest correlation occurs between the lowest X-ray band (3–7 keV) and the highest VHE gamma-ray band (>0.8 TeV), where one finds an approximately linear change in the VHE flux with the X-ray flux. On the other hand, the weakest correlation occurs between the highest X-ray band (30–80 keV) and the lowest VHE gamma-ray band (0.2–0.4 TeV), where the VHE flux changes with the X-ray flux with a slope of ∼0.3 in log–log scale. This indicates that the particle population dominating the emission in the 3–7 keV and >0.8 TeV bands are closely related, while this does not occur for the bands 30–80 keV and 0.2–0.4 TeV. The decrease in the magnitude and significance of the correlation (with respect to the maximum) when increasing the X-ray energy and decreasing the VHE gamma-ray energy are expected if the variability is dominated by the highest-energy electrons; but the rapid change in the correlation pattern with energy is remarkable, and it has never been observed to date neither for Mrk 421 nor for any other TeV blazar.
• 8.  
  Within the context of a synchrotron/SSC model, the approximately linear correlation between VHE and X-ray could be an indication of an expanding blob with a decreasing magnetic field (Katarzyński et al. 2005). It could also be an indication of more than one emitting region (Fossati et al. 2008).
• 9.  
  The temporal and spectral properties of the multiband flares disfavor a single-zone interpretation of the results (see Section 6). A scenario with multiple zones (and possibly with narrow electron energy distributions) is likely needed to explain the achromatic (chromatic) variability of the fast (slow) component of the light curves, as well as the changes of the flux–flux correlation within day-long timescales. We showed that plasmoids forming in the jet of Mrk 421 due to magnetic reconnection might explain some of the main observational results of this campaign. In this scenario, the multi-hour flux variations, which are found to be chromatic for some nights, would be dominated by the combined emission from various plasmoids of different sizes and velocities, while the sub-hour flux variations, which appear to be achromatic, would be dominated by the emission from a single small plasmoid moving with relativistic speeds along the magnetic reconnection layer. Due to the different origins of the fast and slow components of plasmoid-powered light curves, the flux correlations between different energy bands are also expected to differ. In particular, the reconnection model predicts tight correlations (with linear or sub-linear relations between energy bands) for the slow component of the light curve and weaker correlations (with steeper-than-linear relations between fluxes) during fast flares. However, the characteristic loop-like structure in the VHE flux versus X-ray flux plot predicted for fast flares is likely to be missed with the realistic temporal binning.

The accuracy and level of detail in this study, which are unmatched among all VHE gamma-ray emitting blazars (including past observations of Mrk 421), show a large degree of complexity in the variability and correlation patterns. This complexity may be present in other blazars, but it may be difficult to observe owing to insufficient temporal and energy coverage in the observations, especially with the instruments that can resolve the VHE gamma-ray fluxes with high precision. The study presented here for Mrk 421 sheds some light onto this complex behavior and represents a pathfinder to the studies that may be done with the next generation of ground-based VHE gamma-ray instruments like the Cherenkov Telescope Array (CTA), which is expected to resolve many VHE blazars with a level of accuracy comparable to that of MAGIC and VERITAS to resolve Mrk 421 during this period of outstanding activity.

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 and MPG, the Italian INFN and INAF, the Swiss National Fund SNF, the ERDF under the Spanish MINECO (FPA2017-87859-P, FPA2017-85668-P, FPA2017-82729-C6-2-R, FPA2017-82729-C6-6-R, FPA2017-82729-C6-5-R, AYA2015-71042-P, AYA2016-76012-C3-1-P, ESP2017-87055-C2-2-P, and FPA201790566REDC), the Indian Department of Atomic Energy, the Japanese JSPS and MEXT, the Bulgarian Ministry of Education and Science, National RI Roadmap Project DO1-153/28.08.2018 and the Academy of Finland grant No. 320045 is gratefully acknowledged. This work was also supported by the Spanish Centro de Excelencia "Severo Ochoa" SEV-2016-0588 and SEV-2015-0548, and Unidad de Excelencia "María de Maeztu" MDM-2014-0369, by the Croatian Science Foundation (HrZZ) Project IP-2016-06-9782 and the University of Rijeka Project 13.12.1.3.02, by the DFG Collaborative Research Centers SFB823/C4 and SFB876/C3, the Polish National Research Centre grant UMO-2016/22/M/ST9/00382, 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 made use of data from the NuSTAR mission, a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory, and funded by the National Aeronautics and Space Administration. We thank the NuSTAR Operations, Software, and Calibration teams for support with the execution and analysis of these observations. This research has made use of the NuSTAR Data Analysis Software (NuSTARDAS) jointly developed by the ASI Science Data Center (ASDC; Italy) and the California Institute of Technology (USA).

This research has made use of the XRT Data Analysis Software (XRTDAS) developed under the responsibility of the ASI Science Data Center (ASDC), Italy.

D.P. and A.B. are grateful to Amy Furniss and Wystan Benbow for providing the VERITAS VHE gamma-ray fluxes and for useful discussions about them. M.P. acknowledges support from the Lyman Jr. Spitzer Postdoctoral Fellowship and NASA Fermi grant No. 80NSSC18K1745., and J.F. is partially supported by NASA under contract S-15633Y. M.B. acknowledges support from NASA Headquarters under the NASA Earth and Space Science Fellowship Program (grant NNX14AQ07H) and the black hole Initiative at Harvard University, which is funded in part by the Gordon and Betty Moore Foundation (grant GBMF8273) and in part by the John Templeton Foundation.

This research was partially supported by the Bulgarian National Science Fund of the Ministry of Education and Science under grants DN 18-13/2017, DN 18-10/2017, KP-06-H28/3 (2018), and KP-06-PN38/1 (2019), as well as for the Spanish MIMECO (AYA2016-80889-P, RYC-2013-14511) and the IAA-CSIC "Severo Ochoa" program SEV-2017-0709. The St. Petersburg University team acknowledges support from Russian Science Foundation grant 17-12-01029. The Abastumani team acknowledges financial support of the project FR/638/6-320/12 by the Shota Rustaveli National Science Foundation under contract 31/77. T.G. acknowledges support from Istanbul University (Project numbers 49429 and 48285), Bilim Akademisi (BAGEP program), and TUBITAK (project numbers 13AT100-431, 13AT100-466, and 13AT60-430). The Boston University effort was supported in part by NASA grants NNX12AO90G and NNX14AQ58G. Data from the Steward Observatory spectropolarimetric monitoring project were used in this paper. This program is supported by Fermi Guest Investigator grants NNX08AW56G, NNX09AU10G, NNX12AO93G, and NNX15AU81G. The OVRO 40 m monitoring program is supported in part by NASA grants NNX08AW31G and NNX11A043G and NSF grants AST-0808050 and AST-1109911. The Metsähovi team acknowledges the support from the Academy of Finland to our observing projects (numbers 212656, 210338, 121148, and others). W.M. acknowledges support from CONICYT project Basal AFB-170002. I.A. acknowledges financial support from the grant IAA4SKA (R18-RT-3082) funded by the Economic Transformation, Industry, Knowledge and Universities Council of the Regional Government of Andalusia and the European Regional Development Fund from the European Union.

The high accuracy in the VHE gamma-ray flux measurements presented in this work reveals an offset in the simultaneous MAGIC and VERITAS fluxes. This offset, which is energy dependent and can reach up to 30%, is within the flux systematic uncertainties reported by each instrument (on the order of 20%–25%). The factor that dominates the systematic uncertainties in an IACT is the uncertainty of the absolute energy-scale calibration. The energy reconstruction of VHE gamma-ray showers by IACTs relies heavily on Monte Carlo simulations. The main source of systematic uncertainty on the energy scale is the uncertainty on the "light yield" (i.e., the total light throughput) considered in these simulations, associated for instance, with the average transparency of the atmosphere or the light collection efficiency of the telescopes. Each instrument tunes these parameters over long periods of time with stable performance (during years). For this reason, deviations with respect to the considered atmosphere and telescope models are expected for shorter timescales on the order of 15% (as discussed in Madhavan & VERITAS Collaboration 2013; Aleksić et al. 2016).

In order to exploit the excellent coverage of this data set and perform intra-night VHE gamma-ray flux variability and X-ray/VHE flux correlation studies, these energy-dependent offsets need to be corrected. This mismatch shows up prominently when comparing single-night light curves. As an example, the left panel in Figure A1 displays the multiband VHE light curves during the night of April 12, showing an increasing offset with energy reaching about 30% for the VHE gamma-ray fluxes with E > 0.8 TeV. We note that an increasing offset with energy is expected if the main source of this systematic difference between the MAGIC and VERITAS fluxes is the energy scale.

Zoom In Zoom Out Reset image size

The simplest approach to correct for this difference is to scale the light curves of one of the instruments by an energy-dependent factor. As we have been using, throughout this work, three different energy ranges (i.e., 0.2–0.4 TeV, 0.4–0.8 TeV, and >0.8 TeV), three different flux ratios need to be calculated to properly correct for this bias. To calculate these correction factors, two different methods were tested: (a) obtain three scaling factors by normalizing the MAGIC and VERITAS fluxes in the three energy bands; and (b) calculate the energy-scale correction that makes MAGIC and VERITAS spectra compatible and then obtain the three scaling factors for the three energy bands considered. The first approach is purely agnostic, while the second one assumes that the small mismatch between MAGIC and VERITAS VHE gamma-ray measured fluxes is dominated by the systematic uncertainty on the energy scale of each instrument. Both strategies make use of the 2 hr 45 minutes of simultaneous MAGIC and VERITAS observations taken during four different nights (see Table 1).

Inferring the correction factor from available simultaneous data is the most direct method to calculate these factors; however, it cannot be applied to the 0.2–0.4 TeV energy band because the simultaneous MAGIC–VERITAS data relate to MAGIC observations at a zenith angle above 55°, which have an analysis energy threshold above 0.3 TeV (see Section 2). Figure A2 shows a comparison of the VHE gamma-ray fluxes (in the two available energy bands) derived with simultaneous 15 minute time bins, which span throughout four consecutive nights: from April 12 to 15.

Zoom In Zoom Out Reset image size

The figure shows two important characteristics of this offset. First, the ratio at a given energy band is approximately the same for all of the simultaneous observations, hence, indicating a systematic effect between MAGIC and VERITAS, which is not related to the peculiarity of one single day (e.g., bad weather in one of the two telescope sites). Second, the ratio of MAGIC and VERITAS fluxes deviates from 1 by a larger amount in the highest-energy band, hence, confirming the energy dependence of the offset. The correction factors and uncertainties inferred with this method are reported in the first column of Table A1.

As described in Aleksić et al. (2016), the MAGIC analysis allows us to test the impact of a modified light yield over a given reconstructed SED. By applying relative light yield corrections to the MAGIC data that were observed simultaneously with those of VERITAS, we are able to test which value minimizes the MAGIC–VERITAS SED differences. This search determined that a decrease of 20% in the MAGIC light makes the simultaneous MAGIC and VERITAS spectra compatible. ⁷⁶ The flux correction factors for the three energy ranges considered here were then computed through the ratio of the pre- and post-light-yield-corrected MAGIC light curves in the three energy bands. In order to reduce the statistical uncertainties, we used the full MAGIC data set from these four nights (April 12–15) instead of only using the time intervals with simultaneous MAGIC–VERITAS observations. The flux ratios derived with this method are reported in the second column of Table A1.

As shown in Table A1, there is a very good agreement in the flux correction factors calculated with these two different methodologies. Hereafter, the VERITAS light curves are normalized to those of MAGIC using the method derived with the energy-scale correction test (second column in Table A1). We chose this method (instead of applying the flux ratio method) because this one provides a scaling factor for the lowest-energy band, 0.2–0.4 TeV, for which we do not have simultaneous MAGIC–VERITAS data. An example of a light curve after the application of the scaling factors to the VERITAS fluxes is given in the right panel of Figure A1.

In order to test the impact of the scaling factor uncertainty on the VHE/X-ray correlations reported in Section 5 and Tables 3 and D1, we scaled the VERITAS fluxes by a value drawn from a normal distribution with mean and sigma given by the parameters in Table A1. This was repeated 1000 times, and each time, all of the correlations were calculated. The standard deviation of the obtained set of correlation parameters was taken as a measure of the systematic uncertainty associated with the scaling of the VERITAS fluxes. For the nine day correlations in Table 3, the systematic uncertainty was found to be an order of magnitude smaller than the statistical uncertainty. This is expected because the flux corrections are of the order of 20%, while the nine day flaring activity is dominated by flux changes of about one order of magnitude. On the other hand, the impact of this flux scaling is not negligible for quantities derived with the single-night light curves, which are dominated by flux changes of about a factor of two. The uncertainties related to the flux scaling are reported as the second uncertainty in Table D1. Typically, these systematic uncertainties are found to be smaller, or at most comparable, to the statistical uncertainty.

Using the scaling factors from Appendix A, one can now up-scale all of the VERITAS flux measurements and produce single-night light curves with simultaneous VHE and X-ray data spanning about 8–10 hr. The VHE flux measurements can be best compared with the NuSTAR measurements on timescales of 15 minutes, through light curves where the X-ray and VHE gamma-ray fluxes are normalized so that one can compare the relative differences among them.

In order to enable direct (normalized) flux comparisons between the X-ray and VHE bands, the fluxes are normalized by calculating the average flux for each night and each band using only data points that have simultaneous VHE and X-ray observations (filled markers in Figure B1). The resulting normalization factors are reported in Table B1. Subsequently, all of the light-curve fluxes (including the nonsimultaneous data marked with open symbols) are divided by this normalization factor. To be able to visually compare the intra-night flux changes, all VHE and X-ray band combinations for each day are shown in Figure B1. Using the normalization factors from Table B1, the scaling factors from Table A1, and the normalized light-curve flux points from Figure B1, one can retrieve the primal X-ray and VHE light curves from NuSTAR, MAGIC, and VERITAS.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The normalized light curves in Figure B1 all show that, during this few-days-long activity, the X-ray and VHE emissions of Mrk 421 show quite a number of structures on timescales from multiple hours and down to timescales smaller than 1 hr. They also show remarkable correlations in some of the band combinations for selected nights (e.g., the 3–7 keV and >0.8 TeV bands on April 15, or the 3–7 keV and 0.2–0.4 TeV bands on April 11). Looking at single-day band combinations plots, one can also explicitly see different trends for the fast and slow components (e.g., on April 15, the slow change in the 3–7 keV flux is not as perfectly matched in the 0.2–0.4 TeV band as it is with the highest VHE band, while the fast component shows comparable agreement).

This section reports the quantification of the X-ray and VHE gamma-ray multiband flux variability during single nights using the templates defined by Equation (3) in Section 4.2. The best-fit results obtained with the light curves from April 11–14 are reported in Table C1, while the results obtained for April 16 are reported in Table C2. Below, we provide a brief description of the used strategy and the obtained results for each of the nights.

Table C2. Parameters Resulting from the Fit with Equation (C1) to the X-Ray and VHE Multiband Light Curves from 2013 April 16

Band	Flare Amplitude Aa	${\chi }^{2}/{\rm{d}}.{\rm{o}}.{\rm{f}}.$
3–7 keV	0.44 ± 0.02	⋯
7–30 keV	0.38 ± 0.02	⋯
30–80 keV	0.083 ± 0.007	⋯
0.2–0.4 TeV	4.2 ± 0.2	140/26
0.4–0.8 TeV	1.57 ± 0.05	55/26
>0.8 TeV	0.72 ± 0.02	39/26

Note. The three X-ray bands are fit with a single function, yielding t_rise =1.4 ± 0.2 hr, t_fall = 5.0 ± 0.5 hr, break time t₀ = −0.42 ± 0.34 hr from midnight, with ${\chi }^{2}/{\rm{d}}.{\rm{o}}.{\rm{f}}.\,=\,582/57$.

^a The flare amplitude A is given in in 10⁻⁹ erg cm⁻² s⁻¹ for the X-ray bands, and 10⁻¹⁰ ph cm⁻² s⁻¹ for the VHE gamma-ray bands.

Download table as:  ASCIITypeset image

2013 April 11. There are no flares visible in any of the VHE or X-ray bands. The light curves are not consistent with constant flux and are fit with a linear model, i.e., the slow component. The results are shown in Figure C1. For all of the X-ray and VHE gamma-ray bands from these 10 hr–long light curves, the slopes are all consistent within the error bars.

Zoom In Zoom Out Reset image size

2013 April 12. There is a clear flare visible in all X-ray bands around 2 hr past midnight; but unfortunately, there is a gap in the VHE coverage for that time interval. Therefore, VHE light curves are fit with a linear function (slow component) only. The results are shown in Figure C2. All fits are performed only up to 7.7 hr past midnight. There is flickering in VHE bands, with light curves inconsistent with constant flux, and the simple linear model or Slow+Fast component model from Section 4.2 does not describe data well either. Fits with slow component give low slope values, consistent with zero for bands above 0.4 TeV.

Zoom In Zoom Out Reset image size

Slopes of the slow component decrease with increasing energy in X-ray bands as well, while flare amplitude increases with energy.

2013 April 13. There is a flare visible in the two higher VHE bands, with the amplitude and the width of the flare both approximately halved in the 0.4–0.8 TeV band w.r.t. the >0.8 TeV band. In the lowest VHE band, there are flux measurements above the slow component at the location of the flare, but no satisfactory Slow+Fast model fit (Section 4.2) could be obtained. There is a ∼2 hr gap between the MAGIC and VERITAS data, and the fit is performed only with MAGIC data point, denoted by using gray instead of green color for VERITAS data in Figure C3.

Zoom In Zoom Out Reset image size

There are no X-ray data covering the time of the VHE flare. On the other hand, there is a small-amplitude flare around 2.5 hr earlier visible in the X-ray bands, without obvious counterpart in the VHE bands. Being a relatively weak mini-flare, with the highest flux measured at the same time in all three X-ray bands, the fit is performed with Slow+Fast model but with flare time (location) fixed to the time of the highest flux measurement, at midnight—0.24 hr.

2013 April 14. There is no obvious flare component in the VHE bands; therefore, VHE bands are fit with slow components only (Figure C4). In the X-ray bands, the flux goes down, then increases in a broad flare visible in two lower X-ray bands, while no such behavior is seen in VHE bands. Flare amplitudes in X-ray bands are consistent within the error bars, but there are some data missing during the time covered by the fast component; therefore, no definite conclusions about the flaring component can be made.

Zoom In Zoom Out Reset image size

The second part of the light curves (about 1 hour after midnight) has a more simple structure, with a significant decrease in flux. The slow component slopes in all bands are similar.

2013 April 16. There are no characteristic mini-flares in the 10 hr–long light curve from April 16. Instead, there is a prominent rise and then fall of light curves from all energy bands. This behavior cannot be fit by the Slow+Fast model from Section 4.2. Because of that, we modify the model to exclude the slow component and only have a broad "bump" with different rise and fall times,

$\begin{eqnarray}&&{F}_{B}(t)=\displaystyle \frac{2}{{2}^{-\tfrac{t-{t}_{0}}{{t}_{\mathrm{rise}}}}+{2}^{\tfrac{t-{t}_{0}}{{t}_{\mathrm{fall}}}}}\cdot A.\end{eqnarray} \tag{ C1 }$

The fit is performed to all three X-ray bands simultaneously, giving a unique rise time t_rise = 1.38 ± 0.23 hr, fall time t_fall = 4.96 ± 0.47 hr, and break time t₀ = −0.42 ± 0.34 hr from midnight, with ${\chi }^{2}/{\rm{d}}.{\rm{o}}.{\rm{f}}.\,=\,582/57$. Only flare amplitude is allowed to vary in different bands, with results reported in Table C2. After fitting to X-ray data, the function shape is fit to each VHE band, again allowing only the flare amplitude to vary. The results (Figure C5), along with the χ² values for each band are reported in Table C2.

Zoom In Zoom Out Reset image size

In this section, we report the quantification of the VHE versus X-ray correlations for single-night data. In particular, we report the Pearson correlation coefficients and related significances, the values of the DCF, and the slopes and the χ² values from the linear fits to the $\mathrm{log}{F}_{\mathrm{VHE}}\mbox{--}\mathrm{log}{F}_{{\rm{X}} \mbox{-} \mathrm{ray}}$. The numbers in Table D1 are complementary to those reported in Table 3, which report the same information for the nine day data set (instead of single-night data sets).

Table D1. Single-night Correlation Coefficients and Fit Slopes for 2019 April 11–16

VHE Band	X-Ray Band	Pearson coeff.	Nσ (Pearson)	DCF	Linear Fit Slope	χ2/d.o.f.
2013 Apr 11
0.2–0.4 TeV	3–7 keV	${0.83}_{-0.07}^{+0.05}\pm 0.02$	6.1	0.90 ± 0.21 ± 0.007	0.86 ± 0.09 ± 0.05	43/28
	7–30 keV	${0.84}_{-0.07}^{+0.05}\pm 0.01$	6.4	0.92 ± 0.22 ± 0.003	0.83 ± 0.08 ± 0.05	37/28
	30–80 keV	${0.74}_{-0.10}^{+0.08}\pm 0.002$	4.9	1.35 ± 0.37 ± 0.02	1.02 ± 0.15 ± 0.07	28/28
0.4–0.8 TeV	3–7 keV	${0.89}_{-0.05}^{+0.03}\pm 0.04$	7.3	1.06 ± 0.26 ± 0.02	1.01 ± 0.09 ± 0.11	26/28
	7–30 keV	${0.86}_{-0.06}^{+0.04}\pm 0.03$	6.7	1.03 ± 0.26 ± 0.03	0.95 ± 0.09 ± 0.10	30/28
	30–80 keV	${0.64}_{-0.13}^{+0.10}\pm 0.01$	4.0	1.29 ± 0.45 ± 0.06	1.08 ± 0.15 ± 0.13	25/28
>0.8 TeV	3–7 keV	${0.79}_{-0.08}^{+0.06}\pm 0.04$	5.6	1.04 ± 0.28 ± 0.01	1.01 ± 0.13 ± 0.10	28/28
	7–30 keV	${0.78}_{-0.09}^{+0.07}\pm 0.04$	5.4	1.03 ± 0.28 ± 0.02	0.96 ± 0.12 ± 0.09	27/28
	30–80 keV	${0.61}_{-0.14}^{+0.11}\pm 0.02$	3.7	1.35 ± 0.47 ± 0.05	1.2 ± 0.2 ± 0.1	25/28
2013 Apr 12
0.2–0.4 TeV	3–7 keV	${0.26}_{-0.19}^{+0.17}\pm 0.10$	1.4	0.31 ± 0.17 ± 0.11	0.21 ± 0.14 ± 0.09	78/27
	7–30 keV	${0.16}_{-0.19}^{+0.18}\pm 0.09$	0.8	0.19 ± 0.19 ± 0.10	0.15 ± 0.09 ± 0.09	82/27
	30–80 keV	$-{0.19}_{-0.18}^{+0.19}\pm 0.07$	1.0	−0.31 ± 0.24 ± 0.12	−1.03 ± 0.39 ± 14	66/27
0.4–0.8 TeV	3–7 keV	${0.05}_{-0.19}^{+0.18}\pm 0.19$	0.3	0.06 ± 0.16 ± 0.23	0.02 ± 0.16 ± 0.18	90/30
	7–30 keV	$-{0.07}_{-0.18}^{+0.19}\pm 0.17$	0.4	−0.09 ± 0.18 ± 0.21	−0.12 ± 0.18 ± 0.18	89/30
	30–80 keV	$-{0.37}_{-0.15}^{+0.17}\pm 0.12$	2.1	−0.65 ± 0.27 ± 0.23	−0.76 ± 0.23 ± 20	59/30
>0.8 TeV	3–7 keV	${0.19}_{-0.18}^{+0.17}\pm 0.19$	1.0	0.27 ± 0.25 ± 0.24	0.15 ± 0.14 ± 0.17	61/30
	7–30 keV	${0.18}_{-0.18}^{+0.17}\pm 0.16$	1.0	0.26 ± 0.27 ± 0.21	0.18 ± 0.16 ± 0.17	61/30
	30–80 keV	$-{0.04}_{-0.18}^{+0.19}\pm 0.09$	0.2	−0.07 ± 0.36 ± 0.17	0.48 ± 0.20 ± 0.79	61/30
2013 Apr 13 a
0.2–0.4 TeV	3–7 keV	${0.74}_{-0.17}^{+0.11}$	3.2	0.73 ± 0.24	0.98 ± 0.24	52/12
	7–30 keV	${0.70}_{-0.19}^{+0.13}$	2.9	0.69 ± 0.24	0.75 ± 0.21	60/12
	30–80 keV	${0.65}_{-0.21}^{+0.14}$	2.6	0.70 ± 0.28	0.58 ± 0.17	49/12
0.4–0.8 TeV	3–7 keV	${0.86}_{-0.11}^{+0.06}$	4.2	0.89 ± 0.25	0.94 ± 0.16	19/12
	7–30 keV	${0.81}_{-0.13}^{+0.08}$	3.7	0.84 ± 0.25	0.72 ± 0.15	24/12
	30–80 keV	${0.76}_{-0.16}^{+0.10}$	3.3	0.85 ± 0.28	0.50 ± 0.12	23/12
>0.8 TeV	3–7 keV	${0.873}_{-0.094}^{+0.056}$	4.5	0.87 ± 0.22	1.40 ± 0.23	28/12
	7–30 keV	${0.81}_{-0.13}^{+0.08}$	3.8	0.81 ± 0.22	1.07 ± 0.23	41/12
	30–80 keV	${0.70}_{-0.19}^{+0.12}$	2.9	0.76 ± 0.24	0.85 ± 0.21	36/12
2013 Apr 14
0.2–0.4 TeV	3–7 keV	${0.79}_{-0.09}^{+0.07}\pm 0.003$	5.2	0.79 ± 0.17 ± 0.002	0.82 ± 0.15 ± 0.05	130/25
	7–30 keV	${0.78}_{-0.09}^{+0.07}\pm 0.003$	5.2	0.79 ± 0.16 ± 0.002	0.87 ± 0.15 ± 0.06	126/25
	30–80 keV	${0.45}_{-0.18}^{+0.15}\pm 0.002$	2.4	0.55 ± 0.20 ± 0.006	1.5 ± 0.4 ± 0.1	76/25
0.4–0.8 TeV	3–7 keV	${0.78}_{-0.09}^{+0.07}\pm 0.006$	5.3	0.79 ± 0.17 ± 0.006	0.79 ± 0.18 ± 0.11	138/27
	7–30 keV	${0.78}_{-0.09}^{+0.07}\pm 0.009$	5.3	0.80 ± 0.17 ± 0.006	0.85 ± 0.19 ± 0.12	138/27
	30–80 keV	${0.41}_{-0.17}^{+0.15}\pm 0.007$	2.2	0.51 ± 0.19 ± 0.002	1.7 ± 0.5 ± 0.3	94/27
>0.8 TeV	3–7 keV	${0.83}_{-0.07}^{+0.05}\pm 0.003$	6.1	0.85 ± 0.21 ± 0.007	0.94 ± 0.18 ± 0.10	107/27
	7–30 keV	${0.81}_{-0.08}^{+0.06}\pm 0.006$	5.7	0.83 ± 0.19 ± 0.003	0.93 ± 0.19 ± 0.11	112/27
	30–80 keV	${0.42}_{-0.17}^{+0.15}\pm 0.007$	2.3	0.51 ± 0.20 ± 0.003	1.6 ± 0.4 ± 0.2	86/27
2013 Apr 15
0.2–0.4 TeV	3–7 keV	${0.77}_{-0.09}^{+0.07}\pm 0.04$	5.3	0.82 ± 0.20 ± 0.04	0.40 ± 0.05 ± 0.03	68/28
	7–30 keV	${0.75}_{-0.10}^{+0.07}\pm 0.05$	5.1	0.81 ± 0.19 ± 0.04	0.33 ± 0.05 ± 0.02	72/28
	30–80 keV	${0.71}_{-0.11}^{+0.08}\pm 0.05$	4.6	0.77 ± 0.18 ± 0.05	0.29 ± 0.05 ± 0.02	80/28
0.4–0.8 TeV	3–7 keV	${0.80}_{-0.08}^{+0.06}\pm 0.08$	5.8	0.86 ± 0.19 ± 0.08	0.52 ± 0.07 ± 0.04	80/29
	7–30 keV	${0.78}_{-0.09}^{+0.06}\pm 0.08$	5.5	0.84 ± 0.19 ± 0.08	0.41 ± 0.06 ± 0.03	88/29
	30–80 keV	${0.73}_{-0.10}^{+0.08}\pm 0.09$	5.0	0.81 ± 0.18 ± 0.09	0.36 ± 0.06 ± 0.03	93/29
>0.8 TeV	3–7 keV	$0.95\pm 0.02\pm 0.02$	9.5	0.94 ± 0.19 ± 0.02	0.91 ± 0.07 ± 0.05	77/29
	7–30 keV	${0.94}_{-0.03}^{+0.02}\pm 0.02$	9.1	0.94 ± 0.19 ± 0.02	0.74 ± 0.06 ± 0.04	87/29
	30–80 keV	${0.88}_{-0.05}^{+0.04}\pm 0.03$	7.3	0.90 ± 0.15 ± 0.02	0.63 ± 0.06 ± 0.03	90/29
2013 Apr 16
0.2–0.4 TeV	3–7 keV	${0.53}_{-0.20}^{+0.15}\pm 0.05$	2.5	0.62 ± 0.22 ± 0.05	0.56 ± 0.16 ± 0.04	51/18
	7–30 keV	${0.42}_{-0.22}^{+0.18}\pm 0.06$	1.8	0.49 ± 0.24 ± 0.07	0.39 ± 0.16 ± 0.05	62/18
	30–80 keV	${0.21}_{-0.24}^{+0.22}\pm 0.08$	0.9	0.47 ± 0.61 ± 0.17	0.77 ± 0.26 ± 0.02	51/18
0.4–0.8 TeV	3–7 keV	${0.78}_{-0.11}^{+0.08}\pm 0.06$	4.3	1.1 ± 0.3 ± 0.07	0.67 ± 0.13 ± 0.09	19/18
	7–30 keV	${0.72}_{-0.14}^{+0.10}\pm 0.09$	3.7	1.0 ± 0.3 ± 0.09	0.55 ± 0.13 ± 0.10	23/18
	30–80 keV	${0.61}_{-0.17}^{+0.13}\pm 0.12$	2.9	1.6 ± 0.5 ± 0.3	0.68 ± 0.15 ± 0.11	18/18
>0.8 TeV	3–7 keV	${0.92}_{-0.05}^{+0.03}\pm 0.02$	6.5	0.97 ± 0.21 ± 0.01	1.6 ± 0.2 ± 0.09	19/18
	7–30 keV	${0.88}_{-0.07}^{+0.04}\pm 0.03$	5.8	0.94 ± 0.23 ± 0.03	1.4 ± 0.2 ± 0.1	27/18
	30–80 keV	${0.64}_{-0.17}^{+0.12}\pm 0.05$	3.1	1.3 ± 0.4 ± 0.1	1.8 ± 0.4 ± 0.06	32/18

Notes. The statistical and systematic contributions to the total uncertainty are reported separately, and in this order.

^a There are no simultaneous NuSTAR and VERITAS data on April 13, and hence there is no systematic error associated to the usage of the flux-scale factors.

Download table as:  ASCIITypeset images: 1 2