Extreme Hard X-Ray Variability in the TeV Blazar Mrk 421: Insights from NuSTAR Observations

1. Introduction

Blazars constitute the most energetic class of radio-loud active galactic nuclei (AGN) in the universe. These remarkable sources are characterized by relativistic jets aligned closely with the observer’s line of sight, producing intense non-thermal emission across the entire electromagnetic spectrum, from radio to high-energy $\gamma$-rays. The relativistic motion within these jets, aligned close to the line of sight, results in Doppler boosting that significantly enhances the observed emission and drives rapid variability across all wavebands [1,2]. Powered by supermassive black holes, AGN are among the universe’s most potent energy reservoirs. The energy extracted from these black holes drives relativistic jets, accelerating particles to near-light speeds. This results in Doppler-boosted emission with rapid variability across all temporal and spectral scales, making blazars ideal laboratories for studying extreme astrophysical processes and high-energy particle acceleration [3].

The broadband spectral energy distribution (SED) of blazars typically exhibits a characteristic double-peaked structure in the frequency-flux plane. The lower-frequency peak, spanning infra-red to X-ray bands, arises from synchrotron emission generated by relativistic particles accelerated within the jet’s magnetic field. The origin of the higher-frequency peak, extending from X-rays to TeV energies, remains a topic of active debate. Leptonic models attribute this emission to inverse-Compton scattering of soft seed photons by high-energy electrons [4,5]. In the synchrotron self-Compton (SSC) scenario, the same electron population responsible for synchrotron emission upscatters these photons to higher energies. Alternatively, the external Compton (EC) scenario proposes that seed photons originate from external AGN components, such as the accretion disk [6], broad-line region [7], or dusty torus [8]. In contrast, hadronic models suggest that high-energy emission results from interactions between relativistic protons and ambient radiation fields [9,10,11].

Blazars are classified into two main categories: flat-spectrum radio quasars (FSRQs) and BL Lacertae (BL Lac) objects. FSRQs are generally more luminous, displaying prominent optical emission lines and synchrotron peaks at lower frequencies. BL Lac objects, by contrast, are less luminous, with weak or absent emission lines and synchrotron peaks at higher frequencies. BL Lacs represent an extreme subclass, with synchrotron emission extending to soft X-rays (up to ${10}^{16}$ Hz) and inverse-Compton emission spanning hard X-rays to TeV energies. These sources typically lack strong circumstellar photon fields and are thought to accrete matter at relatively low rates [see [12], and references therein]. BL Lac objects are further subdivided based on their synchrotron peak frequency (${\nu }_{\mathrm{sy}}$): low-synchrotron-peaked (LSP, ${\nu }_{\mathrm{sy}}<{10}^{14}$ Hz), intermediate-synchrotron-peaked (ISP, ${10}^{14}$ Hz $<{\nu }_{\mathrm{sy}}<{10}^{15}$ Hz), and high-synchrotron-peaked (HSP, ${\nu }_{\mathrm{sy}}>{10}^{15}$ Hz) sources [13].

Blazars exhibit remarkable variability across all frequencies, from radio to TeV wavebands, spanning timescales from minutes to decades. Multi-wavelength (MWL) variability has emerged as one of the defining properties of blazars–and of AGNs in general–with variability studies playing a significant role in advancing our understanding of these sources. This variability can be characterized as a combination of three primary components: a general red noise-like behavior dominating the broadband power spectral density, occasional high-amplitude flares, and, in certain blazars, quasi-periodic oscillations with diverse characteristic timescales [14,15,16,17,18].

Markarian 421 (Mrk 421), a high-synchrotron-peaked BL Lac object at a redshift of $z=0.031$, is one of the closest and brightest TeV-emitting blazars, making it an ideal target for studying relativistic jet physics in AGN. Its non-thermal emission spans the entire electromagnetic spectrum, from radio to very-high-energy (VHE) $\gamma$-rays, exhibiting significant flux and spectral variability. Consequently, Mrk 421 has been the focus of extensive multi-wavelength (MWL) observational campaigns, covering radio to TeV energies, which have revealed complex variability patterns throughout its observational history [19,20,21,22]. These broadband efforts have produced invaluable datasets, providing crucial constraints for theoretical models of jet emission mechanisms [23,24,25]. Studies utilizing (quasi-)simultaneous X-ray and $\gamma$-ray observations have frequently reported correlated variability between these bands [26,27,28,29,30,31,32,33]. In the X-ray regime, Mrk 421’s variability has been extensively studied using instruments such as XMM-Newton, Swift, NuSTAR, and AstroSAT. These analyses reveal strong correlations between soft and hard X-ray emission, with spectra often described by power-law or curved models [34,35,36,37,38]. Observations with Chandra have also provided evidence for warm-hot intergalactic medium (WHIM) filaments along the line of sight, confirming their role as reservoirs of baryonic matter in the cosmic web [39,40]. Detailed studies of flaring episodes have revealed interesting spectral and temporal features, including spectral hysteresis loops in the flux–hardness ratio plane and evidence of both hard and soft lags [30,41,42]. More recently, observations with the Imaging X-ray Polarimetry Explorer (IXPE) detected an X-ray linear polarization degree of ∼15% during a moderately active flux state [43].

In this paper, we present a comprehensive analysis of all available NuSTAR observations of Mrk 421, focusing on timing and spectral analyses to explore the variability properties of this blazar in the hard X-ray regime. Our study aims to explore the physical processes driving blazar emission by exploring the relationships among variability, mean flux, and spectral properties, revealing insights into the innermost regions of blazar jets, which remain mostly inaccessible to direct observation.

We organize the paper as follows: In Section 2, we describe the 28 NuSTAR observations and the data processing methods. In Section 3, we present the methods used for timing and spectral analyses based on light curves and spectra, and we discuss the results. In Section 4, we discuss key observational findings, including rapid flux and spectral variability, and a correlation between higher flux and harder spectra interpreting these in the context of current blazar models. Finally, in Section 5, we summarize our conclusions.

2. NuSTAR Observations and Data Processing

The Nuclear Spectroscopic Telescope Array (NuSTAR) is a sensitive hard X-ray (3–79 keV) observatory equipped with two focal plane modules, FPMA and FPMB, designed to provide high-resolution imaging and spectroscopy.

The NuSTAR instrument achieves a spectral resolution of ∼400 eV FWHM at 6 keV, increasing to ∼1 keV at 68 keV across its 3–79 keV bandpass. Each telescope provides a field of view of ∼13’ and a half-power diameter of ∼58" for a point source image[44]. This capability makes NuSTAR an ideal tool for studying the hard X-ray variability and spectral properties of blazars like Mrk 421, which often exhibit complex emission in this energy range.

Data reduction was performed using the NuSTAR Data Analysis Software (NuSTARDAS) package version 1.3.1, integrated with HEASOFT1 version 6.33 and the calibration database CALDB version 20220331. All observations were processed with the standard nupipeline script to produce calibrated and cleaned event files. Source emission was extracted from a circular region of 30" radius centered on the source position, while background emission was sampled from a 60" radius circular region in a source-free area of the detector, ensuring minimal contamination while remaining on the same detector chip. Light curves were generated with a uniform time binning of 300 seconds, balancing temporal resolution with a sufficient signal-to-noise ratio for variability analysis. Spectra were extracted and re-binned using the grppha task to ensure a minimum of 20 counts per energy channel, enabling the application of ${\chi }^2$ statistics in spectral fitting. This rigorous data processing approach ensures reliable characterization of the source’s hard X-ray behavior, facilitating detailed comparisons with theoretical emission models.

3. Analysis and Results

To characterize the hard X-ray emission from the blazar Mrk 421, we conducted a comprehensive analysis of both flux variability and spectral properties using 28 NuSTAR observations of the source, selected based on two primary requirements: an exposure time greater than 1000 seconds and an issue flag of zero. The observation date, observation ID, and exposure time for each selected dataset are presented in the 2nd, 3rd, and 4th columns of Table 1, respectively, providing a clear overview of the observational parameters used in this study.

Our investigation can be divided into two main components: timing analysis and spectral analysis, which are presented in the subsequent sections. The timing analysis explores the temporal behavior of the source across multiple timescales, while the spectral analysis examines the underlying emission mechanisms and their evolution. Together, these complementary approaches provide insights into the physical processes driving the hard X-ray emission from this well-studied blazar.

3.1. Timing Analysis

To carry out timing analysis, we processed the data to generate intra-day light curves, as detailed in Section 2, and present them in Figure 1 in chronological order, corresponding to the observations listed in Table 1. The upper panel of Figure 1, displayed on a linear y-scale, reveals a dramatic flux change over the entire period, while the lower panel, using a logarithmic y-scale, emphasizes the intra-day variability, particularly in the low flux states. Building on these light curves, we systematically analyzed the flux variability properties across both intra-day and long-term timescales. The intra-day analysis focused on characterizing rapid fluctuations within individual observing epochs, while the long-term investigation examined variability patterns spanning the entire observational campaign. This dual-timescale approach allows us to probe both the rapid, local processes occurring within the jet and the more sustained evolutionary trends in the emission properties of Mrk 421.

3.1.1. Flux Variability

As the flux state of the source varied across the observation epochs, we computed the mean flux for each of the 28 intra-day hard X-ray light curves, as presented in the fifth column of Table 1. Across the observations, the mean flux exhibits a wide range, from a low of 1.76 counts/s on 2013-01-10 to a high of 133.25 counts/s on 2013-04-14, with an average of 24.28 counts/s over the dataset. This broad span highlights the source’s dynamic behavior in the hard X-ray band, capturing both its baseline and peak activity levels. Notably, the maximum mean flux of 133.25 counts/s on 2013-04-14 indicates a significant flaring event, one of the most intense recorded for this source in this energy regime [45]. The observed flux variations provide valuable insights into the source’s energetic processes over time.

While a light curve’s variance provides an initial measure of variability, it includes both intrinsic source fluctuations and measurement uncertainties. To isolate the source’s true physical variability, we must subtract the variance contribution from measurement uncertainties in the count rates. This corrected metric, known as excess variance, effectively quantifies the intrinsic variability by eliminating statistical noise from the measurement process. For a light curve comprising N measured flux values $X_i$, each with associated finite uncertainties ${\sigma }_{\mathrm{err},\mathrm{i}}$ due to measurement errors, and given $S^2$ as the sample variance of the light curve, the excess variance is computed using the following relation:

$${\sigma }_{XS}^2=S^2-{\overline{\sigma }}_{\mathrm{err}}^2,$$

(1)

where ${\overline{\sigma }}_{\mathrm{err}}^2$ represents square of the mean measurement errors, providing a robust measure of the source’s true variability [46,47]. Fractional variability, which is related to the normalized excess variance, is calculated using the relation,

$$F_{\mathrm{var}}=\sqrt{{\displaystyle \frac{S^2-{\overline{\sigma }}_{\mathrm{err}}^2}{{\overline{X}}^2}}}$$

(2)

[see 48,49]. The associated error in the fractional variability is obtained using,

$${\sigma }_{F_{\mathrm{var}}}=\sqrt{{\left({\displaystyle \frac{1}{\sqrt{2N}}}{\displaystyle \frac{{\overline{\sigma }}_{err}^2}{F_{var}}}{\displaystyle \frac{1}{{\overline{X}}^2}}\right)}^2+{\left(\sqrt{{\displaystyle \frac{{\overline{\sigma }}_{err}^2}{N}}}{\displaystyle \frac{1}{{\overline{X}}^2}}\right)}^2}$$

(3)

[50,51]. We computed the fractional variability (FV) for all 28 intra-day hard X-ray light curves using the relations described above, with results listed in the 6th column of Table 1. The analysis reveals a remarkable range of variability across observations, spanning from a minimal value of approximately 1% (observed on 2024-04-29) to an extraordinary peak of 60% (recorded on 2013-04-11). On average, the intra-day observations exhibit a fractional variability of approximately 18% across the entire dataset. These values highlight the dramatic variability characteristics of this source in the hard X-ray band, demonstrating both extreme fluctuations during its most active states and the wide dynamic range of variability the source experiences over time. The exceptional FV value of 60% during the 2013-04-11 observation represents one of the most violent variability episodes reported for this class of object in the hard X-ray regime. Blazars, including Mrk 421, are well-known for exhibiting rapid flux variations on intra-day timescales across a broad range of electromagnetic wavebands [e.g., [12,37,52,53,54,55]. These short-timescale fluctuations, often detected in the hard X-ray regime, provide a powerful probe of the physical processes within the compact emission regions of blazar jets. Using light-travel-time arguments, the intra-day variability suggests that the emission originates from highly compact regions, likely located close to the central engine. By analyzing the variability in the hard X-ray emission of Mrk 421, we can constrain some of the important properties of these regions, including their size and relativistic kinematics, offering valuable insights into the dynamics and structure of the jet.

3.2. RMS–Flux Correlation

The root-mean-square (RMS)–flux relation in AGN reveals a linear correlation between the RMS variability amplitude and the mean flux, observed across a wide range of timescales [56,57,58]. Initially identified in the X-ray emission from AGN and X-ray binaries, this relation indicates that short-term flux variations scale multiplicatively with longer-term brightness changes, a characteristic of accretion-driven processes. The remarkable linearity persists across diverse timescales and energy bands, suggesting that variability arises from propagating fluctuations in the accretion flow, where turbulent variations in the mass accretion rate at different radii are coupled [56].

Recently, the rms–flux relation has been investigated in several blazars across different observational wavebands, with the aim of probing the disk–jet connection in jetted AGNs. The relation has been found to hold in the gamma-ray band (e.g., [59,60,61]), X-rays [62], and optical wavelengths [63]. However, in long-term variability studies of the blazars 3C 273 and PKS 2155-304 based on XMM-Newton X-ray observations, the relation was not observed [54,55].

Motivated by these findings, we investigated whether the variability amplitude in the hard X-ray emission from Mrk 421 correlates with its mean flux state, aiming to further understand the nature of its variability and the physical processes driving the source’s behavior. We quantified the variability using the root-mean-square (RMS) amplitude, defined as $\mathrm{RMS}=\sqrt{\mathrm{excess}\mathrm{variance}}$, and computed the Pearson’s correlation coefficient between the RMS and the mean count rate across all 28 NuSTAR light curves. The analysis yielded a correlation coefficient of 0.86, indicating a strong positive correlation between the two quantities. This suggests that higher flux states are associated with greater variability amplitudes, a trend commonly observed in blazars and often linked to non-linear, multiplicative processes in the jet [56]. The relationship is visually represented in Figure 2, where a plot of the mean count rate against the RMS clearly illustrates the linear trend, reinforcing the statistical correlation. The relation provides evidence for a connection between the source’s brightness and its variability, potentially reflecting dynamic changes in particle acceleration or energy injection within the jet of Mrk 421.

3.2.1. Flux Distribution and PDF

Understanding the statistical distribution of X-ray flux from blazars like Mrk 421 provides critical insights into the physical processes governing their emission. In particular, the distinction between normal (Gaussian) and lognormal probability density functions (PDFs) serves as a key discriminator between theoretical models of particle acceleration and emission mechanisms. Normal distributions typically arise from additive processes, while lognormal distributions are associated with multiplicative processes, such as turbulent or cascading energy injection in the jet [56]. This analysis is particularly valuable for constraining the underlying physics of blazar variability.

The normal distribution is expressed as:

$$N\left(x\right)={\displaystyle \frac{1}{\sqrt{2\pi {\sigma }^2}}}exp\left(-{\displaystyle \frac{{(x-\mu )}^2}{2{\sigma }^2}}\right),$$

(4)

where $N\left(x\right)$ is the probability density function, x represents the flux value, $\mu$ is the mean, and $\sigma$ is the standard deviation, quantifying the distribution’s spread.

The lognormal distribution is given by:

$$Ln\left(x\right)={\displaystyle \frac{1}{\sqrt{2\pi }sx}}exp\left(-{\displaystyle \frac{{(lnx-m)}^2}{2s^2}}\right),$$

(5)

where x is the positive flux value, m is the mean of the natural logarithm of x, and s is the standard deviation of the logarithm. The $1/x$ factor ensures proper normalization on the logarithmic scale. Unlike the symmetric normal distribution, the lognormal distribution is skewed, with a longer tail toward higher values, a feature often observed in blazar flux distributions.

We combined the all 28 intra-day light curves to construct a histogram of count rate distribution, reflecting the source’s variability across multiple observational epochs. Both normal and lognormal PDFs were fitted to the histogram using a least-squares method, with goodness-of-fit evaluated via the reduced ${\chi }^2$ statistic.

Figure 3 illustrates the histogram of flux variations, with the normal distribution (magenta line) and lognormal distribution (green line) overlaid. The lognormal PDF provides a significantly better fit, yielding a reduced ${\chi }^2$ of $\sim 7$, compared to $\sim 35$ for the normal distribution. An inset text box in the figure reports the ${\chi }^2$ values and degrees of freedom for both fits, allowing a direct comparison. This result indicates that the hard X-ray flux distribution of Mrk 421 is more consistent with lognormal behavior, suggesting multiplicative processes dominate the variability.

The preference for a lognormal distribution has been reported in previous studies using Swift/XRT observations of Mrk 421 [64], and similar lognormal behavior has also been observed in XMM-Newton observations of several other blazars [65]. In contrast, a study based on archival XMM-Newton data spanning a decade found no clear evidence of a distinct probability distribution function (PDF) in the flux distributions of the FSRQ 3C 273 and the BL Lac object PKS 2155–304 [54,55]. Complementing the X-ray results, flux distributions of several blazars—including Mrk 421—have also been shown to be consistent with a lognormal PDF in the gamma-ray [59,60,61] and optical bands [63].

In blazars, the presence of a lognormal PDF in their flux distributions is often attributed to non-linear processes, potentially driven by turbulent energy injection or cascading mechanisms within the relativistic jet [66]. This contrasts with additive processes, often associated with shot-noise-like behavior and characterized by a normal distribution [59]. These findings enhance our understanding of the physical mechanisms driving Mrk 421’s variability and contribute to the broader study of blazar jet dynamics.

3.3. Spectral Analysis

In the following section, we present the spectral analysis of NuSTAR observations of the blazar Mrk 421. This analysis employs a model-independent hardness ratio approach to characterize spectral variability and detailed spectral fitting using various models to describe the NuSTAR spectra.

3.3.1. Hardness Ratio

To investigate the relationship between the variability in different energy bands of Mrk 421, we extracted light curves in two distinct energy ranges—low-energy (LE, 3–10 keV) and high-energy (HE, 10–79 keV)—following the standard data reduction procedures outlined in Section 2. Our objective was to examine how the high-energy light curve responds to variations in the low-energy light curve, thereby probing the connection between these emission components across the entire observation period. To quantify this relationship, we computed the Pearson’s correlation coefficient between the LE and HE light curves, yielding a value of 0.94, which indicates a strong positive correlation between the two energy bands. Several studies, utilizing various X-ray telescopes such as XMM-Newton, Swift/XRT, Chandra, and NuSTAR, have reported a similar positive correlation between the emission in the two bands [67,68,69,71,71].

This high degree of correlation suggests that the variability in the LE and HE bands is closely linked, likely driven by common physical processes within the jet. The correlation is visually represented in Figure 4, where we plot the HE light curve against the LE light curve, with each point corresponding to simultaneous flux measurements in the two bands. Additionally, we performed a linear fit to the data, finding that a slope of 0.6 provides the best fit, indicating a proportional but slightly subdued response of the HE emission to changes in the LE emission. For comparison, a dashed red line representing the $x=y$ relation (slope of 1.0) is also included in the figure, highlighting that the HE variability is less sensitive to LE changes. This slope of 0.6 suggests that while the variability in both bands is strongly correlated, the source’s overall variability is dominated by the soft (LE) band, potentially reflecting differences in the electron populations or cooling mechanisms responsible for the emission in these energy ranges [70].

To further evaluate the relative contributions of different spectral components to the total flux, as measured by count rates, and to examine how these contributions vary across diverse flux states, we define the hardness ratio (HR) as a model-independent metric of the blazar’s spectral shape. The HR is calculated as:

$$\mathrm{HR}={\displaystyle \frac{H}{S}},$$

(6)

where H represents the flux in the hard X-ray energy band (2–10 keV) and S denotes the flux in the soft X-ray energy band (0.3–2 keV). This definition provides a direct measure of the contribution from the hard-band over soft band, facilitating a standardized comparison across observations. The uncertainty in the hardness ratio, ${\sigma }_{\mathrm{HR}}$, is derived using error propagation, given by:

$${\sigma }_{\mathrm{HR}}=\mathrm{HR}·\sqrt{{\left({\displaystyle \frac{{\sigma }_{F_{\mathrm{hard}}}}{F_{\mathrm{hard}}}}\right)}^2+{\left({\displaystyle \frac{{\sigma }_{F_{\mathrm{soft}}}}{F_{\mathrm{soft}}}}\right)}^2},$$

(7)

where ${\sigma }_{F_{\mathrm{hard}}}$ and ${\sigma }_{F_{\mathrm{soft}}}$ are the uncertainties in the hard and soft band fluxes, respectively.

The hardness ratio provides a straightforward, model-independent representation of the spectral properties of Mrk 421 within the specified energy range, with variations in HR directly reflecting spectral changes during the observations. To investigate potential correlations or anti-correlations between the spectral properties and the flux state of Mrk 421, we examined the relationship between the hardness ratio (HR) and the total flux (measured in counts) in the entire band. The right panel of Figure 4 displays the HR values for the entire dataset plotted against the logarithm of the flux in counts, revealing a weak positive correlation characterized by a small slope of 0.05. This shallow trend indicates a nearly negligible dependence between HR and flux, suggesting that the spectral shape remains largely stable despite variations in the source’s brightness. This absence of a correlation between HR and flux suggests that the hard X-ray emission from Mrk 421 likely originates from multiple emission zones, each potentially governed by distinct emission mechanisms, such as synchrotron and inverse Compton processes, or a combination of leptonic and hadronic contributions [71]. These findings highlight the complexity of the source’s spectral behavior and the need for multi-zone models to fully capture the underlying physical processes driving its variability.

3.4. Spectral Modeling

To model the non-thermal emission from the blazar Mrk 421, we conducted a detailed spectral analysis of 28 NuSTAR observations using the XSPEC software package, following standard X-ray spectral analysis procedures [72]. During data processing, we generated source spectra, response matrix files (RMF), ancillary response files (ARF), and background spectra for each observation, binning the spectra to ensure a minimum of 20 counts per energy bin for reliable ${\chi }^2$ statistics. The analysis was performed in the 3–79 keV energy range, simultaneously fitting the spectra from the FPMA and FPMB instruments while incorporating an intercalibration constant to account for minor instrumental differences. The intercalibration constants ranged from 0.97 to 1.04, indicating excellent consistency between the two detectors.

We tested three distinct spectral models to determine the best representation of Mrk 421’s spectral behavior: a power law (PL), a log-parabola (LP), and a broken power law (BPL). The power-law model is defined as:

$${\displaystyle \frac{dN}{dE}}=N{E}^{-\mathrm{\Gamma }},$$

(8)

where N, E, and $\mathrm{\Gamma }$ represent the normalization, photon energy, and photon index, respectively. This model reflects the non-thermal emission from a power-law distribution of high-energy particles in the jet’s magnetic field, with the synchrotron spectral index $\alpha$ related to the particle distribution index p by $\alpha =(p-1)/2$.

The broken power-law model, which allows for a spectral break, is expressed as:

$${\displaystyle \frac{dN}{dE}}=N_0\left\{\begin{array}{cc}{E}^{-{\mathrm{\Gamma }}_1},\hfill & \mathrm{if}E\ge E_b\hfill \\ E^{-{\mathrm{\Gamma }}_2},\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(9)

where ${\mathrm{\Gamma }}_1$ and ${\mathrm{\Gamma }}_2$ are the photon indices below and above the break energy $E_b$, and $N_0$ is the normalization. A broken power-law spectrum often indicates the presence of high-energy relativistic electrons that cool rapidly through radiative processes, leading to a break in the emission spectrum.

The log-parabola model, characterized by a continuous curvature, is given by:

$${\displaystyle \frac{dN}{dE}}=N_0{\left({\displaystyle \frac{E}{E_0}}\right)}^{-\left(\alpha +\beta log\left(E/E_0\right)\right)},$$

(10)

where $N_0$ and $E_0$ (fixed at 1 keV) are the normalization and reference energy, respectively, while $\alpha$ and $\beta$ denote the photon index at $E_0$ and the curvature parameter, respectively [76]. The log-parabola model is often associated with particle acceleration scenarios where the probability of further energy gain decreases with increasing particle energy, resulting in a curved spectrum.

To account for Galactic absorption, we incorporated the tbabs model (Tuebingen-Boulder ISM absorption model; [73]) into each fit, using hydrogen column density values from [74]. For Mrk 421, we fixed the neutral hydrogen column density at $N_{\mathrm{H}}=1.34\times {10}^{20}{\mathrm{cm}}^{-2}$, consistent with previous studies [75]. These spectral models and absorption corrections enable a comprehensive analysis of Mrk 421’s emission properties, providing insights into the physical mechanisms driving its non-thermal radiation in the hard X-ray regime. The results of the spectral fitting are presented in Table 2. The table includes the observation ID and the XSPEC model in the first and second columns, respectively, while the best fit parameters for the PL, BPL, and LP models are provided in the third, fourth, and fifth columns. The final column reports the reduced ${\chi }^2$ statistic for each fit, enabling a direct comparison of the model performance. The table indicates that the LP model provides the best fit for 20 of the 27 observations of Mrk 421, suggesting that the spectra frequently exhibit curvature. This is consistent with particle acceleration scenarios where the probability of energy gain decreases with increasing energy. The BPL model best fits 2 observations, indicating occasional spectral breaks, possibly due to rapid cooling of high-energy electrons. The PL model is the best fit for 5 observations, showing that a simple power-law spectrum adequately describes Mrk 421’s hard X-ray emission in these cases. For one observation (60002023008), fitting convergence could not be achieved, preventing any model from being successfully applied. The preference for LP and BPL models supports the expectation of complex, non-thermal emission from blazar jets [76].

Spectral modeling results for two distinct observation epochs of Mrk 421 are presented in Figure 5, with the left panel corresponding to the lowest flux state and the right panel representing the highest flux state. For both observations, the log-parabola (LP) model provided the best fit, as evidenced by the reduced ${\chi }^2$ values. Notably, the lower panel of the low-state plot reveals a distinct hard excess beginning around 20 keV, indicating a deviation from the LP model at higher energies during this epoch see (also [33]).

To further explore the spectral evolution of Mrk 421 across different flux states, we analyzed the best-fit parameters from Table 2, focusing on their distribution and correlations. The top left panel of Figure 6 presents a plot of the photon index from PL model against the logarithm of the mean count rate for each observation (see Table 1), revealing a clear anti-correlation trend that indicates a slow “harder-when-brighter” behavior; as the flux increases (on a logarithmic scale), the photon index decreases, suggesting a gradual hardening of the spectrum. The top right panel of the same figure displays an anti-correlation between the $\alpha$ and $\beta$ parameters of the LP model, implying that as the spectrum steepens (higher $\alpha$), the curvature parameter $\beta$ decreases, resulting in a flatter spectral shape, with negative $\beta$ values particularly highlighting a hard excess in the two lowest flux states of the source. Similarly, the bottom left panel illustrates a correlation of mean flux level from each observation with the synchrotron peak energy $E_p$, derived using the relation $E_p=E_1\times {10}^{(2-\alpha )/2\beta }$ with $E_1=1\mathrm{keV}$. A positive correlation is observed, with a Pearson’s correlation coefficient of 0.86 but with a shallow slope of 0.02. This subtle shift of $E_p$ toward higher energies as the source brightens further supports the observed spectral hardening during higher flux states. In a similar previous study, log-parabolic spectral fit to the X-ray observations during the flaring state of the source revealed a positive correlation between spectral parameters and flux levels [30,77]. Furthermore, a correlation study examining the relationship between the peak luminosity and peak energy ($L_{\mathrm{p}}$–$E_{\mathrm{p}}$ relation) in a sample of 14 BL Lac objects reveals positive, negative, or no correlations between these parameters [78]. Finally, the bottom right panel of the figure shows no significant correlation between $E_p$ and $\beta$ parameter.

Similarly, to combine the results of timing and spectral analyses, the FV for each observation is plotted against the photon indices derived from PL spectral fitting in Figure 7. The plot reveals a trend of anti-correlation between these variables, suggesting greater variability and a broader spread during the harder spectral states of the source.

4. Discussion

The study of intra-day and long-term hard X-ray variability in the TeV blazar Mrk 421, based on 12 years of NuSTAR observations, reveals flux variations of up to a factor of 22 5. The FV in the X-ray band reaches ∼120%, significantly exceeding the 43% FV in gamma-rays [81] and 30% in the optical band [63] over comparable timescales. The source exhibited its highest flux state on April 14, 2013, and its lowest flux state on January 10, 2013. In addition, we observed a profound flaring event in April 2013, during which the source brightened by a factor of 10 within a period of 15 days. This provides a unique opportunity to probe the physical processes driving blazar variability. By performing comprehensive spectral and timing analyses—including the linear RMS-flux relation, log-normal flux distributions, and the power-law spectral models—we gain deeper insight into the mechanisms responsible for both intra-day and long-term flux and spectral variability, including this remarkable flaring event. These analyses collectively point to multiplicative, non-linear processes within the relativistic jet, likely driven by enhanced particle injection at a shock wavefront, with possible contributions from magnetic reconnection or accretion disk fluctuations.

The hard X-ray emission during the observation period exhibits distinct spectral characteristics. The nearly constant hardness ratio (HR = $F_{\mathrm{HE}}/F_{\mathrm{LE}}$, with a slope of 0.05 in the $logF$ vs. HR relation) and a strong linear correlation between high-energy (HE, 10–79 keV) and low-energy (LE, 3–10 keV) fluxes ($F_{\mathrm{HE}}\approx 0.6F_{\mathrm{LE}}$) suggest that the emission is dominated by leptonic synchrotron radiation from a single electron population with a stable power-law index and a minimal spectral evolution, consistent with self-regulated acceleration balancing radiative cooling ($t_{\mathrm{cool}}\propto {\left(B^2{\gamma }_e\right)}^{-1}$).

The log-parabolic spectral model, expressed by Equation 10, reveals a significant anti-correlation between the photon index $\alpha$ (2.38–3.64) and curvature parameter $\beta$ (-0.33–0.30), in contrast the positive correlation observed in other X-ray observations of the source [e. g., [68,76]. This discrepancy suggests that during the period of hard X-ray emission, particle acceleration is influenced by additional processes, such as stochastic acceleration or energy-dependent cooling, which alter the electron energy distribution and spectral shape. The synchrotron peak energy, calculated as $E_p=E_1\times {10}^{(2-\alpha )/\left(2\beta \right)}$ with $E_1=1\mathrm{keV}$[see [79], shifts from $\approx 0.23\mathrm{keV}$ ($\alpha =2.38$, $\beta =0.30$) to $\approx 302\mathrm{keV}$ ($\alpha =3.64$, $\beta =-0.33$) as the flux increases, aligning with a "harder-when-brighter" trend (photon index $\mathrm{\Gamma }$ decreases from 3.1 to 2.3 as mean count rate rises from 1.76 to 133.25 counts/s).

This spectral evolution is likely driven by an increase in the maximum electron energy (${\gamma }_{e,\max }$) and a modest change in the magnetic field (B). For $E_p\propto B{\gamma }_{e,\max }^2$ to increase by a factor of $(3/0.23)\approx 13$ (ignoring two lowest states), if B remains constant, ${\gamma }_{e,\max }$ must increase by $\approx 3.6$. In particular, to explain the 10-fold flux increase ($F\propto N_e{\gamma }_{e,\max }^2$), the electron density $N_e$ would need to adjust by $\approx 0.77$, suggesting stability or a slight decrease in $N_e$. Alternatively, if B increases by a factor of 2 (e.g., via magnetic reconnection), ${\gamma }_{e,\max }$ need only increase by $\approx 2.55$, with $N_e\approx 0.77$, better aligning with the observed flux and $E_p$ shift. The stable HR and linear flux correlation disfavor significant synchrotron self-Compton (SSC) or hadronic contributions, as SSC would predict a non-linear flux relationship ($F_{\mathrm{HE}}\propto F_{\mathrm{LE}}^2$; in Thomson’s region [80]) and hadronic emission would increase HR with flux (e. g., [81,82]).

A change in the Doppler factor (${\delta }_D$) due to jet reorientation is less likely, as $F\propto {\delta }_D^{3+\alpha }$ (for $\alpha \approx 1.3-2.1$) requires ${\delta }_D$ to increase by $\approx 1.7$, yielding an $E_p\propto {\delta }_D$ increase of only 1.7, far below the observed factor of 13. Thus, enhanced particle injection at a shock wavefront, possibly coupled with a slight increase in B, emerges as the primary driver of the flare, with leptonic synchrotron processes dominating the emission.

The timing properties of Mrk 421’s flare further shed light on the underlying variability mechanisms. The linear trend in RMS-flux relation, where the root-mean-squared (RMS) amplitude of variability generally scales linearly with the mean flux, indicates that higher flux states exhibit greater variability, suggesting a common mechanism linking short- and long-timescale variations. For Mrk 421, this relation in X-ray light curves from AstroSat and XMM-Newton across epochs spanning decades and shorter periods (5 days to ∼10 years), similar to patterns in X-ray binaries and Seyfert galaxies [56,62,83]. A broader study involving a large population of blazars, using 14 years of Fermi-LAT gamma-ray light curves, found that a significant number of sources—particularly BL Lac objects like Mrk 421—exhibited a strong linear RMS-flux relation. This relation was interpreted in terms of internal shock models with variations in shock amplitudes [59]. However, the generality of the relation is debated, as some simulations suggest it appears in only 26–32% of light curves under specific shock scenarios [84].

Complementing the RMS–flux relation, the analysis of NuSTAR light curves reveals a log-normal probability density function (PDF), with a reduced ${\chi }^2\sim 7$, significantly better than that of a normal distribution (${\chi }^2\sim 35$). In the AGN literature, these characteristics are considered indicative of multiplicative, non-linear processes. This log-normality, also reported in AstroSat and XMM-Newton X-ray data for Mrk 421 [62] and in Fermi-LAT gamma-ray data for 46.5% of BL Lacs [59], suggests turbulent or cascading energy injection in the jet [66]. In contrast, flat-spectrum radio quasars (FSRQs) like 3C 273 occasionally show normal distributions, highlighting diverse variability mechanisms [62].

The timing analyses, together with spectral analyses, present a coherent picture of the flare in Mrk 421. The linear RMS-flux relation and log-normal flux distribution point to multiplicative, non-linear processes, likely driven by internal shocks or turbulent energy injection within the jet. These processes align with the spectral findings, where the anti-correlation between $\alpha$ and $\beta$, coupled with the "harder-when-brighter" trend and $E_p$ shift, suggests complex particle acceleration mechanisms, possibly involving stochastic acceleration or energy-dependent cooling. The stable HR and linear flux correlation underscore the dominance of leptonic synchrotron radiation, with internal shocks increasing ${\gamma }_{e,\max }$ and $N_e$, potentially amplified by a modest increase in B.

This variability may originate from accretion rate fluctuations in the disk, which propagate along the jet and carry imprints of disk properties, constraining the disk-jet connection [60,63,85]. The anti-correlation between $\alpha$ and $\beta$, deviating from standard models, further suggests that additional factors, such as magnetic reconnection or spatially varying acceleration regions, modulate the electron energy distribution during the flare. While deviations from log-normality in FSRQs and the non-universal RMS-flux relation indicate complexity, the consistent log-normal PDFs and RMS-flux relation in Mrk 421 reinforce a turbulent jet origin for its variability.

The spectral and timing analyses of the Mrk 421 flare reveal a dynamic interplay of multiplicative, non-linear processes driving its variability. The log-parabolic spectral model, with an anti-correlation between $\alpha$ and $\beta$, points to complex particle acceleration mechanisms, while the "harder-when-brighter" trend and $E_p$ shift suggest enhanced particle injection at a shock wavefront, possibly coupled with magnetic reconnection. The linear RMS-flux relation and log-normal flux distribution support these findings, indicating turbulent energy injection and potential disk-jet coupling. These results highlight the dominance of leptonic synchrotron processes in Mrk 421’s jet and underscore the need for further multiwavelength studies to refine our understanding of blazar jet dynamics and their connections to the central engine.

5. Conclusions

We conducted a comprehensive timing and spectral analysis of the TeV blazar Mrk 421, utilizing 28 observations spanning 12 years from 2012 to 2024, to investigate its hard X-ray variability properties. Our key findings are summarized as follows:

• Mrk 421 exhibited extreme variability, with flux changes of up to a factor of 200 over a 12-year period, including a 15-day flaring event during which the source brightened by more than a factor of 10. This dramatic variability is likely driven by an increase in the maximum Lorentz factor during shock acceleration, highlighting the dynamic nature of the particle acceleration processes in the jet.
• The hard X-ray variability displays a linear RMS-flux relation and a lognormal flux distribution. These characteristics suggest that the variability arises from multiplicative processes, rather than additive ones, pointing to a cascading or interconnected mechanism governing the emission.
• A strong correlation is observed between the count rates in the low-energy (3–10 keV) and high-energy (10–79 keV) bands, indicating that synchrotron emission dominates the X-ray output across these energy ranges.
• Spectral analysis reveals a weak but consistent harder-when-brighter trend, with spectra best described by log-parabolic models exhibiting curvature and occasional breaks. The photon index and curvature parameters of the log-parabola model are anti-correlated, and the spectral peak shifts to higher energies with increasing flux. These features collectively support a scenario of energy-dependent stochastic particle acceleration shaping the emission.

These results provide valuable insights into the physical processes driving Mrk 421’s variability, underscoring the interplay of particle acceleration and synchrotron emission in blazar jets. Our findings pave the way for further studies to refine models of relativistic jet dynamics and their observational signatures.