text:
# Relationship between Electronic and Geometric Structures of the O/Cu(001) System

## I Introduction

Experimental structural studies of Cu(001) with sub-monolayer coverage of oxygen shows clear signs of phase instability of the system (see for example Ref. kittel, ). The recent studies are a continuation of the debate that has circled around two main phases: a $`c(2\times 2)`$ structure on an unreconstructed Cu(001) kono ; holland ; dobler ; sotto  and a $`(2\sqrt{2}\times \sqrt{2})R45^{}`$ structure which involves a reconstructed surface wuttig . Some experiments have indicated the presence of both phases tobin ; lederer . Still other have affirmed the presence of the $`c(2\times 2)`$ overlayer at low coverages, but only in nanometer size domains and not in large well-ordered areas fo96 . Based on these experimental results, it is now widely accepted that for coverage lower than 0.34ML, oxygen forms $`c(2\times 2)`$ islands. If the coverage exceeds this value, the $`(2\sqrt{2}\times \sqrt{2})R45^{}`$ structure is formed kittel . An array of experimental measurements based on surface X-ray diffraction rob , low energy electron diffraction atr90 ; zeng , surface extended X-ray adsorption lederer , and photoelectron diffraction kittel  have provided useful, and at times conflicting, information on $`(2\sqrt{2}\times \sqrt{2})R45^{}`$ superstructure of O on Cu(001). One of the common features of the superstructure, schematically illustrated in Fig. 1, from these studies is the missing of every fourth (010) row of atoms from the top Cu layer. Secondly, all studies point to the oxygen atoms occupying sites located on either side of the missing row between Cu atoms (denoted as Cu2 in Fig. 1). The structure undergoes further reconstruction (some displacement of Cu2 along the (100) direction) and relaxation of mostly the top Cu layer. There is, however, some discrepancy in the exact nature of the relaxation obtained from different experiments. For example, according to Ref. zeng,  the heights $`z_O`$, $`z_{Cu1}`$ and $`z_{Cu2}`$ in Fig. 1, appear in the order $`z_O>z_{Cu2}>z_{Cu1}`$, whereas Ref. kittel,  gives $`z_{Cu1}>z_O>z_{Cu2}`$ (see Table 1 for details).

The competition between the two observed phases of O on Cu(001) and the subsequent reconstruction of the surface with the $`(2\sqrt{2}\times \sqrt{2})R45^{}`$ overlayer naturally beg the question about the origin of this complex behavior. The driving forces for such structural transformation are expected to be related to the chemical binding and electronic structure of the system. Earlier theoretical studies, based on simple models, have correlated the structural properties of O/Cu(001) system to a strong $`p`$O – $`d`$Cu hybridization nors90  and charge transfer col  in the vicinity of the surface. The rationale of these assumptions is supported by calculations, which indicate the presence of the charge transfer from the metal atoms to oxygen madh ; bagus  and a strong $`p`$O – $`d`$Cu hybridization wiell  for O/Cu(001) $`c(2\times 2)`$ system. Our recent first principle calculations for the $`c(2\times 2)`$ O overlayer on Cu(001) stolb01  also confirm the importance of considerations of hybridization and charge transfer in the system and show that these features are eventually traceable to the long-range Coulomb interaction in the system. Clearly, the electronic structure and chemical binding in this system are complex and an adequate description of the relationship between the electronic and atomic structures of O/Cu(001) is yet to be established for the $`(2\sqrt{2}\times \sqrt{2})R45^{}`$ structure. Our goal here is to carry out such a study to understand the driving force for the observed surface reconstruction. For this purpose we take as benchmarks the two proposed surface configurations from Ref. kittel,  and zeng,  and compare their implications. As we shall see, with minor differences in the electronic structure, both geometries provide the same rationale for the ensuing missing row reconstruction. Comparison of the results with those for the $`c(2\times 2)`$ phase provides further evidence for the relative stability of the $`(2\sqrt{2}\times \sqrt{2})R45^{}`$ O overlayer for 0.5 ML coverage.

The calculations in the present paper have been carried out using the real space multiple scattering technique. We have chosen this particular method because apart from being based on first principles, it allows calculations of the charge transfer in the system with ease and provides a good estimate of the nature of the local bonding. This technique also lends itself nicely to the calculation of local densities of electronic states, their projections on the cubic harmonics, and the valence electronic charge densities. We also calculate the lattice Coulomb energies and the perturbation of the Madelung potential induced by the overlayer.

## II Computational Details

Our calculations are based on the density functional theory within the local density approximation (LDA) KSh  and multiple scattering theory. We use the local self-consistent multiple scattering (LSMS) scheme Stocks95 , which is designed for extended heterogeneous solids. Within LSMS a compound is divided into overlapping clusters called local interaction zones (LIZs) centered around atoms of different local environments. For each LIZ, we solve a system of equations for the $`T`$ scattering matrix in the lattice site – angular momentum representation for the muffin-tin (MT) potential Stocks80 :

$$\underset{n_1L_1}{}\left\{\delta _{L^{}L_1}^{n^{}n_1}t_{n^{}l^{}}^1g_{L^{}L_1}^{n^{}n_1}\right\}T_{L_1L}^{n_1n}=\delta _{LL^{}}^{nn^{}},$$
(1)
where $`n`$ denotes the lattice site number, $`L=l,m`$ are angular momentums, $`t_{nl}`$ are the single-site scattering matrix elements, $`g_{LL_1}^{nn_1}`$ denote the free-electron Green’s function matrix elements, and $`T_{LL^{}}^{nn^{}}`$ are scattering path matrix elements. The solutions are used to determine the LIZ Green’s functions, the local densities $`N_l^n(E)`$ of electronic states, their projections $`N_{lm}^n(E)`$ on the cubic harmonics, and local valence charge densities for the atoms located at the center of each LIZ. These charge densities are used further to solve the Poisson’s equation for the entire extended system and hence to build a new potential for the next iteration of the self-consistent process. An application of this method to materials related to the work here has demonstrated its high efficiency and reliability stolb97 ; stolb99 ; stolb01 .

To take into account surface effects on the MT-potential, we follow the approach developed in Ref. ujf, . In this approach the space in the vicinity of the surface is divided into layers belonging to different atomic planes parallel to the surface and the interstitial charge density $`\rho _i`$ is supposed to be layer-dependent reflecting the asymmetry of the system. The detailed description of the technique is given elsewhere stolb01 . Here we only remind our readers that in our calculations the Madelung potential includes monopole $`M_{i\alpha ,j\beta }^{00}`$, dipole $`M_{i\alpha ,j\beta }^{10}`$ and interstitial $`V_{ij}[\rho _j]`$ termsujf :

$$V_{i\alpha }^{Mad}=\underset{j,\beta }{}\left(q_{j\beta }M_{i\alpha ,j\beta }^{00}+d_{j\beta }M_{i\alpha ,j\beta }^{10}+V_{ij}[\rho _j]\right),$$
(2)
where $`i,j`$ denote the layer number, while $`\alpha ,\beta `$ number the atoms in the 2D unit cell, $`d_{j\beta }`$ is a dipole moment of the MT charge density, and the effective charge

$$q_{j\beta }=Z_{j\beta }q_{j\beta }^{MT}+\frac{4\pi }{3}\rho _jR_{j\beta }^3.$$
(3)
The exchange and correlation parts of the potential are determined within LDA using the technique described in Ref. GunLun, .

In this work the above formalism is applied to the case of $`(2\sqrt{2}\times \sqrt{2})R45^{}`$ O overlayer on Cu(001). Fig. 1 is a schematic representation of the system. The parameters for the two proposed configurations from Ref. kittel,  and zeng, , hereafter referred as ”structure I” and ”structure II”, are summarized in Table 1. The second and lower Cu layers are assumed to occupy bulk-like positions. The presence of oxygen atoms and the missing rows reduce the symmetry of the system and create nonequivalent environments for Cu atoms belonging to the same layer. We take this into account by allowing atoms of the first three layers to be nonequivalent as needed. They are marked by different numbers in Fig. 1. Thus we build LIZ’s around oxygen and nine nonequivalent Cu atoms within the top five layers. We take the sites of missing Cu to be occupied by vacancies represented by MT-spheres with zero core electron density and requiring an additional LIZ to be built around the vacancy. These LIZ’s contained 62 to 91 atoms depending on the local configuration, and the requirement that the calculated characteristics of the bulk system closely match those obtained from other reliable methods in the literature.

## III Results and Discussion

Adsorption of gases on surfaces causes both short- and long-range effects on the electronic structure of the system. The short-range effect originates from changes in the local environment of some surface atoms which lead to modifications of their chemical bonds. The long-range effect induced by the Coulomb interaction between the chemisorbed overlayer an the substrate, on the other hand, may cause differing shifts in the local potentials at atomic sites with differing environments. As a result, charge transfer may be induced and the electronic state hybridizations at these sites may be modified.

Here we examine in details both the short- and long-range effects of the chemisorption of an oxygen overlayer on Cu(001) with the view of understanding the driving forces responsible for the observed missing row reconstructed $`(2\sqrt{2}\times \sqrt{2})R45^{}`$ superstructure. We present first a comparative study of the short-range effect in structures I and II through examination of the local densities of electronic states at the oxygen and top Cu sites and their hybridizations. To gain further insight into the nature of the O-Cu bonding and subtle differences emerging from structures I and II, we resort to an analysis of the projections of these densities of states on cubic harmonics. These short-range effects, however, are related to the long-range ones responsible for shifts in the local potentials at each site. We do this next through a systematic evaluation of the effects of the Coulomb interaction on local potentials. This is followed by an analysis of the charge transfer in the system. Finally, we comment on the relative stability of the $`c(2\times 2)`$ phase, in view of the findings here for the $`(2\sqrt{2}\times \sqrt{2})R45^{}`$ phase.

### III.1 Local densities of the pO- and dCu-electronic states

The calculated densities of the $`p`$-states of the oxygen and $`d`$-states of its nearest neighbors Cu1, Cu2 and Cu3 atoms for structures I and II are plotted in Fig. 2. The densities of the $`d`$Cu4 through $`d`$Cu9 states ($`N_d^{Cu}(E)`$) are shown in Fig. 3. Since the difference between $`N_d^{Cu}(E)`$ of Cu6 and Cu7 is found to be negligible, only one is displayed. Figs. 2 and 3 indicate a dramatic difference in $`N_l^n(E)`$ for atomic sites with different environments. The $`N_d^{Cu}(E)`$ of the Cu6 and Cu8 atoms belonging to the third and fourth layers, respectively, differ slightly from bulk values, while that of Cu9 coincides with that of the bulk. The $`d`$Cu4- and $`d`$Cu5-states, on the other hand, consist of a pronounced narrow peak together with some low energy structure. Finally, the densities of the $`p`$O-states and $`d`$-states of oxygen’s nearest neighbors are found to be significantly split and structured. This splitting results in two substructures (”a” and ”b” in Fig. 2), which have the same energetic position and similar shape for all spectra. Such a behavior clearly indicates a strong hybridization and covalent binding between the $`p`$O-states and the $`d`$-states of the neighboring Cu atoms. A relatively high intensity of the low energy peaks ”a” in $`N_d^{Cu3}(E)`$ suggests that the $`d`$Cu3-states are more involved in the $`pd`$ hybridization than the $`d`$Cu1- and $`d`$Cu2-states. To understand the origin of this difference we should remember that the degree of hybridization of electronic states depends on their spatial and energetic overlap. Since, in the present structures, the O-Cu3 interatomic distance is longer than O-Cu1 and O-Cu2 ones zeng ; kittel , the strong $`p`$O-$`d`$Cu3 hybridization does not result from the larger spatial overlap. The energetic overlap, on the other hand, can be quantitatively characterized within the multiple scattering theory by comparison of the single-site resonance energies ($`E_{res}^n`$) for neighboring sites, which reflect energetics of the local potentials and primarily govern the energy of local subbands wein . The calculated $`E_{res}^n`$ for structures I and II, in Table 2, show substantial variation of $`E_{res}^n`$ from site to site. As we shall see in section C, these variations are a result of the long-range interaction and charge transfer in the system. Since from Table 2, the energetic separation between the $`p`$O and $`d`$Cu3 resonances is smallest, we can conclude that the strong covalent $`p`$O – $`d`$Cu3 binding is caused by a large energetic overlap of these states. Also, as $`E_{res}^{Cu4}`$ and $`E_{res}^{Cu5}`$ are much higher than the others in Table 2, there is only a small coupling of the $`d`$Cu4- and $`d`$Cu5-states with those of their neighbors as indicated by the low intensity of the structures seen in the low energy region in Fig. 3.

Having established the strong coupling between the $`p`$O and $`d`$Cu3 electronic states, we proceed next with an analysis of the projections of the local densities of states on cubic harmonics ($`N_{lm}^n(E)`$) which provide more detailed description of covalent binding in the system and help isolate the reasons for the a noticeable difference between the shape and splitting of for structures I and II.

### III.2 Projections of the densities of states on the cubic harmonics

Since in the geometric structures under consideration, the local environment of the oxygen atoms have a reduced symmetry (see Fig. 4) because of the missing rows, the local electronic states are expected to be anisotropic. Therefore, to build a consistent picture of chemical binding in the system, it is important to analyze the projections of electronic state $`N_{lm}^n(E)`$ on the cubic harmonics. Plots in Fig. 5 of such projections of $`p`$O- and the five $`d`$Cu3-states, calculated for structure I, indicate that they differ dramatically from each other. The O-$`p_z`$ and Cu3-$`d_{3z^21}`$ states are found to be strongly split. They form a common subband represented by peaks ”a” and ”b”, which are separated by about 0.5Ry. This strong splitting together with the involvement of major part of these states in hybridization (redistributed into ”a” and ”b” regions) are a clear signature of strong covalent binding. The O-$`p_x`$ and Cu3-$`d_{xz}`$ as well as O-$`p_y`$ and Cu3-$`d_{yz}`$ states also form common split subbands, but not as strongly as the O-$`p_z`$ – Cu3-$`d_{3z^21}`$ states. The cause of this difference in the hybridization becomes clear through a comparison of the angular distribution of corresponding cubic harmonics, which shows that the O-$`p_x`$ – Cu3-$`d_{zx}`$ and O-$`p_y`$ – Cu3-$`d_{zy}`$ spatial overlap is lower than that of O-$`p_z`$ – Cu3-$`d_{3z^21}`$ (Fig. 5). The Cu3-$`d_{xy}`$ and Cu3-$`d_{x^2y^2}`$ states cannot be coupled with the oxygen electronic states by symmetry. Instead they hybridize weakly with the electronic states of Cu4 and Cu5, as reflected by low intensity structure at high energies in Fig. 5. The weakness of this particular hybridization can be traced to the high energetic separation of the single-site resonances (see Table 2).

In Fig. 6, the corresponding splitting of the O-$`p_z`$ – Cu3-$`d_{3z^21}`$ subband for structure II is found to be larger (5.7 Ry versus 5.0 Ry) than that for structure I, implying stronger covalent binding. Although for the O-$`p_x`$ – Cu3-$`d_{zx}`$ and O-$`p_y`$ – Cu3-$`d_{zy}`$ subbands it is harder to quantify the splitting, Figs. 5 and 6 indicate that these subbands are also split stronger for structure II than for structure I. Again, these differences in splitting can be traced to the energetic separation between the $`p`$O and $`d`$Cu3 single-site resonances found for the two structures (see Table 2).

The projections $`N_{dm}^{Cu}(E)`$ calculated for the other nearest neighbors of oxygen, namely, Cu1 and Cu2 are shown in Fig. 7 and 8, respectively. Only a minor part of the $`d`$Cu1- and $`d`$Cu2-electronic states are found to contribute to the peaks responsible for the $`p`$O – $`d`$Cu hybridization. As mentioned above, the single-site $`d`$Cu1- and $`d`$Cu2-resonances are energetically separated from the $`p`$O-resonance that reduces the hybridization. The low energetic structures of $`N_{dm}^{Cu1}(E)`$ and $`N_{dm}^{Cu2}(E)`$ align with ”a” peaks of $`N_{pm}^O(E)`$, which are formed with the $`p`$O-$`d`$Cu3 hybridization. This suggests that the $`d`$Cu1- and $`d`$Cu2-states just admix with the $`p`$O-$`d`$Cu3-subband and the O-Cu1 and O-Cu2 covalent bonds are weaker than the O-Cu3 bond.

### III.3 Long-range Coulomb interaction and charge transfer

For a deeper understanding of the factors controlling the electronic and atomic structures of the O/Cu(001) system, we turn now to considerations of long-range effect of the Coulomb interaction as expressed by the Madelung potential and the subsequent charge transfer in the surface region. To get a simple and intuitive picture we first calculate the Madelung potential for several surface structures of Cu(001) with 0.5 ML of O, within a simple ”frozen” charge density approximation using Eq. 2. Namely, we calculate the effective charges $`q_{i\alpha }`$ introduced in Eq. 3 using bulk-like electronic density for Cu and atom-like density for O. Such an approach gives us a pure geometric effect of the long-range interaction. As expected, the presence of the surface itself causes only small perturbation in the Madelung potential. This is represented by the low values in the column under ”clean” (Cu(001) in Table 3. The admission of 0.5 ML oxygen on bulk terminated Cu(001) causes large perturbation in the Madelung potential at atomic sites in the top few layers, as seen in the column denoted by BT in Table 3. By removing every fourth row of Cu atom from this bulk terminated O/Cu(001) system, we find the perturbation of the Madelung potential to be significantly reduced in the top Cu layer and slightly increased in some of the Cu sites below (see column labeled MR in Table 3). The structures BT and MR are naturally artificial not only in terms of the positions of the top layer Cu atoms, but also that of the O which is taken to occupy the hollow sites in the surface plane. They have been chosen to provide benchmarks for the relative stability of the two observed O superstructures. The calculated perturbation of the Madelung potentials for structure II, i.e. $`(2\sqrt{2}\times \sqrt{2})R45^{}`$ O phase with Cu atoms in relaxed and reconstructed positions, in Table 3, indeed display significant reductions from the values of BT and MR configurations, for all atomic sites. Finally, in Table 3 we have also included the calculated quantities for the $`c(2\times 2)`$ structure (parameters taken from Ref. kittel,  for a coverage less than 0.34 ML) for comparison. Clearly, the $`c(2\times 2)`$ phase, albeit for 0.5 ML coverage, invokes a much stronger perturbation of the Madelung potential than the one in structure II.

An electronic response to the strong perturbation of the Madelung potential is expected to be charge transfer in the vicinity of the surface. The resulting charge distribution would then determine the degree of ionicity of the system, and would influence single-site resonance energies, and consequently, the formation of covalent bonds. Charge transfer is thus a very important characteristic of systems. To characterize it quantitatively, the system has to be divided into spaces belonging to each atom and the charge within this space calculated. In the framework of the MT-approximation, non-overlapping spheres centered on the atomic sites can naturally represent these spaces. A complexity of the present system is that interatomic bond lengths for a number of atoms are different because of the differences in their local environment. Thus different volumes belong to different atomic sites. To make the site charges comparable, we have built the spheres with the fixed radius equal to one-half of the shortest interatomic distance (1.807 Å) and calculated the local valence charges by integrating the self-consistent valence electronic density over these spheres. The results in Table 4 show a significant electronic charge transfer from Cu to O in both structures I and II. The charge transfer pattern is, however, interesting as Cu2, Cu4 and Cu5 donate part of their charge but Cu1 and Cu3 do not. These latter two nearest neighbors of O, in fact, gain a small amount. Such a complex character of the charge transfer reflects the local potential perturbations caused by the missing Cu rows induced by the oxygen overlayer. Indeed, from Table 3 and Table 4 it follows that sites loose electronic charge if the initial Madelung potential is higher and gain it if the potential is lower. Table 4 also reflects subtle differences in the charge transfer to the different sites in the two structures (I and II) considered here. Recall that covalent bonds have also been found stronger in structure II as compared to structure I.

Naturally the Coulomb interaction of the different sublattices plays a very important role in the above considerations. To get further insights into the difference arising from the two structural models chosen in this paper, we compare the Coulomb lattice energies $`E_{lat}`$ for these O/Cu(001) phases. Within the MT approximation schmidt  $`E_{lat}`$ can be expressed as follows:

$$E_{lat}=\frac{1}{2}\underset{i,\alpha }{}\underset{j,\beta }{}M_{i\alpha ,j\beta }q_{i\alpha }q_{j\beta }+\underset{i,\alpha }{}\frac{1}{r_{MT}^{i\alpha }}(1.2(\omega _{i\alpha }\rho _i)^2+3q_{i\alpha }\rho _i\omega _{i\alpha }),$$
(4)
where $`r_{MT}^{i\alpha }`$ and $`\omega _{i\alpha }`$ are the MT-radii and MT-sphere volumes, respectively. The first sum of Eq. 4 represents the Madelung energy and the second comes from the charge neutrality condition. The $`E_{lat}`$ thus includes the Coulomb interaction between the MT charges and that between them and the layer-dependent uniform interstitial charge density in the system. For comparative purposes it is sufficient to include the contribution of the oxygen overlayer and only the top five Cu layers to $`E_{lat}`$. Eq. 4 then yields values of -13.864Ry and -15.420Ry for $`E_{lat}`$, per 2D unit cell, for structure I and structure II, respectively. The lower value of $`E_{lat}`$ along with the finding of stronger O-Cu covalent binding in structure II than in structure I, indicate preference for the former structure.

We have also calculated the Coulomb lattice energy for the $`c(2\times 2)`$ O overlayer on Cu(001) using our self-consistent calculation results from Ref. stolb01,  and obtained the value of $`E_{lat}=`$-8.582Ry. This is much higher than that in the case of both structures I and II. This result is not surprising as we have already shown in Table 3 that the perturbation to the local Madelung potential at the atomic sites is also larger for the $`c(2\times 2)`$ superstructure. Based on these findings we conclude that the long-range Coulomb interaction is the main factor which controls the phase formation in the O/Cu(001) system and makes the $`(2\sqrt{2}\times \sqrt{2})R45^{}`$ phase preferable at oxygen coverage of 0.5ML in agreement with experimental findings.

## IV Conclusions

The results presented above, provide the following microscopic picture for the formation of the $`(2\sqrt{2}\times \sqrt{2})R45^{}`$ O overlayer on Cu(001). Oxygen adsorption in the four-fold hollow sites of Cu(001) induces a strong Coulomb perturbation on the top layer Cu sites. Missing row reconstruction of the Cu surface, lattice relaxation, and charge transfer help significantly reduce this perturbation. In response to the Coulomb perturbation, the electronic charge transfer from the host surface to O is found to have a complex character, involving two of the four nearest neighbors and two next nearest neighbors of O. The perturbation is screened partially and the self-consistent local potentials are such that the energetic separation between the single site $`p`$O- and $`d`$Cu-resonances is different for different Cu sites. This energetic separation essentially determines the degree of the $`p`$O-$`d`$Cu hybridization and the strengths of the covalent bonds formed between oxygen and its nearest neighbors. Thus, the Cu3 atom located right under O is involved in strong covalent binding with O, whereas weak covalent bonds are formed between the O and its other nearest neighbors. The O-Cu3 binding is found to be stronger in structure II than in structure I. Furthermore, as the Coulomb lattice energy is lower for the structure II than for structure I. In extending these consideration to the observed $`c(2\times 2)`$ O overlayer on Cu(001) (for coverage less than 0.34ML), we conclude that the $`c(2\times 2)`$ structure is unstable, for half monolayer coverage, because it induces a strong perturbation at the surface that increases the Coulomb lattice energy of the system. Thus, according to our results, the long-range Coulomb interaction is the main driving force for the O/Cu(001) surface atomic structure formation.

Finally, we note that adsorption processes on metal surfaces are extensively studied by means of such methods as Monte Carlo or molecular dynamics simulations. Many of them are based on the model potentials, which depend only on the local atomic environment, namely on the interatomic distance and the number of the nearest neighbors. Our findings indicate that an application of such approaches to O/Cu(001) can hardly be successful. Indeed, according to our results, the long-range interaction in the system governs not only charge transfer (as usual), but also such ”short-range” effects as the electronic state hybridization and covalent bond formation. Such an implicit impact cannot be taken into account even if a model potential somehow includes the long-range interaction along with a traditional short-range part. This suggests that a proper model potential for O/Cu(001) can be built either by parameterization of the energetic surface obtained from first principle calculations or by development a sophisticated form, which takes into account the interplay between the long- and short-range parts.

###### Acknowledgements.

This work was supported by the National Science Foundation (Grant No CHE9812397). The calculations were performed on the Origin 2000 supercomputer at the National Center for Supercomputing Applications, University of Illinois at Urbana-Campaign under Grant No DMR010001N. The work of TSR was also facilitated by the award of an Alexander von Humboldt Forschungspreis.

Figure Captions

Fig. 1. Schematic illustration of the O/Cu(001) $`(2\sqrt{2}\times \sqrt{2})R45^{}`$ structure.

Fig. 2. Densities of the $`p`$-electronic states of oxygen and $`d`$\- electronic states of its nearest neighbors calculated for structure I (solid line) and structure II (dashed line).

Fig. 3. Densities of the $`d`$-electronic states of the Cu4 through Cu9 atoms calculated for structure I (solid line) and structure II (dashed line).

Fig. 4. Local environment of oxygen atom in the $`(2\sqrt{2}\times \sqrt{2})R45^{}`$ O/Cu(001) structure.

Fig. 5. Projections of the $`p`$O- and $`d`$Cu3-electronic states calculated for structure I.

Fig. 6. Projections of the $`p`$O- and $`d`$Cu3-electronic states calculated for structure II.

Fig. 7. Projections of the $`d`$Cu1-electronic states calculated for structure I (solid line) and the structure II (dashed line).

Fig. 8. Projections of the $`d`$Cu2-electronic states calculated for structure I (solid line) and structure II (dashed line).


format: text

text:
# Megachannel Extraterrestrial Assay Candidates: No Transmissions from Intrinsically Steady Sources

## 1 Introduction

Various searches for extraterrestrial intelligence (SETI) have found signals having all of the expected characteristics of hypothesized ET transmitters, except one. The candidate signals are narrowband, often are in or near a “magic frequency” in a “special” reference frame (e.g., the Galactic barycenter), do not match the characteristics of known, interfering signals of terrestrial or solar system origin, and yet, do not repeat when follow-up observations of the relevant sky positions are performed. Among the relevant programs are the Megachannel Extraterrestrial Assay (META, HS93, 1993, hereinafter HS93), which found several dozen such signals at 1420 and 2840 MHz, and the old Ohio State SETI program (Dixon, 1985), which found the “WOW!” signal at 1420 MHz. (The SERENDIP programs \[Bowyer et al. 1995\] define candidate signals as those that have been detected multiple times at the same sky position, but only modest efforts have been made to analyze these candidates further and it is not clear yet if these candidates can be eliminated as intrinsically steady sources in the way that the META candidates are below.)

By contrast with such an idealized type of signal, it is easy to imagine that received signals might depart significantly from this idealization, particularly with regard to duration and amplitude stability. Cordes, Lazio, & Sagan (1997, hereinafter CLS97) listed (possibly incompletely) reasons for a received signal to appear intermittent:

1. Noise in the receiver electronics, e.g., thermal noise, cosmic ray induced events, and hardware failures;
2. Radio frequency interference (RFI) whose origin is terrestrial, from Earth orbit, or from interplanetary spacecraft;
3. Natural, extrinsic modulation of a (unknown) class of narrowband astrophysical sources such as that caused by interstellar (radio) scintillation (ISS) or gravitational lensing;
4. Natural, extrinsic modulation of ETI sources; and
5. Intrinsic intermittency at the source of an ETI signal due to natural causes, such as planetary rotation or the nature of the transmission (e.g., planetary radar) or, for deliberate reasons, to frustrate detection and decryption by non-target civilizations.

They then proceeded to analyze the impact of Cause 4 upon SETI surveys, with particular emphasis on ISS, motivated by the following considerations. First, it seems to be the simplest mechanism for producing transient signals from otherwise intrinsically steady signals that is amenable to testing. Second, observations of pulsars and active galactic nuclei have demonstrated that ISS is important at centimeter wavelengths (Rickett, 1990) that are commonly used in searches for ETI (Oliver, 1973).

CLS97 were able to show that ISS *by itself* is sufficient to explain the lack of confirmation of any ETI candidate signals in surveys conducted to date. They also predicted the sensitivities required to rule out the scintillation modulation hypothesis.

This paper reports new observations of nine of the eleven META candidates, observations with the requisite sensitivity to allow ISS to be ruled out. In §2 we summarize the extrastatistical candidates from the META and what the ISS hypothesis requires about their signal strengths, in §3 we describe the additional observations we have conducted. We demonstrate that ISS cannot be responsible for modulating an intrinsically steady signal to produce the META candidates in §4 and that gravitational microlensing cannot do the same in §5. We summarize our conclusions in §6.

## 2 The Extrastatistical META Candidates

In this section we summarize briefly the relevant properties of the META extrastatistical signals. For full details of META and the various signals detected, see HS93.

META conducted five surveys of the northern sky, three at the frequency 1420 MHz and two at 2840 MHz, between 1986 and 1991. At each sky position observed, two polarizations were searched in each of three reference frames—the Local Standard of Rest, the Galactic barycenter, and the cosmic microwave background. The vast majority of the signals analyzed were consistent with an exponential distribution as expected for noise from a Fourier transform spectrometer.

Of the roughly $`6\times 10^{13}`$ observations analyzed during the five-year META, only eleven could not be explained as due to either noise or processor failure, e.g., due to cosmic-ray hits on the electronics. It was these eleven extrastatistical detections, four at 1420 MHz and seven at 2840 MHz, that HS93 identified as being candidates for detections of ETI transmitters.

All eleven candidates were narrowband, with the broadest being only four spectrometer channels wide; each channel was 0.05 Hz. Some candidates were unresolved. During the survey, the processing software could halt the survey observations and acquire reobservations on “interesting” signals within 40 s of detection. (The eleven candidates were identified only after the conclusion of META so far more than just the eleven candidates were reobserved in this manner.) In no case was a candidate reobserved during the immediate reobservations. Furthermore, various followup observations were attempted over a five-year period, and no candidate was detected again. (Intense scrutiny of the WOW sky position has also failed to redetect a signal, Gray 1995; Gray & Marvel 2001)

With the exception of the lack of reobservation, the candidates had all of the characteristics expected of ETI transmitters: Narrowband signals in a celestial reference frame. The eleven extrastatistical candidates also all had low Galactic latitudes, consistent with that expected for a Galactic population. HS93 also considered the possibility that these candidates represent signals from a previously unrecognized class of natural sources. Throughout the rest of this paper we shall discuss the candidates as if they are ETI transmitters, but our comments will be equally applicable if they represent a natural class of sources.

## 3 Observations

We observed nine of the eleven META candidates during the interval 1997–1998 as part of Project Phoenix (Cullers, 2000). Table 1 summarizes the observing log. We also report the candidate positions in a common and conventional equinox (J2000).

For a complete discussion of the observational methodology of Project Phoenix, see Cullers (2000). We shall summarize only the relevant details. Coordinated observations were conducted with the 43 m (140 ft) NRAO telescope at Green Bank, WV, and a 30 m telescope in Woodbury, GA. Initial observations were conducted with the NRAO telescope employing a filter to excise terrestrial interference. The detection of a possible signal exceeding a specified threshold (see below) was then transmitted to the telescope in Woodbury. Employing a matched filter and a lower threshold, the telescope at Woodbury re-observed the sky position. A signal had to be detected at both telescopes in order to be considered a genuine celestial signal.

The original META observations used a telescope with a 0$`\stackrel{}{\mathrm{.}}`$5 beamwidth (FWHM) at a frequency of 1420 MHz. In order to cover this area fully with our smaller beam, 0$`\stackrel{}{\mathrm{.}}`$3 at 1420 MHz, a grid of positions centered on the nominal position of the META candidate was observed. The grid included the nominal position of the candidate and pointings separated by 0$`\stackrel{}{\mathrm{.}}`$3 north, east, south, and west of the nominal position.

At each pointing position an integration time of 300 s was used. The bandwidth of the individual spectral channels was 0.68 Hz (vs. 0.05 Hz in META). The equivalent ($`1\sigma `$) sensitivity was $`1.4\times 10^{27}`$ W m-2, and the initial detection threshold for the NRAO telescope was $`7\sigma `$. For comparison the typical META sensitivity ($`1\sigma `$) was $`5.7\times 10^{25}`$ W m-2, and the threshold for the extrastatistical candidates described in §2 was $`31.7\sigma `$.

CLS97 normalized intensities to the mean noise level $`N`$ in the META spectrometer. We shall continue to do so as well. The choice of whether to use the META or Project Phoenix noise level is arbitary, but because the intrinsic signal strength of the (putative) transmitters is also normalized by $`N`$, we must be consistent. Denoting the normalized intensities by $`\eta `$ (as CLS97 did), the various observational thresholds described above become $`\eta =I/N=31.7`$ for META, $`\eta =0.018`$ for the NRAO telescope, and $`\eta =0.0045`$ for the telescope at Woodbury. An important quantity for our analysis will be the dynamic range between the initial candidate detection level $`\eta _1`$ and the subsequent reobservation threshold $`\eta _T`$. For our observations, $`\eta _1/\eta _T10^3`$. For reference, $`\eta _1/\eta _T2`$ for the original META reobservations.

## 4 Interstellar Scintillation

In this section we first summarize key results from CLS97. We then use the formalism developed by CLS97 to show that our new observations exclude the META candidates from being intrinsically steady ETI transmitters.

Compact radio sources observed at centimeter wavelengths (e.g., pulsars) scintillate. Sufficiently distant radio sources ($`0.5`$ kpc) are in the saturated scintillation regime.111We emphasize that like CLS97 we assume that saturated scintillations, in which the rms intensity modulation is 100%, apply. This is unlikely to be the case for the rest of Project Phoenix observations, as most of its observations target nearby stars, so the distance to any potential transmitters is not sufficient for the saturated scintillation regime to obtain. More distant sources in the beam may be in the saturated scintillation regime, though. The signal from an intrinsically steady source of flux density $`S`$ will have an observed flux density of $`gS`$ (in the absence of noise). The probability density function of the ISS gain $`g`$ is $`p(g)=e^gU(g)`$, where $`U(g)`$ is the Heaviside step function. A point discussed at length in Cordes & Lazio (1991) and CLS97 is the ancipital nature of this probability density function: Scintillations can act to render an otherwise undetectable signal detectable, but because the most probable gain is $`g=0`$, ISS will more likely render undetectable an otherwise detectable signal.

A key result of CLS97 was the role of noise in META. They showed that the META candidates could be explained as rare combinations of a gain $`g>1`$ *and* a noise spike. The ISS gain has a decorrelation time scale of minutes to days, depending upon a source’s Galactic longitude and latitude and the velocities of the source, intervening medium, and observer. The noise in the spectrometer decorrelated on the time scales required to compute a fast Fourier transform ($`20`$ s). The failure of the immediate reobservations to detect the source could be understood as due to a reobservation threshold that was too high coupled with the noise decorrelating; later reobservations had the additional complication that the ISS gain had also decorrelated either partially or fully.

CLS97 quantified the probability of redetecting a source with the conditional probability $`P_{2\mathrm{d}}(I_T|I_1;S,\rho )`$ (CLS97, Appendix B). Here a source of intrinsic strength $`S`$ is detected initially at a level $`I_1`$. At a later time, for which the ISS gain correlation coefficient is $`\rho `$ (with $`0\rho 1`$), the source is then reobserved with a reobservation detection threshold of $`I_T`$. Following CLS97 the quantities can be normalized by the noise level at the time of initial detection, $`P_{2\mathrm{d}}(\eta _T|\eta _1;\zeta ,\rho )`$, where $`\zeta S/N`$. However, the expression given by CLS97 (their equation \[B19\]) is incomplete. The correct formulation of $`P_{2\mathrm{d}}(I_T|I_1;S,\rho )`$ is

$$P_{2\mathrm{d}}(I_T|I_1;S,\rho )=\frac{_{I_T}^{\mathrm{}}_{I_1}^{\mathrm{}}𝑑I_2^{}𝑑I_1^{}f_{2I}(I_1^{},I_2^{};S,\rho )}{P_{d,\mathrm{scint}}(I_1;S)}.$$
(1)
The joint intensity probability density function is $`f_{2\eta }(\eta _1,\eta _2;\zeta ,\rho )`$ and is given by equation (B13) of CLS97. The signal detection probability of a scintillating source, above a threshold $`\eta _1`$, is $`P_{d,\mathrm{scint}}(\eta _1,\zeta )`$ and is given by equation (B10) of CLS97.

For completeness Figure 1 shows the corrected $`P_{2\mathrm{d}}(\eta _T|\eta _1;\zeta ,\rho )`$ for the values of $`\eta _T`$ and $`\eta _1`$ relevant to the META observations and with various values of $`\zeta `$. Figure 1 reproduces one of the panels of Figure 3 in CLS97, but with the corrected expression for $`P_{2\mathrm{d}}(\eta _T|\eta _1;\zeta ,\rho )`$, equation (1). Comparison of Figure 1 and Figure 3 in CLS97 shows that the figures in CLS97 that display $`P_{2\mathrm{d}}(\eta _T|\eta _1;\zeta ,\rho )`$ are qualitatively correct and largely quantitatively correct.

Scintillations impact our observations in two ways. First, the original motivation in conducting the observations was to use detection thresholds low enough to rule out the scintillation hypothesis for the META candidates. Second, however, our use of the two telescopes means that a scintillating signal could be undetected in the reobservations with the telescope at Woodbury.

We consider first our coincidence detection scheme, because if this scheme is not robust against scintillations, we can draw no conclusions about the comparisions with the META observations. In the coincidence detection scheme the initial observation threshold is that of the NRAO telescope, $`\eta _1^{}=0.018`$, and the reobservation threshold is that of the telescope at Woodbury, $`\eta _T^{}=0.0045`$. The degree of correlation of the scintillations is unknown, but we consider two possibilities to illustrate the possible range of second detection probabilities.

In the case $`\rho =1`$, there exists no analytic expression for $`P_{2d}(\eta _T|\eta _1;\zeta ,\rho =1)`$, but we can develop a useful approximation. The bivariate intensity density function is

$$f_{2\eta }(\eta _1,\eta _2;\zeta ,\rho =1)=\frac{1}{1+2\zeta }I_0\left(\frac{\zeta \sqrt{\eta _1\eta _2}}{1+2\zeta }\right)\mathrm{exp}\left[\frac{(\eta _1+\eta _2)(1+\zeta )}{1+2\zeta }\right],$$
(2)
where $`I_0(x)`$ is the first order modified Bessel function. We rewrite $`P_{2d}(\eta _T|\eta _1;\zeta ,\rho =1)`$ as

$$P_{2d}(\eta _T|\eta _1;\zeta ,\rho =1)=\frac{1_0^{\eta _T}𝑑\eta _2_0^{\eta _1}𝑑\eta _1f_{2\eta }(\eta _1,\eta _2;\zeta ,\rho =1)}{P_{d,\mathrm{scint}}(\eta _1;\zeta )}$$
(3)
and then, because $`\eta _T^{}`$, $`\eta _1^{}1`$, we expand in powers of $`\eta _T^{}`$ and $`\eta _1^{}`$. We find

$$P_{2d}(\eta _T|\eta _1;\zeta ,\rho =1)\left\{1\frac{\eta _T\eta _1}{1+2\zeta }\left[1\frac{(1+\zeta )(\eta _T+\eta _1)}{2(1+2\zeta )}+\frac{\zeta ^2\eta _T\eta _1}{(1+2\zeta )^2}\right]\right\}e^{\eta _1/(1+\zeta )}.$$
(4)
In order to explain the META candidates as scintillating sources, CLS97 also found that $`\zeta 3`$. For all reasonable values of $`\zeta >0`$, the redetection proabability for the Project Phoenix coincidence detection procedure is $`P_{2d}^{(P)}(\eta _T^{}|\eta _1^{};\zeta ,\rho =1)>0.999`$.

If the scintillations are uncorrelated completely, $`\rho =0`$, then $`P_{2\mathrm{d}}`$ can be evaluated analytically.

$$P_{2\mathrm{d}}(\eta _T|\eta _1;\zeta ,\rho =0)=P_{d,\mathrm{scint}}(\eta _T;\zeta )=e^{\eta _T/(1+\zeta )}$$
(5)
For a source of strength comparable to that preferred by CLS97, $`P_{2\mathrm{d}}^{(P)}(\eta _T^{}|\eta _1^{};\zeta ,\rho =0)0.999`$. A strict *lower* limit to the detection probability is obtained by setting $`\zeta =0`$ (implying that no source is present!) for which $`P_{2\mathrm{d}}^{(P)}0.996`$.

Thus, the coincidence detection scheme is at least 99.6% robust in detecting actual scintillating sources of the source strength estimated by CLS97. Unfortunately, we can place no reasonable limits on $`\rho `$. In order for $`\rho 1`$, the spatial scale length of the scintillations would have to be comparable to or smaller than the baseline between NRAO and Woodbury, $`b500`$ km. At 1.4 GHz ($`\lambda =21`$ cm), such a scintillation scale length is equivalent to a scattering angle $`\theta _d\lambda /b0\stackrel{}{\mathrm{.}}09`$. While large, a scattering angle of this magnitude is not unprecedented. Various low Galactic latitude sources have scattering diameters comparable to or larger than this value (e.g., NGC 6334B (catalog NGC), Moran et al. 1990; Cyg X-3 (catalog ), Molnar et al. 1995). As the META candidates have low Galactic latitudes as well, we can place no constraints on $`\rho `$.

In the roughly 10 yr between the META observations and the observations we report (§3), scintillations would have decorrelated completely. In this case the expression for $`P_{2\mathrm{d}}`$ is given in equation (5). For the META reobservations, the reobservation threshold is that of the NRAO telescope, $`\eta _T=0.018`$. Thus, $`P_{2\mathrm{d}}^{(M)}(\eta _T|\eta _1;\zeta 3,\rho =0)=0.996`$ for a source strength comparable to that estimated by CLS97, and the lower limit is $`P_{2\mathrm{d}}^{(M)}(\eta _T|\eta _1;\zeta =0,\rho =0)=0.982`$.

In order to find the overall probability of detection, we multiply the detection probability from the coincidence scheme with the redetection probability to the original META observations, $`P_{2\mathrm{d}}=P_{2\mathrm{d}}^{(P)}P_{2\mathrm{d}}^{(M)}`$. The strict lower limit to the overall redetection probability ($`\zeta =0`$) is 97.8%, while it is in excess of 99.5% for a source strength comparable to that estimated by CLS97.

These confidence levels are calculated using the nominal META sensitivity. The actual sensitivity varied from $`2.3\times 10^{25}`$ W m-2 to $`4.3\times 10^{24}`$ W m-2. HS93 do not report the appropriate value for each candidate. Using the highest sensitivity reduces the overall redetection probability to $`P_{2\mathrm{d}}(\zeta 3)98.6`$% while using the lowest sensitivity increases the overall redetection probability to $`P_{2\mathrm{d}}(\zeta 3)99.9`$%.

CLS97 also show the number of (uncorrelated) reobservations required to exclude completely the possibility that the META candidates were actual ETI transmitters. In the case that $`\eta _T0`$, the number of reobservations required is 10. Because the maximum number of reobservations for any one candidate is only 3 (at our sensitivity level), we cannot apply this second reality criterion.

## 5 Gravitational Microlensing

As a second extrinsic cause of signal modulation, we consider gravitational microlensing. In this scenario, the continuous transmissions of an ETI transmitter are amplified briefly by a foreground object passing close to the line of sight of the transmitter. This scenario was not covered by CLS97. Gravitational microlensing has an apparent advantage over ISS in that the amplification gain for microlensing $`A`$ is $`A1`$. Bennett & Rhie (1996), Alcock et al. (1996), and Di Stefano (1999) consider the more general problem of detecting terrestrial planets via gravitational lensing.

Gravitational microlensing has been detected in optical monitoring programs of stars toward the Galactic bulge, the Magellanic Clouds, and M 31 (catalog M) (e.g., Udalski et al., 1994; Alcock et al., 1995; Ansari et al., 1996). Because gravitational effects are achromatic, we can make use of the formalism developed to describe these monitoring programs. We shall assume that neither source size nor source blending are important as both can contribute to apparent chromatic effects (Han et al., 2000).

Given the observed signal strengths of the META candidates and the reobservation thresholds in our new observations, the microlensing amplifications required are $`A100`$, somewhat smaller than the ratio $`\eta _1/\eta _T`$. Although we have not conducted extensive tests, a la CLS97, an important conclusion of CLS97 continues to hold: The noise in the META spectrometer would have played a key role in the initial detections, if these are otherwise constant signals. The decorrelation time scale of the noise, approximately 20 s, was a crucial feature of the signal model developed by CLS97. This short decorrelation time explained how the immediate reobservations conducted during META were unable to detect the signal. Like ISS, the decorrelation time of gravitational microlensing, $`\mathrm{\Delta }t_{\mathrm{gl}}>2`$ hr (Paczyński, 1986), is much longer than the time it took for the immediate reobservations to occur. Thus, gravitational microlensing need not account for the full dynamic range $`\eta _1/\eta _T`$. Amplifications this large have been termed extreme gravitational lensing events (EGLEs) or extreme magnification events (EMEs) (Paczyński, 1995; Wang & Turner, 1996; Gould, 1997).

Paczyński (1986) determined the probability for the gravitational microlensing amplification to exceed a fiducial value $`A_0`$, $`P_{\mathrm{gl}}(A>A_0)`$. For $`A_01`$ this function can be approximated as $`P_{\mathrm{gl}}(A>A_0)\tau A_0^2`$, where $`\tau `$ is the optical depth to microlensing. The exact values for $`\tau `$ depend upon the details of the population of lensing objects. Existing estimates for the microlensing optical depth toward the inner Galaxy are sufficient though to place severe constraints on the META candidates. The optical depth toward the inner Galaxy is $`\tau 10^6`$ (Kiraga & Paczyński, 1994; Udalski et al., 1994; Alcock et al., 1995). This estimate depends upon such factors as the mass function and Galactic distribution of the lenses (Gyuk, 1999). In particular, the value can be lower by factors of several, depending upon the Galactic coordinates of the line of sight. For the purposes of this analysis, we ignore these differences. Including these effects would make our conclusions more robust.

In order to explain the META candidates, we require $`A_0100`$ so that $`P_{\mathrm{gl}}(A>A_0)10^{10}`$ toward the inner Galaxy. On average, we expect that in order to obtain an amplification as large as $`A_0`$, the number of *background* objects, i.e., ETI transmitters, must be $`NP_{\mathrm{gl}}1`$. The relevant quantity here is the number of background objects as the number of lenses has already been incorporated through the factor of $`\tau `$ in $`P_{\mathrm{gl}}`$.

We therefore require $`N10^{10}`$ in order to explain the META candidates as being real signals modulated by gravitational microlensing. HS93 estimate that if all of the META candidates represent real signals, the Galactic population of such transmitters is only roughly $`2\times 10^6`$. The observed number of candidates is insufficient, by orders of magnitude, to allow them to be explained as gravitational microlensed transmitters.

We can extend this analysis to future surveys as well. Suppose that we consider lower amplification values. The number of transmitters required within the inner Galaxy must still be $`N\tau ^110^6`$ for gravitational microlensing to be an important factor. This estimate for the gravitational lensing optical depth is for background sources located in the Galactic bulge. We assume that these sources are located within a Galactocentric radius $`R=1`$ kpc. Then the surface density of such transmitters must be $`0.3\mathrm{pc}^2(R/1\mathrm{kpc})^2`$, assuming a disk-like distribution of transmitters. We have assumed a disk-like distribution of transmitters based on the apparent concentration of META candidates to the Galactic plane. Extending this distribution into the Galactic disk, the average distance between the ETI transmitters would be roughly $`0.6\mathrm{pc}(R/1\mathrm{kpc})`$. Given that programs such as the SETI Institute’s Project Phoenix have not detected any transmissions from nearby stars, this surface density is far higher than the actual surface density. We conclude that gravitational microlensing is unlikely to play an important role in modulating the signals from ETI transmitters within the Galaxy, unless there are populations of disk and inner Galaxy ETI transmitters whose densities differ significantly.

## 6 Conclusions

We have reported additional, more sensitive observations of nine of the eleven extrastatistical candidates from the Megachannel Extraterrestrial Assay (Table 1). We have used these observations to evaluate the most simple hypothesis amenable to testing (CLS97): That the extrastatistical candidates were intrinsically steady sources, either extraterrestrial transmitters or an as-yet unknown population of narrowband natural sources, whose intensities were amplified by propagation effects, either interstellar scintillations (§4) or gravitational microlensing (§5).

Using the formalism developed in CLS97 (with a correction for the second detection probability, equation 1), we have shown that these more sensitive limits exclude the scintillation hypothesis: The META candidates do not in general represent a population of scintillating sources. This hypothesis can be excluded at a confidence level of at least 97.8% (for the case of an intrinsically weak source) to a level in excess of 99% (if the source strengths are comparable to that favored by CLS97). Because two of the META candidates were not reobserved with this sensitivity, we cannot rule out the possibility that one or both of them represent ETI transmitters.

We have also shown that the number of transmitters in the inner Galaxy would have to approach $`10^{10}`$ if gravitational microlensing is to be invoked to account for the META candidates. Current constraints on the number of transmitters in the Galaxy are orders of magnitude smaller than this number. Furthermore, existing surveys already suggest that the Galactic population of ETI transmitters is sufficiently small so that gravitational microlensing is unlikely to play a role in future surveys.

HS93 derived various limits on various populations of transmitters. These limits can be revisited in light of our results (see Figure 10 of HS93): There is no more than 1 Kardashev Type II civilization that has constructed an isotropic beacon near the H I line or its second harmonic within the Galaxy or nearest $`10^3`$ galaxies. There are no more than $`10^4`$ Kardashev Type I civilizations broadcasting isotropically within the Galaxy and no more than 1 such civilization with a directed beacon transmitter having 30 dB or more gain anywhere within the Galaxy. More stringent limits on weaker transmissions will have to await future surveys.

We thank J. Cordes for helpful discussions, Jaem for assistance with equation (2), and the referee for suggestions that improved the presentation of our results.


text:
# Heat Conduction and Entropy Production in a One-Dimensional Hard-Particle Gas

## Abstract

We present large scale simulations for a one-dimensional chain of hard-point particles with alternating masses. We correct several claims in the recent literature based on much smaller simulations. Both for boundary conditions with two heat baths at different temperatures at both ends and from heat current autocorrelations in equilibrium we find heat conductivities $`\kappa `$ to diverge with the number $`N`$ of particles. These depended very strongly on the mass ratios, and extrapolation to $`N\mathrm{}`$ resp. $`t\mathrm{}`$ is difficult due to very large finite-size and finite-time corrections. Nevertheless, our data seem compatible with a universal power law $`\kappa N^\alpha `$ with $`\alpha 0.33`$. This suggests a relation to the Kardar-Parisi-Zhang model. We finally show that the hard-point gas with periodic boundary conditions is not chaotic in the usual sense and discuss why the system, when kept out of equilibrium, leads nevertheless to energy dissipation and entropy production.

Low-dimensional systems are special in many ways. Second order phase transitions have anomalous exponents, chemical reactions do not follow the mass action law oz78 , hydrodynamics breaks down due to divergent transport coefficients caused by long time tails aw67 , and electrons in disordered systems are localized a58 . A last item in this list is the divergence of heat conductivity llp01  in $`2`$ spatial dimensions.

For ordered (periodic) harmonic systems it is well known that all transport coefficients are infinite due to the ballistic propagation of modes. Thus one needs either disorder or nonlinear effects in order to have finite conductivity $`\kappa `$. For electric conductivity, disorder in 1-d leads to zero conduction. The main difference between heat and (electronic) charge conduction is that there is no background lattice in the former, i.e. translation invariance is not broken even if the system is disordered. Of course one can study the electronic contribution to heat conduction, but experimentally one never can neglect the ionic contribution foot0 . Thus one has always soft (Goldstone) modes in heat conduction. These soft modes are not localized by disorder mi70 , and they are not affected by nonlinearities. Thus they propagate essentially freely. In high dimensions this has no dramatic consequences. But in low dimensions they become important and lead to the above mentioned divergence. More precisely, one expects a power divergence $`\kappa N^\alpha `$ in $`d=1`$ and a logarithmic divergence in $`d=2`$. Simulations and calculations using the Green-Kubo formula give $`\alpha 0.35`$ to 0.45 llp01 . It is not clear whether the slight discrepancies found between different models are true or artifacts. Theoretically, one would of course prefer a universal value.

There are some exceptions to this general scenario. Apart from models with external potentials and broken translation invariance, the best known 1-d model with finite $`\kappa `$ is the rotor model of glpv00 ; gs00 . Here very fast rotors effectively decouple from their less fast neighbours. Thus very steep jumps of temperature are effectively stabilized and act as barriers to energy transport.

Another model which was recently claimed to have $`\kappa `$ finite ghn01  is the 1-d hard point gas with alternating masses. The same gas with all particles having the same mass is trivial (a collision between particles is indistinguishable from the particles just passing through each other undisturbed), and perturbations just propagate ballistically, leading thus to infinite $`\kappa `$ (i.e., no temperature gradients can build up inside the gas). To break this integrability, it is sufficient to use alternating masses: every even-numbered particle has mass $`m_1`$, and every odd has $`m_2=rm_1`$ with $`r>1`$ ghn01 ; mm83 ; a93 ; acg97 ; h99 ; a01 ; a02 .

The arguments given above for a divergence of $`\kappa `$ with $`N`$ hold also for the hard point gas. It is obviously nonlinear, sound waves dissipate, it is translation invariant, and there is no obvious special feature as in the rotor model. Indeed, prior to ghn01 , heat conduction has been studied in this system by means of simulations in h99 ; a01 . In these papers it was claimed that $`\kappa `$ diverges. But the simulations of h99  had presumably low statistics (very few details were given), while the simulations of a01  are obviously not yet in the scaling regime ($`N`$ is too small) and are compatible with $`\kappa const`$ for $`N\mathrm{}`$. In any case, the exponent $`\alpha `$ suggested by the simulations of a01  is $`0.22`$, much smaller than for all other models.

In view of this confusing situation, and suspecting that the simulations of h99 ; a01 ; ghn01  were not done most efficiently, we decided to make some longer simulations.

We followed a01  in setting $`m_1=1`$ and using Maxwellian heat baths at the ends with $`T_1=2,T_2=8`$ (i.e., after hitting the end, a particle is reflected with a random velocity distributed according to $`P(v)=\mathrm{\Theta }(\pm v)mv/T\mathrm{exp}(mv^2/2T_{1,2})`$. The heat baths sit at $`x=0`$ and $`x=N`$, i.e. the density of the gas is 1. When an even particle with velocity $`v_1`$ collides with an odd one with velocity $`v_2`$, their velocities after the collision are

$$v_1^{}=\frac{(1r)v_1+2rv_2}{(1+r)},v_2^{}=\frac{2v_1+(r1)v_2}{(1+r)}.$$
(1)
Between two collisions, the particles propagate freely. Thus a fast simulation algorithm is event-driven: For each particle $`i`$ we remember its velocity, the time $`t_i`$ of its last collision (initially, $`t_i`$ is put to zero), and the position $`x_i`$ it had at that time. In addition, we maintain a list of future collision times $`\tau _i`$ for each neighboring pair $`(i,i+1)`$ (here the walls are treated formally as particles with $`v_0=v_{N+1}=0,x_0=0,x_{N+1}=N`$). The system is evolved by searching the smallest $`\tau _i`$, updating the triples $`(t_i,x_i,v_i)`$ and $`(t_{i+1},x_{i+1},v_{i+1})`$, and calculating the new future collision times $`\tau _{i1}`$ and $`\tau _{i+1}`$ (the new $`\tau _i`$ is infinite). Since the list $`\{\tau _i\}`$ is essentially a priority queue s90 , we use for it the appropriate data structure of a heap s90 . Using heaps, searching for the next collision takes a CPU time $`O(\mathrm{ln}N)`$. In comparison, a naive search would take $`O(N)`$. This allowed us to make much larger simulations than in previous works. Our largest systems contained 16383 particles and were followed for $`>10^{12}`$ collisions. In spite of this, we had to start with carefully tailored initial configurations to keep transients short. When obtaining statistics one should not forget that measurements should not be made after a fixed number of collisions, but at fixed intervals in real time. The correctness of the results and the absence of transients were checked by verifying that the energy density is constant, as proven in a01 .

Before presenting results, let us discuss the expected dependence on the mass ratio $`r`$. For $`r1`$, equilibration becomes slow (it takes a long time until a fast particle is slowed down to average speed), but perturbations propagate ballistically. Thus a perturbation will be damped out slowly at first, but it will have no long time tails and is damped exponentially. In the other extreme, for $`r\mathrm{}`$ the light particles bounce between pairs of heavy ones which are hardly perturbed. Thus, if a heavy particle is perturbed, we have a situation very similar to the one for $`r1`$. If a light particle is perturbed, its energy is soon given to its two nearest heavy neighbours, which then behave again as for $`r1`$. In contrast, in the intermediate region $`1r\mathrm{}`$ we expect the perturbation to spread non-ballistically for all times. It is in this regime that we expect the fastest convergence to asymptotic behaviour, both with respect to time and to $`N`$.

In Fig.1a we show $`\kappa `$, defined as total energy flux $`J=_im_iv_i^3/2`$ divided by $`\mathrm{\Delta }T`$, versus $`N`$, for four values of $`r`$. The value $`r=1.22`$ is in the small-$`r`$ region and was studied most intensively in a01 . The value $`r=2.62`$ is near the center of the intermediate regime, while $`r=5`$ is clearly above it. Finally, $`r=1.618=(1+\sqrt{5})/2`$ was chosen because it was used in ghn01 , not because of its irrationality (problems with ergodicity related to rational values of $`r`$ exist only for very small $`N`$ a93 ; acg97 ).

A power law would give a straight line with slope $`\alpha `$ in Fig.1a. None of the four curves is really straight. For small $`N`$ the curve for $`r=1.22`$ agrees perfectly with Fig.3 of a01  (which extends only to $`N=1281`$). It shows the strongest curvature (in agreement with the above discussion), and the small-$`N`$ data alone would suggest a cross-over to $`\kappa =const`$. But this curvature stops for large $`N`$ and a closer look shows that the slope increases for $`N>8000`$. The same is true also for the other curves: They all bend down for small $`N`$ but veer up for larger systems (Fig.1b). This is most clearly seen for $`r=1.618`$ and 2.62. It is less clear for $`r=5`$, but the most rational expectation is that also this curve will have the same slope for $`N\mathrm{}`$. Our best estimate $`\alpha =0.32\genfrac{}{}{0pt}{}{+0.03}{0.01}`$ has asymmetric errors because we do not know how much more the curves will bend upwards for very large $`N`$.

The rescaled temperature profiles for $`r=1.22`$ are shown in Fig.2. To verify the claim of a01  that $`T(x)`$ approaches the profile $`T_k(x)`$ predicted by kinetic theory with an inverse power of $`N`$, i.e. $`T(x)T_k(x)N^{0.67}`$, we plot $`T(x)T_k(x)`$ against $`x/N`$. For $`N<2000`$ we see indeed this convergence in perfect agreement with a01 , but not for $`N>2000`$. Instead, it seems that $`T(x)T_k(x)`$ remains different from zero for $`N\mathrm{}`$. The analogous results for $`r=1.618`$ are shown in Fig.3. In that case, the scaling observed in a01  is confined to very small $`N`$, not shown in the figure. The fact that $`T(x)T_k(x)`$ remains finite for $`N\mathrm{}`$ is now obvious. In contrast to a conjecture in ghn01 , the temperature profile also does not become linear for large $`N`$. All results for $`r=1.618`$ are qualitatively also true for $`r=2.62`$ (not shown).

These results are easily understood. For $`r=1.22`$ we are in the small-$`r`$ regime. This explains the slow convergence of $`\alpha `$ with $`N`$ and the weakness of long time tails, manifested in the agreement with kinetic theory. Only at very large $`N`$ we do see the correct asymptotics. For $`r=1.618`$ and 2.62 we are no longer in this regime, the long time tails are stronger, and the disagreement with kinetic theory is more obvious.

In addition to systems driven by thermostats at different temperatures, we also studied systems in equilibrium with periodic boundary conditions. Here the Green-Kubo formula allows $`\kappa `$ to be calculated from an integral over the heat current autocorrelation $`J(t)J(0)`$ llp01 . In Fig.4 we show this, after suitable normalization and after multiplication by a power of $`t`$ which makes it constant for large $`N`$ and $`t`$. We see strong oscillations with periods $`N`$ which reflect the dominance of (damped) sound waves with a fixed velocity of sound (see also a02 ). They were mistaken for statistical fluctuations in ghn01 , showing clearly that the simulations of ghn01  have not reached the asymptotic regime in contrast to what the authors assumed. Our data suggest that $`J(t)J(0)t^{0.66}`$ for large $`N`$ with a cut-off at $`tN`$, giving $`\alpha =0.34`$ in perfect agreement with our previous estimate.

A 1-d hard particle gas should be described macroscopically by hydrodynamics, i.e. by the Burgers resp. Kardar-Parisi-Zhang (KPZ) equation kpz . If we assume that heat diffusion scales like diffusion in KPZ, we might expect $`\alpha =1/3`$ in agreement with our numerics. But we should warn that particle spreading in the 1-d hard particle gas is not superdiffusive (unpublished data; for $`r=1`$ see jepsen ), so the relation with KPZ is not trivial.

Finally, we measured also the propagation of infinitesimal perturbations. Similar simulations were also made in ghn01 , but there the perturbations were added to the ground state ($`E=0`$, all particles are at rest). In contrast, we perturbed equilibrium states, i.e. we performed a standard stability analysis as used e.g. to estimate Lyapunov exponents. Indeed, we only followed the perturbation in velocity space, not in real space. More precisely, after having run the system long enough to have eliminated transients, we chose a tangent vector

$$(\delta v_i(0),\delta x_i(0))=(\delta _{i,0},0)$$
(2)
and iterated, together with the system itself, its linearized variational equations. Notice that $`\delta v_i(t)`$ are independent of the positions perturbations $`\delta x_i(t)`$ for nearly all times (whenever there is no collision), thus it is possible and legitimate to solve the variational equations for velocities only. After the evolution of the $`\delta v_i(t)`$ has been followed for a given time $`t_{\mathrm{max}}`$, $`\delta v_i(0)`$ is again chosen as $`\delta _{i,0}`$, and the integration is continued.

It is easy to prove that that the 1-d hard point gas is not chaotic in the usual sense. Consider the weighted $`L_2`$ norm of the perturbation, $`\delta v(t)_2=[_{i=1}^Nm_i(\delta v_i(t))^2]^{1/2}`$. Since the $`\delta v_i(t)`$ change during a scattering according to the same Eq.(1) as the velocities $`v_i(t)`$ themselves, energy conservation leads to $`\delta v(t)_21`$. Indeed the absence of chaos is quite obvious since there is no local instability. It seems to contradict a widespread folklore that dissipation and entropy production are tightly related to chaos (which sometimes is true; e.g. gd99 , page 231) – although it is also appreciated that this might not be always the case dc00 .

One solution to this puzzle is the observation gs99 , going back to work by Wolfram on cellular automata w84 , that the notion of chaos in systems with infinitely many degrees of freedom is ambiguous and is not necessarily related to any local instability. In a spatially extended system it makes perfect sense to use a norm which, in contrast to the $`L_2`$ norm used above, puts most weight on near-by regions and exponentially little weight on far-away regions. With such a definition, the norm of a perturbation moving towards (away from) the observer with constant velocity will increase (decrease) exponentially. More generally, also perturbations spreading diffusively will lead to an increase of the uncertainty about the local state for a short-sighted observer. For the 1-d hard point gas this means that there is no need for any local instability to generate dissipation, local thermal equilibrium, and mixing. In a non-equilibrium case entropy flow is provided by the stochastic thermostats at the ends, while (coarse-grained) entropy is produced by the diffusive propagation of perturbations.

Indeed, the situation is not quite as simple due to the divergence of $`\kappa `$ with $`N`$ which suggests that perturbations propagate super-diffusively. For a more direct proof, we show in Fig.5 their effective diffusion coefficient, i.e. the average squared distances divided by $`t`$ (notice that $`m_i(\delta v_i(t))^2=1`$),

$$DX^2/t=t^1\underset{i=1}{\overset{N}{}}i^2m_i(\delta v_i(t))^2.$$
(3)
For diffusive spreading, this would be constant. Instead, it increases with $`t`$ for all $`r`$. This increase does not follow a pure power law, but again deviations from a power law are strongest for very small and very large $`r`$. For moderately large $`t`$ they would suggest a crossover to normal diffusion ($`D=const`$), but for all $`r>1.6`$ this is ruled out by the values at very large $`t`$. Again it is hard to estimate the asymptotic behaviour precisely, because the curves bend upward for large $`t`$. If we assume $`Dt^\beta `$, we obtain a precise lower bound $`\beta >0.36\pm 0.01`$, but only a poor upper bound which just excludes ballistic behaviour.

In summary, we have given compelling evidence that heat conduction in the 1-d hard point gas shows the anomalous divergence with system size expected for any generic 1-d system, in spite of strong finite-size and finite-time effects. This “normal” anomalous behaviour holds in spite of the fact that the system is not chaotic in the usual sense, proving again that chaos in the form of local instabilities is not needed for mixing behaviour and dissipation. Finally, we have discussed a possible connection to KPZ scaling.

We are indebted to Roberto Livi and Antonio Politi for very helpful correspondence. W. N. is supported by the DFG, Sonderforschungsbereich 237.


text:
# Spectrum and Duration of Delayed MeV-GeV Emission of Gamma-Ray Bursts in Cosmic Background Radiation Fields

## 1. Introduction

Gamma-ray bursts (GRBs) are the brightest electromagnetic phenomena in the universe. Their radiation features have been well understood even though their origin has been unknown. In the standard fireball shock model (Piran 1999; van Paradijs, Kouveliotou & Wijers 2000; Mészáros 2002), GRBs are explained to be due to the dissipation of kinetic energy of an expanding fireball with an average Lorentz factor more than $`10^2`$ in the internal shocks produced by collisions between different shells in the fireball. Afterglows are also considered to be due to the dissipation of kinetic energy in the external shocks generated by collision of the fireball with its surrounding medium. Such a model has been supported by a variety of observations on GRBs and afterglows (e.g., Wijers, Rees & Mészáros 1997; Waxman 1997; Vietri 1997). Theoretically, it has also predicted high-energy emission from GRBs and afterglows (for a review see Mészáros 2002). The observed spectra of GRBs sometimes extend to GeV-TeV photon energies. For example, the EGRET experiment detected an 18 GeV photon from GRB 940217 (Hurley et al. 1994), and the Milagrito air shower experiment tentatively detected the TeV emission from GRB 970417a (Atkins et al. 2000). These GRBs might have originated at a redshift of $`z<2`$ (Salamon & Stecker 1998; Totani 2000).

The purpose of this paper is to present a new prediction of delayed MeV-GeV emission of GRBs in the cosmic background radiation fields. Such emission arises from inverse-Compton scattering of electron/positron pairs off cosmic microwave background (CMB) photons. In §2 we generally discuss prompt high-energy emission above a few hundreds of GeV due to synchrotron self-Compton scattering in internal shocks based on the standard fireball model. However, such high-energy photons cannot be detected because they may be absorbed in interacting with the cosmic infrared background radiation (Stecker, De Jager & Salamon 1992; Madau & Phinney 1996; MacMinn & Primack 1996; Malkan & Stecker 1998; Salamon & Stecker 1998). In §3 we show that inverse-Compton scattering of the resulting electron/positron pairs off CMB photons will produce delayed MeV-GeV emission by deriving the time-integrated spectrum of such emission and by estimating the emission duration. In §4 we summarize our findings and discuss their implications. Cheng & Cheng (1996) have discussed similar delayed MeV-GeV emission, and here we further calculate the emission spectrum and estimate the emission duration in the standard fireball shock model.

## 2. High-Energy Radiation from Internal Shocks

Let us consider a compact source that produces a wind (i.e., a variable fireball) with an average luminosity of $`L`$ and with an average mass loss rate of $`\dot{M}=L/\eta c^2`$. Initially, the bulk Lorentz factor of the wind, $`\mathrm{\Gamma }`$, increases linearly with radius, until most of the wind energy is converted to kinetic energy and $`\mathrm{\Gamma }`$ saturates at $`\mathrm{\Gamma }\eta `$. After this time, fluctuations of $`\mathrm{\Gamma }`$ due to variabilities in $`L`$ and/or $`\dot{M}`$ on time scale $`\mathrm{\Delta }t`$ could lead to internal shocks in the wind at radius $`R\mathrm{\Gamma }^2c\mathrm{\Delta }t`$. We also assume that the internal shock Lorentz factors $`\mathrm{\Gamma }_i`$ a few in the wind rest frame, so that the internal shocks could reconvert a significant fraction of the kinetic energy to internal energy, which is then emitted as $`\gamma `$-rays by synchrotron and inverse-Compton (IC) radiation of the electrons accelerated by the shocks.

The comoving electron number density of the unshocked wind matter at radius $`R`$ is $`n_e^{}=L/(4\pi \mathrm{\Gamma }^6\mathrm{\Delta }t^2m_pc^5)`$, where $`m_p`$ is the proton mass. The internal energy density of the shocked wind matter is thus approximated by $`e^{}4\mathrm{\Gamma }_i(\mathrm{\Gamma }_i1)n_e^{}m_pc^2`$ (Blandford & McKee 1976). The electron energy distribution behind the shocks is usually a power law: $`dn_e/d\gamma _e\gamma _e^p`$ for $`\gamma _e\gamma _m`$. Assuming that $`ϵ_e`$ and $`ϵ_B`$ are constant fractions of the internal energy density going into the electrons and the magnetic field respectively, we obtain the minimum electron Lorentz factor $`\gamma _m=[(p2)/(p1)](m_p/m_e)ϵ_e\mathrm{\Gamma }_i`$ (where $`m_e`$ is the electron mass) and the magnetic field strength of the shocked wind matter $`B=(8\pi ϵ_Be^{})^{1/2}=6.2\times 10^3ϵ_{B,2}^{1/2}(\mathrm{\Gamma }_i/2)^{1/2}\mathrm{\Gamma }_{500}^3L_{52}^{1/2}\mathrm{\Delta }t_2^1`$ G, where $`ϵ_{B,2}=ϵ_B/10^2`$, $`\mathrm{\Gamma }_{500}=\mathrm{\Gamma }/500`$, $`L_{52}=L/10^{52}\mathrm{erg}\mathrm{s}^1`$, and $`\mathrm{\Delta }t_2=\mathrm{\Delta }t/10^2\mathrm{s}`$. Thus, we estimate the characteristic frequency of synchrotron radiation from the shock-accelerated electrons,

$`\nu _m`$ $`=`$ $`\mathrm{\Gamma }\gamma _m^2{\displaystyle \frac{eB}{2\pi m_ec}}`$ (1)
$`=`$ $`3.0\times 10^{19}[(p2)/(p1)]^2`$
$`\times ϵ_{e,0.5}^2ϵ_{B,2}^{1/2}(\mathrm{\Gamma }_i/2)^{5/2}\mathrm{\Gamma }_{500}^2L_{52}^{1/2}\mathrm{\Delta }t_2^1\mathrm{Hz},`$
where $`ϵ_{e,0.5}=ϵ_e/0.5`$. According to Sari, Piran & Narayan (1998), the cooling Lorentz factor $`\gamma _c\gamma _m`$ for typical parameters of the internal shock model, implying that the shock-accelerated electrons are in the fast cooling regime. Therefore, the synchrotron spectral luminosity is

$$L_\nu \frac{h\nu N_\nu }{\mathrm{\Delta }t}=\{\begin{array}{cc}L_{\nu _m}(\nu /\nu _m)^{1/2}\hfill & \mathrm{if}\nu <\nu _m\hfill \\ L_{\nu _m}(\nu /\nu _m)^{p/2}\hfill & \mathrm{if}\nu \nu _m,\hfill \end{array}$$
(2)
where $`h`$ is Planck’s constant and $`N_\nu `$ is the synchrotron photon number per unit frequency in $`\mathrm{\Delta }t`$.

Although GRBs are optically thin to electron scattering, some synchrotron photons will Compton scatter on the shock-accelerated electrons, producing an additional IC component at higher energies (Papathanassiou & Mészáros 1996; Pilla & Loeb 1998; Panaitescu & Mészáros 2000). The ratio of the IC to synchrotron radiation luminosity can be estimated by

$$Y\frac{L_{\mathrm{IC}}}{L_{\mathrm{syn}}}=\frac{\nu _m^{\mathrm{IC}}L_{\nu _m}^{\mathrm{IC}}}{\nu _mL_{\nu _m}}=\left(\frac{ϵ_e}{ϵ_B}\right)^{1/2},$$
(3)
where we have assumed that $`ϵ_eϵ_B`$ and the radiating electrons are in the fast cooling regime (Panaitescu & Kumar 2000; Sari & Esin 2001; Zhang & Mészáros 2001). Since the characteristic IC frequency

$`\nu _m^{\mathrm{IC}}`$ $``$ $`\gamma _m^2\nu _m`$ (4)
$`=`$ $`1.0\times 10^{26}[(p2)/(p1)]^4`$
$`\times ϵ_{e,0.5}^4ϵ_{B,2}^{1/2}(\mathrm{\Gamma }_i/2)^{9/2}\mathrm{\Gamma }_{500}^2L_{52}^{1/2}\mathrm{\Delta }t_2^1\mathrm{Hz},`$
the ratio of the IC spectral luminosity at $`\nu _m^{\mathrm{IC}}`$ to synchrotron spectral luminosity at $`\nu _m`$ is

$$\frac{L_{\nu _m}^{\mathrm{IC}}}{L_{\nu _m}}=\gamma _m^2Y,$$
(5)
which is derived from equation (3), and the ratio of the IC to synchrotron photon number per unit frequency in $`\mathrm{\Delta }t`$ is further given by

$$\frac{N_\nu ^{\mathrm{IC}}}{N_\nu }=\gamma _m^{p2}Y,\mathrm{for}\nu \nu _m^{\mathrm{IC}}.$$
(6)
The intrinsic cutoff energy in the IC component is defined by $`_{\mathrm{cut}}^{\mathrm{in}}=\mathrm{min}(_{\mathrm{KN}}^{\mathrm{IC}},_0^{\mathrm{in}})`$, where $`_{\mathrm{KN}}^{\mathrm{IC}}`$ is the Klein-Nishina limit and $`_0^{\mathrm{in}}`$ is the energy at which a photon may be attenuated due to pair production through interactions with softer photons (also radiated from the internal shocks) whose frequency is equal to or larger than $`\nu _{\mathrm{an}}=(\mathrm{\Gamma }m_ec^2)^2/(h_0^{\mathrm{in}})`$. We first estimate $`_{\mathrm{KN}}^{\mathrm{IC}}`$. There is the maximum Lorentz factor $`\gamma _M`$ in the electron energy distribution behind the shocks due to the fact that the shock-acceleration time cannot exceed the cooling time. This gives $`\gamma _M=\{3e/[\sigma _TB(1+Y)]\}^{1/2}`$ where $`e`$ is the electron charge and $`\sigma _T`$ is the Thomson cross section. Thus, we obtain the Klein-Nishina limit $`_{\mathrm{KN}}^{\mathrm{IC}}=\gamma _M\mathrm{\Gamma }m_ec^256ϵ_{e,0.5}^{1/4}(\mathrm{\Gamma }_i/2)^{1/4}\mathrm{\Gamma }_{500}^{5/2}L_{52}^{1/4}\mathrm{\Delta }t_2^{1/2}`$TeV. To compute $`_0^{\mathrm{in}}`$, we adopt an analytical approach proposed by Lithwick & Sari (2001). Integrating the synchrotron spectral luminosity over frequency should lead to the total synchrotron luminosity, which further gives the synchrotron photon number per unit frequency at $`\nu _m`$: $`N_{\nu _m}=\{(p2)/[2(p1)]\}ϵ_e(1+Y)^1(h\nu _m^2)^1L\mathrm{\Delta }t`$. The total synchrotron photon number at frequencies larger than $`\nu `$ can be estimated by $`N_{>\nu }_\nu ^{\mathrm{}}N_\nu 𝑑\nu =_\nu ^{\mathrm{}}N_{\nu _m}(\nu /\nu _m)^{(p+2)/2}𝑑\nu =\{(p2)/[p(p1)]\}ϵ_e(1+Y)^1(h\nu _m)^1L\mathrm{\Delta }t(\nu /\nu _m)^{p/2}`$. The attenuation optical depth is

$$\tau _{\gamma \gamma }^{\mathrm{in}}=\frac{(11/180)\sigma _TN_{>\nu _{\mathrm{an}}}}{4\pi (\mathrm{\Gamma }^2c\mathrm{\Delta }t)^2}.$$
(7)
The requirement that $`\tau _{\gamma \gamma }^{\mathrm{in}}=1`$ leads to the upper energy limit

$`_0^{\mathrm{in}}`$ $`=`$ $`160ϵ_{e,0.5}^{(2p3)/p}ϵ_{B,2}^{1/2}(\mathrm{\Gamma }_i/2)^{5(p2)/(2p)}`$ (8)
$`\times \mathrm{\Gamma }_{500}^{4(p+1)/p}L_{52}^{(p+2)/(2p)}\mathrm{\Delta }t_2\mathrm{GeV},`$
where the numerical value is quoted for $`p=2.2`$. It is easy to see that $`_{\mathrm{KN}}^{\mathrm{IC}}_0^{\mathrm{in}}`$ for typical parameters and thus the intrinsic cutoff energy is in fact given by equation (8).

## 3. Pair Production and Reradiation in Cosmic Background Photon Fields

High-energy $`\gamma `$-rays emitted from the internal shocks may not only be intrinsically attenuated due to pair production through interactions with softer photons from the same radiation regions but may also be absorbed in the external background radiation fields when these $`\gamma `$-rays travel towards the observer. In the latter case, the observed cutoff energy, $`_{\mathrm{cut}}^{\mathrm{ob}}`$, is determined by the redshift ($`z`$) and comoving energy, $`\epsilon _\gamma (z)`$, of the background photons that dominate the pair production optical depth ($`\tau _{\gamma \gamma }^{\mathrm{ex}}`$) because the cross section peaks when $`_{\mathrm{cut}}^{\mathrm{ob}}\epsilon _\gamma (z)(m_ec^2)^2/(1+z)`$. It is general that at low redshifts the spectral energy distribution of the starlight background radiation peaks at infrared wavelengths and thus $`_{\mathrm{cut}}^{\mathrm{ob}}`$ is in the TeV energy range (Salamon & Stecker 1998). Detailed predictions of $`_{\mathrm{cut}}^{\mathrm{ob}}`$ differ from model to model, depending on the theoretical treatment adopted for the stellar emissivity. For example, the observed cutoff energy when $`\tau _{\gamma \gamma }^{\mathrm{ex}}=1`$ for $`z=1`$ is $`50\mathrm{GeV}_{\mathrm{cut}}^{\mathrm{ob}}80\mathrm{GeV}`$ for different models used in Salamon & Stecker (1998).

The pair production optical depth, $`\tau _{\gamma \gamma }^{\mathrm{ex}}`$, strongly depends on the $`\gamma `$-ray energy ($`_\gamma ^{\mathrm{ob}}`$). Salamon & Stecker (1998) numerically calculated $`\tau _{\gamma \gamma }^{\mathrm{ex}}`$ as a function of the photon energy for several fixed redshifts by considering the stellar emissivity with and without metallicity correction, shown in their Figures 6 and 7, respectively. It can be seen from their Figure 6 that if $`z=1`$, then $`\tau _{\gamma \gamma }^{\mathrm{ex}}1`$ for $`_\gamma ^{\mathrm{ob}}=50`$GeV, but $`\tau _{\gamma \gamma }^{\mathrm{ex}}10`$ for $`_\gamma ^{\mathrm{ob}}=300`$GeV. In this model of Salamon & Stecker (1998), therefore, a photon with energy of $`_\gamma ^{\mathrm{ob}}300`$GeV must have a pair production optical depth $`\tau _{\gamma \gamma }^{\mathrm{ex}}1`$. This implies that such high-energy $`\gamma `$-rays may be locally attenuated once they are radiated from internal shocks. For simplicity, we will neglect any redshift correction for low-redshift sources discussed in this paper.

The resulting electron/positron pairs have Lorentz factors of $`\gamma _e_\gamma ^{\mathrm{ob}}/(2m_ec^2)\gamma _{\mathrm{min}}_0^{\mathrm{ex}}/(2m_ec^2)=3\times 10^5_{0,300}^{\mathrm{ex}}`$, where $`_0^{\mathrm{ex}}=_{0,300}^{\mathrm{ex}}\times 300\mathrm{GeV}`$ is the minimum energy of photons that are locally attenuated in the external infrared background radiation field. The pairs will Compton scatter the CMB photons. As a result, the initial energy of a microwave photon, $`h\nu _0`$, is boosted by IC scattering up to an average value $`\gamma _e^2h\nu _057(_{0,300}^{\mathrm{ex}})^2`$MeV, where $`h\nu _0=2.7kT`$ is the mean energy of the CMB photons with $`T2.73`$K and $`k`$ is the Boltzmann constant. The IC lifetime (in the local rest frame) of an electron with Lorentz factor of $`\gamma _e`$ reads

$$\tau (\gamma _e)=\frac{3m_ec}{4\gamma _e\sigma _Tu_{\mathrm{cmb}}}2.4\times 10^{14}(_{0,300}^{\mathrm{ex}})^1\mathrm{s},$$
(9)
where $`u_{\mathrm{cmb}}=aT^4`$ is the CMB energy density and $`a`$ is the radiation constant. It is not difficult to find that the typical length in which most of the electron energy is lost, $`c\tau (\gamma _e)`$, is much less than the distance from the source to the observer, implying that energy loss of the electron due to IC scattering is also local. Beyond this length, energy loss of the electron and its emission become insignificant in the absence of any acceleration.

The typical duration of the scattered photons after the GRB trigger in the observer’s frame is estimated by $`\tau ^{\mathrm{ob}}=\mathrm{max}(\tau _1^{\mathrm{ob}},\tau _2^{\mathrm{ob}})`$, where $`\tau _1^{\mathrm{ob}}`$ is the observed IC cooling lifetime,

$$\tau _1^{\mathrm{ob}}\frac{\tau (\gamma _e)}{2\gamma _e^2}1.3\times 10^3(_{0,300}^{\mathrm{ex}})^3\mathrm{s},$$
(10)
and $`\tau _2^{\mathrm{ob}}R_{\mathrm{pair}}/(2\gamma _e^2c)`$ is the angular timescale (Piran 1999). Here $`R_{\mathrm{pair}}=(0.26\sigma _Tn_{\mathrm{IR}})^15.8\times 10^{24}(n_{\mathrm{IR}}/1\mathrm{cm}^3)^1`$ cm is the typical pair-production radius (where $`n_{\mathrm{IR}}1\mathrm{cm}^3`$ is the cosmic infrared photon number density, Protheroe & Stanev 1993). Thus, we have

$$\tau _2^{\mathrm{ob}}1.0\times 10^3(n_{\mathrm{IR}}/1\mathrm{cm}^3)^1(_{0,300}^{\mathrm{ex}})^2\mathrm{s}.$$
(11)
Therefore, most of the scattered photons will reach the observer in $`10^3`$ seconds. This time is much longer than a typical GRB duration.

We next derive the scattered photon spectrum. From equations (2), (3) and (6), we first write the spectrum of the electrons and positrons as follows

$`{\displaystyle \frac{dN_e}{d\gamma _e}}`$ $`=`$ $`2\gamma _m^{p2}Y\nu _mN_{\nu _m}\left({\displaystyle \frac{h\nu _m}{2m_ec^2}}\right)^{p/2}\gamma _e^{(p+2)/2},`$ (12)
$`\mathrm{if}\gamma _{\mathrm{min}}\gamma _e\gamma _{\mathrm{max}},`$
where $`\gamma _{\mathrm{max}}=_0^{\mathrm{in}}/(2m_ec^2)`$. The scattered photon number per unit frequency is given by (Blumenthal & Gould 1970, hereafter BG70)

$$N_\nu ^{\mathrm{sc}}=\left(\frac{dN_e}{d\gamma _e}\right)\left(\frac{dN_{\gamma _e,ϵ}}{dt^{}d\nu }\right)\tau (\gamma _e)𝑑ϵ𝑑\gamma _e,$$
(13)
where $`dN_{\gamma _e,ϵ}/dt^{}d\nu `$, expressed by equation (2.42) of BG70 in the Thomson limit, is the spectrum of photons scattered by an electron with Lorentz factor of $`\gamma _e`$ from a segment of the CMB photon gas of differential number density $`n(ϵ)dϵ`$, and $`t^{}`$ is the time measured in the local rest frame. Please note that the scattered photon number per unit frequency given by equation (13) is time-integrated by multiplying $`\tau (\gamma _e)`$ in equation (2.61) of BG70. After integrating over $`\gamma _e`$ and $`ϵ`$, we obtain the ratio of the scattered to synchrotron photon number per unit frequency

$`{\displaystyle \frac{N_\nu ^{\mathrm{sc}}}{N_\nu }}`$ $`=`$ $`{\displaystyle \frac{9(m_ec^2)^{(p2)/2}}{2^{(p+8)/2}\pi ^2(\mathrm{}c)^3aT^4}}\gamma _m^{p2}Y`$ (14)
$`\times F(q)(kT)^{(p+14)/4}(h\nu )^{(p2)/4}`$
$`=`$ $`2^{p/2}Y{\displaystyle \frac{F(q)}{F(3)}}\left[{\displaystyle \frac{\gamma _m^4(kT)(h\nu )}{(m_ec^2)^2}}\right]^{(p2)/4},`$
where $`F(q)`$ is a function of $`q(p+4)/2`$,

$$F(q)=\frac{2^{q+3}(q^2+4q+11)}{(q+3)^2(q+1)(q+5)}\mathrm{\Gamma }[(q+5)/2]\zeta [(q+5)/2]$$
(15)
with the Gamma function and Riemann zeta function. Therefore, the scattered photon spectrum becomes

$$N_\nu ^{\mathrm{sc}}\nu ^{(p+6)/4}.$$
(16)
The lower energy limit of the scattered photons is

$$h\nu _1\gamma _{\mathrm{min}}^2(2.7kT)=57(_{0,300}^{\mathrm{ex}})^2\mathrm{MeV},$$
(17)
while the upper energy limit is

$`h\nu _2`$ $``$ $`\gamma _{\mathrm{max}}^2(2.7kT)`$ (18)
$``$ $`16ϵ_{e,0.5}^{2(2p3)/p}ϵ_{B,2}^1(\mathrm{\Gamma }_i/2)^{5(p2)/p}`$
$`\times \mathrm{\Gamma }_{500}^{8(p+1)/p}L_{52}^{(p+2)/p}\mathrm{\Delta }t_2^2\mathrm{MeV},`$
which is derived from equation (8).

We now consider a simple example. The typical parameters of the internal shock model are taken: $`ϵ_{e,0.5}=ϵ_{B,2}=L_{52}=\mathrm{\Delta }t_2=\mathrm{\Gamma }_i/2=1`$, $`\mathrm{\Gamma }_{500}=1.4`$ and $`p=2.2`$. Furthermore, $`_{0,300}^{\mathrm{ex}}=1`$ is assumed. Therefore, we obtain $`h\nu _157`$MeV and $`h\nu _20.8`$GeV. From equation (14), we further find $`N_\nu ^{\mathrm{sc}}/N_\nu `$ increases from $`5.0`$ to $`5.8`$ when the observed energy of the scattered photons increases from $`h\nu _1`$ to $`h\nu _2`$. This leads to a delayed MeV-GeV emission component. It is interesting to note that the total energy release of such emission is not only comparable to that of a typical GRB, but also much larger than the total MeV-GeV energy release in an early afterglow as calculated by Dermer, Chiang & Mitman (2000) and Zhang & Mészáros (2001).

## 4. Discussion and Conclusions

We have discussed prompt high-energy emission from the internal shocks, which are produced by collisions between shells with different Lorentz factors in an ultrarelativistic wind ejected from the central engine of GRBs. Such emission may arise from synchrotron self-Compton scattering in the internal shocks. We derived the intrinsic cutoff energy due to pair production through interactions of high-energy photons with softer photons from the same radiation regions. Our intrinsic cutoff energy is consistent with the observed spectra of low-redshift GRBs such as GRB 940217 and GRB 970417a. However, another cutoff energy appears in the $`\gamma `$-ray spectra. This cutoff results from interactions of high-energy photons with external cosmic infrared background photons. The electron/positron pairs produced during such interactions own a significant fraction of the explosion energy and they will Compton scatter the CMB photons. We expanded on the Cheng & Cheng (1996) model by deriving the emission spectrum and duration in the standard fireball shock model. The resulting emission will be able to reach the observer in $`10^3`$ seconds. The time-integrated scattered photon spectrum is $`N_\nu ^{\mathrm{sc}}\nu ^{(p+6)/4}`$. This is slightly harder than the internal shock emission spectrum, $`N_\nu \nu ^{(p+2)/2}`$. For typical parameters of the internal shock model, the observed energies of the scattered photons are in the MeV-GeV range, and the ratio of the scattered to synchrotron photon number per unit frequency, $`N_\nu ^{\mathrm{sc}}/N_\nu `$, is typically a few.

Of course, whether or not there is such delayed emission depends on the parameters of the internal shock model. Some parameters (e.g., $`L`$, $`ϵ_e`$ and $`ϵ_B`$) adopted in this paper are consistent with detailed fits to the multi-wavelength data of a few afterglows (Wijers & Galama 1999; Freedman & Waxman 2001; Panaitescu & Kumar 2002). Two other parameters (e.g., $`p`$ and $`\mathrm{\Delta }t`$) are consistent with the BATSE observations (Preece et al. 2000). The bulk Lorentz factor of the wind, $`\mathrm{\Gamma }`$, is a crucial parameter because $`_0^{\mathrm{in}}`$ is strongly dependent of $`\mathrm{\Gamma }`$ (see equation ). The observed TeV emission from GRB 970417a requires that $`\mathrm{\Gamma }`$ be equal to or larger than 700, provided that the other parameters are given in the simple example of §3. Furthermore, some arguments on the optical flash of GRB 990123 in the reverse shock model (Wang, Dai & Lu 2000; Soderberg & Ramirez-Ruiz 2002) show that $`\mathrm{\Gamma }`$ may be of the order of $`10^3`$ (but also see Sari & Piran , who suggested that $`\mathrm{\Gamma }`$ in this burst could be about 200). If $`\mathrm{\Gamma }`$ is close to or larger than 600, we can indeed see a delayed MeV-GeV emission spectrum. If $`\mathrm{\Gamma }<600`$, however, the intrinsic cutoff energy could be less than the external pair-production minimum energy at which a photon is locally attenuated through interaction with the infrared background radiation, so that there could not be the delayed MeV-GeV emission discussed here.

It is seen from equations (6) and (14) that the Compton parameter $`Y`$ determines the TeV emission from synchrotron self-Compton scattering in internal shocks, and thus the spectral photon number of the delayed MeV-GeV emission. In addition, additional TeV photons produced by other mechanisms, e.g., photo-pion and inelastic proton-neutron collisions (Waxman & Bahcall 1997; Böttcher & Dermer 1998; Derishev, Kocharovsky & Kocharovsky 1999; Bahcall & Mészáros 2000; Dai & Lu 2001), and synchrotron self-Compton scattering both in external reverse shocks (Wang, Dai & Lu 2001a, 2001b) and in early-time forward shocks (Dermer, Chiang & Mitman 2000; Zhang & Mészáros 2001), may be absorbed in the cosmic infrared photon field, leading to electron/positron pairs. The inverse Compton scattering of such pairs off CMB photons may have a non-negligible contribution to the delayed MeV-GeV emission.

The MeV-GeV emission studied in this paper will be detectable by next-generation $`\gamma `$-ray satellites such as the Gamma-Ray Large Area Space Telescope (GLAST). In the Swift-GLAST era, a plenty of GRBs localized by the Swift satellite are expected to be detected by GLAST. Such observations, particularly on the delayed MeV-GeV emission and higher-energy spectral cutoff, may provide a probe of the cosmic infrared background radiation. This in turn may help to constrain some of the most fundamental uncertainties in physical models of the star formation.

We are very grateful to the referee for helpful comments and criticism that enabled us to improve the manuscript significantly, and to Y. F. Huang, S. Kobayashi, X. Y. Wang, D. M. Wei and B. Zhang for valuable discussions. This work was supported by the National Natural Science Foundation of China (grants 19825109 and 19973003) and the National 973 Project (NKBRSF G19990754).