MULTIWAVELENGTH STUDY OF QUIESCENT STATES OF Mrk 421 WITH UNPRECEDENTED HARD X-RAY COVERAGE PROVIDED BY NuSTAR IN 2013

The Owens Valley Radio Observatory (OVRO) 40 m telescope was used for observation at 15 GHz, as a part of a long-term blazar monitoring program. Additional observations were scheduled at times of coordinated observations with X-ray and VHE γ-ray observatories. The data were reduced using standard processing and calibration techniques described in detail in Richards et al. (2011). Radio observations of Mrk 421 between 2.64 and 142 GHz have been obtained within the framework of the F-GAMMA program (Fuhrmann et al. 2007, 2014; Angelakis et al. 2010), a γ-ray blazar monitoring program related to the Fermi Gamma-ray Space Telescope. Observations with the Effelsberg 100 m and Pico Veleta 30 m telescopes are performed roughly once per month. The Effelsberg measurements are conducted with heterodyne receivers at 2.64, 4.85, 8.35, 10.45, 14.60, 23.05, and 32.0 GHz, while the Pico Veleta telescope is used with the EMIR receiver to provide the high-frequency (86.2 and 142.3 GHz) flux measurements. Standard data processing and calibration were performed as described in Angelakis et al. (2008, 2015). The Metsähovi Radio Observatory 14 m telescope also participated in this multi-instrument campaign, providing observations of Mrk 421 at 37 GHz every few days. Details of the observing strategy and data reduction for this monitoring program can be found in Teräsranta et al. (1998).

The coverage at optical frequencies was provided by various telescopes around the world within the GASP-WEBT program (e.g., Villata et al. 2008, 2009). In particular, the following observatories contributed to this campaign: Teide (IAC80), Crimean, Lowell (Perkins telescope), Roque de los Muchachos (KVA and Liverpool telescopes), Abastumani, Pulkovo, St. Petersburg, Belogradchik, Rozhen (50/70 cm, 60 cm, and 200 cm telescopes), Vidojevica, and Lulin. Additionally, many observations were performed with iTelescopes, Bradford Robotic Telescope, ROVOR, and the TUBITAK National Observatory. In this paper, we use only R-band photometry. The calibration stars reported in Villata et al. (1998) were used for calibration, and the Galactic extinction was corrected with the reddening corrections given in Schlafly & Finkbeiner (2011). The flux from the host galaxy was estimated using the R-band flux from Nilsson et al. (2007) for the apertures of 5'' and 75 used by various instruments. We applied an offset of −5 mJy to the fluxes from ROVOR in order to achieve better agreement with the light curves from the other instruments. This difference may be related to the specific spectral response of the filters used, or the different analysis procedures that were employed. Additionally, a point-wise fluctuation of 2% on the measured flux was added in quadrature to the statistical uncertainties in order to account for potential day-to-day differences in observations with any of the instruments.

Polarization measurements are utilized from four observatories: Lowell (Perkins telescope), St. Petersburg, Crimean, and Steward (Bok telescope). The polarization measurements from Lowell and St. Petersburg observatories are derived from R-band imaging polarimetry. The measurements from Steward Observatory are derived from spectropolarimetry between 4000 and 7550 Å with a resolution of ∼15 Å. The reported values are constructed from the median Q/I and U/I in the 5000–7000 Å band. The effective wavelength of this bandpass is not too different from the Kron-Cousins R band, and the wavelength dependence in the polarization of Mrk 421 seen in the spectropolarimetry during this period is not strong enough to significantly affect the variability analysis of the measurements from various telescopes. The observing and data-processing procedures for the polarization measurements are described in Larionov et al. (2008), Smith et al. (2009), and Jorstad et al. (2010).

Swift observations with the UV/Optical Telescope (UVOT; Roming et al. 2005) were performed only with the UV filters (namely, W1, M2, and W2). Observations with the optical filters were not needed because we had organized extensive coverage with ground-based optical telescopes, which have better sensitivity and angular resolution than Swift-UVOT. We performed aperture photometry for all filters in all the observations using the standard UVOT software distributed within the HEAsoft package (version 6.10) and the calibration included in the latest release of the CALDB. Counts were extracted from an aperture of 5'' radius for all filters and converted to fluxes using the standard zero points from Breeveld et al. (2011). The fluxes were then dereddened using the value of (Schlegel et al. 1998; Schlafly & Finkbeiner 2011) with ratios calculated using the mean Galactic interstellar extinction curve from Fitzpatrick (1999). No variability was detected within single exposures in any filter. The results of the processing were carefully verified, checking for possible contaminations from nearby objects within source apertures and from objects falling within background apertures. In almost all observations, Mrk 421 is on the "ghost wings" (Li et al. 2006) from the nearby star 51 UMa, so we estimated the background from two circular apertures of 16'' radius off the source but on the wings, excluding stray light and support structure shadows.

The complete list of Swift X-ray Telescope (XRT; Burrows et al. 2005) and UVOT observations used here is given in Table 1. The observations were organized to be taken simultaneously with (or as close as possible to) the MAGIC/VERITAS and NuSTAR observations, following the fruitful monitoring campaign practice since 2009. Swift observed the source 33 times in 2013 up to the end of March. All Swift-XRT observations were carried out using the Windowed Timing (WT) readout mode. The data set was processed with the XRTDAS software package (version 2.9.0) developed at ASDC and distributed with the HEASoft package (version 6.13). Event files were calibrated and cleaned with standard filtering criteria with the xrtpipeline task using the latest calibration files available in the Swift CALDB. The average spectrum was extracted from the summed cleaned event file. Events for the spectral analysis were selected within a circle of 20 pixel (≃46'') radius, which encloses about 80% of the PSF, centered on the source position. The background was extracted from a nearby circular region of 40 pixel radius. The ancillary response files (ARFs) were generated with the xrtmkarf task applying corrections for PSF losses and CCD defects using the cumulative exposure map. The latest response matrices (version 14) available in the Swift CALDB were used.

Before the spectral fitting, the 0.3–10 keV source spectra were binned using the grppha task to ensure a minimum of 20 counts per bin. Spectra were modeled in Xspec (version 12.8.0) using power-law and log-parabolic models, identical to the modeling presented in detail in Section 3.2. The models include photoelectric absorption by a fixed column density estimated to be cm⁻² (Kalberla et al. 2005). The log-parabolic model was found to fit the data better in each observation (though statistical improvement is marginal in some cases) and was therefore used to compute fluxes in various subbands. Spectral parameters are provided for each observation in Table 2.

Table 2.  Models Fitted to the Swift-XRT Spectra of Each Observation

	Power Law		Log-parabola ( keV)			Time-averaged Fluxa
Start Time (MJD)	Γ	/dof	α	β	/dof	0.3–3 keV	3–7 keV	2–10 keV
56294.7906	2.86 ± 0.03	254/199	2.85 ± 0.03	0.23 ± 0.09	233/198	20.0 ± 0.5	2.7 ± 0.2	5.6 ± 0.4
56296.9370	2.68 ± 0.02	281/227	2.64 ± 0.02	0.28 ± 0.07	235/226	27.0 ± 0.3	4.8 ± 0.2	9.9 ± 0.6
56300.1523	2.75 ± 0.04	171/144	2.74 ± 0.04	0.1 ± 0.1	168/143	10.0 ± 0.4	1.8 ± 0.1	3.6 ± 0.3
56302.1557	2.80 ± 0.03	219/198	2.79 ± 0.03	0.16 ± 0.09	208/197	13.0 ± 0.2	2.0 ± 0.2	4.1 ± 0.3
56302.3418	2.82 ± 0.03	199/183	2.79 ± 0.03	0.2 ± 0.1	187/182	14.0 ± 0.2	2.1 ± 0.2	4.3 ± 0.3
56302.4751	2.89 ± 0.03	175/172	2.87 ± 0.03	0.3 ± 0.1	159/171	13.0 ± 0.4	1.5 ± 0.1	3.3 ± 0.2
56302.6764	2.74 ± 0.02	335/275	2.71 ± 0.02	0.19 ± 0.06	304/274	14.0 ± 0.2	2.4 ± 0.1	5.0 ± 0.2
56302.9601	2.75 ± 0.04	139/155	2.73 ± 0.04	0.2 ± 0.1	132/154	15.0 ± 0.2	2.5 ± 0.3	5.1 ± 0.2
56304.4790	2.60 ± 0.02	252/218	2.56 ± 0.03	0.26 ± 0.08	219/217	21.0 ± 0.6	4.2 ± 0.2	8.4 ± 0.3
56307.0928	2.49 ± 0.01	500/394	2.45 ± 0.01	0.22 ± 0.03	385/393	33.0 ± 0.2	8.3 ± 0.4	16.0 ± 0.4
56307.3519	2.60 ± 0.01	728/431	2.56 ± 0.01	0.26 ± 0.03	477/430	33.0 ± 0.2	6.6 ± 0.2	13.0 ± 0.3
56310.1675	2.85 ± 0.03	198/161	2.83 ± 0.04	0.3 ± 0.1	180/160	12.0 ± 0.3	1.6 ± 0.1	3.5 ± 0.3
56312.2389	2.72 ± 0.01	403/351	2.69 ± 0.01	0.17 ± 0.04	350/350	14.0 ± 0.1	2.5 ± 0.1	5.1 ± 0.1
56314.5070	2.69 ± 0.03	232/197	2.68 ± 0.03	0.07 ± 0.10	229/196	16.0 ± 0.3	3.2 ± 0.2	6.4 ± 0.3
56317.3009	2.71 ± 0.03	165/171	2.69 ± 0.04	0.2 ± 0.1	158/170	17.0 ± 0.3	3.0 ± 0.3	6.2 ± 0.4
56320.3057	2.58 ± 0.02	278/216	2.53 ± 0.03	0.31 ± 0.08	238/215	25.0 ± 0.6	5.1 ± 0.3	10.0 ± 0.5
56324.6601	2.35 ± 0.02	302/261	2.31 ± 0.03	0.21 ± 0.06	270/260	49.0 ± 0.7	15.0 ± 0.5	30 ± 1
56327.1409	2.51 ± 0.02	307/250	2.45 ± 0.02	0.33 ± 0.07	231/249	40.0 ± 0.5	9.1 ± 0.3	18.0 ± 0.9
56329.0586	2.39 ± 0.01	1123/530	2.33 ± 0.01	0.24 ± 0.02	687/529	31.0 ± 0.1	9.0 ± 0.1	18.0 ± 0.3
56333.1279	2.41 ± 0.02	307/256	2.37 ± 0.02	0.21 ± 0.06	277/255	31.0 ± 0.6	8.8 ± 0.4	17.0 ± 0.6
56335.0700	2.51 ± 0.01	528/444	2.49 ± 0.01	0.10 ± 0.03	496/443	34.0 ± 0.3	8.9 ± 0.1	18.0 ± 0.2
56338.0045	2.76 ± 0.03	172/169	2.74 ± 0.04	0.2 ± 0.1	166/168	17.0 ± 0.3	2.8 ± 0.3	5.8 ± 0.6
56340.0047	2.61 ± 0.01	586/419	2.58 ± 0.01	0.17 ± 0.03	493/418	20.0 ± 0.1	4.2 ± 0.1	8.4 ± 0.2
56342.1393	2.53 ± 0.03	243/187	2.51 ± 0.04	0.1 ± 0.1	236/186	19.0 ± 0.3	4.5 ± 0.5	9.0 ± 0.3
56346.3481	2.75 ± 0.03	223/196	2.72 ± 0.03	0.27 ± 0.09	198/195	21.0 ± 0.6	3.3 ± 0.2	6.9 ± 0.4
56350.3573	2.60 ± 0.03	221/196	2.57 ± 0.03	0.20 ± 0.09	206/195	18.0 ± 0.4	3.8 ± 0.3	7.6 ± 0.5
56352.3675	2.46 ± 0.02	269/243	2.44 ± 0.03	0.14 ± 0.07	256/242	24.0 ± 0.5	6.5 ± 0.4	13.0 ± 0.8
56355.9845	2.52 ± 0.02	279/256	2.48 ± 0.02	0.21 ± 0.07	249/255	55.0 ± 0.8	13.0 ± 0.6	26.0 ± 0.6
56356.0538	2.55 ± 0.01	692/444	2.52 ± 0.01	0.21 ± 0.03	518/443	45.0 ± 0.4	10.0 ± 0.2	20.0 ± 0.4
56358.3190	2.39 ± 0.03	230/213	2.41 ± 0.03	−0.07 ± 0.08	228/212	48 ± 1	17 ± 1	33 ± 2
56363.0045	2.57 ± 0.01	837/476	2.53 ± 0.01	0.24 ± 0.02	517/475	34.0 ± 0.2	7.3 ± 0.1	15.0 ± 0.2
56366.2540	2.40 ± 0.02	230/245	2.38 ± 0.03	0.12 ± 0.07	222/244	27.0 ± 0.3	8.4 ± 0.4	16.0 ± 0.9
56368.0609	2.40 ± 0.01	758/516	2.37 ± 0.01	0.14 ± 0.02	643/515	34.0 ± 0.2	10.0 ± 0.2	20.0 ± 0.2
56369.0665	2.37 ± 0.02	352/301	2.33 ± 0.02	0.19 ± 0.05	315/300	44.0 ± 0.5	14.0 ± 0.6	27.0 ± 0.8
56370.0675	2.36 ± 0.02	323/308	2.34 ± 0.02	0.09 ± 0.05	313/307	43.0 ± 0.6	14.0 ± 0.4	28.0 ± 0.8
56374.2797	2.14 ± 0.01	539/425	2.08 ± 0.02	0.18 ± 0.03	457/424	78.0 ± 0.7	35.0 ± 0.6	68 ± 2

Note.

^aFlux calculated from the best-fit model, in units of 10⁻¹¹ erg s⁻¹ cm⁻².

Download table as:  ASCIITypeset image

NuSTAR (Harrison et al. 2013) is a focusing hard X-ray telescope operating in the band from 3 to 79 keV. It is the first X-ray observatory to extend the sensitivity beyond the ≃10 keV cutoff shared by virtually every current focusing X-ray satellite. The inherently low background associated with concentrating target X-rays enables NuSTAR to achieve approximately a 100-fold improvement in sensitivity over the collimated and coded-mask instruments that operate, or have operated, in the same bandpass. All observations are conducted in parallel with two coaligned, independent telescopes called FPMA and FPMB (for Focal Plane Module A and B).

The NuSTAR primary mission includes monitoring of several types of blazars; Mrk 421 has been selected for this program as a representative of the high-peaked BL Lac (HBL) class. In order to maximize the strictly simultaneous overlap of observations by NuSTAR and ground-based VHE γ-ray observatories during the 5-month campaign, three observations per month were scheduled according to visibility of Mrk 421 at the MAGIC and VERITAS sites. A typical NuSTAR observation spanned 10 hr, resulting in 15–20 ks of source exposure after accounting for orbital modulation of visibility and filtering out South Atlantic Anomaly crossings where the background radiation is high. In addition to those observations, Mrk 421 was observed as a bright calibration target in 2012 July and early 2013 January. The total exposure time over 88 orbits of NuSTAR observations in this period is ≃250 ks. A list of all NuSTAR observations considered in this paper is given in Table 3. Analysis of the remainder of the campaign data will be presented elsewhere.

The raw data have been reduced using the NuSTARDAS software version 1.3.1, as a part of the HEAsoft package version 6.12. The spectra of Mrk 421 were extracted from a circular region of 100'' radius centered on the peak of the distribution of cleaned events. Background spectra were extracted from a region encompassing the same detector on which the source was focused, excluding the circular region from which the source counts were extracted. As the background generally differs between different detectors and may be variable on few-orbit timescales, extraction from a region of maximal area on the same detector where the source is present provides the best background estimate over the NuSTAR band. Nevertheless, other background extractions have been attempted and no significant differences have been observed in the results.

The response files were generated using the standard nupipeline and nuproducts scripts and the calibration files from CALDB version 20131223. All flux values reported in this paper have been corrected for the finite extraction aperture by the processing software. The dominant background component above 25 keV is the internal detector background. With good background characterization, the data may be used for spectral modeling up to the high-energy end of the NuSTAR band at 79 keV. The spectra of all NuSTAR observations of Mrk 421 are above the background level at least up to 25 keV and up to ≈40 keV in observations at high flux. For this reason, we quote count rates only up to 30 keV in the remainder of the paper. Three faint serendipitous sources have been found in the NuSTAR field of view (detected only in the deep co-added image using all observations presented in G. B. Lansbury et al., in preparation); however, they do not represent a contamination problem owing to the overwhelming brightness of Mrk 421 in all epochs.

The Large Area Telescope (LAT) on board the Fermi satellite is a pair-conversion telescope with energy coverage from 20 MeV to . LAT has a ∼2.4 sr field of view and provides all-sky monitoring coverage on a ∼3 hr timescale (Atwood et al. 2009). For the analyses presented in this paper, we have selected Source class events with energies in the range 0.1–300 GeV and within 15° of the position of Mrk 421. In order to greatly reduce contamination from Earth limb photons, we have excluded events at zenith angles and any events collected when the spacecraft rocking angle was . The data were analyzed using the P7REP_SOURCE_V15 instrument-response functions and the standard unbinned-likelihood software provided with version 09-33-00 of the Fermi Science Tools.¹⁰⁶

The analyses considered data in day-long and week-long intervals contemporaneous with the NuSTAR observation windows. The likelihood model used for all intervals included the sources from the second Fermi-LAT catalog (Nolan et al. 2012) located within a region of interest centered on Mrk 421, as well as the standard Galactic diffuse, isotropic, and residual instrumental background emission models provided by the Fermi Science Support Center.¹⁰⁷ For all epochs, the spectrum of Mrk 421 was fitted with a power-law model, with both the flux normalization and photon index being left as free parameters in the likelihood fit. We summarize the spectral parameters for four selected epochs (discussed in detail in Section 5) in Table 4. The systematic uncertainty on the flux is estimated as approximately 5% at 560 MeV and under 10% at 10 GeV and above (Ackermann et al. 2012). As variability in the Fermi-LAT band was not significant, these epochs may be considered representative of the entire 2013 January–March period.

MAGIC is a system of two 17 m diameter imaging air-Cherenkov telescopes (IACTs) located at the Roque de los Muchachos Observatory on La Palma, one of the Canary Islands (28°46' N, 17°534 W at 2231 m above sea level). The hardware was substantially upgraded during 2011 and 2012 (Aleksić et al. 2016a), which yielded a performance characterized by a sensitivity of ≃0.7% of the Crab Nebula flux to detect a point-like source above 200 GeV at 5σ in 50 hr of observation. Equivalently, a 1 hr integration yields a detection of a source with approximately 5% Crab flux. The angular resolution is ≲007 (68% containment, GeV), and the energy resolution is 16%. The systematic uncertainties in the spectral measurements for a Crab-like point source were estimated to be 11% in the normalization factor (at ≃200 GeV) and 0.15 in the power-law slope. The systematic uncertainty in the absolute energy determination is estimated to be 15%. Further details about the performance of the MAGIC telescopes after the hardware upgrade in 2011–2012 can be found in Aleksić et al. (2016b).

After data-quality selection, Mrk 421 was observed with MAGIC for a total of 10.8 hr between 2013 January 8 and 2013 March 18. Most of these observations (8 hr in total) were strictly simultaneous with the NuSTAR observations. They were performed in the "false-source tracking" mode (Fomin et al. 1994), where the target source position has an offset of 04 from the camera center, so that both signal and background data are taken simultaneously. Those data were analyzed following the standard procedure described in Aleksić et al. (2016b), using the MAGIC Analysis and Reconstruction Software (MARS; Moralejo et al. 2009). The analysis cuts to extract γ-ray signals from the hadronic background were optimized on the Crab Nebula data and dedicated Monte Carlo simulations of γ-ray-induced showers.

The significance of the source detection, calculated using Equation (17) from Li & Ma (1983), varied between 8.3σ (MJD 56310) and 38.3σ (MJD 56335). Observed intranight variability is not statistically significant, so we used data integrated over complete observations for the spectral analysis. Spectra were modeled with a power-law function with a normalization energy of 300 GeV. The normalization energy was chosen to be 300 GeV for both MAGIC and VERITAS, in order to facilitate a direct comparison of the VHE spectral results. For observations in which the spectrum is not well described by a power-law model (MJD 56302, 56335, and 56363), we additionally fit a log-parabolic model. A summary of the observations and the spectral modeling is given in Table 5. All uncertainties quoted in the table and in the rest of the paper are statistical only.

Table 5.  Summary of the MAGIC Observations of Mrk 421 (2013 January–March)

				Model Fitb,c
Start Time (MJD)	Exposure (minutes)	Zenith Angle (deg)	σa	F0	Γ or α	β	dof	Flux > 200 GeVc (10−11 s−1 cm−2)
56302.1365	122	9–23	16.5	1.7 ± 0.1	2.87 ± 0.07	...	19.7 / 16	...
				1.9 ± 0.1	3.2 ± 0.2	0.8 ± 0.3	9.7 / 15	5.7 ± 0.4
56307.2556	39	21–29	14.9	3.0 ± 0.2	2.48 ± 0.09	...	18.8 / 15	10.4 ± 0.9
56310.2441	54	20–31	8.3	1.3 ± 0.2	2.8 ± 0.1	...	11.6 / 14	3.7 ± 0.5
56312.1718	119	9–33	11.5	1.3 ± 0.1	2.9 ± 0.1	...	8.8 / 20	3.9 ± 0.4
56316.2417	29	24–30	10.9	2.4 ± 0.2	2.3 ± 0.1	...	21.5 / 17	8.1 ± 0.9
56327.0731	25	16–22	25.7	8.6 ± 0.4	2.27 ± 0.05	...	20.5 / 19	34 ± 2
56333.1147	29	9–10	15.8	4.4 ± 0.3	2.34 ± 0.09	...	15.5 / 15	16 ± 1
56335.0795	116	9–24	38.3	5.7 ± 0.2	2.52 ± 0.03	...	25.9 / 19	...
				6.4 ± 0.3	2.53 ± 0.04	0.33 ± 0.09	9.6 / 18	20.5 ± 0.7
56340.1722	29	23–36	13.5	3.0 ± 0.3	2.4 ± 0.1	...	18.7 / 14	10 ± 1
56362.0826	29	15–21	20.1	6.0 ± 0.4	2.36 ± 0.07	...	22.1 / 19	20.7 ± 1.4
56363.1066	56	23–33	29.1	5.8 ± 0.2	2.56 ± 0.04	...	34.0 / 20	...
				6.8 ± 0.4	2.59 ± 0.06	0.56 ± 0.06	18.6 / 19	20 ± 1

Notes.

^aDetection significance, computed according to Equation (17) from Li & Ma (1983) using data above 200 GeV integrated over the exposure time. ^bA power-law model of the form , , is fitted to each observation. A log-parabolic model fit of the form is shown for observations in which it provides a better description of the spectrum than the power-law model. The normalization constant, F₀, is in units of 10⁻¹⁰ s⁻¹ cm⁻² TeV⁻¹. ^cQuoted uncertainties are statistical.

Download table as:  ASCIITypeset image

The Very Energetic Radiation Imaging Telescope Array System (VERITAS) is an array of four 12 m diameter IACTs located in southern Arizona (Weekes et al. 2002; Holder et al. 2006) designed to detect emission from astrophysical objects in the energy range from ∼100 GeV to greater than 30 TeV. A source with 1% of the Crab Nebula flux can be detected in ≃25 hr of observations; equivalently, a source with approximately 5% Crab flux can be detected in a 1 hr integration. VERITAS has an energy resolution of ≃15% and an angular resolution (68% containment) of per event at 1 TeV. The uncertainty on the VERITAS energy calibration is approximately 20%. The systematic uncertainty on reconstructed spectral indices is estimated at 0.2, independent of the source spectral index, according to studies of Madhavan (2013). Details of the sensitivity of the system after the recent hardware upgrade can be found on the VERITAS Web site.¹⁰⁸

VERITAS observations of Mrk 421 were carried out under good weather conditions during the period of the NuSTAR campaign, resulting in a total, quality-selected exposure time of 15.5 hr during the period 2013 January 10 to 2013 March 17, almost all of which is strictly simultaneous with NuSTAR exposures. These observations were taken at 05 offset in each of four cardinal directions from the position of Mrk 421 to enable simultaneous background estimation using the "false-source tracking" method (Fomin et al. 1994). Detected events are parameterized by the principal moments of the elliptical images of the Cherenkov shower in each camera (Hillas 1985). Cosmic-ray background rejection is carried out by discarding events based on a set of selection cuts that have been optimized a priori using VERITAS observations of the blazar 1ES 1218+304 (photon index 3.0) and the Crab Nebula (photon index 2.5). The results were verified using two independent analysis packages (Cogan 2008; Daniel 2008).

The significance of the source detection was calculated using Equation (17) from Li & Ma (1983) and was found to vary between 18.7σ (on MJD 56302) and 40.4σ (on MJD 56368). No significant intranight variability was detected. Since the observations spanned several hours during each night, the energy threshold varied owing to the range of zenith angles observed. Fluxes are therefore quoted at a commonly reached minimum energy of 200 GeV. We modeled the spectra with a power-law function with a normalization energy of 300 GeV. For four observations (MJD 56302, 56356, 56363, and 56368), the spectrum is better described with a log-parabolic model, while for the other observations this model does not provide a significantly better fit than the simpler power-law model. A summary of the observations, VHE flux, and spectral parameters with their statistical uncertainties is given in Table 6.

Table 6.  Summary of the VERITAS Observations of Mrk 421 (2013 January–March)

				Model Fitb,c
Start Time (MJD)	Exposure (minutes)	Zenith Angle (deg)	σa	F0	Γ or α	β	dof	Flux >  200 GeVc (10−11 s−1 cm−2)
56302.3411	131	9–33	18.7	1.7 ± 0.2	3.2 ± 0.2	...	10.2 / 5	...
				1.8 ± 0.2	3.3 ± 0.2	1.1 ± 0.7	7.9 / 4	5.2 ± 0.6
56307.3487	170	8–28	34.5	3.6 ± 0.2	3.1 ± 0.1	...	15.7 / 9	12.2 ± 0.7
56312.3762	197	6–31	32.6	1.5 ± 0.1	3.2 ± 0.2	...	4.9 / 5	5.0 ± 0.4
56329.2864	49	7–28	25.2	4.5 ± 0.5	2.7 ± 0.2	...	8.2 / 5	17 ± 2
56340.3344	89	7–23	22.7	2.5 ± 0.2	3.2 ± 0.2	...	5.1 / 5	8.5 ± 0.8
56356.2352	43	6–21	33.6	6.0 ± 0.5	3.0 ± 0.1	...	10.4 / 8	...
				6.3 ± 0.5	3.0 ± 0.1	0.6 ± 0.4	8.2 / 7	20 ± 2
56363.2355	127	7–17	39.2	3.5 ± 0.2	3.5 ± 0.1	...	16.7 / 6	...
				4.1 ± 0.3	3.9 ± 0.1	1.9 ± 0.4	5.6 / 5	10.9 ± 0.8
56368.1885	123	9–26	40.4	4.3 ± 0.3	3.2 ± 0.1	...	9.2 / 9	...
				4.3 ± 0.3	3.2 ± 0.1	0.7 ± 0.3	4.6 / 8	13 ± 1

Notes.

^aDetection significance, computed according to Equation (17) from Li & Ma (1983) using data above 200 GeV integrated over the exposure time. ^b A power-law model of the form , , is fitted to each observation. A log-parabolic model fit of the form is shown for observations in which it provides a better description of the spectrum than the power-law model. The normalization constant, F₀, is in units of 10⁻¹⁰ s⁻¹ cm⁻² TeV⁻¹. ^cQuoted uncertainties are statistical.

Download table as:  ASCIITypeset image

Figure 1 shows the background-subtracted X-ray light curves extracted from the NuSTAR observations of Mrk 421 listed in Table 3. The split into subbands at 7 keV is based on the spectral analysis and justified in later sections. The count rate above 30 keV is dominated by the background on short timescales and is therefore not shown here. The differences in count rates between observations, as well as the range covered in each particular observation, are entirely dominated by the intrinsic variability of the target. For example, the calibration observation taken in 2012 July (the top panel of Figure 1) includes a possible "flare" in which the count rate increased by a factor of 2.5 over a 12 hr period and dropped by almost a factor of two in 3 hr. Several observations in 2013 have shown a steadily decreasing count rate over the course of ≃12 hr. We did not observe any sharp increases followed by exponential decay typical of flaring events, although we cannot exclude the possibility that the observed count rate decreases are due to such events. All of the observed increases in the count rate (e.g., 2012 July 7 and 8, 2013 February 6, as well as 2013 March 5 and 17 on a shorter timescale) appear rather symmetric with respect to subsequent decreases. The campaign observations up to the end of 2013 March have predominantly covered relatively low flux states of Mrk 421, even though the lowest and the highest observed fluxes span almost an order of magnitude.

Zoom In Zoom Out Reset image size

The observed count rates are not consistent with a constant flux during any of the observations. However, the dominant variations in the count rate can be described as smooth on a timescale of several hours. If a simple exponential decay fit, , is performed on the observations that show significant downward trends (2012 July 8, 2013 January 15, February 12 and 17, March 5 and 12), the typical decay timescale () is found to be between 6 and 12 hr. These rough fits are not meant to describe the light curves fully, but only to provide an estimate of the timescale on which the flux changes significantly. For the remainder of the paper we use hr as our best estimate of this timescale.

In order to characterize the variability on shorter timescales (), we consider the data in individual NuSTAR orbits as this represents a natural, albeit still arbitrary, way of partitioning the data. The NuSTAR orbits are approximately 90 minutes long and contain roughly 50 minutes of source exposure. We treat each orbit independently and fit two simple light-curve models to the observed count rates: a constant flux during the orbit, constant, and a linear trend in time, i.e., . The upper panel of Figure 2 provides an example of both models fitted to the 2012 July data binned into 10-minute time bins, so that each orbit is divided into 5–7 bins per focal-plane module.¹⁰⁹ The lower panel of Figure 2 shows the results of this fitting procedure performed on all 88 orbits.

Zoom In Zoom Out Reset image size

We find that the flux during the majority of orbits is better described by a linear trend than by a constant-flux model. Linear trends account for most of the orbit-to-orbit variability and approximate smooth variations on super-orbital timescales of hr. In 10-minute bins, for example, the variability amplitude typically does not exceed the observed count-rate uncertainty of 3%. Based on the mildly overpopulated tail of the reduced- distribution for the linear-trend fits, we estimate that up to 20% of orbits show excess variance beyond the simple linear trend. Subtracting that trend and comparing the residual scatter to the median rate uncertainty within each orbit, , gives a distribution slightly skewed toward values greater than unity (see lower right panel of Figure 2). This is consistent with intrinsic suborbital variability on a ∼10-minute timescale in ≲20% of orbits, while for the majority of the observations the short-timescale variability can be constrained to have a ≲5% amplitude. These results are independent of the exact choice of the bin size and hold for any subhour binning. Based on a separate analysis of the low- and high-flux data alone, we do not find significant evidence for a change in variability characteristics with flux.

The hardness of the spectrum, defined here as the ratio of count rates in the hard 7–30 keV and in the soft 3–7 keV bands, changes over the course of the observations. In Figure 3 we show the general trend of the spectrum hardening when the count rate is higher. Although the observed range of hardness ratios is relatively large at any specific count rate, the overall trend is clearly present in the binned data shown with thick black lines. There are no apparent circular patterns observed in the count rate versus hardness ratio plane, as previously seen in soft X-ray observations during bright flaring periods (e.g., Takahashi et al. 1996; Ravasio et al. 2004; Tramacere et al. 2009). We note, however, that the circular patterns might not be observable in the NuSTAR data presented here simply because the observations predominantly covered instances of declining flux, whereas the patterns arise from differences in the rising and the declining phases of a flare. The apparent symmetric features observed on 2012 July 7/8, 2013 February 6, and 2013 March 17 do not display enough contrast in flux and hardness ratio to show well-defined circular patterns.

Zoom In Zoom Out Reset image size

We first model the NuSTAR spectra for each observation separately, before examining the clear intra-observation spectral variability in Section 3.3 (see Figures 2 and 3, and Baloković et al. 2013a). All observation-averaged spectra are shown in Figure 4. For spectral analysis we use spectra grouped to a minimum of 20 photons per bin and perform the modeling in Xspec (version 12.8.0). The simplest model for a featureless blazar spectrum is a power-law function:

where Γ is the photon index. The Xspec model is formulated as phabs(zpow), where the phabs component accounts for the Galactic absorption with fixed hydrogen column density of cm⁻² (Kalberla et al. 2005). We first fit each of the 12 observation-averaged spectra with a power-law model and find that this model fits the observations with lower mean flux better than the ones where the mean flux is high (see Table 7). This is likely due to the fact that the higher-flux spectra are somewhat more curved than lower-flux ones, although the curvature is not immediately obvious to the eye, i.e., in Figure 4. The fits confirm that the effective photon index decreases with increasing flux, as suggested by the observed harder-when-brighter behavior shown in Figure 3.

Zoom In Zoom Out Reset image size

Table 7.  Models Fitted to the NuSTAR Spectra of Each Observation

	Power Law		Broken Power Law				Log-parabola ( keV)			Time-averaged Fluxb
Start Time (MJD)	Γ	/dof		Eba		/dof	α	β	/dof	3–7 keV	2–10 keV	7–30 keV
56115.1353	2.82 ± 0.01c	1085/949	2.74 ± 0.02		2.92 ± 0.03	922/947	2.76 ± 0.01	0.21 ± 0.03	906/948	9.47 ± 0.03	18.8 ± 0.1	6.48 ± 0.03
56116.0732	2.87 ± 0.01c	1126/833	2.78 ± 0.03		2.98 ± 0.04	978/831	2.99 ± 0.01	0.24 ± 0.04	967/832	11.08 ± 0.03	22.1 ± 0.1	7.16 ± 0.04
56294.7778	3.10 ± 0.04	415/390			2.9 ± 0.1	399/388	3.16 ± 0.05	−0.2 ± 0.1	403/389	2.87 ± 0.05	6.4 ± 0.1	1.57 ± 0.04
56302.0533	3.07 ± 0.03	512/506	3.08 ± 0.04	7.5 (f)	3.04 ± 0.07	511/505	3.07 ± 0.04	0.0 ± 0.1	512/505	1.91 ± 0.02	4.10 ± 0.05	1.05 ± 0.01
56307.0386	3.02 ± 0.01	865/710	2.89 ± 0.03			742/708	2.92 ± 0.02	0.28 ± 0.05	741/710	6.43 ± 0.03	13.0 ± 0.1	3.55 ± 0.02
56312.0980	3.05 ± 0.02	571/543				568/541	3.03 ± 0.03	0.1 ± 0.1	572/542	2.26 ± 0.02	4.71 ± 0.05	1.22 ± 0.01
56329.0116	2.93 ± 0.01	925/724	2.80 ± 0.02		3.13 ± 0.04	709/722	2.82 ± 0.02	0.39 ± 0.05	709/723	8.27 ± 0.05	16.3 ± 0.1	4.94 ± 0.04
56335.0106	2.73 ± 0.01	839/742	2.66 ± 0.02			738/740	2.64 ± 0.02	0.21 ± 0.05	742/741	9.25 ± 0.05	18.1 ± 0.1	7.00 ± 0.06
56339.9828	3.02 ± 0.02	577/559	2.95 ± 0.03	7.5 (f)	3.09 ± 0.05	559/558	2.96 ± 0.03	0.14 ± 0.09	558/558	3.83 ± 0.03	7.9 ± 0.1	2.14 ± 0.03
56355.9631	3.01 ± 0.01	823/701			3.09 ± 0.04	751/699	2.94 ± 0.02	0.22 ± 0.05	751/700	9.74 ± 0.04	19.9 ± 0.1	5.45 ± 0.05
56362.9690	3.10 ± 0.01	640/640				603/638	3.04 ± 0.02	0.19 ± 0.06	600/639	7.49 ± 0.05	15.6 ± 0.1	3.80 ± 0.04
56368.0210	2.75 ± 0.01	848/760	2.67 ± 0.02		2.91 ± 0.05	756/758	2.68 ± 0.02	0.24 ± 0.05	749/759	9.51 ± 0.05	18.6 ± 0.1	6.93 ± 0.05

Notes.

^aEnergy (in keV) at which the model sharply changes slope from to . For some observations this parameter is unbound, so we fix it to 7.5 keV and mark with (f). ^bFlux calculated from the best-fit model, in units of 10⁻¹¹ erg s⁻¹ cm⁻². Note that the 2–10 keV band requires some extrapolation below the NuSTAR bandpass, but we provide it here for easier comparison with the literature. ^cFormal statistical uncertainty is 0.008; however, the NuSTAR bandpass calibration is limited to 0.01 (Madsen et al. 2015), so we round up the values assuming this lowest uncertainty limit for these cases.

Download table as:  ASCIITypeset image

The fitting results imply that a power-law model with describes the data well for observations with the lowest flux observed in the campaign. The poorer fit of the power-law model for the high-flux observation data may be due to intrinsic curvature, or it may be simply an effect of superposition of different curved or broken-power-law spectra. The latter effect can certainly be expected to be present since the hardness does vary with the flux and the source exhibited significant variability during most of the observations (see Figure 1). We address this issue by examining spectra on a shorter timescale in Section 3.3. In order to better characterize the observation-averaged spectra, we replace the power-law model with two other simple models that allow extra degrees of freedom. The first one is a broken-power-law model, bknpow:

This model provides better fits to highest-flux spectra. However, the broken-power-law form is degenerate at low flux and degrades into the simpler power-law shape discussed above for observations of low mean flux (i.e., the photon indices converge to a single value and E_b becomes unconstrained). Both photon indices depend on the flux, but the dependence of seems to be weaker. The break energy seems to be largely independent of flux and relatively poorly constrained to the range roughly between 5 and 10 keV.

The third Xspec model we use is a simple log-parabolic shape, logpar:

In this model, α and β are free parameters, while E_* is the so-called pivot energy (fixed parameter) at which α equals the local power-law photon index. The β parameter describes deviation of the spectral slope away from E_*. In our analysis we fix the value of E_* to 5 keV, so that α closely approximates the photon index in the 3–7 keV band. This model fits all observations well and also hints at spectral trends outlined earlier. The log-parabolic model does not provide statistically better fits than the broken-power-law model; in most cases they fit equally well (see Table 7). However, the log-parabolic is often used for modeling blazar spectra in the literature and does not contain a rather unphysical sharp break in the spectrum. All relevant parameters of the fits to the observation-averaged spectra are summarized in Table 7. We note that for the NuSTAR observation on 2013 January 2 (MJD 56294), the best-fit parabolic model has marginally significant negative curvature, . As this is the shortest of all NuSTAR observations and the same effect is not observed in any of the other spectra, while negative curvature is a physically possible scenario, it is likely that this anomalous result arises from the fact that the high-energy background is not sufficiently well sampled in such a short exposure.

We next consider spectra integrated over time intervals shorter than the complete NuSTAR observations. Separating the data into individual orbits represents the most natural although still arbitrary way of partitioning. Any particular orbit has a smaller spread in flux compared to a complete observation, since variability amplitude is significantly lower—we have established in Section 3.1 that the dominant flux variations occur on a super-orbital timescale of hr. The shorter exposures significantly reduce the statistical quality but still allow for a basic spectral analysis, such as the one presented in the preceding section, to be performed on spectra from single orbits. The average orbit exposure is 2.8 ks, and the total number of source counts per orbit is between 2000 and 20,000 per focal-plane module.

As with the observation-averaged spectra, we fit power-law, broken-power-law, and log-parabolic models using Xspec. We again find that the broken-power-law model parameter E_b (the break energy) is poorly constrained in general, so we fix it at 7 keV for the remainder of this analysis. Choosing a different value in the interval between 5 and 10 keV does not significantly alter any results; however, break energies outside of that interval cause one of the photon indices to become poorly constrained in a considerable number of orbits. Similarly, the curvature parameter of the log-parabolic model (β) is poorly constrained for the lowest-flux data, likely owing to both the lack of intrinsic curvature and relatively poor photon statistics. In a similar manner to the observation-averaged spectral modeling, the log-parabolic model does not necessarily provide statistically better fits than the broken power-law model. We use it because it provides a smooth spectrum over the modeled energy range, and in order to facilitate comparison to other work in the literature.

With less smearing over different spectral states of the source, the spectral variability is more clearly revealed by this analysis. As shown with the gray data points in Figure 5 (one for each NuSTAR orbit), for any of the three models statistically significant spectral changes occur as the X-ray flux varies. The spectrum becomes harder as the flux increases. Most of the change happens below ≃7 keV, as shown by the substantial variations in the parameters and α, compared to the lower-amplitude variations in and β. In all cases the trends are consistent with the well-established harder-when-brighter behavior, also evident in the more basic representation using hardness ratios in Figure 3. Since for orbits with the lowest count rate the uncertainties on the spectral parameters are relatively large, in the following section we verify that the same spectral variability trends persist for data with higher signal-to-noise ratio (S/N).

Zoom In Zoom Out Reset image size

In order to verify that the spectral parameter trends we identify in the time-resolved spectra are not spuriously produced by relatively low photon statistics at the low-flux and high-energy ends, we proceed to examine stacked single-orbit spectra of similar flux. Stacking provides the highest possible S/N in well-defined flux bins and allows us to use the data up to 70 keV—where the signal is fainter by a factor of a few than the NuSTAR detector background. We combine spectra for each focal-plane module separately due to intrinsic differences in response matrices. Spectra from both modules are fitted simultaneously in Xspec, just like the observation-averaged and the single-orbit ones. Note that this procedure implicitly assumes that the source behaves self-consistently, in the sense that a particular flux level corresponds to a unique spectral shape within the data-taking time interval. The validity of this assumption is further discussed in Section 5.

We first stack the spectra of three complete observations (2013 January 2, 10, and 20), since Mrk 421 displayed a nearly constant low flux during all three (see Figure 1). The combined spectrum is very similar to the spectra from any of the constituent observations but has significantly higher S/N. It can be statistically well described as a simple power law with from 3 to 70 keV (). For completeness, modeling with a broken-power-law model gives the break energy at 6 ± 3 keV, and low- and high-energy photon indices both formally consistent with the Γ value found for the simpler power-law model. Furthermore, the curvature parameter of the log-parabolic model is consistent with zero (), and . These fitting results lead us to conclude that in the lowest-flux state observed in 2013 the hard X-ray spectrum follows a steep power law with . Extrapolating below 3 keV for the sake of comparison with the literature, we derive a 2–10 keV flux of erg s⁻¹ cm⁻² for the power-law model normalized to the lowest observed orbital flux (orbit #6 of the 2013 January 20 observation). A state of such low X-ray flux has not yet been described in the published literature,¹¹⁰ making the results of this analysis new and unique.

We combine the spectra of all 88 orbits according to their 3–7 keV flux in order to obtain higher S/N for the flux states covered in the data set. The choice of flux bins shown with black symbols in Figure 5 is such that relatively uniform uncertainty in spectral parameters is achieved across the flux range; this condition requires stacking of ∼10 orbits of data on the faint end, while a single orbit is sufficient at the bright end. The results, however, are largely independent of the exact choice of which orbits to combine into a particular flux bin. Fitting the stacked spectra with the same simple spectral models as before reveals spectral trends much more clearly than for observation-averaged or time-separated spectra, as shown by black symbols in Figure 5. The spectra of the lowest-flux stack and the highest-flux orbit are displayed in Figure 4 for comparison with the observation-averaged spectra. The analysis performed here describes the spectral changes happening between those two extremes as a smooth function of the X-ray flux.

For each of the parameters of the power-law, broken-power-law, and log-parabolic models we parameterize their dependence on the X-ray flux as

where X stands for any of the spectral parameters (Γ, , E_b, , α, β), s is the slope of the relation, is the 3–7 keV flux, F₀ is a reference flux in erg s⁻¹ (chosen to be the median flux of our data set, ), and X₀ is the vertical offset (parameter value at the reference flux). We find that in all cases this linear function adequately describes the general trends. Since we find that the break energy of the broken power-law model (E_b) is independent of flux within its large uncertainties, we keep it fixed at 7 keV while fitting for the trends in the and parameters. For , the high-energy photon index of the broken-power-law model, the best-fit slope is small but different from zero at the 2σ level. For the rest of the spectral parameters, the trends are statistically more significant. Figure 5 shows the best-fit relations superimposed on the time- and flux-resolved fitting results, clearly matching the trends that the former analysis hinted at. We list the best-fit linear trend parameters s and X₀, with their 1σ uncertainties, in Table 8.

Table 8.  Best-fit Linear Relations Parameterized in Equation (4) for Describing the Change of Spectral Parameters with the X-Ray Flux. The Uncertainties are Quoted at the 1σ Level

Model	Parameter, X	Slope, s	Zero Point, X0
pow	Γ	−0.32 ± 0.07	2.88 ± 0.02
bknpow		−0.46 ± 0.07	2.80 ± 0.02
bknpow	Eb a	1 ± 9	7 ± 2
bknpow		−0.11 ± 0.05	3.00 ± 0.01
logpar	α	−0.42 ± 0.07	2.82 ± 0.02
logpar	β	0.22 ± 0.06	0.24 ± 0.01

Note.

^aOwing to poor constraint on this parameter, it is kept fixed at 7 keV while quantifying the trends in the and parameters.

Download table as:  ASCIITypeset image

Finally, we briefly return to the broken power-law model fits only to make a comparison to the previously observed spectral variability during flares. The RXTE 2–20 keV data analyzed by Giebels et al. (2007) overlap in the 2–10 keV flux only for the highest-flux single-orbit data presented here and extend almost a decade above that. These authors showed that the break energy is essentially independent of flux and (68% confidence interval), which is consistent with the median value of approximately 7 keV found from our nondegenerate fits of the broken-power-law model. The photon indices were found to vary with flux up to approximately 10⁻⁹ erg s⁻¹ cm⁻², above which they saturate at and . The data presented here smoothly connect to those trends (see Figure 6), extending them toward the faint end. Whereas the low-energy photon index ( or ) continues to increase with decreasing flux, reaching at ≲10⁻¹⁰ erg s⁻¹ cm⁻², the high-energy one ( or ) essentially levels off to the same at a factor of a few higher flux. A naive extrapolation of the Giebels et al. (2007) trends is therefore not supported by the new NuSTAR observations. Our analysis reveals a clear low-flux saturation effect that none of the previous studies could have constrained owing to lack of sensitivity.

Zoom In Zoom Out Reset image size

The majority of observations performed in 2013 January through March were coordinated between the participating observatories so as to maximize the strictly simultaneous overlap in the X-ray and VHE bands. In particular, nine 10–12 hr long observations performed by NuSTAR were accompanied with Swift pointings at the beginning, middle, and end, and the ground-based Cherenkov-telescope arrays MAGIC and VERITAS covered roughly half of the NuSTAR exposure each. Approximately 50 hr of simultaneous observations with NuSTAR and either MAGIC or VERITAS resulted in total exposure of 23.5 hr (the remainder being lost owing to poor weather conditions and quality cuts). Figure 7 shows the multiwavelength light curves and highlights the dates of NuSTAR observations taken simultaneously with MAGIC and VERITAS observations with vertical lines. A zoomed-in view of the VHE γ-ray observation times is shown overlaid on the expanded NuSTAR light curves in Figure 1.

Zoom In Zoom Out Reset image size

The VHE flux varied between approximately 0.1 and 2 Crab units, reaching substantially lower and higher than the long-term average of 0.446 ± 0.008 Crab (Acciari et al. 2014), which is considered typical for a nonflaring state of Mrk 421 (Aleksić et al. 2015b). In Figure 7 we show typical fluxes for Fermi-LAT, Swift-XRT, and OVRO bands, represented by medians of the long-term monitoring data that are publicly available. The Fermi-LAT light curve reveals elevated flux with respect to the median, as do the optical and UV data when compared to historical values. Modest soft X-ray flux is apparent from Swift-XRT data in comparison with the long-baseline median (Baloković et al. 2013a; Stroh & Falcone 2013) and the intense flaring episodes covered in the literature (e.g., Acciari et al. 2009; Tramacere et al. 2009; Aleksić et al. 2012). Count rates of s⁻¹ in the 0.3–10 keV band are up to a factor of lower than those observed in quiescent periods during multiwavelength campaigns in 2009 (Aleksić et al. 2015b) and 2010 March (Aleksić et al. 2015c). The radio flux was only slightly elevated above the values that have remained steady for the past 30 years, apart from the exceptional radio flare observed in 2012 October (see Section 4.5 for more details).

Remarkably well correlated flux variability in the X-ray and VHE bands on a timescale of about a week is already apparent from Figure 7 and will be discussed in more detail in Section 4.2.1. The fluxes in the UV and Fermi-LAT bands (to the extent allowed by the limited photon statistics) are consistent with a slow increase in flux between January and March but do not show a clear short-term flux correlation. Further details regarding these bands are presented in Section 4.2.2. The activity observed in the first 3 months of 2013 can be generally described as low. Note in particular that on January 10 and 20, Mrk 421 showed a remarkably low X-ray and VHE flux in comparison to the historical X-ray and VHE fluxes reported in Stroh & Falcone (2013) and Acciari et al. (2014), respectively. Optical polarization, shown in Figure 8, showed random and statistically significant variations around the average polarized fraction of 3%, and the polarization angle also varied significantly without any obvious coherent structure.

Zoom In Zoom Out Reset image size

A general trend observed in the 2013 campaign is a gradual rise in broadband emission between January and March by a factor of ≲10, depending on the band. This was followed by an intense flaring period in 2013 April (not shown in Figure 7), rivaling the brightest flares ever observed for Mrk 421 (Baloković et al. 2013b; Cortina et al. 2013; Paneque et al. 2013; Pian et al. 2014). Analysis of the campaign data from the flaring period and more detailed analysis of the multiwavelength variability properties will be presented in separate publications. In the following sections, we focus on quantifying short-timescale and time-averaged correlations between different spectral bands and on the basic modeling of the Mrk 421 SED in the low-activity state that has not previously been characterized in any detail, except very recently in Aleksić et al. (2015b).

The variability across the electromagnetic spectrum can be described using the fractional variability distribution. Fractional variability, , is mathematically defined in Vaughan et al. (2003), and its uncertainty is calculated following the prescription from Poutanen et al. (2008), as described in Aleksić et al. (2015a). It can be intuitively understood as a measure of the variability amplitude, with uncertainty primarily driven by the uncertainty in the flux measurements and the number of measurements performed. While the systematic uncertainties on the absolute flux measurements¹¹¹ do not directly add to the uncertainty in , it is important to stress that different observing sampling and, more importantly, different instrument sensitivity do influence and its uncertainty: a densely sampled light curve with very small temporal bins and small error bars might allow us to see flux variations that are hidden otherwise, and hence we might obtain a larger . Some practical issues of its application in the context of multiwavelength campaigns are elaborated in Aleksić et al. (2014, 2015b, 2015c).

In this paper we explore two cases, as shown in Figure 9. First, we use the full January–March data set reported in Figure 7 (which has different cadence and different number of observations in each band), and second, we use only data collected simultaneously, in narrow windows centered on observations coordinated between NuSTAR and VHE telescopes. In the latter case the fluxes are averaged over the complete NuSTAR, Swift, and VHE observations, effectively smoothing over any variability on shorter timescales. The optical and radio fluxes are taken from single measurements closest in time to the coordinated observations. In the former case, however, we sample shorter timescales and during a longer time span, which allows us to detect somewhat higher variability, as one can infer by comparing the Swift and MAGIC observations reported by the open/filled markers in Figure 9. for Fermi-LAT is calculated from the weekly binned light curve shown in Figure 7; the relatively low GeV γ-ray flux observed by Fermi-LAT precludes us from using significantly shorter time bins, or dividing the Fermi-LAT band into subbands as we do for Swift-XRT and NuSTAR. Figure 9 shows that determined from our campaign rises significantly from the radio toward the X-ray band (consistent with Giebels et al. 2007), decreases over the Fermi-LAT band (consistent with Abdo et al. 2011), and then rises again in the VHE band. This double-bump structure relates to the two bumps in the broadband SED shape of Mrk 421 and has been recently reported for both low activity (Aleksić et al. 2015b) and high activity (Aleksić et al. 2015c). The less variable energy bands (radio, optical/UV, and GeV γ-ray bands) relate to the rising segments of the SED bumps, while the most variable energy bands (X-rays and VHE γ-ray bands) relate to the falling segments of the SED bumps.

Zoom In Zoom Out Reset image size

4.2.1. X-Ray versus VHE γ-Ray Band

The existence of a correlation between the X-ray and VHE fluxes is well established on certain timescales and in certain activity states of Mrk 421: claims of correlated variability stem from long-term monitoring of fluxes in these bands that include high-activity states (Bartoli et al. 2011; Acciari et al. 2014), as well as observations of particular flaring events that probe correlated variability on timescales as short as 1 hr (Giebels et al. 2007; Fossati et al. 2008; Acciari et al. 2009). The detection of such a correlation in a low state was reported for the first time in Aleksić et al. (2015b), using the X-ray (Swift-XRT, RXTE-PCA) and VHE (MAGIC, VERITAS) data obtained during the 4.5-month multiwavelength campaign in 2009, when Mrk 421 did not show any flaring activity and varied around its typical Swift-XRT 0.3–10 keV count rates of ∼25 s⁻¹ and VHE flux of 0.5 Crab. In this section, we confirm the flux–flux correlation in the X-ray and VHE bands with higher confidence, during a period of even lower activity. We also study the characteristics of such a correlation in different X-ray bands using the strictly simultaneous Swift, NuSTAR, MAGIC, and VERITAS observations.¹¹² We summarize the results in Figure 10 for three nonoverlapping X-ray bands. The flux in each band was calculated from the best-fit broadband model (power-law, or log-parabolic where needed; see Section 3).

Zoom In Zoom Out Reset image size

In the following, we use the discrete correlation function (DCF) and the associated uncertainty as defined in Edelson & Krolik (1988). We carried out the correlation analysis on two timescales: a ∼1 hr timescale, using strictly simultaneous observations, and a ≲1-day timescale, using data averaged over complete 6–10 hr observations (i.e., averaged over one night of observations with each of the VHE observatories). Figure 1 shows the exact overlap of the NuSTAR and VHE observations. For both timescales and for all three X-ray bands we find significant correlations between fluxes in log–log space: DCF in all cases, with typical uncertainty of 0.1–0.2. The DCF is therefore inconsistent with zero with a minimum significance of 3.5σ (nightly averaged fluxes, 7–30 keV band) and maximum significance of 15σ (simultaneous data, 3–7 keV band). As a sanity check, we also compute Pearson's correlation coefficients and find σ significance in all cases. The strongest correlation, at 14σ, is again found for the 3–7 keV band and strictly simultaneous data shown in the upper middle panel of Figure 10. Note that we compute the correlation coefficients and DCF values using the logarithm of flux: because of greater dynamic range, the true flux–flux correlations are even more significant.

The similar slope of the correlation on both timescales may indicate that the correlation is mainly driven by flux variability on a timescale of several days, i.e., between different observations, rather than within single observations spanning several hours. The statistical significance of the correlation on the ∼daily timescale is lower, owing both to the smaller number of data points and to the fact that flux variance is larger because of the presence of strong variability on shorter timescales. For a chosen X-ray band, the best-fit slopes of the relation (a; listed in Figure 10) are statistically consistent with a single value. This is in good agreement with our finding that the dominant X-ray flux variability timescale is 9 hr (see Section 3.1). It could also be indicative of a lag between the bands, which is longer than the binning of the data taken strictly simultaneously; however, such an analysis is outside of the scope of this paper. Results of Aleksić et al. (2015b) point to absence of any lags between the X-ray and VHE bands in a nonflaring state of Mrk 421 in 2009.

An interesting result stems from our ability to broaden the search for the correlation over a very wide band in X-rays, enabled by the simultaneous Swift and NuSTAR coverage. As shown in the upper three panels of Figure 10, the slope of the relation systematically shifts from 1.00 ± 0.08 for the soft 0.3–3 keV band, to 0.80 ± 0.04 for the 3–7 keV band, and to 0.66 ± 0.05 for the hard 7–30 keV band. The same behavior is seen in the nightly averaged data, with somewhat lower significance. The persistence of this trend on both timescales and in all observatory combinations counters the possibility of a systematic bias related to those choices. We interpret it as an indication that the soft X-ray band scales more directly with the VHE flux (which is dominated by soft γ-ray photons on the low-energy end of the VHE band) owing to the emission being produced by the same population of relativistic electrons. The greater relative increase in the hard X-ray flux with respect to the soft band is consistent with both the spectral hardening already revealed by our analysis (see Section 3.4) and the fractional variability distribution determined from our data (Figure 9). Our interpretation would imply that the hard X-ray band scales more directly with the higher-energy VHE flux (e.g., TeV), which we cannot quantify well with the current data.

We find no significant correlation of the simultaneously observed spectral slopes in the hard X-ray and VHE bands. Remarkably, on two dates when the observed X-ray flux was lowest, January 10 and 20 (MJD 56302 and 56312, respectively), steep spectra with were observed by NuSTAR and both VHE observatories. Other simultaneous observations yield in the 3–30 keV band and in the VHE band, with an average photon index of approximately 3 in both bands. In comparison to the previously published results, we note that the observed steepness of the X-ray and VHE spectra presented here is atypical for Mrk 421. In a more typical low state, such as that observed in the 2009 campaign (Abdo et al. 2011), a photon index of has been observed in both bands. Here we compare our NuSTAR spectral slopes to that of the RXTE spectrum (2–20 keV) integrated over several months. Care should be taken in comparing with previously published results, because direct slope measurements in the 3–30 keV band were not available before, especially on short timescales. While the simultaneously observed steep slopes add support to the connection between X-ray and VHE bands, higher-quality data for the quiescent states are clearly needed in order to quantify it further.

4.2.2. UV/Optical versus Other Bands

Despite the low flux observed in the X-ray and γ-ray bands, the range of the UV/optical flux was higher than in some flaring episodes reported in the past (e.g., Aleksić et al. 2012). In this section we present flux correlation analyses with respect to the UV band, as represented by measurements using Swift-UVOT. The choice of band UVW1 ( Å) for this work is arbitrary; results do not change for either of the other two filters, as all of them sample the flux on the opposite side of the extremely variable synchrotron SED peak from the X-ray band. In Figure 11 we show the correlations between the UV and optical, soft X-ray and GeV γ-ray bands, each normalized to the lowest flux observed in the 2013 campaign. As in previous sections, we use the DCF and the associated uncertainty to quantify the correlation significance.

Zoom In Zoom Out Reset image size

A strong correlation is expected between the UV and optical fluxes and is confirmed by the data presented here. Previous work hinted at a possible correlation of the optical flux and the X-ray flux, but over a very narrow dynamic range and with low significance (Lichti et al. 2008). The states of Mrk 421 observed in early 2013 are not consistent with that result, indicating perhaps that a physically different regime was probed. We examine two different timescales in more detail here: for the UV and X-ray measurements taken within 6 hr of each other the DCF is 0.2 ± 0.2, while for weekly averaged values it is 0.3 ± 0.4, i.e., consistent with zero in both cases. Note that the X-ray data require averaging in the latter case, and that the uncertainty in flux is dominated by intrinsic variability. Given the established difference in the variability characteristics (see Figure 9), the lack of a significant correlation, especially on the shorter timescale, is not unexpected.

The most interesting correlation in terms of constraints on physical models is the one between the UV (and hence the optical, given their essentially 1-to-1 correspondence) and the Fermi-LAT band (0.2–100 GeV), shown in the bottom panel of Figure 11. In the framework of the SSC model, emission in these bands is due to electrons of roughly the same energy. Since the Fermi-LAT light curve had to be derived in ≈1-week bins owing to the low photon counts (see Section 2.5), we average the UV flux over the same time periods in order to cross-correlate them. Overall, the DCF is 0.8 ± 0.3, revealing a possible correlation with a 2.7σ significance. In order to examine its robustness against contributions from the outlying data points, we perform the following test: we first remove the highest-flux Fermi-LAT data point and find that it does not change the DCF, and then we remove the lowest-flux Swift-UVOT data point, which increases the DCF to 0.9 ± 0.4 but lowers the significance. We therefore estimate the correlation significance to be ≃2.5σ based on the data presented here. Existence of a real correlation between UV and γ-ray bands cannot be confirmed with the current data, but this may be possible with the Mrk 421 observations at higher γ-ray flux, such as those taken during our multiwavelength campaign in 2013 April.

Previous work on modeling the Mrk 421 SED established that the lower-energy peak of the SED, likely arising from synchrotron processes, is usually located at frequencies between ∼10¹⁷ and ∼10¹⁸ Hz. The peak itself is therefore often directly observable in the Swift-XRT band, as in other similar HBL sources (e.g., Tramacere et al. 2007a; Furniss et al. 2015). Its location in frequency space can be estimated from the UV/optical and X-ray data, using a reasonable smooth interpolation or extrapolation model (e.g., Massaro et al. 2004; Błażejowski et al. 2005; Tramacere et al. 2009; Ushio et al. 2010). For the analysis presented here, we use 30 pairs of Swift-XRT and NuSTAR spectra assembled from data taken simultaneously, together with UV data taken within the same Swift observation and R-band data taken within 24 hr. Since optical variability is significantly lower, especially on short timescales, the nonsimultaneity of the optical flux measurements is not a serious concern. We employ the two most commonly used methods from the literature to localize the synchrotron SED peak. In the top panel of Figure 12 we show examples of both methods applied to two representative sets of data. The two methods are (i) fitting a log-parabolic model to the X-ray data alone, using Xspec model logpar described by Equation (3), and extrapolating to lower energies; and (ii) fitting a log-parabolic model to both optical/UV and X-ray data.

Zoom In Zoom Out Reset image size

The X-ray-based extrapolation underestimates the UV/optical flux by more than an order of magnitude in nearly all cases. The peak frequencies predicted by this method uniformly cover the frequency range from 10¹⁶ to 10¹⁷ Hz and a factor of in peak flux. Interpretation of this simplistic parameterization of the SED would imply that it should consist of two superimposed components in order to match the observed UV/optical flux. However, log-parabolic fits that additionally include the R-band and UVOT fluxes provide a simpler solution that matches the data well on both sides of the SED peak. This is demonstrated by the two examples shown in Figure 12. Both methods are somewhat sensitive to the systematic uncertainties in the cross-normalization between the instruments and to the exact values of the line-of-sight column density and extinction corrections. We conservatively estimate that the combination of these effects results in a factor of uncertainty in the synchrotron peak frequency (), which dominates any statistical uncertainty from the fits. For this reason we do not show the uncertainties for individual estimates in Figure 12.

Both methods consistently show the peak at an atypically low frequency ( Hz), with peak flux comparable to high-activity states (see lower panel of Figure 12 and references listed there). The scatter is found to be larger for the fits using only the X-ray data, which can be easily understood since the curvature is subtle in all but the lowest-energy bins of the Swift-XRT band and the parabola has no constraint at energies below the peak. The optical/UV data provide the leverage to constrain the parabolic curves significantly better. We find an interesting trend using the second method: the flux at the SED peak is approximately proportional to the square root of the peak frequency. This is highlighted with a linear fit shown in the lower panel of Figure 12. The best-fit slope of the relation parameterized as is . The dynamic range over which the relation holds, assuming that this parameterization is valid at all, is unknown. It is clear from the comparison with the peak localization taken from the literature (listed in the lower panel of Figure 12) that the relation is not universal, although its slope is broadly consistent with the slopes of previously identified relations of the same kind. Given the fact that the simple method used to derive it is purely phenomenological and sensitive to systematic uncertainties, we refrain from quantifying and interpreting this correlation further.

For all observations presented here Hz, essentially independent of the choice of model. If only a single peak is assumed to exist, our data unambigously imply that the synchrotron peak frequency was well below the Swift-XRT band (<0.3 keV) during this period of very low X-ray activity. This is atypical for Mrk 421 based on the previously published data. We found that, for the states of lowest X-ray activity, our data show a peak frequency as low as Hz, estimated with a (symmetric) log-parabolic function fitted to both optical/UV and X-ray (Swift-XRT and NuSTAR) data. A peak frequency as low as ∼10¹⁶ Hz has been reported in Ushio et al. (2010), but only when using a specific model that led to a very asymmetric parameterization of the synchrotron bump. When using a symmetric function (comparable to a log-parabola), the fit to the same observational data reported in Ushio et al. (2010) led to a peak frequency of ∼10¹⁷ Hz, which is typical for Mrk 421.

Blazars are classified by the frequency of the peak of the synchrotron emission as LBLs, IBLs, and HBLs (for low-, intermediate-, and high- BL Lacs) if is below 10¹⁴ Hz, in the range 10¹⁴–10¹⁵ Hz, or above 10¹⁵ Hz, respectively. The data presented here show that Mrk 421, which is one of the archetypal TeV HBLs, with a synchrotron peak position well above that of the typical HBL, changed its broadband emission in such a way that it almost became an IBL. This effect may also happen to other HBLs that have not been as extensively observed as Mrk 421 and suggests that the SED classification may denote a temporary characteristic of blazars, rather than a permanent one.

For a better understanding of the empirically observed correlations, we need to consider the complete broadband spectrum. SED snapshots for four selected epochs (marked with dashed vertical lines in Figure 7) are shown in Figure 13. They were selected to show a state of exceptionally low X-ray and VHE flux (January 10 and 20; see Section 4.2.1), in contrast to higher, though not flaring, states (January 15, February 12, for example). In all SED plots we also show data accumulated over 4.5 months in the 2009 multiwavelength campaign (Abdo et al. 2011), which is currently the best-characterized quiescent broadband SED available for Mrk 421 in the literature. For the two epochs of very low X-ray flux, we show for the first time states in which both the synchrotron and inverse-Compton SED peaks are shifted to lower energies by almost an order of magnitude compared to the typical quiescent SED. The accessibility of the low-activity state shows the large scientific potential brought by the improvements in the X-ray and VHE instrumentation in the past several years with the launch of NuSTAR and upgrades to the MAGIC and VERITAS telescopes.

Zoom In Zoom Out Reset image size

We note that the empirical SEDs from the 2013 campaign shown here represent hr of observation in the X-ray and VHE bands, rather than integrations over a time period of weeks or even months. We match the simultaneous UV, X-ray, and VHE data to optical data taken within at most 2 days, radio data taken within at most 2 weeks, and Fermi data integrated over time intervals of 6–10 days centered on the time of the coordinated X-ray and VHE observations. Mrk 421 is a point-like and unresolved source for single-dish radio instruments, which means that the data shown in Figure 13 include emission from spatial scales larger than the jet itself and therefore should be considered as upper limits for the SSC models of jet emission. We further discuss the SED in the context of the SSC model in Sections 5.3 and 5.4.

In addition to the coordinated multiwavelength campaign conducted in 2013, Mrk 421 was observed independently with several instruments in 2012 July–September, including NuSTAR, Fermi, and OVRO. In 2012 July, the flux in the Fermi-LAT band increased above the median level and peaked twice over the following 2 months (see Figure 14). The first peak was reported on July 16 by the Fermi-LAT Collaboration (D'Ammando et al. 2012) and, within the same day, by the ARGO-YBJ collaboration (Bartoli et al. 2012). The daily flux seen by Fermi-LAT increased to s⁻¹cm⁻², a factor of above the average flux reported in the second Fermi-LAT catalog (2FGL; Nolan et al. 2012). Light curves from several observatories monitoring Mrk 421 in 2012 July–September are shown in Figure 14 in order to provide a timeline for this flaring event.

Zoom In Zoom Out Reset image size

An observation of Mrk 421 was performed by NuSTAR on 2012 July 7 and 8, shortly before the start of flaring activity in the γ-ray band. The observation¹¹³ was not originally intended for scientific usage, as the pointing was suboptimal at this early point in the mission (less than a month after launch). However, it represents both the longest and the most variable NuSTAR observation considered in this paper (see Figure 1) and thus represents an important part of the NuSTAR data presented here. The available X-ray and γ-ray data are clearly too sparse to allow for associations to be inferred between any specific features in the light curves. There is indication from the MAXI public monitoring data¹¹⁴ that the X-ray flux in the 4–20 keV band increased further after the NuSTAR observation, peaking between 2 and 5 weeks later (see Figure 14). The NuSTAR observation therefore makes it possible to investigate the hard X-ray spectrum of Mrk 421 at the time when its γ-ray activity was rapidly increasing.

A unique feature of this event is the well-defined rapid radio flare observed at 15 GHz from OVRO (Hovatta et al. 2012, 2015). Approximately 70 days after the first peak of the γ-ray flare, Mrk 421 reached a flux density of Jy, approximately 2.5 times its median flux density. Note that in Figure 14 we show the light curves on a logarithmic scale with a fixed dynamic range in order to facilitate direct comparison with Figure 7; on a linear plot both the γ-ray and the radio flare appear strikingly peaked. Based on the statistical properties of Fermi-LAT and OVRO 15 GHz light curves, Max-Moerbeck et al. (2014) have shown that the γ-ray and the radio flare are likely causally related. During most of the flaring activity in 2012, Mrk 421 was very close to the Sun on the sky, which resulted in relatively poor multiwavelength coverage. We therefore do not attempt a more comprehensive analysis of the sparse data available for this epoch. Hovatta et al. (2015) present the radio and γ-ray data, as well as a physical model for the flaring activity, for both the 2012 flare and the flare observed during the multiwavelength campaign in 2013 April. As the latter event has been covered with numerous instruments (e.g., Baloković et al. 2013b; Cortina et al. 2013; Paneque et al. 2013; Pian et al. 2014), we defer a comprehensive analysis of this period to a separate publication.

The good temporal coverage of the NuSTAR data reveals a typical variability timescale of hr. Significant variability is clearly detectable even in the low-flux states, which is the case for several epochs in early 2013. We find no evidence for strong intrahour variability; on timescales as short as ∼10 minutes the variability amplitude is ≲5%, approximately an order of magnitude lower than the typical flux change over a 10 hr observation. This can be inferred from the fact that in the light curves shown in Figure 1 only a small fraction of adjacent bins differ in count rate by more than 3σ. We summarize this more formally in Figure 2. In contrast to some previous studies (e.g., Takahashi et al. 1996), we observed no clear time-dependent circular patterns in the count rate–hardness plane. The reason for this may be that most of the observations seem to have covered periods of decreasing flux, with no well-defined flare-like events except for a few "mini-flares" of modest % amplitude.

We find that the X-ray spectrum of Mrk 421 cannot generally be described as a simple power law, but instead it gradually steepens between 0.3 and ∼70 keV. For most of the Swift and NuSTAR observations in the 2013 January–March period, we find that the spectra in both bands are better described when a curvature term is added to the basic power law, as in the log-parabolic model available in Xspec. Using this model, we find significant curvature at the highest observed fluxes—still notably lower than in any flaring states—gradually vanishing as flux decreases (see Figure 5). This has a simple explanation because the X-ray band samples the Mrk 421 SED close to the synchrotron peak: when the X-ray flux is low and the SED peak shifts to lower energy (away from the NuSTAR band), the hard X-ray spectra can be well described by a power law. This behavior is consistent with the steady increase of fractional variability with energy through the X-ray band, as shown in Figure 9. The high sensitivity of NuSTAR reveals that the hard X-ray spectrum does not exhibit an exponential cutoff, and it is well described by a power law with a photon index , even during the epochs related to the lowest X-ray fluxes. The NuSTAR data also show no signature of spectral hardening up to ∼80 keV, meaning that the onset of the inverse-Compton bump must be at even higher energies.

The data gathered in the 2013 multiwavelength campaign contribute some unique details to the rich library of blazar phenomena revealed by Mrk 421. The object is highly variable on a wide range of timescales and fluxes, with the fractional variability amplitude highest at the high-energy ends of the synchrotron and inverse-Compton SED bumps (see Figure 9). The well-matched coverage in the X-ray and VHE bands reveals that the steep spectral slope observed in the X-ray band at very low flux occurs simultaneously in time with an atypically steep slope observed in the VHE band. For the first time we observed a simultaneous shift of both the synchrotron and the inverse-Compton SED peaks to lower energies in comparison to the typical quiescent state (see Figure 13), constrained primarily by the X-ray and VHE data. The measurements in those bands do not support the existence of high-energy cutoffs up to and TeV. All of this indicates that the energies of radiating particles must be very high (up to see Section 5.3 below) even when the source is in such a low state.

In Section 4.2.1 we have shown that the X-ray and VHE fluxes are correlated at >3σ significance. Parameterizing the correlation as , the correlation is found to be approximately linear ( 1) on both half-hour and half-day timescales. This is consistent with most previous results considering similar spectral bands: (Tanihata et al. 2004), (Fossati et al. 2008, for averaged and nonflaring periods; also Aleksić et al. 2015b), a = 1 (Acciari et al. 2014; assumed linear), etc. We emphasize the importance of distinguishing between (i) a correlation of count rates versus a correlation of fluxes, since the conversion between them is nonlinear owing to spectral variability, and (ii) a general correlation versus a correlation associated with isolated flares, since those could potentially be produced by different physical mechanisms (as argued by, e.g., Katarzyński et al. 2005). Indeed, for isolated flaring periods Fossati et al. (2008) and Giebels et al. (2007) find and , respectively. Care should be taken with direct comparison of the results in the literature, since the chosen spectral bands are not always the same, and we have shown in Section 4.2.1 that the slope does depend on the band choice as a consequence of the spectral variability.

We note that in the simplest, one-zone SSC model, one expects a close correlation between the X-ray and VHE fluxes. However, if the scattering takes place exclusively in the Thomson regime, the inverse-Compton flux should obey a quadratic (a = 2) relationship, since increasing both the number of electrons and the seed photon flux results in a quadratic increase in the scattering rate. Since we detect a linear relationship, this would argue that the scattering cross section is diminished, possibly because the scattering takes place in the less efficient Klein–Nishina regime. For example, a quadratic relation has recently been observed in a similar HBL object, Mrk 501 (Furniss et al. 2015). It has been shown previously (e.g., Katarzyński et al. 2005) that this implication is valid only if the normalization of the entire electron distribution is changed to produce flux variations. For changes in other parameters of the electron distribution, or in physical conditions within the emission region, this is no longer strictly correct. The linearity of the flux–flux correlation itself does not uniquely indicate Klein–Nishina effects; we therefore combine the broadband SED modeling and variability properties in the following section, in order to further investigate this issue.

In the framework of an SSC model, if the peak energies of the synchrotron and inverse-Compton SED are resolved, then along with constraints from temporal variability and an estimate of the bulk Doppler factor of the emitting material, fairly general and robust estimates can be made of the characteristic particle energies, the magnetic field strength, and the overall size of the emitting region. An estimate for the characteristic electron Lorentz factor (, measured in the comoving frame of the emitting plasma) is given roughly by the square root of the ratio of the energies of the synchrotron and inverse-Compton peaks. The radiation at the peaks is dominated by electrons of intermediate energy, where the Klein–Nishina reduction in the scattering cross section is not expected to be significant. From the SEDs shown in Figure 13, we estimate Hz and Hz. It then follows that . Assuming a bulk Doppler factor of , we can also estimate the magnetic field strength in the plasma comoving frame: G. Finally, an upper limit for the size of the emitting region in the comoving frame is given by the observed variability timescale, , and the bulk Doppler factor δ: cm.

We can constrain the properties of the jet emission region more precisely by directly modeling the multiwavelength SED data shown in Figure 13 with a standard one-zone SSC model. Specifically, we apply an equilibrium version of the SSC model from Böttcher & Chiang (2002), fully described in Böttcher et al. (2013). This model has already been used to represent Mrk 421 in two different states (Acciari et al. 2009). In this model, the emission originates from a spherical region with radius R, containing relativistic electrons that propagate down the jet with a bulk Lorentz factor Γ. In order to decrease the number of free parameters, we assume a value with the jet axis aligned near the line of sight with the critical angle rad , which makes the Doppler factor equal to the jet Lorentz factor (). This simplifying choice is often used in the literature when direct measurements are not available (see, e.g., Abdo et al. 2011 and the discussion therein). A Doppler factor of 25 is higher than the value inferred from VLBA measurements of the blob movement by Piner et al. (2010). This is a common situation in VHE blazars, often referred to as the "bulk Lorentz factor crisis," and requires that the radio and VHE emissions are produced in regions with different Lorentz factors (Georganopoulos & Kazanas 2003; Ghisellini et al. 2005; Henri & Saugé 2006). High Doppler factors (≳10) are required to explain previously reported rapid variations in the VHE band (Gaidos et al. 1996; Celotti et al. 1998; Galante 2011) and are typically used in theoretical scenarios to describe the broadband emission of VHE blazars. Relativistic leptons are injected according to a power-law distribution between and . These particles lose energy through synchrotron and inverse-Compton radiation, leading to an equilibrium between particle injection, radiative cooling, and particle escape. The particle escape is characterized with an escape efficiency factor η, defined so that is the escape time. This results in a particle distribution that propagates along the jet with power L_e. Synchrotron emission results from the interaction of particles with a magnetic field B, generating a Poynting-flux luminosity of L_B. L_e and L_B allow the calculation of the equipartition parameter L_B/L_e. Various other blazars have been represented with this model, with the resulting model parameters summarized in Aliu et al. (2013). In application to the broadband data, the intrinsic source VHE flux from the SSC model is absorbed by the Franceschini et al. (2008) model describing the extragalactic photon field. In Table 9, we list the relevant model parameters that reproduce the observed SED of Mrk 421 for the four selected epochs in 2013.

Table 9.  Model Parameters for the Equilibrium SSC Model for Four Selected Epochs

Parameter	Jan 10	Jan 15	Jan 20	Feb 12
(104)	2.2	3.5	2.4	5.0
(105)	4.0	4.8	5.8	6.8
q	3.5	3.5	3.6	3.0
η	35	35	35	35
B (G)	0.17	0.25	0.16	0.10
Γ	25	25	25	25
R (1016 cm)	0.9	0.55	1.0	1.4
θ (deg)	2.29	2.29	2.29	2.29
τ (hr)	3.4	2.1	3.8	5.3
Le (1043 erg s−1)	3.0	2.6	3.0	3.3
	0.18	0.25	0.20	0.14

Note. The electron energy distribution parameters listed here refer to the electrons injected into the emission region, and the equilibrium distribution is calculated self-consistently within the model, as described in Section 5.3. Model SED curves are shown in Figure 13.

Download table as:  ASCIITypeset image

Since the injected particle distribution in our SSC model follows a single power law, the observed spectral shapes in the GeV and VHE bands imply certain constraints on the model parameters. In the 0.1–100 GeV band, the observed spectra have photon indices in the range Γ ∼ 1.6–1.7, while, by contrast, the VHE spectra have photon indices of Γ ∼ 2.3–3.5 (see Tables 5 and 6). These indices imply spectral breaks of 0.6–1.9, which are moderately to significantly larger than the "cooling" break of that arises from incomplete (or "weak") synchrotron cooling of an injected power-law distribution of electrons. In the strong-cooling regime, i.e., where the synchrotron cooling timescale is shorter than the particle-escape time, the cooled electron distribution has a break at the lower bound of the injected power law, , and has power-law shapes for and for . For the parameters shown in Table 9, this particle distribution implies a synchrotron spectrum with a peak at ≃ (1– Hz and spectral shapes () for and () for (see Table 7 and Figure 13). For and B = 0.2 G, a synchrotron cooling timescale of is obtained in the comoving frame of the emitting plasma; this is slightly larger than the nominal (i.e., in the absence of any scattering) escape time of . The escape efficiency factor, , ensures that the cooled electron distribution extends to sufficiently low energies to model both the optical/UV points and the Fermi-LAT data down to 0.1 GeV. The Larmor radius of the lowest-energy electrons in the modeled magnetic field is small enough that the electrons have sufficient time to cool within the emission region before escaping.

Past SED modeling of HBL-type blazars has often used SSC calculations that have electron distributions assumed to persist in the specified state for the entire duration of the observation. For example, for a given variability timescale, a single, time-averaged, multiply broken power-law electron distribution is used by Abdo et al. (2011) to model the multiwavelength data obtained for Mrk 421 over a 4.5-month period in early 2009. By contrast, the SED calculations we have performed in this paper attempt to model specific flaring or quiescent periods for which most of the data (optical, X-ray, and VHE) were obtained within 12 hr intervals. Our modeling assumes that an initial power-law electron spectrum is injected into the emission region, and we compute the resulting quasi-equilibrium particle distribution for those epochs given the radiative and particle escape timescales. Since the 2009 observations could contain a large number of similarly short flaring and quiescent episodes with a range of physical properties, it would be inappropriate to attempt to model those data with the procedure we have used here. However, as we have indicated above, it should be possible to obtain equivalent time-averaged SED models that have multiply broken power-law electron distributions and that we can compare directly to the Abdo et al. (2011) results.

We have performed such modeling, and in Table 10, we give the parameters for the same four selected epochs appearing in Table 9. The SEDs produced by the two models can be matched very well, as shown in Figure 13 with the blue and dashed red lines. As we noted above, the equivalent time-averaged electron distributions can be represented via a broken power law with a break at and index p_l = 2.0 below and index above the break.

Table 10.  Model Parameters for the Snapshot SSC Model for Four Selected Epochs

Parameter	Jan 10	Jan 15	Jan 20	Feb 12
(103)	2.0	1.3	1.7	3.0
(104)	2.5	3.7	3.0	5.2
(105)	4.0	4.8	5.8	6.8
pl	2.0	2.0	2.0	2.0
ph	4.5	4.5	4.6	4.0
B (G)	0.21	0.28	0.19	0.10
Γ	25	25	25	25
R (1016 cm)	0.93	0.60	1.04	1.69
θ (deg)	2.29	2.29	2.29	2.29
τ (hr)	3.4	2.2	3.9	6.3
Le (1043 erg s−1)	3.3	2.9	3.3	4.1
	0.51	0.47	0.56	0.33

Note. The electron energy distribution parameters listed here refer to the distribution directly responsible for the SSC emission. This simplified model is described in Section 5.3 and used for comparison with the literature. Model SED curves are shown in Figure 13.

Download table as:  ASCIITypeset image

As Abdo et al. (2011) demonstrate and as we discuss above, for data that resolve the shapes of both the synchrotron and self-Compton components of the SED, the model parameters in these sorts of leptonic time-averged models are largely determined once either the variability timescale or the Doppler factor is constrained or set. Therefore, when comparing the current model parameters to those of Abdo et al. (2011), we consider just their results. In terms of the shape of the underlying particle distributions, the value of p_l = 2.0 we find is comparable to their value of p₁ = 2.2, and our values of p_h = 4.0–4.6 are similar to their high-energy index of p₃ = 4.7. While this does not uniquely imply that the same energy-loss mechanisms and acceleration processes are at work in both cases, the consistency is encouraging. The Abdo et al. (2011) modeling does require an additional medium-energy power-law component that is dictated by their generally broader SED peaks (see the gray points in Figure 13). In the context of the quasi-equilibrium modeling, this would arise from a distribution of physical parameters in shorter flaring and quiescent episodes that are averaged over the 4.5-month observation time.

Several physical parameters in the current modeling do differ substantially from those of Abdo et al. (2011). The characteristic electron Lorentz factors are about an order of magnitude lower (– vs. ), while the inferred emitting region radius is about a factor of 3–10 smaller (R = (0.6– cm vs. cm), the inferred magnetic field is substantially higher (B = 0.10–0.28 vs. 0.038), and the resulting jet powers differ by a factor of three to four. As a consequence, the current modeling yields equipartition parameters that are much closer to unity, in the range –0.56 vs. for the Abdo et al. (2011) result. Given the overall similarity in the size and shape of the synchrotron and SSC components among all five data sets (i.e., the four 2013 epochs and the 2009 data shown in Figure 13), the differences in model parameters can be understood as being driven mostly by the combination of the order-of-magnitude larger characteristic electron Lorentz factor and the order-of-magnitude higher peak synchrotron frequency required by the Abdo et al. (2011) data and modeling. Since , the order-of-magnitude higher magnetic fields in the current modeling are readily understood, and those in turn largely account for the equipartition parameters being substantially closer to unity. We note that other authors, such as Aleksić et al. (2015c), have inferred SSC model parameters that are below equipartition by much more than an order of magnitude. However, Aleksić et al. (2015c) considered a flaring state of Mrk 421 that had much higher synchroton peak frequencies, as well as substanially higher fluxes at all wavelengths. Accordingly, their much larger disparity in the partitioning of the jet power compared to the current results is not surprising and roughly fits in with the preceding discussion.

Studying broadband emission of Mrk 421 at different epochs, Mankuzhiyil et al. (2011) found that there were no substantial shifts in the location of the peaks of the synchrotron and the inverse-Compton bumps. They concluded that the variability in the blazar emission was dominated by changes in the parameters related to the environment, namely, the emission-region size, the Lorentz (Doppler) factor, and the magnetic field. The observational results presented here, with substantially broader energy coverage and better instrumental sensitivity due to the advent of new γ-ray and X-ray instruments, differ from those presented in Mankuzhiyil et al. (2011). We show that, besides changes in the magnetic field, the distortions in the broadband emission of Mrk 421 also require changes in the electron energy distribution, which may be due to variations in the mechanism accelerating the electrons to high energies.

Having modeled the broadband SEDs with single-zone SSC calculations, we can test the hypothesis that the VHE emission occurs in the Klein–Nishina regime. The SED modeling yields injected electron Lorentz factors in the range to . Assuming that the target synchrotron photons for inverse-Compton scattering have energies around the synchrotron peak at Hz, the parameter governing the transition between Thomson and Klein-Nishina regimes is (Blumenthal & Gould 1970), which in the observer frame becomes .

When considering photons from the synchrotron peak position ( 40 eV, i.e., about one order of magnitude lower than the typical position of the synchrotron peak in Mrk 421), we obtain 0.4–8, indicating that the inverse-Compton scattering of photons with energy takes place, at least partially, in the Klein–Nishina regime. The X-ray energies probed with Swift-XRT are roughly one order of magnitude above , far above the range where Thomson scattering is relevant, and consistent with the linear () relationship between the soft X-ray and VHE flux.

The electrons responsible for the broadband emission of Mrk 421 lose energy mostly owing to synchrotron cooling, as one can infer from the dominance of the synchrotron bump over the inverse-Compton bump shown in the SEDs from Figure 13. Note that the inefficiency of cooling via the Compton channel is independently implied from the observed slope of the X-ray–VHE flux correlation. The observed variability timescale (measured in a stationary observer's frame) due to synchrotron cooling alone is given by s, where E is the photon energy in keV and B is the comoving frame magnetic field strength in G. Taking as the energy typical for the NuSTAR band, assuming G, as found from our SED modeling, and as before, we arrive at of s.¹¹⁵ This is more than an order of magnitude shorter than the variability timescale hr measured in the observer's frame, as we can derive from the NuSTAR light curves. Since the synchrotron cooling timescales are so short, this requires that the electron acceleration must be happening locally, very close to where the emission takes place.

Considering the disparity between the variability timescale and the synchrotron cooling timescale, along with the similarity of the increases and decreases in flux during the NuSTAR observations (Figure 1), it seems unlikely that the output is dominated by a single shocked region as a site of particle acceleration, such as is often argued to be the case in flaring episodes. Instead, we can interpret the flux changes as a geometrical effect due to a spatially extended region containing multiple particle-acceleration zones contributing comparably to the overall SED. Observation of variability due to geometrical effects of a spatially extended region would lack sharp flux increases in the X-ray band, which might result from sudden particle-acceleration events, because the sharp flux increases and decreases from the different regions (even if partially connected) would not occur at exactly the same time. In this scenario, the shortest variability timescales, comparable to the electron cooling timescales, would be produced only when a single region dominates the overall emission, which is expected to occur during flaring episodes, but not during the relatively low activity reported in this paper. As described in Section 3.1, the observed increases appear at least as smooth and as slow as the observed decreases, consistent with this picture.

One may argue that the X-ray flux variability reported in Section 3.1 is not due to the acceleration/cooling of electrons, but rather produced by variations in the parameters related to the environment (e.g., B, R) or the Dopper factor δ (e.g., due to a change in the viewing angle). In that case, the smooth and relatively slow changes observed in the NuSTAR light curves would not be related to the short electron cooling timescales derived above, but rather to the variations in the above-mentioned parameters. However, such a theoretical scenario is strongly disfavored by the fractional variability as a function of energy reported in Figure 9, as well as by the lack of correlation between optical and X-ray fluxes reported in Figure 11, while there is a correlation between optical and GeV γ-ray fluxes, as well as X-ray and VHE fluxes, reported in Figures 10 and 11, respectively. The only possibility for the parameters R, B, or δ to dominate the measured flux variations would be to have, at least, two distinct emission regions, one dominating the optical and GeV γ-ray bands, and the other dominating the X-ray and VHE bands. Therefore, despite the success of the one-zone SSC scenario in describing the broadband SED (see Section 5.3), we argue that the observed multiwavelength variability and correlations point toward an emission region composed of several distinct zones and dominated by changes in the electron energy distribution. The increase in the fractional variability with energy for both SED bumps and the harder-when-brighter trend that is clearly observed in the X-ray spectra measured with NuSTAR (which is the segment of the broadband SED reconstructed with the highest accuracy) indicate that the changes in the electron energy distribution are generally chromatic,¹¹⁶ with strongest variability in the highest-energy electrons. However, the saturation of the X-ray spectral shape at the lowest and highest X-ray fluxes (see Section 3.4 and Figure 6) suggests that at the times of lowest and highest activity, the variations in the electron energy distribution become achromatic, at least for those electron energies responsible for the X-ray emission. It is possible that at those times the variability is not dominated by acceleration and cooling of the electrons, but rather by variations in the physical parameters of the environment in which particle acceleration occurs. For the periods of very low activity, a possibility would be that the radiation is being produced within a larger region by particles accelerated by Fermi II processes (e.g., stochastic acceleration on magnetic turbulence), as suggested, for instance, by Massaro et al. (2004) and Ushio et al. (2009).

The magnetic field implies a size constraint for the acceleration zones, since electrons cannot attain energies corresponding to a gyroradius significantly larger than the characteristic size of a zone. The NuSTAR data imply no cutoff in the synchrotron SED up to ∼80 keV, so we can estimate the electron gyroradius corresponding to that photon energy using B = 0.2 G and the maximal . Since , we have cm, which is much smaller than the inferred emission-region size of 10¹⁶ cm. Given the large difference of five orders of magnitude between the gyroradius for the highest-energy electrons and the size of the overall emitting region, the electrons cannot travel far from their acceleration site without losing a substantial fraction of their energy, and hence the particle acceleration and the emission need to be essentially co-spatial. We therefore conclude that the set of physical parameters discussed here offers a self-consistent picture in which the observed properties of Mrk 421 in a nonflaring state are consistent with compact zones of particle acceleration distributed within a significantly larger volume that produces the total emission. While detailed characterization of the acceleration process is outside the scope of the paper, one possible scenario involves magnetic reconnection and "mini-jets" formed within a larger emission volume, as suggested, for instance, by Nalewajko et al. (2011) and developed further by Nalewajko et al. (2015). For a recent summary of arguments in favor of magnetic reconnection for powering blazar jets, the reader is referred to Sironi et al. (2015).

Regardless of the exact acceleration mechanism, emitting regions composed of multiple zones, e.g, as in the model proposed by Marscher (2014), would be consistent with other behavior observed in blazars, such as the increase in the degree of polarization of the synchrotron radiation when the polarization electric vector rotates, or curvature in the SED arising from nonuniform particle acceleration and energy losses. In a low-activity state, where no single zone dominates the output, the addition of polarization vectors from individual zones would result in a low overall level of polarization with random fluctuations in both the polarization degree and angle. Our optical polarization measurements, shown in Figure 8, are consistent with that prediction. While multizone scenarios have previously been considered for flaring states (e.g., Massaro et al. 2004; Ushio et al. 2009; Cao & Wang 2013; Aleksić et al. 2015c), it has usually been assumed that the quiescent state can be well described by a simpler single-zone SSC model (e.g., Abdo et al. 2011). The observations presented here, however, show that, even in this state of very low activity, the emission region may have a more complex structure than previously assumed.