Systematic search for VHE gamma-ray emission from X-ray bright high-frequency BL Lac objects

All but three (M87, BL Lac and 3C 279) extragalactic sources detected so far at very high energy (VHE) gamma-rays belong to the class of high-frequency peaked BL Lac (HBL) objects. This suggested to us a systematic scan of candidate sources with the MAGIC telescope, based on the compilation of X-ray blazars by Donato et al. (2001). The observations took place from December 2004 to March 2006 and cover sources on the northern sky visible under small zenith distances zd<30 degrees at culmination. The sensitivity of the search was planned for detecting X-ray bright F(1 keV)>2 uJy) sources emitting at least the same energy flux at 200 GeV as at 1 keV. In order to avoid strong gamma-ray attenuation close to the energy threshold, the redshift of the sources was constrained to values z<0.3. Of the fourteen sources observed, 1ES 1218+304 and 1ES 2344+514 have been detected in addition to the known bright TeV blazars Mrk 421 and Mrk 501. A marginal excess of 3.5 sigma from the position of 1ES 1011+496 was observed and has been confirmed as a source of VHE gamma-rays by a second MAGIC observation triggered by a high optical state (Albert et al. 2007). For the remaining sources, we present here the 99% confidence level upper limits on the integral flux above ~200 GeV. We characterize the sample of HBLs (including all HBLs detected at VHE so far) by looking for correlations between their multi-frequency spectral indices determined from simultaneous optical, archival X-ray, and radio luminosities, finding that the VHE emitting HBLs do not seem to constitute a unique subclass. The absorption corrected gamma-ray luminosities at 200 GeV of the HBLs are generally not higher than their X-ray luminosities at 1 keV.


Introduction
Blazars belong to the most extreme objects in astronomy. Dominated by a non-thermal continuum spectrum, covering up to twenty decades in energy, they show variability on time scales of years down to minutes (Albert et al. 2007f;Aharonian et al. 2007a), and apparent luminosities exceeding 10 49 erg s −1 .
Morphologically, blazars show strongly collimated jets extending from scales not much larger than the event horizon of a supermassive black hole (Biretta et al. 2002) up to megaparsec scales. Superluminal motion of knots in the radio jets indicates relativistic bulk motion (Ghisellini et al. 1993). X-ray knots at distances of more than the radiative cooling length from the nucleus indicate in situ particle acceleration (Biretta et al. 1991), occuring at traveling and stationary shocks in the jet. According to the unified scheme (e.g. Urry & Padovani 1995), blazars are accreting supermassive black holes expelling a relativistic plasma jet at a small angle between the jet axis and the line of sight, with strong boosting of the observed emission due to relativistic bulk motion of the emitting plasma. BL Lacertae objects differ from the generally more luminous quasars by showing only faint or even absent emission lines, the absence of thermal big blue bump emission, and by not showing the otherwise typical luminosity evolution. The Spectral Energy Distribution (SED) of blazars shows two pronounced peaks, the first between IR and hard X-rays, which is commonly believed to be synchrotron radiation of highly relativistic electrons, and the second one at γrays. Depending on the location of the first peak, BL Lac objects are further divided into lowfrequency peaked BL Lacs (LBL, IR to optical) and high-frequency peaked BL Lacs (HBL, UV to X-rays) (Giommi & Padovani 1994). The second peak at high energies can be explained by inverse Compton scattering of low energy photons, produced as synchrotron radiation by the same population of electrons (Synchrotron Self Compton, SSC, Maraschi et al. (1992)), or from ambient thermal photon fields, which could enter directly into the emission region (Dermer & Schlickeiser 1993) or by scattering on material surrounding the jet (Sikora et al. 1994). The origin could also be due to hadronic processes associated with proton and ion acceleration which leads to electromagnetic cascades and proton synchrotron radiation (Mannheim 1993;Aharonian 2000;Muecke & Protheroe 2001). We cautiously remark, that the SED is probed at a sufficient level of sensitivity only in a limited range, there are still large gaps, in particular between 50 keV and 100 MeV where further peaks could show up.
In December 2004, when the regular observations with the Major Atmospheric Gamma-ray Imaging Cherenkov (MAGIC) telescope started, the number of known VHE blazars was six, all of them X-ray bright HBL objects. At the time of writing, the number has increased to 19 1 , includ-ing one LBL object (BL Lacertae, Albert et al. 2007d), one Flat Spectrum Radio Quasar (3C 279, Teshima et al. 2007) and the giant radio galaxy M87 (Aharonian et al. 2006b).
The detection of VHE γ-rays from cosmological distances is made difficult, due to absorption of γ-rays by photon-photon interactions with low energy photons from the evolving metagalactic radiation field (MRF). In the 100 GeV to 10 TeV range, far-infrared to optical photons are most important for the attenuation. It has been realized that this leads to a relation between the γ-ray cutoff energy and the source redshift known as the Fazio-Stecker relation (Fazio & Stecker 1979;Kneiske et al. 2004). The fact that the so-far detected VHE sources have much lower redshifts compared with the EGRET γ-ray sources is in line with the expected effect of γ-ray attenuation, although the lack of curvature in the observed spectra is a source of serious doubts (Aharonian et al. 2006a).
Due to the small field-of-view of an Imaging Air Cherenkov Telescope (IACT) and the limited duty cycle of ∼ 1000 hrs per year, promising canditates for VHE emission have to be selected carefully. All established TeV sources are bright X-ray sources, most of them with comparable luminosities in both regimes, which renders a systematic scan of the X-ray brightest HBL objects a reasonable approach.
Here we report on the results of such an approach pursued with the MAGIC telescope for a sample of 14 HBLs. In Sect. 2 the selection criteria for the sample are discussed, while the description of the observations can be found in Sect. 3. The data analysis technique is described in Sect. 4, and the analysis results are summarized in Sect. 5. A brief explanation of how the SEDs were obtained from the data (using archival radio and X-ray data as well as simultaneous optical data) can be found in Sect. 6. Finally, the resulting properties of the SED of HBLs and inferences on their luminosity function are discussed in Sect. 7.
The selection was made to avoid strong γ-ray attenuation at the energy threshold. At z = 0.3, the expected cut-off energy is still above 200 GeV, where MAGIC has its highest sensitivity. As the energy threshold increases with the zenith distance, all observations were carried out below 40 • , where the analysis energy threshold is around 200 GeV. As most of the established TeV sources show comparable luminosities in X-rays and in γrays, only the X-ray brightest HBLs were selected, leading to a cut at 2 µJy. Assuming the same luminosity at ∼ 200 GeV, it corresponds to ∼ 7% of the flux of the Crab Nebula, which would be detectable for MAGIC within 15 hrs.
The goal was to observe them for at least 15 hrs in order to establish new VHE sources and to put constraints on the SED of HBLs in a systematic fashion. The complete set is listed in Table 1.

Observations
The MAGIC telescope is a single dish IACT, located on the Canary island of La Palma (28.8 • N, 17.8 • W, 2200 m a.s.l.). A 17 m diameter tessellated parabolic mirror with a total surface of 234 m 2 , mounted on a light-weight space frame made from carbon fiber reinforced plastic tubes, focuses Cherenkov light from air showers, initiated by γ-rays or charged cosmic rays, onto a 576-pixel photomultiplier camera with a field-ofview of 3.5 • . The analogue signals are transported via optical fibers to the trigger electronics and each channel is read out by a 300 Msample/s FADC. Further details on the telescope can be found in Baixeras et al. (2004) and Cortina et al. (2005). Note that the readout system has been upgraded to a 2GSample/s FADC in February 2007 (Goebel et al. 2007 b Γ is the spectral index for the differential spectrum (dN/dE) at 1 keV, assuming a power law c Effective observation time after quality selection d Two measurements are above 3 µJy, one below 1 µJy e The earlier reported redshift of 0.200 was recently revised by Albert et al. (2007e) f Results published in Albert et al. (2007c) g Results published in Albert et al. (2006b) h Results published in Albert et al. (2007f) i The earlier reported redshift of 0.018 was recently revised by a lower limit (Sbarufatti et al. 2006) j Proposed but not observed due to bad weather k Results published in Albert et al. (2007b) telescope is pointing to the source (on-mode), the background has to be determined by so-called offdata, where the telescope points to a nearby sky region where no γ-ray source is expected. The off-data cover the same zd range with a similar night-sky background light intensity. The larger fraction of the source sample was observed in the so-called wobble mode, where the pointing of the telescope wobbles every 20 minutes between two symmetric sky locations with an angular distance of 0.4 • to the source. The background in the signal region can be estimated from sky locations placed at the same distance from the camera center as the candidate source. Except 1ES 0927+500 and 1ES 0414+009, all objects were monitored by the KVA telescope (http://users.utu.fi/kani/1m/index.html ) on La Palma in the optical R-band. None of the sources showed flaring activity in the optical during the MAGIC observations. The host galaxy corrected fluxes (Nilsson et al. 2007), taken simultaneously, averaged over the time of the MAGIC observations, are listed in Table 5.

Data analysis technique
The data are processed using the MAGIC Analysis and Reconstruction Software (MARS) (Bretz 2005a). A description of the different analysis steps can be found in Gaug et al. (2005) (including the calibration) and Bretz (2005b).
As the trigger rate strongly depends on the weather condition, only data with a rate above 160 Hz are used to ensure a high data quality.
The moments up to third order of the light distribution are used to characterise each event by a set of image parameters (Hillas 1985). For background suppression, a SIZE-dependent parabolic cut in WIDTH × LENGTH is applied (Riegel 2005). To reconstruct the origin of the shower in the camera plane, the DISP method is employed (Lessard et al. 2001) to estimate the distance between the centre of gravity of the shower and its origin. The third moment determines the direction of the shower development. The constant coefficient ξ from the parametrisation of DISP in the original approach is replaced in this analysis by ξ 0 + ξ 1 · (LEAKAGE) ξ2 , LEAK-AGE being the fraction of light contained in the outermost camera pixels. Thereby the trunca-tion of the shower images at the camera border is taken into account. These coefficients are determined separately for on-off and wobble data using simulated γ-showers, which are produced by CORSIKA, version 6.023 (Heck et al. 1998;Majumdar et al. 2005) for zd below 40 • and energies between 10 GeV and 30 TeV, following a power law with a spectral index -2.6.
The cut coefficients for the background suppression are optimised using Crab Nebula data, taken at similar zd. One set of cut coefficients is derived for data taken in on-off mode and one for wobble mode. The significance of a possible signal is determined from the distribution of the squared angular distance (θ 2 ) between the shower origin and the source position. The signal region is determined as θ < 0.23 • , corresponding to slightly more than two times the γ point-spread function of the MAGIC telescope.
For observation in on-off mode, the off-data have to be scaled to match the on-data. This is done in the region 0.37 • < θ < 0.80 • , where no bias from the source is expected. For wobble observations three regions, located symmetrically on a ring around the camera centre with the same distance from the centre as the source position, are defined as background regions. The scale factor is fixed to 1/3.
For every source the statistical significance according to equation 17 from Li & Ma (1983) is calculated. The upper limit on the excess rate is derived with a confidence level (c.l.) of 99%, using the method from Rolke et al. (2005), which takes also the scaling factor of the background into account. The upper limit for the excess rate is then compared to the excess rate of the Crab Nebula, which leads to an upper limit of the flux in units of the Crab Nebula flux above a certain energy threshold, assuming a Crab like spectrum. The energy threshold is here defined as the energy where the differential distribution (dN/dE vs. E) of simulated γ-ray events, surviving all cuts, peaks. Note that this threshold depends also on the spectral shape.
A large sample of data from the Crab Nebula in on-off as well as in wobble mode is used, spread over the entire observation campaign (see Table 2). This analysis shows, that the excess rate of the Crab Nebula is correlated to the background rate (after γ-hadron separation). Therefore, depend- ing on the background rate of the AGN, a reference value for the excess rate of Crab has to be calculated. This can be understood when taking into account that even after quality selection the rates fluctuate up to 20%, depending on weather conditions. In Fig. 1 and Fig. 2  f. = 13.5/3). As a constant fit gives even a worse result and the fit on the on-data shows a correlation between background-and excess rate, the linear fit for the wobble data is also used to calculate the reference values for the comparison of the excess rates. The Crab units are converted into a flux of photons cm −2 s −1 using the spectrum of the Crab Nebula from Albert et al. (2007g). The systematic error for the flux is estimated to be ∼ 30% (see Albert et al. (2007g) and discus- sion therein). For the u.l. determination there is also the uncertainty of the correct energy threshold (which depends on the source spectrum).

Results of the MAGIC observations
Within this observation program, VHE γ-rays were discovered from 1ES 1218+304 (Albert et al. 2006b) and 1ES 2344+515 was observed in a low flux state with high significance (Albert et al. 2007b). Mrk 421 was observed for more than 25 hrs in 2005. The results are discussed in detail in Albert et al. (2007c). Mrk 501 was observed from May to July 2005 with more than 30 hrs, revealing a high precision lightcurve on a day-by-day basis as well a two exceptionally short-time flares (see Albert et al. (2007f) for more details). For ten sources of the sample, no significant signal is seen. The 2006 observation of 1ES 1218+304 results in a weak signal of 4.6 σ (see Sect. 5.2). A slightly refined analysis of 1ES 1011+496 yields a hint of a signal with a 3.5 σ significance (see Sect. 5.3). The 6 results are listed in Table 3. Observations of 1ES 1727+502 are still pending.

Upper limits
The u.l.s are between 2.3% and 8.6% of the Crab Nebula flux. For a Crab-like spectrum the energy thresholds vary between (190±15) GeV and (230 ± 15) GeV, depending on the zd of the observation. For the threshold calculation the exact zd distribution of every observation is taken into account. As the Crab spectrum at ∼ 200 GeV is quite hard (spectral slope -2.26 for the differential energy spectrum), the u.l.s are also calculated for an -3.0 power law spectrum, which represents quite well the average slope of all HBLs detected at VHE so far. In Table 5 the energy threshold as well as the flux u.l. at 200 GeV are given under the assumption of a -3.0 power law.

1ES 1218+304
The source was observed in 2006 from January the 29th to March the 5th during 15 nights with in total 14.6 hrs. Figure 3 shows the distribution of the squared angular distance between the reconstructed shower origin of each event and the assumed source position. The vertical line indicates the signal region. The background is determined by three off-regions in the camera. The excess has a statistical significance of 4.6 σ (see also Table 3).

1ES 1011+496
The standard analysis performed for the whole source sample yields a significance of 2.5 σ, which gives already a hint for a possible signal. In a more refined analysis, the cut in θ 2 , which determines the signal region, was reduced to θ = 0.20 • . In case of a weak signal the increased signal to background ratio would lead to a higher significance. Also the SIZE-dependent cut for the background suppression is changed to a slightly higher value. The same analysis, performed on a data sample of the Crab Nebula, results in almost the same significance as with the standard coefficients, but with a ∼ 13% lower γ-rate due to the reduced γ acceptance and a 37% lower background rate. This analysis yields 3.5 σ for 1ES 1011+496 which -if interpreted as a detection -corresponds to an integral flux of F (> 180 GeV) = (1.26 ± 0.40) × 10 −11 photons cm −2 s −1 . Further observations with the MAGIC telescope, triggered by an optical outburst in March 2007, shows a clear signal of 6.2 σ within 18.7 hrs of observa- Table 3 Results of the analysis.

1ES 1426+428
The VERITAS collaboration reported a steep spectrum above 300 GeV for their observations in 2001, well fitted by a power law with spectral index −(3.50±0.35) (Petry et al. 2002). Extrapolating the spectral fit to 200 GeV, it yields an integral flux of 0.50 Crab above 200 GeV, which is by a factor of 10 larger than the u.l. presented in this work. Previous measurements yield a marginal detection in 2000 and upper limits for the data taken from 1995 to 1999 (Horan et al. 2002) with the most stringent one of 0.08 Crab above 350 GeV.
The HEGRA collaboration published a much harder spectrum at higher energies (above ∼ 800 GeV) for their combined 1999 and 2000 data, which were well fitted by a power law with spectral index −(2.6 ± 0.6) (Aharonian et al. 2002). An extrapolation of the power law yields an integral flux above 200 GeV of 0.075 Crab. Due to the large extrapolated energy range, combined with the large statistical error of 0.6 for the slope, the uncertainty is a factor of two. Further measurements in 2002 with the HEGRA telescopes showed the source in a 2.5 times lower flux state (Aharonian et al. 2003).
The u.l. for the flux above 200 GeV presented in this work, indicates a lower flux than measured from 1999 to 2001 during several campaigns with different telescopes, whereas it is consistent with the low flux level observed in 2002.

Spectral Energy Distribution
The complete set of sources, described in Sect. 2, amounts to 14 objects (without 1ES 1727+502) and includes the five established TeV sources Mrk 501, 1ES 1218+304, 1ES 1426+428, Mrk 421 and 1ES 2344+514. To better understand the spectral energy distribution of VHE gamma ray emitting blazars, the sample is enlarged by including ten more sources which fullfill the same selection criteria in X-ray flux and redshift, but which deviate only in the zenith angle cut. This means, they can only be observed under large zenith distances from the MAGIC site or not at all, in which case there is still information available from corresponding H.E.S.S. observa-tions. An exception is PG 1553+113 where the redshift is not known. For the total, enlarged sample of 24 HBLs, multi-wavelength data are collected in the following bands: radio (5 GHz, Donato et al. (2001)), optical (R-band, 640 nm, simultaneous data from KVA or, if not available from Donato et al. (2001)), X-rays (1 keV, Donato et al. (2001)), and γ-rays (200 GeV, see Table 4, Table 5 and ref. in Table 1). For 1ES 1426+428, in addition to the upper limit derived in this work, the extrapolation of the spectrum, as measured by HEGRA in 1999/2000, is used to mark the detected flux. The optical data were corrected for galactic extinction, using the coefficients from the NASA Extragalactic Database (NED), which are calculated according to Schlegel et al. (1998).
For the γ-ray flux, sizeable attenuation is expected from current models of the MRF (Hauser & Dwek 2001;Kneiske et al. 2004). Therefore all u.l. at 200 GeV as well as the measured fluxes of the detected HBLs are corrected for the absorption by multiplying with exp(τ (200 GeV, z)) where τ denotes the pair production optical depth. The "best-fit 2006" MRF-model (Kneiske 2007, in preparation) is employed. This model is based on the "best fit" model from Kneiske et al. (2004), but with a lower star formation rate to keep the energy density in the optical band closer to the lower limits, derived from the galaxy number counts. It is also consistent with the u.l. derived by Aharonian et al. (2006a) from the VHE spectrum of 1ES 1101-233 (z = 0.186). The optical depth values for 200 GeV photons for all sources of this sample are listed in Table 5.
All fluxes are K-corrected. Radio spectral indices 2 of ten of the sources can be found in Landt (2003). For the other 14 objects the average value α R = 0.23 of the ten sources is used. For the optical data, the spectral indices of nine sources, calculated at slightly higher wavelengths, are taken from Bersanelli et al. (1992). For the other 15 objects, the average value α O = 0.65 of the nine sources is used. At 1 keV, the spectral indices are taken from Donato et al. (2001), except for 1ES 0229+200, which is not included in this compilation. Instead, the flux is taken from Costamante & Ghisellini (2002) together with the average value for the spectral index α X = 1.36 of all other sources. At 200 GeV the measured spectral indices are used for the detected sources, while for the non-detected ones the average value α γ = 2.0 is used. To take into account the energy dependent attenuation at VHE, which causes a hardening of the spectrum, the measured spectral indices are changed by -0.4 for 0.1 < z < 0.2, -0.8 for 0.2 < z < 0.3 and remain unchanged for z < 0.1.
A special treatment is necessary to derive the flux at 200 GeV from 1ES 0229+200, recently discovered at VHE γ-rays (Aharonian et al. 2007c). The spectrum is measured above 580 GeV, well fitted by a power law with spectral index −(2.51 ± 0.19). As the source is located at z = 0.1396, strong absorption is expected at these energies. Therefore the spectrum is first deabsorbed and afterwards extrapolated to lower energies (Table 4). The resulting intrinsic spectrum is well described by a power law with a spectral index of −(1.09 ± 0.25) (flux normalisation: (4.24 ± 0.81) × 10 −12 cm −2 s −1 TeV −1 at 1 TeV). This result is in good agreement with the results from Stecker & Scully (2007), yielding model dependent intrinsic spectral indices in the range from 1.1 ± 0.3 to 1.5 ± 0.3.
After these corrections, the broad-band spectral indices 3 α 1−2 between the different energy regimes are calculated. Also the luminosities νL ν are calculated assuming isotropic emission and with the use of the following cosmological parameters: H 0 = 71 km s −1 Mpc −1 , Ω Λ = 0.73, and Ω m = 0.27.

Gamma ray emitting HBLs?
One may ask whether the gamma ray emitting HBLs can be distinguished from other HBLs based on their spectral energy distributions. Finding the answer is hampered by a number of problems. The variable peak frequencies, which cannot be easily detected in fixed energy bands, the gamma ray attenuation due to pair production in the metagalactic radiation field, and the flux variability. From Donato et al. (2001), the amplitude of the 3 α 1−2 = − log(F 1 /F 2 )/ log(ν 1 /ν 2 ) , ν 1 < ν 2 flux variability at 1 keV amounts to a factor of six for the sources with multiple entries in the catalogue. Similar or larger amplitudes can be expected at VHE. However, our sample is not triggered by flux variability, the duty cycle of flares generally seems to be rather low, and the observed fluxes or flux upper limits may therefore be characteristic of the quiescent average fluxes. Figure 4 shows the broad-band spectral index α RO vs. α OX for all 24 HBLs as described in the previous section. The distribution is quite homogeneous. As the data are not simultaneously taken, the uncertainties due to flux variations have to be taken into account. In the case of α OX a flux variability of a factor six at 1 keV coresponds to a change in the spectral index of 0.29 (if the optical flux remains the same). This is still below the difference of 0.6 between the lowest and highest values of α OX for the detected VHE sources. Thus there are significant differences in the SED among the HBLs studied here. As the variability in the radio and optical band for HBLs is lower than at X-rays or VHE γ-rays and the change of the spectral index α RO is smaller for different flux ratios, the variation of α RO is much lower than for α OX . However, it is not possible to distinguish between sources detected at VHE and non-detected ones based on spectral indices.

Gamma-to-X-ray spectral index
To unveil the physical state of the emitting plasma, we seek to find the characteristic ratio of the two peaks in the spectral energy distribution of HBLs which is related to the ratio of photon and magnetic field energy densities. Figure 5 shows the broad-band spectral index α Oγ vs. α Xγ . Both indices are distributed in a narrow band around unity for all detected sources, where α = 1 represents the case that the energy output in both frequency bands is the same. The constraints from the u.l.s on the γ-ray flux can not exclude the region which is spanned by the detected TeV sources. In the framework of SSC models, the optical and X-ray band belong to synchrotron emission of highly relativistic electrons (first peak in the SED), while the VHE γ-rays are produced by inverse Compton scattering (second peak in the SED). If the magnetic energy density u B is equal to the photon energy density u ph , the energy output for the synchrotron and the Inverse Compton Table 4 HBLs detected at VHE which do not belong to the sample described in Sect. 2 . emission as well as the peak luminosities are the same. In case of HBLs, the peak frequency of the synchrotron emission is always at higher frequencies than the optical band and below 1 keV for most of them (except the extreme blazars with a hard spectrum at 1 keV). At VHE energies, most of the detected sources show a soft spectrum above ∼ 200 GeV, indicating observed peak energies below 100 GeV. Due to absorption of VHE γ-rays in the MRF, the intrinsic peak energies could also reach several TeV for sources located at higher redshifts. In that context, the scattering of α Oγ and α Xγ around unity could be explained by a continuous distribution of peak frequencies for HBLs, measured with fixed bandwidth. Figure 6 shows the broad-band spectral index α Xγ vs. X-ray luminosity ν X L X . The average energy output at 1 keV never exceeds significantly the one at 200 GeV (α Xγ > 0.97). For half of the detected sources, the energy output in these bands is almost the same (α Xγ ∼ = 1), while for the other ones the energy output at 200 GeV is significantly lower. There is a tendency that this effect shows up at higher X-ray luminosities. For four of the ten u.l., α Xγ = 1 can not be excluded, whereas for the other six sources α Xγ has to be larger than unity. For the non-detected sources of the sample, further observations with longer exposures are needed Table 5 Upper limits on the γ-ray flux at 200 GeV under the assumption of a -3.0 power law spectrum togehter with the optical depth and the simultaneous optical data.  to reach the α Xγ = 1.12-line (corresponding to a nine times lower output at 200 GeV compared to 1 keV), which includes all HBLs detected at VHE so far. Note that if α Xγ = 1 for the peak frequencies is valid, the tendency of increasing α Xγ values with increasing luminosities at 1 keV could be interpreted as a shift of the inverse Compton peak to lower values. Similar studies aiming at a characterization of VHE blazars are also pursued by Wagner (2007).

Source
7.3. Constraints on the luminosity function of γ-ray emitting HBLs Figure 7 shows the luminosity ν γ L γ at 200 GeV vs. the redshift. All detected sources are above or within the line that marks the corresponding luminosity to a flux of 2 µJy at 1 keV. The absorption of γ-rays by the MRF increases with redshift, so that at a redshift of z = 0.3, the emitted luminosity becomes twice larger than the measured one.
The luminosity function at γ-ray energies of HBLs is poorly known, since there has not been a complete survey and the number of detected sources is still rather low. Nevertheless, we can try to constrain the VHE luminosity function from our observations. To this end we derive upper limits on the cumulative omnidirectional flux at 200 GeV from X-ray bright HBLs below z = 0.3, noting that GLAST will measure the diffuse extragalactic background up to 200 GeV. The selection criterion for the declination of the sample corresponds the a patch of the sky with a size of 5.55 sr (or 44% of the sky). The sum over all 14 sources of the sample divided by the 5.55 sr patch results in an upper limit on the total intensity at 200 GeV of I VHE (200 GeV) = ε·2.76×10 −9 GeV cm −2 sr −1 s −1 , where ε accounts for the incompleteness of the sample. Assuming an isotropic distribution of HBLs, ε should be larger than unity. We know about two sources, not included in this calculation (1ES 1727+502 and 1ES 0229+200). There are also 5 HBLs from Donato et al. (2001) that fullfill the criteria of declination and X-ray flux, but are located at higher redshifts. Including these sources the sample would increase to 21 sources. With respect to the ∼ 50% sky coverage of the Einstein Slew Survey, where most of the sources belong to, ε = 3 seems to be a reasonable assumption (factor 1.5 to take the "known" sources into account and a factor 2 for the assumed sky coverage of 50% for the ex- Not that this u.l. is already conservative as also sources are included which are not in the Einstein Slew Survey (would result in a higher sky coverage) and the u.l. is dominated by the two brightest sources Mrk 501 and Mrk 421 (which would mean that three times more sources would not necessarily lead to a three times higher flux). Recently Kneiske & Mannheim (2007) showed that HBLs could contribute up to 30% of the extragalactic background radiation at GeV energies, when including cascade emission from sources at higher redshifts. The luminosity function used for their calculation was derived from the X-ray luminosity function (Beckmann et al. 2003) assuming the same luminosity above 300 GeV as from 0.5 keV to 2 keV. For the HBL contribution of faint point sources at 200 GeV, they obtain (somewhat model-dependent) I

(point) KM
(200 GeV) = (4 − 6) × 10 −8 GeV cm −2 sr −1 s −1 which is within the upper limit obtained here. For the total intensity, including the diffuse component due to electromagnetic cascading, their result is I (diffuse) KM (200 GeV) = 1.0 × 10 −7 GeV cm −2 sr −1 s −1 . This estimate is based on assuming that the emitted VHE spectra generally have peaks at energies well in excess of 200 GeV.

Conclusions
Searching for VHE γ-ray emission from a sample of twelve X-ray bright HBL objects, 1ES 1218+30.4 at a redshift of z = 0.182 was discovered for the first time at VHE (Albert et al. 2006b). The already established VHE source 1ES 2344+514 has been detected with high significance (Albert et al. 2007b), albeit in a state of low activity. For ten sources no significant signal was seen, resulting in upper limits on their integral flux above ∼ 200 GeV between 2.3% and 8.6% of the Crab Nebula flux on a 99% confidence level. There is a hint for a signal from 1ES 1011+496 on a 3.5 σ level, which has been confirmed as a source of VHE γ-rays by a second MAGIC observation campaign, triggered by a high optical state (Albert et al. 2007e).
With fixed-schedule observations a bias to flaring emission states was avoided, tacitly assuming that the duty cycle of the flares is short when compared with the exposure time. As shown in the case of 1ES 2344+514, quiescence does not necessarily preclude the sources from being detected. As a matter of fact, the upper limits obtained for the other sources still lie in a region of parameter space bracketed by the detected sources, and correspond to a VHE energy flux on the level of the X-ray energy flux with an indication of an increasing α Xγ with increasing X-ray luminosity. It thus seems to be a question of time that more sensitive telescopes, and in particular those with a lowered energy threshold such as MAGIC-II, potentially capturing shifting peaks, will eventually lead to a detection of all known bright HBLs. It will be important to obtain a complete catalogue of HBLs from the planned eROSITA and GLAST all-sky surveys to study them on a much larger statistical basis in the future.
The detected sources deviate in no apparent pattern from the so-far non-detected sources. A spectral shape with two equal-height bumps is ex-pected in Synchrotron Self Compton models for the case of balanced energy densities of photons and the magnetic field. The dispersion in the distribution of the X-to-γ-ray luminosity ratio (the lowest X-to-γ-ray luminosity ratio of a detected HBL is 1/9) would then reflect variations of the peak position w.r.t. the observed band or the effect of flux variability (e.g., high state X-ray emission vs. quiescent VHE emission).