A SEARCH FOR SPECTRAL HYSTERESIS AND ENERGY-DEPENDENT TIME LAGS FROM X-RAY AND TeV GAMMA-RAY OBSERVATIONS OF Mrk 421

VERITAS is an array of four 12 m ground-based imaging atmospheric Cherenkov telescopes in southern Arizona, each equipped with a camera consisting of 499 photomultiplier tubes (PMTs) (Holder 2011). It is sensitive to gamma-rays in the energy range from ∼100 GeV to ∼30 TeV with an energy resolution of ∼15%, and covers a 35 field-of-view with an angular resolution (68% containment) of ∼01. It is capable of making a detection at a statistical significance of five standard deviations (5σ) of a point source of 1% C. U. in ∼25 hr. The systematic uncertainty on the energy calibration is estimated at 20%, and that on the spectral index is estimated at 0.2 (Madhavan 2013).

VERITAS has been monitoring Mrk 421 regularly for approximately 20 hr every year, as part of several long-term MWL monitoring campaigns (e.g., Acciari et al. 2011; Aleksić et al. 2015b). The general strategy is to take a 30 minute exposure on every third night when the source is visible, with coordinated, simultaneous X-ray observations (usually with the X-Ray Telescope (XRT) on board the Swift satellite). In contrast, the three long and simultaneous observations with XMM-Newton and VERITAS on 2014 April 29, May 1, and May 3 are specific attempts to catch rapid flares on top of elevated flux states simultaneously in the X-ray and TeV bands.

Due to high atmospheric dust conditions at the VERITAS site on May 1, only data from the ToO observations on April 29 and May 3 have been used. The VERITAS observations on these two nights were taken in "wobble" mode (Fomin et al. 1994) with the source offset 05 from the center of the field-of-view. The zenith angles of the observations were between 10° and 40°. After deadtime correction, the total exposure time from these observations is 6.14 hr. The data were analyzed using the data analysis procedures described in Cogan (2008). Standard gamma-ray selection cuts, previously optimized for sources with a power-law spectrum of a photon index 2.5, were applied to reject cosmic-ray (CR) background events. The reflected-region background model (Berge et al. 2007) was used to estimate the number of CR background events that passed the cuts, and a generalized method from Li & Ma (1983) was used for the calculation of statistical significance. The VERITAS results are shown in Table 2.

To parameterize the curvature in the VERITAS-measured TeV spectra, a power-law model with an exponential cutoff has been used to fit the daily spectra:

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dE}}=K{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{-\alpha }{e}^{-\tfrac{E}{{E}_{\mathrm{cutoff}}}}.\end{eqnarray} \tag{ 1 }$$

However, we used a power-law model in the hysteresis study in Section 3.3, as it adequately describes each 10 minute integrated spectrum without the cutoff energy as an extra degree of freedom.

The Major Atmospheric Gamma-ray Imaging Cherenkov (MAGIC) telescope system consists of two 17 m telescopes, located at the Observatory Roque de los Muchachos, on the Canary Island of La Palma (28.8 N, 17.8 W, 2200 m a.s.l.). Stereoscopic observations provide a sensitivity of detecting a point source at ∼0.7% C. U. above 220 GeV in 50 hr of observation, and allow measurement of photons in the energy range from 50 GeV to above 50 TeV. The night-to-night systematic uncertainty in the VHE flux measurement by MAGIC is estimated to be of the order of 11% (Aleksić et al. 2016).

Mrk 421 was observed by MAGIC for six nights from 2014 April 28 to May 4, as part of a longer MWL observational campaign. The source was observed in "wobble" mode, with 04 offset with respect to the nominal source position (Fomin et al. 1994). After discarding data observed in poor weather conditions, the total analyzed data amounted to 3.3 hr of observations, with exposures per observation ranging from 14 to 38 minutes, and zenith angles spanning from 9° to 42°.

The MAGIC data were analyzed using the standard MAGIC analysis and reconstruction software (Zanin et al. 2013). The integral flux was computed above 560 GeV, the same as the energy threshold found in the VERITAS long-term light curve, in order to use all the observations including those at large zenith angles. The source gamma-ray flux varied between 1.3 and 2.2 C. U. above 560 GeV for different days in this period, with no significant intra-night variability. This flux value is 3–5 times larger than the typical VHE flux of Mrk 421 (Acciari et al. 2014; Aleksić et al. 2015a). These observations are not simultaneous with the XMM-Newton observations. The source is known to change spectral index with flux level (Krennrich et al. 2002), and hence we computed the photon flux above 560 GeV using the measured spectral shape above 400 GeV, which ranged from 2.8 to 3.3.

Fermi-LAT is a pair-conversion high-energy (HE) gamma-ray telescope covering an energy range from about 20 MeV to more than 300 GeV (Atwood et al. 2009). It has a large field-of-view of 2.4 sr that covers the full sky every 3 hr in the nominal survey mode. Thus, Fermi-LAT provides long-term sampling of the entire sky. However, it has a small effective area of ∼8000 cm² for >1 GeV, which is usually not sufficient to resolve variability on timescales of hours or less.

We analyzed the Fermi-LAT Pass 8 data in the week of the VERITAS observations and produced daily averaged spectra and a daily binned light curve. We selected events of class source and type front+back with an energy between 0.1 and 300 GeV in a 10° region of interest (RoI) centered at the location of Mrk 421, and removed events with a zenith angle >90°. The data were processed using the publicly available Fermi-LAT science tools (v10r0p5) with instrument response functions (P8R2_SOURCE_V6). A model with the contributions of all sources within the RoI with a test statistic value greater than 3, a list of 3FGL sources within a source region of a radius of 20° from Mrk 421, and the contribution of the Galactic (using file gll_iem_v06.fit) and isotropic (using file iso_P8R2_SOURCE_V6_v06.txt) diffuse emission was used. This model was fitted to LAT Pass 8 data between 2014 April 1 and June 1 using an unbinned likelihood analysis (gtlike). The test-statistic maps were examined to ensure no unmodeled transient sources were present in the RoI during the period analyzed. All other best-fit parameters in the model were then fixed, except the spectral normalization and the power-law index of Mrk 421, in order to perform a spectral and temporal maximum likelihood analysis.

The XMM-Newton satellite carries the European Photon Imaging Camera (EPIC) pn X-ray CCD camera (Strüder et al. 2001), including two metal-oxide-silicon (MOS) cameras and a pn camera. The reflection grating spectrometers (RGS) with high-energy resolution are installed in front of the MOS detector. The incoming X-ray flux is divided into two portions for the MOS and RGS detectors. The EPIC-pn (EPN) detector receives the unobstructed beam and is capable of observing with very high time resolution. The Optical/ultraviolet (UV) Monitor (OM) onboard the XMM-Newton satellite provides the capability to cover a 17' × 17' square region between 170 and 650 nm (Mason et al. 2001). The OM is equipped with six broad-band filters (U, B, V, UVW1, UVM2, and UVW2).

Three ToO observations were taken simultaneously with the VERITAS observations on 2014 April 29, May 1, and May 3. To fully utilize the high time resolution capability of XMM-Newton in both the X-ray and optical/UV bands, all three ToO observations of Mrk 421 were taken in EPN timing mode and OM fast mode. MOS and RGS were also operated during the observations, but the data were not used due to the relatively low timing resolution and the lack of X-ray spectral lines from the source. The EPN camera covers a spectral range of approximately 0.5–10 keV and, with the UVM2 filter, the OM covers the range of about 200–270 nm.

XMM-Newton EPN and OM data were analyzed using Statistical Analysis System (SAS) software version 13.5 (Gabriel et al. 2004). An X-ray loading correction and a rate-dependent pulse height amplitude (RDPHA) correction were performed using the SAS tool epchain. We ran the SAS task epproc to produce the RDPHA results, which applies calibrations using known spectral lines and is likely more accurate than the alternative charge transfer inefficiency corrections.⁶⁵ Note that even after the RDPHA corrections, residual absorption features (not associated with the source) can still be present in the spectrum (see e.g., Pintore et al. 2014). To account for the source and the residual spectral features, the X-ray spectra were fitted using XSPEC version 12.8.1 with a model including a power law, a photoelectric absorption component representing the Galactic neutral hydrogen absorption, an absorption edge component, and two Gaussian components. The last three components are only associated with the instrument. They account for the oxygen K line at ∼0.54 keV, the silicon K line at ∼1.84 keV, and the gold M line at ∼2.2 keV, respectively. The model can be expressed as follows:

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dE}}=\left\{\begin{array}{l}{e}^{-{n}_{{\rm{H}}}\sigma (E)}\left[{K}_{\mathrm{PL}}{E}^{-\alpha }+\displaystyle \sum _{i}\tfrac{{K}_{G,i}}{\sqrt{2\pi }{\sigma }_{i}}{e}^{-\tfrac{{(E-{E}_{0,i})}^{2}}{2{\sigma }_{i}^{2}}}\right],\\ \qquad E\leqslant {E}_{c};\\ {e}^{-D{(E/{E}_{c})}^{-3}-{n}_{{\rm{H}}}\sigma (E)}\left[{K}_{\mathrm{PL}}{E}^{-\alpha }+\displaystyle \sum _{i}\tfrac{{K}_{G,i}}{\sqrt{2\pi }{\sigma }_{i}}{e}^{-\tfrac{{(E-{E}_{0,i})}^{2}}{2{\sigma }_{i}^{2}}}\right],\\ \qquad E\geqslant {E}_{c};\end{array}\right.\end{eqnarray} \tag{ 2 }$$

where n_H is the column density of neutral hydrogen; K_PL and ${K}_{G,i}$ are the normalization factor for the power-law component and the ith Gaussian component, respectively; E_c, ${E}_{0,i}$, and σ_i are the threshold energy of the absorption edge, the center and the standard deviation of the ith Gaussian component, respectively; D is the absorption depth at the threshold energy E_c; and α is the photon index of the power-law component in the model. The edge component at ∼0.5 keV is replaced by a Gaussian component (the 0th component in Table 4) for data taken on May 3 since the latter provides a better fit. We fix the column density of Galactic neutral hydrogen to N_H ≈ 1.9 × 10²⁰ cm⁻², which was measured by the Leiden/Argentine/Bonn (LAB) survey toward the direction of Mrk 421 (Kalberla et al. 2005). It is worth noting that the best-fit power-law index hardly changes when we set N_H free.

The count rate measured by the EPN camera with a thin filter can be converted to flux using energy conversion factors (ECFs, in units of 10¹¹ cts cm² erg⁻¹), which depend on the filter, the photon index α, the Galactic n_H absorption, and the energy range (Mateos et al. 2009). The flux f, in units of erg cm⁻² s⁻¹, can be obtained by f = rate/ECF, where rate has the units of cts s⁻¹. A similar flux conversion factor is used for the OM UVM2 filter to convert each count at 2310 Å to flux density 2.20 × 10¹⁵ erg cm⁻² s⁻¹ Å⁻¹. A 2% systematic uncertainty error was added to the OM light curve.

The long-term Swift-XRT light curve is produced using an online analysis tool The Swift-XRT data products generator⁶⁶ (Evans et al. 2009). This tool is publicly available and can be used to produce Swift-XRT spectra, light curves, and images for a point source. A light curve of Mrk 421 was made from all Swift-XRT observations available from 2005 March 1 to 2014 April 30, integrated between 0.3 and 10 keV, with a fixed bin width of 50 s. We cut the first 150 s of each WT observation, during which it is possible that the satellite could still be settling, thus causing non-source-related deflections in the light curve.

Regular optical observations of a sample of gamma-ray-bright blazars, including Mrk 421, have been carried out at Steward Observatory since the launch of the Fermi satellite (Smith et al. 2009), and these data are publicly accessible.⁶⁷ For the 2014 April–May MWL observing campaign, the SPOL optical, dual-beam spectropolarimeter (Schmidt et al. 1992b) was used at the Steward Observatory 1.54 m Kuiper Telescope on Mt. Bigelow, Arizona from April 25 to May 4 UTC. When the weather permitted, the usual observing frequency of one observation per night for Mrk 421 was increased to four per night after April 26 so that any rapid changes in linear polarization and optical flux could be better tracked. The spectropolarimeter was configured with a 600 l/mm diffraction grating that yields a dispersion of 4 Å pixel⁻¹, spectral coverage from 4000 to 7550 Å, and resolution of ∼16 Å. The CCD detector is a thinned, anti-reflection coated 1200 × 800 STA device with a quantum efficiency of about 0.9 from 5000 to 7000 Å. All polarization observations of Mrk 421 were made with a 3'' × 50'' slit oriented so that its long (spatial) dimension is east–west on the sky and the CCD was binned by two pixels (∼09) in the spatial direction. An observation of Mrk 421 typically consists of a 30 s exposure at all 16 positions of the λ/2-wave plate, properly sorted into four images with each image containing the two orthogonal polarized beams created by a Wollaston prism in the optical path. Extraction of the sky-subtracted spectra of Mrk 421 is done using a 3'' × 9'' aperture for all polarization observations to keep the contribution of the unpolarized starlight from its host galaxy as constant as possible for all measurements. Medians in the wavelength range of 5000–7000 Å are taken of the resulting spectra of the linear Stokes parameters q and u and used to calculate the observed degree of polarization (P) and the position angle of the polarization on the sky (χ). The instrumental polarization of SPOL has been consistently measured to be ≪0.1% and is ignored. Likewise, Galactic interstellar polarization is negligible in the direction of Mrk 421 based on the amount of reddening estimated in this line of sight (Av ∼ 0.042; Schlafly & Finkbeiner 2011). Given the lack of significant instrumental and Galactic interstellar polarization, the uncertainties in the measurements made by SPOL are dominated by photon statistics and typically σ_p < 0.1% when the spectropolarimetry is binned by 2000 Å. The polarization position angle was calibrated during the campaign by observing the polarization standard stars Hiltner 960 and VI Cyg #12 (Schmidt et al. 1992a).

Optical flux monitoring of Mrk 421 during this period was accomplished using SPOL with a 76 × 50'' slit when conditions were clear. As with the spectropolarimetry, the slit is oriented east–west on the sky and although the larger slit admits more host galaxy starlight, it minimizes slit losses as a function of wavelength. Differential photometry with "Star 1" (Villata et al. 1998) was used to calibrate the V-band magnitude of Mrk 421 within a spectral extraction aperture of 76 × 9''. Generally, a single 30 s exposure at a set wave plate position is obtained for both the blazar and the comparison star. A standard Johnson V filter bandpass transmission curve is multiplied to the extracted spectra and the instrumental fluxes for the objects are compared to derive the brightness of Mrk 421. This measurement is typically performed twice per visit to Mrk 421 to check the consistency of the photometry. The dominant source of uncertainty for the flux measurements is the precision of the V-band calibration of the comparison star (0.02 mag). The flux contribution of the host galaxy in the R-band for a rectangular aperture of 76 × 9'' centered at the Mrk 421 was estimated using the measurements in Nilsson et al. (2007), and converted to V-band using an E galaxy template of age 11 Gyr at a redshift of 0.031 that gives a V − R of 0.686 (Fukugita et al. 1995).

Contemporaneous observations of Mrk 421 were taken with the Owens Valley Radio Observatory (OVRO) at 15 GHz (Richards et al. 2011) and the Combined Array for Research in Millimeter-Wave Astronomy (CARMA) at 95 GHz (Bock et al. 2006).

The OVRO 40 m telescope is equipped with a cryogenic, low-noise high electron mobility transistor amplifier with a 15.0 GHz center frequency and 3 GHz bandwidth. The two off-axis sky beams are Dicke-switched with the source alternating between the two beams, in order to remove atmospheric and ground contamination. The receiver gain is calibrated using a temperature-stable diode noise source. The systematic uncertainty in the flux density scale is estimated to be approximately 5%, and is not included in the error bars. More details of the reduction and calibration procedure can be found in Richards et al. (2011).

The CARMA observations of Mrk 421 were made using the eight 3.5 m telescopes of the array with a central frequency of 95 GHz and a bandwidth of 7.5 GHz. The amplitude and phase gain were self-calibrated on Mrk 421. The absolute flux was calibrated from a temporally nearby observation of the planets Mars, Neptune, or Uranus, or the quasar 3C 273. The absolute systematic uncertainty is estimated to be approximately 10%, and is not included in the error bars.

Figures 1–3 show simultaneous light curves in VHE, X-ray, and UV bands. The VERITAS light curves are binned in 10 minute intervals, and shown with the integrated flux above the highest-energy threshold among all observations taken on the corresponding night: 560 GeV on April 29 (315 GeV for the first ∼3.5 hr) and 225 GeV on May 3. The higher-energy threshold on April 29 is a result of the larger zenith angle at which the source was observed (∼48° during the last 30 minute exposure of the night), which consequently leads to more distant shower maxima from the telescope. The X-ray light curves in the middle panels show the count rates measured by XMM-EPN between 0.5 and 10 keV, binned in 50 s intervals. The bottom panels show UV light curves constructed from the XMM-OM count rate using the UVM2 filter in both Image and Fast modes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The average VERITAS integral fluxes above 0.4 TeV are (1.27 ± 0.03) × 10⁻⁶ photons m⁻² s⁻¹ on April 29 and (1.10 ± 0.04) × 10⁻⁶ photons m⁻² s⁻¹ on May 3. As shown in Table 3, a constant fit to the X-ray light curves yields large reduced χ² values (corresponding p-values <1 × 10⁻⁵, thus rejecting the hypothesis of constant flux), implying the presence of intra-night variability. The corresponding p-values in the VHE band, of 0.003 and 0.07, imply marginally significant intra-night variabilities in the VERITAS light curves >315 GeV on April 29 and >225 GeV on May 3.

Table 3.  Reduced χ² Values for a Constant Fit to the Light Curves and the Corresponding p-values for VERITAS Light Curves

Date	VERITAS		XMM-EPN	XMM-OM
	${\chi }_{\mathrm{red}}^{2}$	p-value	${\chi }_{\mathrm{red}}^{2}$	${\chi }_{\mathrm{red}}^{2}$
2014 Apr 29	2.1 (>315 GeV)	0.003	11.1	0.9
	1.2 (>560 GeV)	0.2
2014 May 01	...	...	48.0	0.9
2014 May 03	1.6	0.07	7.0	0.9

Download table as:  ASCIITypeset image

Another quantity that describes the relative amount of variability is the fractional variability amplitude. Following the descriptions in Vaughan et al. (2003) and Poutanen et al. (2008), the fractional variability F_var and its error ${\sigma }_{{F}_{\mathrm{var}}}$ are calculated as

$$\begin{eqnarray*}&&{F}_{\mathrm{var}}=\sqrt{\displaystyle \frac{{S}^{2}-\langle {\sigma }_{\mathrm{err}}^{2}\rangle }{\langle F{\rangle }^{2}}}\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&{\sigma }_{{F}_{\mathrm{var}}}=\sqrt{{F}_{\mathrm{var}}^{2}+\sqrt{\displaystyle \frac{2\langle {\sigma }_{\mathrm{err}}^{2}{\rangle }^{2}}{N\langle F{\rangle }^{4}}+\displaystyle \frac{4\langle {\sigma }_{\mathrm{err}}^{2}\rangle {F}_{\mathrm{var}}^{2}}{N\langle F{\rangle }^{2}}}}-{F}_{\mathrm{var}},\end{eqnarray*}$$

where S is the standard deviation of the N flux measurements, $\langle {\sigma }_{\mathrm{err}}^{2}\rangle$ is the mean squared error of these flux measurements, and $\langle F\rangle$ is the mean flux.

The fractional variability results from the simultaneous XMM-Newton and VERITAS data, as well as from contemporaneous MWL data, are shown in Figure 4. The VHE fractional variability has been computed using the flux light curves integrated above 315 GeV for data from the first 3.5 hr on April 29, and above 225 GeV on May 3. The VERITAS fractional variability is ∼13% ± 3% on April 29, and ∼8% ± 5% on May 3.

Zoom In Zoom Out Reset image size

The fractional variability of the X-ray flux is low, but significantly above zero in all three energy intervals 0.5–1 keV, 1–3 keV, and 3–10 keV. Comparing these three X-ray energy intervals, a higher fractional variability is observed at higher frequencies on April 29 (navy open diamonds) and on May 3 (cyan open diamonds), in agreement with previous results (e.g., Błażejowski et al. 2005). This may be explained as the manifestation of a different synchrotron-cooling time at different energies, t_cool ∝ E^−1/2. In the slow-cooling regime (the cooling time t_cool is longer than the dynamic timescale R/c), the cooling time is shorter for higher-energy particles, leading to a faster variability in radiation at higher energies. Therefore, more variation at higher energies is observed compared to lower energies on the same timescales, which directly leads to a higher fractional variability for higher-energy emissions. However, the same trend is less obvious on May 1 (blue open diamonds), when only X-ray data are available.

Contemporaneous MWL observations often provide valuable information about the activity of the source, e.g., abrupt changes in the radio and optical polarizations would reveal the emergence of a compact region that may be connected to a flaring event (see e.g., Arlen et al. 2013). MWL light curves of Mrk 421 from MJD 56775 to 56781are shown in Figure 5. The TeV light curve, measured by MAGIC above 560 GeV on each night before the XMM-Newton observations, is shown in blue in the top panel. The fractional variabilities of the VHE flux for a one-week duration with a 30 minute bin width are measured by both VERITAS and MAGIC are shown as the gray open square and the gray cross in Figure 4, respectively. Only statistical uncertainties are taken into account in the calculation of the VHE fractional variability. As mentioned in Section 2.1, the night-to-night systematic uncertainty in the VHE flux measurements from MAGIC in estimated to be ∼11%, the primary contribution to which is the fluctuation of the atmospheric transmission (see Aleksić et al. 2016 for further details). We follow a similar approach and estimate the systematic uncertainty in the VHE flux measured by VERITAS with 10 minute intervals to be less than ∼10% using observations of the Crab Nebula under similar conditions as the observations of Mrk 421 in this work. Adding a 10% systematic uncertainty in quadrature to the statistical uncertainty of the VERITAS-measured VHE flux reduces its one-week F_var value from 24.2% ± 2.5% to 21.9% ± 3.6%, and the F_var value on April 29 from ∼13% ± 3% to ∼8% ± 5%, and the F_var value ∼8% ± 5% on May 3 should be considered as an upper limit.

Zoom In Zoom Out Reset image size

The daily Fermi-LAT light curve (shown in the second panel) does not suggest any significant variability, although there might be a slight drop in GeV flux after May 1. X-ray count rates from the Swift-XRT are also shown, together with those measured by XMM-EPN in the third panel. The Swift-XRT results fill the gaps between the three XMM-Newton observations, and also show significant variability, illustrated by the fractional variability of ∼20% computed from XRT data between April 28 and May 4, shown as the gray diamond in Figure 4. Note that fractional variabilities of ∼20% to ∼40% in the X-ray and ∼30% in VHE during typical non-flaring states were found by Aleksić et al. (2015a), and much higher values, of ∼40% to >60% in both X-ray and VHE during flaring states were found by Błażejowski et al. (2005) and J. Aleksić et al. (2016, in preparation). It is worth noting that fractional variability depends on the bin width, the sampling frequency, and the duration of the light curve (see the discussion section), which makes it more difficult to compare the F_var values measured using different light curves.

The optical photometric variability of Mrk 421 from April 25 to May 4 is mild (see the fourth panel of Figure 5). Mrk 421 varied from 12.8 to 12.9 in V band, with the maximum optical flux observed on MJD 56775. Over a 24 hr period, the object showed a maximum ΔV ∼ 0.1 mag and intra-night variability was generally <0.05 mag. During this period, Mrk 421 was close to the middle of the range of optical brightness it has shown since 2008 (V ∼ 11.9–13.6). The source does not exhibit strong variability in the UV band, nor at 15 GHz or 95 GHz (see the fourth and fifth panels of Figure 5).

In contrast to the flux variations, the optical polarization of Mrk 421 showed more pronounced variability during the dates shown. The observed polarization fraction P peaked at MJD 56779 (4%–5%) with minima of P ∼ 2% two days preceding and one day after the polarization maximum (see the sixth panel of Figure 5). The polarization peak reaches only about half of the highest polarization levels observed for this object (10%–13%) since 2008. From 2008 to 2015, the Steward Observatory blazar data archive identified several periods when the polarization of Mrk 421 is <1%. Like many blazars, the full range of values for the polarization position angle is exhibited by Mrk 421 over a timescale of several years. During the dates shown, the polarization angle χ varied between 20° and 55° with rotations of nearly 20° observed in a day (see the bottom panel of Figure 5). Intra-night variations as large as ∼10° were observed on timescales as short as two hours. During the epoch of the campaign, χ is roughly orthogonal to the position angle of the 43 GHz VLBI jet (∼−35°),⁶⁸ implying that the magnetic field within the region emitting the polarized optical continuum is more-or-less aligned with the jet.

We plot the VHE flux against X-ray count rate in Figure 6. There is no significant evidence for correlation. The Pearson correlation coefficient between X-ray flux and TeV gamma-ray flux >560 GeV is 0.48 (the 90% confidence interval is 0.24–0.67), and that between X-ray flux and TeV gamma-ray flux >315 GeV is 0.60 (the 90% confidence interval is 0.35–0.76). The values of the correlation coefficients only suggest a moderate positive correlation, without considering the uncertainties on the measurements. Therefore, we focus on the flux correlations between three X-ray bands, and between two TeV gamma-ray bands, respectively.

Zoom In Zoom Out Reset image size

3.2.1. Hard/Soft X-Ray Correlation

We further divide XMM-EPN X-ray light curves into three energy bands, 0.5–1 keV, 1–3 keV, and 3–10 keV (as shown in the left panels in Figures 7–9). Z-transformed discrete correlation functions (ZDCFs) between these soft- and hard-X-ray light curves are calculated using a publicly available code, ZDCF v2.2 developed by Alexander (2013), as shown in the right panels in Figures 7, 8, and 9. At least 11 pairs of light curve points in each time delay bin are required to calculate the ZDCF, zero lag is not omitted, and 1000 Monte Carlo runs were used to estimate the measurement error, in addition to the error calculated in the z-space. From the ZDCFs, the corresponding time lags are calculated using PLIKE v4.0 also developed by Alexander (2013). No significant evidence is found for any leads or lags between the soft-X-ray band relative to the hard-X-ray band.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2.2. Gamma-ray Intraband Correlation

The cross-correlations between light curves of blazars at TeV energies are particularly interesting, not only because they can provide insight to the particle acceleration and radiation, but also due to their potential to test Lorentz-invariance violation, which is a manifestation of an energy-dependent speed of light at the Planck scale predicted by foamy structures of spacetime in certain quantum theories (e.g., Aharonian et al. 2008; MAGIC Collaboration et al. 2008; Zitzer for the VERITAS Collaboration 2013).

In a similar fashion as done for the X-ray data (as described in Section 3.2.1), we divide the gamma-ray light curves into two bands, and compute ZDCFs and time lags as shown in Figure 10. The chosen bands are 315–560 GeV and 560 GeV–30 TeV on April 29, and 225–560 GeV and 560 GeV–30 TeV on May 3, so that the event rates are comparable in the higher- and lower-energy bins. ZDCFs are calculated using light curves binned by 10 minutes and 4 minutes, respectively. The 1-σ confidence interval of the time lag of maximum likelihood is calculated between −2000 s and 2000 s using plike_v4.0.

Zoom In Zoom Out Reset image size

To understand the ZDCFs produced by random noise, we simulate flicker-noise (whose power spectral density (PSD) distribution is proportional to 1/f) and Gaussian white noise with 10 minute bin width and similar duration as the data. The 95% confidence regions calculated from ZDCF values between 200 pairs of simulated light curves are plotted in Figure 10, along with ZDCFs calculated from 10 minute binned VERITAS light curves above and below 560 GeV. No evidence for time leads or lags is present in the gamma-ray data; however, given the lack of a strong detection of variability, this lack of evidence for any leads/lags is not surprising since the sensitivity to such leads and lags is dependent on the amplitude of the detected variability in these data.

Besides the possibility of time lags at different energies, the spectral evolution during blazar flares is also informative. A general trend, that the spectrum is harder when the flux is higher, has been observed in blazars in both the X-ray and gamma-ray bands (e.g., Fossati et al. 2008; Acciari et al. 2011). Several possibilities can lead to such a trend if the spectrum is being measured at the HE end of the synchrotron and inverse-Compton SED peaks. For example, this could result from an increase in the maximum electron energy, or a hardening in the electron energy distribution. If, however, the X-ray and gamma-ray observations are sampling the emission near the two peaks of the SED, this "harder-when-brighter" effect could also be the result of an increase of the SED peak frequency, which could arise from an increase in magnetic field strength or Doppler factor.

As well as the "harder-when-brighter" trend in the hardness-flux relation, competition between the acceleration and cooling timescales can lead to spectral hysteresis, since the spectral hardness differs on the rising and falling edges of flares (e.g., Kirk et al. 1998; Li & Kusunose 2000; Sato et al. 2008). If one plots the hardness–flux relation so that the spectrum is harder in the positive y-axis direction, and the flux is higher in the positive x-axis direction, the spectral hysteresis can be seen as a "loop" pattern. Since the spectral hysteresis is driven by the same timescales that determine the time lags at different energies, the direction of the hysteresis loop should be consistent with the sign of the time lag. Specifically, a "hard lag" should correspond to counter-clockwise hysteresis loops (see Section 3.2.1), while a "soft lag" will lead to clockwise hysteresis loops. Therefore, the hardness–flux plot offers an alternative view of the timescales in the system, and can be compared with the time-lag studies presented in the previous section.

The nightly XMM-Newton-EPN X-ray spectra are fitted using a power-law model with absorption from neutral hydrogen, as well as instrumental features from oxygen, silicon, and gold, following the description in Section 2.3. The results are shown in Figure 11 and Table 4. Note that the equivalent width EW_i of each Gaussian component is shown instead of the standard deviation σ_i, because σ_i is small compared to the energy-bin size of the X-ray spectrum and therefore not well-constrained by the fit.

Zoom In Zoom Out Reset image size

Table 4.  XMM-Newton-EPN Spectral Fit Results and ECFs using the Absorbed-power-law Model Plus Three Instrumental Features Described in Equation (2)

Parameter	Unit	Value on 56776	Value on 56778	Value on 56780
Ec	keV	0.540 ± 0.003	0.530 ± 0.003	...
D	⋯	0.118 ± 0.005	0.146 ± 0.009	...
${E}_{\mathrm{0,0}}$	keV	...	...	0.48 ± 0.04
EW0	keV	...	...	<0.11
${K}_{G,0}$	10−3	...	...	0.022 ± 0.015
nH	1020 cm−2	1.9	1.9	1.9
α	⋯	2.649 ± 0.002	2.817 ± 0.004	2.450 ± 0.002
KPL	⋯	0.2714 ± 0.0004	0.1711 ± 0.0004	0.2251 ± 0.0002
${E}_{\mathrm{0,1}}$	keV	1.88 ± 0.02	1.88 ± 0.02	1.87 ± 0.03
EW1	eV	4 ± 2	4 ± 3	2 ± 2
${K}_{G,1}$	10−4	2.2 ± 0.4	1.1 ± 0.4	0.9 ± 0.5
${E}_{\mathrm{0,2}}$	keV	2.26 ± 0.01	2.25 ± 0.01	2.26 ± 0.01
EW2	eV	26 ± 3	19 ± 3	19 ± 7
${K}_{G,2}$	10−4	8.1 ± 0.4	3.3 ± 0.4	5.9 ± 0.6
Reduced χ2	⋯	1.6	1.3	2.3

Note. Note that the fits are insensitive to the width σ_i of the Gaussian component, therefore the equivalent width (EW) values are shown instead. The instrumental feature around 0.5 keV is fitted as an edge component on MJD 56776 and 56778, but as a Gaussian component on 56780.

Download table as:  ASCIITypeset image

We then divide each of the XMM-Newton and VERITAS observations of Mrk 421 on the two nights with good atmospheric conditions into simultaneous 10 minute intervals, and perform spectral fitting for each interval, using an absorbed-power-law model for the X-ray data, and a power-law model for the gamma-ray data. We note that although the X-ray spectral fit is significantly improved by introducing the extra instrumental features, the spectral indices α remain unchanged within the uncertainty. Therefore we are confident that the hardness ratios and the spectral indices derived for each 10 minute interval are robust, and are not severely affected by the instrumental features.

We identify several time intervals with a rise and a subsequent fall of flux in the X-ray light curves, and plot photon index and hardness ratio against flux (or count rate) for these bumps (see Figures 12 and 13). Black arrows indicate the order of time for each point. Measurements taken at different times are also color coded to guide the eye. A "harder-when-brighter" effect can be identified on some individual X-ray branches, (e.g., the blue and green points in the top right panel in Figure 12). The observed "soft lag" on May 3 indicates a harder spectrum when flux rises, and a softer spectrum when flux falls, corresponding to a clockwise loop (in orange) in the top right panel of the spectral hysteresis plot in Figure 13. Similarly, for the "hard lag" scenario on April 29, a counter-clockwise loop is predicted and observed, as shown in Figure 12. It is interesting to note that the time lag and loop direction changes for flares separated by a few days, in spite of similar flux levels.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The same analysis has been carried out for VHE data, and similar plots are shown. The uncertainties in VHE flux, hardness ratio, index, and hence hysteresis direction, are large, and do not allow us to draw any firm conclusions. We note that although the VHE spectrum of Mrk 421 is likely curved, a power-law model describes the data reasonably well for each 10 minute interval (the average reduced χ² value is ∼1.05).

The hardness ratio pattern offers a more crude but less model-dependent estimation of the same signature. However, the hardness–flux diagram of the VERITAS observations also has large uncertainties. At the flux level of roughly 1–2 C. U., such spectral hysteresis studies with the current generation of ground-based gamma-ray instruments is difficult. This offers a reference for the future in defining flux level trigger criteria for ToO observations, aiming for similar goals.

The PSD can be estimated from a periodogram, which is defined as the squared modulus of the Fourier transform of a time series

$$\begin{eqnarray*}P(\nu )=| {F}_{N}(\nu ){| }^{2} & = & {\left[\displaystyle \sum _{i=1}^{N}f({t}_{i})\cos (2\pi \nu {t}_{i})\right]}^{2}\\ & & +{\left[\displaystyle \sum _{i=1}^{N}f({t}_{i})\sin (2\pi \nu {t}_{i})\right]}^{2}.\end{eqnarray*}$$

We apply Leahy normalization to the periodogram to get the PSD, so that Poisson noise produces a constant power at an amplitude of 2.

The top panels of Figure 14 show PSDs calculated from the X-ray light curves of Mrk 421, measured by XMM-Newton-EPN on 2014 April 29, May 1, and May 3, respectively. The light curves are first binned by 50 s intervals, then divided into two equal-length segments, each of 128 bins. A raw power spectrum is calculated for each segment, and averaged over both segments. Then the power spectrum is rebinned geometrically with a step factor of 1.2 (i.e., a bin edge in frequency is the previous bin edge multiplied by a factor of 1.2). To estimate the error of the PSD in each bin, the standard deviation from all PSD points in the bin is used when it contains more than five points, otherwise the theoretical standard deviation of an exponential distribution is used. The PSDs cover a frequency range of 4 × 10⁻⁴ to 10⁻² Hz. At higher frequency, the shape of the PSDs becomes flatter due to Poisson noise, which is shown as the flat line with a power of two under the Leahy normalization. However, we note that the PSD is well above the Poisson noise level up to ∼10⁻³ Hz on all three days, which corresponds to timescales of under an hour. On May 1, the variability is still present, reaching ∼2–3 × 10⁻³ Hz, which is shorter than 10 minutes.

Zoom In Zoom Out Reset image size

We then simulate 1000 light curves at each of the indices ranging from 0.5 to 2.5 in steps of 0.1 following Timmer & Koenig (1995), and compare the data and the simulations to compute success fractions (SuFs) as an estimation of the power-law index α following the method described by Uttley et al. (2002). The results are plotted in the bottom panels of Figure 14. The SuF peaks at the PSD indices of ∼1.5–1.8 on April 29 with a peak value of ∼0.5, ∼1.3–1.6 on May 1 with a peak value of ∼0.8, and ∼1.2–1.6 on May 3 with a peak value of ∼0.7. These values represent the indices of the PSD better than the simple power-law fit to the PSD, as SuFs are less biased by the spectral leakage which is present in both the data and the simulations. The range of PSD indices is consistent with previous studies of the same source at lower frequencies, e.g., a PSD index of 1.35–1.85 was found for frequency range ∼10⁻⁸–2 × 10⁻⁶ Hz (Isobe et al. 2015). However, a break at frequency of ∼9.5 × 10⁻⁶ Hz was found in the X-ray PSD of Mrk 421 by Kataoka et al. (2001), and the PSD index above this frequency was determined to be 2.14. Surprisingly, we do not find evidence of such a steep PSD at ≳4 × 10⁻⁴ Hz.

The SED of simultaneous VERITAS and XMM-Newton data, as well as contemporaneous MWL data, is shown in Figure 15. Daily averaged HE (∼30 MeV–100 GeV) gamma-ray spectra are constructed from Fermi-LAT data between 100 MeV and 300 GeV, and butterfly regions of 95% confidence level are shown. Note that the uncertainty is large because of the scarcity of HE photons in the one-day window. Optical spectra from the Steward Observatory between 400 and 750 nm on May 3, radio data from CARMA at 95 GHz taken on both nights, and from OVRO at 15 GHz on other nights within the week are also shown.

Zoom In Zoom Out Reset image size

The X-ray and VHE spectra are located at the falling slopes of the lower- and higher-energy spectral peaks, respectively. The synchrotron peak is between the UV measurement at ∼10¹⁵ Hz and the soft end of the X-ray spectrum at ∼10¹⁷ Hz. Although the Fermi-LAT spectrum is not very constraining, the HE spectral peak appears to be just below ∼100 GeV, as suggested by the TeV spectrum.

We use a static SSC model described in Krawczynski et al. (2002) to study the observed SEDs. The set of parameters used to generate the solid blue and red curves as shown in Figure 15 are listed in Table 5. The static one-zone SSC model is roughly consistent with the data. The synchrotron peak frequency given by the model is ν_syn ∼ 4 × 10¹⁶ Hz, while the inverse-Compton peak lies at ν_SSC ∼ 5 × 10²⁵ Hz. Using the peak frequencies and the observed spectral indices below and above the synchrotron peak, we can follow Equation (16) given by Tavecchio et al. (1998) assuming the system is in the Klein–Nishina (KN) regime and estimate the strength of the magnetic field to be

$$\begin{eqnarray*}B & \approx & 2.7\times {10}^{-7}\delta \displaystyle \frac{{\nu }_{\mathrm{syn}}}{{\nu }_{\mathrm{SSC}}^{2}}{\left(\displaystyle \frac{{{mc}}^{2}}{h}\right)}^{2}\\ & & \times {\left\{{\rm{\exp }}\left[\displaystyle \frac{1}{{\alpha }_{1}-1}+\displaystyle \frac{1}{2({\alpha }_{2}-{\alpha }_{1})}\right]\right\}}^{2}\displaystyle \frac{1}{1\,+\,z}\approx 0.06\,{\rm{G}},\end{eqnarray*}$$

where $\delta ={[{\rm{\Gamma }}(1-\beta \cos \theta )]}^{-1}\approx 20.3$ is the Doppler factor (Γ and θ taken from the SED modeling results on May 3), α₁ ∼ 0.5 and α₂ ∼ 1.5 are the observed spectral indices below and above the synchrotron peak. We also calculate the Doppler factor limits for the KN and Thomson regimes to be ∼13 and ∼29, respectively, following Equations (13) and (17) in Tavecchio et al. (1998). The Doppler factor value of 20.3 obtained from the SSC model is in between the two limits, therefore the gamma-ray emission is in the Thomson-to-KN transition regime.

Table 5.  Parameters Used for the SSC Model in Figure 15

Data set	Γ	θ	B	R	we	$\mathrm{log}{E}_{\min }$[eV]	$\mathrm{log}{E}_{\max }$[eV]	$\mathrm{log}{E}_{\mathrm{break}}$[eV]	p1	p2
		(deg)	(10−1 G)	(1014 m)	(erg m−3)
Apr 29	23.0	2.7	0.4	4.8	0.001	8.5	11.65	10.9	2.0	4.3
May 3	23.0	2.8	0.4	5.0	0.001	8.5	12.0	10.75	2.0	3.8
B05 low	19.48	5.01	4.05	0.7	0.13777	6.5	11.22	10.34	2.05	3.6
B05 med	20.05	3.2	1.02	1.0	0.03192	6.5	11.55	10.98	2.05	3.4
B05 high	20.0	3.905	2.6	0.7	0.086	6.5	11.6	11.0	2.05	3.4

Note. B05 refers to Błażejowski et al. (2005). Γ, θ, and R are the bulk Lorentz factor, the viewing angle, and the radius of the emitting region, respectively, B is the strength of the magnetic field, w_e is the energy density of the electrons, $\mathrm{log}{E}_{\min }$, $\mathrm{log}{E}_{\max }$, and $\mathrm{log}{E}_{\mathrm{break}}$ are the logarithms of the minimum, maximum, and break electron energy, respectively, and p₁ and p₂ are spectral indices of the electrons below and above the break energy, respectively.

Download table as:  ASCIITypeset image

From April 29 to May 3, the change in SED might be described by an increase in the radius of the emitting region R, along with an increase in the maximum energy E_max, a slight decrease in break energy E_break of the electron distribution, a harder electron spectrum (smaller p₂) after the break energy, and a slight increase in the Doppler factor (see Table 5). Note that a slight increase in magnetic field strength can have a similar effect as the increase in the Doppler factor. If these parameters indeed describe the evolution of the SED, it is consistent with the results of an expansion of the emitting region. The direct result of such an expansion is an increase in the dynamic timescale t_dyn = R/c. Moreover, this will lead to a higher maximum energy of the electrons E_max, since the maximum possible gyro-radius has increased. Also, the synchrotron cooling break, which occurs at the electron energy that satisfies t_syn = t_dyn, decreases since t_syn ∝ γ⁻¹, where γ is the Lorentz factor of the electron.

Under the simple SSC model, the hypothesis of an increase in the maximum energy E_max would be consistent with a counter-clockwise spectral hysteresis pattern on April 29 and a clockwise pattern on May 3. The former is predicted when the system is observed near the maximum particle energy (see e.g., Kirk et al. 1998), in which case the increasing acceleration timescale with higher particle energy results in the higher-energy photons lagging behind the lower-energy ones. Consequently, this leads to a softer spectrum when flux rises because lower-energy photons are emitted faster, and a counter-clockwise spectral hysteresis pattern. If the maximum electron energy E_max is lower on April 29, the observed X-ray frequency (fixed at 0.5–10 keV) is closer to the maximum frequency, therefore the change in flux propagates from low to high frequencies and a counter-clockwise spectral hysteresis pattern should be observed. On the other hand, on May 3, as E_max increases, the observed X-ray frequency covers a relatively lower portion of the entire particle spectrum, therefore the change in flux propagates from high to low energies and a "soft lag" may be observed.