NuSTAR View of TeV Blazar Mrk 501

We report the results of flux and spectral variability studies of all seven {\it Nuclear Spectroscopic Telescope Array (NuSTAR)} observations of TeV $\gamma-$ray emitting blazar Markarian (or Mrk) 501. We found strong evidence of intraday variability in 3-79 keV X-ray light curves (LCs) of Mrk 501 during four out of these seven observations. We examined spectral variability using a model-independent hardness-ratio analysis and found a general"harder-when-brighter"behaviour in two observations. We also investigated the nature of 3-79 keV X-ray spectra of TeV blazar Mrk 501 and found that five out of seven spectra are well described by the curved log-parabola models with photon indices (at 10 keV) $\alpha \sim$ 2.12-2.32 and a curvature $\beta \sim$ 0.15-0.28. The two other spectra are somewhat better represented by simple power-law models with photon indices 2.70 and 2.75. We briefly discuss available physical models to explain our results.


Introduction
Active galactic nuclei (AGN) are highly energetic astrophysical objects that are powered by actively accreting supermassive black holes (SMBH; M BH ≥ 10 6 M sun ) at their centres, and produce unique observational features over the entire electromagnetic (EM) spectrum [1]. About 15% of AGN consist of two well-collimated jets of ultrarelativistic particles extending to kilo-parsec (kpc) or sometimes up to mega-parsec (Mpc) distances and are known as jetted-AGN [2]. The spectra of these jetted-AGN strongly depend on the viewing angles of their jets with respect to the observer's line of sight and this difference forms the base for their classification. Blazars are the jetted-AGN aligned at very small (≤ 15-20 • ) viewing angles [3][4][5]. Their spectra are fully dominated by non-thermal radiation from the relativistic jets. Blazars are usually subdivided into flat-spectrum radio quasars (FSRQs) and BL Lac objects (BLLs) [e.g. 6,7]. The division between these two subclasses could be due to their different accretion regimes with FSRQs having higher accretion rates above 10 −2 of the Eddington rate [8].
Blazars emit detectable radiations at all EM wavelengths and their broadband spectral energy distributions (SEDs) have characteristic double-humped structures. The first or low energy SED hump, which peaks in Infra-red (IR) to X-rays, is dominated by the synchrotron radiation from the jet. The second or high energy SED hump, peaking at γ−rays (GeV to TeV energies), is produced either by the inverse-Compton (IC) radiation (leptonic model; [9]) or via radiation arising from hadronic processes (hadronic model; [10]). Based on the peak frequency of the first SED hump or the synchrotron peak frequency ν syn [11] subdivided blazars into LSP (or low synchrotron peaked blazars; ν syn ≤ 10 14 Hz), ISP (or intermediate synchrotron peaked blazars; 10 14 < ν syn < 10 15 Hz), and HSP (or high synchrotron peaked blazars; ν syn ≥ 10 15 Hz). A slightly modified range of ν syn for these subclasses of blazars was proposed by [12] as: LSP having log ν syn (Hz) ≤ 14.0, ISP with 14.0 < log ν syn (Hz) ≤ 15.3, and HSP having log ν syn (Hz) > 15.3.
Markarian (or Mrk) 501 (α 2000 = 16h53m52s; δ 2000 = +39 • 45 ′ 37 ′′ ) is one of the nearest (z = 0.034; [13]) bright X-rays emitting high-frequency peaked BL Lac object (HBL 1 ). It was first detected at very high energy (VHE) γ− ray (above 300 GeV) in 1995 by the Whipple telescope [14] that made it the second extragalactic object, after Mrk 421 [15], detected at TeV energies. Mrk 501 has been extensively observed in several multiwavelength campaigns to understand its exact nature which is yet not fully understood (e.g., [16,17]). In a recent X-ray variability study of Mrk 501 with Swift, the source showed the most extreme X-ray flare activity in March-October 2014 during its 11.5 yr long monitoring period [18]. During this epoch, several short-term X-ray flares were detected with their amplitude varied by factors of 1.9-4.7 on weekly timescales. They also found a moderate positive correlation between the X-ray and the TeV fluxes, while no significant correlation was detected between 0.3−300 GeV and optical fluxes.
To date, out of about 3561 2 blazars, TeV emissions have been detected only in 68 3 blazars. Most of the TeV blazars belong to the HBL (52) subclass. It has been reported that the X-ray emissions of the TeV γ−ray emitting blazars are usually highly variable on intraday timescales (i.e., in less than a day time interval; (e.g., [19,20], and references therein)). High flux variability on such small timescales has been one of the most puzzling issues in the field of blazar astronomy as it requires very large energy outputs within small regions. These emitting regions often lie near to SMBH. Study of X-ray flux variability on very short (or intraday) timescales is helpful in probing the size of the emitting regions near SMBH as well as the underlying emission mechanisms (leptonic or hadronic; (e.g., [21])). The shape of X-ray spectra of TeV blazars has also been studied for quite some time as it provides information about the distribution of particles at these energies. The X-ray spectra of TeV blazars Mrk 421 and Mrk 501 were examined by [22][23][24] and they described them in terms of the curved log-parabola model. The X-ray spectral shape of 29 TeV blazars observed with Swift/XRT were investigated by [25] and they found that most of them are well described by the curved log-parabola model. In our earlier work [19], we performed the timing analysis of five TeV blazars (including Mrk 501) using Nuclear Spectroscopic Telescope Array (NuSTAR) data. The number of observations of Mrk 501 studied in earlier work was only four, all taken in 2013. In the last few years more (three) NuSTAR observations of the source were carried out. These three new observations were not studied earlier by any author. The motivation of this work is to examine all the NuSTAR light curves of the TeV blazar Mrk 501 for intraday flux and spectral variability, and to investigate the nature of its 3-79 keV X-ray spectra.
The outline of the manuscript is as follows. We briefly describe the NuSTAR observations and data processing in Section 2 as well as discuss the analysis techniques used for examining variability properties in Section 3. The results are given in Section 4. A summary of the work and discussion of our results are presented in Section 5.

NuSTAR Observations and Data Processing
NuSTAR is a space observatory consisting of two hard X-ray focusing telescopes and the two corresponding co-aligned focal plane modules, referred to as FPMA and FPMB [26]. We downloaded all NuSTAR observations of the TeV blazar Mrk 501 from the HEASARC Data archive 4 . The blazar Mrk 501 was observed with NuSTAR on seven occasions between 13 April, 2013 and 19 April, 2018 and the good exposure times (those corrected for the periods of South Atlantic Anomaly (SAA) passage and Earth occultation) ranged from 9.98 ks to 25.76 ks. In 2013, NuSTAR observed Mrk 501 four times as a part of an 1 HBL belongs to HSP subclass. extensive simultaneous multiwavelength campaign [16]. It was observed two times in 2017 in NuSTAR guest observer program cycle 2, while a single observation was made in 2018 in the extragalactic legacy survey as a part of thw NuSTAR extended mission. The observation log of NuSTAR data for Mrk 501 is given in Table 1. The NuSTAR data sets for the TeV blazar Mrk 501 were processed using the NuSTARDAS software within the HEASOFT 5 software package version 6.26.1. We first generated the calibrated and cleaned level 2 event files using the standard nupipeline script and the updated CALDB files version 20200526. For each observation, the source light curve and spectrum were then extracted from a circular region of radius 50 ′′ centred at the source using the nuproducts script. To extract background data for each observation we used a circular region of the similar radius on the same detector module in which the source was located but free from source contamination.
Since the two NuSTAR detectors, FPMA and FPMB, are co-aligned and nearly identical, we summed their background-subtracted light curves and rebinned them using a bin size of 300 s to get the final light curves. We grouped each source spectra using grppha to ensure at least 20 counts per bin.

Fractional Variance
To quantify the amplitude of variability in the blazar light curves (LCs), we employed the fractional variance which is generally used for examining AGN X-ray LCs (e.g., [27,28]). The fractional variance, F var , is explained in detail in our previous papers ( [19,20], and references therein). For a light curve consisting of N data points, the value of F var is given by: The error in F var is calculated as: where S 2 andx are the sample variance and the arithmetic mean of the LC, respectively, while σ 2 err is the mean square error.

Hardness Ratio
To examine spectral variations in the NuSTAR X-ray emission from the TeV blazar Mrk 501, we also extracted LCs in two energy bands: A soft energy band ranging from 3 to 10 keV and a hard energy band from 10 to 79 keV. We then calculated the hardness ratio (HR) as: where S and H are the individual simultaneous fluxes (in the units of counts/sec) in soft and hard energy bands, respectively. The error in HR is computed as: where σ 2 S , and σ 2 H are the errors in soft and hard energy bands, respectively. Since a hardness ratio simply compares the number of counts observed in two different energy bands, it is a model independent technique to study the spectral variability of the source. We examined variations of the HRs with the total flux to study spectral changes with brightness in the 3-79 keV energy range.

X-ray Flux Variability
We plotted all the NuSTAR light curves of the TeV blazar Mrk 501 in Figure 1. Clear intraday variations were seen on a couple of nights of observations. We examined all the light curves of Mrk 501 for intraday flux variations using fractional variability amplitude F var , discussed in Section 3.1. The results of our intraday variability (IDV) analyses are given in Table 2, where dashes "-" indicate that the sample variances < mean square errors, so no real value of F var can be computed for the LC. The source showed significant flux variations on four out of seven nights. The maximum flux variation was detected on MJD 56420 (Obs. ID: 60002024004) with fractional variability amplitude of 15.53%, while the minimum flux variation was observed on MJD 56485 (Obs. ID: 60002024006) with a F var value of 3.90%. The values of F var were 8.45%, and 5.17% on MJD 56486 and MJD 57870, respectively. No significant IDV variation was detected on MJD 56395, MJD 57897, and MJD 58227. We also examined the soft (3-10 keV) energy band LCs and the hard (10-79 keV) energy band LCs for intraday variability. As can be seen from the Table 2, the amplitude of variability is larger in the higher (10-79 keV) energy band than the lower (3-10 keV) energy band during all the observations with significant IDV detection in 3-79 keV X-ray LCs.

X-ray Spectral Variability
We investigated the spectral variability of 3-79 keV X-ray emission from the TeV blazar Mrk 501 using a model-independent hardness-ratio analysis. We plotted the HR against total (3-79 keV) flux in Figure 2. To search for any systematic variation in HR with flux, we fitted a first-order polynomial in each plot of Figure 2. The values of correlation coefficient (r) and the null hypothesis probability (p) are mentioned in each plot. Assuming Gaussian fluctuations and white noise variations, strong positive correlations (r ≥ 0.5 and p < 0.01) between HR and total count rates were found on MJD 56420 (Obs. ID: 60002024004; r = 0.6, p = 3.1 × 10 −10 ) and MJD 56486 (Obs. ID: 60002024008; r = 0.5, p = 4.3 × 10 −3 ) indicating that the spectra become harder with increasing flux during these observations. No statistically significant correlations were detected in the remaining five observations.

NuSTAR Spectra
For each observation, we simultaneously fitted the spectra from the two NuSTAR detectors, FPMA and FPMB, in the XSPEC 6 version 12.10.1f using χ 2 statistics. To account for cross-calibration uncertainties between these two detectors we have included a multiplicative constant factor in each spectral model whose value is kept fixed to 1 for FPMA and free to vary for FPMB. The values of the constant factor ranged from 0.97 to 1.05 which is within the expected values as suggested by [29].
The spectra were modelled with both a simple power-law (PL) and a log-parabolic (LP) model. The power-law model is defined as: where Γ and K are the photon-index and the normalization, respectively, while F(E) is the flux at energy E. The log-parabolic model is given by [22]: where α is the photon-index at fixed energy E pivot = 10 keV, β is the spectral curvature, and K is the normalization.
We multiplied each spectral model by a phabs component with fixed hydrogen column density, taken from [30], to account for the Galactic absorption.
To choose the best-fitted spectral model between PL and LP we performed the F-test 7 using values of χ 2 and the dof 8 for both models. The results of spectral fitting and the F-test are listed in Table 3. The LP model provides a better fit than the simple PL model if the value of F-statistic > 1 and the corresponding 6 https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/XspecManual.pdf. 7 F-test is available in the XSPEC. 8 Degree of freedom. null hypothesis probability, p < 0.01. We found that for five out of seven NuSTAR observations of Mrk 501 the curved LP model provides a better fit over the PL model. For the last two observations, the NuSTAR spectra of Mrk 501 are well described by the PL model as can be seen from the high values (p > 0.01) of null hypothesis probability. The model-fitted spectra and the data-to-model ratios for each observation are plotted in Figure 3. Table 3. Model fits to the NuSTAR spectra of Mrk 501.

Summary and Discussion
In this paper, we analysed seven archival NuSTAR observations of the TeV blazar Mrk 501 performed between 13 April 2013 and 19 April 2018. In particular, we studied hard (3-79 keV) X-ray flux and spectral variability properties of Mrk 501, and investigated the shape of its X-ray spectra in 3-79 keV energy range. We examined all NuSTAR light curves of Mrk 501 for intraday variability using fractional variability analysis. We found significant IDV variations in four out of seven X-ray LCs with the fractional variability amplitudes ranging from 3.90% to 15.53%. No significant variations were detected in the rest three observations. We noticed that during the last two observations (MJD 57897, and MJD 58227) the source was in somewhat lower states, while on MJD 56395 the intrinsic variations in the light curve are probably smaller than the measurement errors.
TeV blazars are known to exhibit large amplitude flux variations at all frequencies on different timescales (e.g., [19][20][21][31][32][33], and references therein). Flux variations in blazar LCs are generally understood to originate from the relativistic jets [34]. At high frequencies (X-rays to γ−rays) large amplitude variations are often observed on very short timescales indicating compact emitting regions (e.g., [35,36]). Such rapid variations are often explained by the interaction of relativistic turbulent plasma with shock within the jets, which accelerates the particles to high energies [34]. The shortest hard X-ray (3-79 keV) variability with doubling time of ∼ 14 minutes was reported from the TeV blazar Mrk 421 by [21] during its 2013 April flare. They suggested that magnetic reconnections taking place in fast-moving emitting regions within the jets could be responsible for such a short timescale variability ("jets-in-a-jet" model; [37]). Long term flux variations in the X-ray light curve of Mrk 501 during the entire NuSTAR monitoring period can be explained by the shock-in-jet model which involves the acceleration of particles to high energies by an internal shock within the jet followed by subsequent cooling by emitting radiations [38]. Other possible explanations for the longer timescale variations include a change in the viewing angle, and hence in the Doppler factor of the emitting region within the jet (e.g., [39,40]), and/or a change in the magnetic field [41].
We examined the hard (3-79 keV) X-ray spectral variability of the TeV HBL Mrk 501 using HR analysis. We found that for observations on MJD 56420 and MJD 56486, the HR increases with increasing flux, that is, the NuSTAR spectra of Mrk 501 became harder with the increasing count rates. Such a "harder-when-brighter" trend has often been noticed in X-ray observations of the HBL-type blazars (e.g., [19,42,43], and references therein). For the other five observations, no strong correlation were seen between HR and count rates. Moreover, as can be seen from Figure 2, the data have uncertainties which were not taken into account during the correlation analysis and thus, the correlation between the HR and the flux may be somewhat lower than those reported here.
We also studied the shape of NuSTAR spectra of TeV blazar Mrk 501 using the simple power-law model and the log-parabola model. We found that five out of seven NuSTAR spectra are well described by the curved LP models with local photon indices α lying in the range of 2.12-2.32 and the spectral curvature β ∼ 0.15-0.28. For the last two observations (Obs. IDs: 60202049004 and 60466006002), when the X-ray fluxes are relatively low the simple PL model provides an equivalently good fit with photon indices Γ ∼ 2.70 and 2.75, respectively. During the NuSTAR monitoring period of Mrk 501, the maximum 3-79 keV unabsorbed flux recorded was 35.31 × 10 −11 erg cm −2 s −1 on MJD 56485, while the minimum flux observed was 1.25 × 10 −11 erg cm −2 s −1 on MJD 58227.
The study of X-ray spectral shape of TeV blazars is helpful in understanding the particle acceleration mechanism and the distribution of emitting particles. The curved LP model was used for the first time to describe the spectral shape of synchrotron emission from the BL Lac type blazars by [44]. Later, [22][23][24] successfully described the X-ray spectra of TeV HBLs Mrk 421 and Mrk 501 using the LP model and suggested that the curved X-ray spectra of TeV blazars can be understood in terms of statistical energy-dependent particle acceleration. The observed curvature of X-ray spectra of TeV HBLs could be either convex or concave. The convex curvature of the X-ray spectra, as found in this work, is likely to be caused by a single-accelerated particle distribution (e.g., [22]), while the concave X-ray spectra observed in some studies (e.g., [45]), maybe a consequence of the spectral upturn at the interaction of the high-energy tail of the synchrotron emission and the low-energy part of the IC emission.