A&A 459, 453-464 (2006)
DOI: 10.1051/0004-6361:20065884

Mrk 421, Mrk 501, and 1ES 1426+428 at 100 GeV with the CELESTE Cherenkov telescope

D. A. Smith¹ - E. Brion^1, - R. Britto² - P. Bruel³ - J. Bussons Gordo² - D. Dumora¹ - E. Durand¹ - P. Eschstruth⁴ - P. Espigat⁵ - J. Holder⁴ - A. Jacholkowska² - J. Lavalle² - R. Le Gallou¹ - B. Lott¹ - H. Manseri³ - F. Münz⁵ - E. Nuss² - F. Piron² - R. C. Rannot^1,6 - T. Reposeur¹ - T. Sako³

1 - Centre d'études nucléaires de Bordeaux Gradignan, CENBG UMR 5797 CNRS/IN2P3, Université Bordeaux 1, Chemin du Solarium, BP 120, 33175 Gradignan Cedex, France
2 - Laboratoire de physique théorique et astroparticules (UMR 5207 CNRS/IN2P3), Université Montpellier 2, 34095, France
3 - Laboratoire Leprince-Ringuet (UMR 7638 CNRS/IN2P3), École Polytechnique, 91128 Palaiseau, France
4 - Laboratoire de l'accélérateur linéaire (UMR 8607 CNRS/IN2P3), Université Paris 11, 91898 Orsay, France
5 - Laboratoire d'astroparticule et cosmologie (UMR 7164 CNRS/IN2P3), Collège de France, 75231 Paris, France
6 - Astrophysical Sciences Division, Bhabha Atomic Research Centre, Mumbai-400 085, India

Received 22 June 2006 / Accepted 11 August 2006

Abstract
We have measured the gamma-ray fluxes of the blazars Mrk 421 and Mrk 501 in the energy range between 50 and 350 GeV (1.2 to $8.3 \times 10^{25}$ Hz). The detector, called CELESTE, used first 40, then 53 heliostats of the former solar facility "Thémis'' in the French Pyrénées to collect Cherenkov light generated in atmospheric particle cascades. The signal from Mrk 421 is often strong. We compare its flux with previously published multi-wavelength studies and infer that we are straddling the high energy peak of the spectral energy distribution. The signal from Mrk 501 in 2000 was weak ($3.4\sigma$). We obtain an upper limit on the flux from 1ES 1426+428 of less than half that of the Crab flux near 100 GeV. The data analysis and understanding of systematic biases have improved compared to previous work, increasing the detector's sensitivity.

Key words: galaxies: BL Lacertae objects: individual: Mrk 421, 1ES 1426+428, Mrk 501 - gamma rays: observations - atmospheric effects

1 Introduction

Measurements of high energy emission from active galactic nuclei (AGN) give insights into a variety of open problems, such as the nature of the AGNs themselves, and the density and evolution of the extragalactic diffuse infrared background. After GLAST is launched in 2007 (Gehrels & Michelson 1999), the large number of high galactic latitude GeV gamma-ray sources to be seen will make AGNs become a background for searches for new classes of emitters - insufficient understanding of high energy AGNs might just limit potential glimpses of a hidden Universe. Several AGNs of the "blazar'' class emit above 250 GeV and have been detected by atmospheric Cherenkov imaging telescopes. They are called HBLs, for high energy BL Lac objects. Of these, Mrk 421 (Bazejowski et al. 2005; Piron et al. 2001a; Aharonian et al. 2003) and Mrk 501 (Dwek & Krennrich 2005; Djannati-Ataï et al. 1997) are the brightest.

The Synchrotron Self Compton model ("SSC'') is consistent with much of the multiwavelength blazar data. For discussion of these models as applied to the blazars discussed in this paper, see for example (Tavecchio et al. 1998) or (Katarzynski et al. 2001). In its simplest variant, relativistic electrons generate synchrotronphotons with a $\nu F_\nu$ peak in the optical to X-ray range, and then up-scatter these same photons to produce the peak seen in the GeV range. The inverse Compton peak position can constrain the physical parameters of the model. Many, if not most, observed characteristics of blazars remain unexplained. Examples are the relation between variations at different wavelengths, as highlighted by the "orphan'' gamma ray flares observed for 1ES1959+650 (Krawczynski et al. 2004), and the role that proton acceleration might play. Progress requires simultaneous measurements across the spectrum.

The CELESTE telescope was one of the first instruments sensitive near the Compton peak for these blazars, beyond the energy range of the relatively large sample of EGRET blazars (Mukherjee et al. 1997; Hartman et al. 1999). Mrk 421 was studied with EGRET, but Mrk 501 was visible only during a flare (Kataoka et al. 1999), while 1ES 1426+428 was undetected by EGRET, although well-measured at TeV energies (e.g. Djannati-Ataï et al. 2002; Horan et al. 2002). A key motivation for our experiment was to provide blazar measurements in an unexplored energy window.

This paper presents the first observations of three blazars below 100 GeV, by a solar heliostat array, and is structured as follows: we describe the detector, the analysis method, and the data samples. We emphasize improvements in our knowledge of the systematic uncertainties. We present light curves for the three blazars as well as the integrated fluxes, and place the results in the context of SSC models.

2 The CELESTE gamma-ray telescope

The CELESTE detector is located in the eastern French Pyrénées (N. $42.50^{\circ}$, E. $1.97^{\circ}$, altitude $1650~{\rm m}$) and described in Giebels et al. (1998), Paré et al. (2002) and in de Naurois et al. (2002). Figure 1 illustrates the principle. CELESTE was dismantled in June, 2004.

CELESTE detected the Crab nebula in 1998 using a solar heliostat array (Smith 1998). STACEE, very similar to CELESTE, reported the flux above 190 GeV (Oser et al. 2001), followed by the flux above 60 GeV from CELESTE (de Naurois et al. 2002). STACEE recently measured the integral flux above 140 GeV from Mrk 421 (Boone et al. 2002) and, assuming power law spectral indices extrapolated from higher energies, bracketed the 50 to 700 GeV differential spectrum. A review of solar tower gamma-ray telescopes can be found in (Smith 2005). Here we recall key features of the CELESTE detector.

	Figure 1: Experimental principle: Heliostats tracking a source reflect Cherenkov light from atmospheric particle cascades to the secondary optics, photomultipliers, and acquisition electronics located high in the 100 meter tall tower.
Open with DEXTER

CELESTE used 53 heliostats (40 until 2001) of the Thémis former solar facility. The CAT Cherenkov imager ran on the same site (Piron et al. 2001a and references therein), as did the THEMISTOCLE (Baillon et al. 1993) and ASGAT (Goret et al. 1993) Cherenkov sampling arrays. Each heliostat has a $54~{\rm m}^2$ mirror on an alt-azimuth mount. Light from the heliostats was reflected to a secondary optical system at the top of the 100 m tower, where one photomultiplier assembly (PMT) viewed each heliostat. Heliostat control and data acquisition systems were also high in the tower.

Data acquisition was triggered when the summed PMT pulse height for at least M of the N groups of heliostats exceeded 4.5 photoelectrons ("p.e.'') per heliostat. For the 40 heliostat array, M=3 of N=5 groups with 8 heliostats each was used most often. For the 53 heliostat array there were N=6 trigger groups of between 6 and 8 heliostats each, and the trigger multiplicity was set at M=3, 4, or 5 for different observations. The deadtime in the electronic delays depends on the rates, which depend on the threshold settings, which are determined in part by the choice of M. We took pains to remove these biases (Manseri 2004), in particular, with an offline group multiplicity cut based on the scaler rates. The acceptance shown in Fig. 5 has an offline multiplicity cut of 4 of 5, for illustration.

Each PMT signal was digitized at 1.06 nanoseconds per sample, with about 3 digital counts (dc) per photoelectron. 100 samples centered around the nominal Cherenkov pulse arrival time were recorded, providing event-by-event pedestal information that we used to determine the night sky background light for each channel. A second 28 sample window contained a timing reference pulse. Scaler rates, PMT anode currents, a GPS time stamp, and some meteorological information rounded out the data record.

3 Monte Carlo simulations

CELESTE's first publication, on the flux from the Crab nebula quoted a $\sim$30 GeV threshold at the trigger level, and 60 GeV after analysis (de Naurois et al. 2002). We estimated the uncertainty on the energy scale to be less than 30%, dominated by disagreement between predictions of Monte Carlo atmospheric cascade generators, and disagreement between estimated and observed PMT illuminations induced by bright stars. Since then we have corrected some small errors in the KASCADE Monte Carlo (Kertzman & Sembroski 1994) and the agreement is improved ( $\delta \epsilon_{\rm MC} = 5$ to 10%). The CORSIKA Monte Carlo atmospheric shower generator (Heck et al. 1998) best reproduces our experimental observables, and we use it for the results presented here.

Cherenkov photons from the shower generator are fed into a detailed simulation of our optics and electronics, to develop our data analysis methods and cuts, as well as to calculate the energy dependent effective area A(E). Here we describe improvements made to the detector simulation since de Naurois et al. (2002). The acceptances will be discussed after the analysis cuts have been explained. Simulations indicate that the trigger rate is dominated by 18 Hz of protons, with 4 Hz of helium nuclei. Photons from the Crab nebula trigger the detector at about 2% of the trigger rate.

3.1 Optical throughput

Previously, the phototube illuminations when viewing stars were lower than predicted by the simulation. We reviewed the ray tracing programs, and found two factors that dominated the discrepancy. First, the combined heliostat and secondary mirror reflectivities are about 20% worse than earlier measurements indicated. Second, the simulated heliostat focussing was sharper than in reality. We compared the image sizes obtained from star "scans'' used for heliostat alignment with simulated sizes, and modified the simulation to improve agreement. (In scans, we record PMT anode currents as the heliostats are pointed successively at a grid of points around the star position: see Fig. 7 in Paré et al. (2002).)

Figure 2 illustrates a consequence of the ray tracing modifications (Brion 2004). When tracking Mrk 421, the star 51 UMa (HD 95934, blue magnitude $M_{\rm B}=6.16$) falls within the $\pm$5 mrad optical field-of-view for heliostats near the center of the array. The On-source PMT anode currents are thus larger ( $i_{\rm On} \sim 13~{\rm\mbox{\usefont{U}{eur}{m}{n}\char22}A}$) than the Off-source currents which are due only to the diffuse night sky light ( $i_{\rm Off}\sim 10~{\rm\mbox{\usefont{U}{eur}{m}{n}\char22}A}$). Using the PMT gains g, the current differences yield the star-induced illuminations $b = (i_{\rm On}-i_{\rm Off})/ge$, in photoelectrons per second (e is the electron charge). The figure shows that measured illuminations vary with the pointing direction, but less than was predicted by the original simulation. The simulations include the stars in the On and Off fields that contribute at least 5% as much as 51 UMa, that is, $M_{\rm B} < 9.4$. We remark that the gamma-ray energy scale predicted by the Monte Carlo simulation shifts less than Fig. 2 suggests, since an atmospheric shower is a light source of angular extent comparable to the detector field-of-view, and so the effect of heliostat "defocussing'' is less violent for showers than it is for the point-like stars.

A few observations constrain the uncertainty in the optical throughput. Photometry using star scans yields 20% channel-to-channel dispersion compared to the simulated values. A census of damaged heliostat mirrors, and the phototube manufacturer's measurements of the photocathode responses, combined together account for half of the dispersion. We attribute the rest to errors in the gains, and focussing anomalies for individual heliostats. However, what matters most is the sums over the 8 heliostats of the trigger groups: the total variation over the 5 groups is $\pm$15%, and the rms is less than 10%. We conclude that the overall effect on the gamma-ray energy scale of the uncertainty in the optical throughput $\delta \epsilon_{\rm opt}$ is less than 10%.

Figure 2: Example of phototube illuminations induced by the star 51 UMa in the field of view of Mrk 421, versus hour angle, for heliostat C07. Squares & crosses: For the 40 pairs of On- and Off-source "double-pointing'' data runs passing current and rate stability criteria, the quantity shown is the difference of the average current during each run divided by the PMT gain g, $(i_{\rm On}-i_{\rm Off})/g$. For the crosses, the trigger rates were stable but lower than usual: optical transmission was low, presumably due to clouds or frost. Circles: original detector simulation. Stars: Simulation results, after "defocussing'' the heliostats and degrading the mirror reflectivities as described in the text. The curves are to guide the eye.

Open with DEXTER

Since de Naurois et al. (2002), the electronics simulation has been completely re-written, leading to improved agreement between predicted and observed signals, in the trigger, and for single channels. The good agreement for the final analysis variable, described below, is the result of the care taken to model the elements of the complex electronics chain.

3.2 Atmospheric transmission

The amount of Cherenkov light that reaches the ground varies between different sites and seasons depending mainly on the aerosol content versus altitude (Bernlöhr 2000). The same also holds for atmospheric extinction in stellar photometry (Hayes & Latham 1975). The effect on CELESTE is particularly striking: the $\sim$25 Hz cosmic ray trigger rates for winter sources (e.g. Crab and Mrk 421) decrease to under 15 Hz every year in spring, affecting sources like Mrk 501 and 1ES 1426+428. Table 1 lists the fraction of the data sets above and below 15 Hz.

The CORSIKA Monte Carlo provides a choice of tabulated wavelength and altitude dependent extinction curves. The default contains some aerosols. The Palomar curve from Hayes & Latham (1975) used for our stellar photometry studies contains fewer aerosols, for an integrated transmission above Thémis at 400 nm 6% greater than for the default CORSIKA curve.

An amateur telescope with filters and a CCD camera was used to measure the total integrated atmospheric transmission during some nights when the trigger rates were above 20 Hz, using the apparent brightness of several stars as a function of their zenith angles (Bouguer method), at three wavelengths (blue: 390-490 nm, green: 500-560 nm, and red: 610-660 nm). The Cherenkov light recorded by CELESTE is mainly in the blue range. The CCD results show that the sky at Thémis corresponds to the Palomar extinction curve. We made a CORSIKA extinction table corresponding to the Palomar curve which we use to simulate good running conditions.

Table 1: Number of data pairs, after applying stability criteria to the PMT anode currents and to the trigger rates. Single Pointing (SP), Doubling Pointing (DP), and Single Pointing with Veto (SPV) are explained in the text. The fraction of the data set with a cosmic ray trigger rate below 15 Hz, which we ascribe to seasonal increases in atmospheric extinction, is also listed.

To calculate the acceptances for data acquired with trigger rates < $15~{\rm Hz}$ we made another CORSIKA extinction table that we call "dusty''. We increased the aerosol contribution to the Palomar curve so that the integrated transmission above Thémis at 400 nm is 60% less. This choice comes from two observations. The first appears in Fig. 2, where the phototube illumination due to the bright star near Mrk 421 in good conditions is halved for data acquired with a low trigger rate. The second is that proton shower simulations using the nominal and "dusty'' extinction tables reproduced the experimental nominal and degraded trigger rates, respectively. Figure 5, discussed below, shows the acceptance curves obtained for the two atmospheres, used to convert our gamma-ray rates to fluxes.

4 Data samples

Maximum Cherenkov emission for 100 GeV $\gamma$-ray showers occurs around 11 km above the site. For CELESTE's oldest data, the heliostats were aimed at this height in the direction of the source, to maximize Cherenkov light collection. This is called Single Pointing data (SP).

SP provides the lowest energy threshold, but also decreases the detector's sensitive area. For one season we took data in a configuration called Double Pointing (DP), where half the heliostats point at 11 km and the rest at 25 km. This increases the area and improves cosmic ray rejection, while raising the energy threshold only slightly.

Our $\pm$5 mrad field-of-view matches the apparent size of a gamma shower, to optimize the ratio of Cherenkov light to night sky background light. Tails of cosmic ray showers are outside the field-of-view, reducing background at the trigger level but weakening offline rejection. In 2001 we added 13 heliostats to CELESTE, and thereafter pointed 41 heliostats at the 11 km point, and the remaining 12 at a ring 150 m around the point. The 12 "veto'' heliostats can see shower tails outside the central field-of-view, enhancing cosmic ray rejection (Manseri 2004). We call this Single Pointing with Veto data (SPV).

Ultimately the performance of the analysis was such that the veto rejection was not used. But the 12 veto heliostats around the edge of the heliostat field led us to restructure the trigger into M=6 groups. The changed trigger topology affects the shape of the recorded showers, leading to different energy dependent acceptances.

Source observations ("On'') last about 20 min and are followed or preceded by an observation at the same declination, offset by 20 min in right ascension. "Off'' data gives a reference for the cosmic ray background - the signal is the difference between On and Off, after analysis cuts. However, changes in sky conditions between On and Off can cause differences having nothing to do with gamma-rays: we select stable On-Off pairs with a cut on the $\chi^2$ value obtained from a least-squares fit to a constant PMT anode current versus time, and to the trigger rate versus time, for both runs. We also rejected pairs with very low data rates (< $10~{\rm Hz}$), or anomalously high trigger dead times. The weather at Thémis is fickle, and we rejected over half of our data.

Table 1 lists the data sets, recorded during clear, moonless nights from December 1999 through May 2004, after data selection.

5 Data analysis

5.1 Event reconstruction and background rejection

Since de Naurois et al. (2002) we have changed our analysis strategy, as described below. The Flash ADC data quality improved, as did the channel-to-channel timing offsets, which both help wavefront reconstruction. Our procedures to remove the biases induced by night sky background differences at the trigger level ("software trigger'') and in the FADC data ("padding'') have been refined. Overall sensitivity has more than doubled.

Just above the trigger threshold, the Cherenkov signal in one channel is comparable to the fluctuations of the night sky background light. Summing the channels greatly improves the signal-to-noise ratio, after correcting for the propagation delay to each heliostat, assuming a spherical Cherenkov wavefront. Lacking knowledge of the shower core position, we scan the plane at 11 km and evaluate the ratio H/W for each point in a grid, where H and W are the height and the width of the summed signal (see Fig. 3). The grid position giving the largest value of the ratio H/W is a measure of the shower core position. The wavefront sphericity correlates with how much H/W decreases 200 m away from the shower core. We measure the relative decrease with a variable called $\xi = (H/W)_{200~{\rm m}}/(H/W)_{\rm max}$, shown in Fig. 4 for an Off observation, Mrk 421 On-Off data, and for a simulation of a $\gamma$-ray spectrum. The Cherenkov wavefront from $\gamma$-ray induced showers is, on average, more spherical than for showers initiated by charged cosmic rays and, as expected, $\gamma$-ray showers have smaller $\xi$ than showers initiated by charged cosmic rays. The method and preliminary results were shown in Bruel (2004) and detailed in Manseri (2004).

The $\xi$ distributions change for different heliostat pointings. This article concerns three blazars with similar declinations, all transiting near zenith at Thémis, so the only changes are the pointing configuration and the hour angle. We studied simulations and the strong signals from Mrk 421 and the Crab to explore the $\xi$ cuts (Brion 2005). Statistical significance (sensitivity) is optimized by cutting a bit above the peak in the gamma ray $\xi$ distribution, e.g. $\xi < 0.35$ in Fig. 4. For other pointings, the peak lies between $\xi \simeq 0.25$ and $\xi \simeq 0.35$. The acceptance depends in part on the efficiency of the $\xi$cut, obtained from the Monte Carlo. This efficiency is related to the integral of the $\xi$ distribution, a steep function of $\xi$ near the peak. Thus we choose to cut beyond the peak, at $\xi < 0.4$, for all data sets, which increases the statistical uncertainty, but decreases the acceptance uncertainty $\delta \epsilon_{acc}$ to less than 10%, and is simpler than having an array of cut values for the various pointings.

5.2 Acceptances and performance

For a gamma-ray source with differential spectrum $\phi(E)$ the number N of detected gammas will be

$$N = T\int_0^\infty A(E) \phi (E) {\rm d}E$$

where T is the On-source time and A(E) is the energy dependent gamma-ray collection area obtained from the Monte Carlo calculations. We have calculated A(E) for each heliostat configuration, and for hour angles of (0, 1.0, 1.5, 2.0) hours, for each of the two atmospheres described above. Figure 5 shows A(E) at the trigger level, and after cuts, for the direction of the blazar Mrk 421 ( $\delta=38\hbox{$^\circ$ }$), one hour after transit. (The few degree discrepancy with the other two blazar's declinations induces a negligible error). For the blazar spectra we assume a power law, $\phi(E)=kE^{-\alpha}$. CELESTE's energy range is near the SED high energy peak for the three blazars, where $\alpha \equiv 2$. We vary $\alpha$ between 1.8 and 2.2 to estimate the uncertainty due to this choice. The energy reconstruction method in Piron et al. (2003) was used to search for a signal in energy bands in Lavalle et al. (2006) but was not applied here.

The combined uncertainties $(\delta\epsilon_{\rm MC} \otimes \delta\epsilon_{\rm opt})< 20$% affect the energy scale. Furthermore, along with $\delta\epsilon_{\rm acc}$ they lead to a flux error. A remaining uncertainty comes from the varying atmospheric transmission. We make a first-order correction by grouping the data into subsets with trigger rates of $\leq$15 and >15 Hz (see Table 1) for use with the nominal or "dusty'' acceptances, A(E), as in Fig. 5. The maximum counting rate after cuts is nearly halved with the "dusty'' atmosphere. The observed extrema of the trigger rates are 12 and 24 Hz, but the distribution clusters around 14 Hz and 22 Hz. We conclude that the overall systematic flux uncertainty is <25%.

To summarize CELESTE's performance: for one hour of On-source Crab data near transit, analysis using $\xi < 0.4$ gives $3.7 \pm 0.2$ gamma/min over the whole data set. The significance is $3.3\sigma$ for SP and DP data, and $5\sigma$ for SPV. The higher sensitivity for the more recent data is due mainly to having six trigger groups instead of five, at the cost of a somewhat higher minimum energy: Table 2 shows that the SP raw data rate is 15% higher than for SPV, but that analysis cuts reject more background for SPV data (90%) than for SP data (83%). The higher sensitivity also stems in part from improvements in FADC quality, timing adjustments, and gain settings over the years. Optimized cuts increase the significances by 20%, and more for Mrk 421 (Brion 2005) but are not used in this paper.

Folding our energy dependent acceptance A(E) (similar to Fig. 5) with the Crab inverse Compton spectrum parametrized in Table 6 of Aharonian et al. (2004) yields a maximum at 90 GeV. The rate falls to 20% of the maximum below 50 GeV and above 350 GeV. Figure 6 shows the acceptance-corrected pair-by-pair Crab fluxes, as a function of date and hour angle. The stability is adequate to search for blazar variability. The power flux at 90 GeV is $\nu F_\nu = 0.89 \pm 0.05 \times 10^{-10}$ erg cm^-2 s^-1, higher than our previous result in de Naurois et al. (2002) but within the uncertainties, and in agreement with the spectrum in Aharonian et al. (2004).

Figure 3: Left: summed Cherenkov pulse - a 100 ns Flash ADC window, centered at the time an ideal spherical wavefront 11 km directly above the site would reach the photomultiplier tubes after reflection from the heliostats, is recorded for each channel. In analysis, we sum the 40 (SP and DP data) or 41 (SPV data) channels repeatedly over a grid of hypothetical shower core positions. Right: the H/W values for each position of the grid, for a single simulated event. The grid position with $(H/W)_{\rm max}$ is a good ($\pm 15$ m) estimator of where the gamma-ray would have hit the ground. The analysis variable $\xi = (H/W)_{200~{\rm m}}/(H/W)_{\rm max}$ uses $(H/W)_{200~{\rm m}}$, which is the value averaged over a ring 200 meters from the grid center. The vertical axis is in ADC counts per nanosecond.

Open with DEXTER

	Figure 4: Left: $\xi$ distributions for Mrk 421 for March 2004 data. Black dots show On-Off data. The broad histogram is Off data (i.e. cosmic ray background). The narrow histogram is simulated gamma-rays. The simulation and Off data are normalized to the number of On-Off events. Right: statistical significance of the excess as a function of the $\xi$ cut value.
Open with DEXTER

Figure 5: Left: energy dependent gamma-ray collection area A(E), at the trigger level (round markers), and after analysis cuts (triangles), for sources culminating near zenith, 1 h after transit, for 11 km single pointing (SP). The curves are splines. The two solid curves used the low-aerosol atmosphere while the two dashed curves used the "dusty'' atmosphere (see text). The analysis trigger group multiplicity cut for these plots was M=4 of the N=5 groups. Right: A(E) multiplied by a typical E-2 flux spectrum, to give the relative counting rate A(E)/E2. Integral fluxes are quoted at the energy corresponding to the maximum rate, while the energy range of 50 to 350 GeV corresponds to 20% of the maximum rate.

Open with DEXTER

Table 2: Number of events before and after cuts, for the Mrk 421 data sets. $N_\sigma$ is calculated after correcting for the $\sim$$80\%$ data acquisition efficiency.

6 Results and discussion

Tables 2, 4, and 5 summarize the data analysis results for the three blazars Mrk 421, Mrk 501, and 1ES 1426+428, respectively. Figures 7-9, show the acceptance-corrected light curves for the blazar Mrk 421 during the years 2000, 2001, and 2004 seasons. Figures 12 and 14 show the light curves for Mrk 501 and 1ES 1426+428, respectively. The curves assume power laws with spectral index $\alpha = 2$. There is an often strong signal from Mrk 421, and no signal for 1ES 1426+428. The excess in the Mrk 501 data is discussed below.

Figure 6: Left: Crab fluxes for single data runs over the 5 years that CELESTE ran. SP pointing was mainly used during the first year, DP for the 2nd year, and SPV for the last. The integral flux above 90 GeV is $0.62 \pm 0.04 \times 10^{-9}$ photons/cm2/s averaged over the three good Crab seasons. Right: same, as a function of the hour angle of the data runs. After acceptance corrections the observed flux is steady.

Open with DEXTER

6.1 Mrk 421

Understanding blazars requires observations at different wavelengths, preferably simultaneous given their variability. Happily, much of our data coincides with previously published results. Table 3 lists the CELESTE, TeV, and X-ray fluxes for the epochs with CELESTE data. CELESTE's richest data taking period was March 2000, with 11 of 14 consecutive nights (epoch (1), from 2000 February 27 to March 12, MJD 51 601.96 to 51 615.09). Figure 7 shows the daily averaged data during this period. CAT (Piron et al. 2001a) and HEGRA (Aharonian et al. 2002) have data for most of the same nights, and the four experiments with data at this epoch are all at a low level. The CELESTE flux averaged over the 2003 season, epoch (3), from February through April (MJD 52 667.12 to 52 762.94) was at an intermediate level, as were the other experiments. Overall, the flux variations observed with CELESTE track the other measurements.

Figure 10 shows the SED for Mrk 421 taken from Bazejowski et al. (2005) to which we have added our result for epoch (2) when TeV data from Cherenkov imagers are available, 4 nights in February 2001 (MJD 51 956.98 to 51 963.11). ASM indicates that the blazar was at the low end of the "medium'' X-ray state as defined in Bazejowski et al. (2005). We have plotted the CAT spectrum for the year 2000 (Piron et al. 2001a), which has the same integral flux as measured on those four nights.

We now look more closely at the light curves. The detailed multi-wavelength study reported in Bazejowski et al. (2005) makes February-March 2004 particularly interesting (epoch 4). CELESTE measured a steady high state of $1.04 \pm 0.1\times 10^{-9}$ cm^-2 s^-1 for the three nights of 2004 February 15 to 18 (MJD 53050.05 to 53053.08, see Table 3 and Fig. 9). VERITAS shows a steady low state, in seeming contrast to CELESTE. But Fig. 1 of (Bazejowski et al. 2005) shows an increasing optical flux on those dates, while PCA and ASM on RXTE show a decreasing X-ray state. This could be a shift of the low energy peak of the SED to lower energy. If the high energy SED peak also shifts during this epoch it could explain that CELESTE sees a high rate when VERITAS does not: VERITAS is above the SED peak while the CELESTE range straddles it.

Similarly, and also in Fig. 9, PCA and ASM show a declining X-ray state followed by a small rebound for 2004 March 16 to 18 (epoch 5). For two nights of data (MJD 53 079.91 to 53 082.02) CELESTE has an average state of $1.4 \pm 0.2\times 10^{-9}$ cm^-2 s^-1, but here VERITAS sees the decline of the flare before returning to the low state. CELESTE's stability during VERITAS' changes is again consistent with Fig. 10, that is, with an SED peak at 100 GeV. Recent spectral work on Mrk 421 by STACEE also seems to indicate that our 50 to 350 GeV energy range matches a high energy SED peak (Carson 2005).

Figure 7: Mrk 421 daily averages in early 2000. The CELESTE data runs from MJD 51572 to 51613 (January 29 to March 10), while epoch (1) is defined in Table 3 and discussed in the text. CAT data is from Piron et al. (2001a, left-hand scale) while the HEGRA points are from Krawczynski et al. (2001), right-hand scale). RXTE ASM (see ASM in references, left-hand scale) and PCA (Krawczynski et al. 2001, right-hand scale) results are also shown.

Open with DEXTER

Figure 8: Mrk 421 daily averages in 2001. Epoch (2) is defined in Table 3 and discussed in the text. CELESTE data runs from MJD 51 957 to 52 015 (February 17 to April 16). STACEE data is from Boone et al. (2002), HEGRA data is from Aharonian et al. (2002), and the VERITAS data is from Holder et al. (2001). The other data points are as in the preceding figure. Note the different left- and right-hand axes.

Open with DEXTER

	Figure 9: Mrk 421 daily averages in February and March 2004. Epochs (4) and (5) are defined in Table 3 and discussed in the text. The RXTE PCA and VERITAS data are taken from Bazejowski et al. (2005), and ASM data from ASM in references.
Open with DEXTER

Table 3: Multiwavelength fluxes from Mrk 421 coincident with CELESTE data. The epoch numbers (1) through (5) are referred to in the text and figures. (a) ASM counts per second (see ASM in references). RXTE PCA 2-10 keV fluxes are in (b) $10^{-10}~{\rm erg~cm^{-2}~s^{-1}}$, 3-20 keV (Krawczynski et al. 2001), and in (c) counts per second (Bazejowski et al. 2005). (d) CAT fluxes are in $10^{-12}~{\rm photons~cm^{-2}~s^{-1}},~E>250~{\rm GeV}$ (Piron et al. 2001a,b). (e) HEGRA fluxes are in $10^{-12}~{\rm photons~cm^{-2}~s^{-1}},~E>1000~{\rm GeV}$ (Krawczynski et al. 2001 for the 2000 data, Aharonian et al. 2003 for 2001). (f) VERITAS fluxes above 300 GeV are in gamma/min, where the Crab gave 2.40 (2.93) gamma/min for 02/03 (03/04) (Holder et al. 2001; Bazejowski et al. 2005). (g) CELESTE fluxes are in $10^{-9}~{\rm photons~cm^{-2}~s^{-1}},~E>100~{\rm GeV}$.

6.2 Mrk 501

The statistical significance of the excess for all the Mrk 501 data is weak ($2.9\sigma$), as shown in Table 4. The same $\xi < 0.4$ cut as for Mrk 421 was applied. The light curve in Fig. 12 shows that the excess comes mostly from the year 2000, with $3.4\sigma$ for those four months. Figure 11 shows that the excess in 2000 ressembles a gamma-ray signal. However, changing the $\xi$ cut values to those that optimize the statistical significance of the Mrk 421 signal do not enhance this excess, suggesting that taking the excess to be a gamma-ray signal is delicate.

CELESTE and the BeppoSAX X-ray satellite have data in common for four nights, from 2000 March 5 to 12, when Beppo measured between $(0.76 \pm 0.08)$ and $(1.18 \pm 0.12) \times 10^{-10}~{\rm erg~cm^{-2}~s^{-1}}$ (Massaro et al. 2004). The ASM rate during this period is a roughly constant $0.4 \pm 0.1$. (ASM rates are unreliable at such low fluxes.) Averaging the ASM count rates for the nights with CELESTE data gives $0.44 \pm 0.06$ and $0.37 \pm 0.09$ for 2000 and 2001, respectively. The dispersion of the points is 0.5 cts/s for both seasons. Averaging ASM over longer periods gives essentially the same result, so that no significant difference in X-ray activity for the two years is apparent.

Figure 6 shows that our experiment was stable. We choose to interpret the excess as due to gamma-rays from this well-known emitter. We use the same acceptances as for Mrk 421, again matched to the conditions of each data run. The flux at 110 GeV thus obtained for the year 2000 is $(0.5\pm 0.15)\times 10^{-10}~{\rm erg~cm^{-2}~s^{-1}}$ (Statistical errors only). Combining both years decreases the result by less than $1\sigma$. Figure 13 shows the flux superimposed on an SED taken from Kataoka et al. (1999). The ASCA X-ray and EGRET gamma ray data were simultaneous with each other. We have added the BeppoSAX X-ray data taken at the same time as CELESTE's. The TeV spectrum for the 1998 season from the CAT Cherenkov imager corresponds to a quiet X-ray period (Piron et al. 2001b), with an average of 0.6 ASM counts per second, close to the year 2000 rate, and no significant flares.

Figure 10: CELESTE average flux measurement (gray bowtie) for Mrk 421 for 2001 February (MJD 51 956.98 to 51 963.11, Epoch (1) in Table 3 and text), superimposed on the spectral energy distribution taken from the multiwavelength study by Bazejowski et al. (2005). The envelope of the bowtie corresponds to the assumptions of $\alpha =1.8$ and 2.2 differential spectral indices. The inner set of lines above and below the bowtie are the statistical errors, while the outer set has in addition a 25% systematic uncertainty added in quadrature. The ASM X-ray activity level in February 2001 corresponds to the "medium'' activity state (the lower of the solid-curved high energy peaks) and thus to the VERITAS TeV measurements represented by the squares. The CAT spectrum in Piron et al. (2001a) for the year 2000 (thick black line) gives the same integral flux as for the nights in February 2001.

Open with DEXTER

6.3 1ES 1426+428

Table 5 and the light curve in Fig. 14 show that CELESTE detected no gamma-rays from 1ES 1426+428, indicating that the $2\sigma$ upper flux limit at 80 GeV is $0.2 \times 10^{-10}~{\rm erg~cm^{-2}~s^{-1}}$. The average ASM X-ray flux during CELESTE's observations was $0.2 \pm 0.1$ counts per second. For the two seasons leading to the detection of this blazar with the Whipple telescope the ASM rate ranged from $0.15 \pm 0.04$ to $0.65 \pm 0.06$, with an average near 0.3, with no discernible correlation between the weak X-ray and TeV signals (Horan et al. 2002).

Prior to the recent HESS result in Aharonian et al. (2006), 1ES 1426+428 was the most distant blazar used to constrain the density of extragalactic diffuse infrared light, possible because TeV photons can be absorbed by low energy photons via electron-positron pair production, $\gamma \gamma \rightarrow \rm e^+e^-$. In particular, Dwek & Krennrich (2005) took the spectra observed at high energies and, for a range of plausible assumptions about infrared spectral shapes and densities, removed absorption effects. Extrapolations of the de-absorbed spectra to 100 GeV suggested fluxes detectable by CELESTE. In the end, our non-detection during a low X-ray state does not constrain the diffuse infrared densities.

7 Conclusions

We have provided some of the first astrophysical flux measurements around 100 GeV, beyond the reach of EGRET on the Compton CGRO and below the range of previous ground-based telescopes, by recording atmospheric Cherenkov light with the reconverted solar tower facility at Thémis. We have increased our detector sensitivity as compared to our earlier work, to a level near that predicted in our original proposal. The gain came mainly from a better analysis procedure but also from refinements in our detector simulation and electronics improvements. 80% of the gamma-rays we detect have energies between 50 and 350 GeV. Bad weather at our site lead to a low duty-cycle, resulting in a modest number of detections in spite of our detector's good performance.

Our results for the blazar Mrk 421 help constrain the shape of the high energy peak of the spectral energy distribution. We have made the first detection of Mrk 501 in this range, although the detection is weak, and we have placed an upper limit on the flux from the distant blazar 1ES 1426+428 during a period of X-ray inactivity.

Table 4: Number of events before and after cuts for the Mrk 501 data sets. $N_\sigma$ is calculated after correcting for the $\sim$$80\%$ data acquisition efficiency.

	Figure 11: Left: $\xi$ distributions for Mrk 501 for the year 2000 data. Black dots show On-Off data. The broad histogram is Off data (i.e. cosmic ray background). The narrow histogram is simulated gamma-rays. The simulation and Off data sets are normalized to the number of On-Off events. Right: statistical significance of the excess as a function of the $\xi$ cut value.
Open with DEXTER

	Figure 12: Gamma- and X-ray light curves for Mrk 501. ${\rm MJD}$ 51600 and 52000 correspond to 2000 February 26 and 2001 April 4. Top: ASM 5-day running averages of the 2 to 12 keV X-ray flux, in cts/s (see ASM in references), and BeppoSax 2 to 12 keV fluxes in $10^{-10}~{\rm erg~cm^{-2}~s^{-1}}$, from Massaro et al. (2004). Bottom: CELESTE daily flux averages.
Open with DEXTER

Table 5: Number of events before and after cuts for the 1ES 1426+428 data sets. $N_\sigma$ is calculated after correcting for the $\sim$80% data acquisition efficiency.

Figure 13: CELESTE flux measurement for Mrk501 from the year 2000 (bowtie), bounded by the assumptions of $\alpha =1.8$ and 2.2 differential spectral indices. The inner set of lines, above and below the bowtie, correspond to the statistical uncertainty, while the outer set have in addition the 25% systematic uncertainty added in quadrature. The Beppo X-ray data is nearly simultaneous (Massaro et al. 2004). The EGRET data and the ASCA X-ray measurements are simultaneous with each other, and correspond to the SSC model curves (Kataoka et al. 1999). The TeV spectra are from the CAT Cherenkov imager (Piron et al. 2001b).

Open with DEXTER

	Figure 14: Light curves for 1ES 1426+428 - CELESTE detected no signal. The CELESTE daily averages run from 2002 March 22 to 2004 April 25. The ASM points are 5 day running averages.
Open with DEXTER

Acknowledgements

This work was supported by the IN2P3 of the French National Center for Scientific Research (CNRS) and by the Région Languedoc-Roussillon. We are grateful for the cooperation of Électricité de France. Antoine Perez, Jacques Maurand, and Stéphane Rivoire kept the Thémis facility operational through hard work and dedication. We thank Wei Cui, Jun Kataoka, and Dieter Horns for providing us with data files and macros.

References

• Aharonian, F. A., et al. 2002, A&A, 393, 89 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Aharonian, F. A., et al. 2003, A&A, 410, 813 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Aharonian, F. A., et al. 2004, ApJ, 614, 897 [NASA ADS] [CrossRef]
• Aharonian, F. A., et al. 2006, Nature, 440, 1018 [NASA ADS]
• ASM (All Sky Monitor) quick-look results are provided by the ASM/RXTE team, http://xte.mit.edu/ASM_lc.html (In the text)
• Baillon, P., et al. 1993, Astropart. Phys., 1, 341 [NASA ADS] [CrossRef] (In the text)
• Bernlöhr, K. 2000, Astropart. Phys., 12, 255 [NASA ADS] [CrossRef] (In the text)
• Blazejowski, M., et al. 2005, ApJ, 630, 130 [NASA ADS]
• Brion, E., et al. 2004, Proc. 35th COSPAR, Paris (In the text)
• Brion, E. 2005, doctoral thesis, http://doc.in2p3.fr/themis/CELESTE/PUB/papers.html (In the text)
• Bruel, P. 2004, Proc. Physics & Astrophysics in Space, Frascati
• Boone, et al. 2002, ApJ, 579, L5 [NASA ADS] [CrossRef] (In the text)
• Carson, J. 2005, doctoral thesis, http://www.astro.ucla.edu/~carson, ApJ in preparation (In the text)
• Djannati-Ataï, A., et al. 1997, A&A, 350, 17 [NASA ADS] (In the text)
• Djannati-Ataï, A., et al. 2002, A&A, 391, (In the text)
• Dwek, E., & Krennrich, F. 2005, ApJ, 618, 657 [NASA ADS] [CrossRef] (In the text)
• Gehrels, N., & Michelson, P. 1999, Astropart. Phys., 11, 277 [NASA ADS] [CrossRef] (In the text)
• Giebels, B., et al. 1998, Nucl. Instr. Meth. A, 412, 329 [CrossRef]
• Goret, P., et al. 1993, A&A, 270, 401 [NASA ADS] (In the text)
• Hartman, R. C., et al. 1999, ApJS, 123, 79 [NASA ADS] [CrossRef]
• Hayes, D. S., & Latham, D. W. 1975, ApJ, 197, 593 [NASA ADS] [CrossRef] (In the text)
• Heck, D., et al. 1998, Report FZKAS 6019, Forschungszentrum Karlsruhe (In the text)
• Holder, J., et al. 2001, 27th ICRC, Hamburg (In the text)
• Horan, D., et al. 2002, ApJ, 571, 753 [NASA ADS] [CrossRef] (In the text)
• Kataoka, J., et al. 1999, ApJ, 514, 138 [NASA ADS] [CrossRef] (In the text)
• Katarzynski, K., Sol, H., & Kus, A. 2001, A&A, 367, 809 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Kertzman, M. P., & Sembroski, G. H. 1994, Nucl. Instrum. Meth. A, 343, 629 [NASA ADS] [CrossRef] (In the text)
• Krawczynski, et al. 2001, ApJ, 559, 187 [NASA ADS] [CrossRef] (In the text)
• Krawczynski, et al. 2004, ApJ, 601, 151 [NASA ADS] [CrossRef] (In the text)
• Lavalle, J., Manseri, H., et al. 2006, A&A, 450 1
• Manseri, H. 2004, doctoral thesis, http://doc.in2p3.fr/themis/CELESTE/PUB/papers.html (In the text)
• Massaro, E., et al. 2004, A&A, 422, 103 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Mukherjee, R., et al. 1997, ApJ, 490, 116 [NASA ADS] [CrossRef] (In the text)
• de Naurois, M., Holder, J., et al. 2002, ApJ, 566, 343 [NASA ADS] [CrossRef] (In the text)
• Oser, S., et al. 2001, ApJ, 547, 949 [NASA ADS] [CrossRef] (In the text)
• Piron, F., et al. 2001a, A&A, 374, 895 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Piron, F. 2001b, in Proc. SF2A-Lyon [arXiv:astro-ph/0109286] (In the text)
• Piron, F., et al., 28th ICRC, Universal Academy Press, 5, 2607 OG 2.3 (Tsukuba 2003)
• Paré, E., et al. 2002, Nucl. Instr. Meth. A, 490, 71 [NASA ADS] [CrossRef]
• Smith, D. A. 1998, Proc. 19th Texas Symposium on Ultrarelativistic Astrophysics, Paris (In the text)
• Smith, D. A. 2005, Proc. "Atmospheric Cherenkov Detectors VII'', Palaiseau (In the text)
• Tavecchio, F., Maraschi, L., & Ghisellini, G. 1998, ApJ, 509, 608 [NASA ADS] [CrossRef] (In the text)