Multiwavelength Study of Quiescent States of Mrk 421 with Unprecedented Hard X-Ray Coverage Provided by NuSTAR in 2013

We present coordinated multiwavelength observations of the bright, nearby BL Lac object Mrk 421 taken in 2013 January-March, involving GASP-WEBT, Swift, NuSTAR, Fermi-LAT, MAGIC, VERITAS, and other collaborations and instruments, providing data from radio to very-high-energy (VHE) gamma-ray bands. NuSTAR yielded previously unattainable sensitivity in the 3-79 keV range, revealing that the spectrum softens when the source is dimmer until the X-ray spectral shape saturates into a steep power law with a photon index of approximately 3, with no evidence for an exponential cutoff or additional hard components up to about 80 keV. For the first time, we observed both the synchrotron and the inverse-Compton peaks of the spectral energy distribution (SED) simultaneously shifted to frequencies below the typical quiescent state by an order of magnitude. The fractional variability as a function of photon energy shows a double-bump structure which relates to the two bumps of the broadband SED. In each bump, the variability increases with energy which, in the framework of the synchrotron self-Compton model, implies that the electrons with higher energies are more variable. The measured multi-band variability, the significant X-ray-to-VHE correlation down to some of the lowest fluxes ever observed in both bands, the lack of correlation between optical/UV and X-ray flux, the low degree of polarization and its significant (random) variations, the short estimated electron cooling time, and the significantly longer variability timescale observed in the NuSTAR light curves point toward in-situ electron acceleration, and suggest that there are multiple compact regions contributing to the broadband emission of Mrk 421 during low-activity states.

variability increases with energy which, in the framework of the synchrotron self-Compton model, implies that the electrons with higher energies are more variable.The measured multi-band variability, the significant X-ray-to-VHE correlation down to some of the lowest fluxes ever observed in both bands, the lack of correlation between optical/UV and X-ray flux, the low degree of polarization and its significant (random) variations, the short estimated electron cooling time, and the significantly longer variability timescale observed in the NuSTAR light curves point toward in-situ electron acceleration, and suggest that there are multiple compact regions contributing to the broadband emission of Mrk 421 during low-activity states.Subject headings: galaxies: active -BL Lacertae objects: individual (Markarian 421) -X-rays: galaxies -gamma rays: general -radiation mechanisms: non-thermal

INTRODUCTION
Markarian 421 (Mrk 421 hereafter) is a nearby active galaxy with a featureless optical spectrum devoid of prominent emission or absorption lines, with strongly polarized variable optical and radio flux, and compact (milli-arcsecond-scale) radio emission.As such, it is classified as a BL Lacertae type.Its spectral energy distribution (SED) is well described by a characteristic twopeak shape (for a review, see, e.g., Urry & Padovani 1995;Ulrich et al. 1997).In the more general context of blazars, Mrk 421 belongs to a subclass of the so-called high-energy-peaked BL Lacertae (HBL) objects, relatively low-luminosity sources with both peaks located at relatively high energies (respectively at ∼1 keV and ∼100 GeV).Mrk 421 is among the closest and most comprehensively studied objects of the HBL class and is also the first extragalactic source detected in the very high energy γ-ray band (E >100 GeV, VHE hereafter; Punch et al. 1992).
The observed properties of Mrk 421, as well as other similar blazars, are best explained as arising from a relativistic jet seen at a small angle to our line of sight (Urry & Padovani 1995).The nonthermal and polarized continuum observed from the radio band to the soft X-ray band suggests that this part of the SED is due to a distribution of relativistic electrons radiating via the synchrotron process.The radiation in the γ-ray band is likely due to inverse-Compton scattering by energetic electrons responsible for the synchrotron radiation, as confirmed by simultaneous, correlated variations in the low-and high-energy SED components (e.g., Giebels et al. 2007;Fossati et al. 2008;Aleksić et al. 2015b).The seed photons are most likely the synchrotron photons internal to the jet.Such "synchrotron self-Compton" (SSC) models, developed by many authors (for early examples see, e.g., Jones et al. 1974;Ghisellini et al. 1985;Marscher & Gear 1985) have been successfully invoked to describe the full SED of HBL objects (e.g., Ulrich et al. 1997;Fossati et al. 2008;Tavecchio et al. 2010).
The range of measured flux of Mrk 421 spans up to two orders of magnitude in some spectral bands, with flares occurring on very short timescales (a day or less; e.g., Gaidos et al. 1996;Tanihata et al. 2003;Fossati et al. 2008).Possibly the best bands to study such variability are the X-ray and VHE γ-ray bands: in the context of the SSC model, they represent radiation from the most energetic electrons, which have the shortest timescales for radiative losses.The cross-correlation of time series measured in various bands provides additional clues to the radiative processes, the acceleration and energy distribution of radiating particles, and the structure and intrinsic power of the relativistic jet.The relative temporal variability in different spectral bands, from radio through VHE γ-rays, provides an important handle on the location of the energy dissipation with respect to the central black hole (Sikora et al. 2009;Janiak et al. 2012).In the context of a specific model for the radiation, the underlying particle distributions may be determined more or less directly from the observed multiwavelength SEDs.Particle-acceleration mechanisms can then be constrained by the shape of those particle distributions.Diffusive shock acceleration, an example of a first-order Fermi (Fermi I) process, is generally associated with single power-law distributions (e.g., Blandford & Eichler 1987;Jones & Ellison 1991).In contrast, log-parabolic distributions are produced in models of stochastic acceleration (e.g., Massaro et al. 2004;Tramacere et al. 2011), which can be considered equivalent to a secondorder Fermi (Fermi II) process.
Mrk 421 and other HBL-type blazars have been extensively studied in the soft X-ray band (e.g., Makino et al. 1987;Takahashi et al. 1996;Ravasio et al. 2004;Tramacere et al. 2007bTramacere et al. , 2009)), revealing a range of spectral slopes in various quiescent and flaring states.Less is known about the hard X-ray ( 10 keV) properties of blazar jet emission: the data are far fewer and available mostly for flaring episodes, or averaged over relatively long timescales (e.g., Guainazzi et al. 1999;Giebels et al. 2007;Fossati et al. 2008;Ushio et al. 2009;Abdo et al. 2011).This energy band probes the most energetic and fastest varying tail of the distribution of synchrotron-radiating particles, and therefore represents an important diagnostic of the content of the jet and the processes responsible for the acceleration of particles to the highest energies.The inverse-Compton component increases with energy, and could potentially contribute significantly to the hard X-ray band.If so, it would also provide a strong constraint on the low-energy part of the electron distribution, which is a significant, if not dominant, part of the total kinetic power of the jet.
Mrk 421 observations were part of the Nuclear Spectroscopic Telescope Array (NuSTAR; Harrison et al. 2013) blazar program, aimed at advancing our understanding of astrophysical jets.The multiwavelength campaign focused on Mrk 421 was carried out between December 2012 and May 2013, with three to four pointings per month, designed to maximize strictly simultaneous overlap with observations by the VHE γ-ray facilities VERI-TAS and MAGIC.We also secured nearly simultaneous soft X-ray, optical and UV observations from the Swift satellite.The γ-ray data from Fermi -LAT which observes Mrk 421 every 3 hours, was also used together with all the coordinated multiwavelength data.Mrk 421 varied in flux throughout the campaign, with a relatively low flux in the X-ray and VHE bands at the beginning, increasing to a major flare toward the end of the campaign.In this paper, we present part of the data collected during the first three months of the campaign with particular emphasis on the detailed shape of the X-ray spectrum, its variability and the correlated variability observed in VHE γ-rays.We also report briefly on the observations of Mrk 421 prior to the start of the campaign, in July 2012, when the object was used for calibration purposes during the in-orbit verification phase of NuSTAR.During this period, Mrk 421 emission was broadly consistent with previously observed quiescent states, which we define here to be characterized by relatively low flux at all frequencies and by the absence of significant flaring (see, e.g., Abdo et al. 2011).The flaring period of the 2013 campaign will be covered in a separate publication.
The outline of the paper is as follows.In § 2 we describe the multiwavelength observations and data used in this paper.We dedicate § 3 to a detailed characterization of the hard X-ray spectrum of Mrk 421 with NuSTAR.The results of the multiwavelength campaign in 2013 January-March are presented in § 4. Discussion of the empirical results and modeling of the broadband properties are given in § 5, and in § 6 we summarize the main results.We adopt a distance of 141 Mpc to Mrk 421, calculated from its measured redshift z = 0.0308 (based on absorption lines in the spectrum of the host galaxy; Ulrich et al. 1975) and the cosmological parameters recently refined by the Planck Collaboration (Ade et al. 2014): h 0 = 0.67, Ω Λ = 0.685.

Radio
The Owens Valley Radio Observatory (OVRO) 40meter 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;Angelakis et al. 2010;Fuhrmann et al. 2014), a γ-ray blazar monitoring program related to the Fermi Gamma-ray Space Telescope.Observations with the Effelsberg 100-meter and Pico Veleta 30-meter 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) and Angelakis et al. (2015).The Metsähovi Radio Observatory 14-meter 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).

Optical
The coverage at optical frequencies was provided by various telescopes around the world within the GASP-WEBT program (e.g., Villata et al. 2008, Villata et al. 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 et al. (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 7.5 ′′ 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 4000-7550 Å band spectropolarimetry 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); Jorstad et al. (2010).

Swift UVOT and XRT
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 E(B − V ) = 0.014 (Schlegel et al. 1998;Schlafly et al. 2011) with A λ /E(B − V ) 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 56358.3190 1 0.9 0.6 27.0 ± 0.9 27.9 ± 0.9 25.3 ± 0.9 25.5 ± 0. 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 logparabolic models, identical to the modeling presented in detail in § 3.2.The models include photoelectric absorption by a fixed column density estimated to be N H = 1.92 × 10 20 cm −2 (Kalberla et al. 2005).The logparabolic 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.
2.4.NuSTAR NuSTAR (Nuclear Spectroscopic Telescope Array; 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 NuS-TAR 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 5month 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 hours, 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 July 2012 and early January 2013.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 NuSTAR-DAS 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 due to the overwhelming brightness of Mrk 421 in all epochs.

Fermi-LAT
The Large Area Telescope (LAT) on board the Fermi satellite is a pair-conversion telescope with energy coverage from 20 MeV to > 300 GeV.The LAT has a ∼ 2.4 sr field of view and provides all-sky monitoring coverage on a ∼ 3 hour time scale (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 > 100 • and any events collected when the spacecraft rocking angle was > 52 • .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 Tools111 .
The analyses considered data in day-long and weeklong 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 15 • regionof-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 Center112 .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 § 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.
2.6.MAGIC 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 • 53.4 ′ 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 hours of observation.Equivalently, a 1-hour integration yields a detection of a source with approximately 5% Crab flux.The angular resolution is 0.07 • (68% containment, > 200 GeV), and the energy resolution is 16%.The systematic uncertainties in the spectral measurements for a Crab-like pointsource 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 h between 2013 January 8 and 2013 March 18.Most of these observations (8 h 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 0.4 • 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 powerlaw 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.

VERITAS
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 hours of observations; equivalently, a source with approximately 5% Crab flux can be detected in a 1-hour integration.VERITAS has an energy resolution of ≃15% and an angular resolution (68% containment) of ∼ 0.1 • 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 website 113 .
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 h during the period 2013 Jan 10 to 2013 Mar 17, almost all of which is strictly simultaneous with NuSTAR exposures.These observations were taken at 0.5 • 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 parametrized 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 113 http://veritas.sao.arizona.edu/specificationsto 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 due 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.

Flux and Hardness Ratio Variability
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, and the range covered in each particular observation, are entirely dominated by the intrinsic variability of the target.For example, the calibration observation taken in July 2012 (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-hour period and dropped by almost a factor of 2 in 3 h.Several observations in 2013 have shown steadily decreasing count rate over the course of ≃12 h.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 March 2013 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.
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, R(t) ∝ e −t/τvar , 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 (τ var ) is found to be between 6 and 12 h.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 τ var = 9 ± 3 hours as our best estimate of this timescale.
In order to characterize the variability on shorter timescales (∆t ≪ τ var ), we consider the data in individual NuSTAR orbits as this represents a natural, albeit still arbitrary, way of partitioning the data.The 6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 7.0 7.1 7.2 7.3 7.4 7.5 7.6  -Count rates for NuSTAR in the 3-7 keV (blue) and 7-30 keV (green) bands.The legend given in the second panel from the top applies to all panels; for both bands FPMA count rates are plotted with a lighter color.The data are binned into 10-minute bins.The time axis of each row starts with the first day of the month and the UTC and MJD dates are printed above the light curves for each particular observation.Note that the data shown in the top panel represent two contiguous observations (broken up near MJD 56115.15).The intervals shaded in grey show times for which simultaneous data from MAGIC and VERITAS are presented in this paper.Both the horizontal and the vertical scales are equal in all panels.in each fit varies slightly due to the different duration of the orbits, but it is typically around 10.The right panel shows the distribution of the residual scatter after subtraction of the best-fit linear trend from the observed count rates in each orbit, in units of the median rate uncertainty within the orbit, (σ R ) orb .The colors reflect the mean count rate of the orbit: the lowest-rate third in red, the midrate third in orange and the highest-rate third in yellow, distributed as shown in the inset.The residual scatter distribution is slightly skewed to values greater than unity, indicating that 20% of orbits show excess variability on suborbital timescales.A Gaussian of approximately matching width is overplotted with a thick black line simply to highlight the asymmetry of the observed distribution.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, R(t) = constant, and a linear trend in time, i.e., R(t) ∝ t.The top panel of Figure 2 provides an example of both models fitted to the July 2012 data binned into 10-minute time bins, so that each orbit is divided into 5-7 bins per focal-plane module. 114The lower panel of Figure 2 shows the results of this fitting procedure performed on all 88 orbits.
We find that the flux during the majority of orbits is better described by a linear trend than by a constantflux model.Linear trends account for most of the orbitto-orbit variability and approximate smooth variations on super-orbital timescales of τ var ≈ 9 hours.In 10minute 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-χ 2 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, (σ R ) orb , gives a distribution slightly skewed towards 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.

Observation-averaged Spectroscopy
114 Because these data were taken in the calibration phase of the mission in suboptimal pointing conditions, a systematic uncertainty of 4% was added to the light curve to reflect the total uncertainty in the true count rate.a Energy (in keV) at which the model sharply changes slope from Γ1 to Γ2.For some observations this parameter is unbound, so we fix it to 7.5 keV and mark with (f).
b Flux calculated from the best-fit model, in units of 10 −11 erg s −1 cm −2 .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.c Formal 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.
We first model the NuSTAR spectra for each observation separately, before examining the clear intraobservation spectral variability in § 3.3 (see Figures 2 and  3, and Baloković et al. 2013b).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 N H = 1.92 × 10 20 cm −2 (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.
The fitting results imply that a power-law model with Γ ≈ 3 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 powerlaw model with two other simple models that allow extra degrees of freedom.The first one is a broken power-law model, bknpow: (2) 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 dependence of Γ 2 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 logparabolic 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, β = −0.2±0.1.
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.
3.3.Time-resolved Spectroscopy 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 § 3.1 that the dominant flux variations occur on a super-orbital timescale of τ var ≈ 9 hours.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 2,000 and 20,000 per focal-plane module.
As with the observation-averaged spectra, we fit powerlaw, 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 due to both the lack of intrinsic curvature and relatively poor photon statistics.In a similar manner to the observationaveraged 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 grey 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 Γ 1 and α, compared to the lower-amplitude variations in Γ 2 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).

Flux-resolved Spectroscopy
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 pro-vides 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 focalplane 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 Γ = 3.05 ± 0.02 from 3 to 70 keV (χ 2 r = 1.05).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 (β = 0.01 ± 0.04), 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 Γ ≈ 3. Extrapolating below 3 keV for the sake of comparison with the literature, we derive a 2-10 keV flux of (3.5 ± 0.2) × 10 −11 erg s −1 cm −2 for the Γ = 3 power-law model normalized to the lowest observed orbital flux (orbit #6 of the 2013 January 20 observation).8.The filled grey circles in the background show spectral modeling results for spectra of 88 individual orbits.
A state of such low X-ray flux has not yet been described in the published literature, 115 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- 115 To the best of our knowledge, the lowest published 2-10 keV flux thus far was (4.1 ± 0.2) × 10 −11 erg s −1 cm −2 (1997 May; Massaro et al. 2004).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 parametrize their dependence on the X-ray flux as where X stands for any of the spectral parameters (Γ, Γ 1 , E b , Γ 2 , α, β), s is the slope of the relation, F 3−7 keV is the 3-7 keV flux, F 0 is a reference flux in erg s −1 (chosen to be the median flux of our dataset, log F 0 = −10.1)and X 0 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 Γ 1 and Γ 2 parameters.For Γ 2 , 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 bestfit X (F 3−7 keV ) 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 0 , with their 1 σ uncertainties, in Table 8.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 analysed by Giebels et al. (2007) overlap in the 2-10 keV flux only for the highest-flux single-orbit data presented here, and extends almost a decade above that.These authors showed that the break energy is essentially independent of flux and E b = 5.9 ± 1.1 keV (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 −9 erg s −1 cm −2 , above which they saturate at Γ 1 ≈ 2.2 and Γ ≈ 2.5.The data presented here smoothly connect to those trends (see Figure 6), extending them towards the faint end.Whereas the low-energy photon index (Γ 1 or Γ L ) continues to increase with decreasing flux, reaching Γ ≈ 3 at 10 −10 erg s −1 cm −2 , the highenergy one (Γ 2 or Γ H ) essentially levels off to the same Γ ≈ 3 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 due to lack of sensitivity.

Multiwavelength Variability
The majority of observations performed in January through March 2013 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-hour long observations performed by NuSTAR were accompanied with Swift pointings at the beginning, middle and the end, and the groundbased Cherenkov-telescope arrays MAGIC and VERI-TAS covered roughly half of the NuSTAR exposure each.Approximately 50 h of simultaneous observations with NuSTAR and either MAGIC or VERITAS resulted in total exposure of 23.5 h (the remainder being lost due 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.
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 Xray flux is apparent from Swift -XRT data in comparison with the long-baseline median (Stroh & Falcone 2013;Baloković et al. 2013b) 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 10 s −1 in the 0.3-10 keV band are up to a factor of ≃ 2 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 October 2012 (see § 4.5 for more details).
Remarkably well-correlated flux variability in the Xray and VHE bands on a timescale of about a week is -Light curves for Mrk 421 from MAGIC , VERITAS (both above 200 GeV, binned in ∼30-minute intervals), Fermi -LAT (0.2-100 GeV, binned weekly), NuSTAR (3-30 keV, binned by orbit), Swift -XRT (0.3-10 keV, complete observations), Swift -UVOT (UVW1, UVM2 and UVW2 bands, complete observations), ground-based optical observatories (R band, intranight cadence), OVRO and Metsähovi (15 and 37 GHz, both with 3-4-day cadence).The host-galaxy contribution in the R band has been subtracted out according to Nilsson et al. (2007).The dynamic range in all panels is 40.Vertical and horizontal error bars show statistical uncertainties and the bin width, respectively, although some of the error bars are too small to be visible in this plot.The vertical lines mark midpoints of the coordinated NuSTAR and VHE observations: dashed lines mark the epochs for which we discuss SED snapshots in § 4.4, while the rest are shown with dotted lines.The horizontal lines in some panels show the long-term median values (see text for details).already apparent from Figure 7, and will be discussed in more detail in § 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 § 4.2.2.The activity observed in the first three months of 2013 can be generally described as low.Note in particular that on January 10 and 20, Mrk 421 showed a remarkably low Xray and VHE flux in comparison to the historical X-ray and VHE fluxes reported in Stroh & Falcone (2013) and Measurement uncertainties are based on photon statistics and are often smaller than the data points plotted.As in Figure 7, the vertical lines mark midpoints of the coordinated NuSTAR and VHE observations: dashed lines mark the epochs for which we discuss SED snapshots in § 4.4, while the rest are shown with dotted lines.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.
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 April 2013 (not shown in Figure 7), rivaling the brightest flares ever observed for Mrk 421 (Baloković et al. 2013a;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, F var , 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 116 do not directly add to the uncertainty in F var , it is important to stress that different observing sampling and, more importantly, different instrument sensitivity, do influence F var and its uncertainty: a densely sampled light curve with very small temporal bins and small error bars might 116 Estimated to be 20% in the VHE band, and around 10% in the optical, X-ray and GeV bands -see § 2 for details.
allow us to see flux variations that are hidden otherwise, and hence we might obtain a larger F var .Some practical issues of its application in the context of multiwavelength campaigns are elaborated in Aleksić et al. (2014); Aleksić et al. (2015b,c).
In this paper we explore two cases, as shown in Figure 9. First, we use the full January-March dataset 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. F var 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 F var determined from our campaign rises significantly from the radio towards 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.

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 which probe correlated variability on timescales as short as 1 hour (Giebels et al. 2007;Fossati et al. 2008;Acciari et al. 2009).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 −1 and VHE flux of 0.5 Crab units.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 Fig. 9.-Fractional variability amplitude, Fvar, as a function of frequency for the period 2013 January-March.The vertical error bars depict the statistical uncertainty, while the horizontal error bars show the energy band covered (with the markers placed in the center of the segments).We show Fvar computed in two ways: using the complete light curves acquired in the campaign (empty symbols), and using only the data taken within narrow windows centered on the coordinated NuSTAR and VHE observations (filled symbols).Note that the points overlap where only the coordinated observations are available.The Fermi -LAT point is based on the weekly-binned light curve shown in Figure 7.
in different X-ray bands using the strictly simultaneous Swift, NuSTAR, MAGIC and VERITAS observations. 117 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 § 3).
In the following, we use the discrete correlation function (DCF) and the associated uncertainty as defined in Edelson & Krolik (1989).We carried out the correlation analysis on two timescales: a ∼1-hour timescale, using strictly simultaneous observations, and a 1-day timescale, using data averaged over complete 6-10-hour 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 0.9 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 > 5 σ 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 log F X−ray − log F VHE correlation on both timescales may indicate that the correlation 117 The MAGIC and VERITAS observations reported in this paper were performed after the extensive hardware upgrades performed on these two facilities in 2011 and 2012.They are therefore much more sensitive than the ones performed in 2009, which allows for a significant detection of lower flux in a single night of observation.
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, due 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 τ var ≈9 hours (see § 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 nightlyaveraged 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) due 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 both with the spectral hardening already revealed by our analysis (see § 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., > 1 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 Xray flux was lowest, January 10 and 20 (MJD 56302 and 56312, respectively), steep spectra with Γ ≈ 3 were observed by NuSTAR and both VHE observatories.Other simultaneous observations yield Γ > 2.6 in the 3-30 keV band and Γ > 2.4 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 ≃ 2.5 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, higherquality data for the quiescent states are clearly needed in order to quantify it further.

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 (λ eff = 2120 Å) 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.
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 hours of each other the DCF is 0.2 ± 0.2, while for weeklyaveraged 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 due 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 In the top and bottom panels we overplot the slope of one-to-one proportionality (not a fit) with a black dotted line.Note that the vertical scale is linear in all panels, but has a different dynamic range in each one.
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, 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 For each of the epochs we show a log-parabolic curve fit to X-ray data alone (yellow) and all data (purple).For each curve, we mark the SED peak with an empty circle of matching color.Lower panel: Results of the SED peak localization based on data from strictly simultaneous Swift and NuSTAR orbits.The colored data points show ν syn.peak and (νFν ) syn. peak , i.e., the frequency of the SED peak and the flux at the peak.The assumption of the log-parabolic model connecting the UV/optical and X-ray data (purple empty circles) reveals a proportionality between log(νFν ) syn. peak and log ν syn.peak ; the dashed black line shows a linear fit best describing that relation.The other method (using only X-ray data) does not show a similar relation.In comparison with the observations published previously, shown here with different hatched grey regions, in 2013 January-March we observed a state in which the peak occurred at atypically low energy and high flux.
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 April 2013.
UV/optical and X-ray data, using a reasonable smooth interpolation or extrapolation model (e.g., Massaro et al. 2004;B laż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 Rband data taken within 24 hours.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; ii) fitting a log-parabolic model to both optical/UV and X-ray data.
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 16 to 10 17 Hz and a factor of ≃ 3 in peak flux.Interpretation of this simplistic parametrization 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 which additionally include the Rband 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-ofsight column density and extinction corrections.We conservatively estimate that the combination of these effects results in a factor of ≃2 uncertainty in the synchrotron peak frequency (ν syn.peak ), which dominates any statistical uncertainty from the fits.For this reason we do not show the uncertainties for individual ν syn.peak estimates in Figure 12.
Both methods consistently show the peak at an atypically low frequency (ν syn.peak < 10 17 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, log(νF ν ) syn. peak ∝ b log (ν syn.peak ), is b = 0.49 ± 0.08 ≈ 0.5.The dynamic range over which the relation holds, assuming that this parametrization 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 ν syn.peak < 10 17 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 ν syn.peak 10 16 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 16 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 17 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-ν syn.peak BL Lacs) if ν syn.peak is below 10 14 Hz, in the range 10 14 − 10 15 Hz, or above 10 15 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.

Broadband SED at Different Epochs
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 § 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 bestcharacterized 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 last several years with the launch of NuSTAR and upgrades to the MAGIC and VERITAS telescopes.
We note that the empirical SEDs from the 2013 campaign shown here represent 12 hours of observation in the X-ray and VHE bands, rather than integrations over -SED snapshots for four selected epochs during the campaign, assembled using simultaneous data from Swift -UVOT, Swift -XRT, NuSTAR, Fermi -LAT, MAGIC, and VERITAS.Most of the data were acquired over a period shorter than 12 hours in each case; the exceptions are the Fermi -LAT data and part of the radio data, which were accumulated over roughly one-week time intervals.The two left panels show low-state SEDs, while the two on the right show elevated states (not flaring, but among the highest presented in this paper).The grey symbols in the background of each panel show the SED of Mrk 421 from Abdo et al. (2011) averaged over a quiescent 4.5-month period.The solid blue lines show a simple one-zone SSC model discussed in § 5.3.To aid comparison, the model curve from the first panel is reproduced in the other panels with a blue dotted line.The dashed red lines show SED models with a time-averaged electron distribution discussed in § 5.3 for comparison with previously published results.
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 § 5.3 and § 5.4.

Brief Summary of the Flaring Activity in 2012
In addition to the coordinated multiwavelength campaign conducted in 2013, Mrk 421 was observed independently with several instruments in July-September 2012, including NuSTAR, Fermi and OVRO.In July 2012, the flux in the Fermi -LAT band increased above the median level and peaked twice over the following two 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 (1.4 ± 0.2) × 10 −6 s −1 cm −2 , a factor of ≃8 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 July-September 2012 are shown in Figure 14 in order to provide a timeline for this flaring event.
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 observation118 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 data119 that the X-ray flux in the 4-20 keV band increased further after the NuSTAR observation, peaking between 2 and  1 for greater time resolution), MAXI (4-20 keV, weekly bins) and OVRO (15 GHz, ∼weekly-daily cadence).Vertical and horizontal error bars show statistical uncertainties and the bin width, respectively, although some of the error bars are too small to be visible in this plot.The colored horizontal lines show the long-term median flux calculated from publicly available monitoring data.The dynamic range in all panels is 40, as in Figure 7, so that the two figures are directly comparable.
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(Hovatta et al. , 2015)).Approximately 70 days after the first peak of the γ-ray flare Mrk 421 reached a flux density of (1.11 ± 0.03) 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, both for the 2012 flare and the flare observed during the multiwavelength campaign in April 2013.As the latter event has been covered with numerous instru-ments (e.g., Baloković et al. 2013a;Cortina et al. 2013;Paneque et al. 2013;Pian et al. 2014), we defer a comprehensive analysis of this period to a separate publication.

Spectral Variability
in the X-ray Band The good temporal coverage of the NuSTAR data reveals a typical variability timescale of τ var ≈ 9 ± 3 hours.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 10hour 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 40% amplitude.We find that the X-ray spectrum of Mrk 421 cannot generally be described as a simple power law, but that 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 Γ ≈ 3, 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.

Correlated Variability in the X-ray
and VHE Spectral Bands 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 ≃ 80 keV and ≃ 1 TeV.All of this indicates that the energies of radiating particles must be very high (up to γ ∼ 10 6 ; see § 5.3 below) even when the source is in such a low state.
In § 4.2.1 we have shown that the X-ray and VHE fluxes are correlated at >3 σ significance.Parametrizing the correlation as log (F X−ray ) = a log (F VHE ) + b, the correlation is found to be approximately linear (a ≈1) both on half-hour and half-day timescales.This is consistent with most previous results considering similar spectral bands: a = 1.7 ± 0.3 (Tanihata et al. 2004), a ≈ 1 (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 due 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 a ≈ 2 and a = 2.9 ± 0.6, 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 § 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.Note.

Interpretation within the Framework
-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 § 5.3.Model SED curves are shown in Figure 13.

of a Single-zone SSC Model
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 (γ c , 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 ν syn.peak ≃ 10 16 Hz and ν IC peak ≃ 10 25 Hz.It then follows that γ c ∼ (ν IC peak /ν syn.peak ) 1/2 ∼ 3 × 10 4 .Assuming a bulk Doppler factor of δ ∼ 25, we can also estimate the magnetic field strength in the plasma comoving frame: B = 4πν syn.peak m e c/(3eδγ 2 c ) ≃ 0.1 G. Finally, an upper limit for the size of the emitting region in the comoving frame is given by the observed variability time scale, τ var ≈ 9 h, and the bulk Doppler factor δ: R cτ var δ ≃ 2 × 10 16 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).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 which propagate down the jet with a bulk Lorentz factor Γ. In order to decrease the number of free parameters, we assume a value Γ = 25 with the jet axis aligned near the line of sight with the critical angle θ = 1/Γ = 0.04 rad = 2.29 • , 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 "bulk Lorentz factor crisis", and requires that the radio and VHE emission 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 dn/dγ ∝ γ −q between γ min and γ max .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 τ esc = ηR/c is the escape time.This results in a particle distribution which propagates along the jet with power L e .Synchrotron emission results from the interaction of particles with a magnetic field B, generating a Poyntingflux 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.
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 Γ = 0.5 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, γ b = γ min , and has power-law shapes dn/dγ ∝ γ −2 for γ < γ b and dn/dγ ∝ γ −(q+1) for γ > γ b .For the parameters shown in Table 9, this particle distribution implies a synchrotron spectrum with a peak at ν syn.peak ≃ 1-2 × 10 16 Hz, and spectral shapes F ν ∝ ν −1/2 (dN/dE ∝ E −1.5 ) for ν < ν syn.peak and F ν ∝ ν −q/2 (dN/dE ∝ E −(2.6−3.2) ) for ν > ν syn.peak (cf.Table 7 and Figure 13).For γ min = 3 × 10 4 and B = 0.2 G, a synchrotron cooling timescale of τ syn = 4 × 10 5 s 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 τ esc, nom = R/c = 3 × 10 5 s.The escape efficiency factor, η = 35, ensures that the cooled electron distribu- Note.
-The electron energy distribution parameters listed here refer to the distribution directly responsible for the SSC emission.This simplified model is described in § 5.3 and used for comparison with the literature.Model SED curves are shown in Figure 13.tion 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 which are 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-hour 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 which 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 γ brk ≃ γ min and index p l = 2.0 below γ brk and index p h = q + 1 above the break.
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 δ = 21 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 1 = 2.2, and our values of p h = 4.0-4.6 are similar to their high-energy index of p 3 = 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 which is dictated by their generally broader SED peaks (see the grey 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 ( γ brk = 2.5-5.2 × 10 4 versus γ brk2 = 3.9 × 10 5 ), while the inferred emitting region radius is about a factor of 3-10 smaller (R = 0.6-1.7 × 10 16 cm versus 5.2 × 10 16 cm), the inferred magnetic field is substantially higher (B = 0.10-0.28versus 0.038), and the resulting jet powers differ by a factor of 3-4.As a consequence, the current modeling yields equipartition parameters that are much closer to unity, in the range ǫ = 0.33-0.56versus ǫ = 0.1 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 datasets (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 orderof-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 ν syn ∝ Bγ 2 , 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-ofmagnitude.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 ∼ 3 × 10 4 to ∼ 6 × 10 5 .Assuming that the target synchrotron photons for inverse-Compton scattering have energies around the synchrotron peak at ν syn.peak ∼ 10 16 Hz, the parameter governing the transition between Thomson and Klein-Nishina regimes is 4hν syn.peak γ/m e c 2 (Blumenthal & Gould 1970), which in the observer frame becomes 4hν syn.peak γ/δm e c 2 .
When considering photons from the synchrotron peak position (E = hν syn.peak ∼40 eV; i.e., about one order of magnitude lower than the typical position of the synchrotron peak in Mrk 421), we obtain 4hν syn.peak γ/δm e c 2 ≃ 0.4−8, indicating that the inverse-Compton scattering of photons with energy hν syn.peak 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 hν syn.peak , far above the range where Thomson scattering is relevant, and consistent with the linear (a ≃ 1) relationship between the soft X-ray and VHE flux.

Toward a multi-zone emission scenario
The electrons responsible for the broadband emission of Mrk 421 lose energy mostly due 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 τ syn = 1.2 × 10 3 B −3/2 E −1/2 δ −1/2 s, where E is photon energy in keV, and B is comoving frame magnetic field strength in G. Taking E ≈ 10 keV as the energy typical for the NuSTAR band, assuming B ≈ 0.2 G, as found from our SED modeling, and δ = 25 as before, we arrive at τ syn of ∼ 10 3 seconds.120This is more than an order of magnitude shorter than the variability timescale τ var ≈ 9 hours 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 particleacceleration 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 § 3.1, the observed increases appear at least as smooth and as slow as the observed decreases, in consistency with this picture.
One may argue that the X-ray flux variability reported in § 3.1 is not due to 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 11 and 10, 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 onezone SSC scenario in describing the broadband SED (see § 5.3), we argue that the observed multiwavelength variability and correlations point towards 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 121 , with strongest variability in the highestenergy electrons.However, the saturation of the X-ray spectral shape at the lowest and highest X-ray fluxes (see § 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 dom-121 In the sense of larger relative increase at higher energies.
inated 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 R G corresponding to that photon energy using B = 0.2 G and the maximal γ ∼ 10 6 .Since R G = γm e c 2 e −1 B −1 , we have R G 10 11 cm, which is much smaller than the inferred emission-region size of 10 16 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 non-uniform 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.

SUMMARY AND CONCLUSIONS
We have observed the blazar Mrk 421 in an intensive multiwavelength campaign in 2013, including GASP-WEBT, Swift, Fermi -LAT, MAGIC, VERITAS, and, for the first time, the new high-sensitivity hard X-ray instrument NuSTAR.In this paper we present part of the data from the campaign between the beginning of January and the end of March 2013, with the focus on the unprecedented coverage of the X-ray part of the broadband spectrum.Another successful aspect of the campaign is the achieved goal of strictly simultaneous observations in the X-ray and VHE γ-ray bands, in order to constrain the correlated variability.During the data-taking period presented in this work, Mrk 421 exhibited relatively low activity, including the lowest-flux state ever investigated with high temporal and broadband spectral coverage.
The rich data set yields a number of important empirical results: • During the first three months of 2013, the X-ray and VHE γ-ray activity of Mrk 421 was among the lowest ever observed.
• NuSTAR performed half-day long observations of Mrk 421 which showed that this source varies predominantly on timescales of several hours, with multiple instances of exponentially varying flux on timescales of 6-12 hours.Mrk 421 also exhibited smaller-amplitude, intrahour variations at the 5 % level.However, only 20 % of the Xray data show any appreciable intrahour variability.Within the dynamic range of our observations, we find no differences in the variability pattern or timescales between the lower and higher flux states.
• We find a systematic model-independent hardening of the X-ray spectrum with increasing X-ray flux.As the X-ray activity decreases, the curvature in the X-ray spectrum decreases and the spectral shape becomes softer.At 2-10 keV fluxes 10 −10 erg s −1 cm −2 , the spectral curvature completely disappears, and the spectral shape saturates into a steep Γ ≈ 3 power law, with no evidence for an exponential cutoff or additional hard components up to ≃ 80 keV.
• For two epochs of extremely low X-ray and VHE flux, in a regime not previously reported in the literature, we observed atypically steep spectral slopes with Γ ≈ 3 in both X-ray and VHE bands.Using a simple steady-state one-zone SSC scenario, we find that in these two epochs the peaks of both the synchrotron and inverse-Compton components of the SED shifted towards lower frequencies by more than an order of magnitude compared to their positions in the typical low states of Mrk 421 observed previously.The peak of the synchrotron bump of Mrk 421 shifted from ∼ 0.5 − 1 keV to ∼ 0.04 keV, which implies that HBLs can move towards becoming IBLs, leading to the conclusion that the SED classification based on the peak of the synchrotron bump may denote only a temporary rather than permanent characteristic of blazars.
• A clear double-bump structure is found in the fractional variability distribution, computed from radio to VHE γ-ray energies.This double-bump structure relates to the two peaks in the broadband SED shape of Mrk 421, and has been recently reported (with less resolution) for both lowactivity (Aleksić et al. 2015b) and high-activity states (Aleksić et al. 2015c).The less variable energy bands (radio, optical/UV and GeV γ-rays) relate to the segments of the SED rising up towards the peaks as a function of photon energy, while the most variable energy bands (X-rays and VHE γrays) sample the SED above the peaks, where it steeply declines with photon energy.
• We find a tight X-ray-VHE flux correlation in three nonoverlapping X-ray bands between 0.3 and 30 keV, with significantly different scaling.These results are consistent with an SSC scenario in which the X-ray and VHE radiation are produced by the same relativistic electrons, and the scattering of Xray photons to VHE energies (∼TeV) occurs in the less-efficient Klein-Nishina regime.From broadband SED modeling with a single-zone SSC model for four epochs, and assuming a constant Doppler factor of 25, we infer a magnetic field B ∼ 0.2 G and electron Lorentz factors as large as γ 6×10 5 .These parameter values, which are typical for describing the broadband SED of HBLs, further support the claim that, in the context of the SSC model, the inverse-Compton scattering responsible for the VHE emission takes place in the Klein-Nishina regime.
• There is tentative evidence for an optical/UV-GeV flux correlation, which is consistent with the emission in these two bands being produced by the same lower-energy electrons within the SSC framework.
• No correlation is found between fluxes in the optical/UV and the soft X-ray bands on either short or long timescales.However, we do find that a simple parametrization of the SED around the synchrotron peak with a log-parabolic function leads to a correlation between the peak flux and the frequency at which it occurs over a limited frequency range.
• The reported increase in the fractional variability with energy (for each of the two SED bumps), and the hardening of the X-ray spectra with increasing flux, suggest that the variability in the emission of Mrk 421 is produced by chromatic changes in the electron energy distribution, with the highestenergy electrons varying the most.The saturation of the X-ray spectral shape at the extremely high and low X-ray fluxes indicates that, for these periods of outstanding activity, the flux variability is instead dominated by other processes that lead to achromatic variations in the X-ray emission.
• The lifetimes of relativistic electrons due to synchrotron losses are estimated to be τ syn 10 3 seconds, which are substantially shorter than the ∼ 3 × 10 4 seconds that dominate the large-amplitude variations in the NuSTAR light curves.Together with the fractional variability distribution and the multiwavelength correlations observed in this campaign, this observation suggests that the broadband emission of Mrk 421 during low activity is produced by multiple emission regions.
• The electron cooling times of τ syn 10 3 seconds are also shorter than the emission-region crossing time ( 10 4 seconds), which points towards in situ electron acceleration.While particle acceleration in shocks is not excluded by our data, the gyroradii of the most energetic electrons (those radiating in the upper part of the NuSTAR band, or the upper part of the VHE band) is 10 11 cm, which is shorter than the cooling (energy-loss) timescales inferred from our modeling.This is suggestive of an electron-acceleration process occurring in relatively compact zones within a larger emission volume.and NNX15AU81G.
The OVRO 40-m monitoring program is supported in part by NASA grants NNX08AW31G and NNX11A043G, and NSF grants AST-0808050 and AST-1109911.
The Metsähovi team acknowledges the support from the Academy of Finland to our observing projects (numbers 212656, 210338, 121148, and others).This research has made use of NASA's Astrophysics Data System, and of Astropy, a communitydeveloped core Python package for astronomy (Astropy Collaboration 2013).
Facilities: NuSTAR, MAGIC, VERITAS, Fermi, Swift Fig.1.-Countrates for NuSTAR in the 3-7 keV (blue) and 7-30 keV (green) bands.The legend given in the second panel from the top applies to all panels; for both bands FPMA count rates are plotted with a lighter color.The data are binned into 10-minute bins.The time axis of each row starts with the first day of the month and the UTC and MJD dates are printed above the light curves for each particular observation.Note that the data shown in the top panel represent two contiguous observations (broken up near MJD 56115.15).The intervals shaded in grey show times for which simultaneous data from MAGIC and VERITAS are presented in this paper.Both the horizontal and the vertical scales are equal in all panels.

R
Fig. 2.-Upper Panel: Light curve of the July 2012 observation (FPMA in black, FPMB in grey) shown as an example for the count rate modeling; the two models fitted to each orbit of data are represented by purple (constant model) and blue lines (linear model).Lower Panels: The colored histograms on the left-hand side (colors matching the upper panel) show the distributions of reduced χ 2 for the two models fitted to every NuSTAR orbit up to the end of March 2013.The number of degrees of freedom (d.o.f.)in each fit varies slightly due to the different duration of the orbits, but it is typically around 10.The right panel shows the distribution of the residual scatter after subtraction of the best-fit linear trend from the observed count rates in each orbit, in units of the median rate uncertainty within the orbit, (σ R ) orb .The colors reflect the mean count rate of the orbit: the lowest-rate third in red, the midrate third in orange and the highest-rate third in yellow, distributed as shown in the inset.The residual scatter distribution is slightly skewed to values greater than unity, indicating that 20% of orbits show excess variability on suborbital timescales.A Gaussian of approximately matching width is overplotted with a thick black line simply to highlight the asymmetry of the observed distribution.

Fig. 3 .
Fig. 3.-Hardness ratio (defined as the ratio of the number of counts in the 7-30 keV band to that in the 3-7 keV band) as a function of the count rate in individual 30-minute bins of NuSTAR data is shown with colored symbols: FPMA plotted with squares, and FPMB with diamonds.The colors distingush different observations and match those in Figure 4.The thick black error bars and symbols show median count rate and hardness ratio in bins (1 ct s −1 width), delimited by the vertical dotted lines.The vertical error bars denote standard deviation within each bin.
Fig.4.-UnfoldedNuSTAR spectra of Mrk 421 in each of the 12 observations.Colors are arranged by the 3-7 keV flux.Also shown are spectra of the single highest-flux orbit (grey symbols) and the stack of three lowest-flux observations (black symbols).Modules FPMA and FPMB have been combined for clearer display and the bin mid-points for each spectrum are shown connected with lines of the same color to guide the eye.The left panel shows the unfolded spectra in the νFν representation, while the right panel shows the same spectra plotted as a ratio to the best-fit model for the lowest-flux stacked spectrum (a power law with Γ = 3.05).Note that the vertical scale is logarithmic in the left, and linear in the right panel.
Fig. 5.-Trends in the hard X-ray spectral parameters as functions of flux for three simple spectral models revealed by time-and flux-resolved analyses of the NuSTAR data.The values fitted to high-S/N stacked spectra separated by flux (see text for an explanation) are shown with black-lined empty circles.A linear function is fitted to each of the trends and the uncertainty region is shaded in grey.Parameters of the fitted linear trends are given in Table8.The filled grey circles in the background show spectral modeling results for spectra of 88 individual orbits.
Fig.6.-Comparison of the spectral trends revealed in our data (black symbols) with the ones published byGiebels et al. (2007) (grey symbols).For this set of fits to the NuSTAR data the break energy (E b ) was kept fixed at 7 keV.The uncertainties are given at the 68% confidence level in order to match the previous results.Note the smoothness of the trends covering nearly two orders of magnitude and the apparent saturation effects at each end.The dotted lines are median photon indices for 2-10 keV flux below 10 −10 erg s −1 cm −2 , representing the apparent low-flux saturation values.
Fig.7.-Light curves for Mrk 421 from MAGIC , VERITAS (both above 200 GeV, binned in ∼30-minute intervals), Fermi -LAT (0.2-100 GeV, binned weekly), NuSTAR (3-30 keV, binned by orbit), Swift -XRT (0.3-10 keV, complete observations), Swift -UVOT (UVW1, UVM2 and UVW2 bands, complete observations), ground-based optical observatories (R band, intranight cadence), OVRO and Metsähovi (15 and 37 GHz, both with 3-4-day cadence).The host-galaxy contribution in the R band has been subtracted out according toNilsson et al. (2007).The dynamic range in all panels is 40.Vertical and horizontal error bars show statistical uncertainties and the bin width, respectively, although some of the error bars are too small to be visible in this plot.The vertical lines mark midpoints of the coordinated NuSTAR and VHE observations: dashed lines mark the epochs for which we discuss SED snapshots in § 4.4, while the rest are shown with dotted lines.The horizontal lines in some panels show the long-term median values (see text for details).
Fig. 8.-Optical polarization of Mrk 421 between 2013 January and March.The degree of polarization is shown in the upper panel and the position angle of polarization is shown in the lower panel.Measurement uncertainties are based on photon statistics and are often smaller than the data points plotted.As in Figure7, the vertical lines mark midpoints of the coordinated NuSTAR and VHE observations: dashed lines mark the epochs for which we discuss SED snapshots in § 4.4, while the rest are shown with dotted lines.
Fig.10.-Flux-fluxcorrelation between the X-ray and VHE (>200 GeV) flux in three different X-ray bands: Swift -XRT 0.3-3 keV in the left panel, Swift -XRT and NuSTAR 3-7 keV in the middle, and NuSTAR 7-30 keV in the right panel.Swift -XRT and NuSTAR measurements are shown with black-filled and white-filled symbols, respectively.Orange symbols mark MAGIC measurements, while dark red mark VERITAS.In the upper panels we show only the data taken essentially simultaneously (within 1.5 hours).The lower panels show data averaged over the nights of simultaneous observations with X-ray and VHE instruments.The N , a and DCF values given in each panel are the number of data points considered, the slope of the log-log relation, and the discrete correlation function.The best fit linear relation (dashed grey line) and its uncertainty region are shown with grey shading.The thin dotted line of slope unity is shown in all panels for comparison.
Fig. 11.-Flux-flux correlations between the UV band (Swift -UVOT filter UVW1) and the optical (R band; top), soft X-ray (Swift -XRT 0.3-3 keV; middle) and GeV γ-ray (Fermi -LAT 0.2-100 GeV; bottom) bands.The green data points are based on flux measurements that are coincident within 6 hours.The blue data points are derived by averaging over 7-day intervals, i.e., integration times used for determination of the flux in the Fermi -LAT band.In the top and bottom panels we overplot the slope of one-to-one proportionality (not a fit) with a black dotted line.Note that the vertical scale is linear in all panels, but has a different dynamic range in each one.
Fig. 12.-Upper panel: Examples of approximate localization of the synchrotron SED component peak for two orbits of simultaneous observations with Swift and NuSTAR.The Swift data are shown as black filled symbols (diamonds for the UVOT and squares for the XRT) and the NuSTAR data are shown as black empty squares.Empty diamonds represent R-band data.For each of the epochs we show a log-parabolic curve fit to X-ray data alone (yellow) and all data (purple).For each curve, we mark the SED peak with an empty circle of matching color.Lower panel: Results of the SED peak localization based on data from strictly simultaneous Swift and NuSTAR orbits.The colored data points show ν syn.peak and (νFν ) syn. peak , i.e., the frequency of the SED peak and the flux at the peak.The assumption of the log-parabolic model connecting the UV/optical and X-ray data (purple empty circles) reveals a proportionality between log(νFν ) syn. peak and log ν syn.peak ; the dashed black line shows a linear fit best describing that relation.The other method (using only X-ray data) does not show a similar relation.In comparison with the observations published previously, shown here with different hatched grey regions, in 2013 January-March we observed a state in which the peak occurred at atypically low energy and high flux.
Fig.13.-SED snapshots for four selected epochs during the campaign, assembled using simultaneous data from Swift -UVOT, Swift -XRT, NuSTAR, Fermi -LAT, MAGIC, and VERITAS.Most of the data were acquired over a period shorter than 12 hours in each case; the exceptions are the Fermi -LAT data and part of the radio data, which were accumulated over roughly one-week time intervals.The two left panels show low-state SEDs, while the two on the right show elevated states (not flaring, but among the highest presented in this paper).The grey symbols in the background of each panel show the SED of Mrk 421 fromAbdo et al. (2011) averaged over a quiescent 4.5-month period.The solid blue lines show a simple one-zone SSC model discussed in § 5.3.To aid comparison, the model curve from the first panel is reproduced in the other panels with a blue dotted line.The dashed red lines show SED models with a time-averaged electron distribution discussed in § 5.3 for comparison with previously published results.
Fig. 14.-Light curves for Mrk 421 around the γ-ray flare detected by Fermi -LAT in 2012, from Fermi -LAT (0.2-100 GeV, binned weekly), NuSTAR (3-30 keV, binned by orbit; see the top panel of Figure1for greater time resolution), MAXI (4-20 keV, weekly bins) and OVRO (15 GHz, ∼weekly-daily cadence).Vertical and horizontal error bars show statistical uncertainties and the bin width, respectively, although some of the error bars are too small to be visible in this plot.The colored horizontal lines show the long-term median flux calculated from publicly available monitoring data.The dynamic range in all panels is 40, as in Figure7, so that the two figures are directly comparable.

TABLE 1
Summary of the Swift observations of Mrk 421 (January-March 2013)

TABLE 2
Models fitted to the Swift -XRT spectra of each observation Flux calculated from the best-fit model, in units of 10 −11 erg s −1 cm −2 . a

TABLE 3
Summary of the NuSTAR observations of Mrk 421 (January-March 2013) a Livetime-corrected sum of all good time intervals comprising the observation.b PSF-corrected source count rate and its uncertainty in the 3-30 keV band averaged over the exposure time.

TABLE 6 Summary of the VERITAS observations of Mrk 421 (January-March 2013)
o.f.Power-law model of the form dN/dE = F0 (E/E ′ ) −Γ , E ′ = 300 GeV, is fitted to each observation.A log-parabolic model fit of the form dN/dE = F0 (E/E ′ ) −α−β log(E/E ′ ) is shown for observations in which it provides a better description of the spectrum than the power-law model.The normalization constant, F0, is in units of 10 −10 s −1 cm −2 TeV −1 .
b c Quoted uncertainties are statistical.

TABLE 7
Models fitted to the NuSTAR spectra of each observation

TABLE 8
Best-fit linear relations parametrized in Equation (4) for describing the change of spectral parameters with the X-ray flux.The uncertainties are quoted at the 1 σ level.Due to poor constraint on this parameter, it is kept fixed at 7 keV while quantifying the trends in the Γ1 and Γ2 parameters. a