text:
# Analytical approximation for the structure of differentially rotating barotropes

## 1 Introduction

Theory of self-gravitating rotating bodies seems to be an unlimited reservoir of difficult problems hardly tractable even under severe simplifications. It is a subject of scientific effort since 1742 when Maclaurin has initiated this field by his studies on incompressible rotating ellipsoids (Maclaurin, 1742). Development of modern numerical calculations resulted in progress in practical applications nowadays such as e.g. 3D hydrodynamical simulations of rotation of complex objects.

Analytical approach has succeeded for constant density, incompressible bodies. Work of Maclaurin, Jacobi, Poincare, Schwarzschild and many others has explained the behaviour of those objects almost completely (Lyttleton, 1953). Behaviour of slowly rotating polytropes has been calculated by Chandrasekhar (1936). By applying the differential equation of hydrostatical equilibrium modified by rotation he reduced the problem to an ordinary differential equation. This method however works only for a uniform rotation. This list would be incomplete without the Roche model. It’s simplicity makes it a very powerful tool for understanding behaviour of rotating objects. Present computational methods allow one to handled numerically two- and three-dimensional problems with complicated governing equations

In this paper we present simple analytical approach which can treat differentially rotating compressible barotropic stars in case of slow or moderately fast rotation. This model could fill a gap between simple analytical methods used for e.g. Maclaurin spheroids or Roche model, and complicated numerical methods such as e.g. HSCF (Hachisu, 1986), or those applying straightforward Newton-Raphson technique (Eriguchi & Müller, 1985).

## 2 Formulation of the problem

We attempt to find a density distribution (iso-density contours) of a single self-gravitating object under the following assumptions:

1. Barotropic EOS $`p=p(\rho )`$
2. Simple rotation $`𝐯=r\mathrm{\Omega }𝐞_\varphi `$ with angular velocity dependent only on the distance from rotation axis $`\mathrm{\Omega }=\mathrm{\Omega }(r)`$
3. Newtonian gravity
4. Axisymmetric density distribution
5. We seek solutions for stationary objects in full mechanical equilibrium, i.e. all quantities are time-independent

With properties (i)–(v) satisfied, the Euler equation becomes, in cylindrical coordinates ($`r`$, $`z`$, $`\varphi `$):

$$r\mathrm{\Omega }(r,z)^2𝐞_r=\frac{1}{\rho }p+\mathrm{\Phi }_g$$
(1)
Continuity equation is then fulfilled automatically. Introducing centrifugal potential:

$$\mathrm{\Phi }_c(r)=_0^r\stackrel{~}{r}\mathrm{\Omega }(\stackrel{~}{r})^2𝑑\stackrel{~}{r}$$
(2)
and enthalpy:

$$h(\rho )=\frac{1}{\rho }𝑑p$$
(3)
we get a simple equation:

$$[h(\rho )+\mathrm{\Phi }_g+\mathrm{\Phi }_c]=0$$
(4)
with a solution

$$h(\rho )+\mathrm{\Phi }_g+\mathrm{\Phi }_c=C=const$$
(5)
Equation (5) is the most important equation in the study of the structure of rotating stars under conditions (i)–(v). We define the integration constant in eq. (3) to be such that the enthalpy satisfies the condition $`h(\rho =0)=0`$. The only term which we haven’t specified yet is the gravitational potential $`\mathrm{\Phi }_g`$. If we use the Poisson equation:

$$\mathrm{\Delta }\mathrm{\Phi }_g=4\pi G\rho $$
(6)
where $`G`$ is the gravitational constant, the equation (5) becomes a non-linear second-order differential equation. This form, however, is very inconvenient, because we have to specify boundary conditions at a surface of the star,111The fact of the surface of non-rotating stars is spherical allows us to specify boundary conditions with one real number – the radius. Generally, in rotating objects we have a surface represented by some function of two variables which has to be determined. which is unknown a priori. More powerful is an integral form of eq. (5) obtained by substitution:

$$\mathrm{\Phi }_g(𝐫)=G\frac{\rho (\stackrel{~}{𝐫})}{|𝐫\stackrel{~}{𝐫}|}d^3\stackrel{~}{𝐫}$$
(7)
This integral form has been used in very successful numerical algorithm developed originally by Ostriker & Mark (1968), and recently improved by Hachisu (1986) and by Eriguchi & Müller (1985). This form will be also used to derive our approximation formula in the next section.

## 3 Approximation for the density distribution

The integral equation form of eq. (5) is:

$$h(\rho )+(\rho )+\mathrm{\Phi }_c=C$$
(8)
where $``$ is the integral operator acting on the density performing the integration on the right-hand side of eq. (7) over entire volume of the star. We define the surface of the star to be manifold consisting of points where $`\rho =0`$, $`h=0`$. Explicit form of the operator $``$ in terms of coordinates will not be needed. For a given EOS and for a fixed rotation law i.e. for given functions $`h(\rho )`$ and $`\mathrm{\Phi }_c`$, the only free parameter is the constant $`C`$. The values of $`C`$ label a family of the stellar models with the same EOS and rotation pattern, which differ in total mass and maximum density222 In differentially rotating stars the central density may be, but generally is not, the maximum density. etc.

Eq. (8) has a form of the Hammerstein non-linear integral equation (Hammerstein, 1930) and can be rewritten in a canonical form:

$$f=\left[F(f)\right]$$
(9)
where:

$$f=C\mathrm{\Phi }_ch(\rho ),F(f)=h^1(f+\mathrm{\Phi }_cC)$$
(10)
In case of linear function $`F`$, eq. (9) could be easily solved by the von Neumann series. This strongly suggests to try the following iteration scheme:

$`f_1`$ $`=[F(f_0)],`$
$`f_2`$ $`=[F(f_1)],`$ (11)
$`\mathrm{}`$
$`f_n`$ $`=[F(f_{n1})]`$
$`\mathrm{}`$
Indeed, an iteration procedure of this type was successfully applied in the so-called self-consistent field method (Ostriker & Mark 1968, Hachisu 1986). We have introduced the canonical form to ensure that the first-order approximation is found in a correct order i.e. by using the first line of the sequence (3). When we go back to non-canonical form (8) the first line of (3) takes the form:

$$C\mathrm{\Phi }_ch(\rho _1)=\left(\rho _0\right)$$
(12)
From eq. (12) above we can find the first-order deviation from sphericity. It seems impossible at first sight to avoid explicit integration in eq. (12). In case of a general $`\rho _0`$, this is true. But let us look at equation (8) in case of vanishing centrifugal potential $`\mathrm{\Phi }_c`$, i.e. with no rotation:

$$h(\rho )+(\rho )=C$$
(13)
When we use a function which satisfies eq. (13) as zero-order approximation:

$$h(\rho _0)+(\rho _0)=C_0$$
(14)
integration in eq. (12) can be easily eliminated:

$$C\mathrm{\Phi }_ch(\rho _1)=(\rho _0)=C_0h(\rho _0)$$
(15)
Finally, our formula takes the form:

$$h(\rho _1)=h(\rho _0)\mathrm{\Phi }_c+CC_0$$
(16)
or simpler, using the enthalpy ($`h(\rho _0)h_0`$, $`h(\rho _1)h_1`$):

$$h_1=h_0\mathrm{\Phi }_c+CC_0$$
(17)
Functions used as zero-order approximation ($`\rho _0`$ or $`h_0`$) are simply density and enthalpy distributions of non-rotating barotropic stars. In case of polytropic EOS, $`p=K\rho ^\gamma `$, these quantities are given by Lane-Emden functions. In more general case we have to find solution of the ordinary differential equation of the hydrostatic equilibrium.

The only unanswered question is what are the value of constants $`C`$ and $`C_0`$333 For a given EOS the function $`\rho _0`$ gives $`C_0`$ and vice versa.. It is an essential part of this work, so we decided to explain it in a separate section.

## 4 Adjusting constants

When we try to find the enthalpy distribution using the formula (17) we have to find the best zero-order function $`h_0`$ and the value of $`C_0`$ given by $`h_0`$. An equivalent problem is to find the equation for which the function $`h_0`$ and the value $`C_0`$ are best zero-order approximations – in this case we seek for $`C`$. We consider the latter case, as we have to find only one real number. Let us denote:

$$\mathrm{\Delta }C=C_0C.$$
(18)
In our approximation, the equation for the first-order enthalpy distribution $`h_1`$, in terms of the initial spherical distribution $`h_0`$ and rotation law is, from (17):

$$h_1=h_0\mathrm{\Phi }_c\mathrm{\Delta }C$$
(19)
where $`\mathrm{\Delta }C`$ is still to be determined.

The terms in eq. (19) behave as follows:

1. $`h_0`$’ is spherical enthalpy distribution, thus it is only a function of the radius, has a maximum at the centre, and goes monotonically to zero, where usually is cut. However, from mathematical point of view, Lane-Emden functions extend beyond the first zero point with negative function values.
2. $`\mathrm{\Phi }_c`$’ is a monotonically increasing function of distance from the rotation axis. It starts with zero at the rotation axis. It does not change the enthalpy along the axis of symmetry. The strongest enthalpy increase takes place along the equatorial plane.
3. $`\mathrm{\Delta }C`$’ shifts the sum of positive functions $`h_0`$ and $`\mathrm{\Phi }_c`$ down.

At first sight, shifting down by $`\mathrm{\Delta }C`$ seems not needed (i.e. one would adopt $`\mathrm{\Delta }C=0`$ ) , because we obtain correct qualitative behaviour – the star is expanded along equator. But often it is enough to introduce slow rotation to get positive value of $`h_0\mathrm{\Phi }_c`$ (i.e. for our approximation to the enthalpy in this case) for any $`r>0,z=0`$ i.e. equatorial radius becomes infinite. It leads directly to physically unacceptable results – infinite volume and mass. So the value $`\mathrm{\Delta }C`$ plays a non-trivial role and has to be found. Fig. 1 shows the behaviour of all terms in eq. (19) along the equator of the star, where the rotation acts most strongly. Horizontal lines show points where enthalpy is cut for a given value of $`\mathrm{\Delta }C`$.

We can distinguish some important values:

1. For $`\mathrm{\Delta }C<C^{(0)}`$ we obtain infinite radius of a star. These values obviously have to be rejected.
2. For $`C^{(0)}<\mathrm{\Delta }C<C^{(1)}=\mathrm{\Phi }_c(R_0)`$, where $`R_0`$ is the radius of a zero-order density distribution, we get finite volume of a star, but we use extension of $`h_0`$ with negative values. This introduces some problems which we discuss later in the article, although the resulting enthalpy and density are positive and physically acceptable.
3. For $`C^{(1)}<\mathrm{\Delta }C<C^{(2)}=h_0(r=0,z=0)`$ we get a density distribution which is topologically equivalent to the ball.
4. For $`C^{(2)}<\mathrm{\Delta }C<C^{(3)}=(h_0\mathrm{\Phi }_c)_{max}`$ we get toroidal density distribution. This case exists only if strong differential rotation is present.
5. For $`\mathrm{\Delta }C>C^{(3)}`$ the star disappears.

We expect to find the solution in the range $`C^{(0)}<\mathrm{\Delta }C<C^{(2)}`$ because we are looking for finite-volume non-toroidal stars.

One can try to find $`\mathrm{\Delta }C`$ both analytically and numerically. To keep the algebraic form and the simplicity of the formula, we now concentrate on the former method.

When we substitute the formula (19) into our basic equation (5) we get:

$$h_0\mathrm{\Phi }_c+CC_0+\mathrm{\Phi }_g(h_1)+\mathrm{\Phi }_c=C$$
(20)
In this formula we have made use of (18). After obvious simplifications, using (14) and denoting $`(\rho _0)=\mathrm{\Phi }_g(h_0)`$ we have:

$$\mathrm{\Phi }_g(h_1)=\mathrm{\Phi }_g(h_0)$$
(21)
This equality is true only if $`\rho _0=\rho _1`$. The same holds for the enthalpy:

$$h_1=h_0.$$
(22)
Using formula (19) again we finally obtain:

$$\mathrm{\Delta }C=\mathrm{\Phi }_c$$
(23)
Left-hand side of eq. (23) is constant, while the right-hand side is a function of distance from the rotation axis, monotonically decreasing from zero. This equality holds only in trivial case $`\mathrm{\Delta }C=0`$ and $`\mathrm{\Phi }_c=0`$ with no rotation at all. In any other case (23) cannot be fulfilled. So instead we try another possibility and require that

$$\mathrm{\Delta }C=\widehat{\mathrm{\Phi }_c}$$
(24)
where ‘hat’ denotes some mean value of the function $`\mathrm{\Phi }_c`$. We have chosen

$$\mathrm{\Delta }C=\widehat{\mathrm{\Phi }_c}=\left(\frac{4}{3}\pi R_{0}^{}{}_{}{}^{3}\right)^1_{V_0}\mathrm{\Phi }_cd^3𝐫$$
(25)
Integration is taken over the entire volume $`V_0=\frac{4}{3}\pi R_0^3`$ of a non-rotating initial star with the radius $`R_0`$. This choice of $`\mathrm{\Delta }C`$ gives good results. But using the mean value theorem:

$$_{V_0}\mathrm{\Phi }_cd^3𝐫=V_0\widehat{\mathrm{\Phi }_c},$$
(26)
where $`\widehat{\mathrm{\Phi }}_c`$ is some value of $`\mathrm{\Phi }_c`$ in the integration area, and taking in account monotonicity of the centrifugal potential $`0<\mathrm{\Phi }_c<\mathrm{\Phi }_c(R_0)`$ we get:

$$\mathrm{\Delta }C=\widehat{\mathrm{\Phi }_c}<\mathrm{\Phi }_c(R_0)$$
(27)
i.e. the value of $`\mathrm{\Delta }C`$ is in the range $`\mathrm{\Delta }C<C^{(1)}`$ from Fig. 1. It forces us to use negative values of non-rotating enthalpy. Moreover, in case of polytropic EOS with fractional polytropic index444 Physically interesting cases like degenerate electron gas in non-relativistic case has fractional polytropic index $`n=3/2`$. Lane-Emden equation (Kippenhahn & Weigert, 1994):

$$\frac{1}{z^2}\frac{d}{dz}\left(z^2\frac{dw}{dz}\right)+w^n=0$$
(28)
has no real negative values, because of fractional power of negative term $`w^n`$. But we can easily write equation, with solution identically equal to solution of Lane-Emden equation for $`w>0`$, and real solution for $`w<0`$ e.g.:

$$\frac{1}{z^2}\frac{d}{dz}\left(z^2\frac{dw}{dz}\right)+|w|^n=0$$
(29)
But, for example, solution of the following equation:

$$\frac{1}{z^2}\frac{d}{dz}\left(z^2\frac{dw}{dz}\right)+\frac{|w|^{n+1}}{w}=0$$
(30)
is again identically equal to solution of Lane-Emden equation for $`w>0`$, but differs from eq. (29) for $`w<0`$. Fortunately, difference between solution of eq. (29) (Fig. 1, below axis, dot-dashed) and eq. (30) (Fig. 1, solid) for $`w<0`$ is small if $`|w|1`$. Example from Fig. 1 (for $`n=1`$) is representative for other values of $`n`$. We will use form (29) instead of the original Lane-Emden equation (28) for calculations in this article.

To avoid problems with negative enthalpy we can put simply:

$$\mathrm{\Delta }C=\mathrm{\Phi }_c(R_0)$$
(31)
which is strictly boundary value $`C^{(1)}`$ from Fig. 1. The great advantage of the eq. (31) is the possibility to analytically perform the integration of the centrifugal potential (2) for most often used forms of $`\mathrm{\Omega }(r)`$. In contrast, in formula (25), not only angular velocity profile (2), but also the centrifugal potential have to be analytically integrable function. In both cases (25, 31) however, possibility of analytical integration depends on the form of $`\mathrm{\Omega }(r)`$. The value of $`\mathrm{\Delta }C`$ from eq. (31) also gives reasonable iso-density contours, cf. Fig. 11 and 8, but global accuracy is poor (Table 1).

As we noticed, the best value of $`\mathrm{\Delta }C`$ in formula (19) could be found numerically. For example, we can use virial theorem formula for rotating stars (cf. Tassoul 2000):

$$2E_k\left|E_g\right|+3pd^3𝐫=0,$$
(32)
where $`E_k`$ and $`E_g`$ is the rotational kinetic energy and the gravitational energy, respectively. We define, so-called virial test parameter $`Z`$:

$$Z=\frac{2E_k\left|E_g\right|+3(\gamma 1)U}{|E_g|}$$
(33)
where we introduced internal energy:

$$U=\frac{1}{\gamma 1}pd^3𝐫$$
(34)
Parameter Z is very common test of the global accuracy for rotating stars models. We may request that our enthalpy satisfy (32), i.e. we choose $`\mathrm{\Delta }C`$ from equation:

$$Z\left(h_0\mathrm{\Phi }_c\mathrm{\Delta }C\right)=0$$
(35)
We can find $`\mathrm{\Delta }C`$ from equation eq. (35) numerically only.

As it is shown on Fig. 2, we can find approximation of the rotating polytrope structure in form (19) satisfying virial theorem (32) up to accuracy limited only by numerical precision.

Values of $`\mathrm{\Delta }C`$ obtained with (31), (25) and from virial test (35) are compared on Fig. 3. Some of the global model properties are very sensitive to value of $`\mathrm{\Delta }C`$ (cf. Figs 2 and 5).

Because virial test is unable to check accuracy of our model, we may also try to compare directly eq. (19) with enthalpy distribution from the numerical calculations $`h_{\mathrm{num}}`$ of e.g. Hachisu (1986), Eriguchi & Müller (1985) and find $`\mathrm{\Delta }C`$ minimizing e.g. the following formula:

$$\left[h_1(\mathrm{\Delta }C)h_{\mathrm{num}}\right]^2d^3𝐫=\mathrm{min}$$
(36)
This method however, requires numerical results (e.g. enthalpy distribution) in machine-readable form.

## 5 Approximate formula accuracy

In the above sections, we tried to be as general as possible. Now we give some examples, and test accuracy of approximation.

In case of polytropic EOS $`p=K\rho ^\gamma `$ the enthalpy is:

$$h(\rho )=\frac{K\gamma }{\gamma 1}\rho ^{\gamma 1}$$
(37)
Zero-order approximation of density (density of non-rotating polytrope, Kippenhahn & Weigert 1994) with n-th555$`\gamma =1+\frac{1}{n}`$ Lane-Emden function $`w_n`$ is:

$$\rho _0=\rho _c\left[w_n(Ar)\right]^n,A^2=\frac{4\pi G}{nK\gamma }\rho _c^{\frac{n1}{n}}$$
(38)
and our formula for density becomes:

$$\rho _1=\left[\rho _{c}^{}{}_{}{}^{1/n}w_n\frac{1}{nK\gamma }(\mathrm{\Phi }_c+\mathrm{\Delta }C)\right]^n$$
(39)
where $`\mathrm{\Delta }C`$ is calculated from (25), (31) or (35).

In certain cases Lane-Emden functions are elementary functions as e.g. $`w_1`$. In cases like this our formula may be expressed even by elementary functions. For example, for $`n=1`$, $`K=1/2`$, $`4\pi G=1`$, $`\rho _c=1`$, $`\mathrm{\Omega }(r)=\mathrm{\Omega }_0/(1+r^2/A^2)`$ and $`\mathrm{\Delta }C`$ from eq. (31) we get a simple formula:

$$\rho _1(r,z)=\frac{\mathrm{sin}\sqrt{r^2+z^2}}{\sqrt{r^2+z^2}}+\frac{1}{2}\frac{\mathrm{\Omega }_0^2A^2r^2}{1+r^2/A^2}\frac{1}{2}\frac{\mathrm{\Omega }_0^2A^2\pi ^2}{1+\pi ^2/A^2}$$
(40)
Functions like this can easily be visualized on a 2D plot. Figure 1 has been made from the formula (40) while figures 8 and 11 from eq. (41).

Now we concentrate on $`n=3/2`$ polytrope. In our calculations and figures we will use $`4\pi G=1`$, $`\rho _c=1`$ and $`K=2/5`$. Now formula (39) becomes:

$$\rho _1=\left(w_n\mathrm{\Phi }_c\mathrm{\Delta }C\right)^{3/2}$$
(41)
Iso-density contours of $`\rho _1`$ from (41) are presented on Fig. 11 and Fig. 8.

To test accuracy of approximation we have calculated axis ratio, total energy, kinetic to gravitational energy ratio, and dimensionless angular momentum. Axis ratio is defined as usual as:

$$\mathrm{Axis}\mathrm{Ratio}=\frac{R_z}{R_{eq}}$$
(42)
where $`R_z`$ is distance from centre to pole and $`R_{eq}`$ is equatorial radius. Total energy $`E_{tot}`$:

$$E_{tot}=(E_k+E_g+U)/E_0$$
(43)
is normalized by:

$$E_0=(4\pi G)^2\frac{M^5}{J^2}$$
(44)
and dimensionless angular momentum is defined as:

$$j^2=\frac{1}{4\pi G}\frac{J^2}{M^{10/3}}\rho _{max}^{1/3}$$
(45)
where $`M`$ and $`J`$ are total mass and angular momentum, respectively; $`\rho _{max}`$ is maximum density. Quantities (42)-(45) are computed numerically from (41), with given angular velocity $`\mathrm{\Omega }(r)`$ and chosen $`\mathrm{\Delta }C`$.

### 5.1 Influence of $`\mathrm{\Delta }C`$

We have made detailed comparison of our $`n=3/2`$ model (41) with $`j`$-const rotation law and $`A=0.2R_0`$ (middle row of Fig. 8) for different values of $`\mathrm{\Delta }C`$ with results of (Eriguchi & Müller 1985, Table 1b). Table 1 show our results for $`\mathrm{\Delta }C`$ from eq. (25). Value of $`\mathrm{\Delta }C`$ from eq. (31) and corresponding virial test parameter $`Z`$ is included here for comparison. Table 2 shows global properties of our approximation with $`\mathrm{\Delta }C`$ equal to the solution of eq. (35), i.e. satisfying virial theorem.

Direct comparison of values from Table 1 and table Table 2 to Table 1b of Eriguchi & Müller (1985) may be difficult, because our driving parameter is central angular velocity $`\mathrm{\Omega }_0`$, while Eriguchi & Müller (1985), following successful approach of Hachisu (1986), use axis ratio (42). More convenient in this case is comparison of figures prepared from data found in Table 1b of Eriguchi & Müller (1985) and our tables. This is especially true, because axis ratio isn’t well predicted by our formula (cf. Fig. 6 and Fig. 7), while global properties ($`E_{tot}`$, $`j^2`$, $`E_k/|E_g|`$, cf. Fig. 45) and virial test $`Z`$ are in good agreement if $`E_k/|E_g|0.1`$.

Fig. 4 shows that our approximation is valid until $`E_k/|E_g|0.05`$, and begins to diverge from numerical results strongly for $`E_k/|E_g|0.1`$. Both values of $`\mathrm{\Delta }C`$ (25,35) give similar behavior here. However, $`\mathrm{\Delta }C`$ from virial test produces better results, and $`E_k/|E_g|`$ values are more sensible for the strongest rotation.

In contrast, total energy (43) is very sensitive to $`\mathrm{\Delta }C`$. Value of $`\mathrm{\Delta }C`$ from eq. (25) produces wrong result. $`E_{tot}`$ begin to increase for $`j^20.007`$, while numerical results give monotonically decreasing $`E_{tot}`$. Use of $`\mathrm{\Delta }C`$ from (35) instead, gives correct result, cf. Fig. 5.

While global properties of our model are in good agreement with numerical results for $`E_k/|E_g|0.1`$, axis ratio tends to be underestimated, even for small values of $`j^2`$. Fig. 6 and Fig. 7 show minor improvements when we use $`\mathrm{\Delta }C`$ from virial test (35) instead of mean value (25).

This subsection clearly show importance of constant value $`\mathrm{\Delta }C`$. Best results are produced with $`\mathrm{\Delta }C`$ from eq. (35), therefore this value will be used in the next subsections to investigate influence of differential rotation parameter $`A`$ and type of rotation law on formula accuracy.

### 5.2 Effects of differential rotation

In addition to the results from previous subsection ($`j`$-const with $`A=0.2R_0`$) we have calculated properties of the almost rigidly ($`A=2R_0`$) and extremely differentially ($`A=0.02R_0`$) rotating model with the same rotation law.

In all three cases we are able to find value of $`\mathrm{\Delta }C`$ satisfying eq. (35). However, this is not enough to find correct solution, because other parameters describing rotating body may be wrong. This is clearly shown on Fig. 9, where $`E_k/E_g`$ versus $`j^2`$ (45) is plotted for three cases of differential rotation. Apparent discrepancy for $`A=2R_0`$ exists. Both $`j`$-const and $`v`$-const angular velocity profiles behaves as rigid rotation in this case. Thus we conclude that our formula is unable to predict correct structure in case of uniform rotation even if rotation is small.

If rotation is concentrated near rotation axis, like in $`A=0.02R_0`$ case, our and numerical results are of the same order of magnitude. Quantitative agreement is achieved only for very small values of $`\mathrm{\Omega }_0`$. Let’s note that in this case $`\mathrm{\Delta }C`$ required by virial theorem (35) is slightly below zero (Table 4). This example shows, that $`\mathrm{\Delta }C`$ may also be negative. All three cases are summarized on Fig. 9.

Results from this section show, that our formula is able to find correct structure of rotating body for differential rotation only. Range of application vary with differential rotation parameters, and best results are obtained in middle range i.e. $`A=0.2R_0`$. With extremal case ($`A=0.02R_0`$) quality of our results is significantly degraded.

In next subsection we examine, if this statement depends on rotation law.

### 5.3 Rotation law effects

In addition to previously described cases, we have calculated global properties of our model in case of $`v`$-const angular velocity profile, with parameter $`A=0.2R_0`$ (Table 5) and $`A=0.02R_0`$ (Table 6). Results with $`A=2R_0`$ aren’t presented, because they are similar to $`j`$-const case (cf. Table 3), where both functions $`\mathrm{\Omega }(r)`$ behave as uniform rotation, and our formula fails in this case.

Figures 10 and 12 show very good agreement of of the global physical quantities ($`E_k/E_g,j^2,E_{tot}`$) with numerical results for entire range of rotation strength covered by both methods. The most extreme case ($`A=0.02`$) also behaves well. Axis ratio (Fig. 13) however, clearly distinguish between approximation and precise solution. Results are quantitatively correct only for small rotation parameters, e.g $`j^20.005`$ i.e. $`E_k/E_g0.1`$.

## 6 Discussion& Conclusions

Comparison of the results obtained with our approximation formula (Fig. 11 – Fig. 13, Table 25) with other (Eriguchi & Müller 1985, Fig. 2–5, Fig. 9, Table 1 and 2) shows a correct qualitative behaviour for even the most simplified version of our approximation formula for a wide range of parameters describing differential rotation and strength of rotation. This make our formula excellent tool for those who are interested in the structure of barotropic, differentially rotating stars, but do not need exact, high precision results. It can be applied for qualitative analysis of structure of rapidly rotating stellar cores (e.g. ‘cusp’ formation, degree of flattening, off-centre maximum density) with arbitrary rotation law, also for initial guess for numerical algorithms. It can also be used as an alternative for high-quality numerical results for use in a more convenient form as long as $`E_k/|E_g|0.1`$ and we are interested mainly in global properties of differentially rotating objects.

ACKNOWLEDGMENTS:

I would like to thank Prof. K. Grotowski and M. Misiaszek for thorough discussion of the problem, and Prof. M. Kutschera for critical reading of the previous version of this article.


format: text

text:
# Folding Mechanism of Small Proteins

## Abstract

Extensive Monte Carlo folding simulations for four proteins of various structural classes are carried out, using a single atomistic potential. In all cases, collapse occurs at a very early stage, and proteins fold into their native-like conformations at appropriate temperatures. The results demonstrate that the folding mechanism is controlled not only by thermodynamic factors but also by kinetic factors: The way a protein folds into its native structure, is also determined by the convergence point of early folding trajectories, which cannot be obtained by the free energy surface.

Understanding how a protein folds is a long-standing challenge in modern science. Computer simulations duan ; skolnick ; dinner ; go ; zhou ; galzitskaya ; shea ; kuss ; alonso ; boczko ; bursulaya ; fernandez ; ghosh  have been carried out to understand the folding mechanism. However, simulation of protein-folding processes with an atomistic model is a very difficult task. The difficulties come from two sources: (i) Direct folding simulation using an all-atom potential requires astronomical amounts of CPU time, and typical simulation times are only about a few nanoseconds. (ii) Atomic pairwise interactions including solvation effects may not be accurate enough. An extensive folding simulation has been carried out for the villin headpiece subdomain (HP-36), where 1-$`\mu s`$ molecular dynamics simulation with an all-atom potential has been performed, producing only candidates for folding intermediates duan . For this reason, direct folding simulations have been mainly focused on simple models, such as lattice models skolnick ; dinner , models where only native interactions are included (Go-type models) go ; zhou ; galzitskaya ; shea , and a model with discrete energy terms whose parameters are optimized separately for each protein kuss . Alternative indirect approaches have also been proposed including unfolding simulations shea ; alonso ; boczko ; bursulaya  starting from the folded state of a protein. However, it is not obvious whether the folding is the reverse of the unfolding dinner ; shea . Moreover, to the best of our knowledge, no one has yet succeeded in folding more than one proteins into their native states using a single potential, even with simplified models kuss .

In this letter, we propose a novel method to fold several of small proteins using Monte Carlo dynamics. This method uses a single atomistic continuous potential which includes all pairwise (native and non-native) interactions, and yet allows us to carry out folding simulations starting from non-native conformations. We observe that all proteins fold into their native-like conformations at appropriate temperatures. We find that the folding mechanism is controlled by both kinetic and thermodynamic factors lee1 : The way a protein folds into its native structure, is determined not only by the free energy surface, but also by the convergence point of early folding trajectories relative to the native state.

We consider a system where an off-lattice potential energy function including non-native interactions is utilized. We study the folding dynamics of proteins using a united-residue (UNRES) liwo  force field where a protein is represented by a sequence of $`\alpha `$-carbons linked by virtual bonds with attached united side chains and united peptide groups. Energy terms are all continuous, and include pairwise electrostatic, van der Waals, and multibody terms (see Ref. liwo  for details). The effect of solvent was indirectly included in the force field. The parameters of the UNRES force field were simultaneously optimized lee2  for four proteins of betanova bursulaya ; kortemme  (20 residues, three-stranded $`\beta `$-sheet), 1fsd dahiyat  (28 residues, one $`\beta `$-hairpin and one $`\alpha `$-helix), HP-36 duan ; fernandez  (36 residues, three-helix bundle), and fragment B of staphylococcal protein A zhou ; shea ; alonso ; boczko ; ghosh ; bai  (46 residues, three-helix bundle). The low-lying local-energy minima for these proteins were found by the conformational space annealing lee3  method. The parameters were modified in such a way that the native-like conformations are energetically more favored than the others. The global minimum-energy conformations found using the optimized force field are of the root-mean-square deviation (RMSD) values 1.5 Å, 1.7 Å, 1.7 Å and 1.9 Å from the experimental structures for betanova, 1fsd, HP-36 and protein A, respectively. After the parameter optimization, one set of the parameters is obtained for four proteins.

In the UNRES force field there are two backbone angles and two side chain angles per residue (no side-chains for glycines). The values of these angles are perturbed one at a time, typically about $`15^{}`$, and the backbone angles are chosen three times more frequently than the side chain angles. The perturbed conformation is accepted according to the change in the potential energy, following Metropolis rule. Since only small angle changes are allowed one at a time, the resulting Monte Carlo dynamics can be viewed as equivalent to the real dynamics. At a fixed temperature, at least ten independent simulations starting from various non-native states of a protein were performed up to $`10^9`$ Monte Carlo steps (MCS) for each run. During simulation the values of RMSD from the native structure and the radius of gyration ($`R_g`$) were calculated using $`C_\alpha `$ coordinates. The lowest RMSD values from folding simulations are 0.78 Å, 1.07 Å, 1.58 Å and 2.07 Å for betanova, 1fsd, HP-36 and protein A, respectively. The fractions of the native contacts ($`Q`$ and $`\rho `$) were also measured during simulations, where $`Q`$ is calculated from the experimental structure. A native contact is defined to exist when two $`C_\alpha `$’s separated more than two residues in sequence are placed within 7.0 Å. To define $`\rho `$ we first characterize the native state conformations by performing simulations starting from the experimental structures, at the same temperatures where folding simulations were performed. We define $`\rho `$ as the fraction of native contacts weighted over their contact probabilities in the native state simulations shea ; boczko ; bursulaya . Due to the fluctuation of the native conformation, the value of $`Q`$ is usually lower than that of $`\rho `$

The time histories of the typical trajectories from the folding simulations of betanova and HP-36 are shown in Fig. 1. The trajectories for 1fsd and protein A are similar to those in the figure. We observe that collapse occurs at a very early stage ($`10^4`$ MCS) for all four proteins, but the details of each folding process appear to depend on the secondary structure contents. Distributions of RMSD, $`Q`$, $`\rho `$ and $`R_g`$ are also accumulated during the whole simulations. The distributions of RMSD are shown in Fig. 2, and those for $`\rho `$ and $`R_g`$ are shown in Fig. 3 as contour plots. The early folding trajectories plotted in Fig. 3 are obtained as follows. We divided the initial $`10^5`$ MCS into 19 intervals (ten $`10^3`$ MCS and subsequent nine $`10^4`$ MCS), and took average over conformations in each interval. These averages were again averaged over 100 independent simulations starting from random conformations. The same procedure was applied to 100 independent simulations starting from a fully extended conformation.

The simulation of betanova at $`T=40`$ (arbitrary unit) starting from an unfolded conformation demonstrates that rapid collapse occurs in about $`10^4`$ MCS; the value of $`R_g`$ decreases, whereas the value of $`Q`$ remains below 0.3 (Fig. 1). During the next $`10^8`$ MCS, more compact states appear, and the value of $`Q`$ moves up as high as 1.0. There are two values of RMSD that are populated (see also Fig. 2(b)): one corresponds to less ordered structures (RMSD$``$5 Å) and the other to native-like structures (RMSD$``$2.5 Å). This demonstrates, although the native-like structure is the most stable one, thermal fluctuation can temporarily kick the protein out of the native structure.

We now analyze the details of the folding behavior of each protein. For betanova at low temperatures ($`T30`$), the probability distributions of various quantities such as RMSD depend on initial conformations, showing its glassy behavior (Fig. 2(a)). At higher temperatures ($`T40`$) this non-ergodic glassy behavior disappears. It should be noted that native-like structures are more easily found from the simulation at $`T=40`$ (Fig. 2(b)) than from the best of ten runs at $`T=30`$. When temperature decreases from $`T=80`$ to 60, the location of the RMSD peak dramatically moves from 8 Å to 3 Å. This demonstrates the cooperative folding characteristics of betanova. For $`40T60`$ betanova folds into its native-like structure. The initial folding trajectories and the distribution of ($`\rho `$, $`R_g`$) at $`T=40`$ are shown in Fig. 3(a) for betanova. Regardless of its initial conformation (either random or fully extended), the average pathways to the folded conformation initially converge to ($`\rho `$, $`R_g`$)$``$(0.35, 8.5 Å), and then they move horizontally to the native structure. This is consistent with the recent folding scenario for proteins with $`\beta `$ structure shea . The populated states bursulaya  around $`\rho `$0.4-0.5 and $`R_g`$11-12 Å are not from initial folding trajectories, but from the fluctuation of native-like structures. This kinetic information is difficult to be captured by free energy calculations alone bursulaya .

For 1fsd the distributions of various quantities for ten independent runs ($`10^9`$ MCS each) show glassy behavior for $`T50`$. Again, the non-ergodic glassy behavior disappears at higher temperatures ($`T70`$). Fig. 2(c) shows the RMSD distributions at various temperatures. The strength of the cooperativity for 1fsd is not as strong as that of betanova. Again after initial collapse to ($`\rho `$, $`R_g`$)$``$(0.3, 9 Å), the average trajectories move horizontally (Fig. 3(b)), although less prominently so compared to the case of betanova.

For HP-36, the distributions of various quantities for ten independent runs ($`10^9`$ MCS each) again show glassy behavior at $`T=60`$. At higher temperatures we have two RMSD peaks (Fig. 2(d)). At $`T=90`$ the peak near the native structure begins to dominate, and as temperature decreases it becomes stronger. This double peak feature demonstrates the cooperative two-state transition. The conformations centered at the higher value of RMSD come from a variety of collapsed states. The conformations from the other peak are native-like. When we examine them, the helix I (residues 4-8) is stably formed fernandez , while the others are fluctuating. The initial folding trajectories (Fig. 3(c)) converge to ($`\rho `$, $`R_g`$)$``$(0.2, 11 Å). These collapsed structures fold into native-like structures ($`\rho `$, $`R_g`$)$``$(0.7, 10 Å). Compared to the case of betanova and 1fsd, the average folding trajectories are more diagonal.

The overall folding characteristics of protein A are similar to those of HP-36. The initial folding trajectories (Fig. 3(d)) converge to ($`\rho `$, $`R_g`$)$``$(0.25, 12 Å). These collapsed states fold into native-like structures in a diagonal fashion similar to the case of HP-36. When we examine the native-like conformations with $`3\mathrm{\AA }\mathrm{RMSD}4\mathrm{\AA }`$, the helix III (residues 42-55) is most stably formed. This is in agreement with recent investigations alonso ; ghosh ; bai . We also observe that a first-order like collapse transition zhou  (from ($`Q`$, $`R_g`$)$``$(0.15, 18 Å) to (0.15, 12 Å)) occurs near $`T=120`$.

By using an atomistic model, we could observe folding processes of four small proteins in realistic settings. In all cases, rapid collapse is followed by subsequent folding process that takes place in a longer time scale. The folding mechanism suggested in this study is as follows: There are two aspects of folding dynamics, (i) non-equilibrium kinetic properties and (ii) equilibrium thermodynamic properties (Fig. 4). The non-equilibrium kinetic properties, relevant to the early folding trajectories (fast process), can be examined only by direct folding simulations. The free energy surface, an equilibrium thermodynamic property, dictates the way an initially collapsed state completes its folding (slow process). The way a protein folds into its native structure, i.e., either horizontally or diagonally in ($`\rho `$, $`R_g`$) plane, is determined by the position of ($`\rho `$, $`R_g`$) where early folding trajectories converge, relative to the native state. It appears that slow folding process of $`\alpha `$-proteins occurs in a diagonal fashion compared to that of proteins containing $`\beta `$-strands shea .

In conclusion, we successfully carried out direct folding simulations of more than one proteins using a single atomistic potential. We also observe that glassy transitions occur at low temperatures. The results provide new insights into the folding mechanism.


text:
# Static impurities in the kagome lattice: dimer freezing and mutual repulsion

## I Introduction

The kagome geometry (Fig. 1) presents one of the most highly frustrated quantum spin systems achievable in two dimensions (2d) with only nearest–neighbor Heisenberg interactions. The $`S`$ = 1/2 kagome lattice has been found to have a spin–liquid ground state,rce  whose ultra–short spin–spin correlation lengths are effectively those of a dimer liquid,rle ; rze  and whose excitation spectrumrlblps  shows a manifold containing an extremely large number of low–lying singlets.rm  These properties are reproduced very well by a short–range resonating–valence–bond (RVB) description,rmm  based in fact only on nearest–neighbor dimer formation (Fig. 1).

On the theoretical level, a similar degree of frustration is found only in quantum dimer models,rrk ; rms  which have also been considered recently in the kagome geometry,rmsp  and in models containing higher–order spin terms and multiple–spin exchange.rmblw  Experimentally, a variety of materials displaying the kagome structure is known to exist, although to date none have been found which contain an ideal, 2d, vacancy–free system of spins $`S=1/2`$. These include the jarosites (H3O)Fe3(OH)6(SO4)2,rwhrs  KFe3(OH)6(SO4)2, and KCr3(OH)6(SO4)2rkea ; rlea  and the volborthite Cu3V2O7(OH)2.H2O.rhhkntkt  The best characterized of these compounds is the magnetoplombite SrCr9pGa12-9pO19 (SCGO),rlmcmombm  in which some of the $`S`$ = 3/2 Cr3+ ions form approximately 2d kagome planes. However, because in this case the ideal stoichiometry ($`p=0`$) remains unachievable, discussion of the influence of static (and also, for Ga3+, spinless) impurities represents an important issue.

Studies of static vacancies in spin liquids have tended to focus on the question of induced free spin degrees of freedom.rlmgmdhlr  In systems such as spin ladders, and also for the 2d square lattice,rbhsl  the total energy appears to decrease as such holes in the spin background are made to approach each other, implying an effective mutual attraction. It is generally believed that in an RVB description of the ground state the minimal disruption of the wave function, and thus the most favorable energy state, occurs when holes occupy adjacent sites, thus generating a simple pairing mechanism for mobile dopants. The issue of impurites in frustrated magnets has received less attention, and appears to be more involved: it has been proposed that an absence of induced free spins corresponds to deconfined, spinon-like excitations.rsv  Consideration of the problem for the 1d frustrated chain,rnm  which provides an example with no free spins localized around the vacancies, suggests the importance of degenerate singlets in the ground–state manifold. This result indicates that the kagome geometry presents the most probable candidate for exotic behaviorrs  as a consequence of strong frustration in 2d.

In this study we investigate the effects of static impurities in the kagome lattice. We find that, despite the very short intrinsic spin–spin correlation lengths, dimer–dimer correlations develop over considerable distances in the presence of dopants. Further, we show that there is no evidence for free local moments induced around impurity sites, and that a highly unusual effective repulsive interaction arises between such holes in the spin background.

In Sec. II we provide an outline of the numerical and anaytical methods to be employed, and discuss the fundamental effects on the excitation spectrum of doping by spinless impurities. In Sec. III we consider the case of one hole in an odd–site cluster, demonstrating “dimer freezing” in the spin–spin correlation functions and dimer–dimer correlations of anomalously long range. The induced magnetizations in the spin sectors $`S`$ = 1 and higher are used to illustrate the consequences for nuclear magnetic resonance (NMR) measurements. Sec. IV presents spin–spin correlation functions which characterize the situation with two impurities, supplemented by total–energy analyses quantifying the novel repulsive interaction between the holes which emerges from this study. Sec. V contains further discussion and our conclusions.

## II Elimination of singlet states and absence of localized spins

### II.1 Exact Diagonalization

By “impurity” we refer henceforth to nonmagnetic impurities. The techniques we employ are based primarily on numerical calculations for small clusters with periodic boundary conditions chosen as symmetrically as possible. Specifically, we perform exact diagonalization (ED) studies of the Heisenberg Hamiltonian

$$H=J\underset{ij}{}𝐒_i𝐒_j,$$
(1)
where $`J`$ is the antiferromagnetic superexchange interaction and $`ij`$ denotes nearest–neighbor sites, for clusters up to 27 sites with one impurity, and up to 24 sites with two impurities for all available configurations. We analyze the total energies, the energies per doped impurity, and also the form of the energy spectra, focusing on the number of singlet states in the low–energy manifold below the first triplet for the finite systems under consideration, as well as on the number of these states removed by impurity doping. The frustrated spin interactions in the presence of impurities are characterized for the entire cluster by the bond spin–spin correlation functions and by the induced site magnetizations in every spin sector.

The density of states integrated over the lowest–lying energy levels is shown in Fig. 2 for systems with and without impurities. Fig. 2 compares ED results from systems of 18 sites, for which in the impurity configuration illustrated 8 of the 12 singlet states occurring below the first triplet, indicated by the vertical lines, are eliminated by the presence of 2 impurities. Although the separation of the triplet from the lowest remaining singlet decreases, this energy gap always remains a finite fraction of the superexchange interaction $`J`$. Thus the most important qualitative result emerging even for the smallest systems is that, in contrast to spin–ladder and dimerized or quadrumerized systems, no low–lying triplet state is induced by the presence of impurities. This implies directlyrnm  that there are no quasi–free spin degrees of freedom, by which is meant spins effectively isolated from their local environment and interacting with each other only on an effective energy scale much smaller than $`J`$. The fact that no localized moments are formed in the vicinity of impurities may be ascribed to the continuum of singlets beyond the ground state, which has the same function as the two–fold degeneracy of the ground state in the frustrated chain,rnm  namely that the rearrangement of dimers incurs no energy penalty. The absence of local spins around nonmagnetic impurities represents one of the most significant differences between frustrated and unfrustrated systems, and has been further interpreted as implying a deconfinement of the spin degrees of freedom.rsv

Table 1 provides more details of the number of states in the singlet manifold which are removed by impurity doping for clusters of 12 and 18 sites, for all possible configurations of 2 impurities. The numbers vary little as a function of the proximity of the impurities, and the average number of low–lying states which are eliminated is in agreement with simple considerations of the number of singlet states in the small system (below). We have also performed analogous calculations for 24-site clusters of different shapes. For all clusters the same behavior is found upon introducing two impurities, in that the number of singlets below the first triplet is reduced, but here the lowest–lying excitation is not a triplet. However, these clusters do not appear to be representative of the generic kagome physics, as the number of low–lying singlets in the absence of impurities is significantly smaller than the expected valuerlblps  of $`1.15^N`$. This situation is presumed to be a consequence of the low symmetry of the 24–site clusters, some of which are essentially 1d systems, and as a result the densities of states are not presented here.

### II.2 Dimer–RVB Basis

We calculate in addition the ground–state properties, namely the energy and spin–spin correlation functions, of clusters doped with one and with two impurities within the dimer–RVB basis. This is the basis of states in which every spin forms a singlet dimer with one of its nearest neighbors in such a way that every spin of the lattice is involved in a dimer (Fig. 1). The use of this basisrmm  naturally implies a truncation of the total number of available states in the manifold of low–lying singlets, but we will show below that it provides a satisfactory description of the ground state of the doped system. Dimer–RVB calculations permit the consideration of somewhat larger cluster sizes, and we will present results up to and including 39 sites with one impurity and 36 sites with two impurities.

Calculations performed in the dimer–RVB basis yield results very similar to the exact ones for the integrated density of states (Fig. 3), giving clear evidence of a strong linear correlation (see also Fig. 9) between the ED and dimer–RVB spectra in both 0– and 2–impurity cases. Demonstration of a quantitative correspondence between the spectra would require shift and scaling factors arising from an energetic renormalization related to the fact that calculations in the dimer–RVB basis are variational in nature, so that the true eigenstates are dressed RVB wave functions, but this procedure is not necessary for the analysis to be presented here.

Given that the RVB basis represents a very significant truncation of the total Hilbert space, the qualitative similarity of the dimer–RVB spectra to those obtained from the exact calculations, both with and without impurities, provides important information about the nature of the exact wave functions of the kagome system. This observation is further confirmed by the semi–quantitative agreement with exact results for spin correlation functions and ground–state energies computed using the dimer–RVB basis. We defer this comparison to Secs. III and IV, and comment here only that our results in the presence of impurities provide additional strong justification for the statement that the ground state of the kagome spin system is described appropriately by a maximally short–ranged RVB wave function, presumably as a consequence of the very strong frustration.

The dimer–RVB basis is a natural framework for understanding the absence of localized spins around an impurity site, as all remaining spins are simply incorporated within rearranged singlet dimer coverings. It also provides a qualitative reflection of the “evaporation” of low–lying singlets from the ground–state manifold (Figs. 2–3 and Table I) which arises on doping a cluster with two impurities, and in fact yields additional insight into the mechanism for this exclusion. It can be shown that the total number of singlet dimer coverings of a doped system is given by

$$N_{\mathrm{coverings}}=2^{N/3N_{\mathrm{imp}}+N_G+N_{\mathrm{et}}},$$
(2)
where $`N`$ is the number of sites of the system under consideration, $`N_{\mathrm{imp}}`$ the number of doped impurities, $`N_G`$ the number of independent clusters into which the system is divided by the dopant distribution, and $`N_{\mathrm{et}}`$ the number of triangles containing 3 impurities. A complete derivation of this result lies beyond the scope of the current analysis, and we present only the following heuristic explanation. By reducing each triangle to an effective $`S=1/2`$ degree of freedom,rm ; rmm  the dimer basis involves a truncation of the $`2^N`$ total states to $`2^{N/3}`$, multiplied by an additional factor of 2 for the overall spin direction on each independent cluster. Each added impurity removes one spin, halving the number of states in every manifold and thus justifying the dependence on $`N_{\mathrm{imp}}`$. The only exception to this statement occurs when the added impurity falls on a triangle which already contains two impurities, as a result of which it is already fully constrained and removal of the third site no longer affects the number of coverings. Eq. (2) can be shown to remain valid for all vacancy concentrations, including most notably values of $`N_{\mathrm{imp}}`$ exceeding the percolation threshold, beyond which $`N_G>1`$, under the condition that the impurity distribution does not cut the system into segments containing odd numbers of spins.

## III Spin correlation functions and induced magnetization

### III.1 Dimer “freezing”

In this section we consider the effects of a single impurity introduced in an odd–site cluster. The bond spin–spin correlation functions $`𝐒_i𝐒_j`$ around a single impurity may be computed by ED for clusters of up to 27 sites. The ED calculations permit in addition the determination of the magnetizations induced on each site by the presence of the impurity, with results for every spin sector corresponding to the situation in all applied magnetic fields, which are presented in Subsec. C below.

The correlation functions obtained by ED are illustrated in Fig. 4 for a 27–site cluster with periodic boundary conditions. The most striking feature of these results is that dimer bonds are effectively frozen in the two depleted triangles, taking a value $`𝐒_i𝐒_j0.69`$ close to the theoretical minimum of $`3/4`$ for a pure singlet state. For comparison, the average value of $`𝐒_i𝐒_j`$ in the pure system, which retains full freedom to resonate, is a factor of 3 smaller ($`0.22`$).rle  The bonds neighboring the frozen dimers then show very weak correlations, while the next neighbors are again enhanced, and the spin correlations oscillate with distance away from the impurity \[also illustrated in Fig. 5(b)\]. Thus although spin–spin correlations in the kagome system are extremely short–ranged,rle  induced dimer–dimer correlations show a considerable spatial extent, reaching in fact well beyond the sizes obtainable by studying small clusters using either ED or dimer–RVB techniques. For the limited data available it is unfortunately not possible to specify the functional form of the decay to equilibrium. Further analysis of this issue, based on data obtained from mutual hole repulsion, is presented in Sec. IV.

### III.2 Dimer–RVB Basis

Fig. 5(a) compares the bond spin correlation values $`𝐒_i𝐒_j`$ obtained from ED and from the dimer–RVB basis for all bonds in clusters of 21 and 27 sites with one impurity. These clusters represent the maximal sizes accessible by ED, and for such systems it is clear that the dimer–RVB basis reproduces the exact values with quantitative accuracy. From the nature of the RVB basis it is perhaps to be expected that nearest–neighbor bond correlations are treated more accurately than any other physical property. We make use of this feature to extend the calculations of spin correlation functions to a system of 39 sites with oie impurity. Because on this cluster the boundary conditions lead to a rather low symmetry, Fig. 5(b) shows the value of $`𝐒_i𝐒_j`$ averaged over all bonds at the same Manhattan distance from the impurity, which illustrates both the oscillatory nature of the frozen dimerization pattern and the very long range over which an impurity acts to disrupt the spin configuration.

### III.3 Induced magnetic moments

The local magnetization pattern induced in the vicinity of a doped nonmagnetic impurity may be computed by exact diagonalization in the sector of total spin $`(S,S_z)`$ = (1,1), and is displayed for the 27–site cluster in Fig. 6. The site magnetizations exhibit a certain anticorrelation with the strength of the dimers formed on each bond, and thus have a tendency to oscillate with increasing distance from the impurity. The induced moments are particularly enhanced on sites far from the impurity, where different paths on the periodic cluster interfere, and the interference in this regime may even result in moments oriented opposite to the effective magnetic field (Figs. 6, 7). The magnetization profile of the surrounding sites is the single most important consequence of impurities which is reflected in local–probe measurements of doped spin systems.

Fig. 7 shows the magnetization per site induced by proximity to a nonmagnetic impurity for each spin sector obtainable by ED for a 27–site cluster. The spin sector may be considered as equivalent to an effective magnetic field which varies between zero ($`S_z=0`$) and saturation ($`S_z=13`$). The most pronounced variation is observed for the nearest–neighbor sites (filled circles in Fig. 7), which have very small magnetization at low fields because they are bound into the strong dimers frozen by the presence of the impurity on the same triangle. However, between $`S_z=4`$ and $`S_z=6`$ the moment shows an abrupt jump to saturation, implying the breaking and full polarization of these dimers at fields in excess of 1/3 of the saturation field. We speculate that this behavior may be due to the presence of a particularly stable state of the remaining spins which corresponds to a 1/3 magnetization plateau of the undoped kagome lattice.rhcs  The moments of the further–neighbor sites show a generally less systematic evolution, usually with a non–monotonic character which implies a very significant rearrangement of the induced local magnetization as a function of field. This evolution is again marked by the most significant changes occurring close to the net magnetization value of 1/3.

NMR has emerged as the technique of choice for analyzing in real space the effects of doped impurities on magnetizations and spin correlations. Information concerning these quantities is extracted from Knight–shifts, linewidth alterations and, for certain cases, from relaxation times. A detailed NMR and NQR (nuclear quadrupole resonance) analysis has been conducted recently for 71Ga impurities in the kagome planes of the $`S`$ = 3/2 spin system SrCr9pGa12-9pO19 (SCGO).rlmcmombm  While a quantitative determination of any shifts or linewidths depends on a detailed Hamiltonian for the hyperfine interactions in a specific material, the results obtained in Fig. 7 may be considered as illustrative for the kagome geometry.

The important qualitative features are the absence of a single spin–1/2 degree of freedom anywhere in the vicinity of the impurity, and the rather large extent of the area of affected sites in spite of the very short spin correlation lengths. The latter effect is contained in Fig. 7 in that the induced site magnetization values on the vertical axis are directly proportional to the peak shifts measured in an NMR field sweep. Impurities in the kagome lattice therefore produce a large number of satellite peaks.

Finally, with regard to the fact that most known kagome systems, including SCGO, are not $`S=1/2`$, there remain rather few studies of this geometry for higher values of the spin.rhcs ; rrsr  On rather general grounds based on the geometry of the lattice, one expects that both triangles and depleted triangles retain the integral or half–integral nature of the spin per site. Thus it is at least plausible that the results for the $`S`$ = 1/2 system are directly relevant to higher half–integral–spin systems such as SCGO.footnote2

## IV Hole repulsion

In this section we present results characterizing the behavior of an even–site cluster doped with two impurities. A comparison with the undisturbed systemrlblps  reveals that the average energy per remaining site not only increases in the presence of impurities, but increases further as the impurities are brought closer together. Both of these features are evident from the ED results of Fig. 8, which shows the energies of all different two–hole states available for the relevant system sizes. The higher energy per site of doped systems confirms, by considering the contributing processes, that the net loss of resonance energy for dimers constrained by the presence of impurities to occupy only the bonds on which they are essentially “frozen” (Sec. III) is not compensated by the gains in dimerization energy for those bonds which are enhanced.

The effective hole repulsion is most unusual in low–dimensional spin systems,rbhsl ; rdrs  and in fact we are unaware of analogous results for any model other than the frustrated chain.rrnm  This repulsion of holes constitutes the central result of our analysis: particularly in the context of short–range RVB wave functions, it contradicts directly the most straightforward expectation for the behavior of doped holes, namely that when mobile these would choose to minimize the net disruption of the system by occupying the same dimer, and thus would exhibit a pairing attraction.ra  The consequences for mobile impurites with repulsive interactions may be expected to include hole crystallization, phase separation and the possibility of transitions between metallic and insulating phases, and would bear considerable further investigation.

### IV.1 Dimer–RVB Basis

We illustrate in Fig. 9 the efficacy of the dimer–RVB basis for the purposes of comparing ground–state energies. The correspondence between the dimer–RVB and ED results is shown by comparing the energies of all two–hole states in clusters of 12, 18 and 24 sites with two impurities, which represent all system sizes accessible by ED. Full details of the extraction of the energy values $`\mathrm{\Delta }_2`$ in Fig. 9 are given below. Although an inspection of the abcissae reveals that the agreement between individual pairs of values is not as convincing as for the spin correlation functions (Sec. III), the main difference between the two sets of data may be reduced to a relative upward shift of all the RVB energies. The direct linear correlation between ED and dimer–RVB results (alluded to in Sec. II) demonstrates unambiguously that comparing the energies of different configurations of the two impurities, which is the primary goal of the energy analysis to follow, is a valid qualitative exercise when performed using the dimer–RVB results.

We show first the spin–spin correlation functions in the presence of a pair of impurities (Fig. 10). When considering the overlap of single–impurity dimerization patterns (Sec. III), it is evident that these are not fully mutually compatible. In fact, one finds by inspection that not only the configuration illustrated in Fig. 10 but all possible configurations of two impurities on the kagome lattice are in some measure incompatible. Comparison with other quantum spin systems suggests that this complete incompatibility is very specific to the kagome geometry, and its extent may be quantified by extracting the effective interaction energy of the holes (Fig. 11). In cluster calculations the frustration is also shown by the presence of relatively large local moments induced on the sites at the mismatch between dimerization regimes (cf. Fig. 6).

### IV.2 Interaction between impurities

Fig. 11 shows twice the energy per hole as a function of separation for all configurations of two impurities on 30– and 36–site clusters, obtained in the dimer–RVB basis. We have analyzed the data in the following manner. $`E_0`$ is the ground–state energy of each cluster in the absence of impurities, and is computed for even clusters then extrapolated for odd clusters from the results for systems of sizes $`N3`$ and $`N+3`$. $`E_1`$ is the ground–state energy in the presence of a single impurity, which is computed for odd clusters and may then be extrapolated for even clusters from the results for systems of sizes $`N3`$ and $`N+3`$. $`E_2`$ is the ground–state energy with two impurities, and again may be computed directly for even clusters and extrapolated for odd clusters. We then define

$$\mathrm{\Delta }_1=E_1E_0$$
(3)
as the energy of one hole,

$$\mathrm{\Delta }_2=E_2E_0$$
(4)
as the energy of a two–hole state, and

$$ϵ=\mathrm{\Delta }_22\mathrm{\Delta }_1=E_22E_1+E_0$$
(5)
as the energy of interaction between two holes. Because of the extrapolation procedure it is important that the data sets for different system sizes be comparable. While certain problems arise for some system sizes due to the relatively large influence of the boundary conditions, we find the process to be valid for systems in excess of 27 sites. In this connection we note that the two–hole energies $`\mathrm{\Delta }_2`$ shown in Fig. 11 are independent of the extrapolation between different system sizes, and it is only the direct comparison with $`\mathrm{\Delta }_1`$ for the 33–site system which requires validation of the procedure.

It is clear that the effective mutual repulsion decays only slowly. Because of the restricted sizes of the acccessible clusters we are unable to speculate on the functional form of this decay, and in particular on the possibility that it may be long–ranged. We note only that for short distances (3–4 lattice constants) it appears to exhibit an approximate $`1/d_r`$ dependence, beyond which a steeper decrease is observed which may be either an intrinsic property of the interaction or a consequence of the system size. The form of the decay appears similar to that of the bond spin correlation functions towards their undisturbed value \[Fig. 5(b)\].

At the phenomenological level, a repulsive interaction between two impurities may be regarded as surprising in that their interference might have been expected to suppress the freezing of individual dimer patterns, leading to a restoration of resonance energy. Our findings to the contrary indicate that the dominant effect of increasing increasing impurity separation lies in the reduction in strength of the partially frozen single–impurity dimerization patterns at the boundaries where they are forced to mesh, which reduces the corresponding energetic penalty. Given the slow decay of the dimer–dimer correlations, the possibility remains that the net repulsive interaction between holes could be long–ranged. Such an unexpected result may indeed be connected with a true deconfinement of spinon–like excitations for the kagome system, which would confirm speculation to this effect in 2d.rsv  We comment further on the issue of deconfinement below.

## V Conclusion

We have illustrated a variety of novel phenomena associated with the effects of static, spinless impurities doped into the $`S=1/2`$ antiferromagnetic Heisenberg model in kagome geometry. By analysis of spin–spin correlation functions and of the local magnetization patterns induced around impurities, we find that, despite the very short spin–spin correlation lengths of the highly frustrated system, dimer–dimer correlations develop over considerable distances due to a hole–induced freezing of their resonance. Furthermore, as a consequence of the continuum of states in the low–lying singlet manifold, there is no evidence for free local moments in the vicinity of impurity sites.

Most importantly, we have shown from the total energies of doped systems that impurities in the kagome lattice have an effective repulsive interaction, a phenomenon unique to our knowledge and in direct contradiction to the conventional expectation for short–range RVB systems of a net pairing attraction. Furthermore, this effective repulsion is found to be rather long–ranged in nature, a most unusual result which is presumably related with the very strong frustration inherent to the kagome geometry.

Our comparisons of doping effects provide additional justification for the validity of the short–range dimer–RVB framework as a suitable description of the ground–state singlet manifold in the kagome antiferromagnet. We have shown that the dimer–RVB basis provides a valuable means of probing the nature of spin and charge degrees of freedom in a further class of quantum magnets. This realization of a highly resonant RVB state, arising as a consequence of the extremely frustrated nature of the system, also presents an important contrast to other low–dimensional spin systems which are known to be described by RVB–type wave functions, particularly with regard to the localization of free spins and to the pairing of holes.

We summarize only briefly the consequences of our analysis for experiment. The majority of known kagome materials do not have complete filling of the lattice sites by a single value of the spin $`S`$, and as such the study of impurities is a significant component of their characterization. Nonmagnetic impurities in the kagome lattice have the important property that they do not generate localized spins at the vacant site, and do generate by the nature of their extended range of influence a significant number of satellites for each peak. On heuristic grounds we believe that our analysis of the qualitative features to be expected in the $`S=1/2`$ system are relevant also for $`S=3/2`$, and perhaps for higher half–integral spins.

Finally, on the issue of deconfinement in highly frustrated systems of dimensionality higher than one, we have shown clearly from the presence of the mutual repulsive interaction between holes that charge degrees of freedom are deconfined. Further, the absence of localized $`S=1/2`$ degrees of freedom around vacancies in the lattice demonstrates that spin and charge are deconfined. However, on the question of whether the spin excitations of the undoped system are appropriately described by deconfined spinon excitations, our results do not provide a proof of this hypothesis.

###### Acknowledgements.

We are grateful to P. Sindzingre for helpful provision of ED spectra for undoped clusters. This work was supported by the Swiss National Science Foundation and by the French Centre National de Recherche Scientifique.


text:
# Galaxy Bias in Quintessence Cosmological Models

## Abstract

We derive the evolution of the linear bias factor, $`b(z)`$, in cosmological models driven by an exotic fluid with an equation of state: $`p_x=w\rho _x`$, where $`1w<0`$ (quintessence). Our aim is to put constrains on different cosmological and biasing models by combining the recent observational clustering results of optical (2dF) galaxies (Hawkings et al.) with those predicted by the models. We find that our bias model when fitted to the 2dF clustering results predicts different bias evolution for different values of $`w`$. The models that provide the weak biasing ($`b_{}1.1`$) of optical galaxies found in many recent observational studies are flat, $`\mathrm{\Omega }_\mathrm{m}=0.3`$ with $`w0.9`$. These models however, predict a weak redshift evolution of $`b(z)`$, not corroborated by N-body simulations.

Keywords: cosmology: theory - large-scale structure of universe

## 1 Introduction

The statistical properties of the distribution of matter on large scales, based on different extragalactic objects, can provide important constrains on models of cosmic structure formation. However, a serious problem that hampers the straight forward use of such an approach is our ignorance of how luminous matter trace the underlying mass distribution. Many authors have claimed that the large scale clustering of different mass tracers (galaxies or clusters) is biased with respect to the matter distribution (cf. Kaiser 1984), which is an essential ingredient for cold dark matter (CDM) models to reproduce the observed galaxy distribution (cf. Davis et al. 1985). Usually, biasing is assumed to be statistical in nature by which galaxies and clusters are identified as high peaks of an underlying, initially Gaussian, random density field.

Furthermore, the bias redshift evolution is very important in order to relate observations with models of structure formation and it has been shown that the bias factor, $`b(z)`$, is a monotonically increasing function of redshift. Indeed, Mo & White (1996) and Matarrese et al. (1997) have developed a model for the evolution of the correlation bias, defined as the ratio of the halo to the mass correlation function, the so called galaxy merging bias model, with $`b(z)(1+z)^{1.8}`$. Nusser & Davis (1994), Fry (1996), Tegmark & Peebles (1998) produced a bias evolution model assuming only that the test particles fluctuation field is proportional to that of the underlying mass and in this case, we have: $`b(z)=1+(b_{}1)(1+z)^1`$, where $`b_{}`$ is the bias factor at the present time. Coles, Melott & Munshi (1999) have developed a bias model within the hierarchical clustering paradigm which gives interesting scaling relations for the galaxy bias. Recently, another approach was proposed by Basilakos & Plionis (2001), in which the bias evolution was described via the solution of a second order differential equation, derived using linear perturbation theory and the basic Cosmological equations.

From the observational point of view, recent advances in cosmology have strongly indicated that we are living in a flat, accelerating universe with low matter (baryonic and dark matter) density (cf. Riess et al. 1998; Perlmutter et al. 1999; de Bernadis et al. 2000; Efstathiou et al. 2002; Persival et al. 2002, Spiergel et al 2003). The available high quality cosmological data (Type Ia supernovae, CMB, etc.) are well fitted by an emerging “standard model”, which contains cold dark matter (CDM), to explain clustering, and an extra component with negative pressure, usually named “quintessence”, which is in agreement with the inflationary flatness prediction ($`\mathrm{\Omega }_{\mathrm{tot}}=1`$) as well.

In this paper we derive the redshift evolution of the linear scale-independent bias parameter for the “quintessence” models, and additionally we attempt to put constraints on the different cosmological models using the recent clustering results (Hawkins et al. 2002) of optical galaxies from the 2dF survey. We compare the observed clustering of 2dF galaxies with that expected in 10 spatially flat cosmological models with dark energy and $`\mathrm{\Omega }_\mathrm{m}=0.3`$. It is worth noting that the galaxy bias has been found to be sensitive to the equation of state (Munshi, Porciani & Wang 2003).

## 2 Linear Bias Evolution Model

In this section we briefly describe our linear bias evolution model (for more details see Basilakos & Plionis 2001) which is based on linear perturbation theory in the matter dominated epoch (cf. Peebles 1993). We assume that the mass tracer population is conserved in time, ie., that the effects of non-linear gravity and hydrodynamics (merging, feedback mechanisms etc) can be ignored (cf. Fry 1996; Tegmark & Peebles 1998; Catelan et al. 1998), and using the linear perturbation theory we can obtain a second order differential equation, which describes the evolution of the linear bias factor, $`b`$, between the background matter and a mass-tracer fluctuation field:

$$\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),$$
(1)
where $`D(t)`$ is the linear growing mode, useful expressions of which can be found for the $`\mathrm{\Lambda }`$CDM model in Peebles et al. (1993) and for the QCDM models in Wang & Steinhardt (1998). The solutions of this second order differential equation provides our bias evolution model. The solution for the different cosmological model enter through the different behaviour of $`D(t)`$ and $`H(t)`$. To see the relevance and difference of our bias evolution model to the Galaxy Conserving one, which is based on similar assumptions, we will derive somewhat differently eq.(1), starting from the continuity equation, which is the starting point of the Galaxy Conserving model. If the galaxies and the underlying matter share the same velocity field then: $`\dot{\delta }+u0`$ and $`\dot{\delta }_g+u0`$, from which we get $`\dot{\delta }\dot{\delta }_g=0`$. Now since we assume linear biasing, we have $`\delta _g=b\delta `$ and using $`y=b1`$, we get that $`d(y\delta )/dt=0`$. Differentiating twice we then have: $`\ddot{y}\delta +2\dot{y}\dot{\delta }+y\ddot{\delta }=0`$. Solving for $`y\ddot{\delta }`$, using the fact that $`y\dot{\delta }=\dot{y}\delta `$ and utilizing the differential time-evolution equation of $`\delta `$ (cf. Peebles 1993) we finally obtain:

$$\ddot{y}\delta +2(\dot{\delta }+H\delta )\dot{y}+4\pi G\rho \delta y=0$$
(2)
which is the corresponding equation 1.

Here we extend our previous bias evolution solutions of Basilakos & Plionis (2001) to flat Friedmann-Robertson Walker Cold Dark Matter (CDM) type models driven by non-relativistic matter and having an exotic fluid (quintessence) with equation of state:

$$p_x=w\rho _x(1w<0).$$
(3)
Following the notation of Basilakos & Plionis (2001), we present the basic steps of our procedure. In order to transform (1) from time to redshift we utilize the following expressions:

$$\frac{dt}{dz}=\frac{1}{H_{}E(z)(1+z)},E(z)=\left[\mathrm{\Omega }_\mathrm{m}(1+z)^3+\mathrm{\Omega }_{x_o}(1+z)^\beta \right]^{1/2}$$
(4)
while $`\mathrm{\Omega }_\mathrm{m}=8\pi G\rho _o/3H_o^2`$ (density parameter), $`\mathrm{\Omega }_{x_o}=8\pi G\rho _{x_o}/3H_o^2`$ (dark energy parameter) at the present time, which satisfy $`\mathrm{\Omega }_\mathrm{m}+\mathrm{\Omega }_{x_o}=1`$, and $`\beta =3(1+w)`$ with $`0\beta <3`$. Taking into account the latter transformations, 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)$$
(5)
with relevant 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}$$
(6)
and

$$Q(z)=\frac{3\mathrm{\Omega }_\mathrm{m}(1+z)}{2E^2(z)}.$$
(7)
Therefore, the general bias solution for all of the flat cosmological models is:

$$b(z)=𝒜E(z)+𝒞E(z)I(z)+1,$$
(8)
where

$$I(z)=\frac{(1+z)^3}{E(z)^3}dz.$$
(9)
The integral of equation (9) is elliptic and therefore its solution, in the redshift range $`[z,+\mathrm{})`$, can be expressed as a hyper-geometric function. Our final general solution is:

$$b(z)1=𝒜E(z)+\frac{2𝒞}{\mathrm{\Omega }_\mathrm{m}^{3/2}}E(z)(1+z)^{\frac{\beta 3}{62\beta }}F[\frac{1}{62\beta },\frac{3}{2},\frac{72\beta }{62\beta },\frac{\mathrm{\Omega }_x}{\mathrm{\Omega }_\mathrm{m}(1+z)^{3\beta }}].$$
(10)
Note that the first term, which is the dominant one, has an approximate redshift dependence $`(1+z)^{3/2}`$ while the second has $`(1+z)`$. Also note that for $`w1`$ the above general bias solution tends to the $`\mathrm{\Lambda }`$CDM case, as it should (eq.39 of Basilakos & Plionis 2001). It is quite instructive to see the solution for the case of $`w=1/3`$, in which case the functional form of the bias evolution is independent of the hyper-geometric function111Where we have used $`\mathrm{\Omega }_x=1\mathrm{\Omega }_\mathrm{m}`$ and $`F[\alpha ,b,b,t]=(1t)^\alpha `$.. We have that:

$$b(z)=𝒜(1+z)(1+\mathrm{\Omega }_\mathrm{m}z)^{1/2}+\frac{2𝒞}{\mathrm{\Omega }_\mathrm{m}}(1+z)+1.$$
(11)
We can easily verify that for $`\mathrm{\Omega }_\mathrm{m}=1`$ this reduces to the Einstein-de Sitter case (eq. 29 of Basilakos & Plionis 2001), as it should.

## 3 Estimating the Bias from Galaxy and Mass correlations

Since our approach gives a family of bias curves, due to the fact that it has two unknown parameters, (the integration constants $`𝒜,𝒞`$) and in order to obtain partial solutions for $`b(z)`$ we need to estimate the values of these constants. In Basilakos & Plionis (2001) we have compared our bias evolution model with the halo merging model (cf. Mo & White 1996; Mataresse et al. 1997) as well as with different N-body results and found a very good consistency, once we fitted the integration constants $`𝒜,𝒞`$ by evaluating our model to two different epochs. We further compare in figure 1 our solution for the $`\mathrm{\Lambda }`$CDM case (see its parameters further below), evaluated at $`z=0`$ and $`z=3`$ using the HDF results (Arnouts et al. 2002; Malioccietti 1999, with the Mataresse et al (1997) model. It is quite evident that our model fits better the $`z`$-dependance of the observational HDF galaxy bias.

Although, this comparison gives consistent results of the functional form of our solution with available data and theoretical models (see Basilakos & Plionis 2001), it does not test directly whether our model, once calibrated observationally, provides a consistent model at high $`z`$’s. Therefore, in order to test it, we evaluate the constants $`𝒜,𝒞`$ by fitting our bias model to the the recent 2dF galaxy clustering results of Hawkins et al. (2002) who used $``$200000 galaxies observed in the 2dF survey to derive their spatial correlation function.

We use the standard theoretical approach to estimate the two point spatial correlation function, using our model for the bias evolution in the different spatially flat cosmological models. We quantify the evolution of clustering with epoch, writing the spatial galaxy correlation function as $`\xi _{\mathrm{model}}(r,z)=\xi _{\mathrm{mass}}(r)R(z)`$, with $`R(z)=D^2(z)b^2(z)`$, where the function $`R(z)`$ characterizes the clustering evolution with epoch (see Basilakos 2001 and references therein). While the $`\xi _{\mathrm{mass}}(r)`$ is the Fourier transform of the spatial power spectrum $`P(k)`$:

$$\xi _{\mathrm{mass}}(r)=\frac{1}{2\pi ^2}_0^{\mathrm{}}k^2P(k)\frac{\mathrm{sin}(kr)}{kr}dk,$$
(12)
where $`k`$ is the comoving wavenumber. Note that we also use the non-linear corrections introduced by Peacock & Dodds (1994). As for the power spectrum of our CDM models, we take $`P(k)k^nT^2(k)`$ with scale-invariant ($`n=1`$) primeval inflationary fluctuations and $`T(k)`$ the CDM transfer function. In particular, we use the transfer function parameterization as in Bardeen et al. (1986), with the corrections given approximately by Sugiyama’s formula (Sugiyama 1995).

In the present analysis we consider flat models with cosmological parameters that fit the majority of observations, ie., $`\mathrm{\Omega }_\mathrm{m}+\mathrm{\Omega }_{x_o}=1`$, $`\mathrm{\Omega }_\mathrm{m}=0.3`$, $`H_{}=100h`$km s-1 Mpc-1 with $`h0.60.7`$ (Freedman et al. 2001; Peebles and Ratra 2002 and references therein) and baryonic density parameter $`\mathrm{\Omega }_\mathrm{b}h^20.02`$ (cf. Olive, Steigman & Walker 2000; Kirkman et al 2003). We investigate such spatially flat cosmological models with negative pressure for different values of $`w`$ (see eq. 3). To this end, all the cosmological models are normalized to have fluctuation amplitude in 8 $`h^1`$Mpc scale of $`\sigma _8=0.50(\pm 0.1)\mathrm{\Omega }_\mathrm{m}^\gamma `$ (Wang & Steinhardt 1998) with $`\gamma =0.210.22w+0.33\mathrm{\Omega }_\mathrm{m}`$.

In order to quantify the integration constants $`(𝒜,𝒞)`$ we perform a standard $`\chi ^2`$ minimization procedure between the measured correlation function for the 2dF galaxies (Hawkins et al. 2002) with those expected in our spatially flat cosmological models,

$$\chi ^2=\underset{i=1}{\overset{n}{}}\left[\frac{\xi _{2df}^i(r)\xi _{\mathrm{model}}^i(r)}{\sigma ^i}\right]^2.$$
(13)
where $`\sigma `$ is the observed correlation function uncertainty.

In figure 2 we present the fit to the data of two models with values of $`w`$ at the opposite end of the range explored. Hawkins et al (2002) found that the 2dF clustering behaviour on small scales is represented by a power law, with correlation length $`r_{}5.05h^1`$Mpc and slope $`\gamma 1.67`$. Due to the interplay between the three unknown ($`𝒜,𝒞,b_{}`$) all of the models (for different values of the constants) fit at a high significance level ($`P_{\chi ^2}0.7`$) the 2dF correlation function. In Table 1 we list the results of the fits for all models, ie., the integration constants and the value of the optical bias, $`b_{}`$, at the present time, derived from our model (eq. 10). It is interesting that for the standard $`\mathrm{\Lambda }`$CDM model ($`w=1`$, $`\mathrm{\Omega }_m=0.3`$ and $`h=0.7`$) the $`𝒞(1+z)`$ factor of the solution is $`0`$, which means that the bias evolution is a very weak function of $`z`$, approximated by: $`b(z)1\frac{1}{33}(1+z)^{3/2}`$.

It is also very interesting that for the $`\mathrm{\Lambda }`$CDM and the $`w=0.9`$ QCDM models (with $`h=0.7`$), the corresponding bias at the present time is $`b_{}1.05`$ and $`b_{}1.11`$, respectively, in very good agreement with the values $`b_{}1.04`$ and 1.11 derived by Verde et al. (2002) and Lahav et al. (2002), respectively. When we use $`h=0.6`$, only the $`\mathrm{\Lambda }`$CDM model provides a value for the present bias which is only roughly consistent with other observational results ($`b_{}1.17`$). The present time bias, $`b_{}`$, as a function of $`w`$, is well fitted by a power law having the form: $`b_{}[1.931.20h](\pm 0.019)\times (w)^{0.411(\pm 0.018)}`$.

Having derived the values of the constants $`𝒜,𝒞`$ and the zero-point bias, $`b_{}`$, from the fit to the 2dF galaxy clustering pattern, we can now investigate the bias evolution in the different cosmological models. In figure 3 we present the function $`b(z)`$ for different values of $`w`$. The biasing is a monotonically increasing function of redshift with its evolution being significantly stronger for lower values of $`w`$. Galaxy clustering in quintessence models evolves more rapidly than in the $`\mathrm{\Lambda }`$CDM model. It is also clear that for the preferred cosmological models (with $`w0.9`$ and $`h=0.7`$) the optical galaxies are only weakly biased, even at high redshifts.

In figure 3 we also compare our evolution model with the galaxy conserving one (Fry 1996; Tegmark & Peebles 1998) by normalizing the latter to the value of $`b_{}`$ provided by our fit. It is interesting that the two models, despite their different form, are consistent with each other. However, we note that this is true for the $`h=0.7`$ case. Had we presented the $`h=0.6`$ case it would have be evident that the consistency is significantly worse. In the same figure we compare our model with the results of a $`\mathrm{\Lambda }`$CDM N-body simulation of Somerville et al (2001) (thick blue line). The simulation galaxy bias evolves significantly more than what our model predicts, which should be attributed to our assumption that the galaxy number density is conserved in time. It is evident that merging processes, not taken into account by our model, are very important in the evolution of clustering.

## 4 Conclusions

In this paper we derived analytically the evolution of linear bias $`b(z)`$ for cosmological models driven by an exotic fluid with an equation of state: $`p_x=w\rho _x`$, where $`1w<0`$. Comparing the optical 2dF galaxy correlation function (Hawkings et al. 2002) with the predictions of various (quintessence) models, we find that the flat cosmological models with $`w0.9`$ provide present bias values in the range $`1.05b_{}1.11`$, which is consistent with observational correlation results. However, we find that for the above cosmological models our bias evolution is too weak to fit N-body results of galaxy clustering, a fact which should be attributed to merging processes, not taken into account in our model.