text:
# Nature of phase transitions in a probabilistic cellular automaton with two absorbing states

## 1 Introduction

Probabilistic cellular automata (PCA) have been widely used to model a variety of systems with local interactions in physics, chemistry, biology and social sciences . Moreover, PCA are simple and interesting models that can be used to investigate fundamental problems in statistical mechanics. Many classical equilibrium spin models can be reformulated as PCA, for example the kinetic Ising model with parallel heat-bath dynamics is strictly equivalent to a PCA with local parallel dynamics . On the other hand, PCA can be mapped to spin models  by expressing the transition probabilities as exponentials of a local energy. PCA can be used to investigate nonequilibrium phenomena, and in particular the problem of phase transitions in the presence of absorbing states. An absorbing state is represented by a set of configurations from which the system cannot escape, equivalent to an infinite energy well in the language of statistical mechanics. A global absorbing state can be originated by one or more local transition probabilities which take the value zero or one, corresponding to some infinite coupling in the local energy .

The Domany-Kinzel (DK) model is a boolean PCA on a tilted square lattice that has been extensively studied . Let us denote the two possible states of each site with the terms “empty” and “occupied”. In this model a site at time $`t`$ is connected to two sites at time $`t1`$, constituting its neighborhood. The control parameters of the model are the local transition probabilities that give the probability of having an occupied site at a certain position once given the state of its neighborhood. The transition probabilities are symmetric for all permutations of the neighborhood, and this property is equivalent to saying that they depend on the sum of “occupied” sites in the neighborhood, whence the term “totalistic” used to denote this class of automata.

In the DK model the transition probability from an empty neighborhood to an occupied state is zero, thus the empty configuration is an absorbing state. For small values of the other transition probabilities, any initial configuration will evolve to the absorbing state. For larger values, a phase transition to an active phase, represented by an ensemble of partially occupied configurations, is found. The order parameter of this transition is the asymptotic average fraction of occupied sites, which we call the density. The critical properties of this phase transition belong to the directed percolation (DP) universality class (except for one extreme point) , and the DK model is often considered the prototype of such a class.

The evolution of this kind of models is the discrete equivalent of the trajectory of a stochastic dynamical system. One can determine the sensitivity with respect to a perturbation, by studying the trajectories originating by two initially different configurations (replicas) evolving with the same realization of the stochasticity, e.g. using the same sequence of random numbers. The order parameter here is the asymptotic difference between the two replicas, which we call the damage. It turns out that, inside the active phase, there is a “chaotic” phase in which the trajectories depend on the initial configurations and the damage is different from zero, and a non chaotic one in which all trajectories eventually synchronize with the vanishing of the damage. In simple models like the DK one, this transition does not depend on the choice of the initial configurations (provided they are different from the absorbing state) and the initial damage .

It has been conjectured that all second-order phase transitions from an “active” phase to a non degenerate, quiescent phase (generally represented by an absorbing state) belong to the DP universality class if the order parameter is a scalar and there are no extra symmetries or conservation laws . This has been verified in a wide class of models, even multi-component, and in the presence of several asymmetric absorbing states . Also the damage phase transition has a similar structure. Once synchronized, the two replicas cannot separate, and thus the synchronized state is absorbing. Indeed, numerical simulations shows that it is in the DP universality class . Moreover, in the DK model, the damage phase transition can be mapped onto the density one .

On the other hand, some models with conserved quantities  or symmetric absorbing states belong to a different universality class called parity conservation (PC) or directed Ising . This universality class appears to be less robust since it is strictly related to the symmetry of the absorbing states; a slight asymmetry is sufficient to bring the model to the usual DP class .

An interesting question concerns the simplest, one-dimensional PCA model with short range interactions exhibiting a first order phase transition. Dickman and Tomé  proposed a contact process with spontaneous annihilation, auto catalytic creation by trimers and hopping. They found a first order transition for high hopping probability, i.e., in the region more similar to mean field (weaker spatial correlations).

Bassler and Browne discussed a model whose phase diagram also presents first and second-order phase transitions . In it, monomers of three different chemical species can be adsorbed on a one-dimensional surface and neighboring monomers belonging to different species annihilate instantaneously. The control parameters of the model are the absorption rates of the monomers. The transition from a saturate to a reactive phase belongs to the DP universality class, while the transition between two saturated phases is discontinuous. The point at which three phase transition lines join does belong to the PC universality class.

Scaling and fluctuations near first order phase transitions are also an interesting subject of study , which can profit from the existence of simple models.

In this paper we study a one-dimensional, one-component, totalistic PCA with two absorbing states. It can be considered as a natural extension of the DK model to a lattice in which the neighborhood of a site at time $`t`$ contains the site itself and its two nearest neighbors at time $`t1`$. This space-time lattice arises naturally in the discretization of one-dimensional reaction-diffusion systems. In our model, the transition probabilities from an empty neighborhood is zero, and that from a completely occupied neighborhood is one. The model has two absorbing states: the completely empty and the completely occupied configurations. The order parameter is again the density; it is zero or one in the two quiescent phases, and assumes other values in the active phase. The system presents a line of symmetry in the phase diagram, over which the two absorbing phases have the same importance. A more detailed illustration of the model can be found in Section 2.

This model can arise as a particular case of a nonequilibrium wetting of a surface. In this framework, only a single layer of particles can be absorbed on the surface. If we assume that particles can be absorbed or desorbed only near the boundaries of a patch of already absorbed particles (when the neighborhood is not homogeneous), then the completely empty and occupied configurations are absorbing states.

This totalistic PCA can also be interpreted as a simple model of opinion formation. It assumes that an individual may change his mind according to himself and his two nearest neighbors. The role of social pressure is twofold. If there is homogeneity of opinions, individuals cannot disagree (absorbing states), otherwise they can agree or disagree with the majority with a certain probability.

The density phase diagram shows two second order phase transition curves separating the quiescent phases from the active one, and a first order transition line between the two quiescent phases, as discussed in Section 3. These curves meet on the line of symmetry in a bicritical point. We use both mean field approximations and direct numerical simulations. The former simple approximation gives a qualitatively correct phase diagram. The numerical experiments are partially based on the fragment method . This is a parallel algorithm that implements directly the evolution rule for different values of the control parameters on the bits of one or more computer words.

In Section 4, we investigate numerically the second order phase transitions and find they belong to the DP universality class. Along the line of symmetry of the model the two absorbing phases are equivalent. In Appendix B we show that on this line one can reformulate the problem in terms of the dynamics of kinks between patches of empty and occupied sites. Since the kinks are created and annihilated in pairs, the dynamics conserves the initial number of kinks modulo two. In this way we can present an exact mapping between a model with symmetric absorbing phases and one with parity conservation. We find that the critical kink dynamics at the bicritical point belongs to the PC universality class.

In Section 5 we study the chaotic phase, using dynamic mean field techniques (reported in Appendix A) and direct numerical simulations. The location of this phase is similar to that of the DK model: it joins the second-order critical curves at the boundary of the phase diagram.

Our model exhibits a first-order phase transition along the line of symmetry in the upper part of the phase diagram. A first-order transition is usually associated to an hysteresis cycle. It is possible to observe such a phenomena by adding a small perturbing field to the absorbing states, as discussed in Section 6.

The DP universality class is equivalent to the Reggeon field theory , which in $`d=0`$ corresponds to a quadratic potential with a logarithmic divergence at the origin. The Langevin description for systems in the PC class yields a similar potential, except for irrelevant terms . It has been shown  that one can reconstruct the potential from the numerical integration of the Langevin equation, which, however, requires special techniques in the presence of absorbing states . In Section 7 we show how the potential is reconstructed from actual simulations of a phenomenological model, such as our original cellular automata or the kink dynamics. In this way we obtain the shape of the potential for a system in the parity conservation universality class.

## 2 The model

We describe here a one-dimensional, totalistic, probabilistic cellular automaton with three inputs. The state of the model at time $`t`$ is given by $`𝒙^t=(x_0^t,\mathrm{},x_{L1}^t)`$ with $`x_i^t\{0,1\}`$; $`t=1,2,\mathrm{}`$ and $`L`$ is the number of sites. All operations on spatial indices are assumed to be modulo $`L`$ (periodic boundary conditions). For simplicity of notation, we write $`x=x_i^t`$, $`x_{}=x_{i1}^t`$, $`x_+=x_{i+1}^t`$ and $`x^{}=x_i^{t+1}`$. We shall indicate by $`\sigma =x_{}+x+x_+`$ number of occupied cells in the neighborhood. The most general three-input totalistic PCA is defined by the quantities $`p_s`$, which are the conditional probabilities that $`x^{}=1`$ if $`\sigma =s`$. The microscopic dynamics of the model is completely specified by

$$x^{}=\underset{s=0}{\overset{3}{}}R_s\delta _{\sigma ,s}.$$
(1)
In this expression $`R_s`$ is a stochastic binary variable that takes the value 1 with probability $`p_s`$ and 0 with probability $`1p_s`$, and $`\delta `$ is the Kronecker delta. In practice, we implement $`R_s`$ by extracting a random number $`r_s`$ uniformly distributed between 0 and 1, and setting $`R_s`$ equal to 1 if $`r_s<p_s`$ and 0 otherwise. Eq. (1) implies the use of four random numbers $`r_s`$ for each site. The evolution of a single trajectory is not affected by the eventual correlations among the $`r_s`$, since only one $`\delta _{\sigma ,s}`$ is different from zero. This is not true when computing the simultaneous evolution of two or more replicas using the same noise. If not otherwise stated, we use only one random number for all the $`r_s`$. More discussions about the choice of random numbers can be found in Section 5 and in Appendix B.

With $`p_0=0`$ and $`p_3=1`$ the model presents two quiescent phases: the phase 0 corresponding to the configuration $`x=(0,\mathrm{},0)`$ and the phase 1 corresponding to the configuration $`x=(1,\mathrm{},1)`$. In this case there are two control parameters, $`p_1`$ and $`p_2`$, and the model is symmetric under the changes $`p_11p_2`$, $`p_21p_1`$ and $`xx1`$ where $``$ is the exclusive disjunction (or the sum modulo 2).

## 3 Phase diagram

In order to have a qualitative idea of the behavior of the model, we first study the mean-field approximation. If $`c`$ and $`c^{}`$ denote the density of occupied sites at times $`t`$ and $`t+1`$ respectively,

$$\begin{array}{cc}\hfill c^{}& =3p_1c\left(1c\right)^2+3p_2c^2\left(1c\right)+c^3.\hfill \end{array}$$
(2)
This map has three fixed points, $`c_0`$, $`c_1`$, and $`c_2`$ given by

$$c_0=0,c_1=\frac{3p_11}{1+3p_13p_2},\text{and}c_2=1.$$

The asymptotic density will assume one of the latter values according to the values of the control parameters and the initial state as we show in Fig. 1. In the square $`1/3<p_11`$, $`0p_2<2/3`$, the only stable fixed point is $`c_1`$. Inside this square, on the segments $`p_22/3=m(p_11/3)`$ with $`m<0`$, $`c_1=1/(1m)`$.The first fixed point $`c_0`$ is stable when $`p_1<1/3`$ and $`c_2`$ is stable when $`p_2>2/3`$. There is a continuous second-order transition from the quiescent phase 0 to the active phase on the segment $`p_1=1/3`$, $`0p_2<2/3`$ and another continuous transition from the active to the quiescent phase 1 on the segment $`1/3<p_11,p_2=2/3`$.

In the hatched region of Fig. 1 $`c_0`$ and $`c_2`$ are both stable. Their basins of attraction are, respectively, the semi-open intervals $`[0,c_1)`$ and $`(c_1,1]`$. Starting from a uniformly distributed random value of $`c`$, as time $`t`$ goes to infinity, $`c`$ tends to $`c_0`$ with probability $`c_1`$, and to $`c_2`$ with probability $`1c_1`$. Since, for $`p_1+p_2=1`$, $`c_1=1/2`$, the segment $`p_1+p_2=1`$, with $`0p_1<1/3`$ and $`2/3<p_21`$, is similar to a first-order transition line between the phase 0 and the phase 1.

In Fig. 2 we show the phase diagram of the model obtained numerically starting from a random initial state with half of the cells occupied. The scenario is qualitatively the same as predicted by the mean-field analysis. In the vicinity of the point $`(p_1,p_2)=(0,1)`$ we observe a discontinuous transition from $`c=0`$ to $`c=1`$. The two second-order phase-transition curves from the active phase to the quiescent phases are symmetric, and the critical behavior of the order parameter, $`c`$ for the lower curve and $`1c`$ for the upper one, is the same.

Due to the symmetry of the model the two second-order phase transition curves meet at a bicritical point $`(p_t,1p_t)`$ where the first-order phase transition line $`p_1+p_2=1`$, $`p_1<p_t`$ ends. Crossing the second-order phase boundaries on a line parallel to the diagonal $`p_1=p_2`$, the density $`c`$ exhibits two critical transitions, as shown in the inset of Fig. 2. Approaching the bicritical point the critical region becomes smaller, and corrections to scaling increase. Finally, at the transition point the two transitions coalesce into a single discontinuous one.

## 4 Critical dynamics and universality classes

We performed standard dynamic Monte Carlo simulations starting from a single site in the origin out of the nearest absorbing state, and measured the average number of active sites $`N(t)`$, the survival probability $`P(t)`$ and the average square distance from origin $`R^2(t)`$ (averaged over surviving runs) defined as

$$\begin{array}{cc}\hfill N(t)& =\frac{1}{K}\underset{k=1}{\overset{K}{}}\underset{i}{}\omega _i^t(k),\hfill \\ \hfill P(t)& =\frac{1}{K}\underset{k=1}{\overset{K}{}}\theta \left(\underset{i}{}\omega _i^t(k)\right),\hfill \\ \hfill R^2(t)& =\frac{1}{KN(t)}\underset{k=1}{\overset{K}{}}\underset{i}{}\omega _i^t(k)i^2.\hfill \end{array}$$
(3)
In these expressions, $`k`$ labels the different runs and $`K`$ is the total number of runs. The quantity $`\omega _i`$ is $`x_i`$ if the nearest absorbing state is $`\mathrm{𝟎}`$ and $`1x_i`$ otherwise; $`\theta `$ is the Heaviside step function, that assumes the value 1 if its argument is greater than 0, and the value 0 if it is smaller than 0.

At the critical point one has

$$N(t)t^\eta ,P(t)t^\delta ,R^2(t)t^z.$$

At the transition point $`(p_1^{}=0.6625(3),p_2=0)`$, we get $`\eta =0.308(5)`$, $`\delta =0.160(2)`$ and $`z=1.265(5)`$, in agreement with the best known values for the directed percolation universality class .

Near the bicritical point, on the line $`p_1+p_2=1`$, the two absorbing states have symmetrical weight. We define a kink $`y_i`$ as $`y_i=x_ix_{i+1}`$. For the computation of the critical properties of the kink dynamics, one has to replace $`\omega _i`$ with $`y_i`$ in Eq. (3). The evolution equation is derived in Appendix B. In the kink dynamics there is only one absorbing state (the empty state), corresponding to one of the two absorbing states $`\mathrm{𝟎}`$ or $`\mathrm{𝟏}`$. For $`p_1<p_t`$ the asymptotic value of the density of kinks is zero and it starts to grow for $`p_1>p_t`$. In models with multiple absorbing states, dynamical exponents may vary with initial conditions. Quantities computed only on survival runs ($`R^2(t)`$) appear to be universal, while others (namely $`P(t)`$ and $`N(t)`$) are not .

We performed dynamic Monte Carlo simulations starting either from one and two kinks. In both cases $`p_t=0.460(2)`$, but the exponents were found to be different. Due to the conservation of the number of kinks modulo two, starting from a single site one cannot observe the relaxation to the absorbing state, and thus $`\delta =0`$. In this case $`\eta =0.292(5)`$, $`z=1.153(5)`$. On the other hand, starting with two neighboring kinks, we find $`\eta =0.00(2)`$, $`\delta =0.285(5)`$, and $`z=1.18(2)`$. These results are consistent with those found by other authors .

## 5 The chaotic phase

Let us now turn to the sensitivity of the model to a variation in the initial configuration, i.e. to the study of damage spreading or, equivalently, to the location of the chaotic phase.

Given two replicas $`x`$ and $`y`$, we define the difference $`w`$ as $`w=xy`$. The damage $`h`$ is defined as the fraction of sites in which $`w=1`$, i.e. as the Hamming distance between the configurations $`x`$ and $`y`$.

The precise location of this phase transition depends on the particular implementation of the stochasticity. Since the sum of occupied cells in the neighborhood of $`x`$ is in general different from that of $`y`$, the evolution equation (1) for the two replicas uses different random numbers $`r_s`$. The correlations among these random numbers affect the location of the chaotic phase boundary .

We limit our investigation to the case of maximal correlations by using just one random number per site, i.e. all $`r_s`$ are the same for all $`s`$ at the same site. This gives the smallest possible chaotic region. Notice that we have to extract a random number for all sites, even if in one or both replicas have a neighborhood configuration for which the evolution rule is deterministic (if $`\sigma =0`$ or $`\sigma =3`$).

One can write the mean field equation for the damage by taking into account all the possible local configurations of two lattices. The evolution equation for the damage depends on the correlations among sites but in the simplest case we can assume that $`h(t+1)`$ depends only on $`h(t)`$ and $`c(t)`$, the density of occupied sites. In Appendix A we find the evolution equation for the damage in the mean-field approximation. In Fig. 3 we show the phase diagram of the chaotic phase in this approximation. There is a qualitative agreement with the mean field phase diagram found for the DK model .

In Fig. 4 we show the phase diagram for the damage found numerically by considering the evolution starting from uncorrelated configurations with initial density equal to 0.5. The damage region is shown in shades of grey. Outside this region there appear small damaged domains on the other phase boundaries. This is due either to the divergence of the relaxation time (second-order transitions) or to the fact that a small difference in the initial configuration can drive the system to a different absorbing state (first-order transitions). The chaotic domain near the point $`(p_1,p_2)=(1,0)`$ is stable regardless of the initial density. On the line $`p_2=0`$ the critical points of the density and the damage coincide at $`p_1^{}`$.

## 6 First-order phase transition and hysteresis cycle

First-order phase transitions are usually associated to a hysteresis cycle due to the coexistence of two phases. In the absence of absorbing states, the coexistence of two stable phases for the same values of the parameters is a transient effect in finite systems, due to the presence of fluctuations. To find the hysteresis loop we modify the model slightly by putting $`p_0=1p_3=\epsilon `$ with $`\epsilon 1`$. In this way the empty and fully occupied configurations are no longer absorbing. This brings the model back into the class of equilibrium models for which there is no phase transition in one dimension but metastable states can nevertheless persist for long times. The mean field equation for the density $`c`$ becomes

$$\begin{array}{cc}\hfill c^{}=& \epsilon (1c)^3+3p_1c\left(1c\right)^2+\hfill \\ & 3p_2c^2\left(1c\right)+(1\epsilon )c^3.\hfill \end{array}$$
(4)
We study the asymptotic density as $`p_1`$ and $`p_2`$ move on a line with slope 1 inside the hatched region of Fig. 1. For $`p_1`$ close to zero, Eq. (4) has only one fixed point, which is stable and close to $`\epsilon `$. As $`p_1`$ increases adiabatically (by taking $`c`$ at $`t=0`$ equal the previous value of the fixed point) the new asymptotic density will still assume this value even when two more fixed points appear one of which is unstable and the other stable and close to one. Eventually the first fixed point disappears, and the asymptotic density jumps to the stable fixed point close to one. Going backwards on the same line, the asymptotic density will be close to one until that fixed point disappears and it will jump back to a small value close to zero. By symmetry, the hysteresis loop is centered around the line $`p_1+p_2=1`$ which we identify as a first-order phase transition line inside the hatched region.

The hysteresis region is found by two methods, the dynamical mean field, which extends the mean field approximation to blocks of $`l`$ sites , and direct numerical experiments. As stated before it is necessary to introduce a small perturbation $`\epsilon =p_0=1p_3`$. We consider lines parallel to the diagonal $`p_1=p_2`$ in the parameter space and increase the value of $`p_1`$ and $`p_2`$ after a given relaxation time $`t_r`$ up to $`p_2=1`$; afterwards the scanning is reverted down to $`p_1=0`$. The hysteresis region for various values of parameters are reported Fig. 5. In numerical simulations, one can estimate the size of the hysteresis region by starting with configurations $`\mathrm{𝟎}`$ and $`\mathrm{𝟏}`$, and measuring the size $`d`$ of the region in which the two simulations disagree.

## 7 Reconstruction of the potential

An important point in the study of systems exhibiting absorbing states is the formulation of a coarse-grained description using a Langevin equation. It is generally accepted that DP universal behavior is represented by

$$\frac{c(x,t)}{t}=ac(x,t)bc^2(x,t)+^2c(x,t)+\sqrt{c(x,t)}\alpha (x,t),$$
(5)
where $`c`$ is the density field, $`a`$ and $`b`$ are control parameters and $`\alpha `$ is a Gaussian noise with correlations $`\alpha (x,t)\alpha (x^{},t^{})=\delta _{x,x^{}}\delta _{t,t^{}}`$. The diffusion coefficient has been absorbed into the parameters $`a`$ and $`b`$ and the time scale. This equation can be obtained by a mean-field approximation of the evolution equation keeping only the relevant terms. The state $`c(x,t)=0`$ is clearly stationary, but its absorbing character is given by the balance between fluctuations, which are of order of the field itself, and the “potential” part $`acbc^2`$ (See also Ref. ).

The role of the absorbing state can be illustrated by taking a sequence of equilibrium models whose energy landscape exhibits the coexistence of an infinitely deep well (the absorbing state) and another broad local minimum (corresponding to the “active”, disordered state), separated by an energy barrier. There is no true stationary active state for such a system (with a finite energy barrier), since there is always a probability of jumping into the absorbing state. However, the system can survive in a metastable active state for time intervals of the order of the inverse of the heigth of the energy barrier. The parameters controlling the height of the energy barrier are the size of the lattice and the length of the simulation: the equilibrium systems are two-dimensional with asymmetric interactions in the time direction . In the limiting case of an infinite system, the height of the energy barrier is finite below the transition point, and infinite above, when the only physically relevant state is the disordered one.

It is possible to introduce a zero-dimensional approximation to the model by averaging over the time and the space, assuming that the system has entered the metastable state. In this approximation, the size of the original systems enters through the renormalized coefficients $`\overline{a}`$, $`\overline{b}`$,

$$\frac{c(x,t)}{t}=\overline{a}c(x,t)\overline{b}c^2(x,t)+\sqrt{c(x,t)}\alpha (x,t),$$

where also the time scale has been renormalized.

The associated Fokker-Planck equation is

$$\frac{P(c,t)}{t}=\frac{}{c}(\overline{a}c\overline{b}c^2)P(c,t)+\frac{1}{2}\frac{^2}{c^2}cP(c,t),$$

where $`P(c,t)`$ is the probability of observing a density $`c`$ at time $`t`$. One possible solution is a $`\delta `$-peak centered at the origin, corresponding to the absorbing state.

By considering only those trajectories that do not enter the absorbing state during the observation time, one can impose a detailed balance condition, whose effective agreement with the actual probability distribution has to be checked *a posteriori*.

A stationary distribution $`P(c)=\mathrm{exp}(V(c))`$ corresponds to an effective potential $`V(c)`$ of the form

$$V(c)=\mathrm{log}(c)2\overline{a}c+\overline{b}c^2.$$

Note that this distribution is not normalizable. One can impose a cutoff for low $`c`$, making $`P(c)`$ normalizable. For finite systems the only stationary solution is the absorbing state. However, by increasing the size of the system, one approximates the limit in which the energy barrier is infinitely hight, the absorbing state unreachable and $`P(c)`$ is the observable distribution.

In order to find the form of the effective potential for spatially extended systems, Muñoz  numerically integrated Eq. (5), using the procedure described by Dickman . It is however possible to obtain the shape of the effective potential from the actual simulations, simply by plotting $`V(c)=\mathrm{log}(P(c))`$ versus $`c`$, where $`c`$ is the density of the configuration.

In Fig. 6 we show the profile of the reconstructed potential $`V`$ for some values of $`p`$ around the critical value or the infinite system $`p_1^{}`$ on the line $`q=0`$. We used rather small systems and followed the evolution for a limited amount of time in order to balance the weight of the $`\delta `$-peak with respect to $`P(c)`$ (which is only metastable). For larger systems the absorbing state is not visible above the transition and dominates below it.

On the line $`q=0`$ the model belongs to the DP universality class. One can observe that the curve becomes broader in the vicinity of the critical point, in correspondence of the divergence of critical fluctuations $`\chi |pp_c|^\gamma ^{}`$, $`\gamma ^{}=0.54`$ . By repeating the same type of simulations for the kink dynamics (random initial condition), we obtain slightly different curves, as shown in Fig. 7. We notice that all curves have roughly the same width. Indeed, the exponent $`\gamma ^{}`$ for systems in the PC universality class is believed to be exactly 0 , as given by the scaling relation  $`\gamma ^{}=d\nu _{}2\beta `$. Clearly, much more informations can be obtained from the knowledge of $`P(c)`$, either by direct numerical simulations or dynamical mean field trough finite scale analysis, as shown for instance in Ref. .

## 8 Discussions and conclusions

We have studied a probabilistic cellular automata with two absorbing states and two control parameters. This is a simple and natural extension of the Domany-Kinzel (DK) model. Despite its simplicity it has a rich phase diagram with two symmetric second-order phase curves that join a first-order line at a bicritical point. The phase diagram and the critical properties of the model were found using several mean field approximations and numerical simulations. The second-order phase transitions belong to the directed percolation universality class except for the bicritical point, which belongs to the parity conservation (PC) universality class. The first-order phase transition line was put in evidence by a modification of the model that allows one to find the hysteresis cycles. The model also presents a chaotic phase analogous to the one present in the DK model. This phase was studied using direct numerical simulations and dynamical mean field.

On the line of symmetry of the model the relevant behavior is given by kink dynamics. We found a closed expression for the kink evolution rule and studied its critical properties, which belong to the PC universality class. The effective potential governing the coarse-grained evolution for the DP and the PC phase was found through direct simulations, confirming that critical fluctuations diverge at most logarithmically in the PC class.

The phase diagram of our model is qualitatively similar to Bassler and Browne’s (BB) one . In both models two critical lines in the DP universality class meet at a bicritical point in the PC universality class, and give origin to a first-order transition line. This suggests that the observed behavior has a certain degree of universality.

An interesting feature of the BB model is that the absorbing states at the bicritical point are indeed symmetric, but the model does not show any conserved quantities. We have shown that the bicritical dynamics of our model can be exactly formulated either in terms of symmetric states or of kinks dynamics, providing an exact correspondence between the presence of conserved quantities and the symmetry of absorbing states.

Furthermore, in order to obtain a qualitatively correct mean-field phase diagram of the BB model, one has to include correlations between triplets, while the mean-field phase diagram of our model is already correct at first approximation. This suggests that we have described a simpler model, which can be used as prototype for multi-critical systems.

## Acknowledgements

Helpful and fruitful discussions with Paolo Palmerini, Antonio Politi and Hernán Larralde are acknowledged. This work benefitted from partial economic support from CNR (Italy), CONACYT (Mexico), project IN-116198 DGAPA-UNAM (Mexico) and project G0044-E CONACYT (Mexico). We wish to thank one of the Referees for having indicated us Ref. .

## Appendix A Damage spreading in the mean field approximation

The minimum damage spreading occurs when the two replicas $`𝒙`$ and $`𝒚`$ evolve using maximally correlated random numbers, i.e. when all $`r_s`$ in Eq. (1) are the same. . Let $`w=xy`$ be the damage at a site $`i`$ and time $`t`$. It is also possible to consider $`w`$ as an independent variable and write $`y=xw`$. We denote $`s=x_{}+x+x_+`$, $`s^{}=y_{}+y+y+=(x_{}w_{})+(xw)+(x_+w_+)`$ and $`s^{\prime \prime }=w_{}+w+w_+`$. The evolution equation for $`h`$, the density of damaged sites $`w`$ at time $`t`$, is obtained by considering all the local configurations $`x_{}xx_+`$ and $`w_{}ww_+`$ of one replica and of the damage

$$h^{}=\underset{\begin{array}{c}x_{}xx_+\\ w_{}ww_+\end{array}}{}\pi (c,s,3)\pi (h,s^{\prime \prime },3)\left|p_sp_s^{}\right|,$$
(A.1)
where

$$\pi (\alpha ,n,m)=\alpha ^n(1\alpha )^{nm}.$$

In this equation all the sums run from zero to one. The value of $`c`$ is given by Eq. (2). The term $`|p_sp_s^{}|`$ is the probability that $`R_sR_s^{}`$ is one using only one random number for the $`r_s`$. The argument of the sum is the probability that $`x^{}y^{}`$. It is possible to rewrite Eq.(A.1) in a different form

$$h^{}=\underset{ss^{}\mathrm{}}{}\left(\genfrac{}{}{0pt}{}{m}{s}\right)\left(\genfrac{}{}{0pt}{}{ms}{\mathrm{}}\right)\left(\genfrac{}{}{0pt}{}{ms}{s^{}\mathrm{}}\right)\pi (c,s,m)\pi (\eta ,s+s^{}2\mathrm{},m)|p_sp_s^{}|,$$
(A.2)
where $`s`$ and $`s^{}`$ are the same as above, and $`\mathrm{}`$ is the overlap between $`x`$ and $`y`$, i.e. $`\mathrm{}=(x_{}y_{})+(xy)+(x_+y_+)`$ ($``$ is the AND operation). Assuming that $`\left(\genfrac{}{}{0pt}{}{a}{b}\right)=0`$ if $`b>a`$, $`a<0`$ or $`b<0`$, the sum of (A.2) can run over all positive integers. This expression is valid for all totalistic rules with a neighborhood of size $`m`$ (here $`m=3`$).

The stationary state of Eq. (A.1) (or Eq. (A.2)) can be found analytically using a symbolic manipulation program. The chaotic transition line is

$$p_2=p_1\frac{1}{9},$$

with $`1/3<p_1<1`$,$`0<p_2<2/3`$.

## Appendix B Kink dynamics

On the segment $`p_1+p_2=1`$, $`p_1<p_t`$ the order parameter is the number of kinks. The dynamics of the kinks $`y_i=x_ix_{i+1}`$ (that for the ease of notation we write $`y=xx_+`$) is obtained by taking the exclusive disjunction of $`x^{}=x_i^{t+1}`$ and $`x_+^{}=x_{i+1}^{t+1}`$ given by Eq. (1). In order to obtain a closed expression for the $`y`$, a little of Boolean algebra is needed.

The totalistic functions $`\delta _{\sigma ,s}`$ where $`s=x_{}+x+x_+`$ can be expressed in terms of the symmetric polynomials $`\xi ^{(j)}`$ of degree $`j`$ . These are

$$\begin{array}{cc}\hfill \xi ^{(1)}& =x_{}xx_+,\hfill \\ \hfill \xi ^{(2)}& =x_{}xx_{}x_+xx_+,\hfill \\ \hfill \xi ^{(3)}& =\xi ^{(1)}\xi ^{(2)}=x_{}xx_+.\hfill \end{array}$$

The totalistic functions are given by

$$\begin{array}{cc}\hfill \delta _{\sigma ,1}& =\xi ^{(1)}\xi ^{(3)},\hfill \\ \hfill \delta _{\sigma ,2}& =\xi ^{(2)}\xi ^{(3)},\hfill \\ \hfill \delta _{\sigma ,3}& =\xi ^{(3)}.\hfill \end{array}$$

In the evolution equation (1), one has $`R_1=1`$ if $`r_1<p_1`$, and $`R_2=1`$ if $`r_2<p_2`$. On the line $`p_1+p_2=1`$ (i.e. $`p_2=1p_1`$), $`R_2`$ takes the value 1 if $`1r_2>p_1`$. Choosing $`1r_2=r_1`$ (this choice does not affect the dynamics of a single replica) we have $`R_2=R_11`$ and Eq. (1) becomes, after some manipulations,

$$x^{}=R(\xi ^{(1)}\xi ^{(2)})\xi ^{(2)},$$

where $`R=R_1`$. One can easily check that

$$\xi ^{(1)}=y_{}yx,$$

and

$$\xi ^{(2)}=y_{}yx.$$

Finally, we obtain the evolution equation for the $`y`$

$$\begin{array}{cc}\hfill y^{}& =x^{}x_+^{}\hfill \\ & =R(y_{}yy_{}y)R_+(yy_+yy_+)y_{}yyy_+y\hfill \\ & =R(y_{}y)R_+(yy_+)y_{}yyy_+y,\hfill \end{array}$$
(B.3)
In this expression $``$ denotes the disjunction operation (OR). The sum modulo 2 (XOR) of all $`y_i`$ over the lattice is invariant with time, since all repeated terms cancel out ($`aa=0`$). Note that the kink dynamics uses correlated noise between neighboring sites.


format: text

text:
# Search for Resonances Decaying to 𝑒⁺-jet in 𝑒⁺𝑝 Interactions at HERA

## 1 Introduction

The $`e^+`$-jet mass spectrum in $`e^+p`$ scattering has been investigated with the ZEUS detector at HERA. An excess of events relative to Standard Model expectations has previously been reported by the H1  and ZEUS  collaborations in neutral current deep inelastic scattering at high $`x`$ and high $`Q^2`$. These events contain high-mass $`e^+`$-jet final states. Several models have been discussed  as possible sources of these events, including leptoquark production  and R-parity-violating squark production . This paper presents an analysis of ZEUS data specifically aimed at searching for high mass states decaying to $`e^+`$-jet.

Candidate events with high transverse energy, an identified final-state positron, and at least one jet are selected. The measured energies ($`E_e^{},E_j`$) and angles of the final-state positron and the jet with highest transverse momentum are used to calculate an invariant mass

$$M_{ej}^2=2E_e^{}E_j(1\mathrm{cos}\xi ),$$
(1)
where $`\xi `$ is the angle between the positron and jet. The angle between the outgoing and incoming positron in the $`e^+`$-jet rest frame, $`\theta ^{}`$, is also reconstructed using the measured energies and angles. No assumptions about the production process are made in the reconstruction of either $`M_{ej}`$ or $`\theta ^{}`$.

The search was performed using $`47.7`$ pb-1 of data collected in the 1994-1997 running periods. In the following, expectations from the Standard Model, leptoquark production and R-parity-violating squark production are summarized. After a discussion of the experimental conditions, the analysis is described and the $`M_{ej}`$ and $`\mathrm{cos}\theta ^{}`$ distributions presented. Since these distributions do not show a clear signal for a narrow resonance, limits on the cross section times branching ratio are extracted for the production of such a state. Limits are also presented in the mass versus coupling plane which can be applied to leptoquark and squark production.

## 2 Model Expectations

High-mass $`e^+`$-jet pairs, produced in the Standard Model (SM) via neutral current (NC) scattering, form the principal background to the search for heavy states. This process is reviewed first. Leptoquark (LQ) production and squark production in R-parity-violating ($`\mathit{}_P`$) supersymmetry are used as examples of physics beyond the SM that could generate the $`e^+`$-jet final state. The diagrams for NC and LQ processes are shown in Fig. 1. The squark production diagrams are similar to the LQ diagrams, but different decay modes are possible, as discussed below.

### 2.1 Standard Model Expectations

The kinematic variables used to describe the deep inelastic scattering (DIS) reaction

$$e^+pe^+X$$

are

$`Q^2`$ $`=`$ $`q^2=(kk^{})^2,`$ (2)
$`y`$ $`=`$ $`{\displaystyle \frac{qP}{kP}}\mathrm{and}`$ (3)
$`x`$ $`=`$ $`{\displaystyle \frac{Q^2}{2qP}},`$ (4)
where $`k`$ and $`k^{}`$ are the four-momenta of the incoming and outgoing positron, respectively, and $`P`$ is the four-momentum of the incoming proton. The center-of-mass energy is given by $`s=(k+P)^2(300\mathrm{GeV})^2`$. The NC interaction occurs between the positron and a parton (quark) inside the proton (see Fig. 1). The production of the large $`e^+`$-jet masses of interest requires high $`x`$ partons, where the valence quarks dominate the proton structure.

In leading-order electroweak theory, the cross section for the NC DIS reaction can be expressed as

$$\frac{d^2\sigma (e^+p)}{dxdy}=\frac{2\pi \alpha ^2}{sx^2y^2}\left[Y_+F_\mathit{2}Y_{}xF_\mathit{3}+y^2F_L\right]$$
(5)
with $`Y_\pm =1\pm (1y)^2`$ and $`\alpha `$ the fine structure constant. The contribution from the longitudinal structure function, $`F_L`$, is expected to be negligible in the kinematic range considered here.

The $`x`$ dependence of the NC cross section is very steep. In addition to the explicit $`1/x^2`$ factor, the structure functions $`F_2`$ and $`xF_3`$ are dominated at large $`x`$ by valence-quark densities that fall quickly for $`x>0.3`$. The $`y`$ dependence of the cross section is dominated by the $`1/y^2`$ term. The structure functions vary slowly with $`y`$ at fixed $`x`$. The uncertainty in the NC cross section predicted by Eq. 5 is dominated by the uncertainty in the structure functions (parton densities), and is small, about $`5`$% at the high-$`x`$ and moderate-$`y`$ ranges of this analysis . The quantity of interest in this paper is the $`e^+`$-jet cross section, which is sensitive to QCD corrections. The uncertainty arising from these corrections has been estimated to be small for this analysis.

For DIS or LQ events produced via the diagrams shown in Fig. 1 (i.e. assuming no QED or QCD radiation), the mass of the $`eq`$ system is related to $`x`$ via

$$M^2=sx$$
(6)
and $`\theta ^{}`$ is related to $`y`$ via

$$\mathrm{cos}\theta ^{}=12y.$$
(7)
The steeply-falling $`x`$ and $`y`$ dependences of DIS events will therefore produce distributions falling sharply with mass and peaking towards $`\mathrm{cos}\theta ^{}=1`$.

### 2.2 Leptoquark Production and Exchange

Leptoquark production is an example of new physics that could generate high-mass $`e^+`$-jet pairs. The set of leptoquarks with $`SU(3)\times SU(2)\times U(1)`$-invariant couplings has been specified . Only LQs with fermion number $`F=L+3B=0`$ are considered here, where $`L`$ and $`B`$ denote the lepton and baryon number, respectively. These leptoquarks are listed in Table 1 together with some of their properties. The $`F=0`$ LQs have higher cross sections in $`e^+p`$ scattering than $`e^{}p`$ scattering since in the $`e^+p`$ case a valence quark can fuse with the positron.

In principle, additional LQ types can be defined  which depend on the generations of the quarks and leptons to which they couple. Only LQs which preserve lepton flavor and which couple to first-generation quarks are considered in this analysis.

As shown in Fig. 1, leptoquark production can generate an s-channel resonance provided $`m_{LQ}<\sqrt{s}`$. Contributions to the $`e^+p`$ cross section would also result from u-channel exchange and interference of LQ diagrams with photon and $`Z^0`$ exchange. The cross section in the presence of a leptoquark can be written as

$$\frac{d^2\sigma (e^+p)}{dxdy}=\frac{d^2\sigma ^{NC}}{dxdy}+\frac{d^2\sigma _{u/NC}^{Int}}{dxdy}+\frac{d^2\sigma _{s/NC}^{Int}}{dxdy}+\frac{d^2\sigma _u^{LQ}}{dxdy}+\frac{d^2\sigma _s^{LQ}}{dxdy}.$$
(8)
The first term on the right-hand side of Eq. 8 represents the SM contribution discussed previously. The second (third) term arises from the interference between the SM and u-channel (s-channel) LQ diagram, and the fourth (fifth) term represents the u-channel (s-channel) LQ diagram alone. The additional contributions to the SM cross section depend on two parameters: $`m_{LQ}`$, and $`\lambda _R`$ or $`\lambda _L`$, the coupling to $`e_{L,R}^+`$ and quark. Leptoquarks of well-defined helicity ($`\lambda _R\lambda _L=0`$) are assumed for simplicity in the limit-setting procedure, and one species of LQ is assumed to dominate the cross section. The $`\mathrm{cos}\theta ^{}`$ dependence varies strongly for the different terms: it is flat for scalar-LQ production in the s-channel, and for vector-LQ exchange in the u-channel, while it varies as $`(1+\mathrm{cos}\theta ^{})^2`$ for vector-LQ production in the s-channel or scalar-LQ exchange in the u-channel. The interference terms produce a $`\mathrm{cos}\theta ^{}`$ dependence which is steeper due to the sharply-peaking $`\mathrm{cos}\theta ^{}`$ distribution in NC DIS.

In general, the s-channel term dominates the additional contributions to the SM cross section if $`m_{LQ}<\sqrt{s}`$, the coupling $`\lambda `$ is small, and the LQ is produced from a quark rather than an antiquark. However, there are conditions for which the other terms can become significant, or even dominant , leading to important consequences for the expected mass spectra and decay angular distributions. The u-channel and interference terms cannot produce a resonance peak in the mass spectrum and the angular distributions from such terms can behave more like those of NC deep inelastic scattering. Limits are presented in this paper for narrow-width LQ and under conditions for which the s-channel term dominates.

The width of a LQ depends on its spin and decay modes, and is proportional to $`m_{LQ}`$ times the square of the coupling. In the narrow-width approximation, the LQ production cross section is given by integrating the s-channel term :

$$\sigma ^{NWA}=(J+1)\frac{\pi }{4s}\lambda ^2q(x_0,\mu )$$
(9)
where $`J`$ represents the spin of the LQ, $`q(x_0,\mu )`$ is the quark density evaluated at $`x_0=m_{LQ}^2/s`$ and with the scale $`\mu =m_{LQ}^2`$. In the limit-setting procedure (Sect. 8), this cross section was corrected for expected QED and QCD (for scalar LQ only) radiative effects. The QCD corrections  enhance the cross section by 20 - 30% for the $`F=0`$ LQ considered here. The effect of QED radiation on the LQ production cross section was calculated and was found to decrease the cross section by 5-25% as $`m_{LQ}`$ increases from $`100290`$ GeV.

### 2.3 R-Parity-Violating Squark Production

In the supersymmetry (SUSY) superpotential, $`R`$-parity-violating terms of the form $`\lambda _{ijk}^{}L_L^iQ_L^j\overline{D}_R^k`$  are of particular interest for lepton-hadron collisions. Here, $`L_L`$, $`Q_L`$, and $`\overline{D}_R`$ denote left-handed lepton and quark doublets and the right-handed down-type quark-singlet chiral superfields, respectively. The indices $`i`$, $`j`$, and $`k`$ label their respective generations.

For $`i=1`$, which is the case for $`ep`$ collisions, these operators can lead to $`\stackrel{~}{u}`$\- and $`\stackrel{~}{d}`$-type squark production. There are 9 possible production couplings probed in $`e^+p`$ collisions, corresponding to the reactions

$`e^++\overline{u}_j`$ $``$ $`\stackrel{~}{\overline{d}}_k,`$ (10)
$`e^++d_k`$ $``$ $`\stackrel{~}{u}_j.`$ (11)
For production and decay via the $`\lambda _{1jk}^{}`$ coupling, squarks behave like scalar leptoquarks and the final state is indistinguishable, event by event, from Standard Model neutral and charged current events. However, as for the scalar leptoquarks, the angular distributions of the final-state lepton and quark will be different and this fact can be exploited in performing searches. Limits derived for scalar LQ production can then be directly related to limits on squark production and decay via $`\lambda _{ijk}^{}`$. In addition to the Yukawa couplings, gauge couplings also exist whereby $`\stackrel{~}{q}`$ can decay by radiating a neutralino or chargino which can subsequently decay. The final-state signature depends on the properties of the neutralino or chargino. The search for such decay topologies from a squark is outside the scope of this analysis.

## 3 Experimental Conditions

In 1994-97, HERA operated with protons of energy $`E_p=820`$ GeV and positrons of energy $`E_e=27.5`$ GeV. The ZEUS detector is described in detail in . The main components used in the present analysis were the central tracking detector (CTD) positioned in a 1.43 T solenoidal magnetic field and the uranium-scintillator sampling calorimeter (CAL). The CTD was used to establish an interaction vertex with a typical resolution along (transverse to) the beam direction of $`0.4(0.1)`$ cm. It was also used in the positron-finding algorithm that associated a charged track with an energy deposit in the calorimeter. The CAL was used to measure the positron and hadronic energies. The CAL consists of a forward part (FCAL), a barrel part (BCAL) and a rear part (RCAL), with depths of $`7,\mathrm{\hspace{0.33em}5}\mathrm{and}\mathrm{\hspace{0.33em}4}`$ interaction lengths, respectively. The FCAL and BCAL are segmented longitudinally into an electromagnetic section (EMC), and two hadronic sections (HAC1,2). The RCAL has one EMC and one HAC section. The cell structure is formed by scintillator tiles; cell sizes range from $`5\times 20`$ cm2 (FEMC) to $`24.4\times 35.2`$ cm2 at the front face of a BCAL HAC2 cell. The light generated in the scintillator is collected on both sides of the module by wavelength-shifter bars, allowing a coordinate measurement based on knowledge of the attenuation length in the scintillator. The light is converted into an electronic signal by photomultiplier tubes. The cells are arranged into towers consisting of $`4`$ EMC cells, a HAC1 cell and a HAC2 cell (in FCAL and BCAL). The transverse dimensions of the towers in FCAL are $`20\times 20`$ cm2. One tower is absent at the center of the FCAL and RCAL to allow space for passage of the beams. The outer boundary of the inner ring of FCAL towers, used to define a fiducial cut for the jet reconstruction, defines a box of $`60\times 60`$ cm2.

Under test beam conditions, the CAL has energy resolutions of $`\sigma /E=0.18/\sqrt{E}`$ for positrons hitting the center of a calorimeter cell and $`\sigma /E=0.35/\sqrt{E}`$ for single hadrons, where energies are in GeV. In the ZEUS detector, the energy measurement is affected by the energy loss in the material between the interaction point and the calorimeter. For the events selected in this analysis, the positrons predominantly strike the BCAL, while the jets hit the FCAL. The in-situ positron-energy resolution in the BCAL has been determined to average $`\sigma /E=0.32/\sqrt{E}0.03`$ while the jet-energy resolution in the FCAL averages $`\sigma /E=0.55/\sqrt{E}0.02`$. The jet-energy resolution was determined by comparing reconstructed jet energies in the calorimeter with the total energy of the particles in the hadronic final-state using Monte Carlo simulation, and therefore includes small contributions from the jet-finding algorithm.

In the reconstruction of the positron and jet energies, corrections were applied for inactive materials located in front of the calorimeter and for non-uniformities in the calorimeter response . For the high energies important in this analysis, the overall energy scale is known to $`1`$% for positrons in BCAL and $`2`$% for hadrons in FCAL and BCAL. The electromagnetic energy scale was determined by a comparison with momentum measurements in the central tracking detector (using lower-energy electrons and positrons). Its linearity was checked with energies reconstructed from the double angle (DA) method . The hadronic-energy scales in the FCAL and BCAL were determined by using transverse-momentum balance in NC DIS events.

The angular reconstruction was performed using a combination of tracking and calorimeter information. From Monte Carlo studies, the polar-angle resolutions were found to be $`2.5`$ mrad for positrons and approximately ($`220/\sqrt{E}4`$) mrad for jets with energies above $`100`$ GeV.

The luminosity was measured from the rate of the bremsstrahlung process $`e^+pe^+p\gamma `$ , and has an uncertainty of $`1.6`$%.

The ZEUS coordinate system is right-handed and centered on the nominal interaction point, with the $`Z`$ axis pointing in the direction of the proton beam (forward) and the $`X`$ axis pointing horizontally toward the center of HERA. The polar angle $`\theta `$ is defined with respect to the $`Z`$ axis.

## 4 Event Selection

The events of interest with large $`e^+`$-jet mass contain a final-state positron at a large angle and of much higher energy than that of the incident positron beam, as well as one or more energetic jets. The only important SM source of such events is NC scattering with large $`Q^2`$. Other potential backgrounds, such as high transverse-energy ($`E_T`$) photoproduction, were determined to be negligible.

The following requirements selected events of the desired topology:

* A reconstructed event vertex was required in the range $`|Z|<50`$ cm.
* The total transverse energy, $`E_T`$, was required to be at least $`60`$ GeV.
* An identified  positron was required with energy $`E_e^{}>25`$ GeV, located either in the FCAL or BCAL. The positron was required to be well-contained in the BCAL or FCAL and not to point to the BCAL/FCAL interface, at approximately $`31^{}<\theta <36^{}`$. Positrons within $`1.5`$ cm of the boundary between adjacent BCAL modules, as determined by tracking information, were also discarded to remove showers developing in the wavelength-shifter bars.
* A hadronic jet with transverse momentum $`P_T^j>15`$ GeV, located in a region of good containment, was required. The jets were reconstructed using the longitudinally-invariant $`k_T`$-clustering algorithm  in the inclusive mode . Only jets with a reconstructed centroid outside the inner ring of FCAL towers were considered. In events where multiple jets were reconstructed, the jet with highest transverse momentum was used. After all cuts, 12% of the events had more than one jet, both in the data and Monte Carlo simulation (see below).

The $`E_T`$ cut, the jet-containment cut and the positron-containment cut define the available kinematic region for further analysis, as shown in Fig. 2. The jet containment cut, in particular, limits the values of $`\mathrm{cos}\theta ^{}`$ that can be measured at the highest $`e^+`$-jet masses. Because most such events have $`\mathrm{cos}\theta ^{}`$ near $`1`$, the acceptance for NC DIS events (with $`E_T>60`$ GeV) falls below $`10`$% for masses beyond 220 GeV. In the region allowed by the cuts shown in Fig. 2, the acceptance is typically $`80`$%.

A total of $`7103`$ events remained after applying all cuts, compared to $`6949\pm 445`$ events predicted by the NC Monte Carlo simulation based on the measured luminosity of $`47.7`$ pb-1 (the sources of uncertainty on the expected number of events are described in Sect. 7.1). The $`E_T`$ distributions for data and NC simulation are compared in Fig. 3a. The positron transverse-momentum ($`P_T^e`$) spectrum, jet transverse-momentum ($`P_T^j`$) spectrum and the ratio $`P_T^e/P_T^h`$, where $`P_T^h`$ is the transverse momentum of the hadronic system, are shown in Figs. 3(b-d), respectively. The missing transverse momentum, $`\mathit{}_T`$, and the longitudinal momentum variable, $`EP_Z`$ , are compared in Figs. 3(e,f). The global properties of the events are well reproduced by the simulation.

## 5 Event Simulation

The SM deep inelastic scattering events were simulated using the HERACLES 4.5.2  program with the DJANGO 6 version 2.4  interface to the hadronization programs. In HERACLES, corrections for initial- and final-state electroweak radiation, vertex and propagator corrections, and two-boson exchange are included. The NC DIS hadronic final state was simulated using the MEPS model of LEPTO 6.5 , which includes order $`\alpha _S`$ matrix elements and models of higher-order QCD radiation. As a systematic check, the NC final state was simulated using the color-dipole model of ARIADNE 4.08 . The CTEQ4 parton-distribution set  was used to evaluate the expected number of events from NC DIS scattering.

The leptoquark events were generated using PYTHIA 6.1 . This program takes into account the finite width of the LQ, but only includes the s-channel diagrams. Initial- and final-state QCD radiation from the quark and the effect of LQ hadronization before decay are taken into account, as are initial- and final-state QED radiation from the positron.

The generated events were input into a GEANT-based  program which simulated the response of the ZEUS detector. The trigger and offline processing requirements applied to the data were applied to the simulated events. The luminosity of the NC Monte Carlo samples ranges from $`46`$ pb-1 at $`Q^2=400`$ GeV2 to $`7.310^6`$ pb-1 at $`Q^2=50000`$ GeV2.

## 6 Mass and $`\theta ^{}`$ Reconstruction

The mass of each $`e^+`$-jet pair was reconstructed from the measured energies and angles of the positron and jet as described by Eq. 1. This formula makes no correction for the finite jet mass. Possible mass shifts and the resolutions for resonant lepton-hadron states were estimated from PYTHIA. Narrow scalar LQ events in the mass range $`150290`$ GeV were simulated. The mean mass for reconstructed events was found to be within 6% of the generated value, while the peak position as determined by a Gaussian fit was typically lower than the generated value by only 1%. The average mass resolution, determined from a Gaussian fit to the peak of the reconstructed mass spectrum, ranged from 5.5% to 3% for masses from 150 to 290 GeV. The RMS of the distribution was typically twice as large.

The positron scattering angle in the $`e^+`$-jet rest frame, $`\theta ^{}`$, was reconstructed as the angle between the incoming and outgoing positron directions in this frame. These directions were determined by performing a Lorentz transformation using the measured positron and jet energies and angles in the laboratory frame. The resolution in $`\mathrm{cos}\theta ^{}`$ near $`|\mathrm{cos}\theta ^{}|=1`$, as determined from a Gaussian fit, was $`0.01`$ degrading to $`0.03`$ as $`|\mathrm{cos}\theta ^{}|`$ decreases. The shift in $`\mathrm{cos}\theta ^{}`$ was less than $`0.01`$ for both the NC MC and the leptoquark MC.

In order to determine limits on leptoquark and squark production, the mass of the electron-hadron system was reconstructed by the constrained-mass method. This method reconstructs the $`e^+`$-hadron mass as

$`M_{CJ}`$ $`=`$ $`\sqrt{2E_e(E+P_Z)}`$ (12)
where $`(E+P_Z)`$ is the sum of the energy and $`P_Z`$ contributions from the positron and all jets satisfying $`P_T^j>15`$ GeV and pseudorapidity $`\eta _j<3`$ (with the highest $`P_T`$ jet required to be outside the FCAL inner ring). The $`\eta _j`$ cut removes contributions from the proton remnant. The constraints $`\mathit{}_T=0`$ and $`EP_Z=2E_e`$, which are satisfied by fully contained events, have been assumed in arriving at this equation. When using this mass-reconstruction method, events with measured $`EP_Z<40`$ GeV were removed to avoid large initial-state QED radiation.

The $`M_{CJ}`$ method gave, on average, improved resolution over the $`M_{ej}`$ method for narrow LQ MC events. The improved resolution occurred at smaller $`\mathrm{cos}\theta ^{}`$ (for $`\mathrm{cos}\theta ^{}0`$ the mass resolution determined from a Gaussian fit to the reconstructed mass distribution for $`m_{LQ}=200`$ GeV was about $`1.5`$% in the $`M_{CJ}`$ method and $`3`$% for the $`M_{ej}`$ method); at the larger $`\mathrm{cos}\theta ^{}`$ values where NC DIS events are concentrated, the resolutions of the two methods were similar (about $`3`$% for $`m_{LQ}=200`$ GeV). The $`M_{CJ}`$ method relies on constraints which do not necessarily apply to a resonant state whose properties cannot be anticipated in detail. We therefore choose to use the $`M_{ej}`$ method as our primary search method. The $`M_{CJ}`$ method is used in the limit-setting procedure.

## 7 $`M_{ej}`$ and $`\mathrm{cos}\theta ^{}`$ Distributions

The reconstructed values of $`M_{ej}`$ are plotted versus $`\mathrm{cos}\theta ^{}`$ for the selected events in Fig. 4. Most of the events are concentrated at large $`\mathrm{cos}\theta ^{}`$ and small mass, as expected from Standard Model NC scattering. The five events indicated as open circles are from data taken in 1994-96, with total luminosity $`20`$ pb-1. They were the subject of a previous publication . In this earlier analysis, the kinematic variables were reconstructed with the DA method. The five events also stand out with the $`M_{ej}`$ reconstruction technique. The average value of $`M_{ej}`$ for these events is $`224`$ GeV, or $`7`$ GeV less than the corresponding mass calculated previously via $`M=\sqrt{sx_{DA}}`$, where $`x_{DA}`$ is the estimator of Bjorken-$`x`$ calculated with the DA method. This mass shift is compatible with expectations based on resolution and initial state radiation effects. With the present luminosity of $`47.7`$ pb-1, 7 events are observed in the region of $`M_{ej}>200`$ GeV and $`\mathrm{cos}\theta ^{}<0.5`$, where $`5.0`$ events are expected.

The $`M_{ej}`$ spectrum for events with $`M_{ej}>100`$ GeV is shown in Fig. 5a on a logarithmic scale. The high-mass part of the spectrum is shown on a linear scale in the inset. The predicted number of events ($`N^{pred}`$) from NC processes is shown as the histogram. The ratio of the measured mass spectrum to the expectation is shown in Fig. 5b. The shaded band indicates the systematic uncertainty on the expectations.

### 7.1 Systematic Uncertainties

The uncertainty on $`N^{pred}`$ varies with mass from $`7`$% at 100 GeV up to 30% at 250 GeV. The most important uncertainties are on the energy scale and the jet position. The NC DIS cross section given in Eq. 5 (neglecting $`F_L`$) can be rewritten in terms of the $`e^+q`$ invariant mass, $`M`$, and the polar angle of the outgoing struck quark in the laboratory frame, $`\gamma `$:

$$\frac{d^2\sigma (e^+p)}{dMd\gamma }=\frac{32\pi \alpha ^2E_e^2\mathrm{sin}\gamma }{M^5(1\mathrm{cos}\gamma )^2}\left[Y_+F_\mathit{2}Y_{}xF_\mathit{3}\right].$$
(13)
The mass dependence is very steep. In addition to the explicit $`M^5`$ dependence, there is also a strong suppression of high masses implicit in the structure functions. An incorrect energy scale will produce a shift in the mass spectrum and potentially a significant error in the number of expected events at a given mass. The dependence on the quark angle is also steep, approximately $`\gamma ^3`$ at small $`\gamma `$. The number of events passing the jet fiducial cut is therefore strongly dependent on the accuracy of the jet position reconstruction. The jet fiducial-volume cut requires the highest-$`P_T`$ jet to point outside the inner ring of FCAL towers. Many distributions from data and MC were compared to search for possible systematic biases.

The dominant sources of uncertainty are itemized below in order of decreasing importance:

1. Knowledge of the calorimeter energy scales:

The scale uncertainties discussed in Sect. 3 are $`1`$% for BCAL positrons and $`2`$% for hadrons, leading to an uncertainty of $`5(18)`$% in $`N^{pred}`$ at $`M_{ej}=100(210)`$ GeV.
2. Uncertainties in the simulation of the hadronic energy flow, including simulation of the proton remnant, the energy flow between the struck quark and proton remnant, and possible detector effects in the innermost calorimeter towers:

Many distributions of data and MC were compared and no important systematic differences were found. Figure 6 shows the fraction of the jet energy in the inner ring of FCAL towers associated with the highest $`P_T`$ jet as a function of $`\eta _j`$. This is shown for all events in Fig. 6a, as well as for those with $`M_{ej}>210`$ GeV in Fig. 6b. For the highest $`\eta _j`$ values considered, this ratio is about $`20`$%. The energy located in the innermost towers of the FCAL and not associated with the highest $`P_T`$ jet is shown in Fig. 6c,d, and compared to the MC simulation. No large differences are seen between data and MC (the lowest $`\eta `$ bin in Fig 6d contains only five data events). The innermost towers of the FCAL have a larger uncertainty in the energy scale than the rest of the FCAL owing to their slightly different construction and proximity to the beam. The energy in these cells has been varied by $`\pm 10`$%. As a test of the simulation of the forward energy flow, the ARIADNE MC has been used instead of the LEPTO MC. These tests yielded variations in $`N^{pred}`$ of 13% at $`M_{ej}=210`$ GeV.
3. Uncertainty in the parton density functions:

The parton density functions were estimated as in , and led to an uncertainty of 5% in $`N^{pred}`$ at $`M_{ej}=210`$ GeV.
4. Uncertainties in the acceptance:

The alignment of the FCAL was determined to better than $`5`$ mm, and various jet position reconstruction algorithms were compared. These studies yielded an uncertainty of 2% in $`N^{pred}`$.
5. Uncertainties in the energy resolution functions.

These were studied by comparing tracking information with calorimeter information for individual events, as well as by comparing different reconstruction methods. The MC energies were smeared by additional amounts to represent these uncertainties, leading to a variation of less than $`5`$% in $`N^{pred}`$.

Other uncertainties include positron finding efficiency, luminosity determination, vertex simulation, multijet production rates, and hadronization simulation. These were found to be small in comparison to the items listed above. The overall systematic uncertainty was obtained by summing the contributions from all these sources in quadrature.

### 7.2 Discussion

The data in Fig. 5 are in good agreement with the SM expectations up to $`M_{ej}210`$ GeV. Some excess is seen at higher masses. For $`M_{ej}>210`$ GeV, $`49`$ events were observed in the data, while $`24.7\pm 5.6`$ events are expected. A careful study of individual events in this mass region uncovered no signs of reconstruction errors. Rather, the events always contain clear examples of a high-energy positron (typically $`70`$ GeV) near $`90^{}`$ and a high-energy jet (typically $`400`$ GeV) in the forward direction (2 events have a second jet, in accord with NC DIS Monte Carlo expectations). The distributions shown in Fig. 3 for all the data are restricted to the events with $`M_{ej}>210`$ GeV in Fig. 7. Whereas the shapes of the distributions are similar, the data lie systematically above the MC, which is normalized to the integrated luminosity.

The events with large $`M_{ej}`$ have characteristics similar on average to NC DIS events. In particular, the $`\mathrm{cos}\theta ^{}`$ projection of the events with $`M_{ej}>210`$ GeV is shown in Fig. 8 and compared to the MC expectations for neutral current DIS (solid histogram). The expectations for narrow s-channel scalar and vector LQ production are also shown for comparison. For $`F=0`$ LQs with $`\lambda <1`$, the u-channel and interference terms would not significantly affect these expectations. The shape of the data and NC MC $`\mathrm{cos}\theta ^{}`$ distributions are qualitatively similar, peaking at high values of $`\mathrm{cos}\theta ^{}`$.

In summary, there is some excess of events with $`M_{ej}>210`$ GeV above the Standard Model predictions. The probability of observing such an excess depends strongly on possible systematic biases. The most important of these are biases in the energy scales. As a test, many MC experiments were generated where the jet energy scale was shifted by +2% and the electron energy scale by +1%. A window of width $`3\sigma (M_{ej})`$, where $`\sigma (M_{ej})`$ is the mass resolution at mass $`M_{ej}`$, was moved over the accessible mass range. For each simulated experiment, the number of observed events within the mass window was compared with the nominal expectations as a function of $`M_{ej}`$, seeking the excess which gave the largest statistical significance. The same procedure was applied to the data. As a result, it was found that 5% of the simulated experiments would observe, somewhere in the mass spectrum, an excess of statistical significance at least as large as the one found in the data. The excess is therefore not statistically compelling. Furthermore, the events have the characteristics of neutral current scattering. Limits are therefore set on the production of narrow scalar or vector states, as discussed below.

## 8 Limits on Narrow Scalar and Vector States

Limits are set on the production cross section times branching ratio into positron+jet, $`\sigma B`$, for a narrow scalar or vector state. For definiteness, limits on coupling strength versus mass for $`F=0`$ leptoquarks are presented, as well as limits on $`\lambda \sqrt{B}`$ versus mass for scalar states coupling to $`u`$ or $`d`$ quarks, such as $`\mathit{}_P`$ squarks. The limits are extracted for $`\lambda 1`$, allowing the use of the narrow-width approximation assumed in Eq. 9.

The $`M_{CJ}`$ mass reconstruction method was used to set limits as described in Sect. 6. The positron fiducial cuts were removed since this method is less sensitive to the positron-energy measurement, while the cut $`EP_Z>40`$ GeV was applied to reduce radiative effects. The mass spectrum reconstructed with this technique is shown in Fig. 9a. In total, 8026 events passed all selection cuts while 7863 events are predicted by the MC.

The leptoquark MC described in Sect. 5 was used to determine the event selection efficiency and the acceptance of the fiducial cuts, as well as to estimate the mass resolution. This MC and the NC background simulation were used to calculate an optimal bin width, $`\mathrm{\Delta }M_{CJ}`$, for each $`M_{CJ}`$, and optimal $`\mathrm{cos}\theta ^{}`$ range, $`\mathrm{cos}\theta ^{}<\mathrm{cos}\theta _{cut}^{}`$, to obtain on average the best limits on LQ couplings. The bin widths were typically 20 GeV. The values of $`\mathrm{cos}\theta _{cut}^{}`$ for setting limits range from $`0.5`$ to $`0.9`$ for vector leptoquarks with masses between $`150290`$ GeV, and from $`0.1`$ to $`0.9`$ for scalar leptoquarks in the same mass range. The mass spectrum after applying the optimal $`\mathrm{cos}\theta ^{}`$ cut for the scalar search is shown in Fig. 9b. No significant deviations from expectations are seen after applying this cut.

The 95% confidence level (CL) limits on $`\sigma B`$ were obtained directly from the observed number of data events with $`\mathrm{cos}\theta ^{}<\mathrm{cos}\theta _{cut}^{}`$ in the particular mass window . The procedure described in  was extended to include the systematic uncertainties in the numbers of predicted events. This was found to have negligible effect on the limits. The limits for a narrow scalar or vector state are shown in Fig. 10. These limits lie between $`1`$ and $`0.1`$ pb as the mass increases from $`150`$ to $`290`$ GeV.

The 95% CL exclusion limits for different species of LQ are given in the coupling versus mass plane in Fig. 11. The limits exclude leptoquarks with coupling strength $`\lambda =\sqrt{4\pi \alpha }0.3`$ for masses up to 280 GeV for specific types of F=0 leptoquarks. The H1 collaboration has recently published similar limits . In Fig. 11, the ZEUS results are compared to recent limits from OPAL. At LEP , sensitivity to a high-mass LQ arises from effects of virtual LQ exchange on the hadronic cross section. The HERA and LEP limits are complementary to Tevatron limits , which are independent of the coupling $`\lambda _{L,R}`$. The limits by D0 (CDF) extend up to $`225(213)`$ GeV for a scalar LQ with 100% branching ratio to $`eq`$. The D0 limits are shown as vertical lines in Fig. 11. The Tevatron limits for vector LQs are model dependent , but are expected to be considerably higher than for scalar LQs.

The ZEUS limits presented in Fig. 11 can also be applied to any narrow state which couples to a positron and a $`u`$ or $`d`$ quark and with unknown branching ratio to $`e^+`$-jet(s). These states correspond to the leptoquark types as labelled in the figure. For these states, the limits are on the quantity $`\lambda \sqrt{B}`$. Examples of scalar states for which these limits apply are $`\mathit{}_P`$-squarks (e.g.,the limit on the $`\stackrel{~}{S}_{1/2}^L(e^+d)`$ LQ can be read as a limit on the $`\lambda _{1j1}^{}`$ R-parity-violating coupling).

## 9 Conclusion

Data from $`47.7`$ pb-1 of $`e^+p`$ collisions at a center-of-mass energy of 300 GeV have been used to search for a resonance decaying into $`e^+`$-jet. The invariant mass of the $`e^+`$-jet pair was calculated directly from the measured energies and angles of the positron and jet. This approach makes no assumptions about the production mechanism of such a state.

The observed mass spectrum is in good agreement with Standard Model expectations up to $`e^+`$-jet masses of about $`210`$ GeV. Above this mass, some excess is seen. The angular distribution of these events is typical of high-$`Q^2`$ neutral current events and does not give convincing evidence for the presence of a narrow scalar or vector state. By applying restrictions on the decay angle to optimize sensitivity to a narrow state in the presence of NC background, limits have been derived on the cross section times decay branching fraction for a scalar or vector state decaying into positron and jet(s). These limits can be interpreted, for example, as limits on leptoquark or R-parity-violating squark production. Limits on the production of leptoquarks and squarks are presented in the coupling strength versus mass plane. At a coupling strength $`\lambda =0.3`$, new states are ruled out at 95% confidence level for masses between $`150`$ and $`280`$ GeV.

## Acknowledgements

We thank the DESY Directorate for their strong support and encouragement, and the HERA machine group for their diligent efforts. We are grateful for the support of the DESY computing and network services. The design, construction and installation of the ZEUS detector have been made possible by the ingenuity and effort of many people from DESY and home institutes who are not listed as authors. It is also a pleasure to thank W. Buchmüller, R. Rückl and M. Spira for useful discussions.


text:
# Untitled Document

UTTG-01-00

A Priori Probability Distribution of the Cosmological Constant

Steven Weinberg***Electronic address: weinberg@physics.utexas.edu

Theory Group, Department of Physics, University of Texas
Austin, TX, 78712

Abstract

In calculations of the probability distribution for the cosmological constant, it has been previously assumed that the a priori probability distribution is essentially constant in the very narrow range that is anthropically allowed. This assumption has recently been challenged. Here we identify large classes of theories in which this assumption is justified.

I. INTRODUCTION

In some theories of inflation1 and of quantum cosmology2 the observed big bang is just one of an ensemble of expanding regions in which the cosmological constant takes various different values. In such theories there is a probability distribution for the cosmological constant: the probability $`d𝒫(\rho _V)`$ that a scientific society in any of the expanding regions will observe a vacuum energy between $`\rho _V`$ and $`\rho _V+\rho _V`$ is given by3,4,5

$$d𝒫(\rho _V)=𝒫_{}(\rho _V)𝒩(\rho _V)d4\rho _V,$$
(1)
where $`𝒫_{}(\rho _V)d\rho _V`$ is the a priori probability that an expanding region will have a vacuum energy between $`\rho _V`$ and $`\rho _V+d\rho _V`$ (to be precise, weighted with the number of baryons in such regions), and $`𝒩(\rho _V)`$ is proportional to the fraction of baryons that wind up in galaxies. (The constant of proportionality in $`𝒩(\rho _V)`$ is independent of $`\rho _V`$, because once a galaxy is formed the subsequent evolution of its stars, planets, and life is essentially unaffected by the vacuum energy.)

The factor $`𝒩(\rho _V)`$ vanishes except for values of $`\rho _V`$ that are very small by the standards of elementary particle physics, because for $`\rho _V`$ large and positive there is a repulsive force that prevents the formation of galaxies6 and hence of stars, while for $`\rho _V`$ large and negative the universe recollapses too fast for galaxies or stars to form.7 The fraction of baryons that form galaxies has been calculated5 for $`\rho _V>0`$ under reasonable astrophysical assumptions. On the other hand, we know little about the a priori probability distribution $`𝒫_{}(\rho _V)`$. However, the range of values of $`\rho _V`$ in which $`𝒩(\rho _V)0`$ is so narrow compared with the scales of energy density typical of particle physics that it had seemed reasonable in earlier work 4,5 to assume that $`𝒫_{}(\rho _V)`$ is constant within this range, so that $`d𝒫(\rho _V)`$ can be calculated as proportional to $`𝒩(\rho _V)d\rho _V`$. In an interesting recent article,8 Garriga and Vilenkin have argued that this assumption (which they call “Weinberg’s conjecture”) is generally not valid. This raises the problem of characterizing those theories in which this assumption is valid and those in which it is not.

It is shown in Section II that this assumption is in fact valid for a broad range of theories, in which the different regions are characterized by different values of a scalar field that couples only to itself and gravitation. The deciding factor is how we impose the flatness conditions on the scalar field potential that are needed to ensure that the vacuum energy is now nearly time-independent. If the potential is flat because the scalar field renormalization constant is very large, then the a priori probability distribution of the vacuum energy is essentially constant within the anthropically allowed range, for scalar potentials of generic form. It is also essentially constant for a large class of other potentials. Section III is a digression, showing that the same flatness conditions ensure tht the vacuum energy has been roughly constant since the end of inflation. Section IV takes up the sharp peaks in the a priori probability found in theories of quantum cosmology and eternal inflation.

II. SLOWLY ROLLING SCALAR FIELD

One of the possibilities considered by Garriga and Vilenkin is a vacuum energy that depends on a homogeneous scalar field $`\varphi (t)`$ whose present value is governed by some smooth probability distribution. The vacuum energy is

$$\rho _V=V(\varphi )+\frac{1}{2}\dot{\varphi }^2,$$
(2)
and the scalar field time-dependence is given by

$$\ddot{\varphi }+3H\dot{\varphi }=V^{}(\varphi ),$$
(3)
where $`H(t)`$ is the Hubble fractional expansion rate, $`V(\varphi )`$ is the scalar field potential, dots denote derivatives with respect to time, and primes denote derivatives with respect to $`\varphi `$. Following Garriga and Vilenkin,8 we assume that at present the scalar field energy appears like a cosmological constant because the field $`\varphi `$ is now nearly constant in time, and that this scalar field energy now dominates the cosmic energy density. For this to make sense it is necessary for the potential $`V(\varphi )`$ to satisfy certain flatness conditions. In the usual treatment of a slowly rolling scalar, one neglects the inertial term $`\ddot{\varphi }`$ in Eq. (3) as well as the kinetic energy term $`\dot{\varphi }^2/2`$ in Eq. (2). With the inertial term neglected, the condition that $`V(\varphi )`$ should change little in a Hubble time $`1/H`$ is that9

$$V^2(\varphi )3H^2|V(\varphi )|.$$
(4)
With the scalar field energy dominating the total cosmic energy density, the Friedmann equation gives

$$|V(\varphi )|\rho _V3H^2/8\pi G,$$
(5)
so Eq. (4) requires

$$\left|V^{}(\varphi )\right|\sqrt{8\pi G}\rho _V.$$
(6)
(The kinetic energy term $`\dot{\varphi }^2/2`$ in Eq. (2) can be neglected under the slightly weaker condition

$$\left|V^{}(\varphi )\right|\sqrt{18H^2\left|V(\varphi )\right|}\sqrt{48\pi G}\rho _V,$$

which is the flatness condition given by Garriga and Vilenkin.) There is also a bound on the second derivative of the potential, needed in order for the inertial term to be neglected. With the scalar field energy dominating the total cosmic energy density, this condition requires that9

$$\left|V^{\prime \prime }(\varphi )\right|8\pi G\rho _V.$$
(7)
As Garriga and Vilenkin correctly pointed out, the smallness of the slope of $`V(\varphi )`$ means that $`\varphi `$ may vary appreciably even when $`\rho _VV(\varphi )`$ is restricted to the very narrow anthropically allowed range of values in which galaxy formation is non-negligible. They concluded that it would be possible for the a priori probability $`𝒫_{}(\rho _V)`$ to vary appreciably in this range. In particular, Garriga and Vilenkin assumed an a priori probability distribution for $`\varphi `$ that is constant in the anthropically allowed range, in which case the a priori probability distribution for $`\rho _V`$ is

$$𝒫_{}(\rho _V)1/|V^{}(\varphi )|$$
(8)
which they said could vary appreciably in the anthropically allowed range.

Though possible, this rapid variation is by no means the generic case. As already mentioned, the second as well as the first derivative of the potential must be small, so that the a priori probability density (8) may change little in the anthropically allowed range. It all depends on how the flatness conditions are satisfied. There are two obvious ways that one might try to make the potential sufficiently flat. Potentials of the first type are of the general form

$$V(\varphi )=V_1f(\lambda \varphi ),$$
(9)
where $`V_1`$ is some large energy density, in the range of $`m_W^4`$ to $`G^2`$; the constant $`\lambda `$ is very small: and $`f(x)`$ is some dimensionless function involving no very large or very small parameters. Potentials of the second type are of the general form

$$V(\varphi )=V_1\left[1ϵg(\lambda \varphi )\right],$$
(10)
where $`V_1`$ is again some large energy density; $`\lambda `$ is here some fixed inverse mass, perhaps of order $`\sqrt{G}`$; now it is $`ϵ`$ instead of $`\lambda `$ that is very small; and $`g(x)`$ is some other dimensionless function involving no very large or very small parameters.

For potentials (9) of the first type, it is always possible to meet all observational conditions by taking $`\lambda `$ sufficiently small, provided that the function $`f(x)`$ has a simple zero at a point $`x=a`$ of order unity, with derivatives at $`a`$ of order unity. Because $`V_1`$ is so large, the present value of $`\lambda \varphi `$ must be very close to the assumed zero $`a`$ of $`f(x)`$. With $`f^{}(a)`$ and $`f^{\prime \prime }(a)`$ of order unity, the flatness conditions (6) and (7) are both satisfied if

$$|\lambda |\left(\frac{\rho _V}{V_1}\right)\sqrt{8\pi G}.$$
(11)
Galaxy formation is only possible for $`|V(\varphi )|`$ less than an upper bound $`V_{\mathrm{max}}`$ of the order of the mass density of the universe at the earliest time of galaxy formation,6 which in the absence of fine tuning of the cosmological constant is very much less than $`V_1`$. The anthropically allowed range of $`\varphi `$ is therefore given by

$$\mathrm{\Delta }\varphi |\varphi a/\lambda |_{\mathrm{max}}=\frac{V_{\mathrm{max}}}{|\lambda f^{}(a)V_1|}.$$
(12)
The fractional change in the a priori probability density $`1/|V^{}(\varphi )|`$ in this range is then

$$\left|\frac{V^{\prime \prime }(\varphi )\mathrm{\Delta }\varphi }{V^{}(\varphi )}\right|=\left|\frac{V_{\mathrm{max}}}{V_1}\right|\left|\frac{f^{\prime \prime }(a)}{f^2(a)}\right|,$$
(13)
with no dependence on $`\lambda `$. As long as the factor $`f^{\prime \prime }(a)/f^2(a)`$ is roughly of order unity the fractional variation (13) in the a priori probability will be very small, as was assumed in references 4 and 5.

This reasoning applies to potentials of the form

$$V(\varphi )=V_1\left[1(\lambda \varphi )^n\right],$$

which, as already noted by Garriga and Vilenkin, lead to an priori probability distribution that is nearly constant in the anthropically allowed range. (In this case $`a=1`$ and $`f^{\prime \prime }(a)/f^2(a)=(1n)/n`$.) But this reasoning also applies to the “washboard potential” that was taken as a counterexample by Garriga and Vilenkin, which with no loss of generality can be put in the form:

$$V(\varphi )=V_1\left[1+\alpha \lambda \varphi +\beta \mathrm{sin}(\lambda \varphi )\right].$$

The zero point $`a`$ is here determined by the condition

$$1+\alpha a+\beta \mathrm{sin}a=0,$$

and the factor $`f^{\prime \prime }(a)/f^2(a)`$ in Eq. (13) is

$$\frac{f^{\prime \prime }(a)}{f^2(a)}=\frac{\beta \mathrm{sin}a}{(\alpha +\beta \mathrm{cos}a)^2}.$$

If the flatness condition is satisfied by taking $`\lambda `$ small, with $`\alpha `$ and $`\beta `$ of order unity, as is assumed for potentials of the first kind, then the factor $`f^{\prime \prime }(a)/f^2(a)`$ in Eq. (13) is of order unity unless $`\alpha `$ and $`\beta `$ happen to be chosen so that

$$\left|1+\alpha \mathrm{cos}^1\left(\frac{\alpha }{\beta }\right)+\beta \sqrt{1\frac{\alpha ^2}{\beta ^2}}\right|1.$$

Of course it would be possible to impose this condition on $`\alpha `$ and $`\beta `$, but this is the kind of fine-tuning that would be upset by adding a constant of order $`V_1`$ to the potential. Aside from this exception, for all $`\alpha `$ and $`\beta `$ of order unity the factor $`f^{\prime \prime }(a)/f^2(a)`$ is of order unity, so the washboard potential also yields an a priori probability distribution for the vacuum energy that is flat in the anthropically allowed range.

In contrast, for potentials (10) of the second kind the flatness conditions are not necessarily satisfied no matter how small we take $`ϵ`$. Because the present vacuum energy is much less than $`V_1`$, the present value of $`\varphi `$ must be very close to a value $`\varphi _ϵ`$, satisfying

$$g(\lambda \varphi _ϵ)=1/ϵ.$$
(14)
This requires $`\lambda \varphi _ϵ`$ to be near a singularity of the function $`g(x)`$, perhaps at infinity, so it is not clear in general that such a potential would have small derivatives at $`\lambda \varphi _ϵ`$ for any value of $`ϵ`$. For instance, for an exponential $`g(x)=\mathrm{exp}(x)`$ we have $`\varphi _ϵ=\mathrm{ln}ϵ/\lambda `$, and $`V^{}(\varphi _ϵ)`$ approaches an $`ϵ`$-independent value proportional to $`\lambda `$, which is not small unless we take $`\lambda `$ very small, in which case have a potential of the first kind, for which as we have seen the a priori probability density (8) is flat in the anthropically allowed range. The flatness conditions are satisfied for small $`ϵ`$ if $`g(x)`$ approaches a power $`x^n`$ for $`x\mathrm{}`$. In this case $`\varphi _ϵ`$ goes as $`ϵ^{1/n}`$, so $`V^{}(\varphi _ϵ)`$ goes as $`ϵ^{1/n}`$ and $`V^{\prime \prime }(\varphi _ϵ)`$ goes as $`ϵ^{2/n}`$, both of which can be made as small as we like by taking $`ϵ`$ sufficiently small.

In particular, if the singularity in $`g(x)`$ at $`x\mathrm{}`$ consists only of poles in $`1/x`$ of various orders up to $`n`$ (as is the case for a polynomial of order $`n`$) then the anthropically allowed range of $`\varphi `$ is

$$\left|\varphi \varphi _ϵ\right|_{\mathrm{max}}\frac{V_m}{V_1ϵ|g^{}(\varphi _ϵ)|}ϵ^{1/n}\left(\frac{V_m}{V_1}\right).$$
(15)
The flatness conditions make this range much greater than the Planck mass, but the fractional change in the a priori probability density (8) in this range is still very small

$$\left|\frac{V^{\prime \prime }(\varphi _ϵ)}{V^{}(\varphi _ϵ)}\right|\left|\varphi \varphi _ϵ\right|_{\mathrm{max}}\frac{V_m}{V_1}1.$$
(16)
To have a large fractional change in the a priori probability distribution in the anthropically allowed range for potentials of the second type that satisfy the flatness conditions, we need a function $`g(x)`$ that goes like a power as $`x\mathrm{}`$, but has a more complicated singularity at $`x=\mathrm{}`$ than just poles in $`1/x`$. An example is provided by the washboard potential with $`\alpha `$ and $`\beta `$ very small and $`\lambda `$ fixed, the case considered by Garriga and Vilenkin, for which $`g(x)`$ has an essential singularity at $`x=\mathrm{}`$.

In summary, the a priori probability is flat in the anthropically allowed range for several large classes of potentials, while it seems to be not flat only in exceptional cases.

It remains to consider whether the small parameters $`\lambda `$ or $`ϵ`$ in potentials respectively of the first or second kind could arise naturally. Garriga and Vilenkin argued that a term in a potential of what we have called the second kind with an over-all factor $`ϵ1`$ could be naturally produced by instanton effects. On the other hand, for potentials of type 1 a small parameter $`\lambda `$ could be naturally produced by the running of a field-renormalization factor. The field $`\varphi `$ has a conventional “canonical” normalization, as shown by the fact that the term $`\dot{\varphi }^2/2`$ in the vacuum energy (2) and the inertial term $`\ddot{\varphi }`$ in the field equation (3) have coefficients unity. Factors dependent on the ultraviolet cutoff will therefore be associated with external $`\varphi `$-lines. In order for the potential $`V(\varphi )`$ to be expressed in a cut-off independent way in terms of coupling parameters $`g_\mu `$ renormalized at a wave-number scale $`\mu `$, the field $`\varphi `$ must be accompanied with a field-renormalization factor $`Z_\mu ^{1/2}`$, which satisfies a differential equation of the form

$$\mu \frac{dZ_\mu }{d\mu }=\gamma (g_\mu )Z_\mu .$$
(17)
At very large distances, the field $`\varphi `$ will therefore be accompanied with a factor

$$\lambda =Z_0^{1/2}=\mathrm{exp}\left\{\frac{1}{2}_0^\mu \frac{d\mu ^{}}{\mu ^{}}\gamma (g_\mu ^{})\right\}Z_\mu ^{1/2}.$$
(18)
The integral here only has to be reasonably large and negative in order for $`\lambda `$ to be extremely small.

III. SLOW ROLLING IN THE EARLY UNIVERSE

When the cosmic energy density is dominated by the vacuum energy, the flatness conditions (6) and (7) insure that the vacuum energy changes little in a Hubble time. But if the vacuum energy density is nearly time-independent, then from the end of inflation until nearly the present it must have been much smaller than the energy density of matter and radiation, and under these conditions we are not able to neglect the inertial term $`\ddot{\varphi }`$ in Eq. (3). A separate argument is needed to show that the vacuum energy is nearly constant at these early times. This is important because, although there is no observational reason to require $`V(\varphi )`$ to be constant at early times, it must have been less than the energy of radiation at the time of nucleosynthesis in order not to interfere with the successful prediction of light element abundances, and therefore at this time must have been very much less than $`V_1`$, which we have supposed to be at least of order $`m_W^4`$. For potentials (9) of the first kind, this means that $`\varphi `$ must have been very close to its present value at the time of helium synthesis. Also, if $`\varphi `$ at the end of inflation were not the same as $`\varphi `$ at the time of galaxy formation, then a flat a priori distribution for the first would not in general imply a flat a priori distribution for the second.

At times between the end of inflation and the recent past the expansion rate behaved as $`H=\eta /t`$, where $`\eta =2/3`$ or $`\eta =1/2`$ during the eras of matter or radiation dominance, respectively. During this period, Eq. (3) takes the form

$$\ddot{\varphi }+\frac{3\eta }{t}\dot{\varphi }=V^{}(\varphi ),$$
(19)
If we tentatively assume that $`\varphi `$ is nearly constant, then Eq. (19) gives for its rate of change

$$\dot{\varphi }\frac{tV^{}(\varphi )}{1+3\eta }.$$
(20)
The change in the vacuum energy from the end of inflation to the present time $`t_0`$ is therefore

$$\mathrm{\Delta }V_0^{t_0}V^{}(\varphi )\dot{\varphi }𝑑t\frac{V^2(\varphi )t_0^2}{2(1+3\eta )}.$$
(21)
The present time is roughly given by $`t_0\eta \sqrt{3/8\pi G\rho _{V0}}`$, so the fractional change in the vacuum energy density since the end of inflation is

$$\left|\frac{\mathrm{\Delta }V}{\rho _{V0}}\right|\left(\frac{3\eta ^2}{2(1+3\eta )}\right)\left(\frac{V^2(\varphi )}{8\pi G\rho _{V0}^2}\right),$$
(22)
a subscript zero as usual denoting the present instant. The factor $`3\eta ^2/2(1+3\eta )`$ is of order unity, so the inequality (6) tells us that the change in the vacuum energy during the time since inflation has indeed been much less than its present value.

III. QUANTUM COSMOLOGY

In some theories of quantum cosmology the wave function of the universe is a superposition of terms, corresponding to universes with different (but time-independent) values for the vacuum energy $`\rho _V`$. It has been argued by Baum2, Hawking2 and Coleman10 that these terms are weighted with a $`\rho _V`$-dependent factor, that gives an a priori probability distribution with an infinite peak at $`\rho _V=0`$, but this claim has been challenged.11 As already acknowledged in references 4 and 5, if this peak at $`\rho _V=0`$ is really present, then anthropic considerations are both inapplicable and unnecessary in solving the problem of the cosmological constant.

Garriga and Vilenkin8 have proposed a different sort of infinite peak, arising from a $`\rho _V`$-dependent rate of nucleation of sub-universes operating over an infinite time. Even granting the existence of such a peak, it is not clear that it really leaves a vanishing normalized probability distribution at all other values of $`\rho _V`$. For instance, the nucleation rate might depend on the population of sub-universes already present, in such a way that the peaks in the probability distribution are kept to a finite size. If $`𝒫_{}(\rho _V)=0`$ except at the peak, then anthropic considerations are irrelevant and the cosmological constant problem is as bad as ever, since there is no known reason why the peak should occur in the very narrow range of $`\rho _V`$ that is anthropically allowed. On the other hand, if there is a smooth background in addition to a peak outside the anthropically allowed range of $`\rho _V`$ then the peak is irrelevant, because no observers would ever measure such values of $`\rho _V`$. In this case the probability distribution of the cosmological constant can be calculated using the methods of references 4 and 5.

ACKNOWLEDGEMENTS

I am grateful for a useful correspondence with Alex Vilenkin. This research was supported in part by the Robert A. Welch Foundation and NSF Grant PHY-9511632.

REFERENCES

1. A. Vilenkin, Phys. Rev. D27, 2848 (1983); A. D. Linde, Phys. Lett. B175, 395 (1986).
2. E. Baum, Phys. Lett. B133, 185 (1984); S. W. Hawking, in Shelter Island II - Proceedings of the 1983 Shelter Island Conference on Quantum Field Theory and the Fundamental Problems of Physics, ed. R. Jackiw et al. (MIT Press, Cambridge, 1995); Phys. Lett. B134, 403 (1984); S. Coleman, Nucl. Phys. B 307, 867 (1988).
3. An equation of this type was given by A. Vilenkin, Phys. Rev. Lett. 74, 846 (1995); and in Cosmological Constant and the Evolution of the Universe, K. Sato, et al., ed. (Universal Academy Press, Tokyo, 1996) (gr-qc/9512031), but it was not used in a calculation of the mean value or probability distribution of $`\rho _V`$.
4. S. Weinberg, in Critical Dialogs in Cosmology, ed. by N. Turok (World Scientific, Singapore, 1997).
5. H. Martel, P. Shapiro, and S. Weinberg, Ap. J. 492, 29 (1998).
6. S. Weinberg, Phys. Rev. Lett. 59, 2607 (1987).
7. J. D. Barrow and F. J. Tipler, The Anthropic Cosmological Principle (Clarendon Press, Oxford, 1986).
8. J. Garriga and A. Vilenkin, Tufts University preprint astro-ph/9908115, to be published.
9. P. J. Steinhardt and M. S. Turner, Phys. Rev. D29, 2162 (1984).
10. S. Coleman, Nucl. Phys. B 310, 643 (1988).
11. W. Fischler, I. Klebanov, J. Polchinski, and L. Susskind, Nucl. Phys. B237, 157 (1989).