© MAGIC Collaboration et al. 2021

1. Introduction

Blazars belong to the group of jetted active galactic nuclei (AGNs) and constitute the most populated class of sources in the extragalactic very-high-energy (VHE; E >  100 GeV) sky¹. They host a super massive black hole (10⁶ − 10⁹ solar masses) surrounded by an accretion disc and display a pair of relativistic plasma jets flowing in opposite directions producing non-thermal radiation. The jet’s axis is oriented at a small angle with the observer’s line of sight, which gives rise to strong relativistic beaming effects of the radiation. Blazars are commonly divided in two families, the flat spectrum radio quasars (FSRQ) and BL Lac objects, depending on their optical spectra (Urry & Padovani 1995). FSRQ exhibit strong emission lines in the optical band, while BL Lac type objects are defined by very weak lines or an absence of such features.

The spectral energy distribution (SED) of BL Lac objects is dominated by non-thermal radiation from the jet and typically shows two continuous components (Ghisellini et al. 2017). Based on polarisation and spectral studies, it is commonly accepted that the low-energy component, peaking in infrared to X-rays, originates from synchrotron radiation by relativistic electrons and/or positrons in a magnetic field. The high-energy component peaks in the GeV-TeV regime and its origin remains under debate. The common scenario involves electron inverse-Compton (IC) scattering off the synchrotron photons emitted by the same population of electrons (see for example Maraschi et al. 1992; Tavecchio et al. 1998; Krawczynski et al. 2004). Such models are labelled as synchrotron self-Compton (SSC) models. In some cases, external Compton models, in which an additional target photon field for IC scattering is introduced, are better suited to describe the SED of BL Lacs (e.g., Madejski et al. 1999; Ghisellini et al. 2005; Böttcher et al. 2013). More complex scenarios invoking hadronic processes can also present a viable explanation (see for example Mannheim 1993; Böttcher et al. 2013; Cerruti et al. 2015). The peak frequency of the low-energy component is frequently used to classify BL Lac type objects in further categories (Padovani & Giommi 1995; Abdo et al. 2010a). Low-frequency BL Lacs (LBL) are defined by a synchrotron peak frequency ν_s <  10¹⁴ Hz and BL Lacs with ν_s >  10¹⁵ Hz are dubbed as high-frequency BL Lacs (HBL). Intermediate-frequency BL Lacs (IBL) have 10¹⁴ Hz < ν_s < 10¹⁵ Hz.

Markarian 421 (Mrk 421; RA = 11^h4′27$\stackrel{″}{.}$31, Dec = 38°12′31$\stackrel{″}{.}$8, J2000) is a HBL at a redshift of 0.031 (de Vaucouleurs et al. 1991). Owing to its brightness and proximity, the source can be well-detected on short timescales (≲1 day) from radio to VHE with current instruments. After its first detection at VHE by the Whipple 10-m Telescope (Punch et al. 1992), numerous multi-wavelength (MWL) campaigns were organised to characterise the broadband emission during individual flares as well as on longer timescales. Similarly to other HBLs, the source shows flux and spectral variability across its entire SED. The variability is most prominent in the X-rays and VHE and was observed down to a sub-hour timescale (Gaidos et al. 1996; Fossati et al. 2008). Using observations over a 14-year time period, Acciari et al. (2014) derived a typical flux above 400 GeV of about 50% that of Crab Nebula flux unit² (C.U.). During low activity, the VHE flux can be as low as 10% C.U., while during strong flaring events, it can reach more than 10 C.U., as observed in an outburst in 2010 (Abeysekara et al. 2020) and 2013 (Acciari et al. 2020).

Several studies reported correlated variability between VHE and X-ray emissions independently from the flux level, which is in agreement with standard SSC models (e.g., Giebels et al. 2007; Fossati et al. 2008; Ahnen et al. 2016; Baloković et al. 2016). The correlation is often parametrised by a simple power law, $F_{\mathrm{V}\mathrm{H}\mathrm{E}}\propto F_{\mathrm{X}-\mathrm{r}\mathrm{a}\mathrm{y}}^x$, where F_VHE and F_{X − ray} are the VHE and X-ray flux, respectively. While the existence of the correlation is well established, the index x may show significant temporal variability: observations revealed sub-linear (x <  1) as well as more-than-quadratic (x >  2) behaviours (Tanihata et al. 2004; Fossati et al. 2008; Baloković et al. 2016). In fact, quadratic or more-than-quadratic trends are typically found during high states. The index x also displays a strong dependency on the exact selection of the spectral bands (Katarzyński et al. 2005; Fossati et al. 2008; Baloković et al. 2016; Acciari et al. 2020). These results underline the complex spectral properties of blazars, but potentially provide insights into the physical processes driving the broadband variability (e.g., Katarzyński et al. 2005). There is currently a lack of detailed investigation of the VHE and X-ray correlation in different energy bands during non-flaring activity, which would provide additional constraints on the emission mechanisms. In order to address this task, sensitive measurements with a dense temporal and wide energy coverage are crucial.

The simplest emission models of blazars (called one-zone models) assume a cospatial particle population responsible for the SED above the infrared (≳10¹³ Hz). In this scenario, a correlation between UV-optical and X-ray photons is generally expected. The observed synchrotron flux is proportional to the product ${\delta }^4{n}_{\mathrm{e}}^{\mathrm{\prime }}{B}^{\mathrm{\prime }2}$, where δ is the Doppler factor, $n_{\mathrm{e}}^{\mathrm{\prime }}$ is the number of electrons, and B′ is the magnetic field (here and in the following, primed quantities refer to quantities in the plasma reference frame). Any change in the latter parameters would simultaneously affect the UV-optical and X-ray emissions. However, only sporadic indications of correlated variability have been reported up to now. Aleksić et al. (2015a) observed a first indication of an anti-correlation during the year 2009. On the other hand, Carnerero et al. (2017) did not report any correlated variability based on long-term observations from 2007 to 2015. It should be noted that these two bands have very different temporal behaviours, the X-ray emission being much more variable, which renders the detection of a correlation challenging.

In this work, we present results from a MWL campaign organised between 2016 December and 2017 June. In order to provide an optimal energy and temporal coverage, we coordinated observations between a large number of instruments covering the emission from radio to VHE. The Florian Goebel Major Atmospheric Gamma Imaging Cherenkov telescopes (MAGIC) and First G-APD Cherenkov Telescope (FACT) carried out observations in the VHE regime. VHE data are complemented by observations from the Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope (Fermi-LAT). Regarding the UV and soft X-ray emission, we organised many observations with the Neil Gehrels Swift Observatory (Swift) satellite in order to obtain a deep temporal coverage. Together with the Swift schedulers team, we coordinated many of these observations to happen simultaneously (or close in time) to the VHE gamma-ray observations performed with MAGIC and FACT in order to be able to properly characterise the temporal evolution of the low- and high-energy SED bumps of Mrk 421. Moreover, we took advantage of four pointings of the Nuclear Spectroscopic Telescope Array (NuSTAR) to obtain a precise characterisation of the hard X-ray (≳10 keV) emission. The four pointings were coordinated to take place strictly simultaneously as the other instruments’ observations. NuSTAR is currently the most sensitive instrument measuring in the hard X-ray regime, corresponding to the high-energy end of the synchrotron spectrum of Mrk 421 (Abdo et al. 2011). Hence, NuSTAR probes the emission from the most energetic particles located in the jet. We investigate in detail the VHE versus X-ray correlation over different spectral bands along the entire MWL campaign. The NuSTAR data significantly widen the energy coverage when combined with Swift data and bring additional constraints to the theoretical model parameters. Furthermore, we study the X-ray versus UV-optical correlated variability.

2. Instruments and analysis

2.1. MAGIC

Mrk 421 belongs to the group of targets that the MAGIC telescopes monitor on a regular basis. The MAGIC telescopes form a system of two 17 m diameter imaging atmospheric Cherenkov telescopes (IACTs). They are situated at an altitude of 2231 m above sea level, on the Canary Island of La Palma at the Roque de los Muchachos Observatory. The integral sensitivity for point sources observations above 220 GeV is (0.66 ± 0.03)% C.U. in 50 h (Aleksić et al. 2016).

The MAGIC dataset presented in this work covers a ≈6-month period from MJD 57727 (2016 December 5) until MJD 57892 (2017 May 19). We analysed the data using the standard analysis tools from the MAGIC Analysis and Reconstruction Software (MARS) package (Zanin et al. 2013; Aleksić et al. 2016). A significant fraction of the observations were carried out with an increased night sky background light contamination (due to the presence of the moon), which directly affects the response of the telescopes (Ahnen et al. 2017a). Due to the varying observing conditions, the data were split into several subsets depending on the level of the moon light contamination. Following Ahnen et al. (2017a), the analysis was then performed by adopting Monte Carlo simulations tuned to match the observing conditions of the different data subsets. The source was observed with zenith angles ranging from 9° to 70°. After data quality selection, ≈70 h of observations were gathered over a total of 48 nights. The light curves were computed in two energy bands: 0.2–1 TeV and > 1 TeV.

The MAGIC spectra were fitted above 100 GeV with a log-parabola spectral model that is defined as follows:

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}E}=f_0{\left(\frac{E}{E_0}\right)}^{-\alpha -\beta log\left(\frac{E}{E_0}\right)}\end{array}$$(1)

where f₀ is the normalisation constant, α the photon index and β the curvature parameter. The normalisation energy E₀ was fixed at 300 GeV. A simple power-law function defined as $\frac{\mathrm{d}N}{\mathrm{d}E}=f_0{\left(\frac{E}{E_0}\right)}^{-\alpha }$ was applied in case a log-parabola was not preferred at a significance above 3σ (based on a likelihood ratio test). The best-fit parameters obtained can be found in Table A.1 together with the fluxes integrated in the two energy bins defined for the light curve. All parameters were evaluated after correcting the spectra for extragalactic background light (EBL) absorption effects using the model of Domínguez et al. (2011).

Four nights had a longer observing time compared to the majority of the observations, which typically lasted between 40 and 60 min. Each of these deep exposures were accompanied by simultaneous NuSTAR and Swift pointings. They took place on MJD 57757, MJD 57785, MJD 57813, and MJD 57840 (2017 January 4, 2017 February 1, 2017 March 1, and 2017 March 28). The MAGIC exposure during these nights varied from 2 h to 6 h depending on the date. For each epoch, light curves with 30-min time bins were produced to study the correlation with the X-ray emission (see Sect. 4).

2.2. FACT

The First G-APD Cherenkov Telescope (FACT) is an imaging atmospheric Cherenkov telescope with a mirror area of 9.5 m² (Anderhub et al. 2013; Biland et al. 2014). It is located close to the MAGIC telescopes, on the Canary Island of La Palma at the Roque de los Muchachos Observatory. FACT measures gamma rays from several hundreds of GeV to about 10 TeV. Operations are performed fully remotely and in an automatic manner. FACT is also pioneering the use of silicon-based photosensors (SiPM aka Geiger-mode Avalanche Photo Diodes or G-APDs), allowing a robust and stable performance (Biland et al. 2014). The observing strategy is specifically tuned to achieve an unbiased long-term monitoring of TeV-emitting blazars.

In this work, Mrk 421 observations between MJD 57720 (2017 November 28) and MJD 57899 (2017 May 26) are presented. The data quality selection was performed based on the cosmic-ray rate (Hildebrand et al. 2017). The effect of the zenith distance on the cosmic-ray rate was taken into account following Mahlke et al. (2017) and Bretz (2019). A total of 279 h of good-quality data spread over 86 nights was obtained.

The analysis was performed following Dorner et al. (2015) using the standard FACT analysis software described in Bretz & Dorner (2010). Using the excess rate from the Crab Nebula, which is a constant source at TeV energies, the excess rate was corrected for its dependence on the zenith distance and the trigger threshold. The analysis threshold of the light curve computed using Monte Carlo simulation is about 0.95 TeV.

2.3. Fermi-LAT

The LAT instrument is a pair-conversion telescope on board the Fermi satellite (Atwood et al. 2009; Ackermann et al. 2012). It monitors the gamma-ray sky in the 20 MeV to > 300 GeV energy range with an all-sky coverage on a ∼3 h timescale. The analysis for this work was performed using the unbinned-likelihood tools from the FERMITOOLS software³ v1.0.10. We used the instrument response function P8R3_SOURCE_V2 and the diffuse background models⁴gll_iem_v07 and iso_P8R3_SOURCE_V2_v1.

We selected all Source class events between 0.2 GeV and 300 GeV in a region of interest (ROI) with a radius of 15° around Mrk 421. The events with a zenith angle > 100° were discarded so that the contribution from limb gamma rays is reduced. A first unbinned analysis was performed over the entire considered time range, between MJD 57720 and MJD 57918. The source model included all sources in the ROI from the fourth Fermi-LAT source catalogue (4FGL; Abdollahi et al. 2020). During the fit, the normalisation and spectral parameters of sources within a radius of 10° of Mrk 421 were left free to vary, while the remaining sources had their parameters fixed to the 4FGL values. Mrk 421 was modelled with a log-parabolic spectral function, $\frac{\mathrm{d}N}{\mathrm{d}E}=f_0{\left(\frac{E}{E_0}\right)}^{-\alpha -\beta log\left(\frac{E}{E_0}\right)}$ with E₀ = 1286.47 GeV as in the 4FGL catalogue. The normalisations of the background components were left free to vary. After this first fit, sources which resulted in a detection test statistic (TS; Mattox et al. 1996) of less than 15 were removed from the model.

Using the simplified model, light curves with a 3-day binning were built in the 0.2–2 GeV and 2–300 GeV bands. In each time bin, the normalisation and the index α of Mrk 421 were left free to vary. The curvature parameter β was fixed to 0.02, the value obtained from the fit over the entire period and is similar to the 4FGL catalogue. For sources within 10°, only the normalisation was left free, while sources further than 10° had all their parameters fixed to their 4FGL values. In case the target was detected with TS < 5 an upper limit at 95% confidence level was computed after fixing the spectral index at α = 1.78, the average value of the entire campaign.

2.4. NuSTAR

The NuSTAR observations in 2017 were performed in four epochs between MJD 57757 and MJD 57840 with a cadence of roughly one month, and in coordination with the observations performed with Swift and the MAGIC telescopes. See Appendices B and C for the specific dates and times, and its relation with the MAGIC, Swift, and optical observations. The raw NuSTAR data were processed following the description in Baloković et al. (2016), but using the updated software and calibration packages NuSTARDAS (version 1.7.1), HEASoft (version 6.21), and CALDB (version 20170222). We used strict event filtering to eliminate high background fluxes during South Atlantic Anomaly passages (using the flags saamode=strict and tenacle=yes for the nupipeline processing script). The resulting exposures for the single-night observations are in the range 16–25 ks. The source counts were selected from a circular extraction region with a radius of 120″, while the background was sampled from the same focal plane detector excluding a circular region with a radius of 180″ centred on Mrk 421, which is a bright X-ray soure and hence clearly detected by NuSTAR. We binned the spectra ensuring that each energy bin exceeds the signal-to-noise ratio of 3. We also verified that a different selection of data processing parameters does not affect the downstream analysis in any significant way.

In addition to the analysis of data for each of the four epochs, we also separated the data per orbit as well as into 30-min bins for which we have strictly simultaneous coverage at VHE gamma-ray energies with the MAGIC telescopes. We find significant variability in the flux and spectral shape in each of the epochs, as shown in detail in Appendices B and C. The spectral analysis was performed in Xspec (Arnaud 1996). The cross-normalisation between the Focal Plane Module A and B was assumed to be a free parameter in spectral fitting. The distribution of the cross-normalisation factor over the various spectral fits is tight with an average of 0.99 and a standard deviation of 0.02, which is firmly within the expectations of NuSTAR (Madsen et al. 2015). In all cases we employed the log-parabolic model with a pivot energy of 1 keV to describe the spectra, finding significantly non-zero curvature parameter β in longer integrations (e.g., see Table 4). For short observations, owing to the lower photon statistics, in some time intervals β is statistically consistent with zero, implying that the spectra are consistent with a power law in those cases (see Appendix B). We provide fluxes calculated based on this model in the 3–7 keV and 7–30 keV bands, as the source is not significantly detected above the background level at energies higher than 30 keV in short time integrations. All uncertainties are given at the 68% significance level.

2.5. Swift-BAT

The BAT light curves were generated using the BAT survey data, which are collected continuously by the spacecraft and are binned into ∼300 s (Markwardt et al. 2007). The survey data were processed using the standard BAT pipeline, batsurvey⁵, which is part of the HEASoft tools and produces count rate data at the source location for each snapshot image in eight energy bands: 14–20, 20–24, 24–35, 35–50, 50–75, 75–100, 100–150, and 150–195 keV.

The count rate for each image was combined into an eight-band light curve file, and the light curve was further binned into the desired time ranges (1 day, 3 days, and 6 days) and energy band (15–50 keV). When performing the binning, errors were calculated with standard error propagation method using the “BKGVAR” column in the light-curve data, which measures the background variation at the source location.

A BAT spectrum was created using the eight-band information produced by batsurvey around MJD 57788.7, which corresponds to the peak activity in the Swift-BAT (see Sects. 3 and 7). For this, we only used data when the source was on-axis relative to the BAT detector plane, in order to avoid complications due to different instrumental sensitivity at different incident angles. Once a spectrum was created, the corresponding instrumental response file was generated using the standard BAT tool, batdrmgen⁶.

2.6. Swift-XRT

The Swift X-ray Telescope (XRT; Burrows et al. 2005) was used to characterise the emission from Mrk 421 in the energy range from 0.3 keV to 10 keV. The Swift-XRT observations were performed in the windowed timing (WT) readout mode, and the data were processed using the XRTDAS software package (v.3.5.0) developed by the ASI Space Science Data Center (SSDC), and released by the NASA High Energy Astrophysics Archive Research Center (HEASARC) in the HEASoft package (v.6.26.1). The calibration files from Swift-XRT CALDB (version 20190910) were used within the xrtpipeline to calibrate and clean the events.

The X-ray spectrum from each observation was extracted from the summed cleaned event file. Events for the spectral analysis were selected within a circle of 20-pixel (∼46 arcsec) radius, which encloses about 90 per cent of the point-spread function (PSF), centred at the source position. The background was extracted from a nearby circular region of 40-pixel radius. The ancillary response files (ARFs) were generated with the xrtmkarf task applying corrections for PSF losses and CCD defects using the cumulative exposure map.

The 0.3 − 10 keV source spectra were binned using the grppha task to ensure a minimum of 20 counts per bin, and then were modelled in XSPEC using power-law and log-parabola models (with a pivot energy fixed at 1 keV) that include a photoelectric absorption by a fixed column density estimated to be N_H = 1.92 × 10²⁰ cm⁻² (Kalberla et al. 2005). The log-parabola model typically fits the data better than the power-law model, and hence it was used to compute the X-ray fluxes in the energy bands 0.3 − 2 keV, and 2 − 10 keV, which are reported in Appendix D. The fluxes were also computed in the 3 − 7 keV range in order to match the low-energy band of NuSTAR and include them in the VHE versus X-ray correlation study in Sect. 4.

2.7. Swift-UVOT

We selected the observations of the Swift UV and Optical Telescope (UVOT, Roming et al. 2005) between 2016 December and 2017 June acquired in the UV filters W1, M2 and W2. We performed photometry on total exposures of each observations available in the official archive with the same apertures for source counts (the standard with 5″ radius) and background (mostly two circles of 16.5″ radii off the source) estimation. We executed the photometry task in the official software version included in the HEAsoft 6.23 package, by the HEASARC, and then applied the official calibrations (Breeveld et al. 2011) included in the more recent CALDB release (20201026). The source is on “ghost wings” (Li et al. 2006) from the near star 51 UMa in most of the observations, so we checked the wing positions and the astrometry very carefully, excluding stray lights and support structure shadows. This quality control removed four UVOT images, leading to a final sample of 95 good-quality observations. The fluxes were dereddened considering a mean Galactic E(B − V) value of 0.0123 mag (Schlafly & Finkbeiner 2011) and using the Galactic interstellar extinction curve from Fitzpatrick (1999).

2.8. Optical

In this paper, we only used photometry in the Cousins’ R band. All the data were provided by the GLAST-AGILE Support Program (GASP, e.g., Villata et al. 2008, 2009; Carnerero et al. 2017) of the Whole Earth Blazar Telescope⁷ (WEBT, e.g., Villata et al. 2002, 2006; Raiteri et al. 2007, 2017). The observations covered the period from MJD 57716 to MJD 57917, and were provided by 24 telescopes in the following 22 observatories spread over the Northern Hemisphere: Abastumani (Georgia), ARIES (India), AstroCamp (Spain), Belogradchik (Bulgaria), Burke-Gaffney (Canada), Calar Alto⁸ (Spain), Crimean (Crimea), Haleakala (US), Hans Haffner (Germany), KVA observatory (Spain), Lulin (Taiwan), McDonald (US), New Mexico Skies (US), Perkins (US), Rozhen (Bulgaria), Siena (Italy), Sirio (Italy), St. Petersburg (Russia), Teide (Spain), Tijarafe (Spain), Vidojevica (Serbia), West Mountain (US).

The data reduction was performed according to standard prescriptions. To minimise offsets among the different data sets due to the presence of the host galaxy, a common aperture radius of 7.5 arcsec was used. Following Nilsson et al. (2007), we subtracted a host galaxy contribution of 8.2 mJy from the observed flux density and then corrected for a Galactic extinction of 0.033 mag according to the NASA/IPAC Extragalactic Database⁹ (NED).

2.9. Radio

The emission at radio frequencies was characterised with the single-dish telescopes at the Metsähovi Radio Observatory, operating at 37 GHz, at the Owens Valley Radio Observatory (OVRO), that operates at 15 GHz, and the Medicina radio telescope at both 8 GHz and 24 GHz. The data from OVRO were retrieved from the instrument team web page¹⁰, while the data from Metsähovi and Medicina were analysed following the prescription from Teraesranta et al. (1998) and Giroletti & Righini (2020), and provided by the instrument teams specifically for this study. For these three single-dish radio instruments, Mrk 421 is a point source and thus the measurements represent an integration of the full source extension. The size of the radio emitting region is expected to be larger than that of the region of the jet that dominates the X-ray and gamma-ray emission, known to vary on much shorter timescales than the radio emission.

3. Multi-wavelength light curves and spectral behaviours

In Fig. 1, we show the MWL light curves from radio to VHE between MJD 57716 and MJD 57918.

Fig. 1. MWL light curves between MJD 57720 and MJD 57918. From top to bottom: MAGIC (0.2–1 TeV & > 1 TeV) and FACT (> 0.95 TeV), Fermi-LAT (0.2–2 GeV & 2–300 GeV), Swift-BAT (15–50 keV), NuSTAR and Swift-XRT (3–7 keV & 7–30 keV), Swift-XRT (0.3–2 keV & 2–10 keV), UV-optical (Swift-UVW1/UVM2/UVW2 & R-band), OVRO (15 GHz), Metsähovi (37 GHz) and Medicina (8 GHz and 24 GHz). The dotted blue horizontal line in the first two panels from the top shows the 1 C.U. flux in the > 1 TeV & 0.2–1 TeV bands, respectively. The Fermi-LAT data have a 3-day binning. The Swift-BAT light curve is computed in 3-day and 6-day binnings, as well as 1-day binning around the flare. The vertical black dashed lines highlight the dates that show an enhanced VHE activity and appear as outliers in the VHE versus X-ray correlation plots (see Sect. 4).

3.1. Gamma rays

In the first two panels from the top, the MAGIC fluxes (0.2–1 TeV & > 1 TeV) are presented. For both energy bands the corresponding C.U. is depicted with a horizontal blue dotted line. The FACT light curve, whose analysis has an energy threshold of about 0.95 TeV, is also shown in the top panel. From MJD 57755 to MJD 57790, the source is more active than usual. The flux is higher than 1 C.U. for several nights, while Acciari et al. (2014) derived a time-averaged flux over 14 years of ≈0.5 C.U. In particular, a bright flare with a duration of about a day is visible on MJD 57788. On that night, the averaged flux measured by MAGIC is ≈3.5 C.U., while the FACT nightly flux is close to 2 C.U. This difference is due to the longer exposure of the FACT observation combined with intra-night variability. The intra-night variability is illustrated in the FACT 20-min binned light curve shown in Fig. 2. The flux decays along the night with a halving time of about 1 h. A peak activity of ≈7 C.U. is measured at the beginning of the observation. Interestingly, the VHE flare does not seem to be accompanied by a comparable burst in the Swift-XRT pass band (0.3–10 keV) as is typically observed for Mrk 421 (Fossati et al. 2008; Aleksić et al. 2015b; Abeysekara et al. 2020; Acciari et al. 2020). During the VHE flare, the 0.3–2 keV flux remains close to the average of the campaign and the 2–10 keV flux is only twice the average. As will be discussed later, this particular flare appears as an outlier in the VHE versus Swift-XRT correlation. The Swift-BAT daily light curve (red points in Fig. 1) reveals a prominent flux increase close to the VHE flare. However, by using a finer binning (3-h), Fig. 2 shows that the Swift-BAT flux simultaneous to the VHE observations is near the typical state. The 15–50 keV flux exhibits a clear increase by a factor ∼3 only a few hours after the MAGIC and FACT observations. After MJD 57789, both the 15–50 keV and VHE fluxes show a drop. Unfortunately, no VHE data are available during the period with the highest state measured by Swift-BAT.

Fig. 2. Zoom on the MWL light curves around the VHE flare on MJD 57788. Two first panels from the top show the FACT (> 0.95 TeV) and the MAGIC (> 1 TeV & 0.2–1 TeV) light curves. The FACT fluxes are computed on 20-min binning while the MAGIC light curve is nightly averaged. On MJD 57788 and MJD 57789 the MAGIC exposures are ∼40 min and ∼100 min respectively, and are depicted with the width of the horizontal bars of the markers. Third panel from the top: the Swift-BAT light curve with 3-h binning. Two panels from the bottom: the Swift-XRT (0.3–2 keV & 3–7 keV & 2–10 keV) and Swift-UVW1/UVM2/UVW2 fluxes. The horizontal orange lines show the average flux over the entire campaign in each energy regime. Top panel: the orange line is the average from the MAGIC > 1 TeV light curve. In the two lowest panels, the orange line depicts the average in the Swift-XRT 0.3–2 keV and Swift-UVW2 bands.

After MJD 57790, the VHE flux is mostly between ≈0.3 and ≈0.6 C.U. The FACT light curve shows nevertheless that during the last 50 days of the campaign, the TeV flux was persistently higher than 0.6 C.U., including several nights with fluxes above 1 C.U.

The MAGIC observations unveil spectral variability on ∼day timescales. The index α lies between 2.80 and 1.95 depending on the night (see Table A.1). The hardest spectrum coincides with the brightest VHE flare on MJD 57788 and the spectrum is best described by a log-parabola shape with α = 1.95 ± 0.04 and β = 0.19 ± 0.05.

Figure 3 shows the hardness ratio (defined as the ratio between the > 1 TeV and the 0.2 − 1 TeV fluxes) obtained from the MAGIC observations. A clear harder-when-brighter behaviour is observed. When considering the 0.2–1 TeV flux, the Pearson’s coefficient is 0.5 ± 0.1 and the correlation significance, computed following the prescription of Press et al. (2007), is 3.6σ. When using the > 1 TeV flux, this behaviour is even more evident. The Pearson’s coefficient is 0.7 ± 0.1 with a corresponding correlation significance of 5.6σ. MAGIC Collaboration (2021) also found a stronger harder-when-brighter behaviour in the > 1 TeV band during the years 2015 and 2016, when Mrk 421 showed very low X-ray and VHE gamma-ray activity. However, unlike the 2015–2016 data, the data from 2017 do not show any saturation at the highest VHE fluxes.

Fig. 3. Hardness ratio (F>1 TeV/F0.2−1 TeV) versus F>1 TeV (top figure) and versus F0.2−1 TeV (lower figure) from the MAGIC measurements. The Pearson’s coefficients and the corresponding significance (following the prescription of Press et al. 2007) are given in both plots.

The Fermi-LAT light curves show no evidence for strong flaring episodes and have moderate flux variability in comparison to that observed at VHE. In the 0.2 − 2 GeV and the 2 − 300 GeV range, the fluxes are ∼10⁻⁷ cm⁻² s⁻¹ and ∼10⁻⁸ cm⁻² s⁻¹ respectively, which are close to the quiescent state described in Abdo et al. (2011).

3.2. X-ray

In the fifth panel from the top of Fig. 1, we show the light curves obtained during the four NuSTAR observations (MJD 57757, MJD 57785, MJD 57813 and MJD 57840). The observations are simultaneous to MAGIC and Swift. The data are binned orbit-wise, and the light curves are computed in two energy ranges, 3–7 keV and 7–30 keV. The first three observations show on average a relatively similar flux state, while during the last observation the flux is roughly four times lower. A flux variation by a factor ≈2–3 is present in each epoch, except on MJD 57757, where the flux varies only by about 30% around 2 × 10⁻¹⁰ erg cm⁻² s⁻¹. We estimated the flux doubling or halving time t_1/2 in the 7-30 keV band (which shows the largest variability) using the prescription of Zhang et al. (1999). The time t_1/2 varies between 4 h and 11 h. It is on the night of MJD 57785 that the source shows the shortest variability. Such variability timescales are similar to those derived by Baloković et al. (2016). More details can be found in Appendix B, where the orbit-wise fluxes along with the exact values of t_1/2 are presented.

The orbit-wise spectra are well described with a log-parabolic model and the results are also listed in Appendix B. Similarly to the Swift-XRT spectral behaviour, the curvature parameter β is not significantly dependent on the flux. In fact, one notes that the curvature on MJD 57740 is comparable to the rest of the observations despite a flux that is about four times lower.

Given the low variability in β, the spectral hardness versus the flux was studied by performing a second series of fits after fixing β to 0.22, which is the mean value of the orbit-wise spectra. The results are shown in Fig. 4. The low state on MJD 57840 is characterised by the softest spectra with α ≈ 2.3 − 2.6. MJD 57757 and MJD 57813 show on average a very similar behaviour with α ≈ 2.0 − 2.1. It is on MJD 57785 that the hardness is the strongest and α ≲ 2.0 for most of the orbits. A linear fit to α versus the 3–7 keV flux gives a slope of −0.46 ± 0.01. This result is in good agreement with Baloković et al. (2016).

Fig. 4. Log-parabola photon index α versus the 3–7 keV flux of the orbit-wise NuSTAR observations. The parameter α is fitted after fixing spectral curvature β = 0.22 in the log-parabolic model. Each day is plotted with a different colour. The black line represents a linear fit, while the grey area is its uncertainty. The resulting slope is −0.46 ± 0.01.

The X-ray observations from Swift-XRT (sixth panel from the top of Fig. 1) display a large variability amplitude. The fluxes lie between F_{0.3−2 keV} ≈ 2 × 10⁻¹⁰ erg cm⁻² s⁻¹ and F_{0.3 − 2 keV} ≈ 10⁻⁹ erg cm⁻² s⁻¹ in the 0.3–2 keV energy range. In the 2–10 keV band, they vary from F_{2 − 10 keV} ≈ 3 × 10⁻¹¹ erg cm⁻² s⁻¹ to F_{2 − 10 keV} ≈ 10⁻⁹ erg cm⁻² s⁻¹, which represents a change by more than a factor 30. The highest X-ray state is registered on MJD 57860 for all energy bands (0.3–2 keV, 3–7 keV and 2–10 keV). A second peak flux is also visible in each band around MJD 57760, after a quasi monotonic increase over more than 40 days (from ≈MJD 57720 to ≈MJD 57760). An indication of an enhanced flux in the Swift-BAT 3-day binned light curve is visible during these high-activity periods and the flux is about three times the campaign average. Unfortunately, no simultaneous VHE observations are available. We note that the flux level observed during those time periods remains moderately high compared to previous published works on Mrk 421 flares (see Hervet et al. 2019 for a 12-year study of the Mrk 421 X-ray flux).

Table D.1 lists for each Swift-XRT observation the spectral parameters of the best-fit log-parabolic and power-law models together with the fluxes in the 0.3–2 keV, 3–7 keV and 2–10 keV bands. For the majority of them, the log-parabolic model provides a better description. The curvature parameter β manifests a very weak dependence on the integral flux, differently from the index α. Therefore, in order to study the spectral hardness as function of the flux, each observation was fitted again with a log-parabolic model with β fixed to 0.16, which is the mean value of the dataset. This procedure removes the correlation between the two parameters of the model, which allows a more illustrative description of the evolution of the hardness versus flux to be obtained. The resulting α versus F_{3 − 7 keV} is shown in Fig. 5. A large variation of α is visible, ranging from ≈1.8 to ≈2.8, with a clear harder-when-brighter behaviour, as is typical in Mrk 421 (Aleksić et al. 2015a; MAGIC Collaboration 2021) and HBLs in general (Pian et al. 1998; Krawczynski et al. 2004; MAGIC Collaboration 2020a). A linear fit results in a slope of −0.64 ± 0.01, and corroborates the harder-when-brighter behaviour visible in the NuSTAR data. For the linear fit we ignored the spectral fits resulting in a p-value¹¹ below 5 × 10⁻² in order to remove the few nights (only 18 out of 107) that are not well described by the spectral model. These nights are depicted with hollow markers in Fig. 5. One can see that they follow closely the general trend and do not show any particular trend or clustering. Thus, ignoring these measurements does not introduce any significant bias to our results.

Fig. 5. Log-parabola photon index α versus the 3–7 keV flux from the Swift-XRT observations. The log-parabolic fits are performed after fixing the spectral curvature β = 0.16. The red cross represents the flare on MJD 57788 seen at VHE energies. The black line represent a linear fit, while the grey area (hardly visible in the plot) is its uncertainty. The slope of α versus the 3–7 keV flux is −0.64 ± 0.01. Hollow markers depict fits with a p-value below 5 × 10−2, which are not considered in the linear fit.

The observation corresponding to the VHE flare (MJD 57788) is characterised by a hard X-ray spectrum despite a moderate flux value (F_{0.3 − 2 keV} and F_{2 − 10 keV}). The latter is shown as a red maker in Fig. 5 and the corresponding α is around 1.8, which is among the hardest spectra measured during the campaign.

In the 15–50 keV band, the Swift-BAT observations do not reveal any particular flaring activity besides the one around MJD 57788 that is already discussed in Sect. 3.1. Although being strongly variable (see Sect. 3.5), the average flux is close to the value derived between 2008 and 2010 by Abdo et al. (2011) when Mrk 421 showed quiescent activity. Using 3-day and 6-day binning, the average flux is at the level of ∼10⁻³ cm⁻² s⁻¹.

3.3. UV-optical

The UV fluxes (from Swift-UVOT) and the R-band fluxes (from GASP-WEBT) are shown in the seventh panel of Fig. 1. They follow a very similar temporal evolution, which is expected given their proximity in energy. Interestingly, they show a continuous flux decay from ≈MJD 57720 to ≈MJD 57760, contrary to the almost monotonic increase in the X-rays during the same period. In addition, the low X-ray activity visible around MJD 57840 is accompanied with high fluxes in both the UV and R-band. This suggests an anti-correlation between X-ray and UV-optical over the campaign. This latter characteristic is investigated in Sect. 5.

3.4. Radio

The light curves from OVRO, Medicina and Metsähovi in the bottom panel unveil no strong variability nor flaring episode. The average flux levels are ∼0.5 Jy at 8 GHz and 15 GHz (Medicina and OVRO) and ∼0.4 Jy at 24 GHz and 37 GHz (Medicina and Metsähovi), respectively. Based on the ∼5-year OVRO and Metsähovi light curves presented in Hovatta et al. (2015), this denotes a typical state during non-flaring activity.

3.5. Multi-wavelength variability

We studied the broadband variability based on the fractional variability, F_var, defined in Vaughan et al. (2003). The uncertainty was computed following the strategy from Poutanen et al. (2008) and the implementation described in Aleksić et al. (2015c). The fractional variability quantifies the variance of the flux normalised to the mean value after subtracting the additional variance caused by the measurement uncertainties. The F_var value is naturally affected by the instrument sensitivities, the binning and the flux sampling. Therefore, great care must be taken when comparing results from different telescopes. We refer the reader to Aleksić et al. (2015a) and Schleicher et al. (2019) for a detailed study of the caveats inherent to the fractional variability method.

The F_var values are shown in Fig. 6 for each of the energy bands. They were computed using a nightly binning for the MAGIC, FACT, Swift-XRT, Swift-UVOT, R-band and radio light curves. Because of the limited sensitivity to detect Mrk 421 on timescales of one day, for Fermi-LAT and Swift-BAT, we adopted a 3-day binning over which to integrate the data and compute the fluxes. Solid markers include all data from Fig. 1. A discrepancy between the MAGIC (> 1 TeV) and FACT F_var is visible and is explained by the different nightly averaged flux measured during the flare on MJD 57788 due to the different integration time, as previously mentioned. When the day of the flare is ignored, the F_var values are fully compatible. All of the MAGIC/FACT/Swift-XRT/BAT/UVOT measurements that are separated from one another by less than 4 h were considered and are shown with hollow markers. Here again, the Swift-BAT fluxes were computed with a 3-day binning.

Fig. 6. Fractional variability Fvar obtained from the light curves shown in Fig. 1. MAGIC, FACT, Swift-XRT, Swift-UVOT, R-band and radio fluxes are nightly binned. Fermi-LAT and Swift-BAT fluxes have a 3-day binning. Results from each instrument are plotted in different colours. The filled markers include all data. The hollow markers include VHE and Swift data lying within a time window of 4 h from each other.

Overall, a clear two-peak structure is visible. The first peak, with F_var ≈ 1, occurs in the hard X-ray band (15–50 keV). The second peak lies in the VHE band (around 1 TeV), also with F_var ≈ 1. The lowest variability is seen in the radio, UV-optical and 0.2–2 GeV band and they all display F_var <  0.2.

In Fig. 7, the fractional variability is shown for the simultaneous MAGIC/NuSTAR observations. Here, again, a significant increase of variability with energy is observed both in the X-rays and in the VHE band.

Fig. 7. Fractional variability Fvar from the simultaneous MAGIC-NuSTAR observations performed on the four nights MJD 57757, MJD 57785, MJD 57813 and MJD 57840 (2017 January 4, 2017 February 1, 2017 March 1 and 2017 March 28). The Fvar values were computed with the fluxes determined on 30-min time bins that are reported in Appendix C.

The patterns reported above are common for Mrk 421, and have been observed on multiple occasions (Aleksić et al. 2015a,b; Baloković et al. 2016). The locations of the two peaks directly correspond to the falling edges of the two bumps noticeable in the SED of Mrk 421 (Abdo et al. 2011). The low F_var values (radio, UV-optical and MeV-GeV) on the other hand match the rising edges of the two SED bumps. In leptonic models, the rising segments of the SED bumps originate from less energetic electrons compared to the falling segments. Because the cooling rate of the particles due to synchrotron radiation is inversely proportional to their Lorentz factor (t_{cool, synch} ∝ 1/γ), the variability is naturally expected to increase with energy around the two SED bumps.

We investigated the flux correlations among the various energy bands that show significant variability, as reported in Fig. 6, finding significant correlations only between the VHE gamma rays and the X-rays and between the X-rays and the UV-optical. The results from these studies are reported in Sect. 4 and 5.2, respectively.

4. Study of the VHE versus X-ray correlation

4.1. Correlation during the MAGIC/NuSTAR/Swift observations

First, we investigated the VHE versus X-ray correlation making use of the strictly simultaneous MAGIC, NuSTAR, and Swift observations performed during four days (MJD 57757, MJD 57785, MJD 57813, and MJD 57840). The combination of NuSTAR and Swift-XRT observations provides an X-ray coverage over two orders of magnitude, from 0.3 keV to 30 keV. Up to now, only a few published works have reported multi-band correlation studies extending into the hard X-rays (i.e., ≳10 keV). We note that in the previous works, the source was either probed during exceptionally low states (Baloković et al. 2016) or during flares (Fossati et al. 2008; Acciari et al. 2020). In this paper, we aim to complete the picture by offering a view of the VHE versus X-ray correlation of Mrk 421 during a close-to-typical activity.

We built light curves for the MAGIC/NuSTAR/Swift strictly simultaneous observations with temporal bins of 30 min for all the instruments. This choice of binning allows sufficiently good statistics in all energy bands, in particular at VHE, where the measurement instruments typically have lower sensitivities compared to the X-ray telescopes. In the X-ray regime, the fluxes were calculated in three energy bands: 0.3–3 keV, 3–7 keV, and 7–30 keV. In the VHE regime, the fluxes were computed in the 0.2–1 TeV and > 1 TeV bands. The resulting light curves are shown in Appendix C. No significant intra-night variability is found in the MAGIC light curves. On the contrary, the fluxes measured by NuSTAR are significantly (> 5σ) incompatible with a constant value for each epoch. For completeness, the light curves in Appendix C also include the densely sampled optical (R-band) data from GASP-WEBT, from which no clear variability pattern is visible.

In Fig. 8, the flux-flux correlations are shown for all energy combinations. For each day, the data are plotted with a different colour. As mentioned above, the X-ray fluxes are quite similar between the first three pointings (MJD 57757, MJD 57785 and MJD 57813). The VHE fluxes, however, differ significantly between those nights, being ≈1.6 C.U. for MJD 57757 and MJD 57785 and ≈0.5 C.U. for MJD 57813. Such a feature points towards significant variability in the Compton dominance of the SED (that is defined as the relative luminosity of the high-energy SED component with respect to the low-energy SED component). The last observation, on MJD 57840, has the lowest VHE flux (≈0.1 − 0.3 C.U.), which coincides with a lower X-ray activity level.

Fig. 8. VHE flux versus X-ray flux during the simultaneous MAGIC/NuSTAR/Swift observations. Flux points are computed over time bins of 30 min. From left to right: the X-ray energy bands in the panels are 0.3–3 keV, 3–7 keV and 7–30 keV. In the upper panels, the MAGIC fluxes are in the > 1 TeV band, while in the lower panels they are in the 0.2–1 TeV band. Data points from Swift-XRT are shown with black-filled markers. Fluxes corresponding to each day are plotted in a different colour. The results of the Pearson coefficients (and its significance in σ), DCF and the slope of the linear fits in the log-log plane (with the χ2/d.o.f. values) are indicated in each of the corresponding subplots.

In order to quantify the correlation, we used the Pearson coefficient and the discrete correlation function (DCF; Edelson & Krolik 1988). Compared to the Pearson coefficient, the DCF has the advantage of naturally including the measurement uncertainties. A linear fit in the log-log plane was used to further study the type of the correlation patterns. All the resulting numbers and the slopes of the linear fits (as well as the χ²/d.o.f. results) are included in the corresponding subplots of Fig. 8.

Each Pearson coefficient shows a value above 0.8 with a significance higher than 5σ. The strong correlation is further confirmed by the DCF, whose values are all inconsistent with 0 at a significance greater than 4σ. At the lowest X-ray energies in the 0.3–3 keV band, the Swift-XRT coverage is much less compared to that of NuSTAR. This prevents a direct comparison of the Pearson coefficients and the DCF obtained at 0.3–3 keV with those obtained in the 3–7 keV and 7–30 keV bands. When restricting ourselves to the 3–7 keV and 7–30 keV panels, the DCF as well as the Pearson coefficients are compatible within statistical uncertainties between the different energy combinations. We stress that the correlations are strongly dominated by the day-by-day flux variations and no significant intra-night VHE versus X-ray correlation is found.

Regarding the correlation slopes, they vary significantly depending on the energies considered. The largest slope is derived for the > 1 TeV versus 0.3–3 keV case and is consistent with a more-than-cubic trend. Conversely, the lowest slope is obtained for 0.2–1 TeV versus 7–30 keV, and indicates a sub-linear trend between these two energies.

It is interesting to note that a significant increase of the slope is visible at all X-ray energies when considering a higher VHE band (i.e., when moving from 0.2–1 TeV to > 1 TeV in a given X-ray band). These trends indicate that the > 1 TeV flux has a stronger scaling with the X-ray flux compared to the 0.2–1 TeV flux.

In a given VHE band, the slopes are compatible within uncertainties between the 3–7 keV and 7–30 keV bands, although the best-fit values seem to indicate lower slopes for the 7–30 keV band. A significant decrease of the slope at higher X-ray energies was reported by Baloković et al. (2016), which we can not confirm in our work with the data at hand. We refrain from a comparison of the 0.3–3 keV band with the other X-ray bands because of the significantly smaller amount of available measurements.

4.2. Correlation over the entire MWL campaign

In this section, we report the MAGIC versus Swift-XRT correlation study over the full MWL campaign. In order to ease the comparison with the results shown in Fig. 8, Swift-XRT fluxes were considered in the 0.3–3 keV and 3–7 keV bands. We correlated MAGIC and Swift-XRT daily averaged fluxes by considering all pairs of observation falling within a time range of δt = 4 h (28 pairs in total fulfil this criterion). The flux-flux plots are shown in Fig. 9.

Fig. 9. VHE flux versus X-ray flux over the MWL campaign using MAGIC and Swift-XRT data. Data are nightly binned, and only pairs of measurement within 4 h are considered. MAGIC fluxes are in the > 1 TeV band (top panels) and in the 0.2–1 TeV band (lower panels). Swift-XRT fluxes are computed in the 0.3–3 keV (left panels) and 3–7 keV bands (right panels). The open red diamond markers highlight measurements on MJD 57757 and MJD 57785 (i.e., the dates of the first two simultaneous MAGIC/NuSTAR/Swift observations). The VHE flare (on MJD 57788) is plotted in open red circle marker. Black lines depict linear functions (in the log-log plane) with slopes x = 1,  2,  3. In the 0.2–1 TeV versus 3–7 keV panel, because of the dynamical range of the data, the linear function with slope x = 3 is replaced with one with slope of x = 0.5. The Pearson coefficients, DCFs and the slopes from linear fits in the log-log plane are shown in Table 1.

The flare on MJD 57788 appears as a clear outlier (plotted with an open red circular marker). The F_>1 TeV and F_{0.2−−1 TeV} are respectively ∼10 times and ∼3 times higher compared to the bulk of the measurements with comparable (and even slightly higher) F_{0.3−3 keV} values. When considering the 3–7 keV band, the VHE flare indeed coincides with the highest X-ray flux among those that lie within δt from a MAGIC observation. This corroborates the fact that in the X-rays the flare is characterised by a hardening of the spectrum without significant increase of the flux amplitude around 1 keV. We stress however that the corresponding F_3−7 keV is only ≈20% greater than the second highest F_3−7 keV (which corresponds to MJD 57763), for which the F_>1 TeV is about seven times lower. We note that for this particular night, the MAGIC and Swift-XRT observations were simultaneous and overlapped, as shown in Fig. 2. This allows us to discard any bias due to the non-simultaneity of the observations. This is particularly important given the rather short duty cycle of the flare at VHE (flux halving time of a few hours, see Fig. 2).

Interestingly, the first two NuSTAR nights on MJD 57757 and MJD 57785, which show an enhanced VHE activity in Fig. 8, also seem to stand out from the distribution. They are plotted with diamond open markers in Fig. 9. The VHE flux is higher than the bulk of the data. Together with the MJD 57788 flare, these measurements suggest a change in the VHE versus X-ray correlation with respect to the typical behaviour.

The correlation was quantified following the same methodology described in the previous section. The Pearson coefficients, the DCF and the slopes from the linear fits are listed in Table 1. The entire dataset (28 measurements) is considered. The Pearson coefficients are consistent within uncertainties between each energy combination and are all greater than 0.7. Also, the significances are all above 4σ. The DCF values are similar to the Pearson coefficients and do not significantly vary between the energy combinations. On the other hand, the slopes show distinct values depending on the energy band considered, similarly to the simultaneous MAGIC/NuSTAR/Swift observations. A roughly cubic trend is obtained for > 1 TeV versus 0.3–3 keV, while an almost linear trend is visible for 0.2–1 TeV versus 3–7 keV. Finally, a close-to-quadratic relation is found for 0.2–1 TeV versus 0.3–3 keV and for > 1 TeV versus 3–7 keV. Within uncertainties, the slopes are consistent with those obtained from the MAGIC/NuSTAR/Swift nights (in the 0.3–3 keV and 3–7 keV bands).

Similarly to the correlation during the MAGIC/NuSTAR/Swift nights, Fig. 9 reveals that the correlation slope increases at higher VHE energies in a given X-ray band. This pattern is also verified in Table 1, where the best-fit slopes in the > 1 TeV band are systematically greater than in the 0.2–1 TeV band.

The χ²/d.o.f. values indicate a poor fit and implies that a $F_{\mathrm{VHE}}\propto F_{\mathrm{X}-\mathrm{ray}}^x$ parametrisation is too simple to describe the relation between these two bands. Nonetheless, the best-fit slopes still provide information about the main trend in the relation between the fluxes from these two energy bands. The above mentioned results suggest steep correlations maintained over monthly timescales.

Katarzyński et al. (2005) and Katarzyński & Walczewska (2010) argued that in one-zone SSC models one would expect in general linear or sub-linear trends between the X-ray and VHE emission in HBLs. One should note that more-than-linear correlations may still occur in some conditions, but demand rather specific evolution of the model parameters (see Katarzyński et al. 2005). This is particularly true when one considers the emission above the peaks of the SED components, which is typically the case for the > 1 TeV and 3–7 keV bands in Mrk 421 (see Sect. 5). From this point of view, under the assumption of a one-zone SSC scenario, it is surprising to derive over the full campaign a slope of 1.80 ± 0.19 for > 1 TeV versus 3–7 keV fluxes. However, from a simple quality check, we found that the previously mentioned three outliers (marked in open red markers in Fig. 9), which show enhanced gamma-ray activity, significantly influence the best-fit slope towards higher values even though they only represent about 10% of the measurements. When ignoring them, the correlation agrees with a linear trend. The mentioned outliers seem to follow a different pattern in the log-log plane, implying that a physical regime different from a more typical behaviour was probed. Therefore, these results suggest a temporal dependence in the correlated variability between the X-ray and VHE bands.

In Fig. 10, we show the results of the correlation between FACT and Swift. As was done previously, only pairs of observations that occurred within a time window of δt = 4 h of one another were considered. The fluxes are nightly binned and only FACT measurements with a signal-to-noise ratio > 2 are taken into account. The results are listed in Table 2. The three outliers MJD 57788 (the VHE flare), MJD 57757 and MJD 57785 are depicted with open blue symbols in Fig. 10. The derived correlation significance is ≳4σ in all energy bands. The DCF and Pearson’s coefficients are in good agreement with those obtained using the MAGIC data.

Fig. 10. VHE flux versus X-ray flux over the MWL campaign using FACT and Swift-XRT data. Data are nightly binned, and only pairs of measurement within 4 h are considered. FACT fluxes are computed above 0.95 TeV and only measurements with a signal-to-noise ratio > 2 are taken into account. Swift-XRT fluxes are computed in the 0.3–3 keV (left panel) and 3–7 keV bands (right panel). The open blue markers highlight measurements on MJD 57757 and MJD 57785 (i.e., the dates of the first two simultaneous MAGIC/NuSTAR/Swift) and the VHE flare (on MJD 57788). Black lines depict linear functions (in the log-log plane) with slopes x = 1 and x = 2. The Pearson coefficients, DCF and slope are shown in Table 1. The Pearson coefficients, DCF and slopes from linear fits in the log-log plane are shown in Table 2.

5. Investigation of the behaviour of the low-energy component in the SED

5.1. Synchrotron peak frequency

By definition, the low-energy component of the SED of HBLs peaks at frequencies above 10¹⁵ Hz. In the case of Mrk 421, the peak frequency generally lies between ∼10¹⁷ Hz and ∼10¹⁸ Hz (Tramacere et al. 2007a), making it an extreme synchrotron blazar (Biteau et al. 2020). However, strong variability has been reported in recent works. For instance, Baloković et al. (2016) derived a peak located at frequencies ≲10¹⁶ Hz during particularly low states, while during flares, values close to 10¹⁹ Hz were reported by Tramacere et al. (2009). Motivated by the large X-ray flux variability amplitude shown in Fig. 1, which is more than a factor 30 in the 2–10 keV band, we carried out an in-depth study of the low-energy component.

From optical to X-rays, the SED is usually well described with a log-parabola shape (Tramacere et al. 2007b):

$$\begin{array}{c}\hfill \nu F_{\nu }(\nu )=f_0·{10}^{-b·{({log}_{10}(\nu /{\nu }_{\mathrm{s}}))}^2}{\mathrm{erg}\mathrm{cm}}^{-2}{\mathrm{s}}^{-1},\end{array}$$(2)

where ν_s is the peak frequency, b the curvature and f₀ the normalisation constant. Equation (2) is equivalent to Eq. (1), with the former now having the peak frequency as a free parameter. In order to locate the peak frequency throughout the campaign, we fitted Eq. (2) to all available pairs of simultaneous Swift-XRT and Swift-UVOT observations. Most of the time, the Swift-XRT and Swift-UVOT instruments were operating together, resulting in a total of 93 pairs of observations. Furthermore, we complemented the SED with Swift-BAT data. For this, we exploited the orbit-wise light curve in the 15–50 keV range and computeed a 6-h averaged photon count rate (in cm⁻² s⁻¹) centred at each of the simultaneous Swift-XRT and Swift-UVOT observations. The Swift-BAT requires a longer integration time because of its limited sensitivity to detect Mrk 421 at hard X-rays in short timescales. The 6-h averaged count rates were converted to a flux point (in erg cm⁻² s⁻¹) using the prescription of Krimm et al. (2013). In Fig. 11, the resulting locations of the peak frequencies, ν_s, are plotted versus time. We show in blue square markers the results from fits that have a p-value above 10⁻³, in violet triangles the results from fits with p-values between 10⁻³ and 10⁻⁴ and in red diamond markers the results from fits with p-values below 10⁻⁴. Around 70% of the fits have a p-value above 10⁻³, and only 17 fits (less than 20%) yielded a p-value below 10⁻⁴, indicating that Eq. (2) provides an acceptable parametrisation of the SED in most of the cases. In Fig. 12, the resulting locations of the peak frequencies, ν_s, are reported in the νF_ν(ν)−ν space. Archival synchrotron peak frequencies are depicted with dashed areas in Fig. 12.

Fig. 11. Evolution of the synchroton peak frequency νs versus time throughout the MWL campaign. The peak frequencies are obtained by fitting a log-parabola to simultaneous Swift-XRT, Swift-UVOT and Swift-BAT observations. Blue square markers indicate fits that result in a p-value above 10−3. Data points in violet triangles correspond to a p-value between 10−3 and 10−4, and red diamond markers have a p-value lower than 10−4.

Fig. 12. Synchrotron peak frequencies throughout the MWL campaign obtained from log-parabola fits to simultaneous Swift-XRT, Swift-UVOT and Swift-BAT observations. Blue square markers indicate fits that result in a p-value above 10−3. Data points in violet triangles correspond to a p-value between 10−3 and 10−4, and red diamond markers have a p-value lower than 10−4. The black dotted line is a linear fit to the data and the obtained slope is indicated on the plot. The dashed areas represent the ranges of archival synchrotron peak frequencies, and are taken from Fig. 12 of Baloković et al. (2016). The respective references are given in the legend.

Figure 12 reveals a strong variability of ν_s. The peak can be as low as ∼10¹⁶ Hz, similarly to the particularly low activity reported by Baloković et al. (2016). The highest ν_s are ∼10¹⁸–10¹⁹ Hz. Such variations roughly cover the entire range of ν_s, which has been observed up to now in the case of Mrk 421. The interesting aspect of our result stems from the fact that this large variability amplitude only happened within a ≈6-month period. The strong flux variability amplitude in the 0.3–2 keV and 2–10 keV bands over the MWL campaign (see Sect. 3) is therefore accompanied by particularly large shifts of the synchrotron bump.

A linear fit was performed in the log-log space of ν_sF_ν(ν_s) versus ν_s. The slope obtained is 0.24 ± 0.02, indicating a significant increase of the X-ray flux with higher peak frequencies. A very similar slope (0.25 ± 0.02) is derived when only the log-parabola fits with a p-value above 10⁻³ (blue data points) are considered. This result is consistent with the results reported in Sect. 3.2, and agrees with the well-known harder-when-brighter behaviour previously reported for Mrk 421 and other HBLs (Pian et al. 1998; Krawczynski et al. 2004; Aleksić et al. 2015a). However, the scatter of the data over the main trend shows that there are many days when the synchrotron emission of Mrk 421 did behave substantially differently. This may be caused by different particle or environment conditions in the main emitting region, or because of the presence of additional components that, during short periods of time, contribute significantly to the overall synchrotron emission of Mrk 421. Either way, the ν_sF_ν(ν_s) versus ν_s data for these days appear as outliers in the main trend with slope of about 0.2, shown in Fig. 12.

Figure 13 shows ν_sF_ν(ν_s) versus the curvature parameter b as defined in Eq. (2). The data tend to indicate a lower curvature with increasing ν_sF_ν(ν_s) values, although this anti-correlation remains weak and some outliers exist. A linear fit yields a slope of −1.17 ± 0.14. A very similar slope is obtained when one considers only the log-parabola fits having a p-value above 10⁻³ (blue data points). A decrease of the curvature with the flux may be indicative of an evolution of the flux driven by stochastic electron acceleration, as demonstrated by Tramacere et al. (2011).

Fig. 13. νsFν(νs) versus the curvature parameter b throughout the MWL campaign obtained from log-parabola fits to simultaneous Swift-XRT, Swift-UVOT and Swift-BAT observations. As in Fig. 12, blue square markers indicate fits that result in a p-value above 10−3. Data points in violet triangles correspond to a p-value between 10−3 and 10−4, and red diamond markers have a p-value lower than 10−4. The black dotted line is a linear fit to the data, and the obtained slope is indicated on the plot.

5.2. UV-optical versus X-ray anti-correlation

Within standard one-zone SSC models, which assume a single emission zone responsible for the SED (from optical to VHE), a positive correlation between the UV-optical and X-ray is generally expected. Changes in parameters describing the emission zone environment (for instance Doppler factor, electron density or the size of the region) would modify the synchrotron emission at each energy, and thus, lead to correlated variability between the UV-optical and X-ray emissions. Nevertheless, the very different temporal behaviour between these two energy bands renders the detection of a correlation difficult. Figure 6 shows that F_var is more than two times higher in the X-rays than in the UV-optical. Hence, for any given flux change in the X-rays, the correlated UV-optical flux variation would be much suppressed. On the other hand, the large variability of the synchrotron peak frequency discussed may facilitate the detection of correlation patterns.

We cross-correlated UV-optical with X-ray fluxes by selecting all possible combinations between the following light curves: optical (R-band), Swift-UVOT, Swift-XRT(0.3–2 keV) and Swift-XRT(2–10 keV). Regarding Swift-UVOT, we considered data from the W2 filter. Given their proximity in frequency, fluxes in the three Swift-UVOT filters are strongly correlated and all of them yield very similar results. The choice of the W2 filter is further motivated by the slightly larger amount of measurements compared to W1 and M2, providing a better coverage of the source. As in Sect. 4, we used the DCF (Edelson & Krolik 1988) to quantify the correlation. It is known that the significance of the correlation assessed directly from the DCF uncertainty given in Edelson & Krolik (1988) can be overestimated (Uttley et al. 2003). This is especially the case for red-noise light curves that are regularly found in blazar observations. In this work, a more careful statistical treatment employing Monte-Carlo light curve simulations was undertaken. The details of the procedure can be found in Appendix E.

The DCF coefficients between Swift-XRT (0.3–2 keV) and Swift-XRT (2–10 keV) on one hand, and Swift-UVOT/W2 and R-band on the other hand, are shown in Fig. 14. The green full lines depict the DCF results obtained with the data, while the blue and red dotted lines are the 2σ and 3σ significance bands (derived by simulations), respectively. Figure 14 reveals a negative peak with DCF ≈ −0.75 and DCF ≈ −1 located between −3 and −5 days. This clear feature is visible in the cross-correlations derived for all the energy bands. The significance is at the level of 3σ for Swift-XRT (0.3–2 keV) versus Swift-UVOT/W2, Swift-XRT (2–10 keV) versus Swift-UVOT/W2 and Swift-XRT (0.3–2 keV) versus R-band. For Swift-XRT (2–10 keV) versus R-band the significance is above 3σ. These results represent evidence for an anti-correlation between UV-optical and the X-rays.

Fig. 14. DCF between X-ray (0.3–2 keV and 2–10 keV) and UV-optical (Swift-UVW2 and R-band) during the MWL campaign. The blue- and red-dashed-lines indicate the 2σ and 3σ confidence intervals estimated from the Monte Carlo simulations, as described in the text.

The multi-instrument LCs from Fig. 1 show that, during the time interval from about MJD 57720 and about MJD 57760, there is an overall decrease in the UV-optical fluxes with time, while the opposite occurs with the X-ray fluxes, that show an overall increase with time. Such a clear trend does not seem to be apparent for time intervals after MJD 57760, and hence it may happen that the above-mentioned time range of about 40 days dominates the marginally significant anti-correlation reported in Fig. 14. To investigate this, we repeated the study and computed the DCF as well as its significance after excluding all data taken prior to MJD 57760. The results are displayed in Fig. F.1 and show that the anti-correlation disappears (significance drops below 3σ) for all the bands probed.

Additionally, we also evaluated the correlation using only data within the time interval MJD 57720–57760. The results are shown in Fig. F.2. This time, the significance is above 3σ for the cross-correlations from all energy bands (the DCF values are also higher), confirming the overall dominance of the time interval MJD 57720–57760 in the cross-correlation results. This observation suggests that these two bands are not persistently related to each other (as occurs between the VHE gamma rays and the X-rays in Mrk 421), but instead show a degree of correlation that changes with time, and is significant only during month timescales.

6. Broadband SED and modelling of the strictly simultaneous MAGIC/NuSTAR/Swift observations

We built four broadband SED snapshots from radio to VHE around the strictly simultaneous MAGIC/NuSTAR/Swift observations. The excellent coverage in the X-rays, combined with the HE-VHE data, offer unprecedented time-resolved SEDs during rather typical states for Mrk 421.

At VHE, the MAGIC spectra were averaged over each observing night. The effects of EBL absorption were corrected for by applying the EBL model of Domínguez et al. (2011). The resulting best-fit parameters of a log-parabola spectral shape are shown in Table 3. A log-parabola model is preferred with respect to a simple power-law model at a significance higher than 5σ (based on a likelihood ratio test) for the four spectra. The NuSTAR and Swift spectra were computed over the exact same time range as the MAGIC observations, limiting possible biases due to non-simultaneity. We note though that the NuSTAR observations consisted of continuous exposures, while for Swift the observations were shorter and scattered over the MAGIC window (see Figs. C.1–C.4). The NuSTAR best-fit spectral parameters assuming log-parabolic models are shown in Table 4. In the high energy (HE) band, owing to the limited ability to significantly detect Mrk 421 in this energy range, the Fermi-LAT spectra were averaged over 3 days, which provides a good compromise between simultaneity and spectral resolution. Optical (R-band) data are not strictly simultaneous, but the closest observation in time was chosen, which resulted in a time difference of at most one day. Because of the low variability in the optical band reported in Figs. 1 and 6, they can be considered as a good proxy for the simultaneous emission in the optical band.

The broadband SEDs are shown with brown data points in Fig. 15. The SED data describing a quiescent state (averaged over 5 months) reported in Abdo et al. (2011) are also plotted, for comparison purposes, with grey full dots. A slight mismatch is visible between Swift-XRT and NuSTAR in their overlapping region (around 10¹⁸ Hz). For the X-ray spectra from MJD 57757, MJD 57785 and MJD 57840, the mismatch is always less than 15%, which is within the systematic uncertainties of these two instruments (Madsen et al. 2017). On MJD 57813 (lower left panel), the mismatch is somewhat larger. The last Swift-XRT point, the one with the largest discrepancy, is 34%±17% off from the lowest-energy point of the NuSTAR spectrum. We note that this remains at the level of 2σ, and thus may well be due to a statistical fluctuation. Moreover, the X-ray flux variability on hour timescales and the longer exposures of NuSTAR (see light curves in Appendix B) may also contribute to this small mismatch between the X-ray fluxes measured by these two instruments.

Fig. 15. SED of the four simultaneous MAGIC/NuSTAR/Swift observations (MJD 57757, MJD 57785, MJD 57813, MJD 57840). The data points and obtained one-zone SSC models are plotted with brown markers and orange lines, respectively. The parameters of the one-zone SSC models are shown in Table 5. For comparison purposes, the one-zone SSC model of the first pointing (MJD 57757) is plotted with black dotted lines in all the other subplots. The Compton dominance parameter AC determined from the model is indicated on each panel. Archival data representing the typical Mrk 421 state from Abdo et al. (2011) are shown in grey.

The SEDs show distinct spectral and flux characteristics across the whole spectrum. While the first three observations (MJD 57757, MJD 57785 and MJD 57813) display a comparable X-ray emission with a synchrotron peak frequency located around 10¹⁷ Hz, the last observation shows a much softer and fainter spectrum as well as a synchrotron peak shifted to lower frequencies by an order of magnitude, at ∼10¹⁶ Hz.

As mentioned in Sect. 4, the first two observations on MJD 57757 and MJD 57785 have an enhanced VHE activity compared to that of MJD 57813, despite the three nights having a comparable synchrotron emission up to the hard X-rays. Figure 15 further shows that the enhanced VHE activity is also reflected in the Fermi-LAT data, revealing an enhancement of the whole IC intensity. The relative luminosity between the two components in the SED can be quantified using the Compton dominance A_C = L_IC, peak/L_{synch, peak} (Finke 2013), where L_IC, peak is the luminosity at the IC peak and L_{synch, peak} the luminosity at the synchrotron peak. The excellent coverage in the X-rays and gamma rays allow these two quantities to be precisely determined. For MJD 57757 and MJD 57785, A_C ≈ 0.6 and A_C ≈ 0.5, respectively. On the other hand, for MJD 57813, A_C ≈ 0.2, which is roughly three times lower. We stress that these three nights are characterised by similar spectral properties both in X-ray and VHE (see Tables 3 and 4). This points towards a simple difference in the relative strengths of both spectral components without spectral changes. For MJD 57840, we find A_C ≈ 0.2.

We adopted a stationary one-zone SSC model (Ghisellini & Maraschi 1996; Tavecchio et al. 1998; Krawczynski et al. 2004) to interpret the four simultaneous SEDs. Within this scenario, the SED from the infrared is dominated by a single emitting zone. The low-energy component is due to synchrotron radiation by relativistic electrons, while the high-energy component originates from IC scattering on the synchrotron photons by the exact same electron population. This simple leptonic scenario was already applied to Mrk 421 and successfully described the SED on a wide range of flux states (Finke et al. 2008; Aleksić et al. 2012; Baloković et al. 2016). Compared to most of the published studies, we benefited from an excellent coverage in the hard X-ray band. The combination of optical, Swift-UVOT, Swift-XRT and NuSTAR observations fully encompasses the synchrotron component over ∼4 orders of magnitude in frequency, which brings stronger constraints on the electron energy distribution (EED).

In this work, we assumed an emitting zone consisting of a spherical blob with radius R′ (primed quantities refer to quantities in the plasma reference frame) that is moving downstream of the jet with a bulk Lorentz factor Γ_b. The blob is homogeneously filled with electrons and is embedded in a homogeneous magnetic field B′. We assume a jet axis at an angle Θ = 1/Γ_b. The advantage of this assumption is to reduce the number of degrees of freedom because, in this configuration, the Doppler factor δ becomes equal to Γ_b. For simplicity, the EED was characterised by a broken power-law (BPL), which is a very common parametrisation in blazar modelling:

$$\begin{array}{c}\hfill N^{\prime }({\gamma }^{\prime })=\{\begin{array}{c}{N}_0^{\prime }{\gamma }^{\prime -{\alpha }_1},{\gamma }_{\min }^{\prime }<{\gamma }^{\prime }<{\gamma }_{\mathrm{br}}^{\prime }\hfill \\ N_0^{\prime }{\gamma }_{\mathrm{br}}^{\prime {\alpha }_2-{\alpha }_1}{\gamma }^{\prime -{\alpha }_2},{\gamma }_{\mathrm{br}}^{\prime }<{\gamma }^{\prime }<{\gamma }_{\max }^{\prime }\hfill \end{array}\end{array}$$(3)

where $N_0^{\prime }$ is a normalisation constant. The corresponding electron energy density is given by $U_{\mathrm{e}}^{\prime }$ (in [erg cm⁻³]). The dimensionless parameters ${\gamma }_{\text{min}}^{\prime }$, ${\gamma }_{\text{br}}^{\prime }$ and ${\gamma }_{\text{max}}^{\prime }$ are defined as the minimum, break and maximum Lorentz factor, respectively. Because of the rather poor constraint on ${\gamma }_{\text{min}}^{\prime }$ provided by the data, we fixed it to 10³ to decrease the numbers of degrees of freedom and to be consistent with previous works on Mrk 421 (Abdo et al. 2011). The parameters α₁ and α₂ are the index below and above the break Lorentz factor, respectively. We note that a pure power-law EED, which has the advantage of having fewer degrees of freedom, clearly fails at describing the SED, especially in the X-ray regime. Consequently, such a simpler functional form is not a viable solution, and the X-ray data require a significant break in the EED. Internal photoabsorption of the VHE photons caused by the interaction with synchrotron radiation was included following Finke et al. (2008). Synchrotron self-absorption occurring in the radio domain was implemented following the prescription in Longair (2011). We employed routines from the open-source software naima to compute the synchrotron and IC interactions¹² (Zabalza 2015). Here, as in all of the models presented in this study, we did not use a minimisation strategy, such as that provided in the naima package, but rather, we performed a simple “eye-ball fit” where we tried to find a parameterisation with sensible model parameters that agreed with the broadband SED.

Prior to fitting the model, we fixed the Doppler factor, δ, to δ = 25. This is representative of the value typically found for HBLs (Tavecchio et al. 2010). The size of the emitting region was constrained by the light crossing time R′≤δt_varc/(1 + z), where t_var is the variability timescale in the observer frame. Depending on the night, t_var ≈ 4 − 11 h (see Appendix B), which is the flux halving time seen in the orbit-wise NuSTAR light curves.

In general, the model is able to describe well the data from optical to VHE for all the observing epochs. We remark that, on MJD 57813, a mismatch in the UV is apparent despite a good match from the X-rays to VHE. The UVOT measurement is under-predicted by a factor ∼2. Such a discrepancy may be attributed to the simplicity of the considered model, and could be better described by an EED with an additional break (two breaks instead of one). Alternatively, a higher γ_min (close to 10⁴) together with a harder α₁ also represents a viable solution to describe the UV-optical emission (not only for MJD 57813, but also for the other three nights), although this would be at the cost of an underestimation of the MeV-GeV SED measured by Fermi-LAT. Finally, one may also consider an additional region that contributes to the emission between the UV and the soft X-ray bands.

Moreover, as is typically found for blazars when considering one-zone SSC models, the radio flux is largely underestimated due to synchrotron self-absorption. Radio emission is most likely coming from regions outside of the inner jet, that have broader and complex morphology not included in our simple model (Giroletti et al. 2006; Lico et al. 2014). In our work, we mostly concentrate on emissions from optical to VHE, whose flux is believed to be dominated by a smaller region inside the jet.

The NuSTAR spectra extending to ≈40 keV allow the parameters of the higher-energy part of the EED (i.e., α₂ and ${\gamma }_{\text{max}}^{\prime }$) to be well-constrained. In HBLs, these two parameters are usually difficult to constrain and show large degeneracy with typical SEDs covering the synchrotron emission only up to the soft and medium X-ray band (see for example Ahnen et al. 2017b). We find a value of ∼10⁶ for ${\gamma }_{\text{max}}^{\prime }$ for each of the four models, with their precise values all lying within a factor of two of each other.

Regarding the break Lorentz factor ${\gamma }_{\text{br}}^{\prime }$, we find that the values are in rough agreement with the expected cooling break ${\gamma }_{\mathrm{b}\mathrm{r},\mathrm{e}\mathrm{x}\mathrm{p}}^{\mathrm{\prime }}$ that is obtained by balancing the synchrotron cooling timescale with the fiducial adiabatic cooling timescale of R′/c related to the expansion of the emitting zone, ${\gamma }_{\mathrm{br},\exp }^{\prime }=\frac{6\pi m_{\mathrm{e}}{c}^2}{{\sigma }_{\mathrm{T}}{B}^{\prime 2}{R}^{\prime }}$ (Tavecchio et al. 1998; Mücke & Protheroe 2001), where σ_T is the Thomson cross section and m_e the electron mass. Indeed, the modelling yields values that are at most a factor ≈3.5 away from ${\gamma }_{\mathrm{b}\mathrm{r},\mathrm{e}\mathrm{x}\mathrm{p}}^{\mathrm{\prime }}$. We note also that for most observing epochs the change of the index at the break, Δα = α₂ − α₁, is higher than the canonical synchrotron cooling break (Δα = 1) expected in models with a homogeneous emitting region (Longair 2011). For instance, we find Δα ≈ 1.6 and Δα ≈ 1.7 for MJD 57757 and MJD 57813, respectively. On MJD 57840, Δα ≈ 2. As mentioned above, the spectral shape of the EED is well-constrained by the data and fixing Δα = 1 during the fit of these particular epochs would result in a worse description of the SED. A large spectral break is in fact a recurrent result found in the modelling of SED and has been already reported for Mrk 421 (Tavecchio et al. 1998; Abdo et al. 2011; Baloković et al. 2016). This points to a more complex and less homogeneous emitting region. The canonical break condition may be loosened and larger values of Δα can be explained by considering velocity gradients in the jet (Ghisellini et al. 2005), for example.

As an alternative to the BPL model for the EED, we also investigated a log-parabolic model with a low-energy power-law branch (LPPL; Massaro et al. 2006; Tramacere et al. 2009). Several works have shown that such curved distributions may be produced through stochastic acceleration (Kardashev 1962; Tramacere et al. 2009, 2011) or via an energy-dependent acceleration probability process (EDAP; Massaro et al. 2004a). We found that a LPPL satisfactorily describes the four simultaneous MAGIC/NuSTAR/Swift observations, with very similar results as those obtained with a BPL EED. The curvature of the LPPL model was also in good agreement with the observed curvature in the synchrotron SED (as derived in Sect. 5). The data do not show a clear preference between the BPL and the LPPL model.

7. Broadband SED and modelling of the intriguing VHE flare on MJD 57788

The flare detected on MJD 57788 deserves special treatment. It is characterised by a strong VHE flux increase: the FACT observation on the day before results in a flux of ≈0.4 C.U., while a peak activity of ≈7 C.U. is measured during the outburst. The VHE flux decays during the night on ∼hour timescale (see Fig. 2). During the MAGIC observations the > 1 TeV flux is ≈3.5 C.U. The X-ray counterpart in the 0.3–10 keV band is much less evident. In fact, Fig. 9 shows that the > 1 TeV flux level is about ten times higher than other nights with a comparable 0.3–10 keV flux. The finely-binned Swift-BAT 15–50 keV light curve displays interesting features in Fig. 2. Simultaneous to the MAGIC highest state, the 15–50 keV flux is around 3.5 × 10⁻³ cm⁻² s⁻¹, while in the following hours (where we lack simultaneous VHE data) the flux is higher by a factor ∼3. This suggests a renewed flaring activity on ∼MJD 57788.7 after the decaying phase around MJD 57788. In the UV-optical and HE, no significant flux enhancement is visible around the outburst (see Fig. 1).

Figure 16 presents SEDs that illustrate the broadband evolution from MJD 57786 (2017 February 2) to MJD 57789 (2017 February 5). The SED on MJD 57786 was built around the closest Swift observation in time before the flare and is dubbed the “pre-flare” state. In the VHE band, we used the FACT data averaged between MJD 57786 and MJD 57787 to improve the statistics since no significant flux variability is visible between these two nights (see Fig. 1). The SED of the flare is plotted in red markers and uses the strictly simultaneous Swift-XRT/Swift-BAT/MAGIC observations. The butterfly from the Swift-BAT best-fit power-law model around MJD 57788.7 is shown, which corresponds to the renewed activity in hard X-rays. As shown in Fig. 2, no strictly simultaneous measurements are available at other wavebands. Finally, in blue markers we show the MAGIC SED measured on MJD 57789 (labelled as “post flare”), for which no simultaneous MWL data are available within less than a day. In each SED, the Fermi-LAT data were averaged over three days centred on the VHE observation. Optical data are not strictly simultaneous but are located less than a day away from the other wavebands. For comparison the typical state of Mrk 421 (Abdo et al. 2011) is shown in grey.

Fig. 16. Simultaneous broadband SEDs of MJD 57786 (pre-flare state), MJD 57788 (flare), and MJD 57789 (post-flare). Fermi-LAT spectral points are integrated over 3 days around the VHE measurements. VHE data with square black-filled markers are obtained from FACT observations, while X-ray data in black-filled markers are from Swift-BAT. The FACT SED for the pre-flare state was averaged from MJD 57786 to MJD 57787. The full black line is the two-zone model for the MJD 57788 flare state. The black dotted line represents the emission from the quiescent zone while the dashed line is the one from the flaring zone. The model parameters are listed in Table 6. Archival data representing the typical Mrk 421 state from Abdo et al. (2011) are shown in grey.

The softness of the pre-flare X-ray spectrum (Swift-XRT power-law index of ≈2.4) suggests a synchrotron peak frequency around 10¹⁷ Hz (≈0.4 keV). An evident hardening occurs on MJD 57788, and the power-law index measured by Swift-XRT is ≈1.9. The combined Swift-XRT and Swift-BAT measurements indicate that the synchrotron peak frequency is located at ∼10¹⁸ Hz (∼4 keV). At VHE, the spectrum is hard. The flux at the highest energies (> 1 TeV) is one order of magnitude higher than the pre-flare level. The IC peak frequency ν_IC is around 0.4–0.5 TeV, which is among the highest values for Mrk 421 (Aleksić et al. 2015b). The SEDs reveal no substantial increase in the UV-optical, where the flux varies by only 15–20%. The ≈2 × 10¹⁷ Hz (≈1 keV) flux remains at the level of the typical Mrk 421 activity (grey points in Fig. 16). Also, the flux at the low-energy peak frequency increases only by a factor ∼2, while that at the high-energy peak frequency seems to exhibit a factor ∼4 enhancement.

The question is what scenario could generate such a sudden change in the Compton dominance with respect to the pre-flare state assuming a one-zone SSC mechanism. In the synchrotron regime, the observed hardening without strong changes in the peak amplitude suggest modifications of spectral parameters from the high-energy part of the EED, i.e., parameters related to electrons emitting synchrotron photons significantly above ν_s. In particular, it suggests a hardening of α₂ and/or an increase of ${\gamma }_{\text{max}}^{\prime }$ or ${\gamma }_{\text{br}}^{\prime }$ to push ν_s to higher energies. The shift by ∼one order of magnitude in ν_s constrains ${\gamma }_{\text{br}}^{\prime }$ to change by a factor ∼3 at most (${\nu }_{\mathrm{s}}\propto {\gamma }_{\mathrm{br}}^{\prime 2}$), while the spectral hardening of ∼0.5 seen by Swift-XRT constrains a hardening of α₂ by ∼1 at most. However, mainly because of the onset of the Klein-Nishina suppression, these modifications are not sufficient to explain the large increase of the VHE flux. This is especially true at ≳1 TeV (i.e., deep in the Klein-Nishina regime) where the spectrum is particularly hard and the flux is about a factor of ten higher than in the pre-flare state. Furthermore, one should note that the high-energy electrons emitting synchrotron photons above ν_s radiate IC emission above ∼1–10 TeV, while the flare is already clearly visible at the low-energy end of the VHE spectra, around a few 100 GeV.

Modifications of the electron density would be needed to enhance the IC luminosity. This would require a rather fine tuning of other parameters (e.g., the magnetic field or the blob radius) in order to keep the UV-optical and ≈1 keV flux at an almost constant level. Since changes in the spectral shape of the EED are not sufficient to explain the hardness at VHE as explained above, a variation of δ would also be required to push the IC component to higher energies. We conclude that a one-zone SSC model demands too much fine-tuning to explain the transition in state from MJD 57786 to MJD 57788 and hence is not a good scenario.

Another scenario would be the appearance of a second emitting zone, in addition to the quiescent zone responsible for the pre-flare state. The variability timescale (∼hours) constrains the second emitting zone to be more compact, i.e., R′≲10¹⁵ cm. The absence of a UV-optical flare also suggests a more energetic and narrower EED. In this way, the second zone would suddenly dominate in the hard X-ray domain and in the VHE regime, and would remain subdominant in the rising segments of the two SED components, which would naturally explain the observations. This two-zone SSC scenario was employed to describe a flaring activity of Mrk 421 in 2010, as reported in Aleksić et al. (2015b).

We tried to perform a similar modelling by assuming two spherical emission zones spatially separated such that their respective radiation fields and particle populations do not interact with each other. In this scenario, the quiescent zone is responsible for the broadband emission during the pre-flare state on MJD 57786 (as well as the previous days), while the flaring zone dominates the hard X-ray and > GeV emission during the flare. The EED for the quiescent zone follows a BPL distribution, while in the case of the flaring zone a simple PL function is adopted. Just as in Sect. 6, the jet axis angle is Θ = 1/Γ_b such that Γ_b = δ. To limit the number of degrees of freedom, the bulk Lorentz factor of the flaring zone is equal to the one of the quiescent zone (Γ_b = 25). The resulting models are shown in Fig. 16. The corresponding parameters are listed in Table 6. Both the pre-flare and the flare are well described by the models, although the UV-optical slope during the flare seems slightly harder than predicted by the model. The data indeed provide an indication of a hardening in the UV-optical during the transition from the pre-flare to the flare state, yet the flux difference is only at the level of 15–20%, which remains a rather mild effect. As mentioned in the previous paragraph, the flaring zone is characterised by a more energetic and narrower EED with respect to the quiescent zone to ensure that it dominates only in the falling edge of two SED components. The lowest electron energies of this second region are constrained by the measured emission at ∼keV and ∼100 GeV, while the highest electron energies are constrained by the multi-TeV emission measured by MAGIC and the lack of strong X-ray emission above 20 keV, as indicated by the Swift-BAT measurement simultaneous to the XRT and MAGIC spectra (see Fig. 16). The model that we use to successfully describe the spectral measurements uses ${\gamma }_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{\prime }}=2\times {10}^4$ and ${\gamma }_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{\prime }}=6\times {10}^5$, and the radius of the flaring zone is R′ = 10¹⁵ cm, consistent with the rapid variability on hour timescales at VHE.

It can be seen from Fig. 2 that the SED on MJD 57788 lies in the decay phase of the flare. Figure 16 includes the Swift-BAT butterfly a few hours later on MJD 57788.7, during which a rapid renewed activity in the hard X-rays (15–50 keV) is observed. The best-fit power-law index is Γ = 1.82 ± 0.42, implying a synchrotron peak frequency roughly in the Swift-BAT pass-band. On MJD 57789, the measured Swift-BAT fluxes are very low (compatible with no signal), and the VHE flux is around 1.5 C.U. Given the fast variability and the lack of strictly simultaneous MWL data during these two epochs, we do not model the emission.

8. Discussion

8.1. VHE versus X-ray correlation

We confirm the existence of a significant correlation between the X-ray and the VHE emission, hence indicating a cospatial origin. Additionally, the fractional variability displays a two-peak shape with the highest variability both in the X-rays and VHE. These results suggest a single population of particles responsible for the emission in those bands. We remark that when using the simultaneous MAGIC/NuSTAR/Swift observations, F_var is significantly higher at VHE gamma rays than at X-rays, as clearly shown in Figs. 6 and 7, as well as in Figs. 8 and 15. These results differ from what was reported during the very low activity from the first months in the year 2013 (Baloković et al. 2016), but they are in agreement with the variability patterns measured during the extremely low activity in the observing campaign from the year 2016 (see Fig. 6 in Acciari et al. 2021). On the other hand, during large X-ray and VHE gamma-ray activity, the fractional variability at VHE gamma rays is equal to or lower than that measured at X-rays (Abeysekara et al. 2020; Acciari et al. 2020).

The correlations show a complex behaviour when different energy bands are compared. In fact, the data point to a correlation slope that depends on the spectral band (see Figs. 8 and 9). Thanks to the sensitive MAGIC measurements, we see a trend of higher slope with increasing VHE energies (for any given X-ray band). Such a behaviour was already noted by Acciari et al. (2020), when the source was strongly flaring. Here, the source was probed at a close-to-typical state. This pattern suggests a more direct relation of the > 1 TeV flux with the X-rays compared to the 0.2–1 TeV flux. Considering a SSC model and assuming a generic magnetic field strength of 0.1 G as well as a Doppler factor of 25, the ≳1 keV flux is dominated by electrons with Lorentz factors of γ′≈2 × 10⁵ (Tavecchio et al. 1998). Through IC processes, these same electrons dominate the ≳ TeV flux. Given the X-ray energies considered in this work, the observed trend is somewhat expected and is consistent with pure leptonic models.

For most of the energy combinations, the correlations are close to linear or even less during the MAGIC/NuSTAR/Swift observations. Nevertheless, we also find steep slopes. The correlation is quadratic or cubic in the VHE versus 0.3 − 3 keV case. Katarzyński et al. (2005) evaluated the VHE versus /X-ray correlation in BL Lac type objects within a one-zone SSC model using different scenarios for the parameters evolution. In most of the cases, the VHE versus X-ray correlation is expected to follow a linear trend rather than a quadratic trend (or cubic). It should be noted that a more-than-linear correlation may still be possible but under rather specific conditions and/or evolution of the parameters. One of the main reasons is that the IC processes responsible for the VHE emission are expected to occur in the Klein-Nishina regime, where the interaction cross section of electrons with seed photons above the synchrotron peak frequency ν_s is strongly suppressed. As a result, the VHE-emitting electrons do not efficiently up-scatter their self-produced synchrotron photons, but photons with lower energies, in the UV-optical. The criteria for having an IC luminosity and peak frequency affected by Klein-Nishina suppression is ${\gamma }_{\mathrm{br}}^{\prime }{\nu }_{\mathrm{s}}\ge \delta \frac{3m_{\mathrm{e}}{c}^2}{4h}$ (Tavecchio et al. 1998). The SEDs in Fig. 15 show synchrotron peak frequencies in the range of ν_s ∼ 10¹⁶ − 10¹⁷ Hz and the break Lorentz factors obtained from the modelling are ${\gamma }_{\mathrm{b}\mathrm{r}}^{\mathrm{\prime }}\approx 0.8\times {10}^4-2\times {10}^5$. Thus, the condition for having significant Klein-Nishina suppression is fulfilled at least for the three first observations. However, Katarzyński et al. (2005) argue that steeper correlation may also be obtained under a specific choice of the spectral bands, even taking into account Klein-Nishina effects. A quadratic and more-than-quadratic trend can occur when selecting X-ray energies close to or below the synchrotron peak, which is the case for the 0.3–3 keV band in Mrk 421. In conclusion, in the specific case of VHE versus 0.3–3 keV, the very steep correlation is possibly caused by the selection of the spectral bands. Moreover, Mrk 421 showed a displacement of the synchrotron peak frequency of two orders of magnitude throughout the campaign (see Fig. 5), and this may also contribute to the modifications in the correlation slopes.

Over the full MWL campaign, the best-fit slopes indicate a roughly cubic relation for > 1 TeV versus 0.3–3 keV, while an almost quadratic trend is seen for 0.2–1 TeV versus 0.3–3 keV and > 1 TeV versus 3–7 keV (see Fig. 9 and Table 1). In the latter situation (> 1 TeV versus 3–7 keV), for the reasons mentioned in the previous paragraph, a long-term quadratic correlation is rather surprising since we are considering synchrotron flux emitted above ν_s (at least for the vast majority of the measurements, see Fig. 12). The steep slope cannot be attributed to the choice of an X-ray band below the peak, as proposed above. Nonetheless, the quadratic slope is not representative of the average behaviour during the campaign. The best-fit slope is biased towards higher values by a small fraction of measurements (10%), including two of the simultaneous MAGIC/NuSTAR/Swift observations and the short flare on MJD 57788. These three observations show a significantly enhanced VHE activity compared to nights with similar X-ray flux, and they appear as clear outliers in Fig. 9. With the exception of these three days with large Compton dominance, the overall trend is well consistent with a ∼linear relation, as can be seen in Fig. 9.

The mentioned outliers suggest sudden emission phases (on timescales of about a few days) that show very different correlation patterns compared to more typical states. This translates into a particularly large scatter of the measurements in the correlation plots. Overall, the simple description of the correlation in the form $F_{\mathrm{V}\mathrm{H}\mathrm{E}}\propto F_{\mathrm{X}-\mathrm{r}\mathrm{a}\mathrm{y}}^x$ is overly simplistic, as indicated by the poor χ²/d.o.f. values. The reason for the scatter, that can also be seen at lower X-ray and VHE fluxes in Fig. 9, may be the consequence of several additional processes contributing to a dominant process driving the overall variability. A large scatter could also be explained by different emitting regions (from different parts of the jet) that dominate the X-rays and VHE bands for a given time period. It is important to stress that, in this study, we are not considering a specific flaring episode (that is most likely caused by the same emitting zone in the jet), but rather the nightly averaged flux spread over monthly timescales that possibly includes several small flaring episodes from different components of the jets. These regions may undergo changes in their physical environments, namely related to B′, R′, δ or $U_{\mathrm{e}}^{\prime }$, which could result in a disparate Compton dominance in the SED (Katarzyński & Walczewska 2010). For instance, a modification of the blob radius and/or electron density would modify the synchrotron photon target density for the SSC processes and increase or decrease the Compton dominance. A different Compton dominance naturally induces a large scatter in the VHE versus X-ray plot, giving rise to an apparent break in the correlation slopes. The modelling of the first three simultaneous MAGIC/NuSTAR/Swift observations (on MJD 57757, MJD 57785 and MJD 57813), which can indeed be explained by the variation of a few parameters related to the environment (see later in Sect. 8.3), is consistent with this hypothesis. On the other hand, while a one-zone SSC model is able to explain the MAGIC/NuSTAR/Swift observations, a deeper investigation of the outlier related to the MJD 57788 flare challenges the simple leptonic scenario (see Sect. 7).

8.2. UV-optical versus X-ray anti-correlation

We find an anti-correlation between UV-optical and X-ray at a significance level of above 3σ over the entire MWL campaign (see Sect. 5.2 for details). We also find that the strength and the significance of the anti-correlation is dominated by the observations during the first 40 days of the multi-instrument campaign (i.e., data taken before MJD 57760). This indicates that the anti-correlation between these two bands is not persistent; it varies over time and may become significant only over month timescales. A first indication of an anti-correlation was reported by Aleksić et al. (2015a), also for Mrk 421. The anti-correlation was marginally significant over a few months timescales, and over time lags ranging from 0 to −20 days, in agreement with what is found here. It was, however, not observed in other multi-instrument campaigns when Mrk 421 did not show strong flaring activity in X-rays (Baloković et al. 2016; Acciari et al. 2021). During a bright outburst of Mrk 421 in 2006 June, which lasted about two weeks, Lichti et al. (2008) reported a hint of positive correlation between the optical and X-ray fluxes. However, the particularly low significance of the correlation did not allow a conclusive claim. In 2010 February, during the brightest VHE gamma-ray activity of Mrk 421 detected to date, a marginally significant positive correlation of about 3σ was observed between the VHE gamma rays and the optical flux in the data taken on 2010 February 17, the day with the highest VHE gamma-ray flux (Abeysekara et al. 2020). Unfortunately, there were no X-ray observations simultaneous to the VHE observations, and hence the X-ray versus optical correlation could not be studied for that night. But owing to the very tight correlation that Mrk 421 always shows between the X-ray and the VHE gamma-ray bands, it is reasonable to assume that, on that night with outstandingly large VHE gamma-ray activity, the X-rays and optical emission may have been positively correlated. On the other hand, when considering the few-week-long dataset from 2010 February, no correlation is observed between the X-ray and optical fluxes, in contrast to what was reported in Lichti et al. (2008). Moreover, in 2013 April, during the second brightest VHE gamma-ray activity from Mrk 421 detected to date, which included several tens of hours of strictly simultaneous optical, X-ray and VHE gamma-ray observations, Mrk 421 did not show any correlation between X-ray and optical, despite showing large variability and a high degree of correlation between VHE and X-rays (Acciari et al. 2020).

To the best of our knowledge, Mrk 421 is the only BL Lac type object where an anti-correlation between these two segments of the SED has been observed to date. This second instance of this characteristic, which we presented and described here with better sampled observations than those presented in Aleksić et al. (2015a), suggests that this anti-correlation pattern, visible over a few month timescales, is a recurrent feature with possibly a real physical origin in the synchrotron emission of BL Lacs.

The first and direct implication is that the synchrotron emission is dominated by a cospatial population of electrons, at least from the UV-optical to the X-ray. One may interpret the anti-correlation as being due to a change in the efficiency of the electron cooling processes. Indeed, assuming a constant acceleration timescale, a stronger cooling shifts the overall EED towards lower energies, resulting in a reduction of the emission in the X-rays, while the UV-optical flux increases. The electron cooling (expected to be dominated by synchrotron radiation in case of HBL where the emitted synchrotron power is larger or at least comparable to the gamma-ray component, see Schlickeiser et al. 2010) may be increased by a stronger magnetic field B′. Alternatively, as suggested by Aleksić et al. (2015a), the EED may also be shifted towards lower or higher energies in case of changes in the acceleration efficiency.

The large variability of the synchrotron peak frequency during the campaign further supports changes in the cooling or acceleration mechanisms. For a given electron population, the synchrotron peak frequency evolves as ${\nu }_{\mathrm{s}}\propto {\gamma }_{\mathrm{br}}^{\prime 2}·B^{\prime }·\delta$ (Tavecchio et al. 1998). The variability of ν_s by about two orders of magnitude, that is reported in Sect. 5, disfavours δ as being the main driving parameter, as it would imply a variation by about two orders of magnitude, leading to unphysical and much larger values usually found in TeV BL Lac type objects (e.g., Tavecchio et al. 2010). Moreover, a change in δ would lead to a positive correlation, in contradiction to the observed anti-correlation pattern. Therefore, the large changes in the peak synchrotron frequency are likely produced by changes in B′ and ${\gamma }_{\text{br}}^{\prime }$, that are in turn parameters linked to acceleration and cooling mechanisms. However, while the shift in the synchrotron peak frequency is typical in Mrk 421, and particularly extreme in the 2017 observing campaign, the correlations between optical and X-rays are extremely rare. It indicates that the ν_s shifts are in most of the cases produced by the appearance and disappearance of emission at hard X-rays without affecting substantially the low-energy emission (e.g., optical and below). This could be produced, for instance, by the acceleration and cooling of the highest-energy electrons without a substantial change in the lowest-energy electrons, or perhaps by the time-variable contribution of an additional component that dominates the emission at the hard X-rays (e.g., Aleksić et al. 2015b). Only in a few cases (e.g., during the first 40 days of this campaign, as shown in Fig. F.2) there seems to be a tight relation between the optical-emitting electrons and the X-ray-emitting electrons, which could be caused by a full shift of the entire EED to higher or lower energies. As shown in Fig. 11, the synchrotron peak frequency ν_s did increase continuously by about two orders of magnitude (from 10¹⁶ Hz to 10¹⁸ Hz) during the first ∼40 days of the campaign, when the optical vs X-ray anti-correlation is most evident. Such a large increase in ν_s cannot be produced by a simple change of the magnetic field, since this would imply an extreme variation by two orders of magnitude of B′. Hence, under the assumption that a single zone is responsible for both the optical and X-ray emission, one must have the contribution from acceleration processes (possibly in shocks and turbulence) and cooling processes (e.g. synchrotron, inverse-Compton, adiabatic cooling) to shape an EED whose effective radiation yields the continuous increase in X-ray emission (by a factor of ∼10) and in ν_s (by a factor of ∼100), while decreasing the overall emission at UV-optical frequencies (by a factor of ∼2).

8.3. Interpretation of the emission during the simultaneous MAGIC/NuSTAR/Swift observations

The broadband emission during the simultaneous MAGIC/NuSTAR/Swift observations reveals intriguing behaviours. The most striking feature is visible when comparing the first three observations: while the synchrotron flux level is very similar, the corresponding Compton dominances A_C are significantly different. Namely, on MJD 57757 and MJD 57785 A_C is about three times higher with respect to MJD 57813. In fact, MJD 57757 and MJD 57785 appear as clear outliers in the correlation plots of Fig. 9. We note that the X-ray and VHE spectral parameters between those nights are very similar as shown in Tables 3 and 4. Within one-zone SSC models, this points towards an electron population with similar spectral characteristics. Hence, the difference in the Compton dominance A_C does not seem to be caused by acceleration and/or cooling mechanisms, as this would inevitably impact the particle distribution and result in a variation of the spectral properties.

We interpreted the simultaneous MAGIC/NuSTAR/Swift observations within a simple SSC model, which successfully describes the data. Benefiting from a wide energy coverage in the X-rays thanks to the combined NuSTAR/Swift-XRT data, the EED is rather well constrained. According to the expectations mentioned above, the EED spectral parameters are comparable between the first three nights. Furthermore, the magnetic field is almost constant between those nights, and B′≈0.06 − 0.07 G. In the Klein-Nishina regime, as is most likely the case here (see Sect. 8.1), the ratio of the peak frequencies of the two emission components relates as $\frac{{\nu }_{\mathrm{s}}}{{\nu }_{\mathrm{IC}}^2}\propto \frac{B^{\prime }}{\delta }$. Under the assumption of a roughly constant Doppler factor, the absence of significant variability in ν_s and ν_IC thus supports the constant behaviour of the magnetic field. The divergence in A_C is reconciled by a change in the emission zone radius R′ and the electron energy density $U_{\mathrm{e}}^{\prime }$: on MJD 57813, R′ is increased by a factor 1.65 and $U_{\mathrm{e}}^{\prime }$ is reduced by a factor ≈4 compared to MJD 57757 and MJD 57785. The lower electron density reduces the target synchrotron photon field, resulting in a reduction of the IC flux. It is important to note that the total number of electrons in the emitting zone, which is proportional to $U_{\mathrm{e}}^{\mathrm{\prime }}{R}^{\mathrm{\prime }3}$, is almost constant between the three nights. In this sense, the difference in the modelling parameters is dominated by a simple increase of the radius, which could happen due to a natural expansion of the emitting region, in the absence of sufficient pressure to constrain it to a given physical size. Under the latter hypothesis of an adiabatic expansion (without significant particle loss) and that the same emitting blob is responsible for the emission during the MAGIC/NuSTAR/Swift observations, it is then difficult to attribute the break in the EED to the cooling of the electrons (estimated by equating the adiabatic expansion timescale with the synchrotron losses, see Sect. 6). Indeed, given that the MAGIC/NuSTAR/Swift observations are separated by ≈30 days, the adiabatic expansion must also happen on similar timescales, which would result in a cooling break located at much lower Lorentz factors than the values of ${\gamma }_{\text{br}}^{\prime }$ found in the modelling. Adiabatic expansion on ∼weekly-monthly timescales may still be possible if one assumes that the emitting blob has a size similar to the cross-sectional radius of the jet and that the blob expansion is driven by its movement downstream in a conical jet structure. In this scenario, the adiabatic expansion timescale (co-moving frame) is given by $t_{\mathrm{ad}}^{\prime }=\frac{3}{2}\frac{R^{\prime }}{\theta {\mathrm{\Gamma }}_{\mathrm{b}}c}$, where θ is the half-opening angle of the jet (Moderski et al. 2003). Assuming Γ_b ∼ 10 and θ ∼ 0.1°, which is in agreement with results of parsec scale studies of blazar jets (Jorstad et al. 2005), one can obtain a value for $t_{\mathrm{a}\mathrm{d}}^{\mathrm{\prime }}$ which matches the weekly-monthly timescales measured in the observer’s frame. The break in the EED may thus be attributed to the intrinsic properties of the acceleration processes. On MJD 57840, the synchrotron component is clearly shifted towards lower energies, which points towards a change in the EED. Accordingly, the break Lorentz factor obtained is lower with respect to the other three MAGIC/NuSTAR/Swift nights, possibly indicating a decrease in the efficiency of the acceleration processes.

In all models we used a minimum Lorentz factor fixed to the fiducial value ${\gamma }_{\text{min}}^{\prime }={10}^3$. At such low energies, the dominant cooling process is most likely the adiabatic one, which may lead to a decrease of the electron energy below ${\gamma }_{\text{min}}^{\prime }={10}^3$. The fact that the same high value of ${\gamma }_{\text{min}}^{\prime }$ is used for the four models, that are separated in time by ≈30 days, suggests a continuous re-acceleration of the electrons in the emitting zone. Abdo et al. (2011) and Baloković et al. (2016) also modelled the broadband emission of Mrk 421 (for timescales of months and days) with ${\gamma }_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{\prime }}\sim {10}^3$, and stated the need for in-situ electron re-acceleration to fully explain these data sets.

The one-zone SSC modelling suggests a significant change of the Compton dominance caused by a modification of only a few parameters related to the jet environment. This provides a natural cause of the rather large scatter in the flux-flux correlation plots in Figs. 8 and 9.

The change in the relative strength of the two SED components could also be interpreted within the framework of external Compton models, in which an additional target photon field for the IC radiation is invoked. A possible approach would be the spine-layer model developed in Ghisellini et al. (2005). In the latter, the jet is assumed to be composed of a central part, the spine, surrounded by a layer. The bulk Lorentz factor of the spine (Γ_S) exceeds that of the layer (Γ_L). The synchrotron and IC emissions of the spine itself dominate over that of the layer. Due to their relative motions, the radiation of one region as observed in the reference frame of the other is seen boosted by a factor ∼Γ^′2, where Γ′=Γ_SΓ_L(1 − β_Sβ_L). This boosted radiation provides additional seed photons to each region and their respective IC luminosities are significantly enhanced, in particular the one from the spine. In this scenario, changes of Γ_L or the radiation density of the layer would therefore qualitatively explain the variability in the Compton dominance without modifying the synchrotron luminosity. One further advantage of the spine-layer model over the SSC model is being able to reach a system that is close to equipartition, i.e., $U_{\mathrm{e}}^{\mathrm{\prime }}\sim U_{\mathrm{B}}^{\mathrm{\prime }}$, as shown by Tavecchio & Ghisellini (2016). We note on the other hand that adding a second region (the layer) doubles the number of free parameters and the data at hand are not sufficient to constrain the model well (differently from the SSC one-zone model).

8.4. Theoretical interpretation of the intriguing VHE gamma-ray flaring activity on MJD 57788

The bright VHE flare on MJD 57788 reaches a peak flux of ≈7 C.U., but it only shows a moderate activity in the X-rays. At lower energies, the UV-optical flux exhibits low variability (≈20%) and no substantial increase with respect to the day prior to the flare. These MWL characteristics are difficult to explain with a one-zone SSC scenario, and lead us to propose a two-zone SSC scenario to describe the broadband behaviour. As described in Sect. 7, the two zones are spatially separated so that they do not interact with each other, with one region dominating the regular (quiescent or slowly varying) broadband emission, and another zone that dominates the flux enhancement observed at hard X-rays and VHE gamma rays.

Within this theoretical scenario, the non-variable UV-optical and MeV-GeV emission could be produced in a shock-in-jet component of relatively typical size dimensions (R′∼10¹⁶ cm), while the X-rays and VHE gamma rays, that show large variability during these days, could be dominated by the emission from a region that is smaller by about one order of magnitude (R′∼10¹⁵cm). This small region could originate in the base of the jet and produce an EED characterised by a very high minimum Lorentz factor ${\gamma }_{\text{min}}^{\prime }$ and with a narrow range of energies (${\gamma }_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{\prime }}/{\gamma }_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{\prime }}=30$). The value ${\gamma }_{\text{min}}^{\prime }$ needs to be above 10⁴ to avoid the overproduction of the keV flux and the 0.1 TeV flux. A similar two-zone SSC scenario was also used in Aleksić et al. (2015b) to describe the temporal evolution of the broadband emission of Mrk 421 during a 2-week flaring activity in 2010 March. This flaring activity contained several days with narrow SED peaks, and the two-zone scenario with a narrow EED (from ${\gamma }_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{\prime }}=3\times {10}^4$ to ${\gamma }_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{\prime }}=6\times {10}^5$) could describe the shape of the (narrow) X-ray and VHE bumps better than the one-zone scenario.

These narrow EEDs may arise through stochastic acceleration by energy exchanges with resonant Alfven waves in a turbulent medium yielding to the production of quasi-Maxwellian distributions of particle energies (Schlickeiser 1985; Stawarz & Petrosian 2008; Asano et al. 2014). An alternative way to produce narrow EEDs is through the emission resulting from an electromagnetic cascade initiated by electrons accelerated to a narrow range of energies in a magnetospheric vacuum gap, as proposed by various authors (Levinson & Rieger 2011; Ptitsyna & Neronov 2016; Wendel et al. 2021), and used successfully to explain an extremely narrow spectral component detected (at marginally significant level of 3–4σ) in the VHE spectrum of Mrk 501 during a large flaring activity in 2014 July (MAGIC Collaboration 2020b). Another scenario that could lead to narrow EEDs is magnetic reconnection, which has been invoked as an efficient particle acceleration process in AGN jets (Romanova & Lovelace 1992; Giannios et al. 2010; Giannios 2013). Blobs of magnetised plasma containing high-energy particles could be formed in the reconnection regions of jets and lead to high-energy emission, as proposed in the semi-analytic model from Petropoulou et al. (2016) and demonstrated in by dedicated particle-in-cell simulations (Christie et al. 2019, 2020). This theoretical framework was recently used to describe the temporal and spectral properties of the multiband flares that Mrk 421 showed in 2013 April (Acciari et al. 2020). Through magnetic reconnection, the dissipated magnetic energy would be converted into non-thermal particle energy, hence leading to a decrease in the magnetic field strength B′ for increasing gamma-ray activity, and hence leading to a ratio $U_B^{\mathrm{\prime }}/U_{\mathrm{e}}^{\mathrm{\prime }}$ as low as 10⁻³, which is needed to explain the measured broadband SED from MJD 57788.

The requirement for a narrow EED is linked to our assumption that the distribution follows a power law with index α₁ = 2 (see Table 6). It is a generic assumption common in blazar modelling that is supported by the prediction of magnetic reconnection (e.g., Sironi & Spitkovsky 2014; Sironi et al. 2015) and also by standard shock wave acceleration mechanisms (Crumley et al. 2019). As an alternative to the narrow EED in the flaring zone, a distribution of electrons spanning a wide range of Lorentz factors but identified by a hard power-law index less than 2 could reproduce the SED. In such a configuration, the flaring zone would remain subdominant in the UV-optical and the MeV-GeV regime even if the EED extends below ${\gamma }_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{\prime }}\sim {10}^3$. Sironi & Spitkovsky (2014) reported clear evidence that magnetic reconnection mechanisms can easily generate electron distributions following a power-law index harder than 2 in the regime where the magnetisation of the blob is larger than 10.

9. Summary

We have reported on a dense MWL observing campaign of Mrk 421 carried out between 2016 December and 2017 June. The MWL dataset comprises more than ten instruments, providing information from radio (with OVRO, Medicina and Metsähovi) to VHE gamma rays (with FACT and MAGIC), and including various instruments covering the optical and UV bands (e.g. GASP-WEBT and Swift-UVOT), X-ray bands (Swift-XRT and Swift-BAT and NuSTAR), and GeV gamma rays (with Fermi-LAT). Owing to the inclusion of NuSTAR data, we obtained a precise characterisation of the hard X-ray emission thought to originate from the high-energy tail of the same population of electrons dominating the VHE gamma-ray emission. This helped us to interpret the measurements within a standard SSC scenario.

The fractional variability versus energy showed the typical double-bump structure, observed in other campaigns. However, in contrast to many other MWL campaigns, when using strictly simultaneous data, the variability in the VHE gamma-ray domain was measured to be larger than that in the X-ray energy range. We found that the VHE gamma rays and X-rays are positively correlated with no time lag, but the strength and characteristics of the correlation change substantially across the various energy bands probed. The multi-instrument light curves and the broadband SEDs showed a large increase in the gamma-ray activity without a clear counterpart in the X-ray range. These orphan VHE gamma-ray activities, present in only a few of the observations (less than 10%), yielded quadratic and cubic dependencies in the VHE versus X-ray flux relations. Removing these few measurements, the relations become linear or sub-linear, which is in agreement with previous observations (Aleksić et al. 2015b; Acciari et al. 2020, 2021), and they are easier to explain with standard SSC models. We have shown that a one-zone SSC scenario with an expanding blob (change of the size of the emission blob without changing the number of electrons) can explain the decrease in the Compton dominance of the broadband SEDs during the MAGIC/NuSTAR/Swift observations, which show orphan gamma-ray activity.

The manuscript reports a substantial harder-when-brighter behaviour in both the gamma rays and X-rays, including displacements in the peak of the synchrotron bump by more than two orders of magnitude in energy. We also show an anti-correlation between UV-optical and X-ray at a significance above 3σ. This is the second time that such a trend is observed in Mrk 421, hence indicating that it is a repeating feature with a real underlying physical mechanism. This might be due to a change in the efficiency of the electron cooling or acceleration processes.

The manuscript also discusses an intriguing VHE flare observed on MJD 57788. In the synchrotron regime, the latter coincides with a clear hardening of the 0.3–10 keV spectrum with respect to the pre-flare state, but the flux around ∼1 keV remains close to the typical Mrk 421 quiescent activity. At VHE, the flare is strong and the flux above 1 TeV is roughly ten times higher than the typical quiescent activity. Within simple one-zone SSC models, this broadband behaviour is difficult to reproduce without fine tuning the parameters. We therefore interpret the flare as being caused by the appearance of a more compact second blob of highly-energetic electrons that span a relatively narrow range of energies (from ${\gamma }_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{\prime }}=3\times {10}^4$ to ${\gamma }_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{\prime }}=6\times {10}^5$) and could have been produced by stochastic acceleration, by magnetic reconnection, or by electron acceleration in the magnetospheric vacuum gap, close to the supermassive black hole.

Acknowledgments

The journal referee is gratefully acknowledged for a constructive list of remarks that helped us improve the contents and clarity of the manuscript. The MAGIC Collaboration would like to thank the Instituto de Astrofísica de Canarias for the excellent working conditions at the Observatorio del Roque de los Muchachos in La Palma. The financial support of the German BMBF, MPG and HGF; the Italian INFN and INAF; the Swiss National Fund SNF; the ERDF under the Spanish Ministerio de Ciencia e Innovación (MICINN) (FPA2017-87859-P, FPA2017-85668-P, FPA2017-82729-C6-5-R, FPA2017-90566-REDC, PID2019-104114RB-C31, PID2019-104114RB-C32, PID2019-105510GB-C31,PID2019-107847RB-C41, PID2019-107847RB-C42, PID2019-107847RB-C44, PID2019-107988GB-C22); the Indian Department of Atomic Energy; the Japanese ICRR, the University of Tokyo, JSPS, and MEXT; the Bulgarian Ministry of Education and Science, National RI Roadmap Project DO1-268/16.12.2019 and the Academy of Finland grant nr. 320045 is gratefully acknowledged. This work was also supported by the Spanish Centro de Excelencia “Severo Ochoa” SEV-2016-0588, SEV-2017-0709 and CEX2019-000920-S, and “María de Maeztu” CEX2019-000918-M, the Unidad de Excelencia “María de Maeztu” MDM-2015-0509-18-2 and the “la Caixa” Foundation (fellowship LCF/BQ/PI18/11630012) and by the CERCA program of the Generalitat de Catalunya; 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 important contributions from ETH Zurich grants ETH-10.08-2 and ETH-27.12-1 as well as the funding by the Swiss SNF and the German BMBF (Verbundforschung Astro- und Astroteilchenphysik) and HAP (Helmoltz Alliance for Astroparticle Physics) are gratefully acknowledged. Part of this work is supported by Deutsche Forschungsgemeinschaft (DFG) within the Collaborative Research Center SFB 876 “Providing Information by Resource-Constrained Analysis”, project C3. We are thankful for the very valuable contributions from E. Lorenz, D. Renker and G. Viertel during the early phase of the project. We thank the Instituto de Astrofísica de Canarias for allowing us to operate the telescope at the Observatorio del Roque de los Muchachos in La Palma, the Max-Planck-Institut für Physik for providing us with the mount of the former HEGRA CT3 telescope, and the MAGIC collaboration for their support. The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l’Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d’Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515. 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 also made use of the XRT Data Analysis Software (XRTDAS) developed under the responsibility of the ASI Science Data Center (ASDC), Italy. A.A.E and D.P acknowledge support from the Deutsche Forschungs gemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. M. B. acknowledges support from the YCAA Prize Postdoctoral Fellowship and from 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 publication makes use of data obtained at the Metsähovi Radio Observatory, operated by Aalto University in Finland. This research has made use of data from the OVRO 40-m monitoring program (Richards et al. 2011) which is supported in part by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G and NSF grants AST-0808050 and AST-1109911. I.A. acknowledges financial support from the Spanish “Ministerio de Ciencia e Innovación” (MCINN) through the “Center of Excellence Severo Ochoa” award for the Instituto de Astrofísica de Andalucía-CSIC (SEV-2017-0709). Acquisition and reduction of the MAPCAT data was supported in part by MICINN through grants AYA2016-80889-P and PID2019-107847RB-C44. The MAPCAT observations were carried out at the German-Spanish Calar Alto Observatory, which is jointly operated by Junta de Andalucía and Consejo Superior de Investigaciones Científicas. C.C. acknowledges support from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program under the grant agreement No 771282. This research was partially supported by the Bulgarian National Science Fund of the Ministry of Education and Science under grants KP-06-H28/3 (2018), KP-06-H38/4 (2019) and KP-06-KITAJ/2 (2020). We acknowledge support by Bulgarian National Science Fund under grant DN18-10/2017 and National RI Roadmap Projects DO1-277/16.12.2019 and DO1-268/16.12.2019 of the Ministry of Education and Science of the Republic of Bulgaria. This research was supported by the Ministry of Education, Science and Technological Development of the Republic of Serbia (contract No 451-03-68/2020-14/200002). G.D. acknowledges observing grant support from the Institute of Astronomy and Rozhen NAO BAS through the bilateral joint research project “Gaia Celestial Reference Frame (CRF) and fast variable astronomical objects” (2020–2022, head – G. Damljanovic). The BU group was supported in part by NASA Fermi guest investigator program grants 80NSSC19K1505 and 80NSSC20K1566. This study was based in part on observations conducted using the 1.8 m Perkins Telescope Observatory (PTO) in Arizona, which is owned and operated by Boston University. This article is partly based on observations made with the LCOGT Telescopes, one of whose nodes is located at the Observatorios de Canarias del IAC on the island of Tenerife in the Observatorio del Teide. This article is also based partly on data obtained with the STELLA robotic telescopes in Tenerife, an AIP facility jointly operated by AIP and IAC. The Abastumani team acknowledges financial support by the Shota Rustaveli National Science Foundation under contract FR-19-6174. Based on observations with the Medicina telescope operated by INAF – Istituto di Radioastronomia.

References

Appendix A: MAGIC analysis results

In this section, we report in detail the results of the MAGIC analysis. Table A.1 lists the nightly fluxes in the 0.2-1 TeV and > 1 TeV bands for each observation. The parameters of the log-parabolic (normalisation energy is 300 GeV) or power-law fits above 100 GeV are given for each day. A log-parabolic fit was applied in cases where it provides a better description of the data at a significance above 3σ in comparison to a power-law function. All fits were performed after correcting for the EBL absorption.

The two nights MJD 57757 and MJD 57789 show a relatively poor χ²/dof in comparison with the other nights, even with a log-parabolic model. Namely, the fits yield χ²/dof = 34.5/12 and χ²/dof = 28.1/14, which corresponds to a p-value of 6 × 10⁻⁴ and 10⁻², respectively. We found that a better fit could be obtained using a power-law with a exponential cutoff:

$$\begin{array}{c}\hfill \frac{\mathit{dN}}{\mathit{dE}}=f_0{\left(\frac{E}{E_0}\right)}^{-\alpha }exp(-E/E_c)\end{array}$$(A.1)

For MJD 57757, the best-fit parameters are f₀ = (1.23 ± 0.06) × 10⁻⁹ cm⁻² s⁻¹ TeV⁻¹, α = 1.86 ± 0.06, and E_C = 1.7 ± 0.3 TeV. The associated χ²/dof is χ²/dof = 25.6/12.

For MJD 57789, the best-fit parameters are f₀ = (1.04 ± 0.04) × 10⁻⁹ cm⁻² s⁻¹ TeV⁻¹, α = 1.87 ± 0.06, and E_C = 2.87 ± 0.6 TeV. The associated χ²/dof is χ²/dof = 23.5/14.

Appendix B: NuSTAR analysis results

In this section, we report in more detail on the orbit-wise light curves from the NuSTAR observations. They are plotted in Fig. B.1. Each light curve is incompatible with a constant flux. In order to quantify the variability timescale, we used the prescription of Zhang et al. (1999), which provides an estimate on the flux doubling or halving time t_1/2. For any consecutive flux measurement i and j, t_1/2 is defined as the minimum of:

$$\begin{array}{c}\hfill t_{1/2}^{i,j}=\left|\frac{F_j+F_i}{2}\frac{T_j-T_i}{F_j-F_i}\right|\end{array}$$(B.1)

where F_i, j is the flux and T_i, j is the time. The uncertainty was obtained by propagating the error of F_i, j. We computed t_1/2 in the 7-30 keV band as it shows the strongest variability. The resulting value for each day is shown in the respective plot from Fig. B.1. It ranges from ∼4 hrs to ∼11 hrs.

Table B.1 list the flux values along with the spectral parameters of the log-parabola fits.

Fig. B.1. NuSTAR orbit-wise 3-7 keV and 7-30 keV light curves of the four observations on MJD 57757, MJD 57785, MJD 57813 and MJD 57840. The flux doubling or halving time is indicated in the legend. It was computed based on the 7-30 keV flux using the prescription of Zhang et al. (1999).

Appendix C: Multi-wavelength light curves during the simultaneous MAGIC/NuSTAR/Swift observations

In this section, the MWL light curves during the four simultaneous MAGIC/NuSTAR/Swift observations are shown in Fig. C.1 to Fig. C.4. The MAGIC, NuSTAR, Swift fluxes are computed in identical 30-min time bins.

Fig. C.1. MWL light curves during the simultaneous MAGIC/NuSTAR/Swift observations on MJD 57757. The first two panels from the top are the MAGIC light curves in the 0.2-1 TeV and > 1 TeV bands with 30 min binning. The third and fourth panels from the top show the NuSTAR and Swift-XRT light curves in the 0.3-3 keV, 3-7 keV and 7-30 keV bands with 30-min binning. The bottom panel shows the R-band observations provided by the WEBT-GASP community.

Fig. C.2. Same description as in Fig. C.1, but for MJD 57785.

Fig. C.3. Same description as in Fig. C.1, but for MJD 57813.

Fig. C.4. Same description as in Fig. C.1, but for MJD 57840.

Appendix D: Swift-XRT analysis results

In this section, we report in detail the results of the Swift-XRT analysis. Table D.1 lists the fluxes in the 0.3-2 keV, 2-10 keV and 3-7 keV bands for each observation. The best-fit power-law index (Γ) as well as the best-fit parameters α and β from the log-parabolic fits (pivot energy of 1 keV) are also given.

Appendix E: Evaluation of the statistical significance of the UV-optical versus X-ray anti-correlation

In this section, we describe the details of the procedure used to assess the statistical significance of the UV-optical versus X-ray anti-correlation. First of all, the power spectral density (PSD) of the light curves were estimated (e.g., Max-Moerbeck et al. 2014). The PSD gives a measure of the strength of the variability as function of the temporal frequency. It is a powerful and widely used tool to characterise the temporal flux behaviour on a large range of timescales. For this study, we adopted the multiple fragments variance function (MFVF) from Kastendieck et al. (2011) to estimate the PSD. The MFVF method has the advantage of not relying on any interpolation and re-binning as required in other widely used PSD estimation method (e.g., the PSRESP method described in Uttley et al. 2002). Because our light curves can show strong variability (especially in the X-rays) and have a very irregular sampling as well as large gaps, applying some interpolation or re-binning might introduce additional systematic effects. We assumed here as a PSD model a simple power-law shape, i.e. P(ν)∝ν^−β. This simple parametrisation was found to describe well Mrk 421, as well as other blazar light curves on a broad energy range (e.g., Uttley et al. 2002; Chatterjee et al. 2008; Abdo et al. 2010b; Aleksić et al. 2015a). The parameter β usually ranges from 1 (referred to as pink noise) to 2 (referred to as red noise). We carried an estimation of β similarly to Nilsson et al. (2018). We summarise below the procedure.

Based on the assumed PSD model, light curves were simulated following the method described in Timmer & Koenig (1995). The light curves were simulated on timescales 100 times longer than the observations in order to take into account red-noise leakage. We then applied the exact same sampling patterns as the real observations to include the distortions of the PSD caused by the sampling. Furthermore, distortions due to aliasing effects were included by generating light curves with a time resolution matching the typical exposure time of the data (Uttley et al. 2002). Finally, the simulated light curves were rescaled to match the flux variance of the observations, and Gaussian noise corresponding to the measurement uncertainties was added. The simulated light curves can now be considered as realistic light curves.

For a fixed β, we simulated 5000 realistic light curves and the MFVF was computed for each of them. The MFVF was characterised down to a minimal temporal frequency f_min = 1/T up to a maximal frequency f_max = 1/Δt₀, where T is the total length of the light curve and Δt₀ is the shortest timescale variability probed. In our case, T lied between ≈150 days and ≈180 days depending on the light curves, while we fixed for all light curves Δt₀ = 1 day, which is the typical shortest cadence of the observations. We then binned the MFVF in seven frequency bins. For each frequency bin f_i, the MFVF probability density function p(β, f_i) was estimated using a Gaussian kernel density estimation. Finally, from the measured light curve, a log likelihood function was computed as follows:

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

where N is the number of frequency bins. The function p(β, f_i) relates to the probability of measuring a MFVF in a frequency bin f_i assuming a power-law index β for the PSD. The best-fit index β_fit maximises ℒ. For this, we scanned a range of β from 0.7 to 2.1, with 0.05 steps. The resulting best-fit indices for the optical (R-band), Swift-UVOT/W2, Swift-XRT(0.3-2 keV) and Swift-XRT(2-10 keV) light curves are summarised in Table E.1. The uncertainty on β_fit was estimated by performing a large number of fits on simulated light curves that have a true β equal to β_fit. The uncertainties were computed from the 68% containment around the mean of the distribution of the resulting best-fit indices.

The R-band and UV light curves yield β_fit values compatible within uncertainties. This is expected given the proximity in energy. Nilsson et al. (2018) found similar results based on an optical light curve spanning ≈10 years. In the two X-ray band, β_fit agree within uncertainties and the values are consistent with what Aleksić et al. (2015a) reported.

The significance of the DCF coefficients between two observed light curves can now be determined. With the same procedure described above, we simulated 2 × 10⁴ realistic and uncorrelated light curves for each energy band assuming the best-fit PSD index β_fit in the PSD model. The DCF coefficients between the two sets of simulated light curves were computed, and at each time lag the distribution of the DCF coefficients was drawn to extract the 2σ and 3σ confidence intervals.

In Fig. 14, the significance bands show a spike at a time lag of ∼30 days in the case of the R-band versus X-ray. We found that this feature is caused by the sampling of the two light curves and is not a statistical artefact. Around this time lag, due to the sampling of the two light curves, the DCF is highly dominated by measurements performed on MJD 57757 and MJD 57785 (corresponding to two of the MAGIC/NuSTAR/Swift simultaneous observations). On those two dates the optical measurements (from GASP-WEBT) have a dense and fine binning on a sub-hour timescale, allowing many pairs contributing to the overall DCF (see Fig. C.1 and Fig. C.3). On such a timescale, the temporal properties of the optical flux lead to an almost constant flux behaviour, but also potentially to an overall decrease or increase of the flux. When correlated to the Swift-XRT data, there is therefore a sizable probability of detecting a spurious (anti-)correlation around this date, which would then dominate the total DCF.

Appendix F: Investigation of the temporal behaviour of the UV-optical versus X-ray anti-correlation

From a visual inspection of Fig. 1, UV-optical versus X-ray anti-correlation is mainly visible for a period of roughly 40 days between MJD 57720 to MJD 57760. In this section, we investigate the impact of this time period on the correlated behaviour derived over the entire MWL campaign.

We repeated the study described in Sect. 5 and Appendix F and computed the DCF significance by ignoring data before MJD 57760. The results are shown in Fig. F.1. While a hint at the level of 2 − 3σ remains in R-band versus X-rays, the significance is below 2σ in the UV versus X-rays. We repeated again the exercise, but this time ignoring data after MJD 57760. The results are shown in Fig F.2. For this time period, the significance is always above 3σ, and the DCF value at the peak is also higher.

Fig. F.1. DCF between X-ray (0.3-2 keV and 2-10 keV) and UV-optical (Swift-UVW2 and R-band) when removing data before MJD 57760. The blue- and red-dashed-lines indicate the 2σ and 3σ confidence intervals estimated from the Monte Carlo simulations, as described in the text.

Fig. F.2. Same description as for Fig. F.1, but this time data after MJD 57760 are removed.

All Tables

All Figures

Fig. 1. MWL light curves between MJD 57720 and MJD 57918. From top to bottom: MAGIC (0.2–1 TeV & > 1 TeV) and FACT (> 0.95 TeV), Fermi-LAT (0.2–2 GeV & 2–300 GeV), Swift-BAT (15–50 keV), NuSTAR and Swift-XRT (3–7 keV & 7–30 keV), Swift-XRT (0.3–2 keV & 2–10 keV), UV-optical (Swift-UVW1/UVM2/UVW2 & R-band), OVRO (15 GHz), Metsähovi (37 GHz) and Medicina (8 GHz and 24 GHz). The dotted blue horizontal line in the first two panels from the top shows the 1 C.U. flux in the > 1 TeV & 0.2–1 TeV bands, respectively. The Fermi-LAT data have a 3-day binning. The Swift-BAT light curve is computed in 3-day and 6-day binnings, as well as 1-day binning around the flare. The vertical black dashed lines highlight the dates that show an enhanced VHE activity and appear as outliers in the VHE versus X-ray correlation plots (see Sect. 4).

In the text

Fig. 2. Zoom on the MWL light curves around the VHE flare on MJD 57788. Two first panels from the top show the FACT (> 0.95 TeV) and the MAGIC (> 1 TeV & 0.2–1 TeV) light curves. The FACT fluxes are computed on 20-min binning while the MAGIC light curve is nightly averaged. On MJD 57788 and MJD 57789 the MAGIC exposures are ∼40 min and ∼100 min respectively, and are depicted with the width of the horizontal bars of the markers. Third panel from the top: the Swift-BAT light curve with 3-h binning. Two panels from the bottom: the Swift-XRT (0.3–2 keV & 3–7 keV & 2–10 keV) and Swift-UVW1/UVM2/UVW2 fluxes. The horizontal orange lines show the average flux over the entire campaign in each energy regime. Top panel: the orange line is the average from the MAGIC > 1 TeV light curve. In the two lowest panels, the orange line depicts the average in the Swift-XRT 0.3–2 keV and Swift-UVW2 bands.

In the text

	Fig. 3. Hardness ratio (F>1 TeV/F0.2−1 TeV) versus F>1 TeV (top figure) and versus F0.2−1 TeV (lower figure) from the MAGIC measurements. The Pearson’s coefficients and the corresponding significance (following the prescription of Press et al. 2007) are given in both plots.
In the text

	Fig. 4. Log-parabola photon index α versus the 3–7 keV flux of the orbit-wise NuSTAR observations. The parameter α is fitted after fixing spectral curvature β = 0.22 in the log-parabolic model. Each day is plotted with a different colour. The black line represents a linear fit, while the grey area is its uncertainty. The resulting slope is −0.46 ± 0.01.
In the text

Fig. 5. Log-parabola photon index α versus the 3–7 keV flux from the Swift-XRT observations. The log-parabolic fits are performed after fixing the spectral curvature β = 0.16. The red cross represents the flare on MJD 57788 seen at VHE energies. The black line represent a linear fit, while the grey area (hardly visible in the plot) is its uncertainty. The slope of α versus the 3–7 keV flux is −0.64 ± 0.01. Hollow markers depict fits with a p-value below 5 × 10−2, which are not considered in the linear fit.

In the text

Fig. 6. Fractional variability Fvar obtained from the light curves shown in Fig. 1. MAGIC, FACT, Swift-XRT, Swift-UVOT, R-band and radio fluxes are nightly binned. Fermi-LAT and Swift-BAT fluxes have a 3-day binning. Results from each instrument are plotted in different colours. The filled markers include all data. The hollow markers include VHE and Swift data lying within a time window of 4 h from each other.

In the text

	Fig. 7. Fractional variability Fvar from the simultaneous MAGIC-NuSTAR observations performed on the four nights MJD 57757, MJD 57785, MJD 57813 and MJD 57840 (2017 January 4, 2017 February 1, 2017 March 1 and 2017 March 28). The Fvar values were computed with the fluxes determined on 30-min time bins that are reported in Appendix C.
In the text

Fig. 8. VHE flux versus X-ray flux during the simultaneous MAGIC/NuSTAR/Swift observations. Flux points are computed over time bins of 30 min. From left to right: the X-ray energy bands in the panels are 0.3–3 keV, 3–7 keV and 7–30 keV. In the upper panels, the MAGIC fluxes are in the > 1 TeV band, while in the lower panels they are in the 0.2–1 TeV band. Data points from Swift-XRT are shown with black-filled markers. Fluxes corresponding to each day are plotted in a different colour. The results of the Pearson coefficients (and its significance in σ), DCF and the slope of the linear fits in the log-log plane (with the χ2/d.o.f. values) are indicated in each of the corresponding subplots.

In the text

Fig. 9. VHE flux versus X-ray flux over the MWL campaign using MAGIC and Swift-XRT data. Data are nightly binned, and only pairs of measurement within 4 h are considered. MAGIC fluxes are in the > 1 TeV band (top panels) and in the 0.2–1 TeV band (lower panels). Swift-XRT fluxes are computed in the 0.3–3 keV (left panels) and 3–7 keV bands (right panels). The open red diamond markers highlight measurements on MJD 57757 and MJD 57785 (i.e., the dates of the first two simultaneous MAGIC/NuSTAR/Swift observations). The VHE flare (on MJD 57788) is plotted in open red circle marker. Black lines depict linear functions (in the log-log plane) with slopes x = 1,  2,  3. In the 0.2–1 TeV versus 3–7 keV panel, because of the dynamical range of the data, the linear function with slope x = 3 is replaced with one with slope of x = 0.5. The Pearson coefficients, DCFs and the slopes from linear fits in the log-log plane are shown in Table 1.

In the text

Fig. 10. VHE flux versus X-ray flux over the MWL campaign using FACT and Swift-XRT data. Data are nightly binned, and only pairs of measurement within 4 h are considered. FACT fluxes are computed above 0.95 TeV and only measurements with a signal-to-noise ratio > 2 are taken into account. Swift-XRT fluxes are computed in the 0.3–3 keV (left panel) and 3–7 keV bands (right panel). The open blue markers highlight measurements on MJD 57757 and MJD 57785 (i.e., the dates of the first two simultaneous MAGIC/NuSTAR/Swift) and the VHE flare (on MJD 57788). Black lines depict linear functions (in the log-log plane) with slopes x = 1 and x = 2. The Pearson coefficients, DCF and slope are shown in Table 1. The Pearson coefficients, DCF and slopes from linear fits in the log-log plane are shown in Table 2.

In the text

Fig. 11. Evolution of the synchroton peak frequency νs versus time throughout the MWL campaign. The peak frequencies are obtained by fitting a log-parabola to simultaneous Swift-XRT, Swift-UVOT and Swift-BAT observations. Blue square markers indicate fits that result in a p-value above 10−3. Data points in violet triangles correspond to a p-value between 10−3 and 10−4, and red diamond markers have a p-value lower than 10−4.

In the text

Fig. 12. Synchrotron peak frequencies throughout the MWL campaign obtained from log-parabola fits to simultaneous Swift-XRT, Swift-UVOT and Swift-BAT observations. Blue square markers indicate fits that result in a p-value above 10−3. Data points in violet triangles correspond to a p-value between 10−3 and 10−4, and red diamond markers have a p-value lower than 10−4. The black dotted line is a linear fit to the data and the obtained slope is indicated on the plot. The dashed areas represent the ranges of archival synchrotron peak frequencies, and are taken from Fig. 12 of Baloković et al. (2016). The respective references are given in the legend.

In the text

Fig. 13. νsFν(νs) versus the curvature parameter b throughout the MWL campaign obtained from log-parabola fits to simultaneous Swift-XRT, Swift-UVOT and Swift-BAT observations. As in Fig. 12, blue square markers indicate fits that result in a p-value above 10−3. Data points in violet triangles correspond to a p-value between 10−3 and 10−4, and red diamond markers have a p-value lower than 10−4. The black dotted line is a linear fit to the data, and the obtained slope is indicated on the plot.

In the text

	Fig. 14. DCF between X-ray (0.3–2 keV and 2–10 keV) and UV-optical (Swift-UVW2 and R-band) during the MWL campaign. The blue- and red-dashed-lines indicate the 2σ and 3σ confidence intervals estimated from the Monte Carlo simulations, as described in the text.
In the text

Fig. 15. SED of the four simultaneous MAGIC/NuSTAR/Swift observations (MJD 57757, MJD 57785, MJD 57813, MJD 57840). The data points and obtained one-zone SSC models are plotted with brown markers and orange lines, respectively. The parameters of the one-zone SSC models are shown in Table 5. For comparison purposes, the one-zone SSC model of the first pointing (MJD 57757) is plotted with black dotted lines in all the other subplots. The Compton dominance parameter AC determined from the model is indicated on each panel. Archival data representing the typical Mrk 421 state from Abdo et al. (2011) are shown in grey.

In the text

Fig. 16. Simultaneous broadband SEDs of MJD 57786 (pre-flare state), MJD 57788 (flare), and MJD 57789 (post-flare). Fermi-LAT spectral points are integrated over 3 days around the VHE measurements. VHE data with square black-filled markers are obtained from FACT observations, while X-ray data in black-filled markers are from Swift-BAT. The FACT SED for the pre-flare state was averaged from MJD 57786 to MJD 57787. The full black line is the two-zone model for the MJD 57788 flare state. The black dotted line represents the emission from the quiescent zone while the dashed line is the one from the flaring zone. The model parameters are listed in Table 6. Archival data representing the typical Mrk 421 state from Abdo et al. (2011) are shown in grey.

In the text

	Fig. B.1. NuSTAR orbit-wise 3-7 keV and 7-30 keV light curves of the four observations on MJD 57757, MJD 57785, MJD 57813 and MJD 57840. The flux doubling or halving time is indicated in the legend. It was computed based on the 7-30 keV flux using the prescription of Zhang et al. (1999).
In the text

Fig. C.1. MWL light curves during the simultaneous MAGIC/NuSTAR/Swift observations on MJD 57757. The first two panels from the top are the MAGIC light curves in the 0.2-1 TeV and > 1 TeV bands with 30 min binning. The third and fourth panels from the top show the NuSTAR and Swift-XRT light curves in the 0.3-3 keV, 3-7 keV and 7-30 keV bands with 30-min binning. The bottom panel shows the R-band observations provided by the WEBT-GASP community.

In the text

	Fig. C.2. Same description as in Fig. C.1, but for MJD 57785.
In the text

	Fig. C.3. Same description as in Fig. C.1, but for MJD 57813.
In the text

	Fig. C.4. Same description as in Fig. C.1, but for MJD 57840.
In the text

	Fig. F.1. DCF between X-ray (0.3-2 keV and 2-10 keV) and UV-optical (Swift-UVW2 and R-band) when removing data before MJD 57760. The blue- and red-dashed-lines indicate the 2σ and 3σ confidence intervals estimated from the Monte Carlo simulations, as described in the text.
In the text

	Fig. F.2. Same description as for Fig. F.1, but this time data after MJD 57760 are removed.
In the text