X-Ray Intraday Variability of Five TeV Blazars with NuSTAR

We have examined the emission of all of the five TeV HBLs observed by NuSTAR: 1ES 0229+200, Markarian (Mrk) 421, Mrk 501, 1ES 1959+650, and PKS 2155−304. This gave us a new opportunity to study the variability nature of these blazars in hard X-ray energies. In order to search for IDV, we only selected observations for which the good exposure times exceeded 5 ks. We downloaded all 43 data sets publicly available from the HEASARC Data Archive⁶ , and our minimum temporal constraint yielded 3 observations for 1ES 0229+200, 22 for Mrk 421, 4 for Mrk 501, 2 for 1ES 1959+650, and 9 for PKS 2155−304. These 40 observations were made between 2012 July 7 and 2014 September 22, and the good exposure times ranged from 5.76 to 57.51 ks. The observing log of NuSTAR data for these five TeV blazars is given in Table 1.

We reduced and analyzed the NuSTAR data using HEASOFT⁷ version 6.17 and CALDB version 20151008. We used the standard nupipeline script to generate calibrated, cleaned, and screened level 2 event files. Each source LC is then extracted, using the nuproducts script, from a circular region centered at the source. We employ a circular background region that is selected to be both relatively close to the source but also far enough away to be free from contamination by the source. The radii of the source and background regions for our five TeV HBLs are listed in Table 2. Since NuSTAR has two co-aligned and nearly identical detectors, FPMA and FPMB, their count rates are background-subtracted and summed to generate the final LCs. We used a bin size of 5 minutes to extract the finer binned LCs. This sampling interval is similar to those employed in most ground-based optical IDV studies of AGNs.

Blazars are known to exhibit rapid and strong flux variations on diverse timescales across the complete electromagnetic spectrum. The strength of this variability is often quantified by calculating the excess variance, which is a measure of a source's intrinsic variance. The excess variance is calculated by subtracting the variance arising from measurement errors from the total variance of the observed LC. If an LC consisting of N measured flux values, x_i, contains corresponding finite uncertainties ${\sigma }_{\mathrm{err},i}$ arising from measurement errors, then the excess variance is calculated (e.g., Vaughan et al. 2003) as

$$\begin{eqnarray}&&{\sigma }_{\mathrm{XS}}^{2}={S}^{2}-\overline{{\sigma }_{\mathrm{err}}^{2}},\end{eqnarray} \tag{ 1 }$$

where $\overline{{\sigma }_{\mathrm{err}}^{2}}$ is the mean square error, given by

$$\begin{eqnarray}&&\overline{{\sigma }_{\mathrm{err}}^{2}}=\displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}{\sigma }_{\mathrm{err},i}^{2}.\end{eqnarray} \tag{ 2 }$$

The quantity S² is the sample variance of the LC, and is given by

$$\begin{eqnarray}&&{S}^{2}=\displaystyle \frac{1}{N-1}\displaystyle \sum _{i=1}^{N}{({x}_{i}-\bar{x})}^{2},\end{eqnarray} \tag{ 3 }$$

where $\bar{x}$ is the arithmetic mean of x_i.

The normalized excess variance is ${\sigma }_{\mathrm{NXS}}^{2}={\sigma }_{\mathrm{XS}}^{2}/{\bar{x}}^{2}$, and the fractional rms variability amplitude, F_var, which is the square root of ${\sigma }_{\mathrm{NXS}}^{2}$, is thus

$$\begin{eqnarray}&&{F}_{\mathrm{var}}=\sqrt{\displaystyle \frac{{S}^{2}-\overline{{\sigma }_{\mathrm{err}}^{2}}}{{\bar{x}}^{2}}}.\end{eqnarray} \tag{ 4 }$$

The uncertainty on F_var is given by (e.g., Vaughan et al. 2003)

$$\begin{eqnarray}&&\mathrm{err}({F}_{\mathrm{var}})=\sqrt{{\left(\sqrt{\displaystyle \frac{1}{2N}}\displaystyle \frac{\overline{{\sigma }_{\mathrm{err}}^{2}}}{{\bar{x}}^{2}{F}_{\mathrm{var}}}\right)}^{2}+{\left(\sqrt{\displaystyle \frac{\overline{{\sigma }_{\mathrm{err}}^{2}}}{N}}\displaystyle \frac{1}{\bar{x}}\right)}^{2}}.\end{eqnarray} \tag{ 5 }$$

The advantage of using a discrete correlation function (DCF) over a classical correlation function is that it can be applied to unevenly sampled data sets, as are typical of astronomical observations and as we have here, without interpolating between data (Edelson & Krolik 1988). The DCF is calculated in the following fashion. First, we calculate the set of unbinned discrete correlations for two discrete data sets x and y (e.g., Gaur et al. 2015) as

$$\begin{eqnarray}&&{\mathrm{UDCF}}_{{ij}}=\displaystyle \frac{({x}_{i}-\bar{x})({y}_{j}-\bar{y})}{\sqrt{{\sigma }_{x}^{2}{\sigma }_{y}^{2}}},\end{eqnarray} \tag{ 6 }$$

where x_i and y_j are the data points, $\bar{x}$ and $\bar{y}$ are their means, and ${\sigma }_{x}$ and ${\sigma }_{y}$ are their standard deviations, respectively. Each of these is associated with the pairwise lag ${\rm{\Delta }}{t}_{{ij}}={t}_{j}-{t}_{i}$. After calculating the UDCF, the correlation function is binned in time. The method does not automatically define a bin size so one must investigate several values for this parameter. If the bin size is too large, information is lost, but if the bin size is too small, we can get spurious correlations, and the timescales may be difficult to interpret. Simulations conducted by Edelson & Krolik (1988) suggest that the results depend only weakly on the choice of bin size.

After binning, the DCF can be calculated by averaging the UDCF values (M in number), for which $\tau -{\rm{\Delta }}\tau /2\leqslant {\rm{\Delta }}{t}_{{ij}}\,\lt \tau +{\rm{\Delta }}\tau /2$, as

$$\begin{eqnarray}&&\mathrm{DCF}(\tau )=\displaystyle \frac{1}{M}\sum {\mathrm{UDCF}}_{{ij}}.\end{eqnarray} \tag{ 7 }$$

Edelson & Krolik (1988) defined the standard error for each bin as

$$\begin{eqnarray}&&{\sigma }_{\mathrm{DCF}}(\tau )=\displaystyle \frac{\sqrt{\sum {[{\mathrm{UDCF}}_{{ij}}-\mathrm{DCF}(\tau )]}^{2}}}{M-1}.\end{eqnarray} \tag{ 8 }$$

In general, a positive DCF peak means that the two data signals are correlated, while a negative DCF peak means that the two data sets are anti-correlated, but no DCF peak, or DCF = 0, means that no correlation exists between the two data sets. When correlating a data series with itself (i.e., x = y), we obtain the ACF with an automatic peak at $\tau =0$, indicating the absence of any time lag. For an ACF, any other strong peak can indicate the presence of periodicity, but strong dips provide an indication of the presence and value of a timescale in the data (Rani et al. 2011).

To examine the spectral variability of the NuSTAR X-ray emission from the five TeV HBLs, we split the LCs into two energy bands: a soft band from 3 to 10 keV and a hard band from 10 to 79 keV. We then computed the hardness ration (HR) as

$$\begin{eqnarray}&&\mathrm{HR}=\displaystyle \frac{(H-S)}{(H+S)},\end{eqnarray} \tag{ 9 }$$

where H and S are the net count rates in the hard (10–79 keV) and soft (3–10 keV) bands, respectively. The HR is a commonly used and simple model-independent method to study spectral variations. We examined variations of the HRs with time to search for spectral changes over this broad X-ray band.

The BL Lac object 1ES 0229+200 (${\alpha }_{2000}$ = 02^h32^m532; ${\delta }_{2000}=+20^\circ 16^{\prime} 21^{\prime\prime}$) is an HBL at $z=0.1396$ (Woo et al. 2005). Recently, Cologna et al. (2015) reported flux variability at VHE on monthly and yearly timescales by studying long-term observations of 1ES 0229+200 (from 2004 to 2013) with H.E.S.S. and also found a hint of correlation between X-ray and VHE emissions.

NuSTAR observed 1ES 0229+200 for 16.26, 20.29, and 18.02 ks on 2013 October 2, 5, and 10, respectively. These LCs and their ACFs are plotted in Figures 1(a) and 3(a), respectively. The count rates are low and the data noisy, so all fractional variances are consistent with no detectable variability on any of those days in the entire 3–79 keV energy band. However, variations in hard and soft bands are nominally seen on 2013 October 2 and 10, respectively. Unsurprisingly, the ACFs show no hint of an IDV timescale in the NuSTAR range.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The soft and hard LCs (left panel), HR plots (middle panel), and the DCF plots between soft and hard band (right panel) of 1ES 0229+200 are shown in Figure 2(a). No significant spectral changes are seen from the HR plots in any of the three observations. Naturally, given the lack of significant variations, the DCF plots are all flat and consistent with 0 throughout.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Markarian 421 (Mrk 421; ${\alpha }_{2000}$ = 11^h04^m19^s; ${\delta }_{2000}=+38^\circ 11^{\prime} 41^{\prime\prime}$) is one of the nearest ($z=0.031$) BL Lac objects and was the first extragalactic object detected at TeV energy (Punch et al. 1992). It is classified as an HBL because its synchrotron peak lies at soft X-rays.

Mrk 421 is highly variable at all timescales over the entire electromagnetic spectrum and has been extensively studied during its flaring states (e.g., Takahashi et al. 1994, 1995; Kerrick et al. 1995; Fan & Lin 1999; Gupta et al. 2004; Teräsranta et al. 2004, 2005; Costa et al. 2008; Pittori et al. 2008; Smith et al. 2008; Gaur et al. 2012c; Lico et al. 2012; Blasi et al. 2013; Racero & de la Calle 2013; Abdo et al. 2014). In the optical region, variations of 4.9 mag were found over the course of several years (Stein et al. 1976). In the spring to summer of 2006, a major flare was recorded with the flux reaching up to ∼8.5 mCrab in the 2.0–10.0 keV energy range (Tramacere et al. 2009; Ushio et al. 2009). In 2010 January and February, strong X-ray flares were detected, with the maximum flux recorded to be 120 ± 10 mCrab and 164 ± 17 mCrab, respectively, with the latter being the largest ever reported from the source (Isobe et al. 2015). In 2012 and 2013, Mrk 421 displayed two flares (Hovatta et al. 2015), and the gamma-ray flare in 2012 was observed without a simultaneous X-ray flare, a behavior called an "orphan flare" (Fraija et al. 2015). Due to its variable nature at all wavelengths and its proximity, it has been studied in several multi-wavelength observational campaigns (e.g., Błażejowski et al. 2005; Fossati et al. 2008; Gupta et al. 2008a; Horan et al. 2009; Gaur et al. 2012c; Aleksić et al. 2015; Li et al. 2016; MAGIC Collaboration et al. 2016). Recently, Sinha et al. (2016) reported on a long-term study of Mrk 421 with the High Altitude Gamma Ray telescope array at Hanle, India, and found strong correlations between the gamma and radio wavelengths and between the optical and gamma wavebands, but saw no correlation between gamma and X-ray emissions. They also saw that the variability depends on energy, being maximum in X-ray and VHE bands.

Mrk 421 was first observed with NuSTAR twice during 2012 July as part of the calibration process for the telescope. It was then observed many times in 2013 January–April as part of an extensive and intensive multi-wavelength campaign that involved simultaneous or quasi-simultaneous data obtained in the radio, optical, soft X-ray, hard X-ray, and GeV through TeV gamma-ray bands (Baloković et al. 2016). The results of the first part of this study, covering the the low-flux state (in 2013 January–March) were reported by Baloković et al. (2016) and Kataoka & Stawarz (2016). Paliya et al. (2015) presented a variability analysis, and Sinha et al. (2015) presented analyses of the spectral variations during the strong flaring state seen in 2013 April using this NuSTAR data.

The IDV LCs of Mrk 421 that we have re-reduced and plotted in Figure 1 indicate that variations appear to seen during every one of the 22 observations. Our overall fractional variability amplitudes (F_var) given in Table 3 confirm these impressions and thus echo the results of Baloković et al. (2016) and Paliya et al. (2015) in this regard. The fractional variability values for soft and hard bands, given in Table 3, indicate that the variability in the hard band is stronger than that in the soft band for 19 out of 22 observations. The complete LC obtained from our analysis of the NuSTAR data from Mrk 421 is shown in Figure 4. Very strong flares were seen while the source was monitored nearly continuously during 2013 April 10–16 (Paliya et al. 2015). The unprecedented outburst during this period is marked with a box in Figure 4, and we have plotted a zoomed-in version of that box as an inset to the same figure. This is a double-peaked outburst in which the first flare appears to have a nearly a Gaussian shape with peak flux at ∼MJD 56395, while the second one, occurring two days later, is even stronger, and evinces a very sharp rise and decay.

Table 3.  X-Ray Variability Parameters

Blazar Name	Obs. ID	${F}_{\mathrm{var}}(\mathrm{percent})$			ACF(ks)	Bin Size(ks)
		Soft (3–10 keV)	Hard (10–79 keV)	Total (3–79 keV)
1ES 0229+200	60002047002	⋯	15.632 ± 4.757	⋯	⋯	2.00
	60002047004	⋯	⋯	⋯	⋯	2.00
	60002047006	8.638 ± 2.557	⋯	4.460 ± 3.411	⋯	2.00
MRK 421	10002015001	21.105 ± 0.341	26.899 ± 1.179	21.217 ± 0.260	⋯	5.00
	10002016001	31.225 ± 0.337	48.503 ± 1.227	32.314 ± 0.305	⋯	3.00
	60002023002	7.531 ±1.183	68.747 ± 8.125	9.159 ±1.029	4.9	0.90
	60002023004	12.173 ± 0.874	75.723 ± 13.099	11.485 ± 0.860	19.9	1.50
	60002023006	20.051 ± 0.402	39.626 ± 1.695	21.595 ± 0.392	11.4	3.00
	60002023008	13.566 ± 0.721	31.496 ± 9.967	10.275 ± 0.757	⋯	3.50
	60002023010	10.371 ± 0.517	8.027 ±3.382	10.264 ± 0.411	19.5	3.00
	60002023012	20.310 ± 0.463	33.744 ± 1.475	21.556 ± 0.424	⋯	3.00
	60002023014	25.203 ± 1.219	21.537 ± 11.947	26.702 ± 0.652	⋯	3.00
	60002023016	12.497 ± 0.550	13.995 ± 2.175	11.815 ± 0.421	⋯	3.00
	60002023018	14.707 ± 0.453	14.930 ± 2.339	14.583 ± 0.420	⋯	3.00
	60002023020	10.217 ± 1.284	⋯	9.164 ±0.374	2.5	0.30
	60002023022	22.052 ± 0.386	40.386 ± 1.243	24.265 ± 0.193	32.8	3.00
	60002023024	15.181 ± 0.431	17.344 ± 1.401	15.458 ± 0.411	⋯	1.00
	60002023025	59.837 ± 0.131	64.919 ± 0.351	60.497 ± 0.123	⋯	6.00
	60002023027	13.406 ± 0.196	17.386 ± 0.486	13.945 ± 0.182	⋯	1.00
	60002023029	25.285 ± 0.184	26.116 ± 0.517	22.266 ± 0.175	13.1	1.00
	60002023031	31.141 ± 0.108	36.316 ± 0.229	32.073 ± 0.098	⋯	2.00
	60002023033	19.344 ± 0.540	25.137 ± 0.458	18.747 ± 0.169	⋯	2.00
	60002023035	39.043 ± 0.178	44.841 ± 0.427	39.943 ± 0.165	⋯	2.00
	60002023037	17.584 ± 0.404	24.000 ± 1.382	18.268 ± 0.387	12.6	1.00
	60002023039	12.310 ± 0.460	17.255 ± 1.798	12.773 ± 0.443	⋯	2.00
MRK 501	60002024002	1.327 ± 1.406	⋯	1.391 ± 1.121	⋯	1.00
	60002024004	14.572 ± 0.627	20.901 ± 1.940	15.540 ± 0.354	⋯	5.00
	60002024006	4.278 ± 0.432	5.138 ± 0.925	3.743 ± 0.351	⋯	2.00
	60002024008	6.338 ± 0.670	9.067 ± 1.400	8.284 ± 0.358	8.0	1.00
1ES 1959+650	60002055002	26.849 ± 1.192	35.449 ± 2.876	26.629 ± 0.451	⋯	3.00
	60002055004	⋯	6.522 ± 2.963	⋯	⋯	1.00
PKS 2155-304	10002010001	12.901 ± 0.915	33.771 ± 5.960	10.601 ± 0.744	29.6	2.00
	60002022002	⋯	21.275 ± 5.478	11.020 ± 1.109	57.4	5.00
	60002022004	12.551 ± 2.437	21.333 ± 7.654	14.931 ± 1.400	⋯	3.00
	60002022006	10.830 ± 4.268	⋯	9.256 ± 2.162	⋯	2.00
	60002022008	⋯	28.233 ± 15.966	13.504 ± 3.000	⋯	3.00
	60002022010	⋯	12.931 ± 15.148	⋯	⋯	3.00
	60002022012	20.963 ± 1.300	24.613 ± 3.958	20.916 ± 1.160	⋯	4.00
	60002022014	16.726 ± 2.082	⋯	16.992 ± 1.554	⋯	4.00
	60002022016	16.810 ± 2.876	⋯	3.298 ± 8.251	⋯	2.00

Download table as:  ASCIITypeset image

The LCs in soft and hard bands (left panel), HR plots (middle panel), and the DCF plots (right panel) are shown in Figures 2(a), (b), (c), and (d). The variations of HR with time show clear features of spectral changes that are stronger during flares. The spectra harden with the increasing flux, providing evidence for a general "harder when brighter" behavior of blazars. Similar X-ray variability behavior was already seen in XMM-Newton observations (e.g., Brinkmann et al. 2003; Ravasio et al. 2004) and in NuSTAR observations (e.g., Paliya et al. 2015; Baloković et al. 2016) of Mrk 421. As can be seen from DCF plots, the soft and hard bands are positively correlated with zero time lag, indicating that their emissions come from the same emitting region at the same time.

The ACFs of Mrk 421 are plotted in Figure 3, and eight of them show structures indicative of timescales. The dates on which observations began for which these structures were seen in the ACFs are 2013 January 2, 10, and 15; 2013 February 6; 2013 March 17; and 2013 April 2, 13, and 18, and the corresponding putative "timescales," which range from 2.5 to 32.8 ks, are given in Table 3. The other ACFs do not show any variability timescales or are too noisy to allow any detections of them.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The TeV blazar Markarian 501 (Mrk 501; ${\alpha }_{2000}$ = 16^h53^m52^s; ${\delta }_{2000}=+39^\circ 45^{\prime} 37^{\prime\prime}$), at $z=0.034$, was the second extragalactic object detected at TeV energies. It was first detected at energy greater than 300 GeV by the Whipple Observatory (Quinn et al. 1996). In 1997, Mrk 501 went into a surprisingly high state with the flux recorded up to 10 Crab at energies >1 TeV, and it displayed strong VHE variability (Aharonian et al. 1997, 1999a, 1999b; Catanese et al. 1997; Samuelson et al. 1998). On 1997 April 16, the highest VHE flux ($F\sim 8.3\times {10}^{-10}$ cm⁻² s⁻¹ at energies > 250 GeV) ever was recorded (Djannati-Atai et al. 1999). Fast VHE flux variations with a flux doubling time of ∼2 minutes were observed from Mrk 501 in 2005 July (Albert et al. 2007). Neronov et al. (2012) reported an orphan VHE gamma-ray flare in 2009 that was not accompanied by an X-ray flare. In the optical band, Mrk 501 has also shown flux variability on different timescales (e.g., Gupta et al. 2008b, 2012, 2016a; Xiong et al. 2016). Recently, in 2014 June, H.E.S.S. observed major flaring activity when the flux reached over 1 Crab and rapid flux variability was recorded at VHE (∼2–20 TeV; Chakraborty et al. 2015).

Mrk 501 was observed with NuSTAR on four occasions between 2013 April 13 and 2013 July 13 as part of an extensive multi-wavelength campaign. Data were collected from radio through optical, UV, X-ray, and several gamma-ray bands, with the last two observations made as part of a target of opportunity program because of elevated states detected by other telescopes (Furniss et al. 2015). The mean fluxes varied significantly from observation to observation, and the LCs in Figure 1(b) show clear IDV during the last three nights of observations. This is confirmed by the substantial values of F_var given in Table 3 for those observations. The F_var values in the soft and hard bands suggest greater variations in the hard band.

The soft and hard LCs of Mrk 501 are plotted in the left panels of Figure 2(d). As seen from the HR plots in the middle panels of that figure, no significant spectral changes were seen during for first observation, while for the other three observations the spectra become harder with increasing flux and softer with decreasing flux. Such spectral behavior was also reported in other X-ray observations (e.g., Pian et al. 1998; Aliu et al. 2016). The DCF plots shown in the right panels of Figure 2(d) show a positive correlation with zero lag between the soft and hard bands, except for the first observation, during which there was negligible variability. Hence, for the first observation, the ACF plot given in Figure 3(b) is essentially noisy. While the second and third observations show clear ACFs, those data do not give any indication of a variability timescale. The ACF of the last observation, however, does indicate a possible timescale of ∼8.0 ks.

The HBL 1ES 1959+650 (${\alpha }_{2000}$ = 19^h59^m598; ${\delta }_{2000}\,=+65^\circ 08^{\prime} 55^{\prime\prime}$) is a BL Lac object at $z=0.48$ (Perlman et al. 1996). This blazar was first detected at X-rays with the Einstein IPC Slew Survey (Elvis et al. 1992) and was further observed with ROSAT in 1996 and with BeppoSAX in 1997 (Beckmann et al. 2002), as well as with RXTE-ARGOS, Swift, and XMM-Newton (Giebels et al. 2002; Tagliaferri et al. 2003; Massaro et al. 2008).

1ES 1959+650 was observed with NuSTAR on 2014 September 17 and 22, with good exposure times of 19.61 ks and 20.34 ks, respectively. The LCs and ACFs are shown in Figures 1(b) and 3(b), respectively. During the first observation, the flux increased dramatically over the span of 30 ks. The count rate was both lower and steadier during the second measurements. The ACF plot for the observation on 2014 September 17 does not show a timescale, while the observation on 2014 September 22 yields a noisy ACF and provides no useful information.

The soft and hard band LCs (left panels), HR plots (middle panels), and the DCF between the two energy bands (right panels) for this blazar are plotted in Figure 2(d). The HR plots reveal no spectral variations. The DCF plot for the observation on 2014 September 17 exhibits a correlation between the soft and hard bands with zero lag, while the lack of variability naturally means that no such correlation can be seen for the observation on 2014 September 22.

The HBL PKS 2155−304 (${\alpha }_{2000}$ = 21^h58^m527; ${\delta }_{2000}=-30^\circ 13^{\prime} 18^{\prime\prime} ;z=0.116$; Falomo et al. 1993; Farina et al. 2016) is the brightest BL Lac object in UV to TeV energies in the southern hemisphere. It was first recognized as a TeV blazar by the Durham MK6 telescopes (Chadwick et al. 1999). VHE flux variability on timescales of minutes was also found (Aharonian et al. 2007). PKS 2155−304 has been observed at all wavelengths, and flux variability on diverse timescales has been reported by many authors (e.g., Urry et al. 1993; Gaur et al. 2010; Gupta 2011; Zhang et al. 2014; Chevalier et al. 2015; Bhagwan et al. 2016; Sandrinelli et al. 2016 and references therein). Lachowicz et al. (2009) reported a possible ∼4.6 hr quasi-periodic oscillation in an X-ray LC using XMM-Newton data.

NuSTAR observed PKS 2155−304 on nine occasions between 2012 July 8 and 2013 September 28, where the exposures ranged from 10.53 to 45.06 ks, as given in Table 1. All these LCs and ACFs are plotted in Figures 1(b) and 3(b), respectively. Although the count rates are low, a visual inspection of the LCs indicates that IDV appears to be present on most dates. The F_var values and their errors in the 3–79 keV range given in Table 3 confirm that IDV was detected during seven of the nine observations. The F_var values in soft and hard bands are also given in Table 3. The ACF plots for observations on 2012 July 8 and 2013 April 23 indicate variability timescales of 29.6 and 57.4 ks, respectively. The other ACF plots are noisy or do not show any timescale of variability.

The soft and hard band LCs (left panel), HR plots (middle panel), and the DCF plots (right panel) of PKS 2155−304 are shown in Figures 2(d) and (e). The HR plots are noisy, providing no useful information. The DCFs for observations on 2012 July 8, 2013 July 16, and 2013 August 26 show zero lag correlations between the bands, while for other observations no significant correlations are observed.

One of the important features of blazar emissions at all measured wavelengths is rapid and strong flux variability on diverse timescales. By doing careful blazar variability studies we can better understand the radiation mechanisms and also get information about the size, location, and structure of the emitting region (e.g., Ciprini et al. 2003). The theoretical models proposed to explain the totality of intrinsic AGN variability can be coarsely classified as the relativistic-jet-based models (e.g., Marscher & Gear 1985; Gopal-Krishna & Wiita 1992; Marscher 2014; Calafut & Wiita 2015) and the accretion-disk-based models (e.g., Chakrabarti & Wiita 1993; Mangalam & Wiita 1993). In blazars, the accretion disk radiation is almost always dominated by the Doppler boosted radiation from the relativistic jets, so the accretion-disk-based models are generally not able to explain variability on any measurable timescales. However, in radio-quiet quasars and for blazars in very low states, it is possible that the IDV and STV can be explained by instabilities or hot spots on the accretion disk (e.g., Chakrabarti & Wiita 1993; Mangalam & Wiita 1993) because the accretion disk flux is not swamped by the jet flux then.

The variations on LTV timescales (months to years) for blazars usually can be explained by shock-in-jet models (e.g., Marscher & Gear 1985; Wagner & Witzel 1995), in which relativistic shocks, assumed to be formed by the disturbances in the inner portion of the jet, propagate outward along the jet and produce major flux changes. Motions of the shock through a helical jet (or helical structures within the jet) can cause variations in the Doppler boosting by changing the effective viewing angle and thus explain some LTV and STV (e.g., Camenzind & Krockenberger 1992; Gopal-Krishna & Wiita 1992). Smaller variations on STV and IDV timescales can be explained by turbulence behind the shock in the relativistic jet (e.g., Marscher 2014; Calafut & Wiita 2015; Pollack et al. 2016).

In this work, we found minimum hard X-ray variability timescales of 2.5, 8.0, and 29.6 ks for Mrk 421, Mrk 501, and PKS 2155−304, respectively. Such rapid X-ray variability timescales have been reported earlier in these sources (Catanese & Sambruna 2000; Kataoka et al. 2000; Cui 2004). Extremely fast TeV variability timescales of a few minutes, as detected in PKS 2155−304 (Aharonian et al. 2007) and Mrk 501 (Albert et al. 2007) are even shorter than the light crossing time of the Schwarzschild radius of the supermassive black holes at the centers of these blazars. Such rapid timescales require the TeV flux emitting region to be very compact; however, the requirement that the TeV photons actually escape such a compact region without being absorbed via pair creation with the synchrotron photons implies the Lorentz factor (Γ) of the emitting region must be $\gtrsim 50$ (e.g., Gopal-Krishna et al. 2006; Begelman et al. 2008). There is very little other evidence for these extremely high bulk Lorentz factors in any AGNs, and they appear to be even more problematical for TeV blazars because much lower values of Γ in these blazars have been inferred from the rather slow apparent motions of their radio knots made using very long baseline arrays (Giroletti et al. 2004; Piner & Edwards 2004).

Models that involve the radio and TeV emission emitting regions having substantially different properties can avoid this apparent contradiction. For example, Ghisellini & Tavecchio (2008) proposed the "needle" model to explain fast TeV variability through the localized magneto-centrifugal acceleration of beams of electrons at small pitch angles along the magnetic field lines. Their model predicts TeV flares without simultaneous X-ray flares, and thus can explain the so-called "orphan" flares. However, during the 2013 April outburst (Figures 1 and 4) very fast (∼14 minutes) hard X-ray variability is seen in the NuSTAR LCs of Mrk 421 (Paliya et al. 2015), accompanied by variations at VHE γ-rays (Cortina & Holder 2013). So the "needle" model fails to explain the spectacular 2013 April outburst of Mrk 421. Those flares can, however, be explained with a "jets-in-a-jet" model proposed by Giannios et al. (2009), which relies on magnetic reconnection. This model predicts rapid TeV variations along with fast X-ray variability. In this model, the relativistic outflow of material from the magnetic reconnection regions gives rise to rapid flares through synchrotron-self-Compton emission. The model also has the flexibility to produce more slowly varying flares due to multiple reconnection regions or the tearing of a large reconnection site. The presence of more slowly varying short-term flares in the complete NuSTAR LC of Mrk 421 (Figure 4) lends some further support to this type of "jets-in-a-jet" model, as discussed by Paliya et al. (2015).

Zoom In Zoom Out Reset image size

Assuming that the hard X-ray emission from the HBLs is mainly due to synchrotron emission, certain parameters can be estimated in a fashion that does not depend on the details of the acceleration models. We initially follow Paliya et al. (2015) and recall that the synchrotron cooling timescale, in the observer's frame, of an electron with energy $E=\gamma {m}_{e}{c}^{2}$ is (Zhang et al. 2002)

$$\begin{eqnarray}&&{t}_{\mathrm{cool}}(\gamma )\simeq 7.74\times {10}^{8}\displaystyle \frac{(1+z)}{\delta }{B}^{-2}{\gamma }^{-1}\,{\rm{s}},\end{eqnarray} \tag{ 10 }$$

where δ is the bulk Doppler factor, B is the magnetic field strength in Gauss, and γ is the electron Lorentz factor. For a given magnetic field strength and electron Lorentz factor the frequency at which synchrotron emission occurs is (e.g., Paliya et al. 2015),

$$\begin{eqnarray}&&\nu \equiv {\nu }_{19}\times {10}^{19}\,\mathrm{Hz}\simeq 4.2\times {10}^{6}\displaystyle \frac{\delta }{1+z}B{\gamma }^{2},\end{eqnarray} \tag{ 11 }$$

where $0.08\lt {\nu }_{19}\lt 2$ for X-rays in the NuSTAR band. Combining these two equations with the physical requirement that the observed minimum variability timescale has to be larger than or equal to the synchrotron cooling timescale, we get for Mrk 421 (with our minimum ${t}_{\mathrm{var}}=2500\,{\rm{s}}$ and $z=0.031$)

$$\begin{eqnarray}&&B\geqslant 0.35\,{\delta }^{-1/3}{\nu }_{19}^{-1/3}\,{\rm{G}}.\end{eqnarray} \tag{ 12 }$$

We note that this expression has different dependences on δ and ${\nu }_{19}$ than does the one (Equation (5)) in Paliya et al. (2015). For $\delta =25$, we find $B\geqslant 0.12{\nu }_{19}^{-1/3}$ G, which is close to the typical value of $B\sim 0.1$ G that is inferred from the SED modeling of Mrk 421. Assuming $\delta =25$ and $B\geqslant 0.12$ G, we can constrain the electron Lorentz factor to

$$\begin{eqnarray}&&\gamma \leqslant 9\times {10}^{5}{\nu }_{19}^{1/2}.\end{eqnarray} \tag{ 13 }$$

We can also estimate the characteristic size of the emitting region as

$$\begin{eqnarray}&&R\leqslant {{ct}}_{\mathrm{var}}\delta /(1+z)\leqslant 1.8\times {10}^{15}\,\mathrm{cm}.\end{eqnarray} \tag{ 14 }$$

Similarly, we can constrain these parameters for Mrk 501 and PKS 2155−304, where we also have reasonable observed timescales from NuSTAR data, where typical values for δ are used (Kataoka et al. 2000; Furniss et al. 2015) and where ${\nu }_{19}=1$ is assumed. These are given in Table 4.

The rapid hard X-ray variability observed by NuSTAR suggests that the relativistic electrons responsible for the hard X-ray emission of TeV HBLs must be repeatedly injected (accelerated), since the high-energy electrons have short cooling timescales. The relativistic electrons are known to be accelerated at shock fronts within jets, although additional acceleration regions are possible. Diffusive shock acceleration (e.g., Blandford & Eichler 1987) is an electron acceleration mechanism that could be responsible for both the flux variations and the observed spectral hardening at high energies.

We examined the X-ray spectral variability of five TeV HBLs in this work using a model-independent HR analysis. Although this method does not provide any direct information about the physical parameters responsible for spectral changes, it is the simplest way to study the spectral variability. Given the steepness of the X-ray spectrum, which means that that the number of hard counts is significantly lower than that of the soft counts (see Figure 2), it is difficult to do a more detailed spectral analysis. For two bright HSP blazars, Mrk 421 and Mrk 501, we found that the HR increases with increasing count rates, that is, their spectra become flatter with increasing flux. Such a "hardening when brightening" trend appears to be a general feature of HSP-type blazars, as it was also noticed in earlier X-ray observations (e.g., Pian et al. 1998; Zhang et al. 2002; Brinkmann et al. 2003; Ravasio et al. 2004). In the case of Mrk 421, we observed that the variations in HR values with time are particularly large during flares.