text:
# Delocalizing effect of the Hubbard repulsion for electrons on a two-dimensional disordered lattice

## I Introduction

The interplay between disorder and electron-electron interactions has been a subject of intense activity in the last few years. Recent experiments have shown the existence of an apparent metal-insulator transition (MIT) in two-dimensional (2D) semiconductor devices . This observation came as a surprise to the community, since the scaling theory of localization, developed for disordered, non-interacting systems predicts insulating states even for infinitesimal disorder strengths . Therefore, these experiments have promoted intense theoretical activity. There is at present no consensus and considerable controversy surrounds this problem . A novel feature of the experimental high mobility samples investigated is their exceptionally low electronic density $`n_s`$. This leads to an unusually high value of the dimensionless parameter $`r_sn_s^{1/2}`$ of up to 80. The parameter $`r_s`$ sets the scale of electron-electron interaction $`E_{ee}`$ as compared to the Fermi energy $`E_F`$ through $`r_sE_{ee}/E_F`$. Thus the question of electron-electron interaction effects becomes important in these disordered systems.

The analytical treatment of the problem of disordered, interacting electrons is possible only in some limiting cases. The effects of weak interactions in disordered systems have been studied in great detail in the metallic (delocalized) regime . The other extreme, corresponding to the strongly localized system can be treated by mean-field methods . Although these analytical approaches have provided many useful physical results, the general treatment of disordered quantum many-body systems remains an unsolved problem. Indeed, one of the most powerful methods developed for the analytical treatment of disordered systems, namely, the supersymmetry approach , cannot handle the effects of electron-electron interactions. Given this context, numerical approaches play a crucial role in the treatment of disordered, interacting systems.

Several numerical approaches have been applied to the study of 2D disordered, strongly correlated systems . Among the approximate approaches, Hartree-Fock calculations with residual interactions similar to the configuration interaction approaches of quantum chemistry have been applied to systems with spinless fermions and electrons with spin . Exact diagonalization approaches have been useful but suffer from severe limitations in the accessible system size and number of particles . From these studies, it has been observed that repulsive electron-electron interactions can have a delocalizing effect in small systems. In addition, it has been shown experimentally  that the spin degrees of freedom play a crucial role in the physics of these strongly interacting systems. The inclusion of the spin degrees of freedom renders the numerical calculations even more difficult and strongly reduces the number of fermions accessible (see e.g. ).

Quantum Monte Carlo (QMC) approaches provide a powerful alternative to the treatment of quantum many-body systems. These methods are in principle exact, apart from statistical errors and allow the treatment of much larger system sizes with many particles. The finite-temperature determinantal QMC approach has been applied to the two-dimensional disordered Hubbard model and signatures of a metal-insulator transition were obtained . However, these calculations were carried out at finite temperature and could not access the ground state of the system. In previous work, we studied the ground state properties of the disordered 2D Hubbard model by the projector quantum Monte Carlo (PQMC) method . We studied the properties of the Green’s function, charge density and inverse participation ratio against model parameters and system size. While we observed some local charge reorganization, we could not detect any significant delocalizing influence of the Hubbard repulsion $`U`$ on the many-body ground state. However, it should be noted that all the physical quantities studied in Ref.  were obtained by integrating over all particles and effectively, the entire energy spectrum. Here we introduce a new approach which allows us to study the properties of a single added particle and therefore emphasizes the physical effects of interactions in the proximity of the Fermi edge.

This approach uses the inverse participation ratio (IPR) extracted from the charge density differences of two many-body ground states, as described below. The IPR, $`\xi `$, for a normalized single-particle wavefunction $`\psi (i)`$ is given by $`\xi =[_i|\psi (i)|^4]^1`$, where $`i`$ is the site index of the system. Clearly, $`|\psi (i)|^2`$ can be identified as the one-particle charge density at the site $`i`$. This definition is usually carried over to many-body systems by renormalizing the total charge density at site $`i`$, $`\rho (i,N_p)`$, which is obtained in the standard way from the ground state many-body wavefunction for $`N_p`$ particles. The IPR $`\xi `$ is then calculated with the effective charge density $`\rho (i,N_p)/N_p`$. We found that the IPR obtained from such a definition showed slight variations with model parameters and practically no variation with system size . We believe that the physical reason for this weak variance is the fact that this procedure effectively integrates over all energies and therefore the dominant contribution comes from the states deeply below the Fermi energy, which remain strongly localized even in the presence of interactions. Therefore, it is necessary to find a method which is more sensitive to the contribution of states in the vicinity of the Fermi energy. Recently, we studied the localization properties of the disordered, attractive Hubbard model in two- and three-dimensions . In Ref. , we showed that an IPR calculated for an added pair of particles is a very sensitive and relevant quantity to study the localization properties of the wavefunction. Based on our results for the attractive Hubbard model, we have introduced a related quantity in the repulsive case, the IPR for a single added particle, which turns out to be much more sensitive to the effects of interaction. This quantity has certain similarities to the single-particle tunneling amplitude. The latter has recently been studied by exact diagonalization methods for small clusters with spin-1/2 fermions . This work provided indications for a delocalizing effect induced by repulsive interactions in disordered systems. However, the number of particles studied was rather restricted due to the limitations of the exact diagonalization methods. With our method and extensive, highly accurate PQMC simulations, we study considerably larger numbers of particles in presence of strong electron-electron interactions. In this study we report a significant delocalizing influence of the Hubbard repulsion on 2D disordered electronic systems. This paper is organized as follows : after this Introduction, we describe the method used and the tests performed in the next section. Our results and discussion are presented in detail in the third section.

## II Model and method

The Hamiltonian studied in this paper is the disordered, repulsive Hubbard model given by:

$$\begin{array}{c}H=H_A+H_I\hfill \\ =\left(t\underset{ij,\sigma }{}c_{i,\sigma }^{}c_{j,\sigma }+\underset{i,\sigma }{}ϵ_ic_{i,\sigma }^{}c_{i,\sigma }\right)+U\underset{i}{}n_in_i\hfill \end{array}$$
(1)
where $`c_{i\sigma }^{}`$ ($`c_{i\sigma }`$) creates (destroys) an electron at site $`i`$ with spin $`\sigma `$ and $`n_{i\sigma }=c_{i\sigma }^{}c_{i\sigma }`$ is the corresponding occupation number operator. The hopping term $`t`$ between nearest neighbor lattice sites characterizes the kinetic energy and the random site energies $`ϵ_i`$ are taken from a box distribution over $`[W/2,W/2]`$. The parameter $`U`$ measures the strength of the screened, repulsive Hubbard interaction ($`U>0`$). We have considered both the one- and two-dimensional cases, with periodic boundary conditions in all directions. In the 2D case, the sites $`i`$ lie on a rectangular lattice of linear dimension $`N_x,N_y`$. The system size $`N=N_x\times N_y`$ then follows accordingly in one- and two-dimensions. In the limit $`U`$ = 0, this Hamiltonian reduces to the Anderson model (given by $`H_A`$) which is a standard model for the study of disordered systems . In the absence of disorder, $`W`$ = 0, this Hamiltonian reduces to the usual Hubbard model, which is one of the best-studied models for correlated electronic systems .

We have studied this Hamiltonian by the projector quantum Monte Carlo (PQMC) method  and exact diagonalization calculations. The PQMC method was initially developed to study the ground state of the Hubbard model (clean limit of Eqn. (1)). The method can be generalized, in principle quite simply, to incorporate disorder via the random site energies. However, the actual implementation and convergence of the algorithm  is highly non-trivial compared to the pure case, as will be discussed in detail below.

The PQMC method consists in filtering out the true ground state $`|\psi _0`$ of the many-body system from an appropriately chosen trial function $`|\varphi `$:

$$|\psi _0=\underset{\mathrm{\Theta }\mathrm{}}{lim}\frac{e^{\mathrm{\Theta }H}|\varphi }{\sqrt{\varphi |e^{2\mathrm{\Theta }H}|\varphi }}.$$
(2)
This method is exact in principle, apart from statistical errors and the sign problem which appears for fermions at $`U>0`$. The Hamiltonian plays the role of the projection operator through the term $`e^{\mathrm{\Theta }H}`$, where $`\mathrm{\Theta }`$ plays the role of the projection parameter. The trial wave-function is usually formed as a product of up and down spin states from the eigenstates of the non-interacting Hamiltonian. In our case, we chose the Fermi sea of $`H_A`$ as the trial wave-function. In the PQMC procedure, the projection operator $`\mathrm{exp}(\mathrm{\Theta }H)`$ is first Trotter decomposed as $`\left(\mathrm{exp}(\mathrm{\Delta }\tau H_A)\mathrm{exp}(\mathrm{\Delta }\tau H_I)\right)^L`$ with $`\mathrm{\Theta }=\mathrm{\Delta }\tau \times L`$. This introduces a systematic error of order $`(\mathrm{\Delta }\tau )^2`$ due to non-commutation of $`H_A`$ and $`H_I`$. We have used the symmetric Trotter decomposition, which introduces a systematic error of $`(\mathrm{\Delta }\tau )^3`$. The interaction is then decoupled by a discrete Hubbard-Stratonovich transformation, by the introduction of $`N\times L`$ Ising-like fields. This Ising model with complicated effective interactions is then treated by a Monte Carlo (MC) procedure to obtain the ground state properties of the system. The quantity calculated during the simulation is the zero-temperature, equal-time Green’s function, $`G_{ij}=_\sigma \psi _0|c_{i,\sigma }^{}c_{j,\sigma }|\psi _0`$. During the simulation this Green’s function can be used to obtain all the other static correlation functions. In the algorithm used $`O(N^2)`$ operations are required to update the Green’s function after a MC step. This procedure introduces cumulative errors and therefore the Green’s function has to be recalculated from scratch regularly during the simulation (every $`L_C`$ steps), which requires $`O(N^3)`$ operations. When the projection parameter $`\mathrm{\Theta }`$ becomes large, which is necessary for good convergence, the different components of the wavefunction tend to become parallel during the projection process. Therefore, it is necessary to reorthogonalise the components of the wavefunction regularly, every $`L_R`$ time steps. In addition, it is well known that quantum simulations of fermionic systems suffer from the sign problem except in some special cases. However, it has been observed that disorder in fact diminishes the magnitude of the sign problem. In our simulations, the sign problem is well under control, with the number of negative signs being less than $`1\%`$ of the total number of steps considered.

We have studied systems with particle number ($`N_p`$) up to 25 fermions on lattice sizes of up to $`6\times 8`$ sites, with particle density always around quarter-filling. The simulations were carried out in the $`S_z=0`$ sector for even numbers of particles and the $`S_z=1/2`$ sector for odd number of particles. The specific quantity studied in this paper is the difference of charge densities $`\delta \rho (i)`$ between systems with $`N_p`$ and $`N_p+1`$ particles (see the next section for more details). Here, we discuss the convergence of the PQMC calculations. The charge densities involved in the calculation $`\delta \rho (i)`$ are obtained from two independent simulations of the same disorder realization with different particle numbers, $`N_p`$ and $`N_p+1`$. Therefore, it becomes necessary to measure very accurately the distribution of this added particle over the lattice. Clearly, it is harder to measure $`\delta \rho (i)`$ accurately than measuring the total charge densities $`\rho (i,N_p)`$ in the ground state with $`N_p`$ particles. We carried out extensive tests to verify the quality of our $`\delta \rho (i)`$ data, by varying the PQMC parameters until convergence was obtained. We have tested for convergence by varying the following PQMC parameters: the projection parameter $`\mathrm{\Theta }`$, the Trotter time-step $`\mathrm{\Delta }\tau `$, the Monte Carlo parameters, the reorthogonalization interval $`L_R`$ and the interval to recalculate the Green’s function $`L_C`$. We have gone up to $`\mathrm{\Theta }=15`$ and $`\mathrm{\Delta }\tau =0.05`$. For the physical parameters used in this paper, we find that $`\mathrm{\Theta }=10`$ with $`\mathrm{\Delta }\tau =0.1`$ are sufficient. This corresponds to a systematic error of $`10^3`$. After varying the parameters $`L_R`$ and $`L_C`$ we decided to recalculate the Green’s function with reorthogonalization of the components, after every 5 time-steps ( $`L_R=L_C=5`$). We have also checked the MC parameters and find that 3000 sweeps are adequate for convergence, with 1000 sweeps for equilibration.

The disorder average was carried out over $`N_R`$ different disorder realizations with $`N_R=16`$ in the PQMC simulations and $`N_R=100`$ for most exact calculations. The site-energies are randomly chosen for $`N_R`$ disorder realizations at $`W/t=1`$ and then scaled proportionally to $`W/t`$ for stronger values of disorder.

While our main motivation is the study of the 2D case, there are many useful reasons to study 1D systems. The 1D case can be studied conveniently by exact diagonalization methods, while changing the system size. Thus, we can study the variation of various properties such as the inverse participation ratio (IPR) $`\xi `$, with system size, exactly. These calculations provide a strong independent check on the PQMC data at small sizes. Further, the localization effect is much stronger in 1D and $`\delta \rho `$ is usually strongly peaked in 1D as compared to 2D. Therefore, if the algorithm is capable of reproducing a localized peak, this provides a strong check on the method, at the given range of physical parameters (examples provided in the next section).

## III Results and Discussion

In order to study the localization properties of the system, we use the charge density distribution for an added particle at the Fermi level, given by: $`\delta \rho (i)=\rho (i,N_p+1)\rho (i,N_p)`$, where $`\rho (i)`$ is the ground state charge density at site $`i`$. The values of $`\rho (i,N_p)`$, $`\rho (i,N_p+1)`$ are obtained from two independent PQMC simulations for the same disorder realization. At $`U`$ = 0, this $`\delta \rho (i)`$ is identically equal to the one-particle probability distribution (added particle) at the Fermi edge. In the interacting case, this charge density distribution ($`\delta \rho (i)`$) is not generally equal to a probability distribution. For example, it is not necessarily positive definite. However, in all the cases studied we find that this quantity is always positive. Therefore, $`\delta \rho (i)`$ can be considered to be similar to a one-particle probability distribution. Thus, we can associate an inverse participation ratio (IPR) for an added particle, for a given disorder realization as: $`\xi =(_i\delta \rho (i)^2)^1`$. The disorder averaged IPR is given by $`\xi `$. At $`U=0`$, the IPR is a standard tool used to obtain the number of sites over which the particle is localized . We have successfully used this approach in previous work on the disordered, attractive Hubbard model for a quantitative description of the localization properties of the ground state .

In Figs. 1-2, we study the behavior of the IPR as a function of $`U/t`$ and $`W/t`$ to try and establish the interesting range of physical parameters in this system. In Fig. 1, we compare the IPR $`\xi `$ for a 1D ring of 12 sites and a $`4\times 3`$ lattice in 2D as a function of interaction strength $`U/t`$. From the figure we see that increasing $`U/t`$ from 0 tends to increase the IPR up to intermediate strengths of the interaction. Thus, the optimal value appears to be around $`U/t=4`$ in both 1D and 2D. Indeed, for $`U/t>4`$, the value of IPR starts to decrease. In Fig. 2, we see the variation of $`\xi `$ for small system sizes for $`5W/t12`$. Since the ratio $`\xi (U/t=4)/\xi (U/t=0)`$ is practically constant, this justifies our choice of parameter range for $`W/t`$. Disorder strength $`W/t<5`$ would lead to states with localization length larger than the accessible system size even at $`U/t=0`$. Convergence of the PQMC calculations becomes progressively more difficult for larger disorder strengths. Therefore, we studied the variation of the IPR for $`W/t=5,7`$ in the 2D case.

Thus, it would seem ideal to study the IPR of 1D and 2D systems for $`5W/t7`$ and $`U/t4`$. In 1D we could access these parameter values conveniently for the system sizes studied by exact calculations. However, this becomes much more difficult in 2D. This is because we study a small difference of large total charge densities in the PQMC simulations. We found that the accuracy of our method for $`\delta \rho (i)`$ is good for for values of the interaction $`U/t2`$. While we are not exactly at the optimal value of $`U/t`$, we are nevertheless in a region with sufficiently strong interactions. Indeed, it can be seen from Fig. 1 that $`U/t=2`$ already has a substantial delocalizing influence on the system. In Fig. 3, we present the PQMC calculation of $`\delta \rho (i)`$ in 1D at $`U/t=2`$ as compared to exact calculations. We note that the $`\delta \rho (U=0)`$ is strongly peaked. Increasing $`U/t`$ to 2 radically changes the picture and shifts the peak completely. The PQMC curve shown reproduces the quantitative picture of charge density difference. It should be noted that the calculation begins with a trial wavefunction corresponding to the $`U/t=0`$ data and changes completely to give the correct physical picture. Quantitatively, the usual errors seen at these values of the physical parameters were around $`35\%`$. The 2D case is more delocalized compared to 1D and we have observed that convergence is better in 2D. For example, we have $`\xi `$ = 7.19 (exact) and 7.43 (QMC) for a $`4\times 3`$ lattice at $`U/t=2`$ and $`W/t=7`$, averaged over 16 disorder realizations. This corresponds to an error of about $`3\%`$ for the most extreme parameter values studied. However, for $`U/t=4`$, the error increases to $`810\%`$ and therefore we restrict our studies to $`U/t2`$ in 2D.

Since we study a model of electrons with spin, it is interesting to consider even-odd effects in the particle number $`N_p`$ and system size $`N`$. In Fig. 4, we consider the effect of progressively adding one and two particles to a 12 site ring with 6 electrons. From Fig. 3 and Fig. 4, we see that the first added electron at $`U/t`$ = 2 and 4, has an entirely different peak as compared to the non-interacting case. In Fig. 4, we see that the second added electron also occupies a different region, as seen from the position of the peak. This leads to the following physical image for the 1D case : each added electron finds its own location in space with optimal energy and it is well localized by external disorder and random distribution of charges of other electrons. This is in contrast with the non-interacting case, where the non-interacting orbital remains the same for odd and next even electron.

The situation is qualitatively different in 2D. Indeed, in Fig. 5, we show the charge density difference for one and two added particles on a $`6\times 6`$ lattice with 18 electrons. This quantity is obtained from PQMC simulations of the system for a particular disorder realization, taken at $`W/t`$ = 7. The data show that the initial $`U/t=0`$ configuration is very clearly localized for the given disorder realization and lattice size. It is seen from the figure that the introduction of a repulsive Hubbard interaction ($`U/t=2`$) leads to a substantial delocalization of the added particle. This is also borne out quantitatively, since the IPR increases practically by a factor of 3, as compared to $`U/t=0`$. In the non-interacting case, both added particles occupy the same orbital. Therefore, the peak for the second added particle, is trivially identical to the first. It is remarkable that in the interacting case, the peak is transformed to a much more extended distribution over the lattice. Furthermore, the second added particle has practically the same distribution $`\delta \rho (i)`$ as the first. Thus, there appears to be almost zero effective repulsion between the two added particles, despite the interaction strength $`U/t=2`$. This is completely in contrast to the scenario in 1D with interactions, where we observe significant repulsion between the two added particles. We note that this is the case for every disorder realization studied.

We now turn to a more quantitative picture of the one and two dimensional systems. For this, it is necessary to average the data over different disorder realizations and to vary the system size. We start the analysis from the data for the 1D case. In Fig. 6, we present the IPR for interacting particles on 1D rings of 4–12 sites at constant density of particles (quarter filling) and compare against the variation of the IPR at $`U/t=0`$. We find a delocalization effect even for 1D rings. We note that the curve for $`U/t=0`$ is well saturated in 1D even at 8–12 sites. The $`U/t=2`$ curve shows signs of saturation. The delocalization effect is most visible for $`U/t=4`$, as expected from Fig. 1, where unfortunately, we do not have access to PQMC data for comparisons. The ratio $`\xi (U/t)/\xi (U/t=0)`$ goes to 1.64 for $`U/t=2`$ and 1.94 for $`U/t=4`$.

In 2D previous exact diagonalization studies have gone up to 6 electrons on a $`4\times 4`$ lattice, or 4 electrons on a $`6\times 6`$ lattice . Therefore, it is of great importance to access larger system sizes with more particles. In Fig. 7, we present the IPR for different lattice sizes in 2D, obtained from PQMC simulations. We have gone up to 48 sites ($`6\times 8`$ lattice) and 24-25 particles, which goes beyond any existing study of the ground state in the literature. We observe a delocalizing effect of the Hubbard repulsion manifested by an increase of the IPR with system size and interaction strength. There is a remarkable difference between the present result and the IPR for the full many-body ground state wavefunction obtained in Ref. . In the previous work, we could notice no size effects over a large range of $`W`$ and lattice sizes. In fact, the curves for different lattice sizes against $`U/t`$ at a given value of $`W/t`$ tended to collapse as seen in Figs. 8 and 9 of Ref. . On the contrary, in the present work, we see from Fig. 7 that effects of lattice size on the IPR for an added particle are considerable. Furthermore, we see that the Hubbard repulsion has a strong delocalizing influence, compared to the data for the non-interacting case. As explained in the previous sections, we attribute the difference between the present results and those presented in Ref.  to the fact that the IPR from $`\delta \rho `$ measures directly the response in the vicinity of the Fermi level.

The ratio $`\xi (U/t)/\xi (U/t=0)`$ can be considered as a quantitative measure of delocalization effect induced by repulsive interaction. From the inset of Fig. 7, it can be seen that this ratio rises with system size for a given density and we do not observe saturation at the system sizes studied. Our data for this ratio for small clusters in 2D are comparable to the enhancement ratio obtained in Ref.  though the quantity studied in Ref.  was slightly different i.e. the one-particle tunneling amplitude. From Figs. 6 and 7 we can compare the values of the ratio $`\xi (U/t)/\xi (U/t=0)`$ in 1D and 2D. At the system sizes studied these values are comparable for $`U/t=2`$. However, it would necessary to analyze the behavior of this ratio for larger system sizes in order to get a clear answer to whether the localization properties are indeed different in 1D and 2D. At the same time, we note that there seem to be qualitative differences based on even-odd effects for added particles in 1D and 2D (see Fig. 4 and Fig. 5 and discussion there). This difference favors the picture of stronger delocalization in 2D as compared to 1D.

Even if our data clearly show a repulsion induced delocalization effect, they do not permit us to draw a definite answer about the existence of a metal-insulator transition for this system in the thermodynamic limit. Indeed, even though we have a significant number of fermions, we cannot go to larger system sizes because of the accuracy constraint of calculating charge differences. It should also be pointed out that our calculations done at $`U/t2`$ are below the optimal value of the interaction strength. We expect that the delocalizing effect is even stronger at $`U/t=4`$.

In conclusion, we have studied the ground state properties of the repulsive Hubbard model with disorder through the powerful PQMC method. Highly accurate simulations permit us to obtain the difference of charge density between two ground states with $`N_p`$ and $`N_p+1`$ fermions respectively. The analysis of this characteristic clearly shows that the Hubbard repulsion has a delocalizing effect in the system. We have observed some qualitative differences between 1D and 2D for this characteristic. However, the restrictions on system size and interaction strength do not permit us to draw a definite conclusion about the existence of a MIT in the thermodynamic limit in 2D.

We thank the IDRIS at Orsay for access to their supercomputers.


format: text

text:
# Comment on “Systematics of the Induced Magnetic Moments in 5d Layers and the Violation of the Third Hund’s Rule”

Recently, Wilhelm et al. Wilhelm  reported x-ray magnetic circular dichroism (XMCD) measurements at the W and Ir $`L_{2,3}`$ edges of W/Fe and Ir/Fe multilayers. The $`5d`$ spin and orbital magnetic moments, $`\mu _S`$ and $`\mu _L`$, of W and Ir were deduced by applying the sum rules. These moments are induced in the $`5d`$ band, through hybridization, by the spin polarization of the Fe, which is reflected in their small values. For Ir, $`\mu _S`$ and $`\mu _L`$ align both parallel to the Fe spin moment, while for W they align both anti-parallel to Fe. As a result, it was suggested that in the case of W, with its less than half filled $`5d`$ band, this implies a violation of Hund’s third rule. It seems a misnomer to call this a violation of Hund’s rule, since Hund’s rules do not generally apply to induced moments, being only strictly valid for isolated atoms. However, on the basis of their results for only two systems, the authors further concluded: “This remarkable finding shows that the induced magnetic behavior of $`5d`$ layers may be radically different than that of impurities and alloys.” Here we show, on the basis of a systematic set of calculations, that the behavior of $`\mu _S`$ and $`\mu _L`$ of $`5d`$ magnetic impurities in an Fe host  is not radically different than that of $`5d`$ interface layers in Fe multilayers, proving that their conclusion, inferred on the basis of two examples, is not correct.

Calculations were performed for the electronic and magnetic properties of periodic multilayers Fe/Z consisting of 5 monolayers (ML) of Fe(100) and 3 ML of the $`5d`$ transition metal, Z, using the relativistic spin-polarized linear muffin-tin orbitals (LMTO) method. Ebert  Figure 1 includes, besides the present results for the $`5d`$ interface layer, two different impurity calculations. These refer to a non-self-consistent fully-relativistic calculation using a self-consistent scalar-relativistic potential Ebert2  and a self-consistent fully-relativistic calculation. Kornherr  The $`\mu _S`$ of $`5d`$ impurities in Fe Kornherr ; Tyer  is aligned anti-parallel for the first part of the $`5d`$ transition metal series, i.e. for less than half filled band. From Os onwards $`\mu _S`$ of the $`5d`$ impurity and Fe are parallel aligned. For $`5d`$ impurities to obey Hund’s third rule, $`\mu _L`$ should be positive throughout the $`5d`$ series: this is not the case for Re, Os, and Ir.

The curves of the calculated $`\mu _S`$ and $`\mu _L`$ for the multilayers show a similar behavior as for the impurities. $`\mu _L`$, being more sensitive to the structure and chemical environment, switches between parallel and anti-parallel at nearly the same location, resulting in a negative value for W, Re, and Os. Therefore, the relative alignment of $`\mu _S`$ and $`\mu _L`$ in the $`5d`$ metals is the same as in the single impurity case, with the exception of W and Ir.

The present calculations agree well with the experimental data which are for $`\mu _S`$ and $`\mu _L`$ $`(0.18,0.016)`$ and $`(0.20,0.019)`$ in W and Ir respectively, while the theory gives $`(0.11,0.031)`$ and $`(0.26,0.003)`$. It just so happens that in the vicinity of W and Ir a sign change of $`\mu _L`$ occurs (Fig.1). Thus it seems obvious why a systematic study should not be based on these two data points alone. Our results indicate that mostly band filling, and to a lesser extent geometrical effects, determine the trend of $`\mu _S`$ and $`\mu _L`$ over the $`5d`$ series. This finding is corroborated by the fact that a similar trend also occurs for the $`4d`$ series. Tyer

R. Tyer, G. van der Laan, W.M. Temmerman, Z. Szotek
Daresbury Laboratory
Daresbury, Warrington WA4 4AD, UK

Received 20 March 2002

PACS numbers: 75.70.Cn, 75.70.-i, 78.70.Dm


text:
# Untitled Document

Standard Clocks, Orbital Precession and
the Cosmological Constant
Andrew W. Kerr, John C. Hauck, and Bahram Mashhoon
Department of Physics and Astronomy
University of Missouri-Columbia

Columbia, Missouri 65211, USA

Abstract

We discuss the influence of the cosmological constant on the gravitomagnetic clock effect and the gravitational time delay of electromagnetic rays. Moreover, we consider the relative motion of a binary system to linear order in the cosmological constant $`\mathrm{\Lambda }`$. The general expression for the effect of $`\mathrm{\Lambda }`$ on pericenter precession is given for arbitrary orbital eccentricity.

PACS numbers: 0420, 9880

1. Introduction

Current cosmological models that take the acceleration of the universe into account involve a cosmological constant $`\mathrm{\Lambda }`$ of magnitude $`\mathrm{\Lambda }H_0^2/c^2`$, where $`H_0`$ is the Hubble parameter. Taking $`H_070\mathrm{km}\mathrm{sec}^1\mathrm{Mpc}^1`$, we find that $`\mathrm{\Lambda }10^{56}\mathrm{cm}^2`$. The local effects of such a cosmological constant are expected to be very small. Nevertheless, great strides are being made in observational techniques in astronomy and it appears worthwhile to investigate some of the observational consequences of the existence of a cosmological constant.

The dependence of the gravitomagnetic clock effect on the cosmological constant is examined in section 2. Section 3 is devoted to the gravitational time delay and its dependence on $`\mathrm{\Lambda }`$. The influence of the cosmological constant on the motion of a binary system is investigated in detail in sections 4 and 5. A bound orbit is shown to precess with a frequency given by $`\frac{1}{2}(\mathrm{\Lambda }c^2/\omega )\sqrt{1e^2}`$, where $`\omega `$ and $`e`$ are the orbital frequency and eccentricity, respectively. Section 6 contains a brief discussion of our results.

2. Clock effects in the Kerr-de Sitter spacetime

The gravitomagnetic clock effect is caused by the net spin of a gravitational source \[1-4\]. In general, the rotation of a massive body introduces a difference in the co- and counter-rotating orbital periods of test masses in motion around the central source. We wish to investigate this gravitomagnetic effect in the equatorial plane of the Kerr-de Sitter spacetime to find what influence the cosmological constant would have on the clock effect. The Kerr-de Sitter metric in the standard Boyer-Lindquist coordinates $`x^\mu =(t,r,\theta ,\varphi )`$ is given by

$`ds^2=`$ $``$ $`\left[1{\displaystyle \frac{2Mr}{\mathrm{\Sigma }}}{\displaystyle \frac{\mathrm{\Lambda }}{3}}(r^2+a^2\mathrm{sin}^2\theta )\right]dt^22a\left[{\displaystyle \frac{2Mr}{\mathrm{\Sigma }}}+{\displaystyle \frac{\mathrm{\Lambda }}{3}}(r^2+a^2)\right]\mathrm{sin}^2\theta dtd\varphi `$ (1)
$`+`$ $`{\displaystyle \frac{\mathrm{\Sigma }}{\mathrm{\Delta }}}dr^2+{\displaystyle \frac{\mathrm{\Sigma }}{\chi }}d\theta ^2+\left[{\displaystyle \frac{2Mr}{\mathrm{\Sigma }}}a^2\mathrm{sin}^2\theta +(1+{\displaystyle \frac{\mathrm{\Lambda }}{3}}a^2)(r^2+a^2)\right]\mathrm{sin}^2\theta d\varphi ^2,`$ (3)
where units are chosen such that $`G=c=1`$, unless specified otherwise, and

$$\mathrm{\Sigma }=r^2+a^2\mathrm{cos}^2\theta ,\chi =1+\frac{\mathrm{\Lambda }}{3}a^2\mathrm{cos}^2\theta ,$$
(4)
$$\mathrm{\Delta }=r^22Mr+a^2\frac{\mathrm{\Lambda }}{3}r^2(r^2+a^2).$$
(5)
The mass of the source is $`M,J=Ma`$ is its angular momentum and $`\mathrm{\Lambda }`$ is the cosmological constant. The metric (1) represents a rotating Kottler spacetime . The Kottler spacetime is also known as the Schwarzschild-de Sitter spacetime. Useful background information on de Sitter spacetime and the cosmological constant can be found in  and , respectively.

It is interesting to note that in the Newtonian limit the physical content of equation (1) reduces to a Newtonian gravitational potential $`\mathrm{\Phi }_N=GMr^1\frac{1}{6}\mathrm{\Lambda }c^2r^2`$, which satisfies the Poisson equation $`^2\mathrm{\Phi }_N=4\pi G\rho _N`$ with $`\rho _N=M\delta (𝐫)+\rho _\mathrm{\Lambda }`$. Here $`\rho _\mathrm{\Lambda }=c^2\mathrm{\Lambda }/(4\pi G)`$ is the effective negative uniform density of “matter” represented by the cosmological constant. The effective repulsive force due to this source is then given by $`\frac{1}{3}\mathrm{\Lambda }c^2𝐫`$ per unit mass. On the other hand, the relativistic interpretation of the cosmological constant in terms of vacuum energy involves a “substance” with equation of state $`\widehat{\rho }c^2+\widehat{p}=0`$ such that the density $`\widehat{\rho }`$ is positive, but the pressure is negative and is given by $`\widehat{p}=c^4\mathrm{\Lambda }/(8\pi G)`$.

Let us now consider circular geodesics of this spacetime. The stability of circular geodesic orbits in the equatorial plane of Kerr-de Sitter geometry has been investigated by Howes . It turns out that stable circular orbits exist out to a radius of $`(3M/\mathrm{\Lambda })^{1/3}`$. On the other hand, the inner stability limits vary depending on the sense of orbital motion as well as the values of $`a`$ and $`\mathrm{\Lambda }`$ . For $`\mathrm{\Lambda }=0`$, we recall that for the corotating case, the inner stability threshold varies from $`r=6M`$ to $`r=M`$ as $`a:0M`$, whereas it varies from $`r=6M`$ to $`r=9M`$ in the counterrotating case.

The geodesic equation for the radial coordinate reduces to

$$g_{tt,r}\left(\frac{dt}{d\varphi }\right)^2+2g_{t\varphi ,r}\left(\frac{dt}{d\varphi }\right)+g_{\varphi \varphi ,r}=0,$$
(6)
which has the solution

$$\frac{dt}{d\varphi }=a\pm \omega _K^1.$$
(7)
Here $`\omega _K`$ is the modified Kepler frequency

$$\omega _K^2=\frac{M}{r^3}\frac{1}{3}\mathrm{\Lambda },$$
(8)
which can be interpreted as being due to a net mass of $`M+(4\pi /3)\rho _\mathrm{\Lambda }r^3`$, where $`\rho _\mathrm{\Lambda }=\mathrm{\Lambda }/4\pi `$ is the constant Newtonian density mentioned before. In connection with equation (6) as well as other modifications of Kepler’s third law to include dark matter, etc., we note that the motion of planets in the solar system can be used to study deviations from Keplerian motion; however, for $`\mathrm{\Lambda }10^{56}\mathrm{cm}^2`$ the corresponding effects are too small to be detectable at present.

There are three distinct clock effects that can be studied using counterrevolving circular geodesic orbits. To find the influence of the cosmological constant on the clock effects, we consider these in turn. Integration of equation (5) over $`2\pi `$ for corotating and $`2\pi `$ for counterrotating test particles implies that $`t_\pm =T_K\pm 2\pi a`$, where $`T_K=2\pi /\omega _K`$ is the Keplerian period. Thus $`t_+t_{}=4\pi a`$, just as in the Kerr case. Consider now an (accelerated) standard clock at rest on the circle of radius $`r`$ in the equatorial plane. A standard clock measures the proper time along its world line. According to this clock, the difference in the time $`\tau ^{}`$ that it would take two free counterrevolving test particles to complete the orbit is

$$\tau _+^{}\tau _{}^{}=4\pi a\sqrt{1\frac{2M}{r}\frac{\mathrm{\Lambda }}{3}(r^2+a^2)},$$
(9)
which is the observer-dependent single-clock clock effect according to the terminology of . For the cosmological constant to have any significant influence on this clock effect, the orbital radius would have to be unreasonably large. In fact, $`\mathrm{\Lambda }10^{56}\mathrm{cm}^2`$ has a totally negligible effect even at a planetary radius of $`r10\mathrm{AU}`$, since $`\mathrm{\Lambda }r^2/310^{28}`$.

Consider now the proper periods of revolution $`\tau _\pm `$ of the two free standard clocks counterrevolving on the circular orbit of radius $`r`$; we find from equations (1) and (5) that

$$\tau _\pm =T_K\sqrt{1\frac{3M}{r}\frac{\mathrm{\Lambda }}{3}a^2\pm 2a\omega _K}.$$
(10)
Expanding this expression in powers of $`ϵ=a/M`$, we find that

$$\tau _+\tau _{}=\frac{4\pi a}{\sqrt{1\frac{3M}{r}}}\left[1+\frac{Ma^2(1\mathrm{\Lambda }r^2)}{2r(r3M)^2}+O(ϵ^4)\right].$$
(11)
Thus the cosmological constant does not affect the observer-dependent two-clock clock effect to first order in $`a/M`$.

Imagine now two standard clocks starting at the event characterized by $`(0,r,\pi /2,0)`$, i.e. $`t=0`$ at $`\varphi =0`$, and moving freely in opposite directions on the circle of radius $`r`$. Since the corotating clock is slower, their first meeting occurs at $`\varphi _1=\pi \alpha `$, where $`\alpha =a\omega _K`$. Let their nth meeting point be the event characterized by $`t_n`$ and $`\varphi _n`$, where $`t_n=\frac{1}{2}nT_K(1\alpha ^2)`$ and $`\varphi _n=n\pi (1\alpha )`$ modulo $`2\pi `$. The observer-independent two-clock clock effect refers to the difference in the proper times $`\tau _n^+\tau _n^{}`$ of the clocks at their nth meeting, where

$$\tau _n^\pm =\frac{1}{2}nT_K(1\alpha )\sqrt{1\frac{3M}{r}\frac{\mathrm{\Lambda }}{3}a^2\pm 2\alpha }.$$
(12)
To first order in $`a/M`$, this clock effect is also independent of $`\mathrm{\Lambda }`$ and is given by

$$\tau _n^+\tau _n^{}\frac{6\pi nJ}{\sqrt{r(r3M)}}.$$
(13)
The diametrical line joining the points of encounter with the origin of the circular orbit undergoes a precession in the opposite sense as the rotation of the source. The precession frequency is given approximately by $`n\pi \alpha /\tau _n^+`$, which to lowest order coincides with the precession frequency of a fixed torque-free test gyroscope in the equatorial plane at radius $`r`$.

The gravitomagnetic clock effect (9) far from the source, $`3M<<r<<(3M/\mathrm{\Lambda })^{1/3}`$, is essentially independent of $`\mathrm{\Lambda }`$, and is basically proportional only to $`a=J/M`$, which is the specific angular momentum of the source. That is, the effect has a topological character as it is practically independent of the radius of the orbit. It is also essentially independent of the gravitational constant $`G`$; therefore, the effect can be “large”. For instance, for space-borne clocks around the Earth the effect is $`10^7`$ sec; however, the detection of the effect requires knowledge of satellite orbits to millimeter accuracy . This seems to be somewhat beyond present capabilities by about an order of magnitude; moreover, the influence of the cosmological constant on the clock effect is too small to be detectable in the foreseeable future. Another aspect of the clock effect is that it takes longer to complete a prograde circular orbit in the equatorial plane than the corresponding retrograde orbit. Thus free motion in the same sense as the rotation of the source is slower than motion in the opposite sense, which is contrary to what one might expect from the “dragging of space” by the rotating source.

There is a close relationship between the clock effect and circular holonomy. In fact, circular holonomy in the Kerr-de Sitter spacetime can be developed in close analogy with that in the Kerr spacetime . This would involve the parallel transport of a vector around a constant-time circular $`\varphi `$-loop in the equatorial plane of the Kerr-de Sitter spacetime. Let us just mention here the result that there is a band of holonomy invariance for $`r_{}<r<(3M/\mathrm{\Lambda })^{1/3}`$ given by $`\stackrel{~}{n}f(r)=\stackrel{~}{m}`$, where $`\stackrel{~}{n}`$ and $`\stackrel{~}{m}`$ are positive integers and

$$f^2(r)=1\frac{2M}{r}\frac{2Ma^2}{r^3}\frac{M^2a^2}{r^4}\frac{1}{3}\mathrm{\Lambda }r^2+\frac{1}{3}\mathrm{\Lambda }a^2\left(1+\frac{2M}{r}\right)$$
(14)
to first order in the cosmological constant. Here $`r_{}`$ is the unique positive root of $`f(r)=0`$ for sufficiently small $`\mathrm{\Lambda }`$.

3. Time delay

Let us consider the time delay in the propagation of electromagnetic rays in the Kerr-de Sitter spacetime. The influence of the cosmological constant on the local propagation effects is expected to be very small; therefore, it is sufficient to consider the linearized Kerr-de Sitter metric in isotropic coordinates on the background Minkowski spacetime. In this case, the time that it takes for the rays to travel from the point $`P_1:(t,𝝆_\mathrm{𝟏})`$ to the point $`P_2:(t,𝝆_\mathrm{𝟐})`$ can be written as

$$t_2t_1=|𝝆_\mathrm{𝟐}𝝆_\mathrm{𝟏}|+\mathrm{\Delta }_{GE}+\mathrm{\Delta }_{GM}+\mathrm{\Delta }_\mathrm{\Lambda },$$
(15)
where $`𝝆`$ is the position vector in the background Minkowski spacetime and $`\mathrm{\Delta }_{GE}`$ and $`\mathrm{\Delta }_{GM}`$ respectively represent the standard gravitoelectric and gravitomagnetic time delays . Moreover, $`\mathrm{\Delta }_\mathrm{\Lambda }`$ is the additional time delay due to the cosmological constant. This delay can be calculated using the approach developed in ; that is,

$$\mathrm{\Delta }=\frac{1}{2c}_{P_1}^{P_2}h_{\alpha \beta }k^\alpha k^\beta 𝑑l,$$
(16)
where $`h_{\alpha \beta }=g_{\alpha \beta }\eta _{\alpha \beta },k^\alpha =(1,\widehat{𝐤}),\widehat{𝐤}`$ is a unit vector along the direction of propagation of the ray in the background Minkowski spactime and $`dl`$ is the element of straight line connecting $`P_1`$ to $`P_2`$.

Concentrating on the cosmological constant, we note that for $`M=a=0`$, equation (1) reduces to the de Sitter metric. Under the coordinate transformation

$$r(\rho )=\frac{\rho }{1+\frac{1}{12}\mathrm{\Lambda }\rho ^2},$$
(17)
the de Sitter metric takes the isotropic form

$$ds^2=\left(\frac{1\psi }{1+\psi }\right)^2dt^2+\frac{1}{(1+\psi )^2}(d\rho ^2+\rho ^2d\mathrm{\Omega }^2),$$
(18)
where $`\psi =\mathrm{\Lambda }\rho ^2/12`$. It follows from equation (16) that for the cosmological constant alone,

$$h_{00}=\frac{1}{3}\mathrm{\Lambda }\rho ^2,h_{0i}=0,h_{ij}=\frac{1}{6}\mathrm{\Lambda }\rho ^2\delta _{ij}.$$
(19)
Using these results in equation (14), we find

$$\mathrm{\Delta }_\mathrm{\Lambda }=\frac{\mathrm{\Lambda }}{12c}_{P_1}^{P_2}\rho ^2𝑑l.$$
(20)
It is straightforward to compute this integral and the result is

$$\mathrm{\Delta }_\mathrm{\Lambda }=\frac{\mathrm{\Lambda }}{36c}|𝝆_\mathrm{𝟏}𝝆_\mathrm{𝟐}|(\rho _1^2+𝝆_1𝝆_\mathrm{𝟐}+\rho _2^2).$$
(21)
We note that for $`\mathrm{\Lambda }10^{56}\mathrm{cm}^2,\mathrm{\Delta }_\mathrm{\Lambda }2`$ sec for a ray of radiation crossing the disk of our Galaxy.

The calculation of first-order gravitational delays in equation (13) assumes that the ray follows a straight line. Therefore, there must be in addition a geometrical delay that takes the actual deflected path into account; for weak lensing, the geometrical delay is of second order in the deflecting potentials. It is important to recognize that the cosmological constant does not participate in the deflection of light rays in the Kottler spacetime ; therefore, the cosmological constant does not contribute to the geometrical delay in the Schwarzschild-de Sitter spacetime.

It follows from equation (19) that locally, i.e. for distances much smaller than $`(3M/\mathrm{\Lambda })^{1/3}`$, the time delay due to the cosmological constant is very small compared to the corresponding Shapiro time delay. On the other hand, equation (19) indicates that for distances approaching $`\mathrm{\Lambda }^{1/2}`$, the time delay could be rather significant, but then our linear approximation scheme may break down. Thus for cosmological observations $`\mathrm{\Delta }_\mathrm{\Lambda }`$ should in general be taken into account, since the present uncertainty in the measurement of gravitational lensing time delay is about $`\frac{1}{2}`$ day .

4. Orbital precession

The solution of the equations of motion of a test particle in the Kerr-de Sitter spacetime is quite complicated. On the other hand, it turns out that one can obtain the main results regarding the gravitomagnetic clock effect, etc., from the linear approximation of general relativity. Therefore, we consider a linear post-Newtonian approach using a Lagrangian of the form $`=mds/dt`$, where $`m`$ is the mass of the test particle and the metric is expressed in $`(t,𝝆)`$ coordinates, where $`𝝆=(x,y,z)`$ indicates a point in space with isotropic (post-Newtonian) Cartesian coordinates. The equations of motion then take the form

$$\frac{d^2𝝆}{dt^2}+\frac{GM𝝆}{\rho ^3}=𝐅,$$
(22)
where $`𝐅`$ is the post-Newtonian perturbing function. The details of this process of reduction of the equations of motion to the form (20) are straightforward and have been described in detail in . Under certain circumstances, the orbital perturbations can be easily characterized in terms of a modified orbit as has been done for the gravitomagnetic clock effect . To illustrate this procedure, it is instructive to digress here and consider the case of bound motion in the Kerr-Taub-NUT spacetime; the results can be compared and contrasted with those of the Kerr-de Sitter spacetime. Linearizing this metric in the angular momentum and the NUT parameter of the source, we find the Schwarzschild metric together with the Lense-Thirring and Taub-NUT off-diagonal terms. Introducing the isotropic Schwarzschild radial coordinate $`\rho `$ and the corresponding Cartesian coordinates, we find that this gravitoelectric background is perturbed by a gravitomagnetic field,

$$𝐁_g=\frac{GJ}{c\rho ^5}\left[3(𝝆\widehat{𝐉})𝝆\rho ^2\widehat{𝐉}\right]\frac{c^2\mathrm{}𝝆}{\rho ^3}$$
(23)
due to the angular momentum (i.e. the gravitomagnetic dipole moment J) and the gravitomagnetic monopole moment $`(c^2\mathrm{})`$ of the source. Here the NUT parameter $`\mathrm{}`$ has the dimension of length. The post-Newtonian gravitomagnetic orbital perturbation can then be obtained from equation (20), where $`𝐅`$ is of the form $`𝐅=2𝐯\times 𝐁_𝐠/c`$. The linear perturbation of the orbit due to the angular momentum of the source has been worked out in detail in  in connection with the gravitomagnetic clock effect. Therefore, we only indicate here the orbital perturbations due to $`𝐅_{\mathrm{}}=2c\mathrm{}𝝆\times 𝐯/\rho ^3`$ using the method developed in .

For the motion of a test mass around a gravitomagnetic monopole, the energy of the particle is conserved and so is its angular momentum, which is however augmented by a contribution from the interaction of the particle with the monopole, i.e. $`𝐣=m𝝆\times 𝐯2mc\mathrm{}𝝆/\rho `$ is the conserved quantity. Following the analysis given in , let us note that the unperturbed orbit is given by a Keplerian ellipse in the $`(X,Y)`$-plane of the $`(X,Y,Z)`$ coordinate system that is related to the $`(x,y,z)`$ system by a rotation,

$`x`$ $`=`$ $`\mathrm{cos}\mathrm{\Omega }X\mathrm{sin}\mathrm{\Omega }\mathrm{cos}iY+\mathrm{sin}\mathrm{\Omega }\mathrm{sin}iZ,`$ (24)
$`y`$ $`=`$ $`\mathrm{sin}\mathrm{\Omega }X+\mathrm{cos}\mathrm{\Omega }\mathrm{cos}iY\mathrm{cos}\mathrm{\Omega }\mathrm{sin}iZ,`$ (25)
$`z`$ $`=`$ $`\mathrm{sin}iY+\mathrm{cos}iZ,`$ (26)
where $`i`$ is the inclination angle and $`\mathrm{\Omega }`$ is the longitude of the line of the ascending nodes. The Keplerian ellipse is then given by $`X=\rho \mathrm{cos}\phi ,Y=\rho \mathrm{sin}\phi `$ and $`Z=0`$, where

$$\rho =\frac{a(1e^2)}{1+e\mathrm{cos}(\phi g)},\frac{d\phi }{dt}=\frac{L_0}{\rho ^2}$$
(27)
and $`L_0=\sqrt{GMa(1e^2)}`$ is the specific orbital angular momentum of the unperturbed orbit. Here $`a,e`$ and $`g`$ are respectively the semimajor axis, eccentricity and argument of the pericenter of the unperturbed elliptical orbit.

Once the gravitomagnetic monopole moment of the source is turned on at $`t=0`$, when observations begin, the only change that occurs in the orbit to linear order is that $`Z`$ is no longer zero but is given by

$$Z=\frac{2}{c}\mu L_0\frac{1\mathrm{cos}(\phi \phi _0)}{1+e\mathrm{cos}(\phi g)},$$
(28)
where $`\mu =c^2\mathrm{}/(GM)`$ is the dimensionless strength of the monopole and $`\phi =\phi _0`$ at $`t=0`$. A detailed examination shows that there is no gravitomagnetic clock effect in this case to first order in $`\mathrm{}`$. In this linear order, the net angular momentum $`𝐣`$ turns out to have componets $`j_X=2mc\mathrm{}\mathrm{cos}\phi _0,j_Y=2mc\mathrm{}\mathrm{sin}\phi _0`$ and $`j_Z=mL_0`$.

The gravitomagnetic monopole moment of the source presumably exists for all time; therefore, the procedure described above has to be re-interpreted in terms of an osculating ellipse. Namely, the position and velocity of the test mass at any given instant of time define an elliptical orbit that is momentarily tangent to the actual perturbed orbit. Therefore, the perturbed motion may be described in terms of the evolution of the osculating ellipse as in the Lagrange planetary equations. In the case under consideration here the observations are assumed to begin at $`t=0`$; hence, the elliptical orbit (25) is simply the osculating ellipse of the perturbed orbit at $`t=0`$.

Let us now return to the Kerr-de Sitter metric and linearize it in the small parameters $`ϵ=a/M`$ and $`\delta =\mathrm{\Lambda }M^2`$ with $`M0`$. We then need to introduce an isotropic radial coordinate $`\rho `$ that puts the Kottler metric in isotropic form. The connection between $`r`$ and $`\rho `$ is given by the differential equation

$$\left(\frac{dr}{d\rho }\right)^2=\frac{r^2}{\rho ^2}\left(1\frac{2M}{r}\frac{\mathrm{\Lambda }}{3}r^2\right).$$
(29)
To first order in $`\delta `$, we find that

$$r(\rho )=\rho \left(1+\frac{M}{2\rho }\right)^2\frac{\delta }{24}\rho \left(1\frac{M^2}{4\rho ^2}\right)I\left(\frac{2\rho }{M}\right),$$
(30)
where $`I(x)`$ is given by

$$I(x)=^x\frac{(1+u)^6}{(1u)^2}\frac{du}{u^3}.$$
(31)
This integral can be evaluated and the result is

$$2I(x)=K+x^2+16x+60\mathrm{ln}x\frac{128}{x1}\frac{16}{x}\frac{1}{x^2},$$
(32)
where $`K`$ is a dimensionless integration constant that we can set equal to zero for the sake of simplicity. A detailed examination of the resulting equation of motion then shows that the dominant effect of the cosmological constant would simply be due to the effective Newtonian perturbation of the form $`𝐅=\lambda 𝝆`$, where $`\lambda =\mathrm{\Lambda }c^2/3`$. Under such a perturbation the orbital angular momentum is conserved and the orbit is thus planar. On the other hand, the orbital energy is augmented by a contribution from the cosmological constant so that the conserved quantity is $`E`$, where $`2E=v^2\lambda \rho ^22GM/\rho `$. Though the perturbing function is simply proportional to $`𝝆`$ in this case, it turns out that the simple method developed in  is not directly applicable in this case for orbits of arbitrary eccentricity, since the perturbing function has the form of an infinite series in powers of the eccentricity (cf. section 5). This is in contrast with the Kerr-Taub-NUT case. It is therefore necessary to adopt a different approach to the motion of the perturbed orbit unless $`e<<1`$, in which case the method of  is adequate. The general solution to first order in the cosmological constant is considered in the next section; in the following, we study the average motion of the perturbed orbit.

It is interesting to compute the average rate of precession of the pericenter under the influence of $`\mathrm{\Lambda }`$. To this end, let us introduce the Runge-Lenz vector of the perturbed orbit

$$𝐐=𝐯\times 𝐋GM𝝆/\rho ,$$
(33)
where $`𝐋=L_0\widehat{𝐙}`$ and we note that for $`𝐅=0,𝐐_\mathrm{𝟎}=GMe(\mathrm{cos}g,\mathrm{sin}g,0)`$ in the $`(X,Y,Z)`$ coordinate system. For $`𝐅=\lambda 𝝆`$, the orbit remains in the $`(X,Y)`$-plane and $`\dot{Q}_X=\lambda L_0Y,\dot{Q}_Y=\lambda L_0X`$ and $`\dot{Q}_Z=0`$. Thus in this case

$$\frac{d𝐐}{dt}=\lambda 𝐋\times 𝝆.$$
(34)
Averaging this relation over the “fast” motion, i.e. the unperturbed orbital motion of the test particle (with period $`T`$) would reveal the “slow” (secular) behavior of the Runge-Lenz vector over a long period of time. Defining the average over a period as

$$<f>=\frac{1}{T}_0^Tf𝑑t=(1e^2)^{3/2}\frac{1}{2\pi }_0^{2\pi }\frac{f(\phi )d\phi }{[1+e\mathrm{cos}(\phi g)]^2},$$
(35)
we obtain

$$<X>=\frac{3}{2}ae\mathrm{cos}g,<Y>=\frac{3}{2}ae\mathrm{sin}g,$$
(36)
where we have used the relation

$$\frac{1}{2\pi }_{\zeta _0}^{\zeta _0+2\pi }\frac{\mathrm{cos}\zeta d\zeta }{(1+e\mathrm{cos}\zeta )^3}=\frac{3}{2}\frac{e}{(1e^2)^{5/2}}.$$
(37)
Thus we find that

$$<\frac{d𝐐}{dt}>=𝛀\times 𝐐_\mathrm{𝟎},$$
(38)
where

$$𝛀=\frac{3}{2}\frac{a\lambda L_0}{GM}\widehat{𝐙}=\frac{1}{2}\frac{\mathrm{\Lambda }c^2}{\omega }\sqrt{1e^2}\widehat{𝐙}.$$
(39)
Here $`\omega =2\pi /T`$ is the Keplerian frequency of the unperturbed orbit. For $`\mathrm{\Lambda }10^{56}\mathrm{cm}^2`$, the precession frequency (37) is too small to be observable at present in the solar system . The precession angle after one period is given by

$$\frac{\pi c^2\mathrm{\Lambda }a^3}{GM}\sqrt{1e^2}.$$
(40)
To first order in eccentricity, this result agrees with previous work , while the exact dependence on eccentricity as $`\sqrt{1e^2}`$ differs from $`(1e^2)^3`$ given in .

To gain a deeper understanding of the average behavior of the orbit, we consider the osculating ellipse in this case and note that the orbital elements of this ellipse vary with time according to the Lagrange planetary equations  as

$`{\displaystyle \frac{da}{dt}}`$ $`=`$ $`{\displaystyle \frac{2\lambda ae\sqrt{1e^2}}{\omega }}{\displaystyle \frac{\mathrm{sin}(\phi g)}{1+e\mathrm{cos}(\phi g)}},`$ (41)
$`{\displaystyle \frac{de}{dt}}`$ $`=`$ $`{\displaystyle \frac{1e^2}{2ea}}{\displaystyle \frac{da}{dt}},`$ (42)
$`{\displaystyle \frac{dg}{dt}}`$ $`=`$ $`{\displaystyle \frac{\lambda \sqrt{(1e^2)^3}}{e\omega }}{\displaystyle \frac{\mathrm{cos}(\phi g)}{1+e\mathrm{cos}(\phi g)}}.`$ (43)
Employing the averaging method once again, we find that

$$<\frac{da}{dt}>=0,<\frac{de}{dt}>=0,$$
(44)
$$<\frac{dg}{dt}>=\frac{3}{2}\frac{\lambda }{\omega }\sqrt{1e^2},$$
(45)
where relation (35) has been used. The rate of advance of the pericenter given by equation (43) is consistent with our previous result. Thus in the presence of a cosmological constant the orbital plane does not change and the orbit keeps its shape on the average while precessing with frequency $`\mathrm{\Lambda }c^2\sqrt{1e^2}/(2\omega )`$.

It follows from equation (37) that the frequency of the pericenter precession caused by the cosmological constant is directly proportional to the period of the binary orbit, which must therefore be very large for the effect to be perceptible. The ratio of this precession frequency to the Einstein pericenter precession frequency is given by $`\mathrm{\Lambda }a^4c^4(1e^2)^{3/2}/(6G^2M^2)`$, which is $`10^{14}`$ for the motion of the Earth around the Sun. We conclude that the influence of the cosmological constant on orbital motion is too small to be measurable in the solar system.

5. Perturbed Keplerian motion

It is interesting to give the general solution of equation (20) for $`𝐅=\lambda 𝝆`$ to first order in the cosmological constant. For a positive cosmological constant, the effective force on the test mass is along the radial direction away from the source. The orbit is thus planar and $`\rho ^2d\phi /dt=L_0`$ is the constant specific orbital angular momentum. In terms of $`u=\rho ^1`$, the radial part of equation (20) reduces in this case to

$$\frac{d^2u}{d\phi ^2}+u=p^1\frac{\lambda }{L_0^2u^3},$$
(46)
where $`p=L_0^2/(GM)`$ is a positive constant. The right side of equation (44) can be expressed as a function of $`\phi `$ if we substitute for $`u`$ the unperturbed solution in accordance with our linear perturbation scheme. The unperturbed solution is given by

$$u_0=\frac{1}{p}\left(1+e\mathrm{cos}\phi \right),$$
(47)
where we have chosen the planar Cartesian coordinates such that $`g=0`$ for the sake of simplicity. Here $`e0`$ is the orbital eccentricity of the conic section. For an ellipse $`e<1`$ and $`p=a(1e^2)`$, where the pericenter distance from the focus is $`\rho _{min}=a(1e)`$ and $`a`$ is the semimajor axis. For a parabola $`e=1`$ and $`\rho _{min}=p/2`$ is the pericenter distance, while for a hyperbola $`e>1`$ and $`p=a(e^21)`$, where the pericenter distance from the focus is $`\rho _{min}=a(e1)`$. Substituting equation (45) for $`u`$ in the right side of equation (44), the resulting equation can be transformed to the form

$$(1\xi ^2)\frac{d^2U}{d\xi ^2}\xi \frac{dU}{d\xi }+U=\frac{q}{(1+e\xi )^3},$$
(48)
where $`\xi =\mathrm{cos}\phi ,U=up^1`$ and $`q=\lambda p^3/L_0^2`$. Equation (46) can be integrated once and we find

$$(1\xi ^2)\frac{dU}{d\xi }+\xi U=\frac{q}{2e(1+e\xi )^2}+C,$$
(49)
where $`C`$ is an integration constant. Equation (47) can be written in the form

$$\frac{d}{d\xi }\left[(1\xi ^2)^{1/2}U\right]=\frac{q}{2e(1\xi ^2)^{3/2}(1+e\xi )^2}+\frac{C}{(1\xi ^2)^{3/2}}$$
(50)
and integrated using formulas 2.264, 2.266 and 2.268 given in . The result for an elliptical orbit is

$$U=\frac{1}{2}\frac{q}{(1e^2)^2}\left[3(1+e\mathrm{sin}\phi )\frac{1e^2}{1+e\mathrm{cos}\phi }(2e+\frac{1}{e})\mathrm{cos}\phi \right]+C\mathrm{cos}\phi +S\mathrm{sin}\phi ,$$
(51)
where $`S`$ is a constant of integration and

$$=\frac{1}{\sqrt{1e^2}}\mathrm{arc}\mathrm{sin}\left(\frac{e+\mathrm{cos}\phi }{1+e\mathrm{cos}\phi }\right).$$
(52)
The integration of equation (48) can be similarly carried out for parabolic and hyperbolic orbits using the relevant formulas of . We expect that $`u`$ reduces to the unperturbed solution $`u_0`$ as $`\lambda 0`$; therefore, it is useful to set $`C=qC^{}+e/p`$ and $`S=qS^{}`$. The constants $`C^{}`$ and $`S^{}`$ can be uniquely determined using the initial conditions, i.e. the initial position and velocity of the test mass. In this way, we find the orbital radius $`\rho =u^1=(p^1+U)^1`$ and then $`\phi (t)`$ by integrating $`d\phi /dt=L_0\rho ^2`$.

6. Discussion

In this paper we have considered some of the local physical consequences of the existence of a cosmological constant. In particular, we have considered the motion of test bodies and the propagation of rays in the Kerr-de Sitter spacetime. While the clock effect is rather weakly affected by the presence of a cosmological constant, there is a (gravitoelectric) time delay in the propagation of rays directly proportional to $`\mathrm{\Lambda }`$. Though this effect is too small to be measurable in the solar system, it should in general be taken into account in the interpretation of gravitational lensing time delay. Moreover, the average behavior of bounded orbits has been examined and the general rate of orbital precession has been determined for an arbitrary eccentric orbit. On the average, the orbital plane remains invariant and the orbit keeps its shape but precesses due to the presence of a cosmological constant; however, the precession rate is too small to be detectable in the solar system.

References

 Cohen JM and Mashhoon B 1993 Phys. Lett. A 181 353

 Mashhoon B, Gronwald F and Theiss DS 1999 Ann. Phys. (Leipzig) 8 135

 Mashhoon B, Gronwald F and Lichtenegger HIM 2000 in: Gyros, Clocks, Interferometers…: Testing Relativistic Gravity in Space, edited by C. Lämmerzahl, C.W.F. Everitt and F.W. Hehl (Berlin : Springer), pp. 83-108

 Mashhoon B, Iorio L and Lichtenegger H 2001 Phys. Lett. A 292 49

 Demiański M 1973 Acta Astron. 23 197

 Carter B 1973 in: Black Holes, edited by C. DeWitt and B.S. DeWitt (New York: Gordon and Breach), pp. 101-103

 Kottler F 1918 Ann. Phys. (Leipzig) 56 410

 Schmidt H-J 1993 Fortschr. Phys. 41 179

 Liu H and Mashhoon B 1995 Ann. Phys. (Leipzig) 4 565

 Peebles PJE and Ratra B 2003 Rev. Mod. Phys., in press; astro-ph/0207347

 Howes RJ 1981 Gen. Rel. Grav. 13 829

 Bini D, Jantzen RT and Mashhoon B 2001 Class. Quantum Grav. 18 653

 Bini D, Jantzen RT and Mashhoon B 2002 Class. Quantum Grav. 19 17

 Maartens R, Mashhoon B and Matravers DR 2002 Class. Quantum Grav. 19 195

 Ciufolini I, Kopeikin S, Mashhoon B and Ricci F 2003 Phys. Lett. A 308 101

 Islam JN 1983 Phys. Lett. A 97 239

 Lake K 2002 Phys. Rev. D 65 087301

 Mashhoon B 2001 in: Reference Frames and Gravitomagnetism, edited by J.-F. Pascual-Sánchez, L. Floría, A. San Miguel and F. Vicente (Singapore: World Scientific), pp. 121-132

 Rindler W 2001 Relativity: Special, General, and Cosmological (Oxford: Oxford University Press), p. 305

 Danby JMA 1988 Fundamentals of Celestial Mechanics, 2nd ed. (Richmond: Willmann-Bell)

 Gradshteyn IS and Rhyzhik M 1980 Table of Integrals, Series and Products (New York: Academic Press)


text:
# EGRET upper limits on the high-energy gamma-ray emission of galaxy clusters

## 1 Introduction

Clusters of galaxies are excellent representatives for the formation and the evolution of structure in the universe. They have been extensively studied at radio, optical and X-ray wavelengths. Within the last decade, radio, extreme UV and hard X-ray observations have revealed emission features that led to the prediction that galaxy clusters might be emitters of high-energy gamma-rays:

* the existence of diffuse radio halos (Giovannini et al., 1993, 1999; Giovannini & Feretti, 2000; Kempner & Sarazin, 2001),
* the rather controversially discussed observations of EUV excess emission in galaxy clusters like A1795, A2199, and the Coma Cluster (Bowyer at al., 1999), Abell 2199 (Lieu at al., 1999), A1367 and A1656 (Coma), A1795 and A2199 (Arabadjis & Bregman, 1999), Virgo (Berghöfer at al., 2000), Virgo and A1795 (Bonamente et al., 2001), the Fornax Cluster (Bowyer at al., 2001), A2199 and A1795 (Berghöfer & Bowyer, 2002), A1795, A2199, A4059, Coma and Virgo (Dureet et al., 2002),
* the observational hint of a distinct non–thermal emission component at hard X-ray wavelengths in the case of the Coma cluster (Fusco–Femiano et al., 1999; Rephaeli et al., 1999), Abell 2199 (Kaastra et al., 1999), Abell 2256 (Fusco–Femiano et al., 2000), and perhaps A754, A119 (Fusco–Femiano et al., 2002).

Various scenarios were suggested to connect and explain the links between these observations and, consequently, to predict a high-energy emission component at gamma-ray wavelengths. Whereas the diffuse radio emission is clearly synchrotron radiation by highly relativistic electrons, the EUV excess emission was first attributed to a second but cooler thermal component. Now a more plausible explanation is Inverse Compton scattering of Cosmic Microwave Background radiation by a non–thermal electron population (Enßlin & Biermann, 1998; Blasi & Colafrancesco, 1999). The hard X-ray excess can be produced by Inverse Compton scattering of the same electron distribution generating the non–thermal radio emission (Giovannini et al., 1993). To avoid the problem of the rather low magnetic field strength in such a scenario, non–thermal bremsstrahlung has been proposed as an alternative emission process (Enßlin et al., 1999). As pointed out by Petrosian (2001), the non–thermal bremsstrahlung cannot be persistently produced on account of the low radiation efficiency of electrons in the 100 keV range. Hadronic particle populations were considered to produce gamma-rays via pp–interactions of high-energy cosmic rays with the intracluster medium (ICM) (Berezinsky et al., 1997), or as the origin of a secondary population of relativistic electrons Atoyan & Völk (2000). Cluster merger systems might offer sufficient cosmic ray injection rates in conjunction with a mechanism for heating the ICM to the observed temperatures (Blasi, 2001; Fujita & Sarazin, 2002).
Gamma-ray radiation from galaxy clusters is also expected as a result of large scale cosmological structure formation scenarios (Dar & Shaviv, 1995, 1996; Colafrancesco & Blasi, 1998; Waxman & Loeb, 2000; Totani & Kitayama, 2000; Kawasaki & Totani, 2002; Miniati, 2002). However, apart from the general prediction of its existence, quantitative estimates range between ’dominant part of the already observed extragalactic diffuse background by EGRET’ to ’magnitudes below the detection threshold of the current gamma-ray instrumentation’ - a range of predictions substantially more uncertain than that for the contribution of unresolved AGN to the extragalactic diffuse gamma-ray background (see i.e. Mücke & Pohl (2000), and references therein). The benefit of dealing with a class of astronomical objects already detected at gamma-ray wavelengths as e.g. AGN is not granted for the galaxy clusters: In contrast to the blazar population well-observed by EGRET, no galaxy cluster has been unambiguously identified at gamma-ray wavelengths to date. Nevertheless, for several individual clusters model predictions exist, which place their gamma-ray fluxes close to or even below the instrumental sensitivity threshold of the EGRET telescope at E $`>100`$ MeV (Dar & Shaviv, 1995; Enßlin et al., 1997; Blasi & Colafrancesco, 1999). Until now, galaxy clusters in gamma-rays have only been analyzed using early EGRET data and preliminary analysis techniques, resulting in non–detections of the Coma cluster (Sreekumar et al., 1996) and several Abell clusters (McGlynn et al., 1994). Therefore galaxy clusters have not been considered as likely counterparts of EGRET sources in the 3EG source catalog (Hartman et al., 1999).
Just recently, claims of an association between galaxy clusters from the Abell catalog and unidentified gamma-ray point sources from the 3EG catalog have been made by Colafrancesco (2001) and Kawasaki & Totani (2002). Likewise, Abell clusters were proposed to be connected with unresolved gamma-ray excesses (Scharf & Mukherjee, 2002). All these detection claims have a statistical significance for association at the 3 $`\sigma `$ level in common.
Here, in order to provide an up–to–date and comprehensive view on the high-energy gamma-ray emission from galaxy clusters, we have expanded a preliminary analysis by Reimer et al. (1999) by considering all relevant EGRET observations between 1991 and 2000. Using the finalized EGRET data, which incorporate the latest instrumental efficiency normalizations, we analyzed individual, nearby X-ray bright galaxy clusters with the likelihood technique. Subsequently, the gamma-ray data from individual galaxy clusters have been coadded in cluster–centered coordinates (Reimer & Sreekumar, 2001). The coadded images were again analyzed using the likelihood technique, however in conjunction with an adequately adapted diffuse gamma-ray foreground model. We also re–examined the statistical associations between unidentified EGRET sources and Abell clusters as a population (Colafrancesco, 2002). For that purpose we measured the cluster autocorrelation and thus derived the correct chance probabilities for the null hypothesis of no correlation between EGRET sources and Abell clusters. Finally, we compare our results, which benefit from the application of the likelihood analysis technique, with the result by Scharf & Mukherjee (2002).

## 2 The flux-limited X-ray bright galaxy cluster sample

For observationally probing the gamma-ray emission of galaxy clusters a sample of X-ray emitting clusters of galaxies has been chosen. This sample consists of the X-ray flux limited cluster catalogs from EINSTEIN (Edge et al., 1990), EXOSAT (Edge & Steward, 1991), and ROSAT surveys (XBACs: Ebeling et al. (1996), BCS north: Ebeling et al. (1998), BCS south: De Grandi et al. (1999)). Cluster selections based on X-ray catalogs currently provide the best way to obtain completeness without introducing biases (i.e. projection effects). Although appearing as extended sources with typical radii of several arcminutes in X-rays, the width of the point spread function of the EGRET instrument ($`\theta =5.85^{}(E_\gamma /100\mathrm{M}\mathrm{e}\mathrm{V})^{0.534}`$, and $`\theta `$ the energy dependent radius for a 68% flux enclosure) does not permit a similar handling of galaxy clusters as extended sources in gamma-rays. Thus, the attempt to analyze clusters of galaxies as point–like excesses at energies above 100 MeV is justified. Here, a total of 58 individual X-ray bright galaxy clusters within z$`<0.14`$ were chosen to represent a reasonable candidate sample for the subsequent analysis at high-energy gamma-ray wavelengths. Although further but similar cluster surveys are on the way or have been completed recently (in particular HIFLUGCS and REFLEX), this sample adequately represents the high-flux end of the log N–log S distribution of X-ray bright galaxy clusters. Almost all clusters that are extensively discussed in the literature for evidence of non–thermal X-ray emission, EUV-excess features and/or characteristic diffuse radio halos are included in this sample. The number of galaxy clusters has been restricted primarily to achieve a manageable amount of analysis work in the gamma-rays, but also because a simple measure like cluster mass M over distance squared D2, as explained in the discussion, should already be major constraint for the detectability of galaxy clusters in gamma-rays. Thus only nearby clusters (z$`<0.14`$) were considered. This choice reflects the expectation that the nearest, most massive galaxy clusters are most likely the ones to be detected as individual sources of gamma-ray emission, whatever their flux. Figure 1 shows the spatial arrangement of the galaxy cluster sample ($`\times `$) and cataloged high-energy gamma-ray sources (Hartman et al., 1999), in galactic coordinates.

EDITOR place figure 1 here

## 3 Analysis of galaxy clusters at high-energy gamma-rays

### 3.1 The study of individual galaxy clusters

Until recently, no positional coincidences between an individual galaxy cluster and gamma-ray point sources in existing EGRET source catalogues have been reported. For the Coma cluster the result of an EGRET analysis has been published, based on observations from CGRO cycle 1 and 2 (Sreekumar et al., 1996). In the analysis described here, EGRET data of individual viewing periods from CGRO observation cycles 1 to 9 were used for the analysis of 58 individual clusters. The latest and presumably final improvements in the efficiency corrections for the instrumental response of the EGRET spark chamber telescope have been fully implemented. Each galaxy cluster has been individually analyzed by means of standard EGRET data reduction techniques (likelihood source finding algorithm and subsequent flux determination at the position of the center of the X-ray emission). This analysis goes beyond the preliminary study presented by Reimer et al. (1999), in which four years of EGRET observations were analyzed in strict congruence with the four years of EGRET observations used for the 3EG catalog of gamma-ray point sources. Coadded images of individual viewing periods, where a cluster has been observed at less the 30° off the pointing axis of the EGRET instrument (standard field-of-view observations) or less then 19° (narrow field-of-view observations), have been searched for gamma-ray excesses after modeling cataloged (and therefore well-known) identified gamma-ray point sources by using the maximum-likelihood technique as described in Mattox et al. (1996). Gamma-ray source fluxes have been determined at the coordinates of the cluster center position known from X-ray observations. Applying the same detection criteria as used and described in the EGRET source catalogs, none of the 58 galaxy clusters are detected in the EGRET data. Special care has been exercised when already cataloged gamma-ray point sources are near the position of an Abell cluster, see footnote of Table 1. Three of these sources are unambiguously identified blazar-class Active Galactic Nuclei: 1633+382 = 3EG J1635+3818 near A2199, 3C279 = 3EG J1255-0549 near A1651, and 1604+159 = 3EG J1605+1553 near A2147. In these cases the particular AGN has been modelled at it’s known radio position and simultaneously taken into account in the deteremination of the gamma-ray flux at the position of the Abell cluster in question. Only one of the remaining three catalog sources shows considerable overlap at the position of an analyzed cluster (A85 with the unidentified source 3EG J0038-0949). Keeping in mind the width of the point spread function of the EGRET telescope, the total number of unidentified gamma-ray sources, and the size of our galaxy cluster sample, this occurrence is perfectly in agreement with chance coincidence: The probability of at least one coincidence between one of the 170 unidentified EGRET sources and one of the 58 considered galaxy clusters at a distance of 0.81° as in case of 3EG J0038-0949 and A85, is 48.1 %. Here, a gamma-ray flux has been obtained at the cluster position without explicitly modeling the known, but unidentified gamma-ray source.
Exemplary for the detailed results given in Tab.1, the intensity and test statistics map of three prominent galaxy clusters (Perseus, Virgo, and Coma) are shown in Fig.2.

EDITOR place figure 2 here

The most significant gamma-ray excess for any individual cluster has been found for A3532, corresponding to a detection significance of 1.6 $`\sigma `$. This is well below the detection threshold for being seriously considered as a source by standards of the EGRET data analysis, which is $`>4\sigma `$ at high-galactic latitudes. Furthermore the achieved detection significance is less than expected from statistical fluctuations at 58 trials alone. Therefore, for each galaxy cluster a 2 $`\sigma `$ upper limit at the position of the cluster center has been determined and is given in Tab.1

EDITOR place table 1 here

### 3.2 The study of the cluster population

Having established that individual clusters are not found in the EGRET data, a further analysis has been performed to study whether or not these galaxy clusters radiate in gamma-rays as a population. For this purpose, the counts, exposure, and intensity maps of the individual galaxy clusters have been used. Each individual map has been transformed into a cluster–centered coordinate system under conservation of the original pixelation. Subsequently, the individual images have been coadded. Three sets of images have been produced in order to assure that the center region of the stacked image is not dominated by already identified point sources or the galactic plane:
(1) a superposition of all 58 galaxy clusters in the sample,
(2) a superposition of 54 galaxy clusters excluding those with unambiguously identified and dominant EGRET sources in the central map bins (A2199: 1633+382, A1650, A1651, and A1689: 3C279), and
(3) a superposition of 50 galaxy clusters excluding those with unambiguously identified and dominant EGRET sources (see above) or the Galactic Plane in the map center (the Oph and Cyg A cluster, 3C129, and 3A 0745-191).

Each of the three sets have been analyzed, however, here we only report from the least contaminated set (3). Due to the wealth of accumulated instrumental exposure the results do not differ dramatically between these selections. The interpretation, however, simplifies considerably since we are not urged to discuss identified and dominant point sources or foreground emission features at or close to the image center anymore. The total exposure in the center region, averaged over the four central 0.5° $`\times `$ 0.5° map bins, is $`3.4\times 10^{10}`$ cm2 s (E $`>100`$ MeV) for the 50 cluster selection, the lowest values at the edge of the 40° by 40° images are about $`1.4\times 10^{10}`$ cm2 s. Figure 3 (upper part) shows the coadded counts, exposure and intensity map for the 50 cluster sample.

EDITOR place figure 3 here

We have analyzed these images using the maximum-likelihood technique as already described. Here, however, we have to provide a customized diffuse galactic gamma-ray foreground model, adapted for the application of our superpositioned cluster sample in the cluster–centered coordinate system. This was achieved by adopting the standard diffuse galactic emission model (GALDIF) (Hunter et al., 1997) used for EGRET likelihood analysis, given on a 0.5° by 0.5° grid, into a specific diffuse foreground model for the galaxy clusters in cluster–centered coordinates: Corresponding in image size and coordinates with the counts, exposure and intensity map of the individual clusters, the appropriate diffuse maps have been taken directly from the GALDIF model. These maps were subsequently transformed into the cluster–centered coordinate system.

In each individual cluster observation we expect the diffuse foreground emission to contribute a certain number of counts, $`c_i`$, that can be calculated as the product of the respective exposure map, $`\epsilon _i`$, and the intensity of the diffuse emission, $`DF_i`$. The appropriate diffuse foreground intensity for the coadded data, $`DF_{tot}`$, is then derived by:

$$DF_{tot}=\frac{1}{\epsilon _{tot}}\underset{i}{}c_i=\frac{1}{\epsilon _{tot}}\underset{i}{}\epsilon _iDF_i,\epsilon _{tot}=\underset{i}{}\epsilon _i.$$

By applying the maximum likelihood procedure in conjunction with the appropriately adapted diffuse model, the 40° by 40° images were searched for excesses. Of interest here is only the map center, corresponding to the emission maximum of the considered galaxy clusters in X-rays. No significant gamma-ray emission excess has been found within a radius of 5° of the origin in the cluster-centered image. With a cumulative exposure of 3.4 $`\times `$ $`10^{10}`$ cm2 s for E $`>`$ 100 MeV the corresponding upper limit is 5.9 $`\times `$ $`10^9`$ cm-2 s-1 (averaged over the four central 0.5° $`\times `$ 0.5° map bins) for the so–constructed average galaxy cluster. Fig.3 (lower part) shows the customized diffuse model and the resulting likelihood test statistics maps. Easily seen in the test statistics image is 3C279, located at $``$ 13° from the map center, and thus far too distant to cause any conflict with the determined upper limit.

## 4 Discussion

### 4.1 Cases of individual galaxy clusters

The negative results from both an analysis of the gamma-ray data from the EGRET instrument at positions of 58 individual galaxy clusters as well as from a superposition of 50 galaxy clusters needs to be critically discussed with respect to underlying systematics. Categorically, the question of an appropriately chosen selection of galaxy clusters might arise. The assumption has been made that the brightest and nearest clusters detected at X-ray wavelengths should be the most likely candidates to emit observable gamma-rays, supported by various models explaining the multifrequency emission properties and the general understanding of confinement and interaction of cosmic rays in the intercluster medium (ICM) of a galaxy clusters (Berezinsky et al., 1997; Völk et al., 1996). Because almost all clusters exhibiting unusual multifrequency emission characteristics (EUV-excess emission, non–thermal X-ray emission and/or a diffuse radio halo) are naturally included here, the above assumption is certainly not artificial. In Tab.2 we compare our results with two different scenarios of gamma-ray emission from cosmic ray interactions in the ICM, quantitatively predicting gamma-ray emission for some individual galaxy clusters.

EDITOR place table 2 here

Our upper limits are consistently below the predictions as given by Enßlin et al. (1997), and especially Dar & Shaviv (1995). Thus, the suggested scenarios are ruled out in the given parameter space. Concerns against the results of Dar & Shaviv (1995) have been already pointed out earlier by Stecker & Salamon (1996) on account of the spectra of the diffuse galactic and extragalactic gamma-ray background. With the apparent conflict between these model predictions and the observational upper limits for the individual galaxy clusters, we clearly disfavor these models against models which predict gamma-ray emission below the sensitivity of the instruments of the Compton gamma-ray observatory era like inverse Compton scenarios from cosmic ray electrons accelerated at accretion shocks by (Colafrancesco & Blasi, 1998; Miniati, 2002).
At present, with quantitative predictions about the gamma-ray emission made only in cases of the most prominent galaxy clusters, and with measured upper limits and predicted fluxes not orders of magnitudes apart, only a moderate conservative interpretation is tenable: The degree of freedom in the parameters in the models in conflict with our measurements, especially if predicting a diffuse extragalactic background component without a single positive detection of the object class in question made does not allow to discriminate between the suggested scenarios solely on the basis of the presented gamma-ray data. Not until the next generation of gamma-ray telescope will be available with good prospects to detect individual galaxy clusters as well as the chance to observationally constraint the constituents of the extragalactic diffuse background, the sensitive measurements to discriminate the various scenarios for gamma-ray emission from galaxy clusters must originate at other wavelengths. Here, especially three wavebands are of particular interest: the new generation of low-threshold ground–based Cherenkov–telescopes (IACTs) presently being built or just starting their observations will provide very sensitive measurements at energies above 100 GeV, thus providing the necessary information to constrain the inverse Compton component at the highest energies; more and better hard X-ray measurements may clearly discriminate between thermal and non–thermal emission components and therefore help to identify the nature of the originating particle population, in particular deciding on the preference of inverse Compton or non–thermal Bremsstrahlung scenarios; finally, more sensitive radio halo measurements may be performed, especially at the high–frequency end of the synchrotron emission component, where currently only the radio halo of the Coma cluster has been sufficiently investigated (Schlickeiser et al., 1987; Thierbach et al., 2003). High-frequency radio measurements from diffuse cluster halos will directly constrain the shape and intensity of the resulting inverse Compton component at high energies.

### 4.2 The case of galaxy clusters as a population

We now discuss the result for the average cluster from our analysis of 50 galaxy clusters population. At the achieved sensitivity there is still no indication of gamma-ray emission from galaxy clusters. The 2 $`\sigma `$ upper limit for the flux of the average cluster is 5.9 $`\times `$ $`10^9`$ cm-2 s-1. This may help to resolve the stark inconsistency between studies performing direct EGRET data analysis on galaxy clusters (i.e. McGlynn et al. (1994), Sreekumar et al. (1996), and this work) and the recently published detection claims originating from information on gamma-ray point sources from the 3EG catalog in conjunction with statistical assessments by Colafrancesco (2001), Colafrancesco (2002) and Kawasaki & Totani (2002). Here, we would like to give a reassessment of these detection claims on a firm spatial–statistical basis. The cluster sample studied by Colafrancesco consists of the entire Abell cluster catalog (Abell et al., 1989) at galactic latitudes $`|b|>20`$° . The corresponding Poissonian probability distribution for spatial association between the cluster sample and gamma-ray point sources form the EGRET catalog is easy to determine and given in Fig.4a. The two-point autocorrelation function, $`\omega _\mathrm{\Theta }`$, for the population of Abell clusters has been intensively studied previously (Hauser & Peebles, 1973; Postman et al., 1986; Olivier et al., 1990; Akylas et al., 2000), but its impact on the correlation analysis can be rigorously calculated only for a very large sample of objects. Here we are concerned with only 170 unidentified EGRET sources, 59 of which are located at $`|b|>20\mathrm{°}`$. Therefore, we have directly determined the chance probability for an association between an arbitrary source and one of the Abell clusters as a function of the radius–of–interest, i.e. the maximum separation of sources considered as associated. Colafrancesco (2002) reports an association of 50 EGRET sources (resp. 18 unidentified EGRET sources) and 70 Abell cluster (resp. 24 Abell clusters) based on a initial sample of 3979 Abell clusters and 128 EGRET sources (resp. 59 unidentified EGRET sources) at $`|b|>20\mathrm{°}`$. Pure Poissonian statistics predicts a total of 34.4 (17.2) single, 8.6 (2.0) double and 1.3 (0.2) triple associations between EGRET sources and Abell associations at 1° roi (resp. 0.81° as the average $`\mathrm{\Theta }_{95}`$). Using the modified chance probabilities (Fig. 4b) to account for the autocorrelation of the Abell clusters, one expects a total of 28.0 (10.6) single, 8.9 (2.5) double and 2.6 (0.5) triple association by chance. Autocorrelation among the unidentified EGRET sources at $``$ 1° scale is per se excluded due to the moderate angular resolution of the EGRET telescope and its relatively inability to discriminate neighboring sources adequately (”source confusion”). Hence the expected chance associations amount to 40.7 EGRET sources and 56.6 Abell clusters for 128 EGRET sources at 1° roi and 13.7 EGRET sources and 17.5 Abell clusters for 59 unidentified EGRET sources at $`\mathrm{\Theta }_{95}`$, respectively. In terms of the cumulative Poisson probability, the significance of the correlation claim by Colafrancesco needs to be reassessed to only 1.36 $`\sigma `$ in case of the 128 EGRET sources (which is anyway meaningless due to the contamination with already identified EGRET sources, the blazars), or 1.03 $`\sigma `$ in case of the 59 unidentified EGRET sources, rigorously indicating the statistical insignificance of the correlation claim.

EDITOR place figure 4 here

Furthermore, the explicitly suggested Abell cluster/unidentified EGRET source associations by Colafrancesco (2002) do not contain the most likely and prominent galaxy clusters as predicted to be the clusters with the best chance to be detected in high-energy gamma rays (Fig.5). Actually, only one of the 18 listed unidentified EGRET sources indeed has a truly X-ray bright Abell cluster counterpart; the remaining 23 clusters are not even included in the most recent and carefully compiled sample of a flux-limited X-ray bright galaxy clusters population (HIFLUGCS: Reiprich & Böhringer (2002)). Thus, the identification sequence as suggested by Colafrancesco (radio halo $``$ X-ray brightness $``$ counterpart in unidentified gamma-ray source) simply does not work this way, neither globally nor just in a few of the referenced cases. Similarly, the predictive power of the Lγ–Lradio– and Lγ–LX–correlation as given in Colafrancesco (2002) must be seriously questioned.

EDITOR place figure 5 here

A further suggestion of possibly merging clusters as counterparts of unidentified gamma-ray sources by Kawasaki & Totani (2002) is based on only seven individual gamma-ray sources, all of which belong to the sample of steady unidentified EGRET sources classified as such by Gehrels et al. (2000). However, more specific variability studies of EGRET sources as performed by Tompkins (1999) and Torres et al. (2001) do not support non-variability for each of these sources. Two sources clearly exhibit gamma-ray flux variability, and a further two sources belong to the group of gamma-ray sources where we could not decide on the basic of the current gamma-ray data whether or not a source is truly non–variable. Thus, a more conservative sample would consist of only three candidates, reducing the statistical significance of the remaining associations to the level of expected chance coincidences. Certainly, those newly-suggested counterparts are worthy of an intensive investigation in order to conclude about their true nature and, subsequently, about their chance to be related to an unidentified gamma-ray source. At present, this study does not leave sufficient confidence for an association of unidentified gamma-ray sources and cluster merger systems. Recently, the problem of such association has been theoretically addressed by Berrington & Dermer (2002), arguing that it is unlikely that more than a few of the isotropically distributed unidentified EGRET sources at high-galactic latitudes can be attributed to radiation from non–thermal particles produced by cluster merger shocks. Similar conclusions were reached by Gabici & Blasi (2003), emphasizing the discrepancy between non–thermal activity in galaxy clusters and the relative ineffectiveness of major shocks to energize the underlying particle population.

Yet another detection claim of galaxy clusters in the gamma-ray data was published recently. A population study by Scharf & Mukherjee (2002) claims a statistical detection at $`3\sigma `$ level from unresolved gamma-ray sources spatially coincident with Abell clusters, in particular those of high optical richness. It quotes a mean cluster flux of $``$ 1.14 $`\times `$ $`10^9`$ cm-2 s-1 at E $`>`$ 100 MeV in a 1° radius aperture, corresponding to a 68 % flux enclosure of the EGRET point-spread function. However, the energy–averaged EGRET flux enclosure at E $`>`$ 100 MeV in a 1° aperture, assuming a source spectrum with a power-law index of -2.0, is in fact much lower, namely 24.1 %. The 68 % flux enclosure is reached at a radius of 3.1° . Therefore the flux quoted for the innermost (1) radial bin in the study, $`1.19\times 10^6`$ ph s-1 cm-2 sr-1, in fact corresponds to a total source flux of 4.7 $`\times `$ $`10^9`$ cm-2 s-1, which is marginally inconsistent with our findings. Considering a 2° aperture by inclusion of the second innermost bin of $``$ 1.2 $`\times 10^6`$ ph s-1 cm-2 sr-1 from Scharf & Mukherjee (2002), the corresponding mean cluster flux will increase to 9.6 $`\times `$ $`10^9`$ cm-2 s-1, which is clearly inconsistent with our upper limit for the average cluster from the sample of nearby, X-ray bright galaxy clusters. Since it is not within the scope of this paper to perform an in–depth study of possible systematic effects, we only state that the results presented by Scharf & Mukherjee (2002), although at comparable sensitivity level, are apparently inconsistent with ours.

In conclusion, we still have to await the first observational evidence for the high-energy gamma-ray emission of galaxy clusters. The last generation of gamma-ray telescopes aboard CGRO was not able to resolve an individual galaxy cluster nor the nearby, X-ray brightest clusters of galaxies as a population. Until the next generation of gamma-ray instruments will challenge this important scientific topic, progress is expected at other wavelengths: from GHz-frequency radio observations of radio halos, from studies of soft and hard X-ray excess features with sufficient statistical significance, and from measurements of the new generation of imaging atmospheric Cherenkov telescopes.

O.R. acknowledges support from a National Academy of Science/National Research Council Associateship at NASA/Goddard Space Flight Center, where part of this study was executed. O.R. and M.P. are supported from DLR QV0002. O.R. would like to thank the many people he had consulted about the topic, in particular D.L. Bertsch, P. Blasi, S. Digel, F. Miniati, G. Madejski, V. Petrosian, D.J. Thompson, R. Schlickeiser and F. Stecker, for the general encouragement to continue this study beyond a conference contribution by reason of the rather unforeseen interest the topic experienced recently.