First Results on Long-Term Multiwavelength Variability Analysis in Mrk 421

Due to the broad-band nature of the emission from blazars, strong efforts are being placed in recent years on simultaneous multiwavelength campaigns as a way to understand the mechanisms responsible for the acceleration of ultra-relativistic particles in these objects. The present work aims at providing a systematic study of the broadband long-term variability properties of the TeV blazar, Mrk 421, by gathering archival multiwavelength data spanning nearly two decades in time over several energy bands, from radio frequencies through x-rays to high and very high γ-ray energies. The lightcurves at these energies have been characterized over long periods (≈months) using several statistical methods, in order to study variability over long-term timescales and the dependence of this variability with flux and energy, aiming ultimately at shedding some light onto the physical mechanisms that drive blazar emission. We present in this work the preliminary results of the variability analysis for Mrk 421, and it will be extended in the future to see if the properties found are particular of this source or if they can be extended to the blazar class.


Introduction
According to the unified scheme of Active Galactic Nuclei (AGN), blazars are radio-loud AGN that display highly variable, beamed, non-thermal emission, covering a broad range from radio to gamma-ray energies [1].The blazar class encompasses BL Lacertae (BL Lac) and flat spectrum radio quasars (FSRQs) objects, whose main differences appear in their emission lines and their Spectral Energy Distribution (SED) properties.Observationally, blazars are characterized by core-dominated emission and rapid variability, and this feature provides limits to the size and the speed of the emitting region.A blazar SED, in a νFν representation, shows two distinctive broad humps arising from different physical processes: synchrotron emission in low energies, and a high energy process of leptonic (e.g.[2] and [3]) or hadronic [4] nature still to be defined, like inverse Compton, π-0 decay or synchrotron emission.Simultaneous multiwavelength campaigns are being carried out in recent years to try differentiate between these two competing models.In the present work we will exploit two decades of archive multiwavelength data for the TeV blazar Mrk 421, to characterize its variability features over long-term timescales and the dependence of this variability with flux and energy.

Observational Data
The seed of our work is a sample of 10 TeV blazars observed by the MAGIC telescopes, two Cherenkov Telea e-mail: eracero@sciops.esa.intb e-mail: icalle@sciops.esa.intscopes with mirrors of 17m diameter, that became available at the end of 2010 into FITS format.Around that time, Tluczykont et al.(2010) combined the Very High Energy (VHE) data of Mrk 421, Mrk 501 and 1ES 1959+650 that had been published by different Imaging Cherenkov Telescopes (ICTs) for the last 20 years, and made them publicly available in FITS format [5].Our main goal was to collect archival data at other wavebands, and in combination with the Cherenkov Telescope data, study in a systematic way the long-term behaviour of these sources in terms of their variability.To achieve our goals, we applied known statistical tests for variability to Mrk 421, the first extragalactic source detected in the Very High Energy sky with the Cherenkov Imaging Telescope at the Whipple Observatory [6], and one of the most monitored TeV blazars at all wavebands with plenty of archival data.F>1TeV (Crab) 1992(Crab) 1993(Crab) 1994(Crab) 1995(Crab) 1996(Crab) 1997(Crab) 1998(Crab) 1999(Crab) 2000(Crab) 2001(Crab) 2002(Crab) 2003(Crab) 2004(Crab) 2005(Crab) 2006(Crab) 2007(Crab) 2008 1 presents a brief summary of the missions/observatories and the websites where the data can be found.All datasets presented in this table were already analyzed by the corresponding mission pipeline or collaboration working group, except for the Large Area Telescope (LAT) data from the Fermi observatory, an waveband we considered relevant for this work due to its strategical spectrum coverage.
All lightcurves used in this work are shown in figure 1, where the data was binned using 14 day bins (the time bin given by UMRAO for the 4.8, 8 and 14.5 GHz datasets) using the same reference time (MJD= 48257 days) for comparison reasons.We computed the weighted average of the data points included in each bin, and associated the uncertainty bars with the error of the weighted mean.To ensure good statistical quality of our data, we impose several filtering criteria on the lightcurves, and only positive flux values (F>0) and 3σ detections were considered for this work.

Data Analysis
This section describes the analysis carried on the individual lightcurves focusing on the long-term variability properties of Mrk 421, and how these properties vary with frequency band.

Lightcurve Characterization
The idea that TeV Blazars show two distinctive states in γ-rays, a baseline or quiescent state, and epochs of flar-ing states superimposed, has being recently discussed [7], where it is proposed a theoretical model that would fit the behaviour observed.Tluczykont et al.(2010) already interpreted from the long-term state distribution that the behaviour of Mrk 421 could be well described by a baseline flux plus a stochastic signal (flares) governed by a multiplicative process.
We derive a mean flux for each lightcurve over the whole time period by defining what we called baseline flux as the quiescent state of the source.In order to derive the baseline flux, we identify the source emission states (flaring and quiescent) and their evolution with time through an iterative process, where the 3σ data points above a linear fit to the lightcurve in each iteration are rejected.The lower limit was set in 15σ to guarantee all the values considered as the baseline state are included.This process ends up when all the points of the baseline flux are within the 3σ and 15σ region.We consider the parameters obtained from the last linear fit as the best estimate of the dependence of the baseline flux with time.From these fits, amongst other things, we derive annual flux variations for all frequency bands, shown in Table 2.The distribution of fluxes of a given lightcurve is an indicator of the nature of the underlying process responsible for the variability observed.To ensure a meaningful representation of the observed fluxes, only good quality data is used, i.e., only flux measurements at the 3σ level are used in the distributions.To account for the different uncertainties on the flux measurements of the different lighcurves, the binning of each one has been chosen such that one bin equals the average statistical error of individual flux measurements.This also ensures that the shape of the histograms is not dominated by the errors of the individual flux measurements.We normalized the flux distributions to their respective means in order to remove the time dependency of the mean flux.We fit each histogram with a Gaussian and Lognormal distributions as well as with a combination of both.Our aim with this exercise is not to extract any meaningful information from these fits, but simply to test both hypothesis.The preliminary results of this exercise can be found in figure 2.

Stationarity
Following Vaughan et al.(2003), we can distinguish weakly-non stationary from strong-non stationary processes studying the variance of AGN lightcurves.This is 04020-p.2based on the premise that, even if every realization of a red-noise process will differ due to the intrinsic random nature of the process, the Power Spectral Density (PSD) will tend to its 'true' value with time and number of data points, and so do the variance, the integral of the PSD.
The intrinsic or excess variance [8] of a time series is defined as the true variance of a sample, taking into account the measurement errors associated σ2 err : where S 2 is the variance and σ2 err is the average mean square errors.The normalized excess variance is thus computed as: Last, the square root of σ 2 NXS is the fractional root mean square (rms) variability amplitude, Fvar, a useful variable to measure the intrinsic variability amplitude in a multiwavelength framework for comparison reasons.Fig. 3 shows the full analysis carried out over the soft x-ray time series as an illustration of the line of work followed.For a red-noise process, the variance variability depends on the power spectrum of the variations and on the time resolution and duration of the lightcurve.This means that measurements of the variance can only be compared in a multiwavelength framework if they have the same duration [9].Thus, we use a common time reference (48257d MJD) and a common time bin to compute the variance estimators described above in all the datasets.
Figure 3 shows as an example the full excess variance analysis carried out on the soft x-ray time series.Red points correspond to the magnitudes (F x , S 2 , σ 2 xs and Fvar) calculated over a fixed time bin of 180 days.This bin was selected as a trade-off between good statistics per segment in all the energy bands to compute S 2 , σ 2 xs and Fvar with good accuracy while maximizing the number of intervals required to extract meaningful results.This compromise was guided by the binning time of the radio datasets (14day bins).To fulfill the minimum condition of N > 10 prescribed in Vaughan et al.(2003) [8] in order to assume a normal distribution of variances, the minimum required time bin is 140 days, and taking into account the gaps exhibited in our time series, we extended the interval up to 6 months.Dashed and dotted black lines correspond respectively to the boundaries of the 90 and 99% confidence regions estimated by Vaughan et al.(2003) [8] for a simulated stationary red-noise process and summarized in Table 1 of the referenced paper.To be conservative, we assume the largest uncertainties associated with a PSD α = −2, N = 10 number of points per segment, which correspond to error values of logS 2 = +0.46−0.78 for the 90% and logS 2 = +0.75−1.16 for the 99% confidence regions.xs , computed for each time bin (180 days).We considered upper limits those cases where σ 2 xs < 0 and σ 2 xs < 1 − σ.Panel 4-Fvar calculated over equally time intervals of 180 days.The variance estimators (panel 2-4) show dashed and dotted lines representing the 90 and 99% confidence levels expected for a red-noise variability [8], and the blue solid lines show the total unweighted mean in all cases.

RMS Spectrum
The 'rms spectrum' measures the dependence of the variability with frequency, where the variability amplitude is computed via Fvar or the fractional root mean square, and it may give us a measure of the spectral variability of a source [10].It has being used for X-ray timing analysis mainly, although it is being applied to other energy bands lately [11].To compute the fractional variability amplitude as a function of energy range in a systematic way for all the time series, the total unweighted average of Fvar 04020-p.3 The Innermost Regions of Relativistic Jets and Their Magnetic Fields with its standard deviation was computed following the steps described in previous Sect.3.2.The rms spectrum of Mrk 421, measured for a 6 month-period, can be seen in figure 4.

Conclusion
The study of the variability features of AGNs in a multiwavelength framework is crucial in the understanding of the relation among the physical properties of these sources observed at different frequencies.To this end we collected archival data spanning almost two decades in time for the TeV blazar Mrk 421 in order to characterize its longterm emission and variability properties in several energy bands, from radio to the highest gamma-ray energy range.The goal of this work is to establish a well-defined set of statistical tools to analyze the aperiodic variability of Mrk 421 systematically to later extend it to a larger population and see if the variability properties derived in this work are intrinsic of the blazar class or particular of the source of study.
The analysis developed here explores several properties of the source, namely the flux evolution with time and its characterization over twenty years of data, and the flux variability by means of variance estimators.The variance estimators help to determine the stationarity of the emission process.More importantly, these properties are determined in all the different wavebands available which let us look into the dependence of the observed variability with energy via the rms spectrum.
The mean baseline flux analysis shows an increase trend of the baseline flux with time for the γ-ray wave-lengths and for both the soft and hard x-ray regimes.However, there seems to be an opposite trend in the radio bands.We will further investigate this behaviour through crosscorrelation analysis.The results from the study of the flux emission show that both the x-ray and γ-ray energy flux distributions can be well described by a log-normal distribution, indicative of underlying multiplicative processes.In the radio regime (4.8, 8 and 15 GHz), the flux distribution can be described by Gaussian distributions, indicative of underlying independent and additive processes.
The variance estimators analysis shows that the emission is consistent with a stationary process for the all radio and Fermi γ-ray time series in timescales from 1 day to 6 months.The soft and hard x-ray time series, and the ICTs γ-ray (> 1T eV) data display a strong non-stationarity for the same timescales.The rms spectrum shows two distinctive peaks with mean Fvar around 50% at the mentioned energy regimes.
Further work will include the verification of these results from the novel analysis of the variance estimators through the study of the Power Spectral Density of the whole set of lightcurves.Due to limitations of the size of the samples, we will implement the logarithmic PSD computation prescribed by Papadakis & Lawrence (1993) to investigate its possible implementation in the tool chain.After the characterization of the variability properties in each time series, the natural step is studying the possible correlations among frequency bands.To do so we intend to make use of the Discrete Correlation Function (DCF) analysis defined by Edelson & Krolik (1998) in order to compute the coherence of the cross-correlation function for discrete time series.

Figure 2 .
Figure 2. Mrk 421 flux characterization at different wavebands.From left to the right: Normalized flux distributions, (G) Gaussian fit, (L) lognormal fit and (G+L) a combination of gaussian plus lognormal fit.Different shading correspond to different events, derived from the baseline characterization of each waveband.

Figure 3 .
Figure 3. From top to bottom: panel 1-initial lightcurve (grey points) with the weighted mean rate (x) of the initial time series binned over 180 days (red points).Blue solid line correspond to the weighted mean of the x values.Panel 2-S 2 calculated over 180-day-interval.Panel 3-σ 2xs , computed for each time bin (180 days).We considered upper limits those cases where σ 2 xs < 0 and σ 2 xs < 1 − σ.Panel 4-Fvar calculated over equally time intervals of 180 days.The variance estimators (panel 2-4) show dashed and dotted lines representing the 90 and 99% confidence levels expected for a red-noise variability[8], and the blue solid lines show the total unweighted mean in all cases.

Figure 4 .
Figure 4. Fvar spectra for ≈6 months of data, averaged over the whole time period for each lightcurve.Vertical bars denote 1-σ uncertainties and the horizontal bars indicate the energy width covered by each instrument.

Table 1 .
Observatory and mission public or private archives where Mrk 421 data was obtained.

Table 2 .
Summary of the mean flux variation per year derived from the linear fit to the baseline flux.