text:
# Can the dark energy equation-of-state parameter 𝑤 be less than -1?

## I Introduction

Cosmological observations strongly indicate that the universe is dominated by a smoothly distributed, slowly varying dark energy component (Riess:1998cb ; Perlmutter:1998np ; for reviews see Carroll:2000fy ; Peebles:2002gy .) The simplest candidate for such a source is vacuum energy, or the cosmological constant, characterized by a pressure equal in magnitude and opposite in sign to the energy density:

$$p_{\mathrm{vac}}=\rho _{\mathrm{vac}}.$$
(1)
While vacuum energy is strictly constant throughout space and time, it is also worthwhile to consider dynamical candidates for the dark energy. A convenient parameterization of the recent behavior of any such candidate comes from generalizing the vacuum-energy equation of state to

$$p=w\rho ,$$
(2)
which should be thought of as a phenomenological relation reflecting the current amount of pressure and energy density in the dark energy. In particular, the equation-of-state parameter $`w=p/\rho `$ is not necessarily constant. However, given that there are an uncountable number of conceivable behaviors for the dark energy, a simple relation such as (2) is a useful way to characterize its current state.

The equation-of-state parameter is connected directly to the evolution of the energy density, and thus to the expansion of the universe. From the conservation-of-energy equation for a component $`\rho _i`$ in a Robertson-Walker cosmology with scale factor $`a(t)`$ and Hubble parameter $`H=\dot{a}/a`$,

$$\dot{\rho }_i=3H(\rho _i+p_i),$$
(3)
it follows that this component evolves with the scale factor as

$$\frac{d\mathrm{ln}\rho _i}{d\mathrm{ln}a}=3(1+w_i).$$
(4)
We notice in particular that vacuum energy remains constant, while the energy density would actually increase as the universe expands if $`w_i<1`$. The Friedmann equations may be written as

$$H^2=\frac{8\pi G}{3}\rho \frac{\kappa }{a^2},$$
(5)
where $`\kappa `$ is the spatial curvature, and

$`{\displaystyle \frac{\ddot{a}}{a}}`$ $`=`$ $`{\displaystyle \frac{4\pi G}{3}}(\rho +3p)`$ (6)
$`=`$ $`{\displaystyle \frac{4\pi G}{3}}(1+3w)\rho .`$
From (6), we see that the universe will accelerate ($`\ddot{a}>0`$) if $`w<1/3`$. (Of course this is the effective $`w`$ of all the energy in the universe; if there is a combination of matter and dark energy, the dark energy will have to have a more negative $`w`$ in order to cause acceleration.) From (5) we see that a flat universe dominated by a component with constant $`w`$ will expand as

$$at^{2/3(1+w)},$$
(7)
unless $`w=1`$, for which the expansion will be exponential. (For $`w<1`$, one should choose $`t<0`$ in this expression.)

What are the possible values $`w`$ may take? It is hard to make sweeping statements about a component of energy about which we know so little. In general relativity, it is conventional to restrict the possible energy-momentum tensors by imposing “energy conditions”. In Garnavich:1998th  it was suggested that a reasonable constraint to impose would be the null dominant energy condition, or NDEC (see Section II for discussion). The physical motivation for a condition such as the NDEC is to prevent instability of the vacuum or propagation of energy outside the light cone. Applied to an equation of state of the form (2), the NDEC implies $`w1`$. Thus, pure vacuum energy is a limiting case; any other allowed component would diminish in energy as the universe expands.

Given our ignorance about the nature of the dark energy, it is worth asking whether this mysterious substance might actually confound our expectations for a well-behaved energy source by violating the NDEC. Given that the dark energy should have positive energy density (to account for the necessary density to make the universe flat) and negative pressure (to explain the acceleration observed in the supernova data), such a violation would imply $`w<1`$. It has been known for some time that such energy components can occur Nilles:1983ge ; Barrow:yc ; Pollock:xe . Their role as possible dark energy candidates was raised by Caldwell Caldwell:1999ew , who referred to NDEC-violating sources as “phantom” components, and has been since investigated by several authors (for some examples see Sahni:1999gb ; Parker:1999td ; Chiba:1999ka ; Boisseau:2000pr ; Schulz:2001yx ; Faraoni:2001tq ; maor2 ; Onemli:2002hr ; Torres:2002pe ; Frampton:2002tu ). Observational limits on $`w`$ Hannestad:2002ur ; Melchiorri:2002ux  are conventionally expressed as allowed regions in the $`w`$-$`\mathrm{\Omega }_\mathrm{M}`$ plane, assuming a flat universe ($`\kappa =0`$, or $`\mathrm{\Omega }_\mathrm{M}+\mathrm{\Omega }_X=1`$, where $`X`$ stands for the dark energy). Current limits Melchiorri:2002ux , obtained by combining results from cosmic microwave background experiments with large scale structure data, the Hubble parameter measurement from the Hubble Space Telescope and luminosity measurements of Type Ia supernovae, give $`1.62<w<0.74`$ at the $`95\%`$ confidence level.

It is straightforward to examine the cosmological consequences of a dark energy component which is strictly constant throughout space, and evolving with any value of $`w`$. Any physical example of such a component, however, will necessarily have fluctuations (so long as $`w1`$). It is therefore important to determine whether these fluctuations can lead to a catastrophic destablization of the vacuum. Unfortunately, because $`p=w\rho `$ is a phenomenological description valid for a certain configuration rather than a true equation of state, specifying $`w`$ is not enough to sensibly discuss the evolution of perturbations, since $`\delta pw\delta \rho `$. We must therefore choose a specific model. In particular, the simplest way to obtain a phantom component ($`w<1`$) is to consider a scalar field $`\varphi `$ with negative kinetic and gradient energy Caldwell:1999ew ,

$$\rho _\varphi =\frac{1}{2}\dot{\varphi }^2\frac{1}{2}(\varphi )^2+V(\varphi ).$$
(8)
Fluctuations in this field have a negative energy, and it may be possible for the vacuum to decay into a collection of positive-energy and negative-energy particles. If the timescale for such an instability is less than the age of the universe, the phantom component would not be a viable candidate for dark energy.

Our goal in this paper is to ask whether phantom components are necessarily plagued by vacuum instability, and hence whether observers should take seriously the possibility that $`w<1`$. We will start with a survey of energy conditions and model-independent considerations in Section II. In Section III, we describe the cosmology of $`w<1`$ models more thoroughly and investigate linear perturbations in a cosmological model, demonstrating that it can be compatible with current observations. To investigate stability beyond linear order, in Section IV we consider a classical toy model of a phantom harmonic oscillator coupled to an ordinary oscillator, and demonstrate numerically that, for small perturbations and sufficiently small values of the coupling constant, there exist both stable and unstable regions of parameter space.

The heart of our paper is in Section V, where we consider the field theory of a phantom component coupled to gravity, and calculate the decay rate of a single phantom particle into several phantoms plus gravitons. The rate is naively infinite, due to the infinite phase space of high-momentum particles. We argue that, considering the phantom Lagrangian as an effective theory, the rate may be rendered finite by imposing a momentum cutoff and demonstrate that, if we restrict attention to couplings in the potential, then, for a momentum cutoff not far below the Planck scale, and for a suitable potential, phantom quanta may be stable against decay into gravitons and other particles over a timescale long compared to the age of the universe. However, when we include derivative couplings of phantom particles to gravitons, we find that such operators can lead to unacceptably short lifetimes for phantom particles.

## II Classical energy conditions

In classical general relativity, without having a specific model for the matter sources, we can nevertheless invoke energy conditions which restrict the form of the energy-momentum tensor $`T_{\mu \nu }`$. In this section we will briefly review these energy conditions, discuss which are relevant for cosmology, and compare them to the condition $`w1`$. A related discussion can be found in McInnes:2002qw .

Each of the energy conditions can be stated in a coordinate-invariant way, in terms of $`T_{\mu \nu }`$ and some vector fields of fixed character (timelike/null/spacelike). For purposes of physical insight, it is often helpful to consider the case of a perfect fluid, for which the energy-momentum tensor takes the form

$$T_{\mu \nu }=(\rho +p)U_\mu U_\nu +pg_{\mu \nu },$$
(9)
where $`\rho `$ is the energy density, $`p`$ the pressure, $`U^\mu `$ the fluid four-velocity, and $`g_{\mu \nu }`$ the metric. (More precisely, $`\rho `$ and $`p`$ are the energy density and pressure as measured in the rest frame of the fluid, but the shorthand designations are standard.) Our metric signature convention throughout this paper is ($``$$`+`$$`+`$$`+`$).

The most common energy conditions are the following:

* The Weak Energy Condition or WEC states that $`T_{\mu \nu }t^\mu t^\nu 0`$ for all timelike vectors $`t^\mu `$, or equivalently that $`\rho 0`$ and $`\rho +p0`$.
* The Null Energy Condition or NEC states that $`T_{\mu \nu }\mathrm{}^\mu \mathrm{}^\nu 0`$ for all null vectors $`\mathrm{}^\mu `$, or equivalently that $`\rho +p0`$.
* The Dominant Energy Condition or DEC includes the WEC ($`T_{\mu \nu }t^\mu t^\nu 0`$ for all timelike vectors $`t^\mu `$), as well as the additional requirement that $`T^{\mu \nu }t_\mu `$ is a non-spacelike vector (i.e., that $`T_{\mu \nu }T^\nu {}_{\lambda }{}^{}t_{}^{\mu }t^\lambda 0`$). For a perfect fluid, these conditions together are equivalent to the simple requirement that $`\rho |p|`$.
* The Null Dominant Energy Condition or NDEC is the DEC for null vectors only: for any null vector $`\mathrm{}^\mu `$, $`T_{\mu \nu }\mathrm{}^\mu \mathrm{}^\nu 0`$ and $`T^{\mu \nu }\mathrm{}_\mu `$ is a non-spacelike vector. The allowed density and pressure are the same as for the DEC, except that negative densities are allowed so long as $`p=\rho `$.
* The Strong Energy Condition or SEC states that $`T_{\mu \nu }t^\mu t^\nu \frac{1}{2}T^\lambda {}_{\lambda }{}^{}t_{}^{\sigma }t_\sigma `$ for all timelike vectors $`t^\mu `$, or equivalently that $`\rho +p0`$ and $`\rho +3p0`$.

In Figure 1 we have plotted these conditions as restrictions on allowed regions of the $`\rho `$-$`p`$ plane. We have also plotted the condition $`w1`$ for comparison. Note that $`w1`$ is not equivalent to any of the energy conditions, although it is implied by the WEC, the DEC, and the NDEC.

The different energy conditions are used in different contexts — for example, the WEC and SEC are used in singularity theorems — and we will not review them in detail here (see hawkingellis  for a discussion). Our present concern is to understand under what conditions a hypothetical dark energy component would be guaranteed to be stable. For this purpose, the relevant result is the “conservation theorem” of Hawking and Ellis hawkingellis ; Carter:2002wz . The conservation theorem invokes the DEC, and uses it to show that energy cannot propagate outside the light cone; in particular, if $`T_{\mu \nu }`$ vanishes on some closed region of a spacelike hypersurface, it will vanish everywhere in the future Cauchy development of that region — energy-momentum cannot spontaneously appear from nothing. A source obeying the DEC is therefore guaranteed to be stable.

For cosmological purposes, however, the DEC is somewhat too restrictive, as it excludes a negative cosmological constant (which is physically perfectly reasonable, even if it is not indicated by the data). This is why Garnavich:1998th  advocated use of the NDEC in cosmology, since the NDEC is equivalent to the DEC except that negative values of $`\rho `$ are allowed so long as $`p=\rho `$. Of course, the less restrictive NDEC invalidates the conservation theorem, as a simple example shows. Consider a theory with a negative vacuum energy $`\rho _{\mathrm{vac}}<0`$, and an ordinary scalar field $`\psi `$ (not a phantom) with potential $`V(\psi )=\frac{1}{2}m^2\psi ^2`$. Then both the vacuum energy and the scalar field obey the NDEC (although their sum may not, since the NDEC is a nonlinear constraint.) Imagine a field configuration with $`\psi =\sqrt{2\rho _{\mathrm{vac}}/m^2}`$ and $`\dot{\psi }=0`$ everywhere throughout space. The energy-momentum tensor of this configuration is exactly zero \[since $`V(\psi )=\rho _{\mathrm{vac}}`$\], but it would instantly begin evolving as $`\psi `$ rolled toward the minimum of its potential. There is nothing unphysical about such a situation, however, since a stable state is achieved once the field reaches $`\psi =0`$. The NDEC, therefore, does not guarantee stability with the confidence that the DEC does. However, it seems legitimate to ask that dynamical fields obey the DEC, while the cosmological constant is allowed as an exception.

The DEC implies that $`w1`$. (Indeed, for cosmological purposes we are interested in a source with $`\rho >0`$; in that case, all of the energy conditions imply $`w1`$.) Therefore, if we allow for pressures which are less than $`\rho `$, we cannot guarantee the stability of the vacuum. The converse, however, is not true; a phantom component will not necessarily allow vacuum decay. In fact, since $`p=w\rho `$ is just a convenient parameterization and not strictly-speaking an equation of state, the issue of stability is somewhat complicated, as we shall see in the next few sections. Nevertheless, it seems very likely that energy sources which violate the DEC generally will be unstable, and this is what we will find in the specific example considered in Section V. Our philosophy, therefore, is neither to dismiss the possibility of DEC violation on the grounds that we cannot prove stability, nor to blithely accept DEC violation on the grounds that we cannot prove instability, but instead to see whether it is plausible that the timescale for instability could be sufficiently long so as to be irrelevant for practical purposes.

## III Cosmological Evolution and Perturbations

Consider a flat Robertson-Walker universe with metric

$$ds^2=dt^2+a^2(t)[dx^2+dy^2+dz^2],$$
(10)
for which the Einstein equations are the Friedmann equations (5) and (6). The cosmology and fate of a universe containing an energy component with constant $`w<1`$ are relatively simple and have been examined in, for example, Caldwell:1999ew . As an example, consider the case in which the universe contains only dust and phantom matter. Then, if the universe ceases to be matter-dominated at cosmological time $`t_m`$, then the solution for the scale factor is

$$a(t)=a(t_m)\left[w+(1+w)\left(\frac{t}{t_m}\right)\right]^{2/3(1+w)}.$$
(11)
From this expression it is easy to see that phantom matter eventually comes to dominate the universe and that, since the Ricci scalar is given by

$$R=\frac{4(13w)}{3(1+w)^2}\left[t\left(\frac{w}{1+w}\right)t_m\right]^2,$$
(12)
there is a future curvature singularity at $`t=wt_m/(1+w)`$. This occurs because, even though the energy densities in ordinary types of matter are redshifting away, the energy density in phantom matter increases in an expanding universe. Thus, the fate of the universe in these models Starobinsky:1999yw  may be very different from that expected Starkman:1999pg ; Krauss:1999hj ; Avelino:2000ix ; Gudmundsson:2001gd ; Huterer:2002wf  in $`w>1`$ dark energy models.

It is, however, simple to construct models of phantom energy in which a future singularity is avoided. As we shall see in this section, scalar field models can yield a period of time in which the expansion proceeds with $`w<1`$ and yet settles back to $`w1`$ (in our case $`w=1`$) at even later times, thus sidestepping the predictions of $`w=`$ constant models. Consider a scalar field theory with action $`S=\sqrt{|g|}d^4x`$, and Lagrange density given by

$$=\frac{1}{2}g^{\mu \nu }(_\mu \varphi )(_\nu \varphi )V(\varphi ).$$
(13)
The notable feature of this model is that the sign of the kinetic term is reversed from its conventional value \[in our conventions the usual expression would be $`\frac{1}{2}(\varphi )^2`$\]. The equation of motion for $`\varphi `$ becomes

$$\ddot{\varphi }+3H\dot{\varphi }a^2^2\varphi V^{}(\varphi )=0,$$
(14)
where $`^2`$ is the spatial Laplacian, $`^2\varphi =_x^2\varphi +_y^2\varphi +_x^2\varphi `$, and a prime denotes differentiation with respect to $`\varphi `$. Here the energy density $`\rho _\varphi `$ and the pressure $`p_\varphi `$ for a homogeneous $`\varphi `$ field are given by

$$\rho _\varphi =\frac{1}{2}\dot{\varphi }^2+V(\varphi ),$$
(15)
$$p_\varphi =\frac{1}{2}\dot{\varphi }^2V(\varphi ),$$
(16)
so that the equation-of-state parameter

$$w=\frac{p}{\rho }=\frac{\frac{1}{2}\dot{\varphi }^2+V(\varphi )}{\frac{1}{2}\dot{\varphi }^2V(\varphi )}$$
(17)
satisfies $`w1`$.

We are interested in the cosmological evolution of this model and in the behavior of linearized perturbations of the phantom scalar field in the resulting cosmological background. As noted in Caldwell:1999ew , the spectrum of fluctuations of a phantom field evolves similarly to that for a quintessence field. Consider metric perturbations in the synchronous gauge,

$$ds^2=dt^2+a^2(\delta _{ij}+h_{ij})dx^idx^j.$$
(18)
A Fourier mode of the phantom field

$$\varphi (t,𝐤)=\frac{1}{\sqrt{2\pi }}\varphi (t,𝐱)e^{i𝐤𝐱}d^3x$$
(19)
satisfies the equation of motion

$$\ddot{\delta \varphi _k}+3H\dot{\delta \varphi _k}+(k^2V^{\prime \prime })\delta \varphi _k=\frac{1}{2}\dot{h}\dot{\varphi },$$
(20)
where $`h`$ is the trace of the synchronous gauge metric perturbation $`h_{ij}`$. The effective mass for the perturbation is $`(k^2V^{\prime \prime })^{1/2}`$. On large scales one may worry that this effective mass could become imaginary for a positive $`V^{\prime \prime }`$. We note, however, that a similar problem exists for a canonical scalar field for negative $`V^{\prime \prime }`$, and is thus easily avoided by the choice of potential; in particular, we may choose a potential with $`V^{\prime \prime }<0`$. An analysis of perturbations in more general non-canonical scalar field models is given in Garriga:1999 . In this paper, we will examine by explicit calculation the evolution of fluctuations in a specific model.

Since scalar fields with negative kinetic terms evolve to the maxima of their classical potential, we consider a gaussian potential,

$$V(\varphi )=V_0e^{(\varphi ^2/\sigma ^2)},$$
(21)
where $`V_0`$ is the overall scale and $`\sigma `$ is a constant describing the width of the gaussian. The potential is represented in Figure 2. We obtained the cosmological evolution by numerically solving the equations of motion (5), (6), and (14). The initial conditions were $`\varphi _{\mathrm{initial}}=M_\mathrm{P}`$ and $`\dot{\varphi }_{\mathrm{initial}}=0`$, where $`M_\mathrm{P}=(8\pi G)^{1/2}`$ is the reduced Planck mass. Since we require dark energy domination at the present epoch in the universe, with $`\mathrm{\Omega }_X0.7`$ we choose $`\sigma =M_\mathrm{P}`$, $`V_0=3M_\mathrm{P}^2H_0^2`$. The results are plotted in figures 3,4 and 5.

Note that during the initial stages of evolution the field $`\varphi `$ is frozen by the expansion and acts as a negligibly small vacuum energy component (with $`w1`$). However, at later times the field begins to evolve more rapidly towards the maximum of its potential, the energy density in the phantom field becomes cosmologically dominant, and during this period the equation of state parameter is much more negative. Finally, in the very late universe, the field comes to rest at the maximum of the potential and a period of $`w=1`$ acceleration begins. Since $`w`$ is no longer less than $`1`$, this ensures that there is no future singularity; rather, the universe eventually settles into a de Sitter phase.

Using the formalism for calculating fluctuations in a general matter field developed in Hu:1998 , one may calculate the effects of the fluctuations in this phantom field on the cosmic microwave background radiation (CMB). The resulting power spectrum of CMB fluctuations is shown in figure 6.

We have not attempted a detailed parameter-fitting to CMB data, but it should be clear that our phantom model with the potential (21) does not predict any significant departures from conventional dark-energy scenarios; in particular, there is no evidence of dramatic instabilities distorting the power spectrum. (A more detailed study of the effect on the CMB power spectrum of phantom fields with a constant $`w`$ can be found in Schulz:2001yx .) However, the formalism for this analysis was based on linear perturbation theory, and it is quite plausible that instabilities only become manifest at higher orders. In the following sections we address this issue, first through numerical investigation of a model with two oscillators, and next by calculating the decay rate of phantom particles into phantoms plus gravitons.

## IV Couplings to Normal Matter and Instability

Excitations of the scalar field from the model in the previous section have negative energy. The existence of negative energy particles may cause the system to be unstable due to interactions involving these particles.

In the next section we examine this possibility by estimating the tree level decay rate of negative energy phantom particles into other phantoms and gravitons. Before delving into that calculation, we will first try to build some intuition about possible instabilities by considering the classical evolution of a simple system, a coupled pair of simple harmonic oscillators. Conservation of energy limits the phase space of a coupled pair of oscillators when both oscillators have positive energy. Though the oscillators may exchange energy, neither oscillator may ever reach an energy greater than the total initial energy of the system. Allowing one of the oscillators to have negative energy removes this limitation on the phase space. In this case the positive energy oscillator may increase its energy to any level so long as the negative energy oscillator decreases its energy by a compensating amount. Hence the positive energy oscillator may reach arbitrarily high energies, and the negative energy oscillator may reach correspondingly large negative energies.

While conservation of energy does not prevent the oscillators in this simple model from reaching arbitrarily large energies, there is no guarantee that the system will be unstable in all regimes. In fact, the following analysis will show that for a weak coupling, the evolution of the energy of each oscillator exhibits a stable oscillation for an arbitrarily long time.

We consider a coupled pair of oscillators, one $`\psi `$ with positive energy, representing normal matter, and the other $`\varphi `$ with negative energy, representing the phantom component.

$$=\left(\frac{1}{2}\dot{\psi }^2m_\psi ^2\psi ^2\right)\left(\frac{1}{2}\dot{\varphi }^2m_\varphi ^2\varphi ^2\right)\lambda \varphi ^2\psi ^2,$$
(22)
where $`\lambda `$ is a dimensionless coupling constant. We are interested in the evolution of the energy of each oscillator,

$`\rho _\psi `$ $`=`$ $`{\displaystyle \frac{1}{2}}\dot{\psi }^2+m_\psi ^2`$
$`\rho _\varphi `$ $`=`$ $`\left({\displaystyle \frac{1}{2}}\dot{\varphi }^2+m_\varphi ^2\right).`$ (23)
The equations of motion are

$`\overline{\psi }^{\prime \prime }`$ $`=`$ $`\left[\left({\displaystyle \frac{m_\psi }{m_\varphi }}\right)^2+\overline{\lambda }\overline{\varphi }^2\right]\overline{\psi }`$
$`\overline{\varphi }^{\prime \prime }`$ $`=`$ $`\left[1\overline{\lambda }\overline{\psi }^2\right]\overline{\varphi },`$ (24)
where we have defined dimensionless rescaled variables

$`\overline{\varphi }`$ $`=`$ $`\varphi /M`$ (25)
$`\overline{\psi }`$ $`=`$ $`\psi /M`$ (26)
$`\overline{\lambda }`$ $`=`$ $`\lambda (M/m_\varphi )^2`$ (27)
$`\tau `$ $`=`$ $`m_\varphi t.`$ (28)
In this section a prime denotes diferentiation with respect to the dimensionless time parameter $`\tau `$, and $`M`$ is the scale of the initial displacements of the oscillators.

We explore the stability of this model by integrating the equations of motion (IV) numerically for a range of parameters and initial conditions. Of particular interest are the simulations with $`m_\psi =0`$, in which one may think of the $`\psi `$ oscillator as analogous to the massless graviton. Three such simulations are shown in Fig. 7. In each of these integrations the oscillators are started at rest, displaced a distance $`M`$ from the origin. These plots are similar to those obtained for other initial conditions, and allowing the $`\psi `$ field to have a mass does not qualitatively change the plots shown here. Note that for $`\overline{\lambda }<1`$, $`\rho _\psi `$ and $`\rho _\varphi `$ exhibit a stable oscillatory behavior, while for $`\overline{\lambda }>1`$, $`\rho _\psi `$ and $`\rho _\varphi `$ rapidly grow. From this analysis we conclude that this simple model is stable for small enough coupling but exhibits an instability for a large coupling. While the analysis of this simple model does not prove that a phantom field in the cosmological context will be stable to excitations and interactions with the gravitational field, it does provide some evidence that stability is possible.

To make some connection to the cosmological model of the previous section, we consider how small the coupling in this simple model should be in order to have a stable solution. If we imagine a period of early universe inflation, then one expects perturbations in the phantom and gravitational fields to be of order $`M10^5M_\mathrm{P}`$. The mass squared of the phantom field is

$$m_\varphi ^2=\frac{d^2V}{d\varphi ^2}\frac{V_0}{M_\mathrm{P}^2}(10^{33}\mathrm{eV})^2.$$
(29)
The restriction $`\overline{\lambda }<1`$ then implies that the coupling in the Lagrangian satisfies $`\lambda <10^{110}`$. While this at first seems to be an absurdly small number, we will see in the next section that, given the cosmological constraints of the previous section, such a coupling naturally results from considering perturbations of the phantom and gravitational fields.

## V Phantom Decay Rate in Field Theory

The toy model of the previous section demonstrated that coupling a phantom oscillator to an ordinary oscillator results in a system which may or may not be unstable, depending on the magnitude of the coupling between them. It is far from clear, however, that this conclusion extends immediately to field theory. Roughly speaking, oscillators with some frequency correspond to field-theory modes of fixed wavelength; even if the model is stable when only certain wavelengths are considered, it does not follow that stability continues to obtain when integrating over all momenta.

One way of stating this concern is to ask about the decay rate of single particles that would conventionally be stable. Because excitations of the phantom field have negative energy, we could imagine a single particle decaying into a large number of phantoms and ordinary particles. The rate for this process can be calculated (at tree level) using ordinary Feynman diagrams. We will find that the rate is infinite when arbitrarily high momenta are included, but can be rendered finite if a cutoff (which might arise, for example, from higher-derivative terms in the action) is introduced.

Since we are interested in the allowed decay modes of phantom particles $`\varphi `$ and their associated decay rates, it is instructive to first analyze the kinematics of reactions involving phantom fields. Let us adopt the convention that ordinary particles have 4-momenta with positive timelike component. Then, since our field has a negative kinetic term, the corresponding 4-momentum will have a negative timelike component. This leads to the following useful dictionary for translating between kinematically allowed reactions for phantom particles and those for ordinary particles.

Denote ordinary particles as $`\psi `$. To ask whether a certain reaction is kinematically allowed, we can just switch the phantom particles from the right side to the left and vice-versa, and ask whether the resulting reaction would normally be allowed. Consider for example the decay of an ordinary particle into another particle plus a phantom:

$$\psi _1\psi _2+\varphi .$$
(30)
This will be allowed if the reaction

$$\psi _1+\varphi \psi _2$$
(31)
would be allowed by conventional kinematics — in particular, if the mass of ordinary particle 2 were greater than the sum of the masses of ordinary particle 1 and the phantom. So it is clear that ordinary particles can decay into heavier particles plus phantoms. For example, if the electron were coupled to the phantom field, processes such as

$$e^{}\mu ^{}+\nu _e+\overline{\nu }_\mu +\varphi $$
(32)
would be allowed. The muon would then decay back into an electron, as part of a potentially disasterous cascade.

Next we turn to the decay of phantoms. First consider decays into ordinary particles

$$\varphi \psi _1+\psi _2$$
(33)
will be allowed if

$$|0\varphi +\psi _1+\psi _2$$
(34)
would be conventionally allowed, which it is not. However, consider a decay into one phantom and one ordinary particle:

$$\varphi _1\varphi _2+\psi $$
(35)
will be allowed if

$$\varphi _2\varphi _1+\psi $$
(36)
would ordinarily be, which requires $`m_{\varphi _2}>m_{\varphi _1}+m_\psi `$. So a phantom can decay into a heavier phantom plus a not-too-heavy ordinary particle.

If there is only one kind of phantom, one might think it would be stable since there would be no heavier phantoms to decay into. However, several lighter particles can mimic the 4-momentum of a heavier particle. Consider the decay of one phantom into two phantoms plus an ordinary particle:

$$\varphi _1\varphi _2+\varphi _3+\psi .$$
(37)
This will be allowed if

$$\varphi _2+\varphi _3\varphi _1+\psi $$
(38)
would ordinarily be, which it aways is for large enough relative velocities of $`\varphi _2`$ and $`\varphi _3`$.

In summary, if there is only one kind of phantom particle (with a unique mass), it can only decay by emitting at least two more phantoms, plus at least one ordinary particle. Ordinary particles, meanwhile, may decay into phantoms plus other ordinary particles with a larger effective mass than the original. A special case to these rules comes from massless particles. Although it would seem kinematically possible for massless particles to decay into massive particles plus phantoms, massless particles cannot decay in flat space (a reasonable approximation in the backgrounds we are considering) simply because there is no rest frame in which to calculate the rate — no proper time elapses along a null path.

Like any other dark-energy scalar field, the phantom should be weakly coupled to ordinary matter (or it would have been detected through fifth-force experiments or variation of the constants of nature Carroll:1998zi ). We therefore restrict our attention to only gravitons and phantoms. With the above rules in mind, the phantom decay channel involving the smallest number of particles is

$$\varphi _ih+\varphi _1+\varphi _2,$$
(39)
where $`h`$ is a graviton and $`\varphi `$ is a phantom particle, illustrated in figure 8.

Consider the specific model investigated in section III, a phantom scalar with potential

$$V(\varphi )=V_0e^{\varphi ^2/M_\mathrm{P}^2}.$$
(40)
We will first consider this potential expanded as a power series around some background value $`\varphi _0M_\mathrm{P}`$. Gravitons, meanwhile, may be represented by transverse-traceless metric perturbations $`g_{\mu \nu }=g_{\mu \nu }^{(0)}+h_{\mu \nu }`$, where $`g_{\mu \nu }^{(0)}`$ is the background Robertson-Walker metric. Strictly speaking, there is no way for a single graviton to couple non-derivatively to the potential, simply because there is no way to construct a scalar from a single traceless $`h_{\mu \nu }`$. But to avoid a surfeit of indices, and because it won’t affect the final answer, we will simply think of the graviton as a dimensionless scalar field $`h`$; to get a canonically-normalized field we multiply by $`M_\mathrm{P}`$.

To study the decay (39), we require the interaction part of the Lagrangian to first order in $`h`$ and third order in $`\varphi `$, which is

$`_I`$ $`=`$ $`{\displaystyle \frac{1}{M_\mathrm{P}}}(M_\mathrm{P}h){\displaystyle \frac{1}{3!}}V^{\prime \prime \prime }(\varphi _0)\varphi ^3`$ (41)
$`=`$ $`\lambda _{\mathrm{eff}}(M_\mathrm{P}h)\varphi ^3,`$ (43)
where

$$\lambda _{\mathrm{eff}}2\varphi _0\frac{V(\varphi _0)}{M_\mathrm{P}^5}\frac{V_0}{M_\mathrm{P}^4}10^{120}.$$
(44)
The decay rate of a phantom particle through this channel is

$$\mathrm{\Gamma }=\frac{1}{m_\varphi }\frac{d^3p_h}{(2\pi )^32E_h}\frac{d^3p_{\varphi _1}}{(2\pi )^32E_{\varphi _1}}\frac{d^3p_{\varphi _2}}{(2\pi )^32E_{\varphi _2}}||^2(2\pi )^4\delta ^{(4)}(p_{\varphi _i}p_{\varphi _1}p_{\varphi _2}p_h),$$
(45)
where the matrix element $`||`$ is just $`\lambda _{\mathrm{eff}}`$ at tree level. To get an upper limit on the reaction rate, and hence a lower limit on the timescale, we assume approximate isotropy, so $`d^3p|p|^2dp`$. Because the relevant momenta are very large (and the masses very small), we can also approximate $`Ep`$. Putting all of this together we obtain

$$\mathrm{\Gamma }\frac{\lambda _{\mathrm{eff}}^2}{m_\varphi }|p_h|𝑑p_h|p_{\varphi _1}|𝑑p_{\varphi _1}|p_{\varphi _2}|𝑑p_{\varphi _2}\delta ^{(4)}(p_{\varphi _i}p_{\varphi _1}p_{\varphi _2}p_h).$$
(46)
If we take the limits on the integrals to be $`\mathrm{}`$, the decay rate is clearly infinite. Hence, the answer to our investigation into instability seems very clear: the theory is dramatically unstable, as individual particles rapidly decay into cascades of phantoms and gravitons.

However, our philosophy has been to think of this model as an effective theory valid at low energies. The reason why the decay rate diverges is because the phase space is infinite, since the phantoms can have arbitrarily large negative energies. Therefore, this result relies on taking the calculation at face value up to infinite momentum transfer. Instead, we should only trust the phase-space integrals up to some cutoff where new physics might enter. For example, we could imagine a higher-derivative term of the form

$$(\varphi )^4,$$
(47)
which would eventually dominate over the negative-kinetic-energy term we have already introduced.

We have not investigated closely the properties of phantom models with higher-derivative terms. Instead, let us crudely approximate the effect of a momentum cutoff at a scale $`\mathrm{\Lambda }`$ by only allowing the phase-space integrals to range up to that cutoff. Looking at (46), we can estimate the truncated decay rate as

$$\mathrm{\Gamma }\lambda _{\mathrm{eff}}^2\frac{\mathrm{\Lambda }^2}{m_\varphi }.$$
(48)
The timescale for decay is $`\tau =\mathrm{\Gamma }^1`$. Despite any fundamental instability, a model will be phenomenologically viable if the lifetime is greater than the Hubble time $`H_0^110^{60}M_\mathrm{P}^1`$. Recalling that $`m_\varphi 10^{60}M_\mathrm{P}`$ and $`\lambda _{\mathrm{eff}}10^{120}`$, the lifetime in units of the Hubble time is

$$H_0\tau 10^{120}\left(\frac{M_\mathrm{P}}{\mathrm{\Lambda }}\right)^2.$$
(49)
In other words, the lifetime from this decay channel exceeds cosmological timescales so long as $`\mathrm{\Lambda }<10^{60}M_\mathrm{P}`$, which is certainly not a stringent constraint.

However, the infinite phase space for this one decay is not the only infinity we have to deal with; a single phantom can decay into arbitrary numbers of gravitons and phantoms, so we must sum over all the channels. For each new final-state particle $`\mathrm{\Gamma }`$ gets multiplied by a factor of $`(\mathrm{\Lambda }/M_\mathrm{P})^2`$. The $`\mathrm{\Lambda }^2`$ comes from the new momentum integral, and the $`M_\mathrm{P}^2`$ is introduced into $`\lambda _{\mathrm{eff}}^2`$. So, if we denote the decay rate (48) by $`\mathrm{\Gamma }_0`$, then if there are $`n`$ additional particles in the final state the total decay rate is given by

$$\mathrm{\Gamma }_n=\left(\frac{\mathrm{\Lambda }}{M_\mathrm{P}}\right)^{2n}\mathrm{\Gamma }_0,$$
(50)
Therefore, the total decay rate is

$$\mathrm{\Gamma }_{\mathrm{total}}=\underset{n=0}{\overset{\mathrm{}}{}}(n+1)\mathrm{\Gamma }_n,$$
(51)
where the factor $`(n+1)`$ comes from the number of different ways the final state can be composed of gravitons and phantoms. This yields

$$\mathrm{\Gamma }_{\mathrm{total}}=\mathrm{\Gamma }_0\left[1\left(\frac{\mathrm{\Lambda }}{M_\mathrm{P}}\right)^2\right]^2.$$
(52)
Thus, the decay rate remains of order $`\mathrm{\Gamma }_0`$ so long as $`\mathrm{\Lambda }`$ is not larger than $`M_\mathrm{P}`$. This still seems like a comfortable result.

There is, however, one more possibility to be accounted for. Part of the reason we obtained reasonable decay rates from interactions originating in the potential $`V(\varphi )`$ was because of the small value of $`V_0(10^3\mathrm{eV})^4`$. This suppression is lost if we consider couplings of gravitons to derivatives of $`\varphi `$. Since we are claiming that our model is simply an effective field theory, we are obligated to include all possible non-renormalizable interactions, suppressed by appropriate powers of the cutoff scale. Consider for example the operator

$$=\frac{\beta }{M_\mathrm{P}\mathrm{\Lambda }}\varphi (M_\mathrm{P}h^{\mu \nu })_\mu \varphi _\nu \varphi ,$$
(53)
where $`\beta `$ is a dimensionless coupling. (We will use the same cutoff $`\mathrm{\Lambda }`$ for our non-renormalizable terms as in our momentum integrals; the $`M_\mathrm{P}`$ in the denominator comes from normalization of $`h^{\mu \nu }`$.) This term will also contribute to the interaction shown in Figure 8. Following similar logic as above (including two powers of the momentum in $``$ from the derivatives), we get

$$\mathrm{\Gamma }\frac{\beta ^2\mathrm{\Lambda }^4}{m_\varphi M_\mathrm{P}^2}.$$
(54)
Using $`H_0m_\varphi 10^{60}M_\mathrm{P}`$, the lifetime in units of the Hubble time is now

$$H_0\tau \beta ^2\left(\frac{M_\mathrm{P}}{10^{30}\mathrm{\Lambda }}\right)^4.$$
(55)
Therefore, if $`\beta `$ is of order unity, to obtain cosmologically viable decay rates ($`H_0\tau >1`$) we require the cutoff to be

$$\mathrm{\Lambda }<10^{30}M_\mathrm{P}10^3\mathrm{eV}.$$
(56)
This is a much smaller cutoff than was required when we considered coupling through the potential, since the small prefactor $`V_0`$ is not around to help us.

The above result is alleviated somewhat by imposing an approximate global symmetry on the theory, that the Lagrangian density be invariant under $`\varphi \varphi +\mathrm{constant}`$. Such a symmetry is quite reasonable, as it is the only known way to ensure both a nearly-flat potential and appropriately small couplings to ordinary matter Carroll:1998zi . In this case, the irrelevant operator of lowest dimension leading to phantom decay into gravitons is

$$=\frac{\gamma }{M_\mathrm{P}^2\mathrm{\Lambda }^4}\left[(M_\mathrm{P}h^{\mu \nu })_\mu \varphi _\nu \varphi \right]^2.$$
(57)
The decay rate is then

$$\mathrm{\Gamma }\frac{\gamma ^2\mathrm{\Lambda }^6}{m_\varphi M_\mathrm{P}^4},$$
(58)
with an associated lifetime

$$H_0\tau \gamma ^2\left(\frac{M_\mathrm{P}}{10^{20}\mathrm{\Lambda }}\right)^6,$$
(59)
which, if $`\gamma 1`$, leads to the relaxed requirement $`\mathrm{\Lambda }<10^{20}M_\mathrm{P}100`$ MeV.

Under our most optimistic assumptions (of an approximate global symmetry), we therefore find that the momentum cutoff characterizing our effective theory must be less than 100 MeV to guarantee that the instability timescale is greater than the age of the universe. We find this value to be uncomfortably low, although not necessarily impossible; dark-energy models are typically characterized by energy scales $`10^3`$ eV (for the vacuum energy today) and $`10^{33}`$ eV (for the mass of the scalar field), so perhaps the scalar field theory is only valid up to relatively low momenta. Alternatively, one might imagine searching for some mechanism which would suppress derivative couplings, leaving only the couplings of gravitons to the potential, which were consistent with a cutoff as high as the Planck scale.

## VI Conclusions

There is no doubt that the discovery of a new component of the energy density of the universe has profound implications for the relationship between particle physics and gravity. Whether this component be a pure cosmological constant (whose magnitude we have no idea how to understand), a dynamical component Sahni:1999gb ; Parker:1999td ; Chiba:1999ka ; Boisseau:2000pr ; Schulz:2001yx ; Faraoni:2001tq ; maor2 ; Onemli:2002hr ; Torres:2002pe ; Frampton:2002tu ; Wetterich:fm ; Ratra:1987rm ; Caldwell:1997ii ; Armendariz-Picon:1999rj ; Armendariz-Picon:2000dh ; Armendariz-Picon:2000ah ; Mersini:2001su  (whose special interactions give rise to the tiny vacuum energy we observe) or an as yet unimagined source, its nature is an outstanding problem of fundamental physics. The data thus far, although pointing conclusively to the existence of this dark energy, do not allow us to distinguish between competing scenarios. In particular, much of the allowed parameter space lies in the region $`w<1`$ in which some cherished notion behind our present theories must be sacrificed.

There are several ways in which one may achieve $`w<1`$, including purely negative kinetic terms, non-minimal kinetic terms, and scalar-tensor theories. In addition, one might imagine that effective superexponential expansion of the universe might be obtained by modifying the Friedmann equation. While this is at odds with a purely four-dimensional general relativistic description of the universe, such an effect might be obtained in the context of brane-world models Kehagias:1999aa . However, any such modification of the Friedmann equation must avoid conflict with the precision predictions of primordial nucleosynthesis Tegmark:2001zc ; Carroll:2001bv ; Zahn:2002rr .

In this paper we have taken seriously the possibility that $`w`$ may be less than $`1`$ and have asked the question: “From a particle-physics and general-relativistic point of view, can such theories be made consistent?” In particular, we have considered a specific toy model in which the null dominant energy condition is violated and hence the resulting space-time may be unstable. In this model the cosmology is well-behaved and the theory may be constructed so that it is stable to small, linear perturbations. When we consider higher order effects, however, the model fails to remain stable. The central result of this paper is a field theory calculation of the decay rate of phantom particles into gravitons. This decay rate would be infinite if the phantom theory was fundamental, valid up to arbitrarily high momenta, and would render the theory useless as a dark energy candidate. We therefore consider the phantom theory to be an effective theory valid below a scale $`\mathrm{\Lambda }`$. Interestingly, couplings of gravitons to an appropriate scalar potential do not lead to decay of phantoms into gravitons and other phantoms on sub-Hubble timescales, so long as the cutoff $`\mathrm{\Lambda }`$ is below the Planck scale. However, in such an effective field theory approach, we are mandated to include in the Lagrangian operators of all possible dimensions, suppressed by suitable powers of the cutoff scale. In particular, we must include couplings of gravitons to derivatives of the phantom field. We find that such operators, even though they may be of high order, can lead to unacceptably short lifetimes for phantom particles unless the cutoff scale is less than $`100`$ MeV, so new physics must appear in the phantom sector at scales lower than this.

Our analysis demonstrates that a model with a well-behaved cosmological evolution and stability to linear perturbations may still exhibit instability due to higher order interactions. We may therefore ask another crucial question: “Should observers seeking to constrain cosmological parameters take seriously the possibility that $`w<1`$?” Unfortunately the answer is somewhat ambiguous. On the one hand, we know of no easy way to construct a viable model of this sort; on the other, it is certainly conceivable that new physics in the dark-energy sector kicks in at low scales to render a phantom model stable, or that phantom behavior is mimicked by even more exotic mechanisms (such as modifications of the Friedmann equation). Therefore, whether or not observations constrain $`w`$ to be greater than or equal to $`1`$ is still an interesting question, although there is a substantial a priori bias against the possibility. For theorists, our conclusion is more straightforward: the onus is squarely on would-be phantom model-builders to show how any specific proposal manages to avoid rapid vacuum decay.

## Acknowledgments

We would like to thank Roger Blandford, Mark Bowick, Dan Chung, David Gross, Don Marolf, Laura Mersini and Emil Martinec for useful conversations. We especially thank Wayne Hu for his assistance with calculating power spectra. MT would like to thank the Kavli Institute for Theoretical Physics for kind hospitality and support during the final stage of this project. The work of SC and MH is supported in part by U.S. Dept. of Energy contract DE-FG02-90ER-40560, National Science Foundation grant PHY-0114422 (CfCP), the Alfred P. Sloan Foundation, and the David and Lucile Packard Foundation. The work of MT is supported in part by the NSF under grant PHY-0094122.


format: text

text:
# Energy landscape of relaxed amorphous silicon

## I Introduction

The dynamics of many complex materials is dominated by activated jumps over energy barriers generally higher than $`k_BT`$. For these systems, the energy-landscape picture, which focuses on the topological relation between locally stable states, has proven very valuable. Disconnectivity graphs, for example, introduced by Czerminski and Elber Czerminksi and Elber (1990) and applied extensively by others, Becker and Karplus (1997); Wales et al. (1998); Doye et al. (1999a, b); Mortenson and Wales (2001) have provided a first classification of the dynamics of complex systems based on the structure of their respective energy landscape.

In parallel with these developments, oriented towards reconstructing the topology of the energy landscape, there has been some efforts in trying to characterize the energetics and the nature of the events in clusters, Wales et al. (1998); Doye et al. (1999a, b); Malek and Mousseau (2000) proteins Mousseau et al. (2001); Mortenson and Wales (2001); G. H. Wei and Derreumaux (2002) and in amorphous materials. Barkema and Mousseau (1998); Mousseau and Barkema (2000); Song et al. (2000); Middleton and Wales (2001) This extensive sampling is essential in order to try to connect the properties of the landscape with the dynamics measured experimentally, it also serves to build a better understanding of the generic properties of various networks: low vs. high connectivity, bulk vs. finite-size systems, covalent vs. ionic bonding, etc.

In this paper, we present an extensive study of the structure of the energy landscape of amorphous silicon around two minima. Less extensive studies of this material were already presented both by our group Barkema and Mousseau (1998); Mousseau and Barkema (2000); Song et al. (2000) and Middleton et al.Middleton and Wales (2001) using two different approaches. Using ART nouveau, we generate more than 42 000 activated events around a well-relaxed minimum and analyze the properties of the reaction paths connecting this initial minimum to a nearby saddle point and new minimum. We find that : (1) the activation mechanisms in a-Si are local, limiting somewhat the usefulness of the configurational energy landscape picture, (2) the activation energy is essentially limited by the energy required to break a single bond, (3) the entropic barrier is almost identical for all events, (4) more than 20 % of all events generated are bond-switches of the Wooten-Winer-Weaire type, and (5) the number of events seems to be of the order 30 to 60 per atom.

## II Details of the simulation

We study here the energy landscape around two well-relaxed 1000-atom configurations of a-Si. We focus mostly on the first one, repeating the simulation on a second configuration only to ensure that the results are generic and not dependent on some specific feature of the configuration.

The energy model used is a modified version (mSW) of the empirical Stillinger-Weber (SW) potential Stillinger and Weber (1985); Vink et al. (2001) developed to ensure that the elastic and structural properties of the amorphous phase of silicon are in agreement with experiment. Vink et al. (2001) With its original parameters, the SW potential describes both the liquid and the crystalline phases with good accuracy but fails to reproduce the experimental structure of the amorphous phase. Ding and Andersen (1986); Cook and Clancy (1993); Barkema and Mousseau (1996) Recently, Vink and collaborators proposed a slightly different set of parameters (see Table 1), which generates the right amorphous structure in addition to providing the correct vibrational properties. Because of its empirical nature and the fact that it is not optimized for dynamics, the energy barriers computed with this potential must be taken with some care. However, since the structural properties of this material are well described by mSW, the qualitative features of the energy landscape are expected to be correct.

### II.1 ART nouveau

While our previous study of relaxation in a–Si used a version of ART that could not identify saddle points, Barkema and Mousseau (1996); Mousseau and Barkema (1998); Song et al. (2000) the results presented here are obtained using ART nouveau, the latest version of the activation-relaxation technique presented in Ref. Malek and Mousseau, 2000; G. H. Wei and Derreumaux, 2002. Using ideas similar to those proposed by Munro and Wales, Munro and Wales (1999) ART nouveau applies a Lanczos scheme to compute directly the lowest second derivative of the energy during the activation phase and ensures convergence to the saddle point to any desired precision.

Events are generated in the following way: In order to leave the harmonic well, one atom is selected at random. This atom and its neighbors, contained in a shell of radius 3.5 Å, are moved iteratively in a randomly chosen direction, while keeping the energy, projected in the perpendicular directions to a minimum. At each step, a Lanczos scheme is used to compute the lowest eigenvalues and eigenvectors associated with the curvature of the energy landscape. We consider that the configuration has left the harmonic well when the lowest eigenvalue falls below $`50`$ eV/Å2.

The activation process per se starts from this point. The configuration is slowly pushed up along the direction of lowest curvature until the modulus of the force falls below 0.5 eV/Å— indicating that the configuration has reached a saddle point— or until the lowest eigenvalue become positive — indicating that the trajectory selected is back in the initial harmonic well region and must be rejected. After reaching the saddle point, the configuration is pushed slightly away from it, and is relaxed into a new local energy minimum, called the “‘final minimum”. This event is stored and a new event is started from the same initial minimum.

### II.2 Properties of the initial configurations

The two initial configurations used here were prepared using ART nouveau and mSW. Starting from a 1000-atom randomly packed unit cell with periodic boundary conditions, ART events were applied until the configurational energy equilibrated. We used a Metropolis accept/reject algorithm, based the energy difference between consecutive local minima, as described in Ref. Barkema and Mousseau, 1996; Mousseau et al., 2000. The first configuration has an energy per atom of $`4.000`$ eV, with 20 five-fold and 26 three-fold defects. The radial distribution function is in good agreement with recent experimental data(Fig. 1). The model is therefore of comparable quality to the models discussed in Refs. Mousseau and Barkema, 1998; Barkema and Mousseau, 2000. This configuration is used as the origin of all events for the first 42 000-event run.

The second configuration was prepared by further applying ART on the initial configuration, at a Metropolis temperature $`T=0.25`$ eV, for a few thousands of events, until the average displacement per atom reached 1 Å. We generated more than 7000 events around this second minimum in order to ensure that the features of the energy landscape were statistically independent of the minimum selected. All the quantities analyzed here are the same for both minimum, confirming that the properties of the landscape are not affected by the details of the topology. In view of the similarities between the properties of these two sets of data, most of the discussion will focus on Minimum 1.

### II.3 Search for activated events

All events presented here are characterized by three configurations: the initial minimum (either Minimum 1 or 2), a first-order saddle point and a second (final) minimum, obtained by relaxing the configurational energy from the saddle point, in a direction opposite to that of the initial minimum. As discussed in Ref. Malek and Mousseau, 2000, these events are reversible both from the saddle point and the final minimum.

Each event is generated starting from the same initial minimum, but with a different random initial direction, the rest of the procedure being deterministic. It takes about 400 force evaluations to generated a single event. About one third of the random launches lead to a trajectory that brings the configuration back in the harmonic basin; these trajectories are simply rejected and a new random direction selected. The whole search took about 4 weeks on a single processor of a Regatta Power4 IBM.

## III Results

### III.1 Distributions on events

We first discuss the distributions of events generated from Minimum 1. 42 581 events were generated following the procedure discussed in the previous section. The distribution of activated ($`E_{\mathrm{saddle}}E_{\mathrm{initial}}`$) and asymmetry ($`E_{\mathrm{final}}E_{\mathrm{initial}}`$) energies for all these events are given in Fig. 2. The activated energies form a continuous spectrum, from 0 to more than 7 eV, with an average barrier height of 3.0 eV and a distribution width of 1.2 eV. Both the width and average barrier energies are much lower than those reported earlier using the previous version of ART. Mousseau and Barkema (2000) This underlines the limit of the original ART, which cannot converge directly to a saddle point. As shown in Ref. Song et al., 2000, the error on the barrier using the initial version of ART could be as high as 1 eV. As with the method of Munro and Wales, Munro and Wales (1999) on the contrary, ART nouveau makes it possible to identify the transition point to any desired accuracy. Comparing with experiment, we find that the average barrier height is in overall agreement with measurements of Shin and Atwater Shin and Atwater (1993) which indicate activation barriers extending from 0.25 eV to about 2.8. The typical activation energies we find correspond also to the isothermal calorimetric data of Roorda et al. which indicate high-activation barriers.

Although there are no experimental information to compare with, it is useful to examine the distribution of the asymmetry energy, i.e., the energy difference between the final and initial minima (see the bottom panel of Fig. 2). The average of the distribution is at 1.7 eV, with a width of 1.4 eV. There are only a few events with a final energy lower than the initial as should be expected, since the initial configuration is already very well-relaxed. As shown below, the distribution around Minimum 1 is essentially identical to that around Minimum 2. This results suggests that the overall shape of the configurational energy landscape is not sensitive to the details of the configuration, contrary to what one could think.

It is also useful to compare our results with those of Middleton and Wales obtained on the same system using the eigenvector following BFGS approach (EF-BFGS); Middleton and Wales (2001) ART nouveau is similar in spirit to this method. In Fig. 3, we plot the activated energy calculated from the initial minimum and the final ($`E_{\mathrm{saddle}}E_{\mathrm{final}}`$) for both sets of events. The barrier distribution calculated from the initial minimum show that the sampling of events differs seriously for both methods; while EF-BFGS seems to favor strongly events with a barrier below 2 eV, ART selects events in a more Gaussian way. Since both methods follow closely the direction corresponding to the lowest eigenvalue to the saddle point, this difference is solely due to the algorithm used to leave the harmonic well. Remarkably, however, this selection has very little impact on the barrier distribution calculated from the final minimum; in this case, both methods find the same barrier distribution, which is heavily skewed towards very lower barriers. We discuss the significance of these dissimilarities in Section IV.

While the energy distribution inform us about the barriers heights, the distribution of Hamming distances between configurations (Fig. 4) provide some insight as to the rearrangements taking place. As was found previously in a-Si and v-SiO2 most saddle points are somewhere at mid-way between the initial and final stage. Mousseau et al. (2000); Mousseau and Barkema (2000) This is somewhat unexpected in the energy landscape picture: in a high-dimensional space, any two randomly selected direction are orthogonal. For truly high-dimensional events, therefore, one would expect to find little correlation between the displacement at the saddle point and that at the final minimum. This is not the case, however, if events are taking place in a much restricted sub-space. The distribution of displacements indicates therefore that even though the simulation system is embedded in a 3000-dimensional space, the effective sub-space in which each event takes place is much smaller and local in nature; one must therefore be cautious when interpreting the dynamics of a material based solely on the energy landscape picture.

The local nature of the dynamics is also found plotting the number of active atoms participating in each event. (Fig. 5). An active atom is defined as one that has moved by more than 0.1 Å from its initial state. As discussed in Refs. Mousseau and Barkema, 1998; Mousseau et al., 2000, the value of the threshold is chosen to be close to the typical atomic displacement due to thermal vibration at room temperature. We see that most events involve a maximum of 25 atoms at the saddle point, and 40 at the final minimum. If the threshold is increased to 0.4 or 0.5 Å, the number of active events drops typically to between 2 and 4 per event. The relatively narrow distribution of active atoms is not an artifact due to some bias in the sampling of the cell. Fig. 6 shows the probability of each atom to participate into an event. The remarkably homogeneous figure confirms that the rearrangements can take place anywhere in the network, with a participation probability varying by at most a factor 3.

### III.2 Properties of the landscape

As the sampling of events is random, the total data set contains redundant events that might bias the analysis of the landscape. It is therefore useful to look at the properties of “unique” events, i.e., to analyze the set of single copies of all events. We identify each event based on the identity of the participating atoms as well as on the energy at the barrier and at the new minimum. In order to decrease the impact of the finite precision in the convergence at the saddle and the minimum, we define the participating atoms as those moving by at least $`\delta r_{\mathrm{threshold}}=0.4`$ and 0.2 Å, respectively, at the saddle point and the final minimum. Two events are considered identical if the active atoms are the same and the energy barrier differs by less than 0.2 eV. We verified that the precise values of the various thresholds do not affect the qualitative results although the precise number of different events obviously depends on the threshold. The number of unique events, as a function of the number of events already sampled, is shown in Fig. 7. We find 6519 different minima, 8799 different saddles and 11014 different events generated within the 42 000+ sequence. On average, each event is therefore visited almost 4 times.

Using the distribution of unique events, it is possible to assess the biases in the ART sampling. Fig. 8 shows the ratio of all saddle points or minima generated over the list of all unique ones. As was the case for the Lennard-Jones clusters, a system with a totally different energy landscape, ART sampling seems to select events with an exponential bias, $`\mathrm{exp}(\mathrm{\Delta }E/4)`$.

These results are obtained using a subset of all events surrounding Minimum 1. Even after 42000 events, the number of different events around minimum is not yet converged. It is possible to obtain some very crude estimate of this number. We first set the full distribution of saddle points identical to that shown in Fig. 7. Supposing a random selection of events, with an exponential bias on the barrier height such as that of the previous paragraph, taken from the distribution of unique events, we find that growth of the total of unique events as a function of trial can be reproduced with an exponential bias $`\mathrm{exp}(\mathrm{\Delta }E/A)`$ with $`A`$ between $`0.40`$ and $`0.60`$, and a total number of events somewhere between 30,000 and 60,000 or 30 to 60 different saddle points per atom.

### III.3 Entropy

Having computed the distribution of activation energy barriers around a minimum, it would be possible to consider associating a time scale with the events. This can be done within the framework provided by the transition state theory Glasstone et al. (1941). For activated mechanisms, the diffusion rate is given by

$$D=\frac{g}{2\alpha }\xi l^2\nu _0e^{\mathrm{\Delta }F/k_bT}$$
(1)
where $`g`$ is the number of equivalent diffusion paths, $`\alpha `$ the spatial dimension, $`l`$ the length of the jump, $`\nu _0`$ a phonon frequency, and $`\mathrm{\Delta }F=\mathrm{\Delta }ET\mathrm{\Delta }S`$, the variation in the Helmholtz free energy between the transition state and the initial state.

In order to use the above equation, we need to measure the change in entropy from the minimum to the transition state. We can do that using the harmonic approximation:

$$\mathrm{\Delta }S=k_b\mathrm{ln}\left[\frac{{\displaystyle \underset{i=1}{\overset{3N}{}}}\nu _i^{(i)}}{\nu _0{\displaystyle \underset{i=1}{\overset{3N1}{}}}\nu _i^{(s)}}\right]$$
(2)
where $`\nu _i^{(i)}`$ and $`\nu _i^{(s)}`$ are the real phonon frequencies at the minimum and at the saddle point respectively. Since there are $`3N1`$ real frequencies at the transition state, we replace the imaginary frequency by a typical phonon frequency, $`\nu _0`$.

Assessing the time scale associated with leaving a given minimum remains a very expensive task even within the harmonic approximation as it would be necessary to compute the entropy associated with each of the 30,000 to 60,000 saddle points. A simpler approach poses that the entropic prefactor is independent of the specific event.

We can ascertain the validity of such assumption by computing the distribution of the entropic difference from the initial minimum to the saddle point. Due to the high cost of diagonalizing a $`3000\times 3000`$ matrix, we have applied this procedure to 50 randomly-selected events. As seen in Fig. 9, the contribution of the entropy to the activation is of the order of 0.0024 eV/K, with a variation of about $`1.5\times 10^4`$ eV/K. The fluctuations in the entropic contributions to the barrier can therefore affect the attempt frequency by about a factor 5. To a first approximation, it is therefore reasonable to consider that the entropy is only a multiplicative constant in the dynamics of the system.

### III.4 Topological classification

In the previous sub-sections, we have analyzed the properties of the energy landscape itself. We now turn to the classification of the events generated around a single minimum. The topological classification used here is described in detail in Ref. Mousseau and Barkema, 1998. Briefly, it only considers changes in the coordination of the atoms, obtained through bond-breaking and binding. Analyzing the more than 40 000 events generated here, we find more than 2000 different types of topological events. Of those, we discuss below only the events that occurred more than 1 % of the time; about 400 times in our sequence. Only 5 types of event satisfy this criterion.

As in previous simulations, the most common event we find, representing about 20 % of all generated events, is the Wooten-Winer-Weaire bond exchange mechanism (or abacbd, in our notation), introduced almost 20 years ago in the sillium model. Since then, this mechanism has been seen both in crystalline and amorphous materials. Wooten et al. (1985); Stillinger and Weber (1985); Motooka (1996); Tang et al. (1997); Bernstein et al. (2000); Goedecker et al. (2002) Figure 10 shows the energy and Hamming distance distribution for this mechanism. Because of the frequency of this mechanism, it might not be surprising to find that both the energy and Hamming distance distributions follow closely those obtained for the whole set of events. These distributions underline once more, nevertheless, the rôle of strain in determining the activation barrier for a given mechanism: the same topological jump can lead to a new configuration with an energy varying by as much as 6 eV.

The topological nature of the other 4 dominant events is also shown in Fig. 11. The class $`abc`$ is associated with a bond jumping from on pair of atoms to the next. $`ab`$ represents a bond breaking at the saddle point and reforming. It does not involve any topological change and is responsible for the low-energy peak in the asymmetry energy distribution. As such it is not very interesting. Finally, the two other classes are modifications on the three most common types.

## IV Discussion and conclusions

The purpose of this study was to characterize in detail the energy landscape around a single minimum in a well-relaxed model of a-Si. Previous work has examined much shorter sequences of events as the configuration relaxed and the atoms diffused. Here, we performed two extensive simulations, always starting from the same two minima. Our results indicate that events are essentially local with a barrier height limited by the cost of breaking one single bond. We also found that there is very little variation in the entropic barrier in spite of the wide spread in the energy part. This is likely related to the fact that the overall shape of the energy landscape is not finely dependent on the specifics of a configuration. Finally, the bias found in ART nouveau indicates that it might be possible to contemplate setting up a kinetic Monte Carlo simulation on this system.

It is interesting to note the strong dissimilarities between the activation barrier distribution as computed from the initial and the final minimum. In particular, the barrier measured from the latter are much lower that those measured from the former configuration. This difference can be explained by the fact that the initial configuration selected is very-well relaxed; most low energy barrier, associated with a very unstable direction, would have been crossed during the relaxation, in previous events. These states are likely not contributing to diffusion or the relaxation of the network since they would rapidly relax back into the initial minimum. The physically relevant barrier distribution for defining the evolution of the system is therefore that measured from the initial minimum.

In spite of the 40,000+ events generated, the total number of events around a minimum remains unknown. We can nevertheless estimate that the number of different paths is somewhere between 30 and 60 paths/atom; as shown here, this number should be independent of the size of the network as all events are local in nature.

Finally the stability of the distribution from one minimum to another, as well as the narrow distribution of entropic barriers suggests that it would be possible to develop an accelerated algorithm for a–Si, similar to that proposed by Henkelman and Jónsson for the diffusion of Cu on Cu Henkelman and Jónsson (2001), and that of Hernandez-Rojas and Wales for LJ glasses. Hernandez-Rojas and Wales (2002) This is made even easier by the exponential bias of ART, also seen in LJ clusters. The origin of this bias is not understood but it simplifies significantly the statistical analysis of an ART sampling.

At this point, however, the main limitation of this type of simulation is the absence of detailed experimental results. Although the overall energy scale of the barriers is in agreement with experimental numbers, it is difficult to assess whether or not the theoretical description is correct: there is no information available experimentally on the shape of the barrier distribution. It is currently impossible to generate a better comparison with experiment. Further experimental work is clearly needed to expand on the current available data and help ensure that the advances in numerical results are on solid grounds.

## V Acknowledgments

We thank F. van Wijland for providing us with the solution do the classical occupancy problem. We also thank D. Wales for useful discussions and for providing their original data. FV is grateful to the Natural Sciences and Engineering Research Council of Canada (NSERC) for a summer scholarship. This work was partially supported by NSERC and the Fonds Nature et Technologie du Québec (NATEQ). NM is a Cottrell Scholar of the Research Corporation.


text:
# Cuts and penalties: a comment on “The clustering of ultra-high energy cosmic rays and their sources”

## I Introduction

Identification of sources of the ultra-high energy cosmic rays (UHECR) is a key to understanding their nature. The observed small-scale clustering of UHECR Hayashida:bc ; Takeda:1999sg ; Tinyakov:2001ic ; Takeda:2001xx  suggests that already existing data Takeda:1999sg ; experiments ; Hayashida:2000zr  may contain information sufficient to make first steps in this direction. In a series of papers Tinyakov:2001nr ; TT\_ICRC ; Tinyakov:st ; Tinyakov:2001ir ; Gorbunov:2002hk ; Tinyakov:2002nm ; Tinyakov:2003nu  we have shown that there exist significant correlations between the arrival directions of UHECR and BL Lacertae objects (BL Lacs), and therefore BL Lacs are likely to be sources of UHECR. Although present small dataset does not allow to determine with certainty the details of UHECR production and propagation, interesting hints may be obtained Tinyakov:2001ir ; Gorbunov:2002hk ; Tinyakov:2002nm ; Tinyakov:2003nu .

The first evidence of connection between UHECR and BL Lacs was found in Ref. Tinyakov:2001nr  where we have shown that the combined set of AGASA events with $`E>4.8\times 10^{19}`$ eV and Yakutsk events with $`E>2.4\times 10^{19}`$ eV (identified previously as the sets with largest clustering) correlates strongly with most powerful BL Lacs. The statistical significance of this correlation was found to be $`6\times 10^5`$ with penalty factor included (for explanation of penalty factor see below).

In the recent paper by Evans, Ferrer and Sarkar (EFS) Evans:2002jy  the validity of this result was called into question. The rational part of the criticism of EFS boils down to the following two issues: i) choice of cuts on BL Lacs and/or calculation of the penalty factor ii) choice of the cosmic ray set. Below we address these issues and demonstrate that criticism of EFS is unjustified. We also clarify some frequently arising questions related to our analysis.

## II Statistical significance of correlations: cuts and penalty factors

### II.1 What is penalty factor

In statistical analysis, one nearly always has to make cuts in order to improve signal-to-noise ratio. The question is how to take them into account correctly. A common wisdom is that cuts should be fixed a priori, i.e., based on theoretical considerations. In that case the cuts simply limit the number of data points, but do not alter the calculation of probabilities.

Not always this is possible. For instance, in the case of UHECRs, their acceleration mechanism is not known. How to impose cuts in a catalogue of astrophysical objects in order to select actual UHECR emitters? In which band — radio, optical, X-ray, $`\gamma `$, TeV? To which flux limit? In Ref. Tinyakov:2001nr  we proposed an approach to this problem which consists in adjusting cut(s) so as to maximize the signal, and then compensating this cut adjustment by a penalty factor.

A penalty factor takes into account the “number of independent attempts” made when searching for the best signal. For instance, if two independent catalogues were tried, the penalty factor is 2 (in the limit of small probabilities). If $`N`$ catalogues tried are not independent, as in the case of cuts in a single catalogue, the penalty factor is smaller than $`N`$ and should be calculated by means of the Monte-Carlo simulation.

In correlation analysis of cosmic rays, the quantity of interest is the probability $`p_{\mathrm{min}}^{\mathrm{data}}`$ that the observed excess of cosmic rays around source positions is the result of a chance coincidence. This probability depends on cuts made in the source catalogue; the latter are adjusted so that the probability is minimum (correlations are maximum). The cut adjustment should be compensated by the penalty factor which is calculated as follows. A random set of cosmic rays is generated and treated exactly as the real data: the same cuts are tried and the resulting minimum probability $`p_{\mathrm{min}}^{\mathrm{MC}}`$ is determined. This procedure is then repeated for a large number of random “cosmic ray” sets, and $`p_{\mathrm{min}}^{\mathrm{MC}}`$ is determined each time. The number of occurrences of a value $`p_{\mathrm{min}}^{\mathrm{MC}}p`$ is then counted as a function of $`p`$. Divided by the total number of sets, this gives the probability $`P(p)`$ that the adjustments of cuts produces $`p_{\mathrm{min}}p`$ for a random set of cosmic rays. $`P(p_{\mathrm{min}}^{\mathrm{data}})`$ is the correct measure of the significance of observed correlations.

When no adjustment of cuts is made one obviously has $`P(p)=p`$, and the quantitative measure of correlation is the probability $`p`$ itself. When cuts are adjusted, small probabilities appear more often by construction, and this relation is modified. The modification can be written in terms of the penalty factor $`F(p)`$,

$$P(p)=F(p)p.$$
(1)
In the presence of cut adjustment, the significance of correlations is determined by the product $`F(p_{\mathrm{min}}^{\mathrm{data}})p_{\mathrm{min}}^{\mathrm{data}}`$.

To summarize, all parameters characterizing UHECR sources and UHECRs themselves can be subject to cuts. Some of these cuts can be decided a priori; they reflect (and depend on) physical assumptions made. These cuts do not require penalty. One can call these cuts fixed. Alternatively, in the absence of physical arguments, the cut should be chosen (adjusted) so as to maximize the signal. These cuts can be called adjustable. These are the adjustable cuts which imply non-trivial penalty factors; the latter can be calculated in a way explained above.

### II.2 UHECR set

In our analysis Tinyakov:2001nr , the cosmic ray set was fixed on physical grounds using results of previous publications Tinyakov:2001ic , and was not adjusted in search for correlations. The motivation was as follows. It was observed earlier Hayashida:bc ; Takeda:1999sg ; Takeda:2001xx  that highest-energy cosmic rays exhibit remarkable property: they auto-correlate on the angular scale consistent with the angular resolution of detectors. The mere existence of these correlations suggests that a relatively small number of point sources of cosmic rays is contributing a considerable fraction of the UHECR flux Dubovsky:2000gv , and that the trace of the original directions to sources is not lost completely. The identification of these sources via cross-correlation analysis may therefore be possible. Within the hypothesis that UHECR clusters are due to point sources, the strongest correlation signal is expected with the UHECR set which has most significant auto-correlations. The latter requirement selects Tinyakov:2001ic  AGASA events with $`E>4.8\times 10^{19}`$ eV and Yakutsk events with $`E>2.4\times 10^{19}`$ eV. The observed angular size of clusters also fixes the angular scale at which correlations with sources are expected.

The fact that AGASA and Yakutsk sets have different energy cuts in this approach is, of course, disturbing. It may, however, be explained by different energy calibration of the two experiments, or smaller number of the Yakutsk events. In either case, within our assumptions it is inconsistent to change energy cuts or to drop Yakutsk events once autocorrelations are found in this set.

The authors of Ref. Evans:2002jy  have rejected the Yakutsk events because Yakutsk array has angular resolution not as good as AGASA. By itself, worse angular resolution does not imply that correlations with sources must be absent in Yakutsk set: even though the angular resolution is worse, the density of UHECR events around actual sources is larger as compared to a random set, and one has an excess in counts even at small angles. The actual presence of autocorrelations in Yakutsk data supports this statement. Therefore, it is not correct that “the correlations found at smaller $`\text{than Yakutsk angular resolution}`$ angles cannot be meaningful” Evans:2002jy .

### II.3 Correlations with BL Lacs

The authors of Ref. Evans:2002jy  have essentially reproduced our Monte-Carlo simulations: solid curve of Fig. 3 of Ref. Evans:2002jy  agrees with our calculations if Yakutsk events are discarded. They concluded, however, that we have underestimated the penalty factor. This conclusion was based on i) a comparison of correlations in the case when one particular set of cuts was chosen to the case when no cuts on BL Lacs were made at all ii) smaller correlations of BL Lacs with a different set of UHECR. Neither such a comparison, nor the dependence on the UHECR set have anything to do with the penalty factor. The UHECR set was fixed as described above and was not adjusted during the calculation, thus requiring no penalty. The penalty associated with the adjustment of BL Lac set should be calculated as described above. This calculation was performed in our paper Tinyakov:2001nr  as follows.

The most complete catalog of QSO veron1  contains 306 confirmed BL Lacs for which the apparent magnitude, redshift (where known) and 6 cm radio flux are listed (for several objects the 11 cm radio flux is also given). In Ref. Tinyakov:2001nr  we considered BL Lacs with $`z>0.1`$ as suggested by statistics of clustering. (We also examined a complimetary small set without further cut adjustment. This “independent attempt” adds one to the penalty factor.) We adjusted cuts on magnitude and 6 cm radio flux in the following way. First, the grid of $`11\times 11`$ cuts was fixed. Namely, the cut on magnitude was varied from $`m<16`$ to $`m<20`$ with the step $`0.4`$; simultaneously the cut on 6 cm radio flux was changed from 0.21 Jy to 0.01 Jy with the step 0.02 Jy. Note that further “refining” of this set of cuts would not change the results as the sets of BL Lacs obtained would overlap almost completely. The lowest probability (best correlation) found Tinyakov:2001nr  in the real data in this way equals to $`p_{\mathrm{min}}^{\mathrm{data}}=4\times 10^6`$ 111EFS quote a different number, $`2\times 10^5`$. This is merely misunderstanding of our paper..

To calculate the penalty factor associated with cut adjustment, a random configuration of 65 UHECR events was generated. Then its correlation with all $`11\times 11`$ subsets of BL Lacs was examined and the minimum probability $`p_{\mathrm{min}}`$ over $`11\times 11`$ cases was found and recorded. The whole procedure was repeated for the next random configuration. A total of $`N_{\mathrm{tot}}=10^5`$ random configurations were treated in this way. The number of occurrences $`N(p)`$ of $`p_{\mathrm{min}}<p`$ was calculated as a function of $`p`$. Divided by $`N_{\mathrm{tot}}`$ it determines the corrected probability $`P(p)`$ and the penalty factor $`F(p)`$ as defined in Eq. (1). The resulting penalty factor $`F(p)`$ is shown in Fig. 1.

The error bars correspond to statistical errors in determination of $`N(p)`$.

Taking the penalty factor at the minimum probability obtained for the real data, $`F(4\times 10^6)14`$ and adding 1 as penalty for BL Lacs with $`z<0.1`$ as explained above, one finds the significance of correlations including penalty,

$$P(4\times 10^6)6\times 10^5.$$

The significance is $`4\sigma `$ in terms of Gaussian distribution. It should be pointed out that, in order to get the correct significance, the penalty factor should be multiplied by the minimum probability obtained for the real data with the same set of cuts.

We stress again that the absence of correlations of cosmic rays with the entire sample of known BL Lacs does not mean anything. Physically, cosmic rays correlate with the brightest BL Lacs. Including dimmer BL Lacs in the analysis dilutes the correlations, and this is the only conclusion one draws from the analysis made in Ref. Evans:2002jy .

### II.4 “Correlations” with GRB

The authors of Ref. Evans:2002jy  have chosen GRBs as a control set where correlations are not expected; it has to be compared with the case of BL Lacs. From the fact that correlations with the whole catalogue are absent in both cases while there are many coincidences, the authors of Ref. Evans:2002jy  concluded that GRBs correlate with UHECR “just as well as do BL Lacs!”, implying that correlations with BL Lacs are due to cut adjustment. This conclusion can be tested by Monte-Carlo simulation: one has to check whether cut adjustment leads to apparent correlations with GRBs. It does not, as we now demonstrate.

The Fourth BATSE Gamma-Ray Burst Catalog Paciesas:1999tp  contains 1637 objects from which we cut out, in sequential order (according to observation date and time), 5 non-overlapping sets of 306 objects each (recall that 306 is the number of confirmed BL Lacs in the catalogue veron1 ). In each set we have looked for maximum correlations by optimizing cuts. We present two different tests: one-dimensional cut which allows transparent graphical representation, and two-dimensional cut which mimics more precisely the case of BL Lacs. We took the same maximally auto-correlated UHECR data set as in Refs. Tinyakov:2001nr ; Evans:2002jy .

In the first test we replace the $`11\times 11`$ cuts on magnitude and radio flux by 121 cut on the total number of objects (we take them in sequential order, so that first few cuts correspond to including 3, 6, 9 … first objects, and the last cut corresponds to including all 306 objects) 222It is not important on which particular parameters the cuts are imposed in the control set. A statement that some parameter is a priori important is already equivalent to the statement that correlations are present.. Like in the case of BL Lacs, we then calculate the probability of the observed excess of UHECRs around GRBs as a function of the cut. The results are presented in Fig. 2. The smallest probability found in all 5 sets is $`>4\%`$ (which by itself is not statistically significant). In view of the similar number of cuts in two cases, Fig. 1 provides an estimate 333The exact calculation of penalty factor is, of course, possible along the lines of Sect. IIA, but rough estimate is sufficient in the case at hand. for the penalty factor, $`F(0.04)5`$, in a single set of 306 GRBs. Multiplying the best probability, the penalty factor and the number of independent GRB sets one gets a number of order 1.

To mimic the case of BL Lacs we perform a similar test with two independent cuts. The first cut is done as before on the sequential number of the event. The second cut is imposed on the time (in seconds of the day, UT) of the event. Note that these two cuts are virtually uncorrelated since BATSE registered roughly one event per day. We have searched for the best correlation signal on the grid of $`11\times 11`$ equidistant cuts from 1 to 306 and from 1 to 86400, respectively. In the above five catalogs consisting of 306 GRBs each, the following lowest probabilities were found: 9.5%, 1.4%, 20%, 1.8%, 24%. Multiplying the lowest of these probabilities 0.014 by the penalty factor $`F(0.01)8`$ and by the number of independent GRB sets, we again obtain a number of order 1.

These examples are a good illustration of compensation between the effects due to cut adjustment and the penalty factor. When compensated by the penalty factor, cut selection does not introduce apparent correlations when in reality they are absent.

To avoid confusion, we note again that we used GRBs here as a control set, and imposed physically unmotivated cuts on purpose.

## III Subtleties of correlation analysis

In this section we discuss some subtleties of correlation analysis which were not mentioned in Ref. Evans:2002jy . We hope that this will answer a number frequently arising questions regarding our work and correlation analysis in general.

### III.1 Why completeness of the BL Lac catalogue is not necessary for establishing correlations with UHECR

A catalogue of astrophysical objects is called complete when it is believed to contain more than, say, 90% of all existing objects of a given type, in a given region of the sky, and down to some fixed level of luminosity in a given waveband(s). The catalogue veron1 , used in our analysis, is not complete simply because it contains all objects known up to date without further selection veron2 .

It is a common wisdom that completeness of the catalogue is crucial for statistical analysis. The reason is obvious: in an incomplete catalogue, the distribution of objects in brightness, position, etc., reflects not only their actual abundance, but also observational bias. When one studies, for instance, the evolution of abundance of some objects with time, one has to make sure that logN/logS dependence reflects actual change in spatial density rather than the fact that remote (and therefore dim) objects are more difficult to observe. Indeed, as it is, the catalogue veron1  is not suitable for any analysis where statistical properties of BL Lacs are of interest.

In the correlation analysis of UHECR the question is different: given a set of candidate sources (BL Lacs in the case at hand), is the distribution of cosmic rays random, or does it peak around BL Lac positions? Clearly, these are the statistical properties of cosmic rays which are of interest in that case, not the statistical properties of BL Lacs. Technically speaking, when calculating the correlation function by the algorithm of Ref. Tinyakov:2001nr , the cosmic ray directions are simulated, while source positions are held fixed. For this reason the method is applicable without change to studying, for instance, a correlation of UHECR with one particular direction (say, the Galactic center), i.e. in the case when the notion of completeness does not apply.

Completeness is not necessary for establishing the fact of correlations: if statistically significant correlations are found with an incomplete catalogue, this is a real signal. Incompleteness of the catalogue of BL Lacs cannot be a source of correlations with UHECR. Indeed, the objects which have to be added (removed) from a catalogue to make it complete are absent from (present in) the catalogue for reasons not related to cosmic rays. Therefore, correlations with UHECR can only be weakened by such an incompleteness.

### III.2 UHECR autocorrelations

The UHECR set used in our calculations is known to contain event clusters — in fact, it was chosen in such a way that auto-correlations are maximum. When studying cross-correlations of such a set with potential sources, one has to be careful to take clusters into account when calculating the probability of chance coincidence. Namely, one has to make sure that the signal observed in the data is not due to chance coincidence of clusters of cosmic ray events with candidate sources. For a given cluster, such a coincidence is roughly as probable as for a given single event, but contributes more into correlation function. This could produce artificial enhancement of correlations if not taken into account in the Monte-Carlo simulation.

In our calculations Tinyakov:2001nr  this problem is solved by introducing in each Monte-Carlo cosmic ray set the same number of doublets and triplets as there are in the real data. Then chance coincidences between clusters and candidate sources happen in simulated sets as often as in the real data, and thus are correctly accounted for. This is confirmed by calculations presented in Sect. II.4 where correlations with GRB catalogs do not appear despite auto-correlations in the UHECR dataset.

### III.3 Choice of angular scale

The observed angular size of clusters suggests the angular scale at which correlations with sources are expected. In our analysis Tinyakov:2001nr  we have fixed it to the previously published Hayashida:bc ; Takeda:1999sg  AGASA value of $`2.5^{}`$. This value is treated as not adjustable.

When the angular resolution of the experiment is known, there exists a preferred choice of the angular scale — the scale at which correlations are expected to be maximum. This angle can be determined by means of the Monte-Carlo simulation as follows. One has to generate cosmic ray sets which are correlated with sources (taking into account the experimental angular resolution) and then measure angle at which the correlation signal is maximum. In the case of AGASA events this procedure gives the result very close to $`2.5^{}`$ (see dotted curve on Fig. 3 of Ref. Tinyakov:2001nr ), so in fact the choice of Ref. Tinyakov:2001nr  was close to optimal.

Alternatively, one may choose not to fix the angular scale and treat it as a free parameter. Then one should adjust it to maximize the correlation and calculate corresponding penalty factor. A serious problem of this approach is that the result would depend on the limits within which the angular scale is varied. So finally one would have to input in one or the other way the information about expected angular scale of correlations in order to obtain a definite answer. We do not follow this approach in our calculations.

### III.4 Hidden penalty

We believe that we have accounted for all contributions to the penalty factor (“effective number of tries”) directly related to our work. Still, the obtained significance may be somewhat overestimated. The point is that we are not the first who are looking for correlations between UHECR and astrophysical objects. Some of these attempts have been published in the literature; the others may have never been reported. All these attempts should contribute, in principle, into the penalty factor. However, it does not seem possible to account correctly for all such contributions.

The way around this problem is obvious. First, only very low values of $`P`$ should be interpreted as a signal (in our calculations we considered $`P<10^4`$ to be sufficiently low to report our results). Second, and more important, the results have to be confirmed with a new independent data set.

## IV Further evidence

The correlations found in Ref. Tinyakov:st  for a particular set of UHECR and BL Lacs suggest, strictly speaking, only that UHECR and BL Lacs are connected to each other. They say little about the acceleration mechanism, particle nature and other relevant physical parameters. These questions can be fully addressed only with much larger dataset than available now. Some attempts, however, can be made. The final purpose is to arrive at the end at a consistent picture which incorporates all known features of UHECR. Here is a sketch of these attempts, each bringing additional evidence of the connection between UHECR and BL Lacs.

In Ref. Gorbunov:2002hk  it was noted that cuts in BL Lac catalogue, chosen in Ref. Tinyakov:2001nr  so as to maximize correlations with UHECR, select automatically $`\gamma `$-ray loud BL Lacs. When this observation is consistently elaborated, i.e. a sub-sample of the catalogue veron2  is selected on the basis of a single criterium, namely, the cross-correlation with the EGRET sources, the resulting subset of 14 $`\gamma `$-ray loud BL Lacs correlates with UHECR at the level of $`10^7`$ of chance coincidence. This number cannot, of course, be interpreted as the significance of correlations between UHECR and BL Lacs because of a posteriori selection; rather the conclusion is that $`\gamma `$-ray loudness may be a distinctive feature of those BL Lacs which are UHECR accelerators.

In Ref. Tinyakov:2001ir  an attempt was made to determine the charge composition of UHECR by reconstructing actual arrival directions of UHECR particles bent in the Galactic magnetic field. The idea was that such a reconstruction should improve correlations of UHECR with BL Lacs if the latter are the sources. Substantial improvement was indeed observed for particles with the charge +1, which is an indication of the presence of protons.

In Ref. Tinyakov:2002nm  it was observed that if correlations of UHECR and BL Lacs were due to chance coincidence, the coinciding rays would be distributed over the sky randomly, reflecting only the local density of BL Lacs and exposure of a cosmic ray experiment. Thus, any significant deviation in the distribution of correlating rays over the sky from this expectation speaks in favor of real physical connection between cosmic rays and BL Lacs. In fact, the UHECRs correlating with BL Lacs form two “spots”, with low probability to occur by chance Tinyakov:2002nm . This non-uniformity of the distribution of correlating rays may be due to several factors: (1) anisotropy of extragalactic magnetic fields at scales of order $`500`$ Mpc; (2) poor knowledge of the Galactic magnetic field in some areas of the sky; (3) fluctuations in the space distribution of the nearest sources.

## V Conclusions

According to textbooks, any successful statistical analysis consists in formulation of a “null hypothesis” and its subsequent falsification, at some confidence level, by comparing to the experimental data. In the case at hand the “null hypothesis” which is being tested is that BL Lacs (and, in particular, any subset of them) and UHECR are uncorrelated. It takes only one counter-example to disprove a hypothesis, while “pro-examples” do not prove its validity. An illustration of this general rule has been discussed in Sect. II: the absence of significant correlations of UHECR with the whole BL Lac catalogue does not prove that UHECR and BL Lacs are uncorrelated (i.e., that there is no subset of BL Lacs which are sources of UHECR and thus correlate with them).

Having found such counter-example (i.e. the case when correlation is significant), the only thing which can be concluded, to a certain confidence level, is that UHECR and a particular subset of BL Lacs are correlated. The nature and physical implications of these correlations have to be studied separately by formulating and testing different “null hypothesis”. Refs.Tinyakov:2001ir ; Gorbunov:2002hk ; Tinyakov:2002nm  are first attempts in this direction.

###### Acknowledgements.

We are grateful to V. Rubakov, D. Semikoz, M. Shaposhnikov and P. Veron for reading the manuscript and valuable comments. The work of P.T. is supported in part by the Swiss Science Foundation, grant 20-67958.02.

## Appendix A Statistics of clustering and the number of UHECR sources.

The authors of Ref. Evans:2002jy  have misinterpreted our earlier paper Dubovsky:2000gv  which concerns the statistics of clustering of UHECR and the number of their sources. In view of the growing confusion we would like to clarify this issue.

In Ref. Evans:2002jy  one reads: “$``$observed occurrence of clusters of UHECR$``$ …was used to estimate the spatial density of sources to be $`6\times 10^3\mathrm{Mpc}^3`$ Dubovsky:2000gv . This would obviously place stringent constraints on candidate astrophysical sources, e.g. $`\gamma `$-ray bursts (GRBs) have a spatial density of only $`10^5\mathrm{Mpc}^3`$. However, a more careful analysis Fodor:2001yi  shows that the uncertainties in this estimate are very large. The true number is $`n=2.77_{2.53(2.70)}^{+96.1(916)}\times 10^3\mathrm{Mpc}^3`$ at the 68% (95%) c.l.; moreover relaxing the assumptions made, viz. that the sources all have the same luminosity and a spectrum $`E^2`$, increases the allowed range even further, e.g. to $`n=180_{165(174)}^{+2730(8817)}\times 10^3\mathrm{Mpc}^3`$ …”.

The actual situation is different from what is described in this extract. First, intrinsic inaccuracy of the estimate was fully acknowledged in Ref. Dubovsky:2000gv . Its physical reason is obvious: clustering is not sensitive to the number of dim sources. The estimates of Ref. Fodor:2001yi  cited above are a good illustration of the point: upper limits are huge and strongly model-dependent.

Second, the main point of Ref. Dubovsky:2000gv  was to show that there exists a model-independent lower bound on the number of sources. This bound is presented in Table 1 of Ref. Dubovsky:2000gv : at 90% and 99% confidence levels the number of sources is $`n>2.3\times 10^4\mathrm{Mpc}^3`$ and $`n>3.2\times 10^5\mathrm{Mpc}^3`$, respectively. This is consistent with the calculations of Ref. Fodor:2001yi  performed in particular models.