text:
# 5@5 - a 5 GeV energy threshold array of imaging atmospheric Cherenkov telescopes at 5 km altitude

## 1 Introduction

The high detection rate, the ability of effective separation of electromagnetic and hadronic showers, and the good accuracy of reconstruction of the direction of primary $`\gamma `$-rays are three remarkable features of the imaging atmospheric Cherenkov telescope (IACT) technique (see e.g. ). The recent detections of TeV $`\gamma `$-rays from several galactic and extragalactic objects (see e.g. ) provide the basis for the further development of ground-based gamma-ray astronomy.

Qualitative improvements of the IACT technique in the next few years will most probably be linked to stereoscopic observations of $`100\mathrm{GeV}`$ $`\gamma `$-rays . This approach not only allows an unambiguous determination of the energy and of the arrival direction of $`\gamma `$-ray primaries on an event-by-event basis, but also significantly improves the efficiency of rejection of hadronic showers produced by cosmic rays , as it was recently demonstrated by the HEGRA system of 5 imaging telescopes operating in the energy region from 500 GeV to 20 TeV .

One of the important issues for future detectors is the choice of the energy region based on two principal arguments: (a) astrophysical significance (goals) and (b) the experimental feasibility/reliability (cost). If one limits the energy region from above to a relatively modest energy threshold around 100 GeV, then the performance of IACT arrays and their practical implementation can be predicted with confidence. In practice, an energy threshold of $`100`$ GeV can be achieved by a stereoscopic system of IACTs consisting of 10 m diameter class optical mirrors and equipped with high resolution cameras, based on conventional photo-multipliers (PMTs). Currently three such arrays are under development/construction in Australia (CANGAROO-3 ), in Namibia (H.E.S.S. ), and in Arizona (VERITAS ). With their superior angular resolution of several arcminutes and energy flux sensitivity close to $`10^{13}\mathrm{erg}/\mathrm{cm}^2\mathrm{s}`$, these projects perfectly suit to the energy range from 50 GeV to 10 TeV, which, from the point of view of scientific motivations and the potential astronomical targets, can be considered as a spectral domain in its own right.

On the other hand, it is expected that the next generation major satellite $`\gamma `$-ray mission GLAST (see e.g. ), the successor of the EGRET instrument aboard the Compton Gamma Ray Observatory, will extend the exploration of the $`\gamma `$-ray sky up to 100 GeV. Thus the gap between space-based and ground-based $`\gamma `$-ray instruments will eventually disappear. It should be noticed, however, that in many cases this statement has a rather conditional or even symbolic character. Although at GeV energies GLAST will improve the EGRET sensitivity by almost two orders of magnitude, the capability of GLAST (and likely that of any post-GLAST space-based project) at energies well beyond 10 GeV will be quite limited because of the limited detection area.

This circumstance justifies recent activities to reduce the energy threshold of atmospheric Cherenkov detectors below 100 GeV. Currently two low threshold projects, CELESTE in France  and STACEE in USA , based on the concept of conversion of existing solar power plants into atmospheric Cherenkov telescopes, are in their final stage of realization. Although this technique allows in principle a reduction of the energy threshold down to 20-30 GeV (see e.g. Ref. ), the challenge remains to reach an adequate detection efficiency of $`\gamma `$-rays at such low energies. Another approach has been proposed by the MAGIC collaboration with a single imaging Cherenkov telescope having a large, 17 m diameter reflector. Equipped with a standard PMT-based high resolution camera this telescope is expected to allow effective detection of $`\gamma `$-rays at energies above 30 GeV .

With some exceptions, the GeV $`\gamma `$-ray sources ($`E0.1\mathrm{GeV}`$) are expected to be quite different from TeV sources ($`E0.1\mathrm{TeV}`$). The proximity of the intermediate domain below 100 GeV to the energy range covered by EGRET suggests that many objects established as GeV emitters have a good chance to be detected also by the above mentioned ground-based instruments. This argument, however, cannot yet guarantee definite success. Indeed, although the two largest $`\gamma `$-ray source populations identified by EGRET, radiopulsars and distant AGN, do not show a significant steepening or cutoff up to 10 GeV, the theoretical studies of $`\gamma `$-ray production and absorption conditions in these objects, as well as rather general phenomenological considerations predict cutoffs in the energy spectra around 10 GeV or less. In addition, for any reasonable model of the diffuse extragalactic cosmic background radiation, we should expect sharp cutoffs in the spectra of distant extragalactic objects with redshift $`z1`$ at energies as low as 50 GeV (see e.g. ). This implies that for the study of cosmologically distant sources, like the GeV blazars discovered by EGRET, or Gamma Ray Bursts (GRBs), the energy threshold of the detectors should be less than 10 GeV at which energy the Universe is most likely transparent up to at least $`z3`$. An instrument like GLAST, operating effectively in the 0.1 to 10 GeV energy region, nicely suits this task. In particular, it is expected that the number of AGN that GLAST will detect could exceed several thousands . At the same time, the relatively small detection area of GLAST, $`A_{\mathrm{eff}}0.8\mathrm{m}^2`$, limits the potential of this instrument for detailed studies of the temporal and spectral characteristics of highly variable $`\gamma `$-ray sources like blazars, which have variability timescales less than a few hours, or of solitary events like GRBs with a duration of $`10^2`$ to $`10^3`$ seconds. In this regard, GLAST can hardly match the performance of current X-ray detectors that have similar detection areas but operate in a regime of photon fluxes that exceed the fluxes of MeV/GeV $`\gamma `$-rays by many orders of magnitude.

The idea of a “Gamma-ray timing explorer” to study transient $`\gamma `$-ray phenomena with an adequate photon detection rate motivated, to a large extent, the present investigation. It concerns the possible extension of the domain of ground-based Cherenkov technique with its huge detection area of $`10^4`$ to $`10^5\mathrm{m}^2`$ down to energies of several GeV. We shall argue that such a goal could be best achieved by stereoscopic systems of several large, 20 m class imaging atmospheric Cherenkov telescopes located at very high, $`H5\mathrm{km}`$, mountain altitude.

## 2 Concept of an IACT array with 5 GeV threshold

The concept of stereo imaging is based on the simultaneous detection of a single air shower in different projections by at least two telescopes, separated at a distance comparable with the “effective radius” $`R_\mathrm{C}100\mathrm{m}`$ of the Cherenkov light pool. The stereoscopic approach allows (i) unambiguous and precise reconstruction of shower parameters on an event-by-event basis, (ii) superior rejection of hadronic showers, and (iii) effective suppression of the background light from different sources - the nigh sky background (N.S.B.), local muons, etc. . All these three advantages over single IACTs have been convincingly demonstrated at TeV energies by the HEGRA IACT system .

Compared with single (“stand alone”) telescopes, which can adequately measure the shower inclination in the direction perpendicular to the plane containing the telescope axis, but poorly in the in-plane direction, the stereoscopic approach allows full reconstruction of the arrival direction of individual $`\gamma `$-ray showers. Apart from the good directional information, stereoscopic IACT systems make use of the fact that the Cherenkov images of a shower detected by several different, spatially separated telescopes, are only partially correlated. Therefore the stereoscopic measurements significantly improve the efficiency of rejection of hadronic (background) showers at both the hardware (trigger) and the software levels. The only disadvantage of the stereoscopic approach is a non-negligible loss in the detection rate because of the overlap of the shower detection areas of individual telescopes located from each other at distances $`2R_\mathrm{C}`$. However, this loss of statistics is compensated, especially for steep spectra of primary $`\gamma `$-rays, by a significant reduction of the energy threshold of the telescopes operating in coincidence mode.

At $`\gamma `$-ray energies above 100 GeV stereoscopic IACT arrays do provide an excellent angular resolution of about $`0.1^{}`$ or less, and a “gamma/hadron” separation efficiency (including the hadron rejection at the trigger level) of 1000:1. This improves the flux sensitivity dramatically compared with the sensitivity of single telescopes. The efficiency of the imaging technique is somewhat lower at energies below 100 GeV. In particular, the Cherenkov images become less elongated and less regular. In practice this introduces significant uncertainties in the reconstruction of image parameters. Even so, below it will be shown that the performance of the stereoscopic imaging remains adequately high even below 10 GeV.

The effective energy threshold of IACTs is basically determined by two conditions: (i) the number of photoelectrons in the image should be sufficient for an appropriate image analysis; typically, “good imaging” requires $`n_{\mathrm{ph}.\mathrm{e}.}^{\mathrm{min}}50100`$ electrons111For example, in the case of the HEGRA stereoscopic IACT system the minimum number of photoelectrons corresponding to the showers classified as “high quality events”, i.e. the showers which are accepted for the further image analysis, is close to 40 electrons per telescope ., (ii) the accidental trigger rate introduced by the N.S.B. should not exceed the detection rate of $`\gamma `$-rays. For an ideal IACT the first condition in principle should dominate over the second, technical condition.

The number of photoelectrons detected by the imager - a multi-pixel camera placed in the focal plane of the mirror - depends on (i) the telescope photo-electron response (or aperture), $`S_{\mathrm{ph}.\mathrm{e}.}=S_{\mathrm{mir}}\xi _{\mathrm{ph}\mathrm{e}}`$ (where $`S_{\mathrm{mir}}`$ is the geometrical area of the mirror, and $`\xi _{\mathrm{ph}\mathrm{e}}`$ is the photon-to-photoelectron conversion factor), and (ii) the density of optical Cherenkov photons $`\rho (R,E)`$ produced at the typical distance $`R100\mathrm{m}`$ from the axis of the shower initiated by a primary $`\gamma `$-ray of energy $`E`$.

The density of the Cherenkov light at an elevation of $`H=5`$ km a.s.l. within 100 m of the shower core produced by a primary $`\gamma `$-ray photon of energy $`E=5\mathrm{GeV}`$ is close to $`1\mathrm{photon}/\mathrm{m}^2`$. Correspondingly, the number of photoelectrons detected by a telescope with aperture $`S_{\mathrm{ph}.\mathrm{e}.}`$ is equal to $`n_{\mathrm{min}}1S_{\mathrm{ph}.\mathrm{e}.}`$. Thus, for reduction of the energy threshold down to 5 GeV, the telescope aperture $`S_{\mathrm{ph}.\mathrm{e}.}`$ should be as large as $`50\mathrm{m}^2`$, assuming that the minimum number of electrons required for an image analysis $`n_{\mathrm{ph}.\mathrm{e}.}^{\mathrm{min}}50`$. Detailed Monte-Carlo calculations presented below confirm that the energy threshold of such an instrument, determined as the energy at which the differential $`\gamma `$-ray detection rate reaches to its maximum, could indeed be as low as 5 GeV. For conventional aluminized optical mirrors and PMT-based cameras with typical conversion factor $`\xi _{\mathrm{ph}\mathrm{e}}0.150.2`$, $`S_{\mathrm{ph}.\mathrm{e}}=50\mathrm{m}^2`$ would require a large optical reflector of approximately 20 m diameter. With the successful development of novel, fast (nsec) detectors of optical radiation with a quantum efficiency exceeding $`50\%`$, the energy threshold of the telescopes could be pushed further down to 2 or 3 GeV, which is an absolute limit determined by the minimum energy of secondary electrons capable of producing Cherenkov light in the upper atmosphere.

The Cherenkov light density increases monotonically with elevation. Thus, for the given telescope configuration, the installation of an IACT array at very high altitude would allow a straightforward and unconditional reduction of the energy threshold. At the same time the increase of the telescope aperture leads not only to a proportional increase of number of photoelectrons registered per shower, but also to an increase of the accidental rate caused by the N.S.B. Thus, a significant gain from large aperture telescopes, with the ultimate goal to operate the telescopes in the linear regime, $`E_{\mathrm{th}}1/S_{\mathrm{ph}.\mathrm{e}}`$, can be achieved only through effective suppression of the N.S.B. Otherwise the reduction of the energy threshold would be rather slow, $`E_{\mathrm{th}}1/S_{\mathrm{ph}.\mathrm{e}}^{1/2}`$, and therefore would be difficult to justify economically.

In the imaging technique, the accidental events introduced by the N.S.B are suppressed by a trigger condition that requires signals above a threshold $`q_0`$ (photoelectrons) in $`m2`$ adjacent pixels. A smaller pixel size not only reduces significantly the noise level due to the N.S.B., but it allows higher trigger multiplicity as well. At the same time the optimal size of the pixel is determined by the condition of “good imaging”, which implies that the area which is covered by the minimum number of pixels (typically 10) used in the image analysis, should not exceed the characteristic size of the image. At $`E10\mathrm{GeV}`$, the image area of electromagnetic showers is less than $`0.2`$ square degree. Therefore at such low energies high resolution cameras with a pixel size close to $`0.1^{}`$ provide an adequate imaging quality. Although such small pixels allow also significant reduction of the N.S.B., it appears that for a single IACT with aperture $`50\mathrm{m}^2`$, even a pixel size of $`0.1^{}`$ is still not sufficient to operate the telescope at the minimum possible energy threshold, i.e. in the regime when the detection threshold is determined by the Cherenkov light amplitudes rather than by the N.S.B. noise. Further reduction of the pixel size makes the design of the imager with many thousands of channels technically very difficult. A solution of this dilemma is the stereoscopic mode of observations which offers a more feasible and economic approach. It requires simultaneous detection of a shower by at least two telescopes. Because of the flat lateral distribution of the Cherenkov radiation from electromagnetic showers, this requirement does not effect the $`\gamma `$-ray detection efficiency (if the distance between telescopes does not significantly exceed 100 m), but significantly reduces the accidental rate caused by the N.S.B. Below we will indeed show that the suppression of the N.S.B. by a stereoscopic system of IACTs with camera pixel size $`0.1^{}`$ is sufficient to operate the $`50\mathrm{m}^2`$ aperture telescopes in the energy regime below 10 GeV.

The arrangement of an IACT array, in particular the number of telescopes, and the spacing between them, can be understood from the following simple considerations. The Cherenkov light pool on the ground produced by primary $`\gamma `$-rays of energies $`E100\mathrm{GeV}`$ has a flat radial distribution with a pool radius of approximately 100 m. Hence the optimal spacing between the telescopes should be of the order of 100 m. A significantly smaller spacing reduces not only the detection area, but also the quality of images. A spacing of significantly more than 100 m reduces the coincidence rate dramatically, especially at low energies, and thus increases the energy threshold. The special investigation carried out for 100 GeV class telescopes , generally confirms this simple conclusion. We expect that this should be the case also for the sub-10 GeV array, although detailed Monte-Carlo simulations are needed for optimization of the arrangement of reflectors in this array. Here we consider a reasonable arrangement with a spacing of 100 m. The question of optimization of the spacing will be discussed in a separate paper. The minimum number of telescopes is determined from the consideration that at least 3 stereoscopic views are needed to reconstruct the shower parameters reliably . For clarity, we assume here a square baseline with telescopes at the four corners, placed at an altitude $`H`$ of 5 km a.s.l. In order to have more homogeneous coverage of distances of impact parameters, especially for showers with cores outside the square, we assume in addition a further (5th) telescope in the center of the square. The main objective of this paper is to study the basic performance of a 5 GeV IACT array placed at very high altitude. The optimization of the arrangement of the array is outside of the scope of this paper, but we believe that the suggested layout is not far from the optimum design.

## 3 Basic characteristics of Cherenkov radiation induced by sub-10 GeV gamma-rays

The results discussed in this section have been obtained by using the ALTAI code  that simulates electromagnetic and hadronic showers and their Cherenkov radiation in the Earth’s atmosphere. This code has been used before for calculations of the performance of the HEGRA system of IACTs , as well as for studies of the expected characteristics of the new generation “100 GeV - threshold” IACT arrays  that are the basis of the H.E.S.S. (High Energy Stereoscopic System) project in Namibia . The predictions for the HEGRA IACT system have been thoroughly checked using the detected hadronic (background) showers produced by cosmic rays with a relatively well known energy spectrum and mass composition . More importantly, the calculations of the characteristics of the instrument for primary $`\gamma `$-rays, in particular the lateral distribution of the Cherenkov light, which determines the $`\gamma `$-ray detection area, the angular resolution, and the gamma/hadron separation efficiency, have been experimentally confirmed by the HEGRA collaboration using $`\gamma `$-ray data obtained during the active state of Mkn 501 in 1997 . This extraordinary high state of the source, with a duration of several months, resulted in approximately 40,000 $`\gamma `$-rays in the energy interval from 500 GeV to 20 TeV, detected by the HEGRA IACT system under almost background-free conditions .

In this paper we are interested in the possibility of detecting primary $`\gamma `$-rays over an interval of very low energies from several GeV to 100 GeV, by a system of IACTs located at $`H5\mathrm{km}`$ a.s.l. Both conditions are rather extreme, and exhibit features which differ significantly from traditional studies of air showers. In particular, at energies of primary $`\gamma `$-rays below 10 GeV, we deal with Cherenkov radiation from only a handful of first generation electrons, while at TeV energies the Cherenkov radiation is contributed by a large number of electrons produced during the full cascade development. At such low energies we have therefore to expect more fluctuations of the parameters characterizing the showers and their Cherenkov radiation. In particular, the Cherenkov images of sub-10 GeV showers are expected to have a less regular shape, as well as to be more strongly affected by the geomagnetic field, as compared with the showers in the $`100\mathrm{GeV}`$ energy region. Also, both the very low energy domain and the very high elevation of the location of the suggested IACT array imply not only better transparency conditions for the Cherenkov radiation, but also a non-negligible reduction of the light produced above the telescopes from hadronic showers. All these effects add new features to the characteristics of the Cherenkov radiation of air showers.

### 3.1 Lateral and longitudinal development of showers

The lateral distribution of Cherenkov radiation from air showers, induced by $`\gamma `$-rays of energy 10 GeV and by protons of energy 100 GeV, are shown in Fig. 1 at 3 different observational levels - $`H=`$2.2 km, 3.5 km and 5 km. It is seen that within the radius of 100 m from the shower core which typically determines the detection area of IACTs, the density of Cherenkov radiation from $`\gamma `$-ray induced showers has (i) a rather flat lateral distribution, and (ii) increases significantly with elevation above sea level. In particular, the rise from 2.2 km to 5 km results in an increase of the density of the Cherenkov light by a factor of 2 to 3 (Fig. 1a). This elevation effect is less pronounced for proton-induced showers ($`50`$ per cent at $`r=100\mathrm{m}`$; see Fig. 1b). This implies that the choice of very high elevation for the operation of telescopes would allow a significant reduction of the energy threshold, by a factor of $`2.5`$, as well as a noticeable improvement of the background rejection of hadronic showers.

The slow increase of the Cherenkov light density from hadronic showers with the rise of elevation is explained by the effect of deep penetration of the “$`\pi ^\pm `$-branches” (sub-showers) of the cascade into the atmosphere (Fig. 2), which results in the production of a non-negligible amount of Cherenkov light below the observation level at $`H5\mathrm{km}`$. This effect is quantitatively demonstrated in Fig. 3 where the longitudinal distributions of the Cherenkov light, produced by a 10 GeV $`\gamma `$-ray photon and a 100 GeV proton, are shown.

The low energy threshold itself already leads, in fact more strongly than the elevation effect, to a significant improvement of the gamma/hadron ratio. This effect, caused by the reduction of the yield of Cherenkov light by sub-100 GeV hadronic showers (see e.g. ) is demonstrated in Fig. 4 for $`H=`$5 km a.s.l. The density of Cherenkov light in the electromagnetic showers is approximately proportional to the $`\gamma `$-ray energy down to $`E10\mathrm{GeV}`$, the dependence becoming a bit stronger below 10 GeV (Fig. 4a). At the same time, the density of

Cherenkov light from hadronic showers drops significantly faster with reduction of the primary energy (Fig. 4b). In Fig. 5 we show the ratio of densities of the Cherenkov light $`\rho _{100}^{(\gamma )}(E_\gamma )`$ and $`\rho _{100}^{(\mathrm{p})}(E_\mathrm{p})`$ produced by $`\gamma `$-rays and protons at 100 m distance from the shower core.

Dramatic rise of this ratio below 100 GeV implies a monotonic increase of the gap between the energies of electromagnetic and hadronic showers that produce the same amount of Cherenkov light. This effect is demonstrated in Fig. 6. The solid line corresponds to the relation between energies of $`\gamma `$-rays $`E_\gamma `$ and protons $`E_\mathrm{p}`$ defined from equation $`\rho _{100}^{(\gamma )}(E_\gamma )=\rho _{100}^{(\mathrm{p})}(E_\mathrm{p})`$. For comparison, the line “$`E_\gamma =E_\mathrm{p}`$” is also shown. It is seen, in particular, that the gap between $`E_\gamma `$ and $`E_\mathrm{p}`$, which is close to $`2`$ at high energies (e.g. 1 TeV electromagnetic showers versus 2 TeV hadronic showers), becomes significantly larger at very low energies below 10 GeV (e.g. 2 GeV $`\gamma `$-ray events versus 30 GeV hadronic showers).

### 3.2 Field of view and pixel size

One of the principal parameters of the imaging camera is its field of view (FoV). At TeV energies the image centroids of showers detected at distances $`R100\mathrm{m}`$ or beyond are shifted in the focal plane from the center of the camera by $`1^{}`$. Therefore an inner region with diameter $`3^{}`$ provides high detection efficiency for $`\gamma `$-rays, and can be treated as optimal zone for the hardware trigger. On the other hand, the IACT technique requires that the camera should be larger than the trigger zone by about one degree, in order to avoid a distortion of the Cherenkov images because of a limited FoV; for lower $`\gamma `$-ray energies the shift of the image centroid decreases . This effect becomes significant especially at very low energies.

In particular at $`E=10\mathrm{GeV}`$ the shift does not exceed $`0.5^{}`$ (Fig. 7). Thus in the energy range $`E100\mathrm{GeV}`$, the camera can be more compact compared with the conventional (FoV$`4^{}`$) cameras designed for high energies. In fact, a FoV$`3^{}`$ seems a reasonable for point-like or even moderately extended ($`0.5^{}`$) $`\gamma `$-ray sources. This is an important circumstance which keeps the number of the channels of the imager within reasonable ($`10^3`$) limits, especially taking into account that at such low energies the pixel size must be small, $`0.1^{}`$ or less, dictated by two the equally important conditions: effective suppression of the N.S.B., and high image quality (see Sec. 2).

### 3.3 Trigger integration gate

Truncation of the Cherenkov pulse integration gate is another effective way to suppress the N.S.B. Ideally, the integration gate should be comparable with the duration of the time impulse of the Cherenkov light which for 10 GeV $`\gamma `$-rays, detected at 5 km a.s.l. altitude, is less than 5 ns (Fig. 8). The compression of the integration gate from conventional 20 ns to 5 ns reduces the average noise level - the number of photoelectrons produced by the N.S.B. - by a factor of 4, and correspondingly lowers the energy threshold significantly.

Note that both techniques, fast timing and small pixel size, have already been successfully employed by the CAT collaboration . A further suppression of the N.S.B. can be effectively achieved in the stereoscopic approach by requiring simultaneous detection of showers by two or more telescopes. This technique has been convincingly demonstrated by the HEGRA collaboration . Below we will show that, even for telescopes with an aperture as large as $`50\mathrm{m}^2`$, the combination of all three techniques - small pixel size ($`0.1^{}`$), fast timing ($`5\mathrm{ns}`$), and stereoscopy \- allows effective operation of the IACT system in the energy regime as low as several GeV.

### 3.4 Effect of geomagnetic field

An important issue concerning the quality of Cherenkov images at very low energies

is connected with the deflection of the secondary (cascade) electrons in the geomagnetic field (see e.g. Ref. ). Generally, our results agree well with the conclusion of

Patterson and Hillas  that down to $`\gamma `$-ray energies of 100 GeV the geomagnetic field has a rather small effect on the lateral distribution of the Cherenkov radiation as long as the perpendicular component of the magnetic field $`B_{}`$ is less than 0.3 G. One cannot a priory exclude, however, that at very low energies the effect of the geomagnetic field is stronger. The main parameter characterizing the effect of the geomagnetic field is the average energy $`\overline{E_\mathrm{e}}`$ of the electrons radiating Cherenkov light at the shower maximum, rather than the energy of primary particles. In Fig. 9 we compare the parameter $`\overline{E_\mathrm{e}}`$ for showers produced by primary $`\gamma `$-rays at three energies - 10 GeV, 100 GeV, and 1 TeV. It is seen that, at any fixed depth, the average energy of electrons producing Cherenkov light is larger for cascades induced by higher energy $`\gamma `$-rays. However, the average electron energies become quite similar when we compare $`\overline{E_\mathrm{e}}`$ at the depths corresponding to the shower maxima for the given energy of the primary $`\gamma `$-ray photon. This interesting effect is explained by the fact that at low primary energies the maximum of an electromagnetic shower occurs at high altitudes where the density is low, and therefore the electrons need to be more energetic to produce Cherenkov light. In practice this implies that lowering the energy of $`\gamma `$-rays should not result in significant amplification of the geomagnetic effect. In Fig. 10 the density of the Cherenkov light from 10 GeV showers at the observation level are shown for three different values of the perpendicular component of the geomagnetic field - $`B_{}`$=0, 0.5 G, 1 G. The split of the patterns (azimuthal asymmetry) produced by electrons and positrons due to their deflections in opposite directions becomes noticeable only for a very large perpendicular component of the field, $`B_{}`$=1 G. More quantitatively, the effect of the geomagnetic field on the lateral distribution of Cherenkov photons from 10 GeV $`\gamma `$-rays is shown in Fig. 11. It is seen that for a reasonable field $`B_{}0.5\mathrm{G}`$, the effect is less than 15 per cent.

## 4 Detection areas

For a given telescope configuration and arrangement of the IACT array, the energy threshold and the effective detection area of primary $`\gamma `$-rays are determined by the lateral and angular distributions of the Cherenkov light and the hardware trigger conditions.

The configuration of the IACT array assumed in this paper is the following: the system consisting of 5 IACTs is installed at the altitude $`H=5\mathrm{km}`$ a.s.l., four telescopes are located at the corners, and one in the center of a square with a linear size $`d=100\mathrm{m}`$; each telescope has an aperture $`S=50\mathrm{m}^2`$, and is equipped with a $`n=721`$-channel camera of individual pixel size $`0.12^{}`$. This comprises an effective field of view FoV$`=3.2^{}`$. The hardware trigger is organized in a way that requires (1) signals above some critical threshold $`q_0`$ in $`m`$ adjacent pixels (“$`m/nq_0`$”) in each individual telescope (the “local trigger”); (2) detection of a shower by at least 2 telescopes (the “system trigger”). For calculations of the threshold $`q_0`$ we assume a 5 ns trigger gate, and require that the accidental rate caused by the N.S.B. is less than 1 Hz, i.e. less than 1 per cent of the detection rate of cosmic ray electrons (see below).

In Fig. 12 we show the detection areas for $`\gamma `$-rays from 1.5 GeV to 100 GeV calculated for 3 local trigger conditions: $`2/721`$ 6 ph.-e., $`3/721`$ 5 ph.-e. $`4/721`$ 5 ph.-e., assuming a standard N.S.B. flux $`F_{N.S.B.}1.5\times 10^{12}\mathrm{ph}/\mathrm{m}^2\mathrm{sr}\mathrm{s}`$ (see e.g. ). It is seen that for a chosen pixel size of $`0.12^{}`$, a multiplicity greater than two ($`m>2`$) actually reduces the detection area, even though the requirement of higher multiplicity allows a lower trigger threshold $`q_0`$. Therefore all calculations presented below were performed for $`m=2`$.

In Fig. 13 the detection area of the IACT array for a somewhat higher threshold, $`q_0=8`$ ph.-e., is shown in order to demonstrate how sensitive the detection area is to the choice of $`q_0`$, especially at low, sub-10 GeV energies (compare curves 1 in Fig. 12 and Fig. 13). For comparison, we show also the overall detection area of 5 telescopes if they would be located at large, $`100\mathrm{m}`$ distances from each other and operating independently in a stand alone mode, i.e. being not integrated in the system trigger. It is seen that at energies $`E3\mathrm{GeV}`$ the array of independent telescopes has a larger detection area, but the difference compared with the array operating in the stereoscopic mode at energies $`E5\mathrm{GeV}`$ does not exceed a factor of two. On the other hand, the analysis of images obtained in the stereoscopic mode provides much better suppression of the background caused by the cosmic ray protons and electrons which effectively results in a significant improvement of not only the quality of the data, but also of the flux sensitivity of the instrument.

At energies below 100 GeV, the cosmic ray background detected by the system of IACTs operating in the stereoscopic mode is essentially dominated by electromagnetic showers produced by cosmic ray electrons. This is seen in Fig. 14 where the detection rates of both cosmic ray protons and electrons are also shown. The strong dominance of the electronic component is explained by the combination of several effects, in particular (i) the large, up to a factor of 10, difference between the energies of electrons and protons producing the same amount of Cherenkov light, (ii) the high altitude of observations, (iii) the compact Cherenkov images of electromagnetic showers compared with hadronic showers, (ii) the noticeable increase (approximately $`E^{0.5}`$) of the electron-to-proton ratio of cosmic rays down to $`E10\mathrm{GeV}`$.

In Fig.14 we show also the detection rate of $`\gamma `$-rays from a point source with a power-law spectrum $`\mathrm{d}J_\gamma /\mathrm{d}EE^{2.5}`$, and integral flux $`J_\gamma (1\mathrm{GeV})=3\times 10^7\mathrm{ph}/\mathrm{cm}^2\mathrm{s}`$. The latter is somewhat larger than the fluxes of most of the EGRET sources . This implies that the $`\gamma `$-ray detection rate shown in Fig. 14 should be considered as an upper limit for “standard” EGRET sources. This curve still lies significantly below the rate of detection of cosmic ray electrons. However, for a point-like source the electron background can be reduced significantly if we select showers arriving from the direction of the $`\gamma `$-ray source.

The detection rates of electromagnetic showers within the detector’s energy-dependent Point Spread Function (PSF), given by the angle $`\varphi `$, which is the half-angle of the cone around the source direction containing $`67\%`$ of events at energy $`E`$, are determined as follows:

$$R_\gamma =\frac{\mathrm{d}J}{\mathrm{d}E}A_{\mathrm{eff}}(E)\kappa _\gamma ,$$
(1)
and

$$R_\mathrm{e}=\frac{\mathrm{d}J_\mathrm{e}}{\mathrm{d}E\mathrm{d}\mathrm{\Omega }}A_{\mathrm{eff}}(E)\mathrm{\Omega },$$
(2)
for $`\gamma `$-rays and cosmic ray electrons, respectively, where $`\mathrm{d}J/\mathrm{d}E`$ is the differential flux of $`\gamma `$-rays from a point source, $`\mathrm{d}J_\mathrm{e}/\mathrm{d}E\mathrm{d}\mathrm{\Omega }`$ is the differential flux of cosmic ray electrons per solid angle, $`A_{\mathrm{eff}}`$ is the detection area for electromagnetic showers, $`\mathrm{\Omega }=2\pi (1\mathrm{cos}\varphi )`$, and $`\kappa _\gamma =0.67`$ is (by definition) the $`\gamma `$-ray acceptance.

For the chosen configuration of the IACT array, the Monte-Carlo calculations of the effective detection area $`A_{\mathrm{eff}}`$ (Fig. 15) in the interval from 1.5 GeV to 100 GeV can be presented in the following form

$$A_{\mathrm{eff}}=8.5E^{5.2}[1+(E/5\mathrm{GeV})^{4.7}]^1\mathrm{m}^2,$$
(3)
where E is the energy of a $`\gamma `$-ray or electron in units of GeV.

This presentation shows a strong energy-dependence of the detection area at energies below 10 GeV, $`A_{\mathrm{eff}}(E)E^{5.2}`$, but at higher energies it gradually turns into a slow increase with energy, $`A_{\mathrm{eff}}(E)E^{0.5}`$ Close to 1 TeV the detection area actually becomes constant, in essence because of the limited field of the view ($`3^{}`$) of the camera.

The stereoscopic approach allows the determination of the arrival direction of primary $`\gamma `$-rays on an event-by-event basis. The determination of the arrival direction of primary $`\gamma `$-rays is described in ref., and its practical implementation in the case of the HEGRA IACT system can be found in ref. . In Fig. 16 we show the Monte-Carlo calculations of the PSF which can be approximated in the simple form

$$\varphi =0.8(E/1\mathrm{GeV})^{0.4}\mathrm{degree}.$$
(4)
Despite the small pixel size $`0.1^{}`$, the angular resolution of the 5 GeV IACT array at energies $`E100\mathrm{GeV}`$ is significantly poorer than the resolution of “100 GeV” threshold instruments, which could be as good as $`0.1^{}`$, even for a larger pixel-size of about $`0.25^{}`$ . Since for both categories of instruments the number of photoelectrons, or the so-called size of the image, from a detected $`\gamma `$-ray photon are comparable ($`100(E/E_{\mathrm{th}})`$ photoelectrons), the lower performance of the IACT technique is rather an intrinsic feature of Cherenkov images at very low energies (the images are less elongated and less regular). This effect can be seen also in the so-called alpha-distribution of images in a single telescope, where the alpha parameter is indicative of the orientation of the image in the camera (see e.g. ref.): while at TeV energies most of $`\gamma `$-rays from a point source have an angle alpha less than 5-8 deg , the alpha distribution is significantly broader in the energy region around 10 GeV (see Fig. 17). The suppression of the cosmic ray background at such low energies becomes correspondingly less effective. The characteristic values of the acceptance of the isotropic cosmic ray showers $`\kappa _{\mathrm{CR}}`$, and the point-source $`\gamma `$-ray showers $`\kappa _\gamma `$, as well as the so-called Q-factor, $`Q=\kappa _\gamma /\kappa _{\mathrm{CR}}^{1/2}`$, which characterizes the improvement of the signal-to-noise ratio after application of the image cuts, are shown in Table 1. It is seen that the best improvement of the signal-to-noise ratio is achieved for an alpha cut at $`2025^{}`$ which gives a rather modest Q-factor, $`Q1.4`$. In contrast, at 10 GeV the stereoscopic measurements allow the determination of the shower direction with an accuracy of $`0.3^{}`$ (Fig. 16). Therefore the $`\gamma `$-ray signal could be improved by a factor of $`Q=\kappa _\gamma /\kappa _{\mathrm{CR}}^{1/2}=\kappa _\gamma (\mathrm{\Psi }/2\varphi )3.35`$, where the efficiency of the rejection showers from cosmic ray electrons $`\kappa _{\mathrm{CR}}(2\varphi /\mathrm{\Psi })^2`$ is determined by the FoV of the imager, $`\mathrm{\Psi }3^{}`$, and the angular resolution $`\varphi 0.3^{}`$. Thus, despite the smaller (by a factor of two) detection area of the array operating in the stereoscopic mode compared with the overall area of 5 independent IACTs (see Fig. 13), the stereoscopic array would have at least by a factor of 2 better sensitivity, even disregarding other advantages of the stereoscopic approach, in particular, the complete removal of the hadronic background.

## 5 Detection rates and the energy threshold

The differential detection rates of $`\gamma `$-rays from a point source, calculated for the collection area given by Eq. (3) and assuming pure power-law spectrum in the form

$$\frac{\mathrm{d}J}{\mathrm{d}E}=10^7(E/1\mathrm{GeV})^\alpha \mathrm{cm}^2\mathrm{s}^1\mathrm{GeV}^1,$$
(5)
are shown in Fig. 18a. Eq. (5) implies that, independent of the spectral index $`\alpha `$, the differential $`\gamma `$-ray flux is normalized at 1 GeV to $`10^7\mathrm{cm}^2\mathrm{s}^1\mathrm{GeV}^1`$. The latter corresponds to the typical flux of “standard” EGRET sources.

In Fig. 18b we show the differential detection rates of $`\gamma `$-rays with a hard power-law spectrum with $`\alpha =2`$ and an exponential cutoff at $`E_0`$:

$$\frac{\mathrm{d}J}{\mathrm{d}E}E^2\mathrm{exp}(E/E_0),$$
(6)
for 4 different values of $`E_0`$, and assuming the same absolute flux normalization at 1 GeV as in Fig. 18a.

It is seen from Figs.18a and 18b that for a large variety of $`\gamma `$-ray spectra the peak of the detection rate appears in a rather narrow band between 4 and 6 GeV. Now, defining the energy threshold as the energy at which the differential $`\gamma `$-ray detection rate reaches to its maximum, we may conclude that the suggested IACT array has an effective energy threshold of about 5 GeV.

In Fig. 18a,b we show also the differential detection rate of cosmic ray electrons within the cone limited by the PSF of the instrument given by Eq. (4). The energy spectrum of cosmic ray electrons is shown in Fig. 19. At energies above 10 GeV the differential spectrum is very steep with a power-law index $`\alpha _\mathrm{e}3.2`$. Below 10 GeV it becomes flatter. Within the uncertainties of the measured fluxes, the electron spectrum can be approximated in the entire region from several GeV to 1 TeV by the following function shown by the solid line in Fig. 19 :

$$\frac{\mathrm{d}J_\mathrm{e}}{\mathrm{d}E\mathrm{d}\mathrm{\Omega }}=1.36\times 10^7E^1[1+(E/5\mathrm{GeV})^{2.2}]^1\mathrm{m}^2\mathrm{s}^1\mathrm{sr}^1\mathrm{GeV}^1.$$
(7)
In Fig. 20 we show integral detection rates $`R_\gamma (\mathrm{E})`$ for power-law $`\gamma `$-ray spectra represented by Eq.(5). It is seen that for a relatively flat $`\gamma `$-ray spectrum with $`\alpha _\gamma 2`$ and for an integral flux above 1 GeV of $`10^7\mathrm{ph}/\mathrm{cm}^2\mathrm{s}`$ (this approximately corresponds to the total, i.e. pulsed plus unpulsed, flux from the Crab), the detection rate of $`\gamma `$-rays from the EGRET sources can be as high as 6 events per second, against the cosmic-ray background rate of about 25 events per sec caused by cosmic-ray electrons. This implies that an observation time of approximately 20-30 sec would be sufficient to detect a statistically significant signal from such a source. Remarkably, for the brightest persistent $`\gamma `$-ray source, the Vela pulsar with photon index $`\alpha _\gamma 1.7`$ from 100 MeV to 10 GeV, and the integral flux $`J_\gamma (1\mathrm{GeV})1.5\times 10^6\mathrm{ph}/\mathrm{cm}^2\mathrm{s}`$, the detection rates would exceed 100 events per 1 sec. Thus a statistically significant signal from the source could be obtained during an observation time less than 1 sec ! For the given normalization of the differential flux at 1 GeV, $`10^7\mathrm{ph}/\mathrm{cm}^2\mathrm{s}\mathrm{GeV}`$, the detection of sources with steep $`\gamma `$-ray spectra would require significantly longer exposure. Even so, the time needed for detection of sources with very steep power-law spectra with an index $`\alpha _\gamma 3`$ (like curve 4 in Fig. 18a), or with a sharp, e.g. exponential cutoff at a few GeV (like the curve 4 in Fig. 18b), does not significantly exceed 1 h.

## 6 Flux Sensitivity

The curves in Fig. 14 correspond to the detection rates before the image analysis. Remarkably, even after such effective rejection of hadronic showers at the trigger level, there still remains room for further suppression of the background from cosmic ray protons and nuclei by analyzing the shapes of the Cherenkov images of the detected showers. In Fig. 21 we show the so-called mean-scaled width parameter distribution of showers which have already passed the hardware trigger condition. This parameter represents the mean value of the width parameter measured by all telescopes and normalized to the impact distances and the image amplitudes . It is seen that the distributions of the electromagnetic and hadronic showers are rather well separated. The efficiencies of the acceptance of both type of showers for different mean scaled width cuts are presented in Table 2. We see that even very loose cuts at the level of $`<w>=1.3`$ provide suppression of the hadronic showers by a factor of 5, while the $`\gamma `$-ray acceptance can be as high as 90 per cent. This implies that after such a loose cut, which practically does not reduce the $`\gamma `$-ray (or CR electron) statistics, we may push the detection rates of cosmic ray protons and nuclei further down, and thus make $`\gamma `$-ray detection in the entire energy region below 100 GeV essentially free from hadronic background.

This makes the calculations for the differential flux sensitivity of the instrument straightforward and simple and, more importantly, there is no need to specify the spectrum of primary $`\gamma `$-rays. Indeed, because the showers produced by electrons are very similar to $`\gamma `$-ray showers222Actually there are some differences. In particular the primary electrons start to produce Cherenkov light earlier, but in this paper we will ignore these effects., the condition of detection of a $`\gamma `$-ray signal with statistical significance $`m`$-sigma in the energy interval $`[E\mathrm{\Delta }E,E+\mathrm{\Delta }E]`$, $`N(E)\mathrm{\Delta }E=m\sqrt{N_\mathrm{e}\mathrm{\Delta }E}`$, provided that the number of detected $`\gamma `$-rays $`N_{\mathrm{min}}=J(E)2\mathrm{\Delta }EA_{\mathrm{eff}}(E)T10`$, gives the minimum detectable differential flux for the observation time $`T`$:

$$J_{\mathrm{min}}(E)\frac{m}{\kappa _\gamma }\frac{\varphi }{(A_{\mathrm{eff}}(E)T)^{1/2}}\sqrt{\frac{2\pi J_\mathrm{e}(E)}{E}}$$
(8)
where it is assumed that $`\mathrm{\Delta }E=E/4`$, which corresponds to a rather conservative 25 per cent accuracy of reconstruction of the energy of the primary electron or $`\gamma `$-ray photon. The results of calculations for the so-called spectral energy distribution (SED), $`E^2J(E)`$, based on Eqs.(3),(4) and (7) requiring 3-sigma detection ($`m=3`$) at each energy $`E`$, are presented in Fig. 22 for 2 different observation times, $`T=`$1 h and 25 h.

In the same figure we present also the power-low fluxes of $`\gamma `$-rays represented by Eq.(6). It is seen that 1 hour observations by the IACT array would be sufficient to detect a statistically significant signal from a “standard” EGRET source even at the presence of an exponential cutoff in the $`\gamma `$-ray spectrum as low as 3 GeV. In the case of a cutoff at 10 GeV or higher energies, the detection time ($`t1/J^2`$) could be reduced to $`1`$ min.

For comparison in Fig. 22 we show the expected sensitivity of GLAST for an observation time $`T=30\mathrm{days}`$. Since at energies above several GeV GLAST will operate at almost background free conditions (for point-like sources), the flux sensitivity is determined by the photon statistics, $`N_{\mathrm{min}}=10`$ (see e.g. ). Note that at energies above 5-10 GeV, 1 hour observation time by the IACT array could provide better sensitivity that the minimum detectable fluxes achievable by GLAST during 1 month of continuous observations. Moreover, even very short observations by the IACT array can give unusually rich (for $`\gamma `$-ray astronomical standards) photon statistics over the whole energy region from few GeV to 100 GeV; the number of detected $`\gamma `$-rays exceeds 100 at each energy interval $`E\pm E/4`$ (see Fig. 23). This would guarantee an appropriate $`\gamma `$-ray spectroscopy with energy resolution of about 20-25 per cent below 10 GeV, and better than 15 per cent at higher energies.

The flux sensitivity shown in Fig. 22 is obtained under a very robust and to a large extent non-standard condition which requires detection of a signal with at least 3-sigma significance in each energy band with width $`E/2`$ centered on $`E`$, provided that the number of detected $`\gamma `$-rays in this band exceeds 10. Note that this definition of sensitivity does not require any knowledge about the shape of the spectrum of the primary $`\gamma `$-rays.

In Fig.24 we present the flux sensitivities determined in a more traditional way, namely requiring 5-sigma detection of $`\gamma `$-rays above the given energy, $`J_{\mathrm{min}}(\mathrm{E})`$. This definition of the integral flux sensitivity obviously requires an assumption about the shape of the energy spectrum. The curves shown in Fig. 24 are calculated for $`T=25\mathrm{h}`$ observation time, assuming power-law spectra of $`\gamma `$-rays with photon indices $`\alpha _\gamma =`$1.5 (dotted curve), 2 (dot-dashed curve), and 2.5 (dashed curve). The expected GLAST sensitivity shown by the solid curve corresponds to 1 year of continuous observations of the source.

## 7 Discussion

Results presented in this paper show that a stereoscopic array of large, 20 m diameter class imaging Cherenkov telescopes installed at very high mountain altitudes could effectively enter into the domain of satellite-borne $`\gamma `$-ray astronomy. A 5 GeV energy threshold array of IACTs at 5 km a.s.l. - hereafter 5@5 - could provide a deeper probe of $`\gamma `$-ray sources compared with GLAST - the most powerful current satellite-borne $`\gamma `$-ray project. However, the scientific goals of these instruments are essentially different. While GLAST with its almost $`2\pi `$ steradian field of view can provide very effective simultaneous monitoring of a very large number (hundreds or even thousands) quasi-stable $`\gamma `$-ray sources, as well as a study the galactic and extragalactic components of the diffuse $`\gamma `$-ray background radiation, 5@5 has an obvious advantage for the search and study of highly variable or transient $`\gamma `$-ray sources. The flux sensitivity of this instrument at 5 GeV of about $`2\times 10^{11}\mathrm{erg}/\mathrm{cm}^2\mathrm{s}`$ (see Fig. 22) would allow the detection of any $`\gamma `$-ray flare with apparent luminosity $`2\times 10^{39}(\mathrm{d}/1\mathrm{Mpc})^2`$, lasting only 1 h, where $`d`$ is the distance to the source. Of special interest are the gamma-ray blazars detected by EGRET (see e.g. ). A detailed study of the time structure of $`\gamma `$-radiation for these highly variable objects on timescales of several minutes by 5@5 would provide unique information about the relativistic non-thermal processes in astrophysical jets. An effective operation of this instrument in the sub-10 GeV regime guarantees detection of $`\gamma `$-rays arriving from cosmological distances up to $`z3`$ or so, for which the intergalactic medium becomes almost transparent. The dynamic range from several GeV to 100 GeV would allow important cosmological measurements, in particular a study of the diffuse ultraviolet extragalactic background by detecting intergalactic $`\gamma \gamma `$ absorption features in the spectra of $`\gamma `$-rays below 100 GeV. The confusion problem (spectral cutoff due to the internal or extragalactic absorption ?) at such redshifts can probably be overcome by simultaneous observations at optical and X-ray wavelengths.

5@5 can also be effectively used for the study of galactic transient sources, in particular for the detection of short time ($`1\mathrm{day}`$) $`\gamma `$-ray activity expected during synchrotron radio flares of microquasars.

And finally, 5@5 can serve as a very powerful instrument for a study of the phenomenon of Gamma Ray Bursts (GRBs). If the spectra of GRBs extend to high energies without abrupt cutoffs up to several GeV, which is the case at least for some of GRBs , then the sensitivity of 5@5 would allow very detailed studies of the spectral and temporal features of GRBs in this extremely important energy region. In Fig. 25 we show the minimum time $`t_{\mathrm{min}}`$ required for detection of GeV $`\gamma `$-ray flares with a given energy flux $`f(E)`$ at 4 different $`\gamma `$-ray energies: 2 GeV, 5 GeV, 10 GeV, and 50 GeV. The calculations correspond to the 3-sigma signal at each energy $`E`$ within the interval $`[E0.25E,E+0.25E]`$, provided that the number of detected $`\gamma `$-rays exceeds 10. In the regime of low fluxes, typically $`f10^9\mathrm{erg}/\mathrm{cm}^2\mathrm{s}`$, the $`\gamma `$-rays are detected in the presence of the heavy background induced by cosmic ray electrons. Therefore $`t_{\mathrm{min}}f^2`$. For fluxes larger than $`10^9\mathrm{erg}/\mathrm{cm}^2\mathrm{s}`$, the detection occurs under almost background-free conditions, and therefore $`t_{\mathrm{min}}f^1`$.

The results shown in Fig. 25 demonstrate the capability of 5@5 for detection of GeV counterparts of GRBs. The detection of $`5`$ GeV episodic events with typical GRB fluxes between $`10^8`$ and $`10^6\mathrm{erg}/\mathrm{cm}^2\mathrm{s}`$ would require only $`0.1\mathrm{s}`$ observation time. Thus it would be possible to monitor the spectral evolution of the source with a typical GRB duration from several second to 100 seconds. Remarkably, even for fluxes as low as $`10^{10}\mathrm{erg}/\mathrm{cm}^2\mathrm{s}`$, the required exposure time does not exceed 100 sec. This implies that 5@5 could serve as a unique tool to study GRBs in the late stages of evolution, i.e. during the afterglows.

5@5 is a detector with a small field of view. Therefore it requires special strategies for the search and study of multi-GeV $`\gamma `$-ray emitters. The proximity of this energy region to the energy range covered by EGRET suggests that almost all (more than 300) EGRET sources should also be detected by 5@5. The typical observation time for detection of a “standard” EGRET source would not exceed 1 hour, even if the spectrum of $`\gamma `$-rays cuts off at energies of several GeV. The full overlap of the energy range of this instrument with the energy domain of GLAST would make the latter a “the best guide” for developing a strategy for the study of persistent galactic and extragalactic objects. Generally, all sources seen by GLAST can be potential targets for observations with the 5@5. These observations with very large photon statistics - not achievable by GLAST - could provide detailed studies of the spectral and temporal features of $`\gamma `$-ray sources in the multi-GeV region. For highly variable objects like blazars or galactic sources with relativistic jets, a multi-wavelength approach including observations with radio, optical and X-ray detectors would be very important. These observations would not only inform about the pre-flaring or flaring states of the sources, but also would provide complementary information for understanding and comprehensive modeling of the physical processes in these objects.

A special strategy should be developed for the search for GeV radiation from GRBs during and after the main event. Apparently, prompt information (within 10 sec or so) from the new generation GRB detectors like SWIFT (and possibly also from GLAST) would be needed, containing the angular coordinates of an event with an accuracy better than $`1^{}`$. It would be very worthwhile to have also a nearby ground-based prompt optical telescope like ROTSE . In their turn, the telescopes of 5@5 should be rather fast in order to be directed to the source not later than 1 minute after receiving the alarm from these detectors.

Because of effective rejection of hadronic showers by 5@5, the cosmic-ray background below 100 GeV is dominated by the showers from cosmic-ray electrons. This component of electromagnetic showers remains a part of the background which can be hardly removed, and thus it is the most serious limiting factor of flux sensitivities, especially for extended sources. On the other hand, these electromagnetic showers with a known flux and spectrum of cosmic ray electrons, measured up to energies 1 TeV, can be used for absolute energy calibration of the instrument.

Energy calibration with the aid of the cosmic-ray electrons provides a unique tool for the continuous (on-line) control of the characteristics of the detector (e.g. the energy threshold, the detection area, etc.) during the observations. For example, a 10 min exposure will be enough for the detection of hundreds of electrons at any energy $`E`$ within $`\pm E/4`$ in the entire dynamical region from few GeV to 100 GeV (see Fig. 23). It is difficult to overestimate the significance of such calibration and control, especially for the study of the spectral characteristics of highly variable $`\gamma `$-ray sources on sub-hour timescales.

The basic elements of the suggested 5@5 detector are the large optical reflectors and the multichannel high resolution cameras. Presently, 20 m diameter alt-azimuth mounts with the required precision of about 1 arcminute could be designed and built by many companies specialized in the construction of large radio dishes. The area of the optical reflector could be composed of several hundreds to thousand $`1\mathrm{m}`$ diameter glass mirrors with protective quartz coating, quite similar to the mirrors used in the current or planned imaging Cherenkov telescope projects. The general requirements on the imagers are a relatively large ($`3^{}`$ or so) field of view with a pixel size of about $`0.1^{}`$. This implies less than 1 thousand fast channels. Such a camera with similar parameters has been already built and successfully operated as part of the CAT imaging Cherenkov telescope .

The decrease of the energy threshold of the imaging atmospheric Cherenkov technique to several GeV would depend to a large extent on the availability of exceptional sites with a dry and transparent atmosphere at an altitude as high as 5 km. Nature does provide us with such an extraordinary site - The Llano de Chajnantor in the Atacama desert in Northern Chile. This site with its very arid atmosphere was recently chosen for the installation of one of the most powerful future astronomical instruments - the Atacama Large Millimeter Array (ALMA), a funded US-European project which will consist of 64 12-meter radio antennas with spacing from approximately 150 meters to 10 km (see http://www.alma.nrao.edu/). The large flat area on that site could certainly accommodate an additional Cherenkov telescope array as well which requires a relatively compact area with a radius of about 100 m. Another attractive feature of this site seems to be an adequate infrastructure which will be built up during the next several years for the ALMA project. The foreseen technological developments of ALMA, concerning in particular the construction of very large antennas operating in robotic or semi-robotic mode, could help very much in the design of the telescopes of 5@5. Moreover, the neighboring Cerro Toco site  offers suitable areas at even higher altitudes, $`H5.2`$ and $`5.6`$ km a .s.l., with the same infrastructure advantages.

Acknowledgments

We thank the anonymous referee for her/his critical comments and remarks which helped us to improve the paper significantly. HQ is grateful for the award of a Presidential Chair in Science (Chile). His research is partially funded by FONDECYT Grant No. 8970009.


format: text

text:
# Universal gap fluctuations in the superconductor proximity effect

A normal metal in the proximity of a superconductor acquires characteristics that are typical of the superconducting state . One of those characteristics is that the quasiparticle density of states vanishes at the Fermi energy. This superconductor proximity effect is most pronounced in a confined geometry, such as a thin metal film or metal grain, or a semiconductor quantum dot. In that case, provided the scattering in the normal metal is chaotic, no excitations exist within an energy gap $`E_\mathrm{g}\mathrm{}/\tau `$, where $`\tau `$ is the typical time between collisions with the superconductor .

If the coupling to the superconductor is weak (as for the point contact coupling of Fig. 1), the functional form of the density of states becomes independent of microscopic properties of the normal metal, such as the shape, dimensionality, or mean free path. Weak coupling means that $`\tau `$ is much bigger than the time $`\tau _{\mathrm{erg}}`$ needed for ergodic exploration of the phase space in the normal region. For a point contact with $`N1`$ propagating modes at the Fermi level $`\epsilon =0`$, the density of states has a square root singularity at the excitation gap ,

$$\rho _{\mathrm{mf}}(\epsilon )=\frac{1}{\pi }\sqrt{\frac{\epsilon E_\mathrm{g}}{\mathrm{\Delta }_\mathrm{g}^3}}.$$
(1)
For a ballistic point contact and in the absence of a magnetic field, $`E_\mathrm{g}=cN\delta `$ and $`\mathrm{\Delta }_\mathrm{g}=c^{}N^{1/3}\delta `$, where $`c=0.048`$ and $`c^{}=0.068`$ are numerical constants and $`\delta `$ is the mean level spacing in the normal metal when it is decoupled from the superconductor.

Equation (1) was obtained in a self-consistent diagrammatic perturbation theory that uses $`\tau \delta /\mathrm{}N^1`$ as a small parameter. Such a mean-field theory provides a smoothed density of states for which energies can only be resolved on the scale of $`\mathrm{}/\tau N\delta `$, not on smaller energy scales, and is unable to deal with mesoscopic sample-to-sample fluctuations of the excitation gap. Mesoscopic fluctuations arise, e.g., upon varying the shape of a quantum dot or the impurity configuration in a metal grain. The lowest excited state $`\epsilon _1`$ fluctuates from sample to sample around the mean field value $`E_\mathrm{g}`$, with a probability distribution $`P(\epsilon _1)`$. It is the purpose of this paper to go beyond mean field theory and to study the mesoscopic fluctuations of the excitation spectrum close to $`E_\mathrm{g}`$. Our main result is that the gap distribution $`P(\epsilon _1)`$ is a universal function of the rescaled energy $`x=(\epsilon _1E_\mathrm{g})/\mathrm{\Delta }_\mathrm{g}`$, in a broad range $`|x|N^{2/3}`$. The Fermi level itself ($`\epsilon =0`$) falls outside this range, which is why the universal gap distribution was not found in a recent related study . Our main findings are illustrated in Fig. 2.

We first consider the gap distribution in the absence of a magnetic field, and then include a time-reversal symmetry breaking magnetic field. Starting point of our calculation is the effective Hamiltonian

$$=\left(\begin{array}{cc}H& \pi WW^{}\\ \pi WW^{}& H^{}\end{array}\right).$$
(2)
Here $`H`$ is an $`M\times M`$ Hermitian matrix representing the Hamiltonian of the isolated quantum dot, and $`W`$ is an $`M\times N`$ matrix that describes the coupling to the superconductor via an $`N`$-mode point contact. For a ballistic point contact, $`W_{mn}=\pi ^1\delta _{mn}(M\delta )^{1/2}`$ . The number $`M`$ is sent to infinity at the end of the calculation . The effective Hamiltonian is a valid description of the low-lying excitations if the Thouless energy $`N\delta `$ is much smaller than the order parameter $`\mathrm{\Delta }`$ of the bulk superconductor. In the absence of a magnetic field, the matrix $`H`$ is symmetric. To describe an ensemble of chaotic quantum dots (or disordered metal grains), we take $`H`$ from the Gaussian orthogonal ensemble (GOE) of random-matrix theory ,

$$𝒫(H)\mathrm{exp}\left(\frac{\pi ^2}{4\delta ^2M}\mathrm{Tr}H^2\right).$$
(3)
Calculation of the density of states of $``$ using perturbation theory in $`N^1`$ yields the result (1) discussed in the introduction. Our problem is to go beyond perturbation theory and find the probability distribution $`P(\epsilon _1)`$ of the lowest positive eigenvalue $`\epsilon _1`$ of the Hamiltonian (2).

We have solved this problem numerically by exact diagonalization of the effective Hamiltonian $``$. Before presenting these results, we first describe an entirely different approach, that leads to an analytical prediction for the gap distribution. We invoke the universality hypothesis of random-matrix theory, that the local spectral statistics of a chaotic system depends only on the symmetry properties of the Hamiltonian, and not on microscopic properties. This universality hypothesis has been proven for a broad class of Hamiltonians in the bulk of the spectrum , but is believed to be valid near the edge of the spectrum as well. A proof exists for so-called trace ensembles, having $`𝒫(H)\mathrm{exp}[\text{tr}f(H)]`$, with $`f`$ an arbitrary polynomial function .

The mean-field density of states near the edge can be written in the form

$$\rho _{\mathrm{mf}}(\epsilon )=\frac{1}{a}\left(\frac{\epsilon b}{a}\right)^p,\epsilon >b.$$
(4)
According to the universality hypothesis, the spectral statistics near the edge, in rescaled variables $`(\epsilon b)/a`$ depends only on the exponent $`p`$ and on the symmetry index $`\beta `$ \[$`\beta =1`$ ($`2`$) in the presence (absence) of time-reversal symmetry\]. Generically, $`p`$ is either $`1/2`$ (soft edge) or $`1/2`$ (hard edge). For our problem, we have $`\beta =1`$, $`p=1/2`$, $`a=\pi ^{2/3}\mathrm{\Delta }_\mathrm{g}`$, $`b=E_\mathrm{g}`$, cf. Eq. (1). The corresponding gap distribution is given by

$`P(\epsilon )`$ $`=`$ $`{\displaystyle \frac{d}{d\epsilon }}F_1\left[(\epsilon E_\mathrm{g})/\mathrm{\Delta }_\mathrm{g}\right],`$ (5)
$`F_1(x)`$ $`=`$ $`\mathrm{exp}\left({\displaystyle \frac{1}{2}}{\displaystyle _{\mathrm{}}^x}[q(x^{})+(xx^{})q^2(x^{})]𝑑x^{}\right).`$ (6)
The function $`q(x)`$ is the solution of

$$q^{\prime \prime }(x)=xq(x)+2q^3(x),$$
(7)
with asymptotic behavior $`q(x)\mathrm{Ai}(x)`$ as $`x\mathrm{}`$ \[$`\mathrm{Ai}(x)`$ being the Airy function\].

The distribution (5) is shown in Fig. 3 (solid curve). It is centered at a positive value of $`x=(\epsilon _1E_\mathrm{g})/\mathrm{\Delta }_\mathrm{g}`$, meaning that the average gapsize $`\epsilon _1`$ is about $`\mathrm{\Delta }_\mathrm{g}`$ bigger than the mean-field gap $`E_\mathrm{g}`$. For small $`x`$ there is a tail of the form

$$P(x)\frac{1}{4\sqrt{\pi }|x|^{1/4}}\mathrm{exp}\left(\frac{2}{3}|x|^{3/2}\right),x1.$$
(8)
Non-universal corrections to the distribution (5) become important for energy differences $`|\epsilon E_\mathrm{g}|E_\mathrm{g}`$, hence for $`|x|N^{2/3}`$. Since the width of the gap distribution is of order unity in the variable $`x`$, the probability to find a sample with an excitation gap in the non-universal regime is exponentially small.

In order to verify our universality hypothesis, we compare Eq. (5) with the results of an exact diagonalization of the Hamiltonian (2). As one can see in Fig. 3, the numerical data are in good agreement with the analytical prediction. The small deviations can be attributed to the finiteness of $`N`$ and $`M`$ in the numerics.

Let us now consider the effect of a weak magnetic field on the gap distribution. In the effective Hamiltonian, the presence of a magnetic field is described by replacing $`H`$ by

$$H(\alpha )=H+i\alpha A.$$
(9)
Here $`A`$ is an $`M\times M`$ real antisymmetric matrix, whose off-diagonal elements have the same variance as those of $`H`$. The parameter $`\alpha `$ is proportional to the magnetic field,

$$M\alpha ^2=\eta \left(\frac{\mathrm{\Phi }}{\mathrm{\Phi }_0}\right)^2\frac{\mathrm{}}{\tau _{\mathrm{erg}}\delta },$$
(10)
where $`\mathrm{\Phi }`$ is the magnetic flux through the quantum dot, $`\mathrm{\Phi }_0=h/e`$ is the flux quantum, and $`\eta `$ is a non-universal numerical constant . The case $`\alpha =0`$ corresponds to the GOE that we considered above; the case $`\alpha =1`$ corresponds to the Gaussian unitary ensemble (GUE) of fully broken time-reversal symmetry.

The effect of a magnetic field on the density of states in mean-field theory is known . The square-root singularity (1) near the gap still holds, but the magnitude of the gap is reduced. The critical flux $`\mathrm{\Phi }_\mathrm{c}`$ at which $`E_\mathrm{g}=0`$ and hence the proximity effect is fully suppressed is given by

$$M\alpha ^2N\mathrm{\Phi }_\mathrm{c}\mathrm{\Phi }_0\sqrt{\frac{N\tau _{\mathrm{erg}}\delta }{\mathrm{}}}.$$
(11)
This is a much larger flux than the flux $`\mathrm{\Phi }_{\mathrm{bulk}}`$ at which the spectral statistics in the bulk of the spectrum crosses over from GOE to GUE, which is given by

$$M\alpha ^21\mathrm{\Phi }_{\mathrm{bulk}}\mathrm{\Phi }_0\sqrt{\frac{\tau _{\mathrm{erg}}\delta }{\mathrm{}}}.$$
(12)
We will now argue that the characteristic flux $`\mathrm{\Phi }_{\mathrm{edge}}`$ for the spectral statistics at the edge of the spectrum is intermediate between $`\mathrm{\Phi }_\mathrm{c}`$ and $`\mathrm{\Phi }_{\mathrm{bulk}}`$. We consider the effect of the magnetic field on the lowest eigenvalue $`\epsilon _1`$ of $``$ to second order in perturbation theory,

$$\delta \epsilon _1=\underset{j1}{}\alpha ^2\frac{|1|𝒜|j|^2}{\epsilon _1\epsilon _j},𝒜=i\left(\begin{array}{cc}A& 0\\ 0& A\end{array}\right).$$
(13)
Since typically $`|1|𝒜|2|^2M\delta ^2/\pi ^2`$ and $`\epsilon _2\epsilon _1\mathrm{\Delta }_\mathrm{g}`$, we see that the effect of level repulsion from the neighboring level $`\epsilon _2`$ on the lowest level $`\epsilon _1`$ becomes comparable to $`\mathrm{\Delta }_\mathrm{g}N^{1/3}\delta `$ if

$$M\alpha ^2N^{2/3}\mathrm{\Phi }_{\mathrm{edge}}\mathrm{\Phi }_0\sqrt{\frac{N^{2/3}\tau _{\mathrm{erg}}\delta }{\mathrm{}}}.$$
(14)
The terms in Eq. (13) with $`j1`$ give a uniform shift of all low-lying levels, and hence do not affect the fluctuations. For $`N1`$ the flux scale (14) for breaking time-reversal symmetry at the edge of the spectrum is much smaller than the critical flux $`\mathrm{\Phi }_\mathrm{c}`$ needed to suppress the proximity effect. What is needed is $`N^{2/3}N`$. This condition is difficult to satisfy in a numerical calculation. The analytical prediction for fully broken time-reversal symmetry is

$`P(\epsilon )`$ $`=`$ $`{\displaystyle \frac{d}{d\epsilon }}F_2[(\epsilon E_\mathrm{g})/\mathrm{\Delta }_\mathrm{g}],`$ (15)
$`F_2(x)`$ $`=`$ $`\mathrm{exp}\left({\displaystyle _{\mathrm{}}^x}(xx^{})q^2(x^{})𝑑x^{}\right).`$ (16)
This curve is shown dashed in Fig. 3. The tail for small $`x`$ is now given by

$$P(x)\frac{1}{8\pi |x|}\mathrm{exp}\left(\frac{4}{3}|x|^{3/2}\right),x1.$$
(17)
To make contact with Ref.  we briefly discuss the implications of our results for the ensemble averaged density of states $`\rho (\epsilon )`$ in the sub-gap regime. The tail of $`P(x)`$ for $`x1`$ is the same as the tail of $`\rho `$, cf. Fig. 2. We conclude that

$$\rho (x)\mathrm{exp}\left(\frac{2\beta }{3}x^{3/2}\right)$$
(18)
over a broad range $`\mathrm{\Delta }_\mathrm{g}E_\mathrm{g}\epsilon E_\mathrm{g}`$ inside the mean–field gap. A different exponential decay (with a power $`2`$ instead of $`3/2`$ in the exponent) was predicted recently by Beloborodov, Narozhny, and Aleiner , for the sub-gap density of states of an ensemble of superconducting grains in a weak magnetic field. Since the mean-field density of states in that problem is also of the form (1), the universal GUE edge statistics should apply. The reason that the universal decay (18) was not obtained in Ref.  is that their theory applies to the non-universal energy range $`\epsilon E_\mathrm{g}`$ near the Fermi level. To emphasize the significance of the universal energy range we note that the probability to have the lowest energy level in that range is larger than in the non-universal range by an exponentially large factor $`\mathrm{exp}[(E_\mathrm{g}/\mathrm{\Delta }_\mathrm{g})^{3/2}]`$.

In conclusion, we have argued that the proximity effect in a mesoscopic system has a gap distribution which is universal once energy is measured in units of the energy scale $`\mathrm{\Delta }_\mathrm{g}(E_\mathrm{g}\delta ^2)^{1/3}`$ defined from the mean-field density of states $`\rho (\epsilon )=[(\epsilon E_\mathrm{g})/\mathrm{\Delta }_\mathrm{g}^3]^{1/2}/\pi `$. This universal distribution is the same as the distribution of the smallest eigenvalue of the Gaussian orthogonal or unitary ensembles from random-matrix theory, depending on whether time-reversal symmetry is broken or not. We have identified the magnetic field scale for breaking time-reversal symmetry and verified our results by exact diagonalization of an effective Hamiltonian. Characteristic energy and magnetic field scales are summarized in Table 1. The universality of our prediction should offer ample opportunities for experimental observation.

We thank I. Aleiner, I. Beloborodov, E. Mishchenko, and B. Narozhny for useful discussions. This work was supported by the Cornell Center for Materials Research under NSF grant No. DMR–9632275 and by the Dutch Science Foundation NWO/FOM.


text:
# Acknowledgments

## Acknowledgments

This work was partially financed by Programa de Apoio a Núcleos de Excelência (PRONEX) and CNPq, BRAZIL.

## Figure Captions

Fig. 1: First order contribution to the unitarity corrections. In (a) these corrections are controlled by the proton radius, while in (b) by the constituent radius.

Fig. 2: The $`F_2`$ structure function as a function of the variable $`ln(\frac{1}{x})`$ for different virtualities and radii. Only the unitarity corrections in the quark sector are considered. Data from H1 .

Fig. 3: The $`F_2`$ slope as a function of the variable $`x`$ for different radii. (a) Only the unitarity corrections in the quark sector are considered. (b) The unitarity corrections in the gluon-quark sector are considered. Data from ZEUS . The data points correspond to a different $`x`$ and $`Q^2`$ value.

Fig. 4: Comparison between the DGLAP and Glauber-Mueller (GM) predictions for the behavior of the $`F_2`$ slope using as input in the calculations the GRV(94) or GRV(98) parameterizations. Data from ZEUS . The data points correspond to a different $`x`$ and $`Q^2`$ value.

Fig. 5: Comparison between the Glauber-Mueller (GM) prediction and DGLAP using as input in the calculations the GRV(94) or GRV(98) parameterizations, for the behavior of the $`F_2`$ slope. Preliminary data from H1 .


text:
# Montecarlo simulation of the role of defects as a melting mechanism

## I Introduction

Melting of crystals has always been an exciting subject in condensed matter physics. The Lindemann criterion allows to estimate the melting temperature in a simple way. However, the mechanism responsible for the melting is still debated. It was widely admitted that the melting in three dimensions (3D) occurs when one of the phonon modes is softened by the temperature $`T`$ so that instability of the crystalline phase takes place leading to the liquid phase. In 2D, this scenario is not valid. Mermin  has shown in 1968 that long-range crystalline order is destroyed by $`T`$ if the elastic interaction is power-law-decayed with distance. Nelson and Halperin have shown that the 2D melting is due to defects, analogous to the case of Kosterlitz-Thouless transition for XY spins in 2D . Inspired by this defect-mediated melting, several workers have attempted to prove that in 3D the melting may be also due to dislocations and defects . The soft-mode scenario, though theoretically possible, is not the one frequently observed in numerical calculations. Works on defect theory up to 1989 have been summarized and developed in the book by Kleinert. Recently, several numerical investigations and analytic approximations have shown that dislocations and defects are responsible for melting. In particular, Burakowsky, Prestoni and Silbar  have shown that melting properties of most of elements of the Periodic Table can be explained by excitations of linear dislocations in the crystal using a polymer theory.

Katsnelson and Trefilov  have also recently developed a route toward a melting theory based in defects. They have emphasized the concept of the geometrical frustration , i. e. the 3D euclidean space cannot be filled by the closest packing structure, the tetrahedral one. Therefore they proposed that superdense regions and voids around them are to be formed near melting by thermal activation. These regions were suggested to be the precursor of the destruction of the crystalline order.

In this paper, we investigate by Monte Carlo (MC) simulation the melting of a 3D crystal where atoms interact with each other via the Lennard-Jones (LJ) potential. We show that defects which occur in the solid phase when the transition temperature is approached from below are responsible for the melting. Section II is devoted to a description of our model and MC method. Results are shown in Section III. Concluding remarks are given in Section IV.

## II Model and Monte Carlo Method

We consider a crystal of face-centered cubic (fcc) structure which is described by the following Hamiltonian

$$=\underset{(ij)}{}U(r_{ij})$$
(1)
where the interaction between atoms at $`r_i`$ and $`r_j`$ is described by the potential $`U(r_{ij})`$. For simplicity, the distance dependence is supposed to be given by the LJ potential $`U(r_{ij})=4ϵ[(r_0/r_{ij})^{12}(r_0/r_{ij})^6]`$ where $`r_{ij}=r_ir_j`$, $`r_0`$ being a characteristic length of the system chosen in such a way that the NN distance in the fcc lattice is equal to $`\sqrt{2}/2`$ when only NN interactions are taken into account. The fcc lattice constant is therefore equal to 1 in the ground state. Different potentials for cohesive interactions are envisageable, for example the so-called Gupta many-body potential which has been recently used to study the melting process of a fcc crystal. We have shown that defects created near melting play an important role in the melting mechanism. In the present paper we would like to clarify how defects are topologically distributed and how they destroy the crystal order. To this end, we choose the LJ potential which is simpler for implementing a volume-variable algorithm as will be described below. However, most of the conclusions concerning the melting mechanism are similar for these two potentials. Details of a precise comparison will be given elsewhere.

In this work, we studied a fcc lattice with different system sizes , using periodic boundary conditions, at constant pressure. Interactions up to a cutoff distance $`r_c=1.57`$ have been taken into account. This is about the fourth-nearest-neighbor distance in the perfect fcc crystal.

The following algorithm was used. Starting from the solid state where atoms are on the fcc lattice sites, we heat the system up to a temperature $`T`$. We equilibrate the system first locally at constant volume and then globally at variable volume, as explained hereafter. The local equilibration at constant volume is done as follows: we take an atom and move it to a nearby random position. This position is accepted if it lowers the atom energy. Otherwise it is accepted with a probability according to the Metropolis algorithm. We repeat this for all atoms: we say we achieve one MC step/atom. Next, we change the system volume by a random amount: all atom positions are thus rescaled with the volume variation. We recalculate the energy and accept or reject this new volume using a constant pressure Metropolis algorithm. For a general method, see Ref.. In the following we work at zero pressure. We find that the equilibrium is reached very fast with alternately 10 consecutive local MC steps/atom followed by one volume variation step, and so forth. In all, we performed about $`10^5`$ MC sweepings at each $`T`$. Physical quantities such as averaged internal energy per atom $`U`$ and radial distribution function $`g(r)`$ are averaged over the next $`10^5`$ MC steps/atom.

From the plot of $`E`$ versus $`T`$, we can identify the transition temperature from solid to liquid state. Furthermore, to investigate the melting mechanism we have computed the following quantities: angular distribution of neighbors, distribution of defects, number of nearest-neighbors (NN) etc. As we discuss below, due to the particular fcc structure we cannot use the Voronoi method which would give confusing results if care is not taken.

## III Results

Let us show first in Fig. 1(a) $`U`$ versus $`T`$. One observes a discontinuity of $`U`$ indicating clearly that the transition is of first-order as expected for a 3D melting. This is confirmed by the jump of the average fcc lattice constant shown in Fig. 1(b). These figures show results for three different system sizes. Since the transition is of first order, the transition temperature $`T_m`$ cannot be defined with precision (due to hysteresis). Moreover for sizes greater than 500 particles the melting temperatures fall within the same interval. Therefore, in the following we will show results for the system of N=1372 particles. We take $`T_m`$ as the lower limit of the transition temperature region. From Fig. 1, we take $`T_m=0.76`$. Note that as our simulation cell does not contains surfaces, this transition temperature, corresponds to the metastability limit of the crystal’s ability to superheat rather than to the thermodynamical melting point . In the following, we will call for simplicity, Tm at this transition temperature and associate the phenomena with an apparent melting point since the structure loss all local order at this temperature.

Fig. 2 shows the radial distribution function between atoms at a low $`T`$ ($`T=0.25`$) in the solid phase. The first (second, third,…) peak corresponds to the NN (NNN, third neighbors, …) distance of the fcc structure.

To calculate the coordination number $`c`$ of a site, we integrate the radial distribution function up to the first minimum. For example in Fig. 2, integrating $`g(r)`$ over $`r`$ up to about $`r_1=0.85`$ (first minimum) will give $`c=12`$. Integrating $`g(r)`$ from $`r_1`$ up to the next minimum at $`r_2=1.1`$ gives $`c_2=6`$ and so on. These coordination numbers verify the fcc structure. We use this method instead of the frequently used Voronoi construction because the Voronoi polyhedra in a slightly distorted fcc structure have small faces and short edges and the coordination number could not be precisely calculated because these small faces appear and disappear frequently due to thermal excitations .

Furthermore, to complete our structure determination, we also calculate $`g(\theta )`$, the angular distribution of neighbors around a site within a plane up to a certain distance. The function $`g(\theta )`$ gives the probability to have two neighbors of an atom forming an angle $`\theta `$. We show in Fig. 3 an example of angular distribution taken at $`T=0.25`$ . The peak at 60 (90, 120) degrees corresponds to the angular distribution of the NN (NNN, 3rd NN) atoms in the fcc structure.

Let us increase $`T`$ near the transition. The radial distribution $`g(r)`$ at $`T=0.7`$ is shown in Fig. 4. The main differences with the low $`T`$ structure shown in Fig. 2 are:

* All the minima between the peaks are raised. This is evidence that the crystal order is reduced at $`T=0.7`$. In particular, the minimum between the first and second peaks does not correspond to $`g(r)=0`$. This implies that a diffusional dynamics is activated near melting. This is the origin of the formation of defect clusters we will discuss below.
* Another feature is the strong reduction of the second and fourth peaks. The cubic symmetry is therefore weakened and the remaining correlated part of the system is within the (111) planes.

Integrating $`g(r)`$ over $`r`$ between two consecutive minima, one obtains the coordination number $`c=12`$ only at $`74\%`$ of the total crystal sites. $`15.5\%`$ of the sites have 11 NN, $`7.3\%`$ have 13 NN and $`2\%`$ have 10 NN. At this stage, it is worthwhile to stress that though the temperature is still below $`T_m`$, such an important percentage of sites ($`26\%`$)exhibit defects. The crystal structure cannot be therefore considered as a defect-free one which is used in the search for a phonon soft-mode responsible for the elasticity instability at the melting.

The angular distribution for $`T=0.7`$ is shown in Fig. 5. One observes a strong deviation from the low-$`T`$ one shown in Fig. 3. The cubic symmetry is reduced and the persistence of the (111) structure is seen in this figure: the peak at 90 degrees characteristic of the cubic symmetry is broadened and reduced whereas the other two peaks of the (111) one still have some structure.

For $`T=0.79`$, i.e. a temperature just above the transition, we find the following striking result: $`30\%`$ of sites with 11 NN, $`25\%`$ with 10 NN, $`20\%`$ with 12 NN, $`12\%`$ with 9 NN, $`4\%`$ with 8 NN, $`6\%`$ with 13 NN, $`1\%`$ with 7 NN and with 14 NN. Note that the percentage of 12-NN atoms has dropped from $`74\%`$ at $`T=0.7`$ to $`20\%`$ at this temperature. The average NN number is therefore 10.71.

Let us analyze now the structure of the observed defects. For this purpose we calculated the radial and angular distributions between defects as follows: whenever we are at a defect, i.e. atom with a coordination number different from 12, we search for defects around it up to a certain distance and realize histograms. Figure 6 shows the radial distribution of defects at $`T=0.7`$ around defects with $`c=13`$. One observes that i) the defects are surrounded by other defects at a distance (0.696) slightly shorter than the equilibrium NN distance indicating that at this temperature defects form clusters ii)there is a larger number of defects at a distance 0.876 very close to $`\sqrt{3}/2`$, indicating that defects occupy the sites at the middle of the fcc cube which are normally vacant in the equilibrium configuration. This finding is important since it shows that defects forms clusters with two shells: the inner are defects at 0.696 and the outer at 0.876 from the center defect. At this temperature, this kind of defects are uncorrelated at larger distance (see Fig. 6). The situation for defects with $`c=11`$ is shown in Fig. 7. In this case the second peak is less intensive and order at larger distance appear as we discuss below. Angular distributions between defects at $`T=0.7`$ shows a shift of a few degrees of the peaks toward smaller angles for defects with $`c=13`$ as expected for $`c`$ larger than 12 (see Fig. 8). However for defects with $`c=11`$, only the first peak is lowered, while the other two peaks are shifted a few degrees higher.

Let us analyze more carefully the distribution of the clusters formed by the defects. We have computed the number of defects inside each cluster. As we previously stated we have seen that defects are not isolated but they are arranged in separated inner and outer groups within a defect cluster. The outer groups appear and dissapear while the simulation evolves, making the cluster size vary. To quantify the sizes of these clusters we have obtained a histogram representing the probability to have a cluster with a given number of particles. The results are plotted in Figs. 9 at low, intermediate and near the melting point temperatures. At low temperature \[Fig. 9(a)\], where the density of defects is very low, they appear mainly in pairs. These pairs of defects are completely isolated and not important for the thermodynamical properties. As the temperature increases, larger clusters of defects start to be created \[Fig. 9(b)\]. Note that at temperatures near melting, clusters of all sizes are present with the same probability up to a large given number of particles \[Fig. 9(c)\]. When the cluster size reaches a given critical value the system melts. We can interpret these defect clusters as a set of dislocation arrays or alternatively as liquid zones inside a solid bulk. To enforce this interpretation we show the radial distribution between defects in Fig. 10 for defects with $`c=13`$ at $`T=0.79`$. One observes a large double peak indicating defects at distances between NN and NNN equilibrium distances and pronounced peaks at around 3rd neighbor distance ($`\sqrt{2}`$) and at 2. This suggests that defects are linked over large distances near the transition. Note that this correlation is in fact insinuated below $`T_m`$ as is seen in Fig. 7 for defects of c=11 at T=0.7. The scenario proposed by Burakowsky, Prestoni and Silbar  is somewhat verified. We recall that this theory states that melting appears as a result of dislocation generation. When the density of the dislocation array reach a given critical value, the entropy of these arrays compensates the increase of energy produced by the local breakdown of the perfect crystalline order.

Finally, we emphasize that our results were obtained with simulations using a combination of local and volume updates. However the use of only collective updates with variable volume as in our previous work does not alter our conclusion though the equilibrating time is somewhat longer.

## IV Concluding Remarks

We have studied the melting mechanism of a fcc solid with the LJ potential. The results show that defects created in the solid phase become numerous enough to cause the crystal to melt. At the transition, about one fourth of the total sites do not have the coordination number 12 of the perfect structure. Just above the transition, only about 20$`\%`$ have 12 NN. Our analysis of the structure of defects shows that a defect is surrounded by atoms at $`r`$ about $`\sqrt{2}/2`$ (NN distance) and by atoms at $`r`$ about $`\sqrt{3}/2`$ which is not the NNN distance (=1) of the perfect structure. A closer examination shows that these atoms are themselves defects created by dislocating atoms to somewhere between NN and NNN distances. These dislocated positions correspond to ’bridge’ positions in the potential landscape. We note that statistics taken between defects shown above indicate that defects are linked together at the transition. In other words, the assumption of linear defects by Burakowsky, Prestoni and Silbar is somewhat verified here. By linear defects, one should understand ’strings’ of defects which are not necessarily straight lines of defects. Of course, other aspects of their theory should be further checked. But this is out of the scope of the present work.

Acknowledgements: We are grateful to E. Jagla for useful discussions. L.G. and A.D. thank the University of Cergy-Pontoise for hospitality.

’Laboratoire de Physique Théorique et Modélisation’ is associated with CNRS (ESA 8089).