VHE $\gamma$-ray emission of PKS 2155$-$304:
spectral and temporal variability

The integral flux for a given period of observations is determined in a standard way. For subsequent discussion purposes, the formula applied is given here :

where $T$ represents the corresponding live-time, $A(E)$ and $R(E,E^{\prime})$ are, respectively, the collection area at the true energy $E$ and the energy resolution function between $E$ and the measured energy $E^{\prime}$, and $S(E)$ the shape of the differential energy spectrum as defined in Eq. 1 and 2 . Finally, $N_{\mathrm{exp}}$ is the number of measured events in the energy range [E_min,E_max].

In the case that $S(E)$ is a power law, an important source of systematic error in the determination of the integral flux variation with time comes from the value chosen for the index $\Gamma$. The average 2005–2007 energy spectrum yields a very well determined power law index¹¹1 The average 2005–2007 energy spectrum yields a power law with a photon index $\Gamma=3.37\pm 0.02_{\mathrm{stat}}$. One should note that some curvature is observed at higher energies, resulting in a better spectral determination when the alternative hypothesis shown in (Eq. 2) is used, yielding a harder index ($\Gamma=3.05\pm 0.05_{\mathrm{stat}}$) with an exponential cut-off at energy $E_{\mathrm{cut}}=1.76\pm 0.27_{\mathrm{stat}}\leavevmode\nobreak\ \mathrm{TeV}$. This curved model is prefered at a $8.4\,\sigma$ level as compared to the power law hypothesis. However, the choice of the model has little effect on the determination of the integral flux values above $200\,{\rm GeV}$, the integral being dominated by the low-energy part of the energy spectrum.. However, in Section 4 it will be shown that this index varies depending on the flux level of the source. Moreover, in some cases the energy spectrum of the source shows some curvature in the TeV region, giving slight variations in the fitted power law index depending on the energy range used.

For runs whose energy threshold is lower than $E_{\mathrm{min}}$, a simulation performed under the observation conditions corresponding to the data shows that an index variation of $\Delta\Gamma=0.1$ implies a flux error at the level of $\Delta\Phi\sim 1\%$, this relation being quite linear up to $\Delta\Gamma\sim 0.5$. However, this relation no longer holds when the energy threshold is above $E_{\mathrm{min}}$, as the determination of $\Phi$ becomes much more dependent on the choice of $\Gamma$. For this reason, only runs whose energy threshold is lower than $E_{\mathrm{min}}$ will be kept for the subsequent light curves. The value of $E_{\mathrm{min}}$ is chosen as $200\,{\rm GeV}$, which is a compromise between a low value which maximises the excess numbers used for the flux determinations and a high value which maximises the number of runs whose energy threshold is lower than $E_{\mathrm{min}}$.

From 2005 to 2007, PKS 2155$-$304 is almost always detected when observed (except for two nights for which the exposure was very low), indicating the existence, at least during these observations, of a minimal level of activity of the source. Focussing on data set $D_{QS}$ (which excludes the July 2006 data where the source is in a high state), the distribution of the integral fluxes of the individual runs above $200\,{\rm GeV}$ has been determined for the 115 runs, using a spectral index $\Gamma=3.53$ (the best value for this data set, as shown in 3.4). This distribution has an asymmetric shape, with mean value $(4.32\pm 0.09_{\mathrm{stat}})\times 10^{-11}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$ and root mean square (r.m.s.) $(2.48\pm 0.11_{\mathrm{stat}})\times 10^{-11}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$, and is very well described with a lognormal function. Such a behaviour implies that the logarithm of fluxes follows a normal distribution, centered on the logarithm of $(3.75\pm 0.11_{\mathrm{stat}})\times 10^{-11}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$. This is shown in the left panel of Fig. 3, where the solid line represents the best fit obtained with a maximum-likelihood method, yielding results independent of the choice of the intervals in the histogram. It is interesting to note that this result can be compared to the fluxes measured by H.E.S.S. from PKS 2155$-$304 during its construction phase, in 2002 and 2003 (see Aharonian et al.  2005b and Aharonian et al.  2005c). As shown in Table 3, these flux levels extrapolated down to $200\,{\rm GeV}$ were close to the value corresponding to the peak shown in the left panel of Fig. 3.

The right panel of Fig. 3 shows how the flux distribution is modified when the July 2006 data are taken into account (data set $D$ in Table 2): the histogram can be accounted for by the superposition of two Gaussian distributions (solid curve). The results, summarized in Table 4, are also independent of the choice of the intervals in the histogram. Remarkably enough, the characteristics of the first Gaussian obtained in the first step (left panel) remain quite stable in the double Gaussian fit.

This leads to two conclusions. First, the flux distribution of PKS 2155$-$304 is well described considering a low state and a high state, for each of which the distribution of the logarithms of the fluxes follows a Gaussian distribution. The characteristics of the lognormal flux distribution for the high state are given in Sections 5, 6 and 7. Secondly, PKS 2155$-$304 has a level of minimal activity which seems to be stable on a several-year time-scale. This state will henceforth be referred to as the “quiescent state” of the source.

In order to determine if the measured width of the flux distribution (left panel of Fig. 3) can be explained as statistical fluctuations from the measurement process a simulation has been carried out considering a source which emits an integral flux above $200\,{\rm GeV}$ of $4.32\times 10^{-11}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$ with a power law photon spectrum index $\Gamma=3.53$ (as determined in the next section). For each run of the data set $D_{QS}$ the number $n_{\gamma}$ expected by convolving the assumed differential energy spectrum with the instrument response corresponding to the observation conditions is determined. A random smearing around this value allows statistical fluctuations to be taken into account. The number of events in the off-source region and also the number of background events in the source region are derived from the measured values $n_{\mathrm{off}}$ in the data set. These are also smeared in order to take into account the expected statistical fluctuations.

10000 such flux distributions have been simulated, and for each one its mean value and r.m.s. (which will be called below RMSD) are determined. The distribution of RMSD thus obtained, shown in Fig. 4, is well described by a Gaussian centred on $0.98\times 10^{-11}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$ (which represents a relative flux dispersion of 23%) and with a $\sigma_{RMSD}$ of $0.07\times 10^{-11}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$.

It should be noted that here the effect of atmospheric fluctuations in the determination of the flux is only taken into account at the level of the off-source events, as these numbers are taken from the measured data. But the effect of the corresponding level of fluctuations on the source signal is very difficult to determine. If a conservative value of 20% is considered ²²2 A similar procedure has been carried out on the Crab Nebula observations. Considering this source to be perfectly stable it allows us to determine an upper limit to the fluctuations of the Crab signal due to the atmosphere, and a value of $\sim 10\%$ was derived. Nonetheless, this value is linked to the observations’ epoch and zenith angles, and to the source spectral shape. , which is added in the simulations as a supplementary fluctuation factor for the number of events expected from the source, a RMSD distribution centred on $1.30\times 10^{-11}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$ with a $\sigma_{RMSD}$ of $0.09\times 10^{-11}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$ is obtained. Even in this conservative case, the measured value for the flux distribution r.m.s. ($(2.48\pm 0.11_{\mathrm{stat}})\times 10^{-11}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$) is very far (more than 8 standard deviations) from the simulated value. All these elements strongly suggest the existence of an intrinsic variability associated with the quiescent state of PKS 2155$-$304.

The energy spectrum associated with the data set $D_{QS}$, shown in Fig. 5, is well described by a power law with a differential flux at 1 TeV of $\phi_{0}=(1.81\pm 0.13_{\mathrm{stat}})\times 10^{-12}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}\mathrm{TeV}^{-1}$ and an index of $\Gamma=3.53\pm 0.06_{\mathrm{stat}}$. The stability of these values for spectra measured separately for 2005, 2006 (excluding July) and 2007 is presented in Table 5. The corresponding average integral flux is $(4.23\pm 0.09_{\mathrm{stat}})\times 10^{-11}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$, which is as expected in very good agreement with the mean value of the distribution shown in the left panel of Fig. 3.

Bins above 2 TeV correspond to $\gamma-$ray excesses lower than 20$\,\gamma$ and significances lower than $2\,\sigma$. Above 5 TeV excesses are even less significant ($\sim 1\sigma$ or less) and 99% upper-limits are used. There is no improvement of the fit when a curvature is taken into account.

The spectral state of PKS 2155$-$304 has been monitored since 2002. The first set of observations (Aharonian et al.  2005b), from July 2002 to September 2003, shows an average energy spectrum well described by a power law with an index of $\Gamma=3.32\pm 0.06_{\mathrm{stat}}$, for an integral flux (extrapolated down to $200\,{\rm GeV}$ ) of $(4.39\pm 0.40_{\mathrm{stat}})\times 10^{-11}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$. No clear indication of spectral variability was seen. Consecutive observations in October and November 2003 (Aharonian et al.  2005c) gave a similar value for the index, $\Gamma=3.37\pm 0.07_{\mathrm{stat}}$, for a slightly higher flux of $(5.22\pm 0.54_{\mathrm{stat}})\times 10^{-11}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$. Later, during H.E.S.S. observations of the first (MJD 53944, Aharonian et al.  2007a) and second (MJD 53946, Aharonian et al.  2009) exceptional flares of July 2006, the source reached much higher average fluxes, corresponding to $(1.72\pm 0.05_{\mathrm{stat}})\times 10^{-9}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$ and $(1.24\pm 0.02_{\mathrm{stat}})\times 10^{-9}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$ ³³3corresponding to data set T200 in Aharonian et al.  (2009) respectively. In the first case, no strong indications for spectral variability were found and the average index $\Gamma=3.19\pm 0.02_{\mathrm{stat}}$ was close to those associated with the 2002 and 2003 observations. In the second case, clear evidence of spectral hardening with increasing flux was found.

The observations of PKS 2155$-$304 presented in this paper also include the subsequent flares of 2006 and the data of 2005 and 2007. Hence, the evolution of the spectral index is studied for the first time for a flux level varying over two orders of magnitude. This spectral study has been carried out over the fixed energy range 0.2–1 TeV in order to minimize both systematic effects due to the energy threshold variation and the effect of the curvature observed at high energy in the flaring states. The maximal energy has been chosen to be at the limit where the spectral curvature seen in high flux states begins to render the power law or exponential curvature hypotheses distinguishable. As flux levels observed in July 2006 are significantly higher than in the rest of the data set (see Fig. 6), the flux–index behaviour is determined separately first for the July 2006 data set itself ($D_{JULY06}$) and secondly for the 2005-2007 data excluding this data set ($D_{QS}$).

On both data sets, the following method was applied. The integral flux was determined for each run assuming a power law shape with an index of $\Gamma=3.37$ (average spectral index for the whole data set), and runs were sorted by increasing flux. The set of ordered runs was then divided into subsets containing at least an excess of 1500 $\gamma$ above $200\,{\rm GeV}$ and the energy spectrum of each subset was determined⁴⁴4even for lower fluxes, the significance associated with each subset is always higher than 20 standard deviations.

The left panel of Fig. 7 shows the photon index versus integral flux for data sets $D_{QS}$ (grey crosses) and $D_{JULY06}$ (black points). Corresponding numbers are summarized in Appendix B. While a clear hardening is observed for integral fluxes above a few $10^{-10}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$, a break in this behaviour is observed for lower fluxes. Indeed, for the data set $D_{JULY06}$ (black points) a linear fit yields a slope $\rm d\Gamma$/$\rm d\Phi=(3.0\pm 0.3_{\mathrm{stat}})\times 10^{8}\mathrm{cm}^{2}\mathrm{s}$, whereas the same fit for data set $D_{QS}$ (grey crosses) yields a slope $\rm d\Gamma$/$\rm d\Phi=(-3.4\pm 1.9_{\mathrm{stat}})\times 10^{9}\mathrm{cm}^{2}\mathrm{s}$. The latter corresponds to a $\chi^{2}$ probability $\rm P(\chi^{2})=71\%$; a fit to a constant yields $\rm P(\chi^{2})=33\%$ but with a constant fitted index incompatible with a linear extrapolation from higher flux states at a 3$\,\sigma$ level. This is compatible with conclusions obtained either with an independant analysis or when these spectra are processed following a different prescription. In this prescription the runs were sorted as a function of time in contiguous subsets with similar photon statistics, rather than as a function of increasing flux.

The form of the relation between the index versus integral flux is unprecedented in the TeV regime. Prior to the results presented here, spectral variability has been detected only in two other blazars, Mrk 421 and Mrk 501. For Mrk 421, a clear hardening with increasing flux appeared during the 1999/2000 and 2000/2001 observations performed with HEGRA (Aharonian et al. (2002)) and also during the 2004 observations performed with H.E.S.S. (Aharonian et al. 2005a ). In addition, the Mrk 501 observations carried out with CAT during the strong flares of 1997 (Djannati-Ataï et al. (1999)) and also the recent observation performed by MAGIC in 2005 (Albert et al.  2007) have shown a similar hardening. In both studies, the VHE peak has been observed in the $\nu F_{\nu}$ distributions of the flaring states of Mrk 501.

In this section, the spectral variability during the flares of July 2006 is described in more detail. A zoom on the variation of the integral flux (4-minute binning) for the four nights containing the flares (nights MJD 53944, 53945, 53946, and 53947, called the “flaring period”) is presented in the top panel of Fig. 8. This figure shows two exceptional peaks on MJD 53944 and MJD 53946 which climax respectively at fluxes higher than $2.5\times 10^{-9}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$ and $3.5\times 10^{-9}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$ ($\sim 9$ and $\sim 12$ times the Crab Nebula level above the same energy), both about two orders of magnitude above the quiescent state level.

The variation with time of the photon index is shown in the bottom panel of Fig. 8. To obtain these values, the $\gamma$ excess above $200\,{\rm GeV}$ has been determined for each 4-minute bin. Then, successive bins have been grouped in order to reach a global excess higher than 600 $\gamma$. Finally, the energy spectrum of each data set has been determined in the 0.2–1 TeV energy range, as before (corresponding numbers are summarized in Appendix Table 14). There is no clear indication of spectral variability within each night, except for MJD 53946 as shown in Aharonian et al.  (2009). However, a variability can be seen from night to night, and the spectral hardening with increasing flux level already shown in Fig. 7 is also seen very clearly in this manner.

It is certainly interesting to directly compare the spectral behaviour seen during the flaring period with the hardness of the energy spectrum associated with the quiescent state. This is shown in the right panel of Fig. 7, where black points correspond to the four flaring nights; these were determined in the same manner as for the left panel (see 4.1 for details). A linear fit here yields a slope $\rm d\Gamma$/$\rm d\Phi=(2.8\pm 0.3_{\mathrm{stat}})\times 10^{8}\mathrm{cm}^{2}\mathrm{s}$. The grey cross corresponds to the integral flux and the photon index associated with the quiescent state (derived in a consistent way in the energy range from 0.2–1 TeV), showing a clear rupture with the tendancy at higher fluxes (typically above $10^{-10}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$).

These four nights were further examined to search for differences in the spectral behaviour between periods in which the source flux was clearly increasing and periods in which it was decreasing. For this, the first 16 minutes of the first flare (MJD 53944) are of special interest because they present a very symmetric situation: the flux increases during the first half, and then decreases to its initial level, the averaged fluxes are similar in both parts ($\sim 1.8\times 10^{-9}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$), and the observation conditions (and thus the instrument response) are almost constant — the mean zenith angle of each part being respectively $7.2$ and $7.8$ degrees. Again, the spectra have been determined in the 0.2–1 TeV energy range, giving indices of $\Gamma=3.27\pm 0.12_{\mathrm{stat}}$ and $\Gamma=3.27\pm 0.09_{\mathrm{stat}}$ respectively. To further investigate this question and avoid potential systematic errors from the spectral method determination, the hardness ratios were derived (defined as the ratio of the excesses in different energy bands), using for this the energy (TeV) bands [0.2–0.35], [0.35–0.6] and [0.6–5.0]. For any combination, no differences were found beyond the $1\,\sigma$ level between the increasing and decreasing parts. A similar approach has been applied — when possible — for the rest of the flaring period. No clear dependence has been found within the statistical error limit of the determined indices, which is distributed between 0.09 and 0.20.

Finally, the persistence of the energy cut-off in the differential energy spectrum along the flaring period has been examined. For this purpose, runs were sorted again by increasing flux and grouped into subsets containing at least an excess of 3000 $\gamma$ above $200\,{\rm GeV}$⁵⁵5To be significant, the determination of an energy cut-off needs a greater number of $\gamma$ than for a power law fit.. For the seven subsets found, the energy spectrum has been determined in the 0.2–10 TeV energy range both for a simple power law and a power law with an exponential cut-off. This last hypothesis was found to be favoured systematically at a level varying from 1.8 to 4.6 $\sigma$ compared to the simple power-law and is always compatible with a cut-off in the 1–2 TeV range.

All fluxes in the energy bands of Fig. 8 show a strong variability which is quantified through their fractional r.m.s. $F_{\rm var}$ (which depends on observation durations and their sampling). Measurement errors $\sigma_{i,{\rm err}}$ on each of the $N$ fluxes $\phi_{i}$ of the light curve are taken into account in the definition of $F_{\rm var}$:

where $S^{2}$ is the variance

and where $\sigma^{2}_{\mathrm{err}}$ is the mean square error and $\overline{\phi}$ is the mean flux.

The energy-dependent variability $F_{\rm var}(E)$ has been calculated for the flaring period according to Eq. 4 in all three energy bands. The uncertainties on $F_{\rm var}(E)$ have been estimated according to the parametrization derived by Vaughan et al.  (2003), using a Monte Carlo approach which accounts for the measurement errors on the simulated light curves.

Fig. 9 shows the energy dependence of $F_{\rm var}$ over the four nights for a sampling of 4 minutes where only fluxes with a significance of at least 2 standard deviations were considered. There is a clear energy-dependence of the variability (a null probability of $\sim 10^{-16}$). The points in Fig. 9 are fitted according to a power law showing that the variability follows $F_{\rm var}(E)\propto E^{0.19\pm 0.01}$.

This energy dependence of $F_{\rm var}$ is also perceptible within each individual night. In Table 6 the values of $F_{\rm var}$, the relative mean flux, and the observation duration, are reported night by night for the flaring period. Because of the steeply falling spectra, the low energy events dominate in the light curves. This lack of statistics for high energy prevents to have a high fraction of points with a significance more than 2 standard deviation in light curves night by night for the three energy bands previously considered. On the other hand, the error contribution dominates and does not allow to estimate the $F_{\rm var}$ in all these three energy bands. Therefore, only two energy bands were considered: low (0.2–0.5 TeV) and high (0.5–5 TeV). As can be seen in Table 6 also night by night the high-energy fluxes seem to be more variable than those at lower energies.

While $F_{\rm var}$ characterizes the mean variability of a source, the shortest doubling/halving time (Zhang et al.  1999) is an important parameter in view of finding an upper limit on a possible physical shortest time scale of the blazar.

If $\Phi_{i}$ represents the light curve flux at a time $T_{i}$, for each pair of $\Phi_{i}$ one may calculate $T_{2}^{i,j}=|\Phi\Delta T/\Delta\Phi|$, where $\Delta T$=$T_{j}$-$T_{i}$, $\Delta\Phi$=$\Phi_{j}$-$\Phi_{i}$ and $\Phi=(\Phi_{j}+\Phi_{i})/2$. Two possible definitions of the doubling/halving are proposed by Zhang et al.  (1999): the smallest doubling time of all data pairs in a light curve ($T_{2}$), or the mean of the 5 smallest $T_{2}^{i,j}$ (in the following indicated as $\tilde{T_{2}}$). One should keep in mind that, according to Zhang et al.  (1999), these quantities are ill defined and strongly depend on the length of the sampling intervals and on the signal-to-noise ratio in the observation.
This quantity was calculated for the two nights with the highest fluxes, MJD 53944 and MJD 53946, considering light curves with two different binnings (1 and 2 minutes). Bins with flux significances more than $2\,\sigma$ and flux ratios with an uncertainty smaller than 30% were required to estimate the doubling time scale. The uncertainty on $T_{2}$ was estimated by propagating the errors on the $\Phi_{i}$, and a dispersion of the 5 smallest values was included in the error for $\tilde{T_{2}}$.

In Table 7, the values of ${T_{2}}$ and $\tilde{T_{2}}$ for the two nights are shown. The dependence with respect to the binning is clearly visible for both observables. In this table, the last column shows that the fraction of pairs in the light curves which are kept in order to estimate the doubling times is on average $\sim$45%. Moreover, doubling times $T_{2}$ and $\tilde{T_{2}}$ have been estimated for two sets of pairs in the light curves where $\Delta\Phi$=$\Phi_{j}$-$\Phi_{i}$ is increasing or decreasing respectively. The values of the doubling time for the two cases are compatible within $1,\sigma$, therefore no significant asymmetry has been found.

It should be noted that these values are strongly dependent on the time binning and on the experiment’s sensitivity, so that the typical fastest doubling timescale should be conservatively estimated as being less than $\sim 2\,{\rm min}$, which is compatible with the values reported in Aharonian et al.  (2007a) and in Albert et al.  (2007), the latter concerning the blazar Mrk 501.

Time lags between light curves at different energies can provide insights into acceleration, cooling and propagation effects of the radiative particles.

The Discrete Correlation Function (DCF) as a function of the delay (White & Peterson (1984); Edelson & Krolik (1988)) is used here to search for possible time lags between the energy-resolved light curves. The uncertainty on the DCF has been estimated using simulations. For each delay, $10^{5}$ light curves (in both energy bands) have been generated within their errors, assuming a Gaussian probability distribution. A probability distribution function (PDF) of the correlation coefficients between the two energy bands has been estimated for each set of simulated light curves. The r.m.s. of these PDF are the errors related to the DCF at each delay. Fig. 10 shows the DCF between the high and low energy bands for the four-night flaring period (with 4 minute bins) and for the second flaring night (with 2 minute bins). The gaps between each 28 minute run have been taken into account in the DCF estimation.

The position of the maximum of the DCF has been estimated by a Gaussian fit, which shows no time lag between low and high energies for either the 4 or 2 minute binned light curves. This sets a limit of $14\pm 41\,{\rm s}$ from the observation of MJD 53946. A detailed study on the limit on the energy scale on which quantum gravity effects could become important, using the same data set, are reported in Aharonian et al.  (2008a).

Having defined the shortest variability time scales, the nature of the process which generates the fluctuations is investigated, using another estimator: the excess r.m.s.. It is defined as the variance of a light curve (Eq. 5) after subtracting the measurement error:

Fig. 11 shows the correlation between the excess r.m.s. of the light curve and the flux, where the flux here considered are selected with an energy threshold of $200\,{\rm GeV}$. The excess variance is estimated for 1- and 4-minute binned light curves, using 20 consecutive flux points $\Phi_{i}$ which are at least at the $2\,\sigma$ significance level (81% of the 1 minute binned sample). The correlation factors are $r_{1}=0.60^{+0.21}_{-0.25}$ and $r_{4}=0.87^{+0.10}_{-0.24}$ for the 1 and 4 minute binning, excluding an absence of correlation at the $2\,\sigma$ and $4\,\sigma$ levels respectively, implying that fluctuations in the flux are probably proportional to the flux itself which is a characteristic of lognormal distributions (Aitchinson & Brown (1963)). This correlation has also been investigated extending the analysis to a statistically more significant data set including observations with a higher energy threshold in which the determination of the flux above $200\,{\rm GeV}$ requires an extrapolation (grey points in the top panel in Fig. 8). In this case the correlations found are compatible ($r_{1}=0.78^{+0.12}_{-0.14}$ and $r_{4}=0.93^{+0.05}_{-0.15}$ for the 1 and 4 minute binning, respectively) and also exclude an absence of correlation with a higher significance ($4\,\sigma$ and $7\,\sigma$, respectively).

Such a correlation has already been observed for X-rays in the Seyfert class AGN (Edelson et al. (2002), Vaughan, Fabian & Nandra (2003), Vaughan et al.  2003, M^cHardy et al. (2004)) and in X-ray binaries (Uttley & M^cHardy (2001), Uttley (2004), Gleissner et al. (2004)), where it is considered as evidence for an underlying stochastic multiplicative process (Uttley, M^c Hardy & Vaughan (2005)), as opposed to an additive process. In additive processes, light curves are considered as the sum of individual flares “shots” contributing from several zones (multi-zone models) and the relevant variable which has a Gaussian distribution (namely Gaussian variable) is the flux. For multiplicative (or cascade) models the Gaussian variable is the logarithm of the flux. Hence, this first observation of a strong r.m.s.-flux correlation in the VHE domain fully confirms the log-normality of the flux distribution presented in Section 3.2.

For practical reasons, simulated values of $\ln\Phi$ were calculated from Fourier series, thus with a discrete PSD. The fundamental frequency $\nu_{0}=1/T_{0}$, which corresponds to an elementary bin $\delta\nu\equiv\nu_{0}$ in frequency, must be much lower than $1/T$ if $T$ is the duration of the observation. The ratio $T_{0}/T$ was chosen to be of the order of 100, in such a way that the influence of a finite value of $T_{0}$ on the average variance of a light curve of duration $T$ would be less than about 2%. Taking $T_{0}=9\times 10^{5}\,\mathrm{s}$, this condition is fulfilled for the following studies. With this approximation, the simulated flux logarithms are given by:

where $N_{\mathrm{max}}$ is chosen in such a way that $T_{0}/N_{\mathrm{max}}$ is less than the time interval between consecutive measurements (i.e., the sampling interval). According to the definition of a Gaussian random process, the phases $\varphi_{n}$ are uniformly distributed between 0 and 2$\pi$ and the Fourier coefficients $a_{n}$ are normally distributed with mean 0 and variances given by $P(\nu)\,\delta\nu/2$ with $\nu=n\,\delta\nu=n\,\nu_{0}$.
From the long simulated time-series, light curve segments were extracted with the same durations as the periods of continuous data taking and with the same gaps between them. The simulated fluxes were further smeared according to measurement errors, according to the observing conditions (essentially zenith angle and background rate effects) in the corresponding data set.

For a fixed PSD, characterized by a set of parameters $\left\{\alpha,K\right\}$, 500 light curves were simulated, reproducing the observing conditions of the flaring period (MJD 53944–53947), according to the procedure explained in Section 6.1.

For each set of simulated light curves, segments of 20 minutes duration sampled every minute (and alternatively segments of 80 minutes duration sampled every 4 minutes) were extracted and, for each of them, the excess r.m.s. $\sigma_{xs}$ and the mean flux $\bar{\Phi}$ were calculated as explained in Section 5.4. For a large range of values of $\alpha$ and $K$, simulated light curves reproduce well the high level of correlation found in the measured light curves. On the other hand, the fractional variability $F_{\mathrm{var}}$ and $\bar{\Phi}$ are essentially uncorrelated and will be used in the following. A likelihood function of $\alpha$ and $K$ was obtained by comparing the simulated distributions of $F_{\mathrm{var}}$ and $\bar{\Phi}$ to the experimental ones. An additional observable which is sensitive to $\alpha$ and $K$ is the fraction of those light curve segments for which $F_{\mathrm{var}}$ cannot be calculated because large measurement errors lead to a negative value for the excess variance. The comparison between the measured value of this fraction and those obtained from simulations is also taken into account in the likelihood function. The $95\%$ confidence contours for the two parameters $\alpha$ and $K$ obtained from the maximum likelihood method are shown in Fig. 12 for both kinds of light curve segments. The two selected domains in the $\left\{\alpha,K\right\}$ plane have a large overlap which restricts the values of $\alpha$ to the interval (1.9, 2.4).

Another method for determining $\alpha$ and $K$ is based on Kolmogorov structure functions (SF). For a signal $X(t)$, measured at $N$ pairs of times separated by a delay $\tau$, $\left\{t_{i},t_{i}+\tau\right\}\,(i=1,...,N)$, the first-order structure function is defined as (Simonetti et al. (1985)):

In the present analysis, $X(t)$ represents the variable whose PSD is being estimated, namely $\ln\Phi$. The structure function is a powerful tool for studying aperiodic signals (Rutman (1978), Simonetti et al. (1985)), such as the luminosity of blazars at various wavelengths. Compared to the direct Fourier analysis, the SF has the advantage of being less affected by “windowing effects” caused by large gaps between short observation periods in VHE observations. The first-order structure function is adapted to those PSDs whose variability spectral index is less than 3 Rutman (1978), which is the case here, according to the results of the preceding section.

Fig. 13 shows the first-order SF estimated for the flaring period (circles) for $\tau<60$ hr. At fixed $\tau$, the distribution of $\mathrm{SF}(\tau)$ expected for a given set of parameters $\left\{\alpha,K\right\}$ is obtained from 500 simulated light curves. As an example, taking $\alpha=2$ and $\log_{10}(K/{\rm Hz}^{-1})=2.8$, values of $\mathrm{SF}(\tau)$ are found to lie at 68% confidence level within the shaded region in Fig. 13.

In the case of a power law PSD with index $\alpha$, the SF averaged over an ensemble of light curves is expected show a variation as $\tau^{\alpha-1}$ (Kataoka et al. (2001)). However, this property does not take into account the effect of measurement errors, nor of the limited sensitivity of Cherenkov telescopes at lower fluxes. For the present study, it was preferable to use the distributions of $\mathrm{SF}(\tau)$ obtained from realistic simulations including all experimental effects. Using such distributions expected for a given set of parameters $\left\{\alpha,K\right\}$, a likelihood function can be obtained from the experimental SF and further maximized with respect to these two parameters. Furthermore, the likelihood fit was restricted to values of $\tau$ lower than $10^{4}\,\mathrm{s}$, for which the expected fluctuations are not too large and are well-controlled. The 95% confidence region in the $\left\{\alpha,K\right\}$ plane thus obtained is indicated by the dotted line in Fig. 12. It is in very good agreement with those based on the excess r.m.s.–flux correlation and give the best values for $\alpha$ and $K$:

The variability index $\alpha$ at VHE energies is found to be remarkably close to those measured in the X-ray domain on PKS 2155$-$304, Mrk 421, and Mrk 501 (Kataoka et al. (2001)).

Simulations were also used to investigate if the estimator $T_{2}$ can be used to constrain the values of $\alpha$ and $K$. However, for $\alpha$ less than 2, no significant constraints on those parameters are obtained from the values of $T_{2}$. For higher values of $\alpha$, doubling times only provide loose confidence intervals on $K$ which are compatible with the values reported above. This can be seen in Table 8, showing the 68% confidence intervals predicted for $T_{2}$ and $\tilde{T_{2}}$ for a lognormal process with $\alpha$=2 and $\log_{10}(K/(\rm Hz^{-1}))=2.8$, as obtained from simulation. Therefore, the variability of PKS 2155$-$304 during the flaring period can be consistently described by the lognormal random process whose PSD is characterized by the parameters given by Eq. 9.