Rapid Variability of Mrk 421 During Extreme Flaring as Seen Through the Eyes of XMM-Newton

By studying the variability of blazars across the electromagnetic spectrum, it is possible to resolve the underlying processes responsible for rapid flux increases, so-called flares. We report on an extremely bright X-ray flare in the high-peaked BL Lacertae object Mrk 421 that occurred simultaneously with enhanced 𝛾 -ray activity detected at very high energies (VHE) by FACT on 2019 June 9. We triggered an observation with XMM-Newton , which observed the source quasi-continuously for 25 hours. We find that the source was in the brightest state ever observed using XMM-Newton , reaching a flux of 2 . 8 × 10 − 9 erg cm − 2 s − 1 over an energy range of 0.3 – 10 keV. We perform a spectral and timing analysis to reveal the mechanisms of particle acceleration and to search for the shortest source-intrinsic time scales. Mrk 421 exhibits the typical harder-when-brighter behaviour throughout the observation and shows a clock-wise hysteresis pattern, which indicates that the cooling dominates over the acceleration process. While the X-ray emission in different sub-bands is highly correlated, we can exclude large time lags as the computed zDCFs are consistent with a zero lag. We find rapid variability on time scales of 1 ks for the 0.3 – 10 keV band and down to 300 s in the hard X-ray band (4 – 10 keV). Taking these time scales into account, we discuss different models to explain the observed X-ray flare, and find that a plasmoid-dominated magnetic reconnection process is able to describe our observation best


INTRODUCTION
Among active galactic nuclei (AGN), blazars show the strongest and the most rapid variability (e.g., Stein et al. 1976;Wagner & Witzel 1995;Ulrich et al. 1997).Their luminosity and extreme properties ★ E-mail: gokus@wustl.edu,dr.andrea.gokus@gmail.comoriginate from a collimated outflow, the so-called jet, of relativistically moving particles close to the line of sight (Antonucci 1993;Urry & Padovani 1995), and their emission is Doppler boosted towards the observer.Their characteristic broadband spectral energy distribution has a double-hump structure and sub-classes of blazars are broadly defined by their overall luminosity and the peak position of the first hump, which lies somewhere between near-IR to X-ray energies.The emission of the low-energy hump is typically strongly dominated by leptonic synchrotron radiation, while the high-energy hump, which peaks at -ray energies, can be explained via leptonic Inverse Compton emission (e.g., Ghisellini et al. 1985;Maraschi et al. 1992) as well as hadronic processes (e.g., Mannheim 1993;Atoyan & Dermer 2001;Mücke & Protheroe 2001;Mücke et al. 2003;Kelner & Aharonian 2008;Böttcher et al. 2013).
Blazar variability is observed across the entire electromagnetic spectrum and from sub-hour time scales to years (e.g., Stein et al. 1976;Wagner & Witzel 1995;Ulrich et al. 1997;Pian 2002;Bhatta & Dhital 2020).Many coordinated multi-wavelength studies, partially with a focus on high-amplitude flux increases, that is flares, have been conducted to shed light on the emission regions within the jet as well as the underlying particle acceleration processes.The fastest variability is typically observed during bright flares and can happen on sub-hour time scales.The most rapid flux variations, on timescales below 10 minutes, have been observed at -ray energies for five blazars so far: PKS 2155−304 (Aharonian et al. 2007), Mrk 501 (Albert et al. 2007b), IC 301 (Aleksić et al. 2014), 3C 279 (Ackermann et al. 2016), and CTA 102 (Meyer et al. 2019).These fluctuations challenge established models such as the shock-in-jet model (Marscher & Gear 1985) or the spine-sheath-model (Henri & Pelletier 1991;Ghisellini et al. 2005), and motivated the development of new scenarios, such as the minĳet-in-a-jet model (e.g., Ghisellini & Tavecchio 2008;Giannios et al. 2009;Narayan & Piran 2012;Zhang et al. 2021), with slightly different underlying processes, or the jet interacting with a star or a gas cloud (Barkov et al. 2012;Heil & Zacharias 2020).Alternatively, stochastic perturbations in the particle acceleration process are able to explain the observed, lognormal flux distributions, which are characteristic for blazars (e.g., Sinha et al. 2018;Burd et al. 2021).
Due to its proximity ( = 0.031; Ulrich et al. 1975) and brightness, which enables obtaining high signal-to-noise data at all wavelengths, the high-synchrotron peaked (HSP) source Markarian 421 (Mrk 421) has been extensively monitored.In 1992, it was detected as an extragalactic TeV-emitter with the Whipple telescope (Punch et al. 1992), and since then has exhibited frequent and bright -ray flares with the flux exceeding that of the Crab nebula at TeV energies by a factor of three or more (e.g., Albert et al. 2007a; Bartoli et al. 2016), even over ten Crab flux units in three instances (Acciari et al. 2014).While a large amount of data has been gathered during different activity states of the source, as well as transitioning states, the behaviour of Mrk 421 has not been fully understood because of its complexity.The source shows signatures of rapid variability, i.e., on sub-hour time scales at VHE -ray energies (e.g., Gaidos et al. 1996;Abeysekara et al. 2020;Acciari et al. 2020;MAGIC Collaboration et al. 2021) as well as at X-ray energies (e.g., Abeysekara et al. 2017;Chatterjee et al. 2021;Noel et al. 2022).Its TeV -ray and X-ray emission have been found to be usually correlated during flares (e.g., Macomb et al. 1995;Buckley et al. 1996;Albert et al. 2007a;Fossati et al. 2008;Acciari et al. 2011;Cao & Wang 2013), but also during quiescent phases (Horan et al. 2009;Aleksić et al. 2015;Acciari et al. 2021).However, Mrk 421 has also shown indications for so-called orphan flares at VHE -rays (Błażejowski et al. 2005;Acciari et al. 2011;Fraĳa et al. 2015), which lack correlated X-ray flares and whose origins are unclear.In a study covering a time span of 5.5 years and densely sampled data in the radio, X-ray, and -ray (GeV and TeV) energy range, Sliusar et al. (2019) report 29 flares that occurred simultaneously in the X-ray and TeV band, and two orphan TeV flares.On average, this would result in a -ray flare occurring roughly every 65 days.A 14-year long study with the Whipple telescope, covering a time span from 1995 to 2009, determined an average flux of 0.446±0.008Crab flux units for Mrk 421 (Acciari et al. 2014).(E th ∼730 GeV) was obtained through the long-term monitoring campaign by FACT, while monitoring of the source in the X-ray band (0.3-10 keV) was done with Swift -XRT (dark blue), which is taken from Gokus et al. (2022).The XMM-Newton data (teal) is included as well and covers the majority of the X-ray flare.The color scale for the -ray light curve bins displays the detection significance for Mrk 421 for each observation.The inlay of the color scale does not cover any -ray data since there was no visibility of Mrk 421 after MJD 58649 for FACT.Note that the uncertainties for the X-ray flux are so well constrained that error bars are not visible.
In this work, we present an X-ray observation of Mrk 421 that was performed with XMM-Newton (Jansen et al. 2001), and coincident with high VHE -ray activity in June 2019.In order to obtain such a simultaneous observation, we use the long-term monitoring campaign of nightly observations with the First G-APD Cherenkov Telescope (FACT; Anderhub et al. 2013) for selected blazars (Dorner et al. 2022;Arbet-Engels et al. 2021).Our campaign, which has been designed to trigger several multi-wavelength instruments during a VHE -ray flare is described in Kreikenbohm (2019), and the details on the follow-up observations during June 2019 are reported in Gokus et al. (2022) and Gokus (2022).A light curve displaying the monitoring of Mrk 421 in the -ray and X-ray band in June 2019, the time of the extreme X-ray flare presented in this work, is shown in Fig. 1.The goal was to trigger multi-wavelength follow-up observations for the occurrence of a VHE -ray flux above 2 Crab Units (CU).This condition for sending the requests for target-of-opportunity observations was fulfilled on 2019 June 9 based on the quick-look analysis, which is described in Dorner et al. (2015).However, this analysis does not include a correction for the dependence of the threshold on the zenith distance and ambient light (see Arbet-Engels et al. 2021, for details), and a subsequent analysis found that the -ray flux (E th ∼730 GeV) was slightly lower.The data were analysed with the Modular Analysis and Reconstruction Software (MARS; Bretz & Dorner 2010), using the background suppression method described in Beck et al. (2019).Good quality data were selected based on the cosmic-ray rate (Beck et al. 2019;Hildebrand et al. 2017) calculating the artificial trigger rate R750 above a threshold of 750 DAC1 -counts.A corrected rate, 750 cor , is calculated for an evolving zenith distance (Mahlke et al. 2017;Bretz 2019).Because the cosmic-ray rate changes from season to season, a monthly reference value 750 ref was derived, and a cut for good quality data for 0.93 < 750 cor /750 ref < 1.3 was applied.The highest VHE -ray flux is measured on MJD 58645.9, which is coincident with the X-ray flare, and our analysis yields 1.63 ± 0.35 CU, or (5.1 ± 1.1) × 10 −11 ph cm −2 s −1 .Since the VHE -ray flux did not quite reach the trigger threshold of 2 CU (see Fig. 1), we do not refer to the -ray activity as a flare, but as high VHE -ray activity.
Here, we focus on the variability analysis of the X-ray data obtained simultaneously with the -ray activity.With a peak 0.3-3 keV flux of 1.6 × 10 −9 erg cm −2 s −1 , the flux during our observation was the highest that has ever been detected with XMM-Newton, and is close to the brightest 0.3-3 keV fluxes ever measured for this source during a flare in 2013 April (2.2 × 10 −9 erg cm −2 s −1 ; Acciari et al. 2020) and a flaring episode in 2018 January (∼ 2.8 × 10 −9 erg cm −2 s −1 , estimated from the reported flux of 5 × 10 −9 erg cm −2 s −1 in the 0.3-10 keV band; Kapanadze et al. 2020) 2 .
This paper is structured as follows: In Section 2, we describe the data extraction of the XMM-Newton observation.We elaborate on the spectral analysis of these data in Section 3, and perform a timing analysis in Section 4. We discuss the implications of our findings and conclude in Section 5.

OBSERVATION WITH XMM-Newton
The XMM-Newton observation presented in this work was performed between 2019 June 10 19:11:27 UTC and 2019 June 11 20:26:42 UTC (ObsID: 0845000901) as a target-of-opportunity observation after enhanced VHE -ray activity was observed by FACT at 2019 June 9 22:30 UTC.We concentrate on the observations with the European Photon Imaging Camera (EPIC) pn-camera (Strüder et al. 2001).Observations done with EPIC MOS-camera (Turner et al. 2001) were heavily piled-up and could not be used for scientific analysis.Two changes of observing mode split the observation in three sections.The first 10 ks were observed in the EPIC-pn Burst mode, followed by 35 ks in Timing mode, and lastly 38 ks in the Burst mode again.Between these sections, gaps of roughly 3 ks length exist.The reason for the second switch of observation mode, i.e., from Timing to Burst mode, was the rise in flux by the source, which meant that the observation was no longer feasible in Timing mode.We extract the EPIC-pn data following the standard methods of the XMM-Newton Science Analysis System (SAS, Version 20.0).The part of the observation taken in Timing mode was affected by pileup.We mitigate the pile-up effects by ignoring the central five CCD rows around the peak of the photon distribution.The extraction region for the Timing mode data is therefore 20 ≤ RAWX ≤ 35 && 41 ≤ RAWX ≤ 56.The extraction region for the Burst mode data covers 40 rows, 20 ≤ RAWX ≤ 60.For these data, we apply an additional cut at RAWY≤ 140 in order to avoid pile-up that can occur during the readout (see Kirsch et al. 2006).
For the full exposure time, we generate light curves with different time binnings (100 s and 1 ks) and for different energy ranges for the variability analysis.The 100 s-binned light curve for both the full energy range (0.3-10 keV) as well as the soft (0.3-1 keV), medium (1-4 keV) and hard (4-10 keV) band are shown in Fig. 2. In addition, we extract a light curve with time binnings of 10 ms for the data that were taken in the Timing mode, which are used for computing the power spectral density (see Sect. 4.2).In the first 10 ks of the observation, strong variability is visible in the hardest energy band at around 5 ks and 7 ks from the start of the observation.We extract the light curve for energies above 10 keV in the background region in order to check for the same pattern being present at hard X-rays, which is a hint for background flaring occurring during the observation.We find that an excess of up to five times the average background count rate is present at the same time as the peaks in the 4-10 keV light curve within the first section of the hard X-ray light curve (see lightgray data in Fig. 2).Hence, we exclude these data from our study.In addition to light curves, we extract time-resolved spectra by splitting the data into sections of 2 ks exposure each.
Because data were obtained both in the EPIC-pn Timing and Burst modes, the uncertainties on the count rates and spectral parameters (flux and photon index) are different between the two modes.The reason for this is the different readout methods, and therefore live time, of the detector in these science modes.While it is 99.5% for the Timing mode, in the Burst mode the detector collects data only during 3% of the time.Accordingly, spectral parameters and count rates are less constrained when data were taken for the Burst mode.

SPECTRAL VARIABILITY
Using this XMM-Newton data set, it is possible to probe quasicontinuous X-ray emission on time scales longer than 20 hrs and resolve the processes in the jet.Ignoring the first 10 ks due to potential background flaring, our observation covers about 74 ks including one gap of ∼ 3 ks length, which restricts our dataset from being strictly continuous.Nonetheless, the overall duration of the observation enables us to reveal significant spectral changes on time scales of a few hours.As a first step, we want to use a model-independent approach to broadly investigate how the distribution of photons in different energy bands changes without any assumption of the underlying process with which the data can be modelled.A hardnessintensity diagram (HID) is an ideal tool to quantify spectral changes during the observation, and it is defined as (e.g., Park et al. 2006) where  and  are the number of counts in the high-energy (harder) and low-energy (softer) band, respectively.By definition, the hardness value can range from −1, meaning that the low-energy band has a overwhelming amount of photons, to +1, for which significantly more photons are found in the high-energy band.Naturally, more photons are detected at lower energies, and the effective area of the EPIC favors the detection of soft photons, too, but changing emission processes or absorption can affect this order.While this quantity does not rely on any physical assumptions, it is biased by the chosen energy ranges that are compared to one another.The uncertainty on the hardness is computed via Gaussian error propagation.
In Figure 3, we show the light curve for the 0.3-10 keV band in 1 ks binning, which is the same binning we choose also to create the HIDs. Figure 4 displays the HIDs for the medium (1-4 keV) versus soft (0.3-1 keV) band on the left, and the hard (4-10 keV) versus soft band on the right.As expected, the values in both HIDs are negative due to the larger amount of photon counts in the soft band compared to the medium and hard band.With a typical photon index > 2 during an average flux state (e.g., Abeysekara et al. 2017;Markowitz et al. 2022), the X-ray emission from Mrk 421 is rather soft, similar to that of most HSP sources.In the first part of the light curve, the flux constantly rises and the harder-when-brighter trend is visible in both HIDs.During the transition from the first to the second section of the light curve (Timing to Burst mode), the count rate seems to increase without a change of the ratio between soft and medium band (left HID), while the hard vs. soft HID indicates a slight decrease of photons in the 4-10 keV band compared to those in the 0.3-1 keV band.The harder-when-brighter trend continues for the second part of the light curve, and above 1000 cts s −1 , the behaviour in the medium vs. soft band appears to be following the same trend for both increasing and decreasing count rates.However, the trend seen in the HID for the hard vs. soft band seems to be slightly different: above 1000 cts s −1 , the hardness ratio stays at the same level during the increase of the count rate.Then, the hardness decreases when the count rate falls back to the level of ∼1050 counts s −1 after the peak flux.This behaviour creates a counter-clockwise rotation in the HID, which indicates a change to an improved cooling efficiency since the ratio of hard to soft photons stays the same during the flux increase, and becomes softer than before while the flux drops after reaching its maximum.
As a next step, we want to apply a model to the extracted spectra and determine physical parameters, that is, the slope of the underlying power law and the flux of the source.Given the extreme brightness of Mrk 421 during the observation, spectra with a relatively short exposure of 2 ks suffice to constrain both parameters well.We rebin the 2 ks-exposure spectra such that each bin has a minimum signal-tonoise ratio of 5, excluding bins below 0.3 keV and above 10 keV.Each spectrum is individually fitted with an absorbed pegged power law (tbabs*pegpwrlw).For the absorption in the interstellar medium, we adopt the vern cross sections (Verner et al. 1996) and wilm abundances (Wilms et al. 2000).The absorption parameter is set to the Galactic hydrogren column value,  H,Gal = 1.33 × 10 20 cm −2  The grey data points are the hardness ratios of the 100-s binned light curve and display the scatter of hardness ratios on very short time scales, while the colored data points are obtained from the 1 ks-binned time series.The colors match the time bins of the light curve displayed in Fig. 3.The sections that were taken in Timing and Burst mode, respectively, are separated by the dashed gray line.(HI4PI Collaboration, et al. 2016), and kept frozen to avoid correlation with the photon index.The best-fit results are listed in Table 1.The value of reduced  2 for the power law model is very close to unity, indicating that the model is sufficient to describe the data sets.However, sometimes a log-parabola is found to provide a better model of the X-ray spectrum of Mrk 421 (e.g., Kapanadze et al. 2016).Here, we choose a simple power law on purpose: the curvature of the parabola () is correlated with the photon index (), and in order to describe the spectrum accurately, usually both parameters are kept free in the modelling process.In this work, we want to track the change of the spectral slope through the photon index, which is why a power law model is our favoured option to avoid cross-correlation of parameters.We choose to visualise the time-resolved spectral changes with a so-called hysteresis curve (Kirk et al. 1998).In order to create this hysteresis curve, we use the best-fit results of the extracted spectra, which are able to verify sub-hour changes while at the same time providing adequately constrained flux and photon index values (see Table 1).During this 2019 June observation, Mrk 421 presents a hard spectrum compared to spectra taken earlier during low flux states (∼ 1-5×10 −10 erg cm −2 s −1 with power law photon indices between 2.4 and 2.8, e.g., Baloković et al. 2016).This comes as no surprise as Table 1.Best fit results for the absorbed power law fit to the 2 ks spectra of this XMM-Newton observation of Mrk 421.The flux is given for the energy range from 0.3 to 10 keV.The results are used to display the hysteresis curve in Fig. 5. the 'harder-when-brighter' behaviour is a well established correlation for HSP sources, and in particular Mrk 421 (e.g., Zhang et al. 2019).
A distribution of spectral indices obtained via Swift observations from 2015 until 2018 peaks at roughly 2; however, this distribution has wide tails extending to 1.6 and 2.9 (Kapanadze et al. 2020) demonstrating the strong and frequent changes in the steepness of the X-ray spectrum exhibited by this source.Compared to the spectral parameters obtained by this study for similar flux levels as reported here, we find a slightly harder spectrum.However, since the data analysed by Kapanadze et al. (2020) was taken with Swift and our data set was obtained through XMM-Newton, we cannot exclude the possibility that potential systematic differences between the two instruments may contribute to the discrepancy.The resulting hysteresis curve and an accompanying light curve in 2 ks-binning are shown in Fig. 5.During the majority of the observing time, the X-ray flux of Mrk 421 increases.In addition to the increasing flux, the spectrum becomes significantly harder, hence the movement from the lower right corner to the upper left corner in the hysteresis plot.This harder-when-brighter behaviour is also seen in the model-independent HID.The overall spectral hardening takes place within less than 1 day, with the continuous increase being roughly 65 ks long.For the first section of the light curve, the flux steadily increases most of the time, while the photon index initially decreases, before behaving in the opposite way: increasing while the flux gets slightly lower (e.g., the small, intermediate peak seen at ∼35 ks).Here, no loop can be observed.In the second section of the light curve, the flux peaks at ∼70 ks, before it visibly drops.For this behaviour, a clear clockwise rotation can be observed in the hysteresis curve, which indicates that the cooling time scale dominates the underlying physical processes, as the high-energy particles cool efficiently and the accelerated particle population interacts with the low-energy photons later, i.e., showing a 'soft lag' (c.f.Kirk et al. 1998).The 'soft lag'-behaviour is a common occurrence in the X-ray emission of HSP-type blazars (e.g., Falcone et al. 2004;Wang et al. 2018), and has also been seen in the optical (Agarwal et al. 2021).For Mrk 421, both clockwise and anti-clockwise rotating hysteresis curves have been found in the past: data taken with XMM-Newton (0.5-10 keV) as well as RXTE/PCA (2-60 keV) revealed opposing behaviour, in some cases for observations spaced just a few days apart (Ravasio et al. 2004;Cui 2004;Abeysekara et al. 2017).The dissimilar evolution of flux and spectral shapes indicates that we might see different types of flares, or different stages within flaring events.However, observing times often only cover a relatively short period of time (usually hours) compared to the length of an activity phase, and since strictly continuous monitoring over several days with an instrument like XMM-Newton is currently not available, it is only possible to work with data that covers only a small segment, which is biased towards flaring states due to target-of-opportunity based observation strategies.

TIME SERIES ANALYSIS
While we are able to track down spectral variations down to a few kilo-seconds for this data set, the time resolution of the EPIC-pn onboard XMM-Newton enables us to study the variability properties of sources down to minutes or even seconds.We utilize the light curves to derive energy-resolved values for the fractional variability amplitude (Sect.4.1), to search the shortest source-intrinsic variability signatures (Sect.4.2), and to determine whether or not a lag is present between the different X-ray energy bands (Sect.4.3).

Fractional variability amplitude
The fractional variability amplitude ( var ) is a quantity commonly used to derive a measure of the variability in a time series normalized to a source's average flux. var is a useful tool for unevenly sampled data and has therefore often been used for the analysis of data obtained in multi-wavelength campaigns.Here, we use it to provide values that can be compared to other work that has been done on similar data sets.It was defined by Vaughan et al. (2003) as where the excess variance  nxs (Nandra et al. 1997;Edelson et al. 2002)  Figure 5. XMM-Newton light curve (left) in 2 ks binning and corresponding hysteresis curve (right).The flux is given for the 0.3-10 keV band.The color-coding is used to enable an easier connection of the light curve to the data points of the hysteresis curve.The sections that were taken in Timing and Burst mode, respectively, are separated by the dashed gray line, and the smaller uncertainties for flux and photon indices measured in the Timing mode are due to the significantly higher live-time and therefore photon collecting ability of this observing mode.During both the Timing and the Burst mode, the increase and decrease of a (localized) peak in the light curve can be observed at 35 ks and 70 ks, respectively.While for the peak in the Timing mode section of the light curve the path on the hysteresis curve is identical for rise and decay, a hysteresis loop pattern is visible for the peak in the Burst mode section.

err(𝐹 var
where  is the number of all measurements.We compute  var for the 100s-binned light curve and list the results in Table 2.For Mrk 421, the fractional variability amplitude at X-ray energies ranges typically from 40% to 80% when covering time ranges of several days to weeks or months (e.g., Giebels et al. 2007;Horan et al. 2009;Baloković et al. 2016;Ahnen et al. 2016;Acciari et al. 2020), however, our data cover a time range of less than a day, resulting in a lower fractional variability amplitude.A systematic study by Schleicher et al. (2019), which is based on FACT monitoring data of blazars, has shown that the size of the analysed data set and the time span which a data set covers can strongly influence the fractional variability amplitude obtained.Therefore, we compare our results to a study of 25 observations (average duration of ∼ 38 ks) of Mrk 421 with XMM-Newton, which covered all observations between 2000 and 2017 with a minimum duration of 10 ks (Noel et al. 2022).The data set we use to determine  var contains a total of 726 data points.The value of  var for this observation is significantly above the average value of  var derived from the systematic study.While Noel et al. (2022) have potentially found a slight positive correlation between the variability of a source and its flux exists, the length of an observation also influences fractional variability towards higher values.This connection indicates the presence and dominance of red noise in the light curves of Mrk 421.
In addition, the values for  var exhibit a clearly visible energy dependence, with the fractional variability amplitude being almost double for the hard vs. soft band.Prior studies have found the same trend that the value of  var increases with X-ray energy for Mrk 421

Search for shortest time scales
By searching for the shortest time scales on which a source exhibits variability, we can obtain knowledge of the underlying processes and potential structure within the emission region.Our first approach to search for short-time variability is computing the flux variability time scale, also known as e-folding time.Burbidge et al. (1974) defined an estimate of the flux-normalized (or weighted) variability time scale via with Δ being the time interval between two selected flux measurements,  1 and  2 .An uncertainty for Δ var can be derived via standard error propagation (see also Bhatta et al. 2018), which yields Δ. (5) We compute  var for all combinations of two data points in the 100s-binned light curve.From those results, we select the smallest value that fulfils the condition of | 1 −  2 | > Δ 1 + Δ 2 , which ensures that the variation lies outside of the flux uncertainties.The resulting values are listed in Table 2.
We find a value of  var of ∼ 1 ks for the full energy range, which translates to variability on a time scale of roughly 16 minutes.The value of  var for the sub-bands is even smaller, with a fastest variability time scale of ∼5 min at the highest energies.
In previous studies of Mrk 421 that have used  var as a tool, similar time scales for the full energy band have been found (Yan et al. 2018;Chatterjee et al. 2021;Noel et al. 2022).Noel et al. (2022) do not find any clear correlation of  var with the flux state of Mrk 421.Another study which did not include data from Mrk 421, however, reported a tentative correlation of  var with the mean flux, for which the fastest variability occurs during the lowest flux state of a source (Bhatta et al. 2018).This study analysed NuSTAR data of 13 blazars, both of the FSRQ and the BL Lac type.Since the observation of Mrk 421 presented in this work reveals rapid variability during a very high X-ray flux state, our findings do not match their proposed correlation.
A more refined method to study variability in any given time series is using a discrete Fourier transformation   in order to compute a power spectral density (PSD), which is defined as with  *  being the complex conjugate.Since the Timing mode is able to provide a large live-time in combination with a high time resolution of less than 1 ms, we will use the data taken in this mode, and use light curves with different binning to create a combined PSD that spans a large frequency range from 10 −4 to > 10 Hz.We use a light curve with a binning of 10 ms, and in order to stay consistent, we create light curves with a binning of 1 s and 10 s from this initial light curve.We ignore all values for which either the rate is nan, or the fractional exposure is lower than 90%.The count rate is steadily rising in the light curve, and therefore does not initially fulfill the condition for stationarity, which is necessary for a timing analysis via a PSD.In order to account for this behaviour, we fit the 10 s light curve with a linear function,  + , and subtract a straight line with parameters from the best fit results from each light curve.Then, we add the average count rate determined by a linear fit of each light curve, respectively, in order to get the appropriate normalization.
We compute periodograms; to reduce the scatter inherent in single periodograms, we perform averaging across multiple periodograms derived from multiple light curve segments.We optimise the segment length to cover the maximum possible time range that doesn't include any gaps, while also producing a reasonable amount of segments.For the 10 ms, 1 s, and 10 s binned light curves we produce 4084, 419, and 6 segments that we create periodograms for and average over, respectively.In addition, we apply a frequency rebinning of Δ  /  = 0.4.We perform a Monte Carlo simulation for each PSD to estimate the level of the Poisson noise contribution to the PSD based on the measured light curves by taking into account the statistical properties, the linear increase, and occurrences of gaps, and compute PSDs in the same way as has been done for the real data.We use the results from the simulation to subtract the Poisson noise from the PSDs in order to determine up to which frequency, and therefore time scales, variability is present in the light curve of Mrk 421.In Figure 6, we show the computed PSDs in the top panel, and the Poisson noise-subtracted PSDs in the other two panels, with the bottom one showing the PSDs multiplied by frequency.The upper 1 uncertainty envelope is indicated by the gray line, and makes visible where power remains above the Poisson noise level.As can be seen clearly in the middle and lower panel in Fig. 6, no variability remains at frequencies above ∼ 10 −3 Hz.This translates to a time scale of about 1 ks, which is in agreement with our calculation of  var , as well as other studies (Abeysekara et al. 2017;Yan et al. 2018;Chatterjee et al. 2021;Noel et al. 2022), and seems to be a characteristic time scale for Mrk 421.

Search for energy-related time lags
The detection of a presence, or a lack, of time lags between two energy bands can reveal information about acceleration vs. cooling time scales during flares.Studies searching for lags between the soft and hard X-ray band have yielded different results as Mrk 421 seems to present an inconsistent behaviour with showing no lag during some times, but also exhibiting either a positive or negative lag that can be as large as a few kiloseconds at other times (e.g., Takahashi et al. 1996;Fossati et al. 2000;Takahashi et al. 2000;Tanihata et al. 2001;Ravasio et al. 2004;Lichti et al. 2008;Arbet-Engels et al. 2021;Markowitz et al. 2022).A tool that is commonly used for lag determination is the discrete correlation function (DCF; Edelson & 4.4 ± 5.9 Krolik 1988).It is ideal for irregularly sampled data sets and data that show gaps, such as our XMM-Newton data set.In this work, we are using the z-transformed DCF by Alexander (1997Alexander ( , 2013)), which alleviates biases associated with the classic DCF.We use the 100sbinned light curves to resolve potential short-time lags, and compute the zDCFs for the combinations soft vs. medium band, soft vs. hard band, and medium vs. hard band.Our maximum lag is the length of the light curve (75 ks) with a binning of 100 s for the zDCF curve and we ensure that the binning is centered around a lag value of zero seconds.The resulting three zDCFs are shown in Fig. 7.The correlation for the two combinations, 0.3-1 keV vs. 4-10 keV (Soft/Hard) and 1-4 keV vs. 4-10 keV (Medium/Hard), appears more asymmetric towards positive lag times than the correlation of the 0.3-1 keV vs. 1-4 keV (Soft/Medium) bands.In order to determine the lag, we use two techniques: (1) determining the lag  peak at the position of the maximum correlation factor  max , and ( 2) averaging over all points with correlation coefficients > 0.8  max in order to obtain a central time lag value  cent .The results for both approaches are listed in Table 3.The error on  peak is calculated via the zDCF using a maximum likelihood method.We note that for the Soft/Hard and Medium/Hard correlation, the curves reach briefly above the threshold of > 0.8  max at lags of ∼ 23 ks and ∼ 37 ks, respectively.We exclude these parts in order to have a continuous portion of the zDCF to fit.We consistently find that peak lag values are slightly positive (suggesting harder photons leading softer ones), although values are formally consistent with zero lag.Across all light curves, we can safely exclude that soft photons lead hard photons by ≥ 0.3-0.4ks.We also exclude that hard photons lead soft ones by values ≥ 1.0, 2.1, and 1.7 ks in the Medium/Soft, Hard/Soft, and Hard/Medium correlations, respectively.The tendency towards a soft lag in relation to the hard X-ray band is also visible in the hysteresis curve in Fig. 5, where a clock-wise loop appears that describes the last 30 ks of the observation.
In general, Mrk 421 seems to exhibit both soft, hard, and also no lags in the X-ray band.While the X-ray emission in different energy bands is typically correlated, this correlation is often found to be consistent with a time lag of zero, as was reported based on data obtained by NuSTAR (3-10 keV vs. 10-79 keV; Pandey et al. 2017), XMM-Newton (0.3-2 keV vs. 2-10 keV Noel et al. 2022), Chandra (0.3-2 keV vs. 2-10 keV; Aggrawal et al. 2018), and ASTROSAT (0.7-7 keV vs. 7-20 keV; Das & Chatterjee 2023).In other cases it was found that the lower-energy emission leads the higher-energy emission (i.e., a hard lag) based on data obtained with BeppoSAX (0.1-1.5 keV vs. 3.5-10 keV / 0.1-2 keV vs. 2-10 keV; Fossati et al. 2000;Zhang 2002, respectively), ASTROSAT (0.6-0.8 keV vs. several sub-bands reaching up to 7 keV; Markowitz et al. 2022), and INTEGRAL (20-40 keV vs. 40-100 keV;Lichti et al. 2008), with the tendency to be observed during flaring events.A soft lag was reported based on data from an XMM-Newton observation (0.6-0.8 keV vs. 4-10 keV Zhang et al. 2010), as well as an ASCA observation (0.5-1 keV vs. 2-7.5 keV Takahashi et al. 1996).Some studies found a time lag being present only for a part of the analysed light curve (Brinkmann et al. 2005;Markowitz et al. 2022), indicating that the involved emission processes can change relatively quickly.Due to the uncertainties on our derived values, which are determined by the standard deviation of  cent , our results are consistent with a zero time lag, even though they tentatively hint at a soft lag.

DISCUSSION & CONCLUSIONS
We analysed a target-of-opportunity XMM-Newton observation of Mrk 421 during a very bright X-ray flare that occurred simultaneously to VHE -ray activity in June 2019.Our spectral analysis reveals typical harder-when-brighter behaviour.From the clock-wise hysteresis loop and a tentative soft lag, which is found by employing the z-transformed discrete correlation function, we conclude that the electron cooling dominated over the acceleration process during the flare.Our timing analysis determines a fractional variability amplitude above the average level of single pointed observations and rapid variability down to time scales of ∼ 300 s at the highest energies, and 1 ks in the full energy range.The fastest variability timescales that we find can give as an estimation of the size of the emission region for such rapid variability.In our calculations, we use  = 0.031, and  BH = 1.9 × 10 8 M ⊙ (Woo & Urry 2002), which yields a Schwarzschild radius of 1  S = 2.9 × 10 13 cm.If we naively assume that the emission region spreads over the entire width of the jet, and that the jet employs a conical shape (which has been found to be true for most HSP and radio-loud quasars up to several kpc, e.g., Pushkarev et al. 2017), we can compute the distance of the emission region from the jet base.For an opening angle , this distance of the jet base is given by with  r being the upper limit for the size of the emission region (based on the light crossing time of the fastest variability time scale), and  being the opening angle.The intrinsic opening angles of blazar jets have been found to lie between 0.1 • and 9.4 • , with the median being 1.2 • for sources detected by the Large Area Telescope onboard the Fermi satellite (Pushkarev et al. 2017), which we choose for Mrk 421.
The Doppler factors of blazars can be estimated via several different methods (see, e.g., Liodakis & Pavlidou 2015, for a review), but often these yield only lower limits, or even contradictory results for the same sources.Studies modelling broadband SEDs, in particular those that contain a data set representing a flaring state, usually require large values of  in order to explain the high luminosity at -ray energies.For example, to model the high-activity states of Mrk 421, Doppler factors in a range from 20 to 48 have been derived (Tavecchio et al. 1998;Donnarumma et al. 2009;Kapanadze et al. 2016;Banerjee et al. 2019;Zheng et al. 2021;Markowitz et al. 2022).However, other work relying on radio data and measuring the variability Doppler factor of Mrk 421 found  = 1.7 (Liodakis et al. 2017), and  = 2.03 (Liodakis et al. 2018), and a lower limit of  > 2.47 based on -ray opacity was determined by Pei et al. (2020).In Table 4, we therefore list the computed values for the time scales in the jet frame ( jet ), the size of the emission region ( r ), and the distance of the emission region to the jet base ( jb ) only in relation to .We derive these parameters for the fastest variability time scales in the full and the highest energy band, as well as based on the long-term variability defined by the duration of the X-ray flux increase.
In a similar way, we can estimate the magnetic field present in the emission region, both for the entire slowly-varying envelope as well as for the region where the rapid variability is created.Given that no lag is observed in the light curve, we know that cooling alone does not dominate the observed variability.Still, we can assume that a significant amount of variability is created through synchrotron cooling, for which the cooling time scale for an electron in the comoving frame is given as (e.g., Tavecchio et al. 1998) where   is the electron rest mass,  T is the Thomson cross section,  is the Lorentz factor of the ultra-relativistic electron, and  is the magnetic field in the co-moving frame of the jet.We can rearrange Eq. 8 to solve for the magnetic field .By applying  cool =  cool,obs /(1 + ), we obtain Given that we have no means to measure the Lorentz factor , we replace the Lorentz factor through its relation with the peak synchrotron emission  sync energy instead.Using the expression by Nalewajko et al. (2011), which is given for an averaged magnetic pitch angle, we can compute an estimate for the magnetic field  via We use  sync = 5 keV as the average energy of the 0.3-10 keV band.
For the large region that shows the slowly varying envelope and an increase time of roughly 57 ks in the full energy band (0.3-10 keV), we obtain a lower limit on the magnetic field of  ≥ 0.09 −1/3 G, which is similar to the findings of Markowitz et al. (2022).Using the shortest variability time for the 0.3-10 keV band, we estimate  ≥ 1.4 −1/3 G, and for 4-10 keV (using  sync = 7 keV for the average energy of that band)  ≥ 2.7 −1/3 G. Considering the range of values for  obtained through multiwavelength studies of the source, we can test if the peak luminosity of the high-energy component, which is created through synchrotronself Compton scattering, is at the expected value.For this, we make the simple assumption that the VHE -ray flux obtained with FACT is covering the decline of the SSC component close to the peak and use equation (2.4) from Ghisellini et al. (1996) to solve for a lower limit of the expected luminosity of the Inverse Compton component,  C .For the following calculations, we adopt a flat cosmology with  0 = 67.8km s −1 Mpc −1 , Ω  = 0.692, and Ω M = 0.308 (Planck Collaboration et al. 2016).Using a spectral index of  = 3.35 for the -ray spectrum, which is the average slope of FACT spectra reported by Arbet-Engels et al. (2021), and a flux of 5.1 × 10 −11 ph cm −2 s −1 , we obtain   ≥ 2.5 × 10 44 erg s −1 .To estimate the luminosity of the synchrotron peak, we use the spectral parameters during the highest X-ray flux measured with XMM-Newton (see Table 1), which yields L syn ∼ 6.9 × 10 45 erg s −1 .For a range of Doppler factors from 20 to 48, we find 1.6 × 10 44 erg s −1 ≤ L C ≤ 1.7 × 10 46 erg s −1 .Since the luminosity of the SSC component during our observation is at the lower end of that range, we adopt a Doppler factor of  = 40 for the continuing discussion regarding the origin of the X-ray emission in the jet.
If we use the shortest time scales to constrain the full size of the emission region responsible for the overall flaring activity, we find that it has to be located relatively close to the central engine, since  jb ≤ 280 R S and  jb ≤ 920 R S using  var measured for the 4-10 keV band and 0.3-10 keV band, respectively.However, in addition to the rapid variability embedded in the X-ray light curve, Mrk 421 exhibits a steady flux increase from  = 13 ks to  = 70 ks (see, e.g., Fig. 2).Using a time scale of ∼ 57 ± 1 ks for the light crossing time of the emission region in the observers frame, we conclude that the emission region responsible for the steady increase, if covering the entire width of the jet, would likely lie at a distance of ∼1 pc from the central engine.The size of this larger emission region is comparable to what has been found by Tramacere et al. (2022) for Mrk 421, based on a model involving an expanding adiabatic blob, and which is able to produce the frequently observed time lags between GeV -ray and radio emission from blazars.However, neither this model nor the standard shock acceleration model are able to explain the observed rapid variations, requiring us to turn to other, complementary models, such as a 'minĳets-in-jet'(e.g., Giannios 2013; Shukla & Mannheim 2020) or a 'spine-sheath' / 'current sheet' scenario, in which (kink) instabilities lead to magnetic reconnection events (e.g., Uzdensky et al. 2010;Sironi et al. 2015;Petropoulou et al. 2016;Bodo et al. 2021;Zhang et al. 2022).Both of these models were developed to explain rapid -ray variability on time scales of a few minutes, which was observed for a handful of sources during extremely bright flares (Aharonian et al. 2007;Albert et al. 2007b;Aleksić et al. 2014;Ackermann et al. 2016;Meyer et al. 2019); but similarly rapid variability is predicted for the X-ray regime as well (Christie et al. 2019, e.g.,).In the context of the 'mini-jets' scenario, the extremely rapid -ray variability events are distinctly observed as additional flares on top of a more slowly evolving (-ray) light curve.Similarly, we would expect a coherent short-term rise and decay as explicit structures in the light curve while a mini-jet crosses our line of sight and synchrotron emission appears enhanced due to additional Doppler boosting, but we do not observe such a signature.Therefore, we deem such a scenario unlikely for this observed flare in Mrk 421.Instead of distinct features, the variability time scale and PSDs capture fluctuations in the light curve, which could either originate from 1/   noise, or a multitude of small acceleration sites that are not resolvable in our X-ray light curve such as plasmoids in a reconnection layer in the jet that power flares when exiting the current sheet or 'spine' (Petropoulou et al. 2016).When comparing our observables to the predictions, we find that an emission region with a light-travelling distance equal to the rapid variability time scale in the 0.3-10 keV band has a size of roughly 5.7 × 10 14 cm, which falls into the range of sizes for small and fast plasmoids derived by Petropoulou et al. (2016).While we are not able to resolve individual flares produced by these small but highly relativistic plasmoids, we might see their signatures as variability on top of the slowly increasing envelope of the flare.Plasmoid-dominated reconnection is also able to explain the full X-ray flare that we observe over the duration of > 15 hours when considering large, or so-called 'monster' plasmoids (Christie et al. 2019).In addition, the merging of two plasmoids can be seen as an instant injection of particles that is biased towards high energies, which would result in, e.g., spectral hardening (Christie et al. 2019), which we clearly observe in the X-ray spectra presented in this work.We conclude that a scenario based on plasmoids undergoing magnetic reconnection in spine layer embedded in the jet is able to explain the observed X-ray flare in Mrk 421.
The observation presented in this paper reveals a strong correlation between different X-ray bands consistent with zero lag as shown in the zDCF results.The hysteresis curve shows a clockwise rotation, which implies the electron cooling being dominant over the electron acceleration.Indeed, this energy dependence, also visible in the higher fractional variability amplitude at higher energies, implies that the rapid variability likely originates from acceleration/cooling mechanisms instead of a change in the magnetic field or the Doppler factor.Recent X-ray polarization observations of Mrk 421 with the Imaging X-ray Polarimetry Explorer (IXPE) yield a detailed look into the jet region where X-rays are produced.During the first observation of the source with IXPE in May 2022, the polarization degree and angle stayed constant while the flux varied, which points towards shock acceleration being the dominant mechanism responsible for the X-ray emission (Di Gesu et al. 2022).However, two observations in June 2022 revealed a significant rotation of the X-ray polarization angle in combination with a harder-when-brighter behaviour, while no angle rotations were observed at optical or radio frequencies (Di Gesu et al. 2023).From their observations, Di Gesu et al. (2023) conclude that the X-ray emission site is a shock region moving along a helical magnetic field of the jet, similar to a spine-sheath scenario.This X-ray observation of Mrk 421 was taken during an enhanced state of VHE -ray activity.While the -ray data lacks the time resolution obtained with XMM-Newton, the -ray flux measured at the beginning of the XMM-Newton observation and again shortly after, reveals in an increase in amplitude by a factor of 2 (Gokus et al. 2022).In 2019, no instrument existed to take X-ray polarization measurements, hence we cannot know what the polarization properties during our XMM-Newton observation were.Rapid polarization swings can be associated with -ray flares (Zhang et al. 2022), however given the existing FACT data, we observe, if any, only a moderate flare at VHE -rays compared to past flaring activity.
While this XMM-Newton observation provides temporally very high-resolution X-ray data, no comparable data for the VHE -ray energy range exists for the duration of this flare, hence it is not clear if Mrk 421 exhibited variability on similar time scales at TeV energies as well.A simultaneous coverage of X-ray and TeV energies holds the potential of detecting contemporaneous rapid variability, such as predicted by plasmoid-dominated reconnection (Christie et al. 2019).Testing these particle acceleration models will be feasible via multi-wavelength campaigns with the Cherenkov Telescope Array and existing as well as future X-ray missions, such as the High-Energy X-ray Probe (HEX-P; Madsen et al. 2019), which is designed to provide very high resolution timing data across a broad X-ray energy range from 0.2 to ∼150 keV, and therefore would be ideally suited for high-energy timing studies of blazars (Marcotulli et al. 2023).

Figure 1 .
Figure1.VHE -ray (top) and X-ray (bottom) light curve covering the time span from 2019 June 2 to 2019 June 18.Data of Mrk 421 at -ray energies (E th ∼730 GeV) was obtained through the long-term monitoring campaign by FACT, while monitoring of the source in the X-ray band (0.3-10 keV) was done with Swift -XRT (dark blue), which is taken fromGokus et al. (2022).The XMM-Newton data (teal) is included as well and covers the majority of the X-ray flare.The color scale for the -ray light curve bins displays the detection significance for Mrk 421 for each observation.The inlay of the color scale does not cover any -ray data since there was no visibility of Mrk 421 after MJD 58649 for FACT.Note that the uncertainties for the X-ray flux are so well constrained that error bars are not visible.

Figure 2 .
Figure 2. 100s-binned light curve of Mrk 421 taken with XMM-Newton.The upper panel displays the full energy range (0.3-10 keV), while the subsequent panels show the varying flux in the soft (0.3-1 keV), medium (1-4 keV) and hard (4-10 keV) band, respectively.The data taken in the first 10 ks (light gray) are not taken further into account due to a high risk of being contaminated by strong background flaring.

Figure 3 .
Figure3.XMM-Newton light curve of Mrk 421 in 1 ks binning, displaying the count rate for the full energy range (0.3-10 keV).The rainbow color-coding provides an indicator of progressing time and has been applied to connect the hardness ratios in Fig.4to their respective occurrence in the time frame of the observation.

Figure 4 .
Figure 4. HIDs for medium (1-4 keV) vs. soft (0.3-1 keV) band (left) and the hard (4-10 keV) vs. soft band (right).The grey data points are the hardness ratios of the 100-s binned light curve and display the scatter of hardness ratios on very short time scales, while the colored data points are obtained from the 1 ks-binned time series.The colors match the time bins of the light curve displayed in Fig.3.The sections that were taken in Timing and Burst mode, respectively, are separated by the dashed gray line.

]Figure 6 .
Figure 6.Combined power spectral densities for Timing mode light curves with 10 s, 1 s, and 10 ms binning (top panel), the Poisson noise subtracted version of these PSDs (middle panel), and the Poisson-reduced PSDs multiplied by frequency (bottom panel).The gray line shows the upper 1 envelope for the uncertainty of the Poisson reduction.

Figure 7 .
Figure 7. zDCFs computed for the energy-resolved XMM-Newton light curves of Mrk 421.From top to bottom, the panels display the correlation between the soft and medium, soft and hard, and medium and hard energy band.The dashed line marks the zero lag, and the threshold for 0.8 ×  max is marked by the gray dotted horizontal line.Values that fulfil > 0.8 ×  max and are taken into account for computing  cent are highlighted by the red-shaded area in the background.
is divided by the average flux of the time span taken into account. 2 describes the general variance, and  2 err is the expected measurement uncertainty.The uncertainty of  var can be calculated via

Table 2 .
Fractional variability amplitude,  var , and minimum variability time scale,  var , computed for the full energy range as well as the three energy sub-bands.

Table 3 .
Results from applying the discrete correlation function to the 100sbinned light curves.Lags are defined such that positive values indicate harder photons leading softer photons.

Table 4 .
Transferred time scales in the jet frame ( jet ), size of the emission region ( r ), and distance of the emission region to the jet base ( jb ) under the assumption that the region covers the full width of the jet.We assume  BH = 1.9 × 10 8 M ⊙ ,  = 1.2 • , and  = 0.031.Using those parameters, 1 R S = 2.9 • 10 13 cm.