text:
# 1 Introduction

## 1 Introduction

If one looks at textbooks to learn about the phase diagram of the two dimensional ($`2D`$) $`O(N)`$ models the situation seems clear: for $`N=2`$, there is a transition to a low temperature phase with only power law decay of correlations, whereas for the nonabelian case $`N>2`$ there is exponential decay at all temperatures. But while the first statement has been proven rigorously a long time ago , the second one remains an open mathematical question . The standard belief is rooted in the perturbative asymptotic freedom of the models for $`N>2`$; but over the years we have brought forth many reasons why we think it is unfounded . The absence of a mathematical proof together with ambiguous numerical results left the issue wide open.

In this paper we would like to present what we regard as convincing numerical evidence that in fact the $`2D`$ $`O(3)`$ model possesses a massless phase for sufficiently large $`\beta `$ and give a rigorous proof that this is incompatible with asymptotic freedom in the massive phase. We will also give a heuristic explanation of why and where the phase transition happens.

The models we are considering consist of classical spins $`s`$ taking values on the unit sphere $`S^{N1}`$, placed at the sites of a $`2D`$ regular lattice. These spins interact ferromagnetically with their nearest neighbors. Let $`ij`$ denote a pair of neighboring sites. We will consider two types of interactions between neighbouring spins:

* Standard action (s.a.): $`H_{ij}=s(i)s(j)`$
* Constrained action (c.a.): $`H_{ij}=s(i)s(j)`$ for $`s(i)s(j)c`$ and $`H_{ij}=\mathrm{}`$ for $`s(i)s(j)<c`$ for some $`c[1,1)`$.

The corresponding Gibbs measures are (for a finite lattice) given by

$$d\mu _{s.a.}=\frac{1}{Z}\underset{ij}{}e^{\beta H_{ij}}\underset{i}{}d\nu \left(s(i)\right)$$
(1)
for the standard action and

$$d\mu _{c.a.}=\frac{1}{Z}\underset{ij}{}\left[e^{\beta H_{ij}}\theta \left(s(i)s(j)c\right)\right]\underset{i}{}d\nu \left(s(i)\right)$$
(2)
for the constrained action, where $`d\nu `$ is the standard measure on the two sphere $`S^2`$ and the product $`_{ij}`$ is over nearest neighbors.

Almost a decade ago we showed  that one can rephrase the question of the existence of a soft phase in these models as a percolation problem and in fact this is the reason we introduced the c.a. model. It should be noted that the c.a. model shares with the s.a. model not only invariance under $`O(N)`$, but has also the same perturbative (= low temperature) expansion and the same ‘smooth’ classical solutions as the s.a. model. It is therefore to be expected that the s.a. and c.a. models fall in the same universality class (possess the same continuum limit) and, as we shall show shortly, the numerical evidence supports this expectation. The advantage of studying the c.a. model stems from the following fact: let $`ϵ_c=\sqrt{2(1c)}`$ and $`S_{ϵ_c}`$ the set of sites such that $`|sn|<ϵ_c/2`$ for some given unit vector $`n`$. Our rigorous result  was that if for some $`ϵ>ϵ_c`$ the set $`S_ϵ`$ on the triangular (T) lattice does not contain a percolating cluster, then the $`O(N)`$ model must be massless at that $`c`$. For the abelian $`O(2)`$ model we could prove the absence of percolation of this equatorial set $`S_{ϵ_c}`$ for $`c`$ sufficiently large  (modulo certain technical assumptions which were later eliminated by M. Aizenman ). For the nonabelian cases we could not give a rigorous proof. We did however present certain arguments  explaining why the percolating scenario seemed unlikely.

In this paper we will present numerical evidence that there is an $`ϵ_)`$ such that for $`ϵϵ_0`$ $`S_ϵ`$ does not percolate for any $`c`$; for sufficiently large $`c`$ this $`ϵ_0`$ will be larger than $`ϵ_c`$ and the model will thus be massless. We will also show that due to a rigorous inequality derived by us in the past , the existence of a finite $`\beta _{crt}`$ in the s.a. model on the square (S) lattice is incompatible with the presence of asymptotic freedom in the massive continuum limit of the model.

## 2 Percolation and masslessness

In this section we briefly review the special features of percolation in two dimensions and give a brief sketch of our argument relating percolation properties to the absence of a mass gap. We restrict the discussion to the T lattice; this keeps the arguments simpler because the T lattice is self-matching and no distinction has to be made between connectedness and $``$-connectedness (where points are also considered connected along diagonals).

The following two facts special to $`2D`$ are relevant for our discussion:

1. Noncoexistence of disjoint percolating sets: Let $`A`$ be the subset of the lattice defined by the spin lying in some subset $`𝒜S^{N1}`$ and $`\stackrel{~}{A}`$ its complement. Then with probability 1 $`A`$ and $`\stackrel{~}{A}`$ do not percolate at the same time. This has been proven rigorously only for special cases like Bernoulli percolation and the $`+`$ and $``$ clusters of the Ising model, but is believed to hold quite generally. (Aizenman  showed that in the case of $`O(2)`$ one does not need to invoke this principle).

2. Russo’s lemma : If neither $`A`$ nor its complement $`\stackrel{~}{A}`$ percolate, then the expected size of the cluster of $`A`$ attached to the origin, denoted by $`A`$, diverges; the same holds for its complement $`\stackrel{~}{A}`$. (In this simple form the lemma only holds for a self matching lattice like the T lattice). If $`\stackrel{~}{A}`$ percolates, then $`A`$ is expected to be finite.

The subsets of the sphere $`S^2`$ that interest us here are the following:

* ‘equatorial strip’ $`𝒮_ϵ`$, defined by $`|sn|<ϵ/2`$ for some fixed unit vector $`n`$.
* ‘upper polar cap’ $`𝒫_ϵ^+`$, defined by $`snϵ/2`$,
* ‘lower polar cap’ $`𝒫_ϵ^{}`$, defined by $`snϵ/2`$.
* ‘union of polar caps’ $`𝒫_ϵ=𝒫_ϵ^+𝒫_ϵ^{}`$.

The subsets of the lattice defined by these subsets of the sphere we denote by the corresponding roman letters $`S_ϵ`$ etc. and for brevity we say ‘a certain subset of the sphere percolates’ instead of ‘the subset of the lattice induced by a certain subset of the sphere percolates’ etc..

According to the discussion above, there are the following possibilities: either $`S_ϵ`$ percolates, or $`P_ϵ`$ percolates, or neither $`S_ϵ`$ nor $`P_ϵ`$ percolates and then both have divergent mean size (we shall call this third possibility in short formation of rings).

Let us now briefly review our argument  that relates percolation properties to the absence of a mass gap. Our statement was that if there was an equatorial strip $`𝒮_ϵ`$ that did not percolate for a certain $`c>1ϵ^2/2`$, there could be no mass gap in the system.

The argument is based on the imbedded Ising variables $`\sigma _i\mathrm{sgn}(s(i))`$. Using these variables, the s.a. Hamiltonian becomes:

$$H_{ij}=\sigma _i\sigma _j|s_{}(i)s_{}(j)|s_{}(i)s_{}(j)$$
(3)
where $`s_{}(i)=s(i)n`$ and $`s_{}(i)=n\times (s(i)\times n)`$. The c.a. model can be similarly described in terms of the variables $`\sigma _i`$, $`|s_{}(i)|`$ and $`s_{}(i)`$. In both models one thus obtains an induced Ising model for which the Fortuin-Kastleyn (FK) representation  is applicable. In this representation the Ising system is mapped into a bond percolation problem: In the s.a. model a bond is placed between any like neighboring Ising spins with probability $`p=1\mathrm{exp}(2\beta s_{}(i)s_{}(j))`$. For the c.a. model a bond is also placed if after flipping one of the two neighboring Ising spins the constraint $`s(i)s(j)c`$ is violated. From the FK representation is follows that the mean cluster size of the site clusters joined by occupied bonds (FK-clusters) is equal to the Ising magnetic susceptibility. In a massive phase the latter must remain finite. Hence, if the FK-clusters have divergent mean size, the original $`O(3)`$ ferromagnet must be massless (the Ising variables $`\sigma `$ are local functions of the originally spin variables $`s`$).

Now notice that by construction for the c.a. model the FK-clusters with, say, $`\sigma =+1`$ must contain all sites with $`s(i)n>\sqrt{(1c)/2}`$. Therefore the model must be massless if clusters obeying this condition have divergent mean size. But the polar set $`P_ϵ`$ consists of the two disjoint components $`P_ϵ^+`$ and $`P_ϵ^{}`$. For $`c>1ϵ^2/2`$ there are no clusters containing elements of both $`P_ϵ^+`$ and $`P_ϵ^{}`$. Hence if for such values of $`c`$ clusters of $`P_ϵ`$ form rings, so do clusters of $`P_ϵ^+`$ separately and the $`O(3)`$ model must be massless by the argument just given.

If we want to study the percolation or absence of percolation of the set $`A`$ corresponding to a subset $`𝒜S^{N1}`$ of positive measure numerically, we have to consider a sequence of tori of increasing size $`L`$. On these tori we measure the mean cluster size of of $`A`$. If $`A`$ percolates in the thermodynamic limit, by translation invariance $`A=O(L^2)`$; if its complement percolates $`A`$ should approach a finite nonzero value, and if $`A`$ forms rings we expect $`A=O(L^{2\eta })`$ for some $`\eta >0`$. Therefore, if we define the ratio

$$r=P_ϵ/S_ϵ,$$
(4)
for $`L\mathrm{}`$ it should either go to 0 if $`S_ϵ`$ percolates or to $`\mathrm{}`$ if $`P_ϵ`$ percolates; if neither $`S_ϵ`$ nor $`P_ϵ`$ percolates, then both form rings and the ratio $`r`$ could diverge, go to 0 or approach some finite, nonzero value depending upon the value of the critical index $`\eta `$ for the two types of clusters.

In the next section we will describe what our numerical simulations tell us about the percolation properties of equatorial strips and polar caps.

## 3 Main numerical results

Our results were obtained from a Monte Carlo (MC) investigation using an $`O(3)`$ version of the Swendsen-Wang cluster algorithm  and consist of a minimum of 20,000 lattice configurations used for taking measurements. For each value of $`ϵ`$ we studied $`L=160`$, 320 and 640 (for $`ϵ=.78`$ we also studied $`L=1280`$).

In Fig.1 we show the numerical value of the ratio $`r`$ as function of $`c`$ for $`\beta =0`$ for four values of $`ϵ`$ for the c.a. model on a T lattice. Three distinct regimes are manifest for each of the four values of $`ϵ`$ investigated:

* For small $`c`$, $`r`$ is increasing with $`L`$, presumably diverging to $`\mathrm{}`$ (region 1).
* For intermediate $`c`$, $`r`$ is decreasing with $`L`$, presumably converging to 0 (region 2).
* For $`c`$ sufficiently large depending upon $`ϵ`$, $`r`$ shows a very mild dependence upon $`L`$.

This can only mean that for these values of $`ϵ`$ for small $`c`$ $`P_ϵ`$ percolates, for intermediate $`c`$ $`S_ϵ`$ percolates and for sufficiently large $`c`$ both $`P_ϵ`$ and $`S_ϵ`$ form rings with quite similar (possibly equal) values of $`\eta `$. Our data allow to deduce a semiquantitative ‘phase diagram’ in the $`(c,ϵ)`$-plane of the percolation problem induced by the c.a. model on the T lattice for $`\beta =0`$.

This is shown in Fig.2. The solid line C is the curve $`c=1ϵ^2/2`$; above that line the two polar caps cannot touch and therefore their union cannot percolate. The dashed line D represents the minimal equatorial width above which $`S_ϵ`$ percolates. The point T at the intersection of the curves D and C gives an upper bound for $`c_{crt}`$, the value of $`c`$ above which the c.a. model is massless.

Let us explain how this picture was obtained: since for $`c=1`$ (no constraint) the model reduces to independent site percolation, for which the percolation threshold is known rigorously to be $`ϵ=1`$, curve D has to start at $`ϵ=1`$. With increasing $`c`$ that threshold shifts to smaller values of $`ϵ`$. Four points of the dashed line are in fact determined by the data displayed in Fig.1: the clearly identifiable four points where the lines for different lattice sizes $`L`$ cross and the ratio $`r`$ becomes independent of $`L`$ determine the percolation threshold for the chosen value of $`ϵ`$ and so determine a point of the dashed line D.

Two features of this diagram are worth emphasizing:
1. An equatorial strip of width less than approximately $`ϵ=.76`$ never percolates.
2. In Fig.1 for approximately $`c>0.4`$ a range of $`ϵ`$’s appears such that the ratio $`r`$ again becomes approximately independent of $`L`$. This is signalling the appearance of a new ‘phase’ in which both $`S_ϵ`$ and $`P_ϵ`$ form rings (the dotted line separates it from the region of percolation of $`P_ϵ`$).

This regime of ring formation of both $`S_ϵ`$ and $`P_ϵ`$ is lying between the dotted and the dashed lines in Fig.2. Our data give strong evidence of its existence, but they do not determine in detail where the boundaries are. The dotted line has to run to $`ϵ=0`$ for $`c=1`$ because below it there is percolation of $`P_ϵ`$, and this is not possible above the solid line C $`ϵ=ϵ_c`$ (because it would conflict with the principle of non-coexistence of disjoint percolating sets). We drew the dashed line into upper right corner because we expect that eventually, for $`c`$ approaching 1, any polar cap will start forming rings, thereby preventing percolation of the corresponding equatorial strips.

But what it is essential for our conclusion that there is a massless phase is only that there is a regime below the dashed line D and above the solid line C, in which $`S_ϵ`$ does not percolate and the two polar caps do not touch. In other words, the lines C and D have to cross (the crossing point is denoted by T in Fig.2). Since we found that for $`ϵ<ϵ_0=0.76`$ $`𝒮_ϵ`$ never percolates, this means that for $`cc_{ϵ_0}=0.71`$ the model is massless. In fact the massless phase must start earlier, and for instance based on our data we estimate that at c=0.61 the c.a. model is already massless.

In the next section we will further corroborate the fact that for $`ϵ<0.76`$ the equatorial strip does not percolate for any $`c`$.

We would like to comment briefly on another recent paper dealing with percolation properties of equatorial strips in the $`O(3)`$ model: Allès et al  published a study showing that for $`ϵ=1.05`$ and $`\beta =2.0`$ in the s.a. model $`S_ϵ`$ percolates. Although strictly speaking our percolation argument applies only to the c.a. model, the result of Allès et al is not surprising since at $`\beta =2.0`$ the s.a. model is clearly in its massive phase , hence, by analogy with what happens in the c.a. model, one would expect that clusters of a sufficiently wide equatorial strip percolate (see ). The real issue, which the authors of  did not seem to appreciate, is whether in the c.a. model clusters of the equatorial strip $`S_{ϵ_c}`$ continue to percolate for $`c`$ sufficiently close to 1. The numerics presented in Fig.1 suggest that that is not the case.

## 4 Corroborating numerics

To corroborate our most important result, namely that for approximately $`ϵ<.76`$ $`S_ϵ`$ does not percolate for any value of $`c`$, we also measured (at $`\beta =0`$) the ratio of the mean cluster size of the set $`P_ϵ^{}^+`$ with $`ϵ^{}=.5`$ to that of the set $`S_ϵ`$ with $`ϵ=0.75`$ ($`ϵ^{}`$ was chosen so that $`P_ϵ^{}^+`$ has equal density with $`S_ϵ`$). The results are shown in (Fig.3). This figure shows that for $`c`$ less than about 0.4 (and greater than 0) the ratio grows very rapidly with $`L`$, indicating that $`P_ϵ^{}^+`$ forms rings while $`S_ϵ`$ has finite mean size; this region terminates around $`c=0.4`$, where presumably also $`S_ϵ`$ starts forming rings, and the dependence of the ratio upon $`L`$ becomes much milder. Since for $`c>0.4`$ the ratio continues to grow with $`L`$, at equal density, clusters of the polar cap are larger than those of the equatorial strip. The larger average cluster size of the polar cap compared to the strip of the same area is probably due to the fact that the strip has a larger boundary than the polar cap. This is in agreement with a general conjecture stated in , namely that for $`c`$ sufficiently large, if two sets have equal area but different perimeters, the one with the smaller perimeter will eventually, for $`c`$ approaching 1, have larger average cluster size. For the case at hand, this is apparently true for all values of $`c`$.

## 5 Comparison to the $`O(2)`$ model

The general belief, which we criticized in ref., is that there is a fundamental difference between abelian and nonabelian models. To test this belief we also studied the ratio $`r`$ for the c.a. $`O(2)`$ model on the T lattice. The phase diagram is shown in Fig.4. Since in the $`O(2)`$ model the set $`𝒫_ϵ`$ can also be regarded as a set $`𝒮_{\stackrel{~}{ϵ}}`$ where $`\stackrel{~}{ϵ}=\sqrt{4ϵ^2}`$, certain features of that diagram follow from rigorous arguments. For instance it is clear that in the c.a. model there exist two intersecting curves C and $`\stackrel{~}{\mathrm{C}}`$ and in the region to their right the model must be massless . The precise location of the curves $`D`$ (or $`\stackrel{~}{D}`$) must be determined numerically, something which we did not do. We did verify though that the ring formation region begins around $`c=0.5`$.

## 6 Universality between the s.a. and c.a. models

In our opinion the arguments and numerical evidence provided so far give strong indications that the c.a. $`O(3)`$ model on the T lattice has a massless phase. Universality would suggest that a similar situation must exist for the s.a. models on the T and S lattices. To test universality we measured on the S lattice the renormalized coupling both on thermodynamic lattices in the massive phase and in finite volume in the presumed critical regime (as in ). Our data for the c.a. model on the S lattice only determine an interval (about $`.5`$ to $`.7`$) in which the massless phase of the model sets in; we tried to see if we could get a similar $`L`$ dependence for the renormalized coupling in the s.a. model at a suitable $`\beta `$ as for $`c=.61`$ in the c.a. model at $`\beta =0`$. This seems to be indeed the case for $`\beta `$ roughly $`3.4`$. We went only up to $`L=640`$, hence this equivalence between $`c`$ and $`\beta `$ should be considered only as a rough approximation, but there seems to be no doubt that the two models have the same continuum limit.

We also compared the step scaling curve in the s.a. and c.a. models. The step scaling curve is obtained as follows: on a periodic lattice of size $`L\times L`$ we define an apparent correlation length $`\xi (L)`$. Namely let $`P=(p,0)`$, $`p=\frac{2n\pi }{L}`$, $`n=0,1,2,\mathrm{},L1`$. Then define

$$\xi (L)=\frac{\sqrt{3}}{4\mathrm{sin}(\pi /L)}\sqrt{G(0)/G(1)1}$$
(5)
where

$$G(p)=\frac{1}{L^2}|\widehat{s}(P)|^2;\widehat{s}(P)=\underset{x}{}e^{iPx}s(x)$$
(6)
Leaving $`\beta `$ respectively $`c`$ fixed, one doubles $`L`$ and measures also $`\xi (2L)`$. The step scaling function gives the ratio $`2\xi (L)/\xi (2L)`$ versus $`L/\xi (L)`$. In the continuum limit ($`L\mathrm{}`$ and $`\beta \beta _{\mathrm{𝑐𝑟𝑡}}`$ at $`L/\xi (L)`$ fixed), this procedure produces a unique curve characterizing the universality class of the model. In Fig.5 we present the step scaling function for the c.a. and s.a. models. The data were produced by adjusting $`\beta `$ and $`c`$ so that at $`L=20`$ we obtain roughly the same $`\xi (L)`$ in the two models. After that, leaving $`\beta `$ respectively $`c`$ fixed, $`L`$ was doubled until $`L=320`$. As can be seen, the two step scaling curves agree reasonably well; the slight disagreement is probably due to the fact that the two curves have to agree only in the continuum limit ($`\xi \mathrm{}`$), i.e. they have different lattice artefacts, whereas the largest value of $`\xi `$ reached was only approximately 35 lattice units.

## 7 Heuristic explanation of the transition

It is intersting to note that there is a heuristic explanation for both the existence of a massless phase in the s.a. $`O(3)`$ model and for the value of $`\beta _{\mathrm{𝑐𝑟𝑡}}`$. Indeed it is known rigorously that in $`2D`$ a continuous symmetry cannot be broken at any finite $`\beta `$. In a previuos paper  we showed that the dominant configurations at large $`\beta `$ are not instantons but superinstantons (s.i.). In principle both instantons and s.i. could enforce the $`O(3)`$ symmetry. In a box of diameter $`R`$ the former have a minimal energy $`E_{\mathrm{𝑖𝑛𝑠𝑡}}=4\pi `$  while the latter $`E_{s.i.}=\delta ^2\pi /\mathrm{ln}R`$, where $`\delta `$ is the angle by which the spin has rotated over the distance $`R`$. Now suppose that $`\beta _{\mathrm{𝑐𝑟𝑡}}`$ is sufficiently large for classical configurations to be dominant. Then let us choose $`\delta =2\pi `$ (restoration of symmetry) and ask how large must $`R`$ be so that the superinstanton configuration has the same energy as one instanton. One finds $`\mathrm{ln}R=\pi ^2`$. But in the Gaussian approximation

$$s(0)s(x)1\frac{1}{\beta \pi }\mathrm{ln}x$$
(7)
Thus restoration of symmetry occurs for $`\mathrm{ln}x\pi \beta `$. This simpleminded argument suggests that for $`\beta \pi `$ instantons become less important than s.i.. Now in a gas of s.i. the image of any small patch of the sphere forms rings, hence the system is massless. While this is not a quantitative argument, we believe it captures qualitatively what happens: a transition from localized defects (instantons) to super-instantons.

## 8 Absence of asymptotic freedom

Next let us discuss the connection between a finite $`\beta _{crt}`$ and the absence of asymptotic freedom. It follows from our earlier work concerning the conformal properties of the critical $`O(2)`$ model . We refer the reader for details to that paper and give only an outline of the argument. The s.a. lattice $`O(N)`$ model possesses a conserved isospin current. This currrent can be decomposed into a transverse and longitudinal part. Let $`F^T(p)`$ and $`F^L(p)`$ denote the thermodynamic values of the 2-point functions of the transverse and longitudinal parts at momentum $`p`$, respectively. Using reflection positivity and a Ward identity we proved that in the massive continuum limit the following inequalities must hold for any $`p0`$:

$$0F^T(p)F^T(0)=F^L(0)F^L(p)=2\beta E/N$$
(8)
Here $`E`$ is the expectation value of the energy

$$E=s(i)s(j)$$

at inverse temperature $`\beta `$. Since $`E1`$ it follws that if $`\beta _{\mathrm{𝑐𝑟𝑡}}<\mathrm{}`$ $`F^T(0)F^T(p)`$ cannot diverge for $`p\mathrm{}`$ as required by perturbative asymptotic freedom . In fact, for $`\beta _{crt}=3.4`$ (which is a reasonable guess) $`F^T(p)`$ must be less than 2.27, in violation of the form factor computation giving $`F^T(0)F^T(\mathrm{})>3.651`$ .

## 9 Concluding remarks

Since the implications of our result, that for the c.a. model a sufficiently narrow equatorial strip never percolates, are so dramatic, the reader may wonder how credible are the numerics. The only debatable point is whether our results represent the true thermodynamic behaviour for $`L\mathrm{}`$ or are merely small volume artefacts. While we cannot rule out rigorously the latter possibility, certain features of the data make it highly implausible:

* Small volume effects should set in gradually, while the data in Fig.1 indicate a rather abrupt change from a region where $`r`$ is decreasing with $`L`$ to one where $`r`$ is essentially independent of $`L`$.
* For $`c1`$ at fixed $`L`$, $`r`$ must approach the ‘geometric’ value $`r=2/ϵ1`$. As can be seen, in all the cases studied, throughout the ‘ring’ region $`r`$ is clearly larger than this value, while it should go to 0 if $`S_ϵ`$ percolated.
* In Fig.3 there is no indication of the ratio going to 0 for increasing $`L`$. Moreover the dramatic change in slope around $`c=.4`$ indicates that the polar cap $`P_ϵ^{}`$ starts forming rings at a smaller value of $`c`$ than the equatorial strip $`S_ϵ`$.

We have additional numerical evidence that clusters of a polar cap $`P_ϵ^+`$ smaller than a hemisphere ($`sn>ϵ/2>0`$) form rings for some $`c<1`$. Namely we investigated the case $`ϵ=0.1`$. For the the case $`c=1`$ (Bernoulli percolation) it is known rigorously that clusters of this set have finite mean size. As can be seen from Fig.6 our numerical values at $`c=1`$ corroborate this fact. In the same figure we show the mean cluster size of cluster of $`P_{0.1}^+`$ at $`c=0`$, where the correlation length is approximately 53 lattice units. Even though we increased $`L`$ up to 1280, the mean cluster size shows no sign of leveling off, growing in fact like some power of $`L`$, consistent with the formation of rings.

Therefore there is good numerical evidence that for $`c=0`$, where we can reach the thermodynamic limit, clusters of this polar cap form rings. The natural expectation would be that the mean cluster size of a subset of a hemisphere is a nondecreasing function of $`c`$. This is borne out by the numerics, as shown in Fig.7. There we represent the mean cluster size of $`P_{0.1}^+`$ at fixed $`L=640`$ function of $`c`$. The data support the assertion that for any $`c>0`$ the mean size of the clusters of $`P_{0.1}^+`$ diverges, which, via our argument, implies that the c.a. model must cease being massive for some $`c<1`$.

Thus we doubt very much that the effects we are seeing represent small volume artefacts. Moreover, if $`s_z`$, the $`z`$-component of the spin $`s`$ remained massive at low temperature and in fact an arbitrarily narrow equatorial strip percolated, one would have to explain away our old paradox : if such a narrow strip percolated, an even larger strip would percolate and on it one would have an induced $`O(2)`$ model in its massless phase, in contradiction to the Mermin-Wagner theorem.

Consequently it seems unavoidable to conclude that the phase diagram in Fig.2 represents the truth, that a soft phase exists both in the s.a. and the c.a. model and that the massive continuum limit of the $`O(3)`$ model is not asymptotically free. In a previous paper  we have already shown that in nonabelian models even at short distances perturbation theory produces ambiguous answers. The present result sharpens that result by eliminating the possibility of asymptotic freedom in the massive continuum limit.

Acknowledgement: AP is grateful to the Humboldt Foundation for a Senior US Scientist Award and to the Werner-Heisenberg-Institut for its hospitality.


format: text

text:
# Nuclear matter properties in relativistic mean field model with 𝝈-𝝎 coupling

## Abstract

The possibility of extending the linear $`\sigma `$-$`\omega `$ model by introducing a $`\sigma `$-$`\omega `$ coupling phenomenologically is explored. It is shown that, in contrast to the usual Walecka model, not only the effective nucleon mass $`M^{}`$ but also the effective $`\sigma `$ meson mass $`m_\sigma ^{}`$ and the effective $`\omega `$ meson mass $`m_\omega ^{}`$ are nucleon density dependent. When the model parameters are fitted to the nuclear saturation point (the nuclear radius constant $`r_0=1.14`$fm and volume energy $`a_1=16.0`$MeV) as well as to the effective nucleon mass $`M^{}=0.85M`$, the model yields $`m_\sigma ^{}=1.09m_\sigma `$ and $`m_\omega ^{}=0.90m_\omega `$ at the saturation point, and the nuclear incompressibility $`K_0=501`$MeV. The lowest value of $`K_0`$ given by this model by adjusting the model parameters is around $`227`$MeV.

Keywords: relativistic mean field theory, $`\sigma `$-$`\omega `$ coupling, nuclear matter properties

As the starting point for the microscopic relativistic description of nuclear many-body system, within the framework of quantum hadrodynamics, the well-studied linear $`\sigma `$-$`\omega `$ model has been proved to be able to describe the saturation and other properties of nuclear matter. However, this model yields the nuclear incompressibility $`K_0`$ around $`550`$MeV which is unacceptably high, and also the effective nucleon mass $`M^{}`$ around $`0.54M`$ which seems uncomfortably low. In order to remedy this situation, many works have been done on the extension or generalization of this model. Further inclusion of nonlinear self-interaction of $`\sigma `$ mesons has been proved to be successful in reducing the nuclear incompressibility $`K_0`$ significantly, in spite of some controversy about the stability of the system specified by the fitted parameters. On the other hand, the density-dependent N-$`\sigma `$ interaction given by the Zimanyi-Moszkowski model is shown to be able to give reasonable nuclear incompressibility $`K_0`$ and effective nucleon mass $`M^{}`$ as well. Since the Zimanyi-Moszkowski model contains a coupling between $`\sigma `$ and $`\omega `$ mesons, it is interesting to investigate what is the direct effect of the coupling between $`\sigma `$ and $`\omega `$ mesons, when detached from the density-dependence of the N-$`\sigma `$ coupling. The purpose of this letter is to explore the possibility of extending the linear $`\sigma `$-$`\omega `$ model by introducing a $`\sigma `$-$`\omega `$ coupling phenomenologically, and to investigate whether this coupling alone is able to reduce the nuclear incompressibility $`K_0`$ and in the same time to increase the effective nucleon mass $`M^{}`$.

Let us start with the following Lagrangian density:

$$=\overline{\psi }[\gamma _\mu (i^\mu g_\omega \omega ^\mu )(Mg_\sigma \varphi )]\psi $$

$$+\frac{1}{2}\left[(_\mu \eta g_\omega ^{}\omega _\mu )\varphi (^\mu +\eta g_\omega ^{}\omega ^\mu )\varphi m_\sigma ^2\varphi ^2\right]$$

$$\frac{1}{4}F_{\mu \nu }F^{\mu \nu }+\frac{1}{2}m_\omega ^2\omega _\mu \omega ^\mu ,\eta =i\text{and}\mathrm{\hspace{0.17em}1},$$
(1)
where $`F^{\mu \nu }=^\mu \omega ^\nu ^\nu \omega ^\mu ,`$ $`\psi `$, $`\varphi `$ and $`\omega `$ are the nucleon, $`\sigma `$ and $`\omega `$ meson fields with masses $`M`$, $`m_\sigma `$ and $`m_\omega `$, respectively, while $`g_\sigma `$, $`g_\omega `$ are the respective coupling constants, and $`g_\omega ^{}`$ is the $`\sigma `$-$`\omega `$ coupling constant. $`\eta =i`$ gives a $`\sigma `$-$`\omega `$ coupling similar to that introduced by usual covariant derivatives, while $`\eta =1`$ gives a gauge-like virtual coupling similar to the imaginary coupling. Mathematically, it is equivalent to add a term $`\eta ^2g_\omega ^2\omega _\mu \omega ^\mu `$ to $`m_\sigma ^2`$, which is the same as the $`\eta _2`$ term introduced in the nuclear effective field theory. The present model reduces to the original Walecka model when the $`\sigma `$-$`\omega `$ coupling constant is zero, $`g_\omega ^{}=0`$.

For static nuclear matter, the field equations derived from this Lagrangian are reduced to the following equations in the mean field approximation:

$$(i\gamma ^\mu _\mu g_\omega \gamma ^0\omega _0M^{})\psi =0,$$
(2)
$$m_\sigma ^2\varphi =g_\sigma \rho _s,$$
(3)
$$m_\omega ^2\omega _0=g_\omega \rho _N,$$
(4)
where the effective nucleon mass $`M^{}`$, effective $`\sigma `$ meson mass $`m_\sigma ^{}`$ and effective $`\omega `$ meson mass $`m_\omega ^{}`$ are defined respectively as

$$M^{}=Mg_\sigma \varphi ,$$
(5)
$$m_\sigma ^2=m_\sigma ^2+\eta ^2g_\omega ^2\omega _0^2,$$
(6)
$$m_\omega ^2=m_\omega ^2\eta ^2g_\omega ^2\varphi ^2,$$
(7)
and $`\rho _s=\overline{\psi }\psi `$ is the scalar density, $`\rho _N=\overline{\psi }\gamma ^0\psi `$ the baryon density. For calculating these effective masses, the above three equations can be rewritten as the following self-consistent equations:

$$\xi =\frac{M^{}}{M}=1\frac{\alpha }{s}\frac{\rho _s}{\rho _0},$$
(8)
$$s=\xi _\sigma ^2=\frac{m_\sigma ^2}{m_\sigma ^2}=1+\frac{\alpha _\sigma }{v^2}\frac{\rho _N^2}{\rho _0^2},$$
(9)
$$v=\xi _\omega ^2=\frac{m_\omega ^2}{m_\omega ^2}=1\frac{\alpha _\omega }{s^2}\frac{\rho _s^2}{\rho _0^2},$$
(10)
where $`\rho _0`$ is the standard nucleon number density, $`\alpha `$, $`\alpha _\sigma `$ and $`\alpha _\omega `$ are the dimensionless composite parameters defined respectively as

$$\alpha =\frac{g_\sigma ^2\rho _0}{m_\sigma ^2M},\alpha _\sigma =\frac{\eta ^2g_\omega ^2g_\omega ^2\rho _0^2}{m_\omega ^4m_\sigma ^2},\alpha _\omega =\frac{\eta ^2g_\omega ^2g_\sigma ^2\rho _0^2}{m_\sigma ^4m_\omega ^2}.$$
(11)
It should be noted that $`\alpha _\sigma `$ and $`\alpha _\omega `$ are positive for $`\eta =1`$ while negative for $`\eta =i`$.

The nuclear matter equation of state derived from Lagrangian density (1) can be expressed in terms of the nuclear energy density $``$ as $`e=/\rho _NM`$, and

$$=_k+_\sigma +_\omega ,$$
(12)
$$_k=\frac{M^4\xi ^4}{\pi ^2}\underset{i=p,n}{}F_1(k_i/\xi M),$$
(13)
$$_\sigma =\frac{1}{2}(1\xi )M\rho _s,$$
(14)
$$_\omega =\frac{1}{2}(1\xi )yM\rho _s,y=\frac{(s1)(2v1)}{s(1v)},$$
(15)
where $`k_p`$ and $`k_n`$ are the proton and neutron Fermi momenta respectively, and the function $`F_m(x)`$ is defined as

$$F_m(x)=_0^x𝑑xx^{2m}\sqrt{1+x^2}.$$
(16)
The baryon density $`\rho _N`$ and scalar density $`\rho _s`$ can be expressed as

$$\rho _N=\frac{1}{3\pi ^2}\underset{i=p,n}{}k_i^3,$$
(17)
$$\rho _s=\frac{M^3\xi ^3}{\pi ^2}\underset{i=p,n}{}f_1(k_i/\xi M),$$
(18)
where the function $`f_m(x)`$ is defined as

$$f_m(x)=_0^x𝑑x\frac{x^{2m}}{\sqrt{1+x^2}}.$$
(19)
Having the equation of state, the pressure $`p`$ can be derived as

$$p=+\rho _N\frac{}{\rho _N}=\frac{1}{3}(_k_\sigma M\rho _s)+_\omega \frac{2(1s)}{s}_\sigma .$$
(20)
Instead of Fermi momenta $`k_p`$ and $`k_n`$, we will use the nucleon density $`\rho _N=\rho _n+\rho _p`$ and relative neutron excess $`\delta =(\rho _n\rho _p)/\rho _N`$ as independent variables for describing the nuclear matter. At the standard state $`\rho _N=\rho _0`$, $`\delta =0`$, the pressure should be zero,

$$p(\rho _0,0)=0.$$
(21)
In addition, the depth of the equation of state is related to the nuclear volume energy $`a_1`$ as

$$e(\rho _0,0)=a_1.$$
(22)
The standard nucleon number density $`\rho _0`$ defined by Eq.(21) is related to nuclear radius constant $`r_0`$ and Fermi momentum $`k_F`$ of standard nuclear matter as $`\rho _0=3/4\pi r_0^3=2k_F^3/3\pi ^2`$.

The generalized nuclear incompressibility $`K(\rho _N,\delta )`$ can be defined and derived as

$$K9\frac{p}{\rho _N}=\frac{3(p+)+12_\omega }{\rho _N}\left[3M\rho _s+\frac{12(32s)}{s}_\sigma \right]\frac{1}{\rho _s}\frac{\rho _s}{\rho _N}+\frac{12(3s)}{s}\frac{_\sigma }{s}\frac{s}{\rho _N}\frac{12(1v)}{12v}\frac{_\omega }{v}\frac{v}{\rho _N},$$
(23)
$$\rho _N\frac{s}{\rho _N}=2(1s)\left[1\frac{2(1v)}{(1\xi )v}\rho _N\frac{\xi }{\rho _N}\right],$$
(24)
$$\rho _N\frac{v}{\rho _N}=\frac{2(1v)}{1\xi }\rho _N\frac{\xi }{\rho _N},$$
(25)
$$\rho _N\frac{\rho _s}{\rho _N}=Q+3(\rho _sQ)\frac{\rho _N}{\xi }\frac{\xi }{\rho _N},$$
(26)
$$\rho _N\frac{\xi }{\rho _N}=\frac{1}{s}\left[2(1\xi )(1s)+\frac{\alpha Q}{\rho _0}\right]\left[1\frac{4(1s)(1v)}{sv}+\frac{3(1\xi )(\rho _sQ)}{\xi \rho _s}\right]^1,$$
(27)
$$Q=\frac{M^3\xi ^3}{3\pi ^2}\underset{i=p.n}{}\frac{k_i}{\xi M}f_1^{}(k_i/\xi M).$$
(28)
The usual nuclear incompressibility $`K_0`$ can be obtained from $`K(\rho _N,\delta )`$ as

$$K_0=K(\rho _0,0)=9\left(\rho _N^2\frac{^2e}{\rho _N^2}\right)_0,$$
(29)
the subscript $`0`$ stands for the standard state $`\rho _N=\rho _0`$ and $`\delta =0`$. In addition, the following formula for symmetry energy $`J`$ can be derived:

$$J\frac{1}{2}\frac{^2e}{\delta ^2}|_0=\frac{1}{6}\frac{k_F^2}{\sqrt{k_F^2+M^2\xi _0^2}}.$$
(30)
Around the standard state $`\rho _N=\rho _0,\delta =0`$, the nuclear equation of state and thus the nuclear matter properties are specified essentially by the standard density $`\rho _0`$, volume energy $`a_1`$, symmetry energy $`J`$, incompressibility $`K_0`$, density symmetry $`L`$ and symmetry incompressibility $`K_s`$ , where $`L`$ and $`K_s`$ are, respectively,

$$L\frac{3}{2}\left(\rho _N\frac{^3e}{\rho _N\delta ^2}\right)_0=\frac{3}{2\rho _0}\frac{^2p}{\delta ^2}|_0,$$
(31)
$$K_s\frac{9}{2}\left(\rho _N^2\frac{^4e}{\rho _N^2\delta ^2}\right)_0=6L+\frac{1}{2}\frac{^2K}{\delta ^2}|_0.$$
(32)
For the description of nuclear matter, there are three independent parameters $`\alpha `$, $`\alpha _\sigma `$ and $`\alpha _\omega `$ in the present model. The saturation point of nuclear matter ($`\rho _0`$,$`a_1`$) defined by Eqs.(21) and (22) can be used as input to determine $`\alpha _\sigma `$ and $`\alpha _\omega `$ as function of $`\alpha `$. The way is as follows. At the standard state ($`\rho _N=\rho _0`$, $`\delta =0`$), these two equations can be solved for $`s`$ and $`y`$,

$$s=\frac{3(1\xi )M\rho _s}{3\rho _0(Ma_1)2_k\xi M\rho _s},$$
(33)
$$y=\frac{2\rho _0(Ma_1)2_k(1\xi )M\rho _s}{(1\xi )M\rho _s},$$
(34)
and thus from Eq.(15)

$$v=\frac{sy+s1}{sy+2s2}.$$
(35)
For given $`\rho _0`$ and $`a_1`$, $`s`$ and $`v`$ can be calculated by Eqs.(33)-(35) for a chosen $`\xi `$. Then $`\alpha `$, $`\alpha _\sigma `$ and $`\alpha _\omega `$ can be calculated by Eqs.(8)-(10) for a chosen $`\xi `$. In our calculation, the nuclear radius constant $`r_0=1.14`$fm which corresponds to the standard nucleon number density $`\rho _0=0.161`$fm<sup>-3</sup>, the nuclear volume energy $`a_1=16.0`$MeV, nucleon mass $`M=938.9`$MeV and constant $`\mathrm{}c=197.327053`$MeV$``$fm are used.

Fig.1 plots $`\alpha _\sigma `$ and $`\alpha _\omega `$ as function of $`\alpha `$. At the point $`\alpha =\alpha _W=0.4908`$ (indicated by an arrow), we have $`\alpha _\sigma =\alpha _\omega =0`$, and the present model reduces to Walecka model. For $`\alpha >\alpha _W`$, $`\alpha _\sigma `$ and $`\alpha _\omega `$ are negative, which corresponds to $`\eta =i`$. For $`\alpha <\alpha _W`$, $`\alpha _\sigma `$ and $`\alpha _\omega `$ are positive, which corresponds to $`\eta =1`$. There is a lower limit $`\alpha _{min}=0.1037`$, below this point there is no physical solution, as the effective $`\omega `$ meson mass becomes imaginary, $`\xi _\omega ^2<0`$. Fig.2 shows the dimensionless effective nucleon mass $`\xi `$, the dimensionless effective $`\sigma `$ meson mass $`\xi _\sigma `$ and the dimensionless effective $`\omega `$ meson mass $`\xi _\omega `$ as function of $`\alpha `$. At the point $`\alpha _W`$, $`\xi =0.5437`$ and $`\xi _\sigma =\xi _\omega =1`$. From $`\alpha _W`$ to the right, $`\xi `$ decreases as $`\alpha `$ increases, and we have $`\xi _\sigma <1`$ and $`\xi _\omega >1`$. From $`\alpha _W`$ to the left, $`\xi `$ increases as $`\alpha `$ decreases, and we have $`\xi _\sigma >1`$ and $`\xi _\omega <1`$. Fig.3 gives $`K_0`$, $`J`$, $`L`$, and $`K_s`$ as function of $`\alpha `$. At $`\alpha _W`$, $`K_0=553`$MeV, $`J=20.2`$MeV, $`L=70.6`$MeV, and $`K_s=88`$MeV. As $`\alpha `$ decreases, $`K_0`$ decreases at first and then increases slowly, passing through a very flat plateau, and finally decreases rapidly to a lower limit $`227`$MeV. $`J`$ and $`L`$ as well as $`K_s`$ decrease slowly as $`\alpha `$ decreases. It is worthwhile to note that $`K_s`$ is negative in the low $`\alpha `$ side, in opposition to what is obtained in the usual $`\sigma `$-$`\omega `$ model. Experimentally, $`K_s`$ obtained from the isoscalar giant monopole resonance energy is between $`566\pm 1350`$ to $`34\pm 159`$MeV. On the other hand, $`K_0`$ increases to very high values rapidly as $`\alpha `$ increases. In this aspect, the case of $`\eta =i`$ is not acceptable.

The dimensionless effective nucleon mass $`\xi `$ at the saturation point can be used as the third input to fix the parameter $`\alpha `$. For example, if $`\xi =0.85`$ at the saturation point, we have $`\alpha =0.1822`$. Correspondingly we have $`\alpha _\sigma =0.1134`$, $`\alpha _\omega =0.2938`$, $`\xi _\sigma =1.085`$, $`\xi _\omega =0.895`$, $`K_0=501`$MeV, $`J=13.8`$MeV, $`L=30.1`$MeV, and $`K_s=38`$MeV, where the values of $`\xi _\sigma `$ and $`\xi _\omega `$ are those at the saturation point. Instead of $`\xi `$, the nuclear incompressibility $`K_0`$ could be also used to fix $`\alpha `$. If $`K_0=250`$MeV is chosen, it yields $`\alpha =0.1050`$, $`\alpha _\sigma =0.00783`$, $`\alpha _\omega =2.642`$, $`\xi =0.939`$, $`\xi _\sigma =1.295`$, $`\xi _\omega =0.328`$, and $`K_s=31.2`$MeV, where the values of $`\xi `$, $`\xi _\sigma `$ and $`\xi _\omega `$ are those at the saturation point.

The symmetry energy $`J`$ can be increased by including the $`\rho `$ meson field in the present model. In this case, the additional term in the Lagrangian density is

$$_\rho =\frac{1}{4}𝐁_{\mu \nu }\mathbf{}𝐁^{\mu \nu }+\frac{1}{2}m_\rho ^2𝐛_\mu \mathbf{}𝐛^\mu g_\rho \overline{\psi }\gamma _\mu 𝝉\mathbf{}𝐛^\mu \psi ,$$
(36)
where $`𝐁^{\mu \nu }=^\mu 𝐛^\nu ^\nu 𝐛^\mu `$, $`𝐛^\mu `$ is the $`\rho `$ meson field with mass $`m_\rho `$ and coupling constant $`g_\rho `$, and $`𝝉`$ are the isospin matrices. The resulting nucleon field equation, in the mean field approximation for static nuclear matter, has an additional term $`g_\rho \tau _3\gamma ^0b_0`$:

$$(i\gamma ^\mu _\mu g_\omega \gamma ^0\omega _0g_\rho \tau _3\gamma ^0b_0M^{})\psi =0.$$
(37)
The equation for the $`\rho `$ meson field is

$$b_0=\frac{g_\rho \rho _N\delta }{m_\rho ^2},$$
(38)
while Eqs.(3) and (4) are the same. The additional terms to the nuclear energy density, pressure, generalized incompressibility, symmetry and density symmetry energies are respectively

$$_\rho =\frac{1}{2}\alpha _\rho M\rho _N^2\delta ^2/\rho _0,p_\rho =_\rho ,K_\rho =\frac{18_\rho }{\rho _N},J_\rho =\frac{1}{2}\alpha _\rho M,L_\rho =3J_\rho ,$$
(39)
where

$$\alpha _\rho =\frac{g_\rho ^2\rho _0}{m_\rho ^2M}$$
(40)
is the $`\rho `$ meson dimensionless composite parameter. It is easy to prove that there is no contribution to $`K_s`$ from the $`\rho `$ meson field.

It can be seen that, as $`_\rho =0`$ at the standard state where $`\delta =0`$, the inclusion of $`\rho `$ meson produces no change in nuclear energy density $``$, pressure $`p`$ and generalized incompressibility $`K`$ at the standard state. Therefore, the parameters $`\alpha `$, $`\alpha _\sigma `$ and $`\alpha _\omega `$ have no change, and thus the dimensionless effective masses $`\xi `$, $`\xi _\sigma `$ and $`\xi _\omega `$ are the same even the $`\rho `$ meson field is included. The parameter $`\alpha _\rho `$ can be fixed by using the symmetry energy $`J`$ as the fourth input, if parameters $`\alpha `$, $`\alpha _\sigma `$ and $`\alpha _\omega `$ have been fixed.

In conclusion, the possibility of extending the linear $`\sigma `$-$`\omega `$ model by introducing a $`\sigma `$-$`\omega `$ coupling phenomenologically is explored. It is shown that, in contrast to the usual Walecka model, not only the effective nucleon mass $`M^{}`$ but also the effective $`\sigma `$ meson mass $`m_\sigma ^{}`$ and the effective $`\omega `$ meson mass $`m_\omega ^{}`$ are nucleon density dependent. When the model parameters are fitted to the nuclear radius constant $`r_0=1.14`$fm and volume energy $`a_1=16.0`$MeV as well as to the effective nucleon mass $`M^{}=0.85M`$, the model yields $`m_\sigma ^{}=1.09m_\sigma `$ and $`m_\omega ^{}=0.90m_\omega `$ at the same nuclear saturation point, and the nuclear incompressibility $`K_0=501`$MeV. This incompressibility seems too high. On the other hand, if the model parameters are fitted to $`K_0=250`$MeV, in addition to $`r_0=1.14`$fm and $`a_1=16.0`$MeV, it yields $`M^{}=0.938M`$, $`m_\sigma ^{}=1.295m_\sigma `$ and $`m_\omega ^{}=0.328m_\omega `$ at the nuclear saturation point. This effective $`\omega `$ meson mass seems too low. In addition, the lower limit of incompressibility $`K_0=227`$MeV is not low enough. Therefore, even this model is able to reduce the nuclear incompressibility and also to increase the effective nucleon mass simultaneously, there is still not enough degrees of freedom to adjust the parameters in order to give more reasonable result. In this aspect, additional physics should be included for improving this model.


text:
# Detecting X-ray Synchrotron Emission in Supernova Remnants: Implications for Abundances and Cosmic Rays

## 1. Introduction

A growing number of SNRs (listed in Table 1) are thought to show evidence for synchrotron X-rays, for a variety of reasons. Most are suspected of nonthermal X-rays because they show hard energy excesses that can be fit with powerlaws. A few (SN 1006, G347.3-0.5) have nearly featureless X-ray spectra which are clearly dominated by nonthermal emission. Some remnants, while not dominated by nonthermal emission in the ASCA band, show other evidence pointing to its presence. OSSE observations of Cassiopeia A at 400-1250 keV (The 1996) are well described by a broken power law, steepening to higher energies (Allen 1997). RCW 86 is an especially interesting case. It shows anomalously low abundances when fit with thermal models (Vink 1997), but the addition of a synchrotron model bring the abundances to within expected ranges (Borkowski 2001).

The most straightforward explanation for nonthermal X-ray emission (non-plerionic in these SNRs) is synchrotron emission from the same mechanism that is responsible for the radio emission. It is straightforward to determine the synchrotron spectrum from the electron energy population However, the extrapolation from the radio spectrum measured by single dish instruments, over-predicts the X-ray flux (see Figure 1) – there must be a mechanism that causes the electron population to roll off at the highest energies. Reynolds & Chevalier (1981) attributed this rolloff to synchrotron losses.

The sharpest physically plausible spectral rolloff would be the emission resulting from an exponential cutoff of the electron power-law energy distribution, radiating in a uniform magnetic field. This was embodied in the XSPEC model SRCUT and used by Reynolds & Keohane (1999) to set limits on the energies to which SNRs could accelerate electrons (with implications for cosmic ray ions). A more sophisticated treatment requires addressing the mechanism limiting electrons at the highest energies. The cutoff could be due to the age of the SNR, radiative losses from synchrotron or inverse-Compton emission, or electron escape(Reynolds 1998).

## 2. The Escape-limited Synchrotron Model

The escape-limited model is a single parameter model that is particularly well suited to X-ray analysis programs such as XSPEC, and has been released as SRESC in XSPEC 11 (it should also be available for Chandra data analysis in CIAO 2.1). This model uses the Sedov solution to specify the electron distribution and assumes a uniform upstream magnetic field, with a shock compression ratio of 4. The model takes the radio spectral index and 1 GHz flux as inputs (available from Green 2000 for galactic SNRs and now for some Magellanic Cloud SNRs in Williams 2000) and fits one parameter, the rolloff frequency – proportional to the product of the maximum MHD wavelength scattering the electrons and the magnetic field strength. We have applied SRESC to archival observations of the integrated ASCA spectrum of SN 1006 as a test case. As well as describing the physics behind the emission the SRESC model is much better constrained for fitting purposes than a powerlaw model. The synchrotron model allows more accurate abundance determinations in the thermal emission – these abundances are close to predictions for Type Ia SN (Iwamoto et al. 1999) and include half a solar mass of iron. Since TeV detections can be used to determine the relativistic electron density in SN 1006, fitting the rolloff frequency with the SRESC model uniquely determines the magnetic field. This suggests that SN 1006 is in fact far from equipartition, violating a common assumption for astrophysical plasmas. The model also makes spatial predictions and is capable of reproducing the morphology of SN 1006 at energies at which it is dominated by nonthermal emission, shown in Figure 2.

However, a much more stringent test of the model is to apply it to spatially resolved spectra from SN 1006. As shown in Figure 2 we have developed subsets of the escape-limited synchrotron model for different regions of the supernova remnant, describing the polar caps, the limbs and the center. These models predict slightly different curvature for the nonthermal emission in each region for a particular rolloff energy. The models should allow us to make more accurate determinations of the abundances as well as mapping the thermal and nonthermal distribution across the remnant.

Preliminary results have been very promising. Contrary to expectation, all regions of SN 1006 have a significant amount of nonthermal emission, including the polar caps and center which have been fit in the past with thermal-only models. The limbs are dominated by nonthermal emission but require a thermal component while fits to the thermally dominated center will tolerate a significant amount of nonthermal emission. In addition we find that the north limb and south limb will not support the identical thermal and nonthermal models – we find the spectra to be detectably different.

## 3. Conclusions

We have found, both in integrated and spatially resolved models, that the measured abundances and temperature change with the addition of accurate synchrotron models. In the future, accurate accounting for a possible nonthermal component will be a necessary precondition to using nucleosynthesis model predictions to interpret abundances in any SNR. SRCUT and SRESC, new models in XSPEC 11, accurately describe this emission and are a substantial improvement over powerlaw models both in terms of physical description and fitting.

All regions of SN 1006, even the polar caps which show obvious thermal lines, seem to have significant nonthermal emission. We also find that the spectra of the north and the south limb are not identical. While most of the emission from the center does seem to be thermal, the fit can tolerate a significant amount of nonthermal emission.

Technology is surpassing theory – Chandra arcsecond-resolution observations of SNRs including SN 1006 will soon be available, requiring models capable of describing extremely small scale variations. These synchrotron models will be useful both for analyzing higher resolution observations of SN 1006, and for interpreting other SNRs with nonthermal X-ray emission.

### Acknowledgments.

This work is supported by NASA’s Graduate Student Research Program through GSFC. http://education.nasa.gov/gsrp/

## References

Allen, G. E. et al. 1997, ApJ, 487, L97

Bamba, A., Koyama, K., Tomida, H. 2001 Accepted for publication in PAS astro-ph/0007365

Borkowski, K.J., Rho, J., Reynolds, S.P., Dyer, K.K. 2001 Accepted for the Astrophysical Journal March 10, 2001

Decourchelle, A. & Petre, R. 1999, Astronomische Nachrichten, 320, 203

Dyer, K.K., Reynolds, S.P., Borkowski, K.J., Petre, R., Allen, G.E.2001 astro-ph/0010424 Accepted for the Astrophysical Journal 549 March 10, 2001

Favata, F. et al. 1997, A&A, 324, L49

Green D.A., 2000, http://www.mrao.cam.ac.uk/surveys/snrs/

Koyama, K., Kinugasa, K., Matsuzaki, K., Nishiuchi, M., Sugizaki, M., Torii, K., Yamauchi, S. & Aschenbach, B. 1997, PASJ, 49, L7

Koyama, K., Petre, R., Gotthelf, E. V., Hwang, U., Matsura, M., Ozaki, M. & Holt, S. S. 1995, Nature, 378, 255

Iwamoto, K., Brachwitz, F., Nomoto, K., Kishimoto, N., Umeda, H., Hix, W. R. & Thielemann, F. 1999, ApJS, 125, 439

Reynolds, S.P. 1998 ApJ 493, 357

Reynolds, S.P. & Chevalier, R. A. 1981, ApJ, 245, 912

Reynolds, S. P. & Keohane, J. W. 1999, ApJ, 525, 368

Slane, P., Gaensler, B.M., Dame, T.M., Hughes, J.P., Plucinsky, P.P. & Green, A. 1999, ApJ, 525, 357

Sun, M., Wang, Z. & Chen, Y. 1999, ApJ, 511, 274

The, L. -., Leising, M. D., Kurfess, J. D., Johnson, W. N., Hartmann, D. H., Gehrels, N., Grove, J. E. & Purcell, W. R. 1996, A&AS, 120, C357

Vink, J., Kaastra, J. S. & Bleeker, J. A. M. 1997, A&A, 328, 628

Williams, R. 2000 Magellanic Catalog http://www.renaissoft.com/snrc/


text:
# 1 Introduction

## 1 Introduction

The Standard Model allows $`B^0\overline{B^0}`$ oscillations to occur via second order weak interactions. The frequency of oscillation is determined by the mass differences, $`\mathrm{\Delta }m`$, between the mass eigenstates in the $`B^0`$ system. The mass difference in the $`B_s^0`$ system ($`\mathrm{\Delta }m_s`$) and in the $`B_d^0`$ system ($`\mathrm{\Delta }m_d`$) are proportional to the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements $`\left|V_{ts}\right|^2`$ and $`\left|V_{td}\right|^2`$, respectively. A measurement of $`\mathrm{\Delta }m_d`$ can in principle be used to extract the CKM matrix element $`\left|V_{td}\right|`$. However, the extraction of $`\left|V_{td}\right|`$ from $`\mathrm{\Delta }m_d`$ is complicated by a large theoretical uncertainty on the hadronic matrix elements. The complication can be circumvented by taking the ratio of $`\mathrm{\Delta }m_s`$ and $`\mathrm{\Delta }m_d`$. In the ratio, the theoretical uncertainty is reduced to the 5% level . Therefore, a direct measurement of $`\mathrm{\Delta }m_s`$, combined with the current measurement of $`\mathrm{\Delta }m_d`$, can be translated to a precise value of $`\left|V_{td}\right|`$.

## 2 Experimental Technique

The results presented in this paper are based on 400,000 hadronic $`Z^0`$ decays collected during the 1996-98 run with an average electron beam polarization of 73%. A detailed description of the experimental apparatus can be found elsewhere . The three main ingredients for measuring the time dependent $`B_s^0\overline{B_s^0}`$ oscillations are: (1) determination of the flavor at production (initial state tag), (2) determination of the flavor at the decay vertex (final state tag), and (3) reconstruction of the proper decay time of the $`B_s^0`$. Three methods were explored at SLD for studying the $`B_s^0`$ oscillations and the methods are referred to as “$`D_s`$+Tracks”, “Lepton+D”, and “Charge Dipole” . All three analyses share the common initial state tag as well as the B energy reconstruction algorithms but differ in the event selection, decay length reconstruction, and the final state tag.

Several techniques are used to determine the initial state of the $`B_s^0`$. The most powerful method, unique to the SLD, is the polarization tag. In a polarized $`Z^0b\overline{b}`$ decay, the outgoing quark is produced preferentially along the direction opposite to the spin of the $`Z^0`$ boson. Therefore by knowing the helicity of the electron beam and the direction of the jet, the flavor of the primary quark in the jet can be determined. To further enhance the initial state tag, information in the opposite hemisphere (e.g. momentum weighted jet charge, vertex charge, lepton and kaon tracks) is used. Combining all available tags, the average initial state correct tag probability is about 78%.

An ideal mixing analysis requires high efficiency, high $`B_s^0`$ purity, clean initial and final state tags, and excellent proper time resolution. In practice, experimental constraints necessitate trade-offs between the four key elements. The “$`D_s`$+tracks” is the most exclusive analysis at SLD. The analysis partially reconstructs the $`B_s^0`$ by combining a fully reconstructed $`D_s`$ (via $`\varphi \pi `$ and $`K^0K`$ modes) with other secondary B decay tracks. By taking the exclusive approach, the analysis is able to achieve a high average $`B_s^0`$ purity of 38% and an excellent decay length resolution of 48$`\mu m`$ (60% core resolution) for the $`B_s^0`$ events. However, the analysis suffers from low efficiency and only 361 candidates are selected in the final sample. The “Lepton+D” takes a slightly more inclusive approach by selecting events with a partially reconstructed D vertex and a lepton in the same hemisphere. The identified lepton not only enhances the b-hadron fraction but also provides a clean final state tag (final state mistag $`<`$ 10%). The estimated B energy resolution is comparable to the other two analyses and is about 7% (60% core fraction) and 20% (tail) for the $`B_s^0(bl)`$ events. The “Lepton+D” has 2087 candidates in the final sample with an average $`B_s^0`$ purity of about 16%. The most inclusive method and also the analysis with the highest sensitivity at SLD is the “Charge Dipole”. The “Charge Dipole” selects events that contain both a secondary and a tertiary vertex. To enhance the $`B_s^0`$ fraction, the total track charge (from secondary and tertiary) is required to be zero. The final state tag is based on the dipole value which is defined as the sign of the charge difference weighted by the distance between the secondary and the tertiary vertices. The final state correct tag probability using the dipole method is highly dependent on the decay toplogy and ranges from 53% for $`B_s^0D\overline{D}X`$ to 91% for $`B_s^0D_sX`$ decays. A total of 8556 decays is selected in the “Charge Dipole” analysis with a $`B_s^0`$ purity of 15%.

## 3 Results

The study of the $`B_s^0`$ oscillations is performed using the amplitude fit method . In this method, the probability for mixing, which is proportional to $`1cos(\mathrm{\Delta }m_s\tau )`$, is modified by introducing the amplitude A in front of the cosine (same modification for the unmixed expression). The amplitude plot is generated by scanning the $`\mathrm{\Delta }m_s`$ value over a specified range, and for each $`\mathrm{\Delta }m_s`$ value, fitting for the parameter A. We expect the fitted value of A to be consistent with zero when the chosen $`\mathrm{\Delta }m_s`$ is away from the true value and the amplitude A to reach the value 1 near the true $`\mathrm{\Delta }m_s`$. If no signal is seen, a 95$`\%`$ C.L. lower limit can be set for frequencies at which $`𝐀+1.645\sigma _A<1`$. The 95$`\%`$ C.L. sensitivity is defined as the value of $`\mathrm{\Delta }m_s`$ at which $`1.645\sigma _A=1`$.

The amplitude plot for the three SLD analyses combined is shown in Fig. 1. The combined plot takes into account the correlated systematic uncertainties. Furthermore, the samples were selected such as to remove any statistical overlap between analyses. No evidence of a signal is observed up to $`\mathrm{\Delta }m_s`$ of 25 $`ps^1`$. The preliminary SLD results exclude the following values of $`\mathrm{\Delta }m_s`$ at the 95% C.L.: $`\mathrm{\Delta }m_s`$ $`<`$ 7.6 $`ps^1`$ and 11.8 $`<`$ $`\mathrm{\Delta }m_s`$ $`<`$ 14.8 $`ps^1`$. The SLD combined sensitivity at the 95% C.L. is 13.0 $`ps^1`$.

## Acknowledgments

We thank the personnel of the SLAC accelerator department and the technical staffs of our collaborating institutions for their outstanding efforts. This work was supported by the Department of Energy, the National Science Foundataion, the Instituto Nazionale di Fisica of Italy, the Japan-US Cooperative Research Project on High Energy Physics, and the Science and Engineering Research Council of the United Kingdom.

## ∗∗ List of Authors

Kenji Abe,(15) Koya Abe,(24) T. Abe,(21) I. Adam,(21) H. Akimoto,(21) D. Aston,(21) K.G. Baird,(11) C. Baltay,(30) H.R. Band,(29) T.L. Barklow,(21) J.M. Bauer,(12) G. Bellodi,(17) R. Berger,(21) G. Blaylock,(11) J.R. Bogart,(21) G.R. Bower,(21) J.E. Brau,(16) M. Breidenbach,(21) W.M. Bugg,(23) D. Burke,(21) T.H. Burnett,(28) P.N. Burrows,(17) A. Calcaterra,(8) R. Cassell,(21) A. Chou,(21) H.O. Cohn,(23) J.A. Coller,(4) M.R. Convery,(21) V. Cook,(28) R.F. Cowan,(13) G. Crawford,(21) C.J.S. Damerell,(19) M. Daoudi,(21) S. Dasu,(29) N. de Groot,(2) R. de Sangro,(8) D.N. Dong,(21) M. Doser,(21) R. Dubois,(21) I. Erofeeva,(14) V. Eschenburg,(12) S. Fahey,(5) D. Falciai,(8) J.P. Fernandez,(26) K. Flood,(11) R. Frey,(16) E.L. Hart,(23) K. Hasuko,(24) S.S. Hertzbach,(11) M.E. Huffer,(21) X. Huynh,(21) M. Iwasaki,(16) D.J. Jackson,(19) P. Jacques,(20) J.A. Jaros,(21) Z.Y. Jiang,(21) A.S. Johnson,(21) J.R. Johnson,(29) R. Kajikawa,(15) M. Kalelkar,(20) H.J. Kang,(20) R.R. Kofler,(11) R.S. Kroeger,(12) M. Langston,(16) D.W.G. Leith,(21) V. Lia,(13) C. Lin,(11) G. Mancinelli,(20) S. Manly,(30) G. Mantovani,(18) T.W. Markiewicz,(21) T. Maruyama,(21) A.K. McKemey,(3) R. Messner,(21) K.C. Moffeit,(21) T.B. Moore,(30) M. Morii,(21) D. Muller,(21) V. Murzin,(14) S. Narita,(24) U. Nauenberg,(5) H. Neal,(30) G. Nesom,(17) N. Oishi,(15) D. Onoprienko,(23) L.S. Osborne,(13) R.S. Panvini,(27) C.H. Park,(22) I. Peruzzi,(8) M. Piccolo,(8) L. Piemontese,(7) R.J. Plano,(20) R. Prepost,(29) C.Y. Prescott,(21) B.N. Ratcliff,(21) J. Reidy,(12) P.L. Reinertsen,(26) L.S. Rochester,(21) P.C. Rowson,(21) J.J. Russell,(21) O.H. Saxton,(21) T. Schalk,(26) B.A. Schumm,(26) J. Schwiening,(21) V.V. Serbo,(21) G. Shapiro,(10) N.B. Sinev,(16) J.A. Snyder,(30) H. Staengle,(6) A. Stahl,(21) P. Stamer,(20) H. Steiner,(10) D. Su,(21) F. Suekane,(24) A. Sugiyama,(15) S. Suzuki,(15) M. Swartz,(9) F.E. Taylor,(13) J. Thom,(21) E. Torrence,(13) T. Usher,(21) J. Va’vra,(21) R. Verdier,(13) D.L. Wagner,(5) A.P. Waite,(21) S. Walston,(16) A.W. Weidemann,(23) E.R. Weiss,(28) J.S. Whitaker,(4) S.H. Williams,(21) S. Willocq,(11) R.J. Wilson,(6) W.J. Wisniewski,(21) J.L. Wittlin,(11) M. Woods,(21) T.R. Wright,(29) R.K. Yamamoto,(13) J. Yashima,(24) S.J. Yellin,(25) C.C. Young,(21) H. Yuta.(1)

(The SLD Collaboration)

(1)Aomori University, Aomori , 030 Japan, (2)University of Bristol, Bristol, United Kingdom, (3)Brunel University, Uxbridge, Middlesex, UB8 3PH United Kingdom, (4)Boston University, Boston, Massachusetts 02215, (5)University of Colorado, Boulder, Colorado 80309, (6)Colorado State University, Ft. Collins, Colorado 80523, (7)INFN Sezione di Ferrara and Universita di Ferrara, I-44100 Ferrara, Italy, (8)INFN Lab. Nazionali di Frascati, I-00044 Frascati, Italy, (9)Johns Hopkins University, Baltimore, Maryland 21218-2686, (10)Lawrence Berkeley Laboratory, University of California, Berkeley, California 94720, (11)University of Massachusetts, Amherst, Massachusetts 01003, (12)University of Mississippi, University, Mississippi 38677, (13)Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, (14)Institute of Nuclear Physics, Moscow State University, 119899, Moscow Russia, (15)Nagoya University, Chikusa-ku, Nagoya, 464 Japan, (16)University of Oregon, Eugene, Oregon 97403, (17)Oxford University, Oxford, OX1 3RH, United Kingdom, (18)INFN Sezione di Perugia and Universita di Perugia, I-06100 Perugia, Italy, (19)Rutherford Appleton Laboratory, Chilton, Didcot, Oxon OX11 0QX United Kingdom, (20)Rutgers University, Piscataway, New Jersey 08855, (21)Stanford Linear Accelerator Center, Stanford University, Stanford, California 94309, (22)Soongsil University, Seoul, Korea 156-743, (23)University of Tennessee, Knoxville, Tennessee 37996, (24)Tohoku University, Sendai 980, Japan, (25)University of California at Santa Barbara, Santa Barbara, California 93106, (26)University of California at Santa Cruz, Santa Cruz, California 95064, (27)Vanderbilt University, Nashville,Tennessee 37235, (28)University of Washington, Seattle, Washington 98105, (29)University of Wisconsin, Madison,Wisconsin 53706, (30)Yale University, New Haven, Connecticut 06511.


text:
# Cosmological Evolution of Linear Bias

## Abstract

Using linear perturbation theory and the Friedmann-Lemaitre solutions of the cosmological field equations, we derive analytically a second-order differential equation for the evolution of the linear bias factor, $`b(z)`$, between the background matter and a mass-tracer fluctuation field. We find $`b(z)`$ to be a strongly dependent function of redshift in all cosmological models. Comparing our analytical solution with the semi-analytic model of Mo & White, which utilises the Press-Schechter formalism and the gravitationally induced evolution of clustering, we find an extremely good agreement even at large redshifts, once we normalize to the same bias value at two different epochs, one of which is the present. Furthermore, our analytic $`b(z)`$ function agrees well with the outcome of N-body simulations.

Keywords Cosmology: theory - large-scale structure of universe

## 1 Introduction

The concept of biasing between different classes of extragalactic objects and the background matter distribution was put forward by Kaiser (1984) and Bardeen et al. (1986) in order to explain the higher amplitude of the 2-point correlation function of clusters of galaxies with respect to that of galaxies themselves.

In this framework biasing is assumed to be statistical in nature; galaxies and clusters are identified as high peaks of an underlying initially Gaussian random density field. Biasing of galaxies with respect to the dark matter distribution was also found to be an essential ingredient of CDM models of galaxy formation in order to reproduce the observed galaxy distribution (cf. Davies et al. 1985; Benson et al. 2000).

The classical approach to study the redshift evolution of bias utilises the ratio of the correlation functions of objects and dark matter, which are assumed to be related via the square of a scale independent bias factor. However, in this study we will use the definition by which the extragalactic mass tracer (galaxies, halos, clusters) fluctuation field, $`\delta _{tr}`$, is related to that of the underlying mass, $`\delta _\mathrm{m}`$, by

$$\delta _{\mathrm{tr}}=b\delta _\mathrm{m},$$
(1)
where $`b`$ is the linear bias factor. Note that the former definition results from the latter but the opposite is not necessarily true. The bias factor may have many dependencies; even assuming that it is scale independent, it necessarily depends on the type of the mass tracer as well as on the epoch $`z`$, since the fluctuations evolve with time as gravity draws together galaxies and mass. It is evident, therefore, that the bias factor should also depend on the different cosmological models and dark matter content of the Universe (for a recent overview see Klypin 2000).

More realistic biasing schemes have been proposed in the literature. Coles (1993) introduced the idea of biased galaxy formation in which galaxies form with a probability given by an arbitrary function of the local mass density. Mann, Peacock & Heavens (1998) investigated the properties of different bias models of galaxy distributions that results from local transformations of the present-day density field. The deterministic and linear nature of eq.(1) has been challenged (cf. Dekel & Lahav 1999; Tegmark & Bromley 1999) and indeed some non-linearity of the biasing relation is necessary to reconcile high biasing with deep voids. Despite the above, the linear biasing assumption is still a useful first order approximation which, due to its simplicity, it is used in most studies of large scale (linear) dynamics (cf. Strauss & Willick 1995 and references therein; Branchini et al. 1999; Schmoldt et al. 1999; Plionis et al. 2000).

Different studies have indeed shown that the bias factor is a monotonically decreasing function of redshift. An important advancement in the analytical treatment of the bias evolution was the work of Mo & White (1996) in which they used the Press-Schechter (1974) formalism and found that in an Einstein-de Sitter universe the linear bias factor evolves strongly with redshift. Using a similar formalism, Matarrese et al. (1997) extended the Mo & White results to include the effects of different mass scales (see also Catelan et al 1998).

Steidel et al. (1998) confirmed that the Lyman-break galaxies are very strongly biased tracers of mass and they found that $`b(z=3.4)>\text{ }6,4,2`$, for SCDM, $`\mathrm{\Lambda }`$CDM $`(\mathrm{\Omega }=0.3)`$ and OCDM $`(\mathrm{\Omega }=0.2)`$, respectively (see also Giavalisko et al 1998). A similar value for the $`\mathrm{\Lambda }`$CDM model was obtained by Cen & Ostriker (2000) using high resolution Nbody/hydro simulations in which they treated DM, gas as well as star formation. The use of high resolution N-body simulations (cf. Klypin et al. 1996; 1999, Cole et al. 1997 and references therein) have shown that anti-biasing ($`b<1`$) should exist at scales $`r38h^1`$ Mpc, for the open and flat low-$`\mathrm{\Omega }`$ models, in contrast with $`\mathrm{\Omega }=1`$ models, where $`b>1`$. Colin et al (1999), using high-resolution N-body simulations of SCDM, $`\mathrm{\Lambda }`$CDM, OCDM and $`\tau `$CDM models, which avoid the so called “overmerging” problem, found that indeed biasing evolves rapidly with redshift, while Kauffmann et al. (1999) combining semi-analytic models of galaxy formation and N-body simulations has also studied the evolution of clustering in different cosmologies.

In this paper we will not indulge in such aspects of the problem but rather, working within the paradigm of linear and scale-independent bias, we will derive the functional form of its redshift evolution in the matter dominated epoch and in all cosmological models. The Einstein de-Sitter case has been studied in the past (cf. Nusser & Davis 1994; Fry 1996; Bagla 1998) using the continuity equation, which is a first order differential equation, to derive a solution, $`(1+z)`$, valid only for low $`z`$’s. Our approach is to use the perturbation evolution equation which combines the continuity, the Euler and the Poisson equations and which is a second order differential equation. We should therefore expect to find a further component to the known solution.

The paper is organised as follows: in section 2 we discuss the basic models for the linear bias evolution, in section 3 we derive the basic differential equation describing the evolution of the linear bias factor, while in section 4 we present its analytical solution for the different cosmological models and a comparison with previous models and N-body simulation results. Finally, in section 5 we summarise our main results.

## 2 Models for bias evolution

Theoretical expectations regarding the cosmological evolution of bias have been investigated using analytical calculations, semi-analytical approximations and N-body simulations. In this section we shortly describe some of these models in order to compare them with our results.

### 2.1 Test Particle or Galaxy Conserving Bias (M1):

This model, proposed by Nusser & Davis (1994), Fry (1996), Tegmark & Peebles (1998), predicts the evolution of bias, independent of the mass and the origin of halos, assuming only that the test particles fluctuation field is related proportionally to that of the underlying mass. Thus, the bias factor as a function of redshift can be written:

$$b(z)=1+(b_{}1)D(z)^1$$
(2)
where $`b_{}`$ is the bias factor at the present time. Bagla (1998) found that for SCDM model and in the range $`0z1`$ the above formula describes well the evolution of bias.

### 2.2 Halo Models (M2):

Mo & White (1996) using the Press-Schecter formalism, have developed a model for the evolution of the correlation bias, which depends on halo mass, and found, in an Einstein-de Sitter Universe, that:

$$b(z)=1+\frac{1}{\delta _c}\left[\left(\frac{D(z_{})}{D(z)}\right)^21\right],$$
(3)
with $`z_{}`$ a reference redshift 111Consider the distribution of halos of mass $`M`$, or larger, before typical halos of this mass have collapsed. Therefore, to quantify this, the parameter $`\nu (M,z)=\delta _c/\sigma (M,z)`$, which related directly with the bias, $`b=1+\frac{1}{\delta _c}(\nu ^21)`$, is usually utilized (where $`\delta _c=1.69`$ is the critical overdensity for spherical collapse at $`z=0`$, and $`\sigma (M,z)`$ is the rms linear mass fluctuation on the scale $`M`$ of halos linearly extrapolated to redshift $`z`$) . Thus, $`z_{}`$ is fixed by requiring $`\nu (M,z_{})=1`$ (cf. Bagla 1998)., $`\delta _c`$ the critical overdensity for a spherical top-hat collapse model and $`D(z)=(1+z)^1`$ is the linear growth rate of clustering. Parametrising this equation to the present epoch one gets:

$$b(z)=0.41+(b(0)0.41)D(z)^2$$
(4)
Similarly, Matarrese et al. (1997) parametrising the evolution of bias for halos above a certain mass $`M`$, obtain a similar expression for an Einstein-de Sitter Universe:

$$b(z)=0.41+(b_{\mathrm{eff}}(0)0.41)D(z)^\beta $$
(5)
with $`b_{\mathrm{eff}}`$ the bias of a sample of halos with a range of masses and $`\beta `$ depending on the minimum mass scale that contributes to the halo correlation function (with $`\beta <2`$).

## 3 Basic Equations

The central issue here is to derive the basic differential equation which describes the evolution of bias. The present analysis is based on linear perturbation theory in the matter dominated epoch (cf. Peebles 1993) and it is an extension of the M1 model.

The time evolution equation for the mass density contrast, $`\delta _\mathrm{m}=(\delta \rho /\rho )_\mathrm{m}`$, modelled as a pressureless fluid with general solution of the growing mode: $`\delta _\mathrm{m}A(x)D(t)`$, is (cf. Padmanabhan 1993):

$$\frac{^2\delta _\mathrm{m}}{t^2}+2H(t)\frac{\delta _\mathrm{m}}{t}=4\pi G\rho _\mathrm{m}\delta _\mathrm{m},$$
(6)
Assuming for simplicity that the mass tracer population is conserved in time, ie., that the effects of non-linear gravity and hydrodynamics (merging, feedback mechanisms etc) do not significantly alter the population mean, then a similar evolution equation, containing in the right hand side the gravitational contributions of all the perturbed matter, should be satisfied for $`\delta _{tr}`$ (see also Fry 1996; Catelan et al. 1998):

$$\frac{^2\delta _{tr}}{t^2}+2H(t)\frac{\delta _{tr}}{t}=4\pi G\rho _\mathrm{m}\delta _\mathrm{m}.$$
(7)
Differentiating twice eq.(1) and using eq.(7) and eq.(6) we obtain:

$$\ddot{b}\delta _\mathrm{m}+2\left[\dot{\delta }_\mathrm{m}+H(t)\delta _\mathrm{m}\right]\dot{b}+\left[\ddot{\delta }_\mathrm{m}+2H(t)\dot{\delta }_\mathrm{m}\right]b=4\pi G\rho _\mathrm{m}\delta _\mathrm{m}$$
(8)
Then from equation eq.(8), eq.(6) and $`\delta _\mathrm{m}A(x)D(t)`$ we have:

$$\ddot{b}D(t)+2\left[\dot{D(t)}+H(t)D(t)\right]\dot{b}+4\pi G\rho _\mathrm{m}D(t)b=4\pi G\rho _\mathrm{m}D(t).$$
(9)
In order to transform eq.(9) from time to redshift we use the following expression:

$$\frac{dt}{dz}=\frac{1}{H(z)(1+z)},$$
(10)
where the Hubble parameter is given by:

$$H(z)=H_{}E(z)$$
(11)
with

$$E(z)=\left[\mathrm{\Omega }(1+z)^3+\mathrm{\Omega }_R(1+z)^2+\mathrm{\Omega }_\mathrm{\Lambda }\right]^{1/2}$$
(12)
and $`\mathrm{\Omega }=8\pi G\rho _o/3H_o^2`$ (density parameter), $`\mathrm{\Omega }_R=(H_o\alpha _oR)^2`$ (curvature parameter), $`\mathrm{\Omega }_\mathrm{\Lambda }=\mathrm{\Lambda }/3H_o^2`$ (cosmological constant parameter) at the present time which satisfy $`\mathrm{\Omega }+\mathrm{\Omega }_R+\mathrm{\Omega }_\mathrm{\Lambda }=1`$ and $`H_{}`$ is the Hubble constant.

Finally, the growing solution (cf. Peebles 1993) as a function of redshift is:

$$D(z)=\frac{5\mathrm{\Omega }E(z)}{2}_z^{\mathrm{}}\frac{(1+x)}{E^3(x)}dx.$$
(13)
Therefore, as the time evolves with redshift, utilised eq.(10), eq.(11) eq.(12) and the relation

$$4\pi G\rho _\mathrm{m}=4\pi G\rho _{}(1+z)^3=\frac{3H_{}^2}{2}\mathrm{\Omega }(1+z)^3,$$
(14)
then the basic differential equation for the evolution of the linear bias parameter takes the following form:

$$\frac{\mathrm{d}^2b}{\mathrm{d}z^2}P(z)\frac{\mathrm{d}b}{\mathrm{d}z}+Q(z)b=Q(z)$$
(15)
with basic factors,

$$P(z)=\frac{1}{1+z}\frac{1}{E(z)}\frac{\mathrm{d}E(z)}{\mathrm{d}z}\frac{2}{D(z)}\frac{\mathrm{d}D(z)}{\mathrm{d}z}$$
(16)
and

$$Q(z)=\frac{3\mathrm{\Omega }(1+z)}{2E^2(z)}.$$
(17)
It is obvious that the above generic form depends on the choice of the background cosmology. Thus, the functional form which satisfies the general bias solution for all of the cosmological models is:

$$b(z)=y(z;\mathrm{\Omega };\mathrm{\Omega }_\mathrm{\Lambda })+1$$
(18)
where $`y`$ is the general solution of the homogeneous differential equation:

$$\frac{\mathrm{d}^2y}{\mathrm{d}z^2}P(z)\frac{\mathrm{d}y}{\mathrm{d}z}+Q(z)y=0.$$
(19)
Whereas it is obvious that the present theoretical approach takes into account the gravity field, it does not interact directly with the nature of the DM particles.

## 4 Bias Evolution in different Cosmological Models

In this section using both eq.(19) and Friedmann-Lemaitre solutions of the cosmological field equations we present the analytical solution of bias evolution for the Einstein-de Sitter, the open and low-density flat cosmological models.

### 4.1 Elements of the Differential Equation Theory

Without wanting to appear too pedagogical, we remind the reader some basic elements of differential equation theory (cf. Bronson 1973). If one is able to find any solution $`y_1`$ of eq.(19), then a second linearly independent solution be found very easily. Let the second solution can be written as

$$y_2(z)=y_1(z)u(z)$$
(20)
where $`u(z)`$ is to be determined. Inserting eq.(20) into eq.(19) and remembering that $`y_1`$ satisfies the same equation, we find the following equation for $`u(z)`$:

$$\frac{\mathrm{d}^2u}{\mathrm{d}z^2}+\left[\frac{2(\mathrm{d}y_1/\mathrm{d}z)P(z)y_1(z)}{y_1}\right]\frac{\mathrm{d}u}{\mathrm{d}z}=0.$$
(21)
Integrating the above equation we have

$$\frac{\mathrm{d}u}{\mathrm{d}z}=\frac{\mathrm{const}}{y_1^2(z)}\mathrm{exp}\left(_{z_0}^zP(x)dx\right)$$
(22)
where $`z_0`$ is an arbitrary initial point. A further integration of eq.(22) yields $`u(z)`$, and inserting this value into eq.(20), we obtain the second solution

$$y_2(z)=\mathrm{const}y_1(z)_{z_0}^z\frac{\mathrm{d}x}{y_1^2(x)}\mathrm{exp}\left(_{z_0}^xP(t)dt\right).$$
(23)
The Wronskian of the two solution $`y_1`$ and $`y_2`$ is

$$W(y_1,y_2)=y_1^2\frac{\mathrm{d}u}{\mathrm{d}z}=\mathrm{const}\mathrm{exp}\left(_{z_0}^zP(x)dx\right).$$
(24)
Thus the Wronskian never vanishes which implies that any general solution of eq.(19) is a linear combination ($`y=c_1y_1+c_2y_2`$) of the fundamental set of solutions $`y_1`$ and $`y_2`$.

### 4.2 Einstein - de Sitter Model

In this case the basic cosmological equations are the following:

$$E(z)=(1+z)^{3/2},$$
(25)
while the growth factor of the linear density contrast is

$$D(z)=\frac{1}{(1+z)}$$
(26)
and thus

$$P(z)=\frac{3}{2(1+z)}$$
(27)
and

$$Q(z)=\frac{3}{2(1+z)^2}.$$
(28)
It can be found that the function $`y_1(z)=(1+z)=D^1(z)`$ is a solution of the eq.(19) which is to be expected from the M1 model. Therefore, we are looking for the second independent solution of the eq.(19). Thus according to the procedure, described before we can calculate the second solution directly from eq.(20) $`y_2(z)=(1+z)^{3/2}=D^{3/2}(z)`$. The general solution of the second order differential eq.(19) is the following:

$$y(z;1;0)=𝒜D^1(z)+D^{3/2}(z)=𝒜(1+z)+(1+z)^{3/2}$$
(29)
with general bias solution $`b(z)=y(z;1;0)+1`$. To this end this analysis generalise the M1 model in the sense that the added function $`y_2(z)`$ dominates the functional form of the bias evolution. Of course in order to obtain partial solutions for $`b(z)`$ we need to estimate the values of the constants $`𝒜`$ and $``$, which means that we need to calibrate the $`b(z)`$ relation using two different epochs: $`b(z=0)=b_{}`$ and $`b(z=z_1)=b_1`$. Therefore, utilised both the above general bias solution and the latter parameters, we can give the expressions for the above constants as a function of $`b_{}`$ and $`b_1`$:

$$𝒜=\frac{(b_{}1)D^{3/2}(z_1)(b_11)}{D^{3/2}(z_1)D^1(z_1)},$$
(30)
$$=\frac{(b_11)(b_{}1)D^1(z_1)}{D^{3/2}(z_1)D^1(z_1)}.$$
(31)
For $`=0`$ (M1 model) we obtain, as we should, $`𝒜=b_{}1`$.

Our generalised solution does not suffer from limitations in the value of $`b`$ (as does the M1 solution); $`b`$ can take values $`>1`$ and $`<1`$. It is interesting to compare our generalised test-particle bias with the more elaborate halo and merging models. Since our approach gives a family of bias curves, due to the fact that it has two unknown parameters, (the integration constants $`𝒜,`$), we evaluate the latter by using Steidel et al. (1998) value of the bias for Lyman break galaxies which gives for $`\mathrm{\Omega }=1`$, $`b(3.4)7`$. Inserting this into Mo & White (1996) model we obtain $`b(0)=0.75`$. In figure 1 we compare our solution with the Mo & White model and to our surprise we find an excellent agreement. This implies that the complete test particle bias solution is an extremely good approximation to the more elaborate halo solutions which takes into account, via the Press-Schechter formalism, the collapse of different mass halos at the different epochs.

We further compare our analytic solution with N-body estimates provided by Colin et al. (1999) and Kauffmann et al. (1999). In figure 2 our model is represented by lines while the numerical results with the different symbols. It is evident that our analytic function, normalized to two different epochs of the numerical results, fits extremely well the behaviour of the N-body derived bias evolution.

### 4.3 Low Density Universes

The basic cosmological equations in low-density Universes become more complicated than in an $`\mathrm{\Omega }=1`$ Universe and eq.(19) does not have simple analytical solutions. We therefore present approximate analytical solutions which are valid in the high-redshift regime. In order to do so we consider that (i) for a low density open Universe, the Einstein-de Sitter growing mode is a good approximation for $`z>\text{ }\mathrm{\Omega }^11`$ and (ii) for a low density flat Universe the growing mode is well approximated by the Einstein de-Sitter case for $`z>\text{ }\mathrm{\Omega }^{1/3}1`$ (cf. Peebles 1984b; Carrol, Press & Turner 1992).

#### 4.3.1 Analytical Approximation for the Open Universe

Using eqs (16), (17) and $`z>\text{ }\mathrm{\Omega }^11`$ we obtain the following basic factors of the differential equation (19):

$$P(z)=\frac{2}{1+z}\frac{\mathrm{\Omega }}{2(1+\mathrm{\Omega }z)}$$
(32)
and

$$Q(z)=\frac{3\mathrm{\Omega }}{2(1+z)(1+\mathrm{\Omega }z)}$$
(33)
where we have used $`E(z)=(1+z)(1+\mathrm{\Omega }z)^{1/2}`$ and eq.(26). In this case it can be easily found that the function $`y_1(z)=(1+z)+4(1\mathrm{\Omega })/3\mathrm{\Omega }`$ is a solution of the eq.(19). Thus, after some calculations we can obtain the general bias solution:

$$b(z)1=y(z>\text{ }(\mathrm{\Omega }^11);\mathrm{\Omega };0)=𝒜\left[(1+z)+4\frac{(1\mathrm{\Omega })}{3\mathrm{\Omega }}\right]+𝒞\left[(1+z)+4\frac{(1\mathrm{\Omega })}{3\mathrm{\Omega }}\right]u(z)$$
(34)
with

$$u(z)=\frac{(1+z)^2\mathrm{d}z}{\left[(1+z)+4\frac{(1\mathrm{\Omega })}{3\mathrm{\Omega }}\right]^2(1+\mathrm{\Omega }z)^{1/2}}.$$
(35)
Performing the latter integration one finds that the bias evolution is given by:

$$b(z)1=𝒜\left[(1+z)+4\frac{(1\mathrm{\Omega })}{3\mathrm{\Omega }}\right]+\frac{2𝒞}{\mathrm{\Omega }}\left[(1+z)+4\frac{(1\mathrm{\Omega })}{3\mathrm{\Omega }}\right]\left[(1+\mathrm{\Omega }z)^{1/2}+\frac{8(1\mathrm{\Omega })(1+\mathrm{\Omega }z)^{1/2}}{(1\mathrm{\Omega })+3(1+\mathrm{\Omega }z)}\right]$$
(36)
It is obvious that for $`\mathrm{\Omega }1`$ the above solution tends to the Einstein-de Sitter case, as it should.

#### 4.3.2 Analytical Approximations for the $`\mathrm{\Lambda }`$ Universe

In a flat universe with non zero cosmological constant the growing mode approximation leads to the following basic factors:

$$P(z)=\frac{3}{1+z}\frac{3\mathrm{\Omega }(1+z)^2}{2[\mathrm{\Omega }(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }]}$$
(37)
and

$$Q(z)=\frac{3\mathrm{\Omega }(1+z)}{2[\mathrm{\Omega }(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }]}$$
(38)
where we have used $`E(z)=[\mathrm{\Omega }(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }]^{1/2}`$ and eq.(26). It is obvious that $`y_1(z)=E(z)`$ is a solution of eq.(19) and following a similar procedure to that of the previous subsection we obtain the general solution:

$$b(z)1=y(z>\text{ }(\mathrm{\Omega }^{1/3}1);\mathrm{\Omega };\mathrm{\Omega }_\mathrm{\Lambda })=𝒜\left[\mathrm{\Omega }(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }\right]^{1/2}+𝒞\left[\mathrm{\Omega }(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }\right]^{1/2}u(z)$$
(39)
where

$$u(z)=\frac{(1+z)^3\mathrm{d}z}{[\mathrm{\Omega }(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }]^{3/2}}.$$
(40)
The integral of equation (40) is elliptic and therefore its solution, in the redshift range $`[z,+\mathrm{})`$, can be expressed as a hyper-geometric function. We finally obtain:

$$b(z)1=𝒜[\mathrm{\Omega }(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }]^{1/2}+\frac{2𝒞}{\mathrm{\Omega }^{3/2}}[\mathrm{\Omega }(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }]^{1/2}(1+z)^{1/2}F[\frac{1}{6},\frac{3}{2},\frac{7}{6},\frac{\mathrm{\Omega }_\mathrm{\Lambda }}{\mathrm{\Omega }(1+z)^3}].$$
(41)
If $`\mathrm{\Omega }1`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }0`$ the above bias solution tends to the Einstein-de Sitter case, as it should. Note that for $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ our solution is valid even at low redshifts, since the growing mode of the fluctuations evolution is well approximated by the Einstein de Sitter solution. Therefore we present in figure 3 a comparison between the results of the high-resolution N-body simulations of Colin et al. (1999) and our solution, parametrised to $`b(0)=0.75`$ and the bias value of Colin et al. (1999) at $`z=3`$. As it is evident the agreement is excellent.

#### 4.3.3 Analytic Solution for two Limiting Cases

Although, due to the complex form of the $`P(z)`$ and $`Q(z)`$ functions in the case of $`\mathrm{\Omega }<1`$ Universes, we cannot solve the bias evolution problem analytically for all redshifts, we can produce, however, a complete analytical bias evolution solution for the limiting case of $`\mathrm{\Omega }=0`$ (Milne Universe) or $`\mathrm{\Omega }=0`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=1`$ (de Sitter Universe).

In the former case we have from eqs (16) and (17) that $`Q(z)=0`$ and $`P(z)=0`$. Applying these to eq.(19), we find the general bias solution for this special open cosmological model to be:

$$b(z)=𝒜(1+z)+$$
(42)
Interestingly, also the density fluctuations, $`\delta `$, have such a $`z`$-dependence in the $`\mathrm{\Omega }=0`$ Universe.

In the de Sitter Universe, which is dominated only by vacuum energy, putting $`\mathrm{\Omega }0`$ into eq.(16) and eq.(17) we obtain again $`Q(z)=0`$, while the $`P(z)`$ factor is given by:

$$P(z)=\frac{1}{1+z},$$
(43)
with solution

$$b(z)=𝒜(1+z)dz+.$$
(44)
Performing the latter integration one finds:

$$b(z)=𝒜(1+z)^2+,$$
(45)
with initial condition $`b_{}=𝒜+`$.

We know that the two special solutions (42) and (45) do not correspond to realistic Universes. Nevertheless, these solutions can operate as limiting cases of the generic problem. For example, using the general solution of the bias (eq.18) and putting $`\mathrm{\Omega }_{}=0`$, or $`\mathrm{\Omega }_\mathrm{\Lambda }=1`$ and $`\mathrm{\Omega }_{}=0`$, the factors $`(1+z)`$ and $`(1+z)^2`$ should survive in the general case, as well. In Figure 4, we compare the bias evolution of these special low-density models with the Einstein-de Sitter one, normalising them to the same $`b(0)`$ and using the results of Steidel et al. (1998). Of course, by no means do we imply that these models predict the same value of $`b(0)`$, we only parameterise our solution in order to compare the behaviour of $`b(z)`$ in these limiting cases.

## 5 Summary

We have introduced analytical arguments and approximations based on linear perturbation theory and a linear, scale-independent bias between a mass tracer and its underlying matter fluctuation field in order to investigate the cosmological evolution of such a bias. We derive a second order differential equation, the solution of which provides the functional form of the of bias evolution in any Universe. For the case of an Einstein-de Sitter Universe, we find an exact solution which is a linear combination of the known solution $`(1+z)`$ (cf. Bagla 1998 and references therein), derived from the continuity equation, and a second term $`(1+z)^{3/2}`$ which dominates. This solution once parametrised at two different epochs, compares extremely well with the more sophisticated halo models (cf. Mo & White 1996) and with N-body simulations.

For the two low-density cosmological models we find exact solutions, albeit only in the high-redshift approximation (where the growing mode of perturbations can be approximated by the Einstein-de Sitter solution). We also derive analytical solutions for two limiting low-density Universes (ie., $`\mathrm{\Omega }=0`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$ and $`\mathrm{\Omega }=0`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=1`$).

We thank Stefano Borgani and Peter Coles for their useful comments and suggestions. We also thank the anonymous referee for the critical comments and useful suggestions. M.Plionis acknowledges the hospitality of the Astrophysics Group of the Imperial College, were this work was completed.