text:
# A CLASSIFICATION OF SPHERICALLY SYMMETRIC SELF-SIMILAR DUST MODELS

## 1 Introduction

Spherically symmetric self-similar solutions to Einstein’s equations have the feature that every dimensionless variable is a function of some dimensionless combination of the cosmic time coordinate $`t`$ and the comoving radial coordinate $`r`$. In the simplest situation, a similarity solution is invariant under the transformation $`rat,tat`$ for any constant $`a`$, so the similarity variable is $`z=r/t`$. Geometrically this corresponds to the existence of a homothetic Killing vector. Such solutions have been the focus of much attention in General Relativity because the field equations simplify to ordinary differential equations (Cahill & Taub 1971).

Even greater simplification is afforded if one focusses on the situation in which the source of the gravitational field is pressureless dust since, in this case, the solutions can often be expressed analytically and are just a special subclass of the more general spherically symmetric Tolman-Bondi solutions (Tolman 1934, Bondi 1947, Bonnor 1956). A number of people have studied such solutions: for example, Gurovich (1967), Carr & Hawking (1974), Dyer (1979), Wu (1981), Maharaj (1988), Carr & Yahil (1990), Ori & Piran (1990), Henriksen & Patel (1991) and Sintes (1996). These papers are described in some detail by Krasinski (1997). As examples of applications of particular physical intererst, self-similar dust solutions have also been used in studying naked singularities by Eardley & Smarr (1984), Christodoulou (1984), Lake (1992), Joshi & Dwivedi (1993), Joshi & Singh (1995), Dwivedi & Joshi (1997) and in modelling cosmic voids by Tomita (1995, 1997a, 1997b) and Carr & Whinnett (1999).

In this paper we will present a classification of spherically symmetric homothetic dust models which is “complete” subject to certain specified restrictions. This serves as a first step in the more general analysis, presented in an accompanying paper (Carr & Coley 1999a), of spherically symmetric self-similar perfect fluid models with equation of state $`p=\alpha \mu `$. Despite the simplifications entailed in dropping pressure, we will find that many of the features of the dust ($`\alpha =0`$) solutions carry over to the $`\alpha 0`$ case, at least in the supersonic regime. In particular, many (though not all) of the types of solutions with pressure have direct analogues in the dust case. Our aim is therefore to use the exact analytic dust solutions to derive qualitative features of self-similar solutions which will also turn out to pertain when there is pressure but which can then only be demonstrated numerically. It should, of course, be stressed that the introduction of pressure is also associated with many new features - especially in the subsonic regime - which cannot be understood in this way.

The plan of the paper is as follows: In Section 2 we will show how the dust equations can be regarded as a special case of the equations with pressure providing one adopts a specific prescription in taking the limit $`\alpha 0`$. In Section 3 we discuss the general family of self-similar solutions, showing that they can be conveniently categorized as asymptotically Friedmann and what we term asymptotically quasi-static. In each case, we will focus on the important physical features of the solutions, such as the presence of apparent horizons, event horizons and singularities. We will also clarify the connection between models with positive and negative $`z`$. We will specify the sense in which our classification is “complete” in Section 4.

## 2 Spherically Symmetric Similarity Solutions

In the spherically symmetric situation one can introduce a time coordinate $`t`$ such that surfaces of constant $`t`$ are orthogonal to fluid flow lines and comoving coordinates ($`r,\theta ,\varphi `$) which are constant along each flow line. The metric can then be written in the form

$$ds^2=e^{2\nu }dt^2e^{2\lambda }dr^2R^2d\mathrm{\Omega }^2,d\mathrm{\Omega }^2d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2$$
(2.1)
where $`\nu `$, $`\lambda `$ and $`R`$ are functions of $`r`$ and $`t`$. The Einstein equations have a first integral

$$m(r)=\frac{1}{2}R\left[1+e^{2\nu }\left(\frac{R}{t}\right)^2e^{2\lambda }\left(\frac{R}{r}\right)^2\right].$$
(2.2)
This can be interpreted as the mass within comoving radius $`r`$ at time $`t`$:

$$m(r)=4\pi _0^r\mu R^2\frac{R}{r^{}}𝑑r^{}.$$
(2.3)
where $`\mu (r,t)`$ is the energy density and we choose units in which $`c=G=1`$. This is constant if $`p=0`$. Eqn (2.2) can be written as an equation for the energy per unit mass of the shell with comoving coordinate $`r`$:

$$E\frac{1}{2}(\mathrm{\Gamma }^21)=\frac{1}{2}U^2\frac{m}{R},Ue^\nu \left(\frac{R}{t}\right),\mathrm{\Gamma }e^\lambda \left(\frac{R}{r}\right).$$
(2.4)
This can be interpreted as the sum of the kinetic and potential energies per unit mass. $`E`$ and $`\mathrm{\Gamma }`$ are conserved along fluid flow lines in the $`p=0`$ case.

By a spherically symmetric similarity solution we shall mean one in which the spacetime admits a homothetic Killing vector $`𝝃`$ that satisfies

$$\xi _{\mu ;\nu }+\xi _{\nu ;\mu }=2g_{\mu \nu }.$$
(2.5)
This means that the solution is unchanged by a transformation of the form $`tat`$, $`rar`$ for any constant $`a`$. Solutions of this sort were first investigated by Cahill & Taub (1971), who showed that by a suitable coordinate transformation they can be put into a form in which all dimensionless quantities such as $`\nu `$, $`\lambda `$, $`E`$ and

$$S\frac{R}{r},M\frac{m}{R},PpR^2,W\mu R^2$$
(2.6)
are functions only of the dimensionless variable $`zr/t`$. This means that the field equations reduce to a set of ordinary differential equations in $`z`$. Another important quantity is the function

$$V(z)=e^{\lambda \nu }z,$$
(2.7)
which represents the velocity of the surfaces of constant $`z`$ relative to the fluid. These surfaces have the equation $`r=zt`$ and therefore represent a family of spheres moving through the fluid. The spheres contract relative to the fluid for $`z<0`$ and expand for $`z>0`$. This is to be distinguished from the velocity of the spheres of constant $`R`$ relative to the fluid:

$$V_R=\frac{U}{\mathrm{\Gamma }}=e^{\lambda \nu }\left(\frac{R/t}{R/r}\right).$$
(2.8)
This is positive if the fluid is collapsing and negative if it is expanding. Special significance is attached to values of $`z`$ for which $`|V|=1`$ and $`|V_R|=1`$. The first corresponds to a Cauchy horizon (either a black hole event horizon or a cosmological particle horizon) and the second to a black hole or cosmological apparent horizon. One can show that the existence of an apparent horizon is also equivalent to the condition $`M=1/2`$.

Although our main focus in this paper is (pressureless) dust solutions, it is elucidating to start by considering the equations for a fluid with pressure. The only barotropic equation of state compatible with the similarity ansatz is one of the form $`p=\alpha \mu `$ ($`1\alpha 1`$). If one introduces a dimensionless function $`x(z)`$ defined by

$$x(z)(4\pi \mu r^2)^{\alpha /(1+\alpha )},$$
(2.9)
then the conservation equations $`T_{;\nu }^{\mu \nu }=0`$ can be integrated to give

$$e^\nu =\beta xz^{2\alpha /(1+\alpha )}$$
(2.10)
$$e^\lambda =\gamma x^{1/\alpha }S^2$$
(2.11)
where $`\beta `$ and $`\gamma `$ are integration constants. The remaining field equations reduce to a set of ordinary differential equations in $`x`$ and $`S`$:

$$\ddot{S}+\dot{S}+\left(\frac{2}{1+\alpha }\frac{\dot{S}}{S}\frac{1}{\alpha }\frac{\dot{x}}{x}\right)[S+(1+\alpha )\dot{S}]=0,$$
(2.12)
$$\left(\frac{2\alpha \gamma ^2}{1+\alpha }\right)S^4+\frac{2}{\beta ^2}\frac{\dot{S}}{S}x^{(22\alpha )/\alpha }z^{(22\alpha )/(1+\alpha )}\gamma ^2S^4\frac{\dot{x}}{x}\left(\frac{V^2}{\alpha }1\right)=(1+\alpha )x^{(1\alpha )/\alpha },$$
(2.13)
$$M=S^2x^{(1+\alpha )/\alpha }\left[1+(1+\alpha )\frac{\dot{S}}{S}\right],$$
(2.14)
$$M=\frac{1}{2}+\frac{1}{2\beta ^2}x^2z^{2(1\alpha )/(1+\alpha )}\dot{S}^2\frac{1}{2}\gamma ^2x^{(2/\alpha )}S^6\left(1+\frac{\dot{S}}{S}\right)^2,$$
(2.15)
where the velocity function is given by

$$V=(\beta \gamma )^1x^{(1\alpha )/\alpha }S^2z^{(1\alpha )/(1+\alpha )}$$
(2.16)
and an overdot denotes $`zd/dz`$. The other velocity function is

$$V_R=\frac{V\dot{S}}{S+\dot{S}}$$
(2.17)
and, from eqns (2.4) and (2.15), the energy function is

$$E=\frac{1}{2}\gamma ^2x^{(2/\alpha )}S^6\left(1+\frac{\dot{S}}{S}\right)^2\frac{1}{2},$$
(2.18)
which necessarily exceeds $`1/2`$.

As discussed by Carr & Yahil (1990), we can best envisage how these equations generate solutions by working in the $`3`$-dimensional $`(x,S,\dot{S})`$ space. At any point in this space, for a fixed value of $`\alpha `$, eqns (2.14) and (2.15) give the value of z; eqn (2.13) then gives the value of $`\dot{x}`$ unless $`|V|=\sqrt{\alpha }`$ and eqn (2.12) gives the value of $`\ddot{S}`$. Thus the equations generate a vector field $`(\dot{x},\dot{S},\ddot{S})`$ and this specifies an integral curve at each point of the 3-dimensional space. Each curve is parametrized by $`z`$ and represents one particular similarity solution. This shows that, for a given equation of state parameter $`\alpha `$, there is a $`2`$-parameter family of spherically symmetric similarity solutions. In general there would be a sonic point, with possible associated discontinuities, at $`|V|=\sqrt{\alpha }`$. This corresponds to crossing a 2-dimensional surface in the solution space but we need not discuss this complication here.

We now focus on the dust solutions. Although eqns (2.10) to (2.16) break down when $`\alpha =0`$, in that some of the terms appearing there disappear or diverge, we will show that they are still formally applicable providing the function $`x`$ defined by eqn (2.9) is set to $`1`$ whenever it does not appear with the exponent $`1/\alpha `$. Otherwise one must make the substitution

$$x^{1/\alpha }(4\pi \mu r^2)^1,\frac{1}{\alpha }\frac{\dot{x}}{x}\frac{dln(\mu r^2)}{dlnz},$$
(2.19)
as suggested by eqn (2.9). In fact, the relevant equations are most simply obtained by noting that both the energy and mass within comoving radius $`r`$ are conserved, so that $`E`$ and $`m/r=MS`$ are constant. If we put $`m/r=\kappa `$, then eqn (2.3) implies

$$4\pi \mu R^2\frac{R}{r}=\frac{dm}{dr}=\kappa $$
(2.20)
and this can be combined with eqns (2.9) and (2.11) to give

$$e^\lambda =\gamma ^1(4\pi \mu r^2S^2)^1=\kappa ^1\gamma ^1\frac{R}{r}.$$
(2.21)
On the other hand, eqn (2.4) implies $`e^\lambda =\mathrm{\Gamma }^1(R/r)`$ where $`\mathrm{\Gamma }=\pm \sqrt{1+2E}`$ and so the constant $`\kappa `$ is just $`\mathrm{\Gamma }/\gamma `$. The mass function is therefore

$$M=\frac{\mathrm{\Gamma }}{\gamma S}=\frac{\sqrt{1+2E}}{\gamma S},$$
(2.22)
where we have taken the positive square root for $`\mathrm{\Gamma }`$ to ensure that the mass if positive. (We discuss the negative mass case later.) Putting $`x=1`$ in eqn (2.10) also gives $`e^\nu =\beta `$, so eqn (2.4) can be rewritten as

$$E=\frac{1}{2\beta ^2}z^4\left[\frac{dS}{dz}\right]^2\frac{\sqrt{1+2E}}{\gamma S}.$$
(2.23)
Once this equation has been integrated to give $`S(z)`$, all the other functions can be obtained, so one has solved the problem completely.

We now show that eqns (2.10) to (2.16) are all formally satisfied if one uses the prescription given by eqn (2.19). Eqn (2.14) can be written as

$$MS=m/R=S^2(4\pi \mu r^2)(S+\dot{S})=4\pi \mu R^2\frac{R}{r}$$
(2.24)
and this is just equivalent to eqn (2.20). Since eqn (2.21) implies

$$(4\pi \mu r^2)^1=\gamma S^2(S+\dot{S})/\mathrm{\Gamma },$$
(2.25)
we also have

$$\frac{dln(\mu r^2)}{dlnz}=\frac{2\dot{S}}{S}\frac{\dot{S}+\ddot{S}}{S+\dot{S}}.$$
(2.26)
Eqn (2.12) is then automatically satisfied and from eqn (2.23) formally corresponds to $`\dot{E}=0`$. Finally we can substitute for the $`x`$ terms in eqn (2.13), using eqns (2.19) and (2.25), to obtain

$$\ddot{S}+\dot{S}=\frac{\mathrm{\Gamma }}{S^2z^2}$$
(2.27)
and this can be integrated to give eqn (2.23).

It is now convenient to scale the $`r`$ and $`t`$ coordinates so that $`\beta =\gamma =1`$. Eqn (2.23) then implies

$$\frac{dS}{dz}=\pm \frac{\sqrt{2E+2\mathrm{\Gamma }/S}}{z^2}$$
(2.28)
and this can be integrated to give

$$D\frac{1}{z}=\{\begin{array}{cc}\frac{\sqrt{ES^2+\mathrm{\Gamma }S}}{\sqrt{2}E}\frac{2\mathrm{\Gamma }}{(2E)^{3/2}}\mathrm{sinh}^1\sqrt{\frac{ES}{\mathrm{\Gamma }}}\hfill & (E>0)\hfill \\ & \\ \frac{\sqrt{2}}{3}S^{3/2}\hfill & (E=0)\hfill \\ & \\ \frac{2\mathrm{\Gamma }}{(2E)^{3/2}}\mathrm{sin}^1\sqrt{\frac{ES}{\mathrm{\Gamma }}}\pm \frac{\sqrt{ES^2+\mathrm{\Gamma }S}}{\sqrt{2}E}\hfill & (1/2<E<0)\hfill \end{array}$$
(2.29)
where D is an integration constant and $`\mathrm{sin}^1`$ is taken to lie between $`0`$ and $`\pi `$. In the first two cases, the upper and lower signs on the left apply for $`dS/dz`$ positive and negative, respectively, this sign being constant for any particular solution. In the third case, the sign on the left is fixed for a particular solution, even though $`dS/dz`$ may switch sign, but the sign of the last term is plus if $`\mathrm{sin}^1`$ lies between $`0`$ and $`\pi /2`$ and minus if it lies between $`\pi /2`$ and $`\pi `$. If we took the negative square root in eqn (2.22), corresponding to $`M`$ and $`\mathrm{\Gamma }`$ being negative, there would be another solution for $`E>0`$ given by

$$D\frac{1}{z}=\frac{\sqrt{ES^2|\mathrm{\Gamma }|S}}{\sqrt{2}E}+\frac{2|\mathrm{\Gamma }|}{(2E)^{3/2}}\mathrm{cosh}^1\sqrt{\frac{ES}{|\mathrm{\Gamma }|}}(E>0).$$
(2.30)
This solution is unphysical, since the mass is negative, but it is of interest for comparison with the solutions with pressure.

Eqns (2.16), (2.17), (2.19), (2.25) and (2.28) give the velocity functions as

$$V=\frac{Sz\pm \sqrt{2E+2\mathrm{\Gamma }/S}}{\mathrm{\Gamma }}$$
(2.31)
and

$$V_R=\pm \frac{\sqrt{2E+2\mathrm{\Gamma }/S}}{\mathrm{\Gamma }},$$
(2.32)
while eqns (2.16) and (2.31) imply that the density is given by

$$4\pi \mu t^2=\frac{1}{zS^2V}=\frac{\mathrm{\Gamma }}{zS^2(Sz\pm \sqrt{2E+2\mathrm{\Gamma }/S})},$$
(2.33)
where the upper and lower signs again correspond to $`dS/dz`$ being positive and negative, respectively. Note that V and $`\mu `$ is negative (corresponding to tachyonic models) for the solution given by eqn (2.30) and $`\mu `$ is also negative. In all cases the metric can be written as

$$ds^2=dt^2\frac{(S+\dot{S})^2}{1+2E}dr^2r^2S^2d\mathrm{\Omega }^2,$$
(2.34)
which is the standard Tolman-Bondi form with constant energy function $`E(r)`$.

## 3 Classification of Solutions

Eqn (2.29) implies that there is a 2-parameter family of similarity solutions (as in the general $`\alpha `$ case). In this section we will provide a complete description of these solutions. We will start by considering the simplest one: the flat Friedmann solution with $`D=E=0`$. We will then consider the one-parameter family of solutions with $`E0,D=0`$. Finally we will consider the full two-parameter family of solutions with $`E0,D0`$. In each case, we will show the form of the physically interesting quantities $`S`$, $`V`$ and $`\mu t^2`$ as functions of $`z`$. In obtaining the full family of solutions, it is crucial that we allow $`z`$ to be either positive or negative. Our analysis will also cover the (presumably unphysical) solutions with negative mass because they relate to some of the solutions with pressure.

### 3.1 $`E=D=0`$ solution

In this case eqns (2.21), (2.29), (2.31) and (2.33) give

$$S=(\sqrt{2}z/3)^{2/3}=M^1,V=(z/6)^{1/3},\mu =(6\pi t^2)^1.$$
(3.1)
This corresponds to the standard dust Friedmann model with zero curvature constant. The metric can be put in the usual form by making the substitution $`\widehat{r}=(9r/2)^{1/3}`$, which gives

$$ds^2=dt^2t^{4/3}[d\widehat{r}^2+\widehat{r}^2d\mathrm{\Omega }^2].$$
(3.2)
Note that the curvature constant must be zero because otherwise there would be an intrinsic scale, which would contradict the similarity assumption.

### 3.2 $`D=0`$ solutions

Solutions with $`D=0`$ are asymptotically Friedmann as $`|z|\mathrm{}`$ and are specified entirely by the energy parameter $`E`$ and were studied by Carr & Hawking (1974) and Carr & Yahil (1990). The form of $`S(z)`$ in these solutions is shown in Figure (1a), the arrows always corresponding to the direction of increasing time. The solutions with $`z>0`$ correspond to initially expanding Big Bang models: they start from an initial Big Bang singularity ($`S=0`$) at $`t=0`$ ($`z=\mathrm{}`$) and then either expand indefinitely ($`S\mathrm{}`$) as $`t\mathrm{}`$ ($`z0`$) for $`E0`$ or recollapse to a black hole singularity ($`S=0`$) at

$$z_S=\frac{(2E)^{3/2}}{2\pi \sqrt{1+2E}}$$
(3.3)
for $`E<0`$. Note that $`z_S`$ corresponds to the physical origin since $`R=rS=0`$ there. The form of $`V(z)`$ in the $`z>0`$ solutions is shown in Figure (1b). In the first case, $`V`$ decreases monotonically from $`\mathrm{}`$ to 0. In the second case, it reaches a minimum before rising to $`\mathrm{}`$ at $`z_S`$. One can show that the values of $`z`$ and $`V`$ at the minimum both decrease as $`E`$ increases; the minimum will exceed 1 (in which case the whole Universe is inside the black hole) if $`E`$ is less than some critical negative value $`E_{}`$ and it will be less than 1 (in which case there is a black hole event horizon and a cosmological particle horizon) if $`E`$ exceeds $`E_{}`$. The solutions with $`z<0`$ are the time-reverse of the $`z>0`$ ones and the sign of $`V`$ is also reversed: as $`t`$ increases from $`\mathrm{}`$ to 0 (i.e. as $`z`$ decreases from 0 to $`\mathrm{}`$), the $`E0`$ models collapse from an infinitely dispersed state ($`S=\mathrm{}`$) to a Big Crunch singularity ($`S=0`$); the $`E<0`$ models also collapse to a Big Crunch singularity but they emerge from a white hole and are never infinitely dispersed.

Both $`S`$ and $`V`$ have the same $`z`$-dependence as in the $`E=0`$ Friedmann solution as $`|z|\mathrm{}`$:

$$S[9\sqrt{1+2E}/2]^{1/3}|z|^{2/3},V[6(1+2E)]^{1/3}z^{1/3}.$$
(3.4)
However, the $`E0`$ solutions deviate from the $`E=0`$ solution at small values of $`|z|`$. The $`E<0`$ solutions never reach $`z=0`$ at all, while the $`E>0`$ ones have

$$S(2E)^{1/2}|z|^1,V(1+2E)^{1/2}E^1z\mathrm{ln}[(2E)^{3/2}(1+2E)^{1/2}|z|]$$
(3.5)
as $`|z|0`$. The first relation implies that the circumference function $`R(r,t)=Sr`$ is non-zero in limit $`r0`$ unless $`E=0`$ since

$$R(0,t)=\sqrt{2E}t.$$
(3.6)
This means that the “coordinate” origin ($`r=0`$) is an expanding 2-sphere. \[This feature is specific to the dust case and does not arise if there is pressure.\] This has a natural physical interpretation since the forms of $`S`$ and $`V`$ are similar to those in the Kantowski-Sachs (1966) solution, in which all the matter is localized on a shell \[cf. eqns (3.19) and (3.23) in Carr & Coley (1999a)\]. However, we note that there is no exact self-similar Kantowski-Sachs solution in the dust case. To obtain a complete solution, one must therefore match the self-similar solution onto a (non-self-similar) part inside $`R(t,0)`$. In the $`E<0`$ case, we have seen that the physical origin is the black hole singularity $`z_S`$, so only for $`E=0`$ can one identify $`z=0`$ with the physical origin.

The form of the density function $`\mu t^2`$ can be derived from eqn (2.33) and is shown in Figure (1c). For a given fluid element, this specifies the density as a function of time $`\mu (t)`$. For a given time, it also specifies the density profile $`\mu (r)`$ and this illustrates that a non-zero value of $`E`$ necessarily introduces an inhomogeneity into the model. Solutions with $`E>0`$ are everywhere underdense relative to the Friedmann model, with eqns (2.33) and (3.5) implying that $`\mu t^2`$ goes to 0 as $`(\mathrm{ln}|z|)^1`$ as $`z0`$. (This suggests that the interior non-self-similar region should be a vacuum.) Solutions with $`E<0`$ are everywhere overdense relative to Friedmann, with $`\mu `$ diverging at the singularity. Note that eqns (2.33) and (3.4) imply that $`\mu t^2`$ is independent of $`E`$ to 1st order as $`|z|\mathrm{}`$.

The form of the mass function $`M(z)`$ in the $`D=0`$ solutions is not shown explicitly but can be immediately deduced from the expression for $`S`$ since eqn (2.21) gives $`M=\mathrm{\Gamma }/S`$. In the $`E0`$ case, there is always a single point where $`M=1/2`$ and this corresponds to the cosmological apparent horizon. In the $`E<0`$ case, eqn (2.29) implies that $`S`$ has a maximum of $`\mathrm{\Gamma }/|E|`$ and so eqn (2.21) shows that $`M`$ has a minimum of $`|E|`$. Since this is less than 1/2, there are always two points where $`M=1/2`$, one corresponding to the black hole apparent horizon and the other to the cosmological apparent horizon. Note that a black hole’s apparent horizon always lies within or coincides with its event horizon (Hawking & Ellis 1973), which is why the first can exist without the second. Eqns (2.21) and (3.3) imply that the mass associated with this singularity is

$$m_S=(MSz)_St=(2E)^{3/2}t/(2\pi ).$$
(3.7)
It therefore starts off zero when the singularity first forms at $`t=0`$ but then grows as $`t`$. The mass of the black hole is given by a similar formula but with $`z`$ having the value appropriate for the event horizon or apparent horizon. Since the former may not exist, it is more appropriate to use the latter.

Finally, we consider the $`E>0`$ negative-mass solutions given by eqn (2.30). Their form is indicated by the dotted curves in Figures (1). $`S`$, $`V`$ and $`\mu `$ have the same form as in the positive mass solutions for small values of $`|z|`$ except that V and $`\mu `$ reverse their signs. However, the solutions are very different at large values of $`|z|`$ since eqn (2.30) shows that $`S`$ must always exceed $`|\mathrm{\Gamma }|/E`$. Indeed it tends to this value asymptotically, so we have

$$S\sqrt{1+2E}/E,Vz/E,\mu r^2E^3/\sqrt{1+2E}$$
(3.8)
as $`|z|\mathrm{}`$. The form of this solution is closely related to that of the $`\alpha <<1`$ static solution \[cf. eqn (3.29) of Carr & Coley (1999a)\], although there is no static solution in the $`\alpha =0`$ case itself.

### 3.3 $`E=0`$ solutions

We now put $`E=0`$ and consider the effect of introducing a non-zero value for the constant $`D`$. \[In this case, eqn (2.28) does not permit $`\mathrm{\Gamma }<0`$, so there are no negative-mass solutions.\] The form of $`S(z)`$ for the $`D>0`$ solutions is shown in Figure (2a). There are two types of solutions in this case, one expanding and the other collapsing. For the expanding solutions (solid lines), $`S=0`$ at $`z=1/D`$ and so the Big Bang occurs before $`t=0`$ (i.e. it is “advanced”). As $`t`$ increases to 0 (i.e. as $`z`$ decreases to $`\mathrm{}`$), $`S`$ tends to the finite value

$$S_{\mathrm{}}(D)=(3D/\sqrt{2})^{2/3}.$$
(3.9)
As $`t`$ further increases from 0 to $`+\mathrm{}`$ (i.e. as $`z`$ jumps to $`+\mathrm{}`$ and then decreases to 0), $`S`$ increases monotonically to $`\mathrm{}`$. For the contracting solutions (broken lines), $`S`$ starts infinite at $`t=\mathrm{}`$ ($`z=0`$) and then decreases to $`S_{\mathrm{}}(D)`$ as $`t`$ increases to 0 ($`z\mathrm{}`$). As $`t`$ further increases (i.e. as $`z`$ jumps to $`+\mathrm{}`$ and then decreases), $`S`$ continues to decrease until it reaches 0 at the Big Crunch singularity at $`z=1/D`$. Both types of solutions are characterized by the fact that they have just one singularity and span both positive and negative values of $`z`$. Note that for each value of $`S_{\mathrm{}}`$ the two asymptotic solutions just correspond to the plus and minus signs in eqn (2.28).

The form of $`V(z)`$ in the $`D>0`$ solutions is shown in Figure (2b). For the expanding solutions (solid lines), it starts off at $`\mathrm{}`$ at the Big Bang ($`z=1/D`$), reaches a negative maximum and then, from eqn (2.31), tends to

$$V=S_{\mathrm{}}(D)z=(3D/\sqrt{2})^{2/3}z$$
(3.10)
as $`z\mathrm{}`$. When z jumps $`+\mathrm{}`$, $`V`$ becomes positive but eqn (3.10) still applies. As $`z`$ decreases from $`+\mathrm{}`$ to 0, $`V`$ decreases monotonically to 0. For the contracting models (broken lines), $`V`$ starts from zero at $`z=0`$ and monotonically decreases as $`z`$ goes to $`\mathrm{}`$, being again given by eqn (3.10) asymptotically. When $`z`$ jumps to $`+\mathrm{}`$, $`V`$ jumps to $`+\mathrm{}`$ and then decreases to a minimum before rising to infinity at the Big Crunch singularity. Note that the maximum value of $`V`$ for the expanding solutions will exceed $`1`$ and the minimum value for the contracting ones will be less than $`+1`$ (i.e. the minimum value of $`|V|`$ is less than $`1`$) if $`D`$ exceeds some critical value $`D_+`$. A simple calculation (see later) shows that the values of $`|z|`$ and $`|V|`$ at the stationary point are given by

$$|z|_{min}=\left(\frac{2+\sqrt{3}}{3D}\right),|V|_{min}=\left(\frac{26+15\sqrt{3}}{3D}\right)^{1/3},$$
(3.11)
so the stationary point moves towards the origin as $`D`$ increases and $`D_+=26/3+5\sqrt{3}17`$. For $`D>D_+`$, the condition $`|V|=1`$ will be satisfied at three values of $`z`$. As illustrated by Figure 14 of Ori & Piran (1990), this means that the contracting solutions will form a black hole in which the central singularity is naked for a while.

The crucial feature of these solutions is that, while the form of $`V(z)`$ is like that in the $`D=0`$ case for small $`|z|`$, $`V`$ scales as $`z`$ rather than $`z^{1/3}`$ \[cf. eqn (3.4)\] for large $`|z|`$. This is because any solution with finite $`S`$ at infinity must be “nearly” static in the sense that $`dS/dz`$ tends to zero. However, the solutions are not asymptotic to an exact static solution (indeed this does not exist in the $`\alpha =0`$ case) because eqn (2.32) implies that $`V_R`$ tends to a non-zero value:

$$V_R^{\mathrm{}}=\pm \left(\frac{4}{3D}\right)^{1/3}.$$
(3.12)
We therefore term these solutions asymptotically “quasi-static”. If $`V_R^{\mathrm{}}`$ is positive, the fluid is collapsing at infinity; if $`V_R^{\mathrm{}}`$ is negative, it is expanding. Note that eqn (2.28) implies that both $`dS/dz`$ and $`zdS/dz`$ tend to zero at large $`|z|`$ but $`z^2dS/dz`$ \[which directly relates to $`V_R^{\mathrm{}}`$ from eqns (2.28) and (2.32)\] tends to a non-zero value except in the limit $`D\mathrm{}`$.

The form of the density function $`\mu t^2`$ in the $`D>0`$ solutions is also interesting and is illustrated in Figure (3c). From eqn (2.33) the density parameter is given by

$$\mathrm{\Omega }6\pi \mu t^2=\frac{1}{(13Dz)(1Dz)}$$
(3.13)
where the upper and lower signs apply for positive and negative values of $`dS/dz`$, respectively. \[The inverse of the factor $`6\pi t^2`$ corresponds to the density in a flat Friedmann dust universe, as indicated by eqn (3.1).\] For a given fluid element, this describes how the density evolves as a function of time and it has the expected form. However, at a given time it also prescribes the density profile and one sees immediately that a non-zero value of $`D`$ (like a non-zero value of $`E`$) introduces an inhomogeneity. This inhomogeneity has a particularly interesting form. In the $`z<0`$ regime, the profile for the collapsing solutions is homogeneous for $`|z|<<1/D`$ but has $`\mu r^2`$ for $`|z|>>1/D`$ (i.e. it resembles an isothermal sphere with a uniform core). In the $`z>0`$ regime, the collapsing solutions again have $`\mu r^2`$ for $`|z|>>1/D`$ but the density diverges at $`z=1/D`$ (i.e. one has a density singularity at the centre of an isothermal sphere). For the expanding solutions, the signs of $`z`$ are reversed. These features are illustrated in Figure (2c) and have an obvious physical interpretation.

Although the asymptotically quasi-static solutions have an natural cosmological interpretation when $`z`$ is allowed to span both positive and negative values, we see that the $`z>0`$ and $`z<0`$ solutions also have a non-cosmological interpretation when considered separately: they just represent collapsing and expanding self-similar models which evolve from an initially isothermal distribution. It is interesting that the isothermal model (which is usually associated with a static solution) features prominently in both regimes, despite the fact that there is no exact static solution in the dust case. Note that the behaviour at $`z=0`$ is different from the $`D=0`$ case, in that $`\mathrm{\Omega }`$ tends to $`1`$ rather than $`0`$ for $`E>0`$.

The mass function is $`M=S^1`$ in this case and therefore decreases or increases monotonically. There is just one value of $`z`$ at which $`M=1/2`$ (corresponding to a black hole or cosmological apparent horizon) but this may be in either the positive or negative $`z`$ region. Eqn (3.9) implies that the asymptotic value of $`M`$ as $`|z|\mathrm{}`$ is less than $`1/2`$ for $`D>4/3`$. In this case, the collapsing solutions have their apparent horizon in $`z>0`$, whereas the expanding ones have it in $`z<0`$. The mass of the (possibly naked) singularity in these solutions is

$$m_S=(MSz)_St=t/D$$
(3.14)
from eqn (2.21). As in the asymptotically Friedmann case, it starts off zero at $`t=0`$ but then grows as $`t`$. In the limit $`D\mathrm{}`$, one gets a naked singularity of zero mass at the origin (cf. the static solution).

Finally, we consider the $`D<0`$ solutions. The form of $`S(z)`$ in this case is shown by the dotted curves in Figure (2a). Such solutions are confined to $`|z|<1/D`$, with $`S`$ either decreasing monotonically for $`z<0`$ (i.e. as $`t`$ increases from $`\mathrm{}`$) or increasing monotonically for $`z>0`$ (i.e. as $`t`$ increases to $`+\mathrm{}`$). However, these solutions break down when $`S`$ is too small. This is because eqn (2.31) implies that $`|V|`$ increases to some maximum value and then falls to zero at $`|z|=1/(3D)`$; this is indicated by the dotted curve in Figure (2b). From eqn (2.33) this means that the density diverges there, as shown by the dotted curve in Figure (2c). This divergence is associated with the formation of a shell-crossing singularity since the model resembles the Kantowski-Sachs solution at this point. For $`1/D>|z|>1/(3D)`$, the density and velocity functions become negative but this is presumably unphysical.

### 3.4 $`D0,E0`$ solutions

The forms of $`S(z)`$ for the $`(D>0,E0)`$ solutions are indicated in Figure (3a). The figure assumes that $`D`$ is fixed but allows $`E`$ to vary. The $`(D>0,E>0)`$ solutions are qualitatively similar to the $`(D>0,E=0)`$ ones in that they are monotonically expanding or collapsing and span both positive and negative $`z`$, as illustrated by the upper solid and broken curves, respectively. They are also asymptotically quasi-static, in the sense that $`S`$ tends to a finite value as $`|z|\mathrm{}`$, even though $`V_R`$ is non-zero there from eqn (2.32). The form of the solutions near $`z=0`$ is still given by eqns (3.5) for $`E>0`$, so the behaviour is like that in the $`(E>0,D=0)`$ case here. In particular, $`z=0`$ no longer corresponds to the physical origin, so one again has to attach the solution to a non-self-similar central region.

The $`(D>0,E<0)`$ solutions are qualitatively different from the $`(D>0,E=0`$) ones in that the models no longer collapse from or expand to infinity. This is clear from eqn (2.28), which implies that $`S`$ has a maximum value of $`\mathrm{\Gamma }/|E|`$, so all the solutions start off expanding and then recollapse. Note that there is no exact static solution since that would be incompatible with eqn (2.13), the term on the right-hand-side being non-zero for $`\alpha =0`$. Figure (3a) shows that there are two types of $`(D>0,E<0)`$ solutions. One type (illustrated by the lower solid curves) expands from the Big Bang singularity at $`z=1/D`$ and then recollapses to a black hole singularity at

$$z_S=\left[\frac{2\pi \sqrt{1+2E}}{(2E)^{3/2}}D\right]^1.$$
(3.15)
This reduces to the value given by eqn (3.3) if $`D=0`$. The other type (illustrated by the lower broken curves) expands from a white hole singularity at $`z_S`$ and then recollapses to a Big Crunch singularity at $`z=1/D`$.

Eqn (3.15) implies that $`z_S=1/D`$, so that the solution is symmetric in $`z`$, if $`D`$ has the value

$$D_{sym}\frac{\pi \sqrt{1+2E}}{(2E)^{3/2}}.$$
(3.16)
One can invert this condition to obtain the associated value of $`E`$ in terms of $`D`$:

$$E_{sym}\frac{4\pi }{\sqrt{3}D}\mathrm{sinh}\left[\frac{1}{3}\mathrm{sinh}^1\left(\frac{3\sqrt{3}D}{8\pi }\right)\right]$$
(3.17)
and this specifies a $`1`$-parameter family of solutions with $`V_R^{\mathrm{}}=0`$. If one considers the limit of the symmetric solution as $`D0`$, one finds $`E_{sym}1/2`$ and $`z_S\mathrm{}`$. On the other hand, if one considers the limit as $`D\mathrm{}`$, one finds $`E_{sym}0`$ and $`z_S0`$, so that both singularities go the origin. This is the closest one can get to a static solution in the dust case. Note also that $`z_S0`$ in the limit $`E0^{}`$ whatever the value of $`D`$; the sudden transition as one goes from $`E=0^{}`$ to $`E=0^+`$ is illustrated in Figure (3a).

This value of $`E`$ given by eqn (3.17) has a special physical significance in that it prescribes the minimum value of $`E`$ allowed for given $`D`$, as indicated by the lower boundary in Figure (4). The proof of this is as follows. If one takes the limit of eqns (2.14) and (2.15) with $`\alpha =0`$ as $`z\mathrm{}`$, using eqns (2.19) and (2.32) and the fact that $`\dot{S}0`$ from eqn (2.28), one obtains

$$z^2\dot{S}^2=(1+2E)(V_R^{\mathrm{}})^2=(4\pi \mu r^2)_{\mathrm{}}^2S_{\mathrm{}}^6+(8\pi \mu r^2)_{\mathrm{}}S_{\mathrm{}}^21$$
(3.18)
where one can regard $`(\mu r^2)_{\mathrm{}}`$ and $`S_{\mathrm{}}`$ as independent asymptotic parameters. If one now fixes $`(\mu r^2)_{\mathrm{}}`$ and assumes that it is positive, then the right-hand-side of eqn (3.18) decreases monotonically with decreasing $`S_{\mathrm{}}`$. One therefore gets a real solution for $`V_R^{\mathrm{}}`$ only if $`S_{\mathrm{}}`$ exceeds a certain value and - by monotonicity - this must be the value associated with the symmetric solution. Thus a real solution (with positive density) exists only for $`E>E_{sym}`$ or, equivalently, $`D<D_{sym}`$. This means that $`z_S`$ is always positive and less than $`1/D`$ and that the maximum of $`S`$ will always occur at the opposite sign of $`z`$ as the $`|z|=1/D`$ singularity. These features are illustrated by the curves in Figure (3a).

The form of $`V(z)`$ in the these solutions is shown in Figure (3b). For $`(D>0,E>0)`$ it is similar to that in the $`(D>0,E=0)`$ case. As $`z`$ decreases from $`0`$ to $`\mathrm{}`$, $`V`$ decreases monotonically from $`0`$ to $`\mathrm{}`$; it then jumps to $`z=+\mathrm{}`$ and, as $`z`$ continues to decrease, it falls to a minimum and rises to infinity at the Big Crunch singularity at $`z=1/D`$. Also as in the $`(D>0,E=0)`$ case, this minimum will fall below 1, corresponding to a naked singularity, providing $`D`$ exceeds some value $`D_+(E)`$. We derive an implicit expression for $`D_+(E)`$ later but, for the present, note that it increases with increasing $`E`$ and reduces to the value $`D_+`$ which arose in Section 3.3 when $`E=0`$. The condition for a naked singularity can also be expressed as the requirement that $`E`$ exceed some critical value $`E_+(D)`$.

For $`(D>0,E<0)`$ the form of $`V`$ is similar to that in the $`(D=0,E<0)`$ case. As $`z`$ decreases from $`1/D`$, $`V`$ rises from $`\mathrm{}`$ until some maximum and then falls to $`\mathrm{}`$ quasi-statically as $`z\mathrm{}`$. It then jumps to $`z=+\mathrm{}`$ and falls to a minimum before rising to $`\mathrm{}`$ at $`z_S`$. As in the $`(D=0,E>0)`$ case, there will be a black hole event horizon if the minimum value of $`V`$ is less than $`1`$; this requires that $`E`$ exceed some value $`E_{}(D)`$, which must reduce to the value $`E_{}`$ given in Section 3.2 when $`D=0`$. We derive an implicit expression for $`E_{}(D)`$ below. It should be stressed that the value $`E_{}(D)`$ is associated with the minimum of $`|V|`$ near the singularity at $`|z|=z_S`$ and is different from the value $`E_+(D)`$ associated with the minimum near $`|z|=1/D`$. The relationship between the values $`E_+(D)`$ and $`E_{}(D)`$ is discussed below.

In order to understand the form of the curves in Figures (3a) and (3b) more precisely, it is useful to specify their asymptotic behaviour. As $`|z|\mathrm{}`$, eqns (2.29) and (2.31) with $`E>0`$ imply that $`S(z)`$ and $`V(z)`$ have the following asymptotic forms:

$$S_{\mathrm{}}D\sqrt{2E},V_{\mathrm{}}\left(\frac{2E}{1+2E}\right)^{1/2}Dz$$
(3.19)
for $`D>>[(1+2E)/E^3]^{1/2}`$ and

$$S_{\mathrm{}}\left(\frac{1+2E}{4}\right)^{1/6}(3D)^{2/3},V_{\mathrm{}}\frac{(3D)^{2/3}z}{2^{1/3}(1+2E)^{1/3}}$$
(3.20)
for $`D<<[(1+2E)/E^3]^{1/2}`$. The transition value for $`D`$ between these two regimes is just an extrapolation of the expression for $`D_{sym}`$ given by eqn (3.16) into the $`E>0`$ regime; it scales as $`E^1`$ for $`E>>1`$ and $`E^{3/2}`$ for $`E<<1`$. Note that eqn (3.20) agrees with eqns (3.9) and (3.10) in the limit $`E=0`$. For $`E<0`$, eqn (3.20) still applies if $`D<<D_{sym}`$ but one has

$$S_{\mathrm{}}\sqrt{1+2E}/|E|,V_{\mathrm{}}z/|E|$$
(3.21)
for $`DD_{sym}`$ (i.e. $`S_{\mathrm{}}`$ tends to the value associated with the symmetric solution).

These equations prescribe the asymptotic forms for $`S_{\mathrm{}}`$ and $`V_{\mathrm{}}`$ in the different $`(E,D)`$ regimes of Figure 4. For fixed $`D`$, eqns (3.19) to (3.21) show that $`S_{\mathrm{}}`$ always increases with $`E`$ but is roughly constant for $`|E|<<1`$. The behaviour of $`V_{\mathrm{}}`$ is more complicated: for $`D<<1`$, it first decreases with increasing $`E`$ and then flattens off; for $`D>>1`$, it first increases with increasing $`E`$ before flattening off. These features are indicated in the Figures (3a) and (3b). Note that all these solutions are quasi-static and not exactly static asymptotically since eqn (2.32) gives

$$V_R^{\mathrm{}}\pm \frac{4^{1/3}}{(3D)^{1/3}(1+2E)^{1/3}},V_R^{\mathrm{}}=\pm \left(\frac{2|E|}{1+2E}\right)^{1/2}$$
(3.22)
for $`D>>[(1+2E)/E^3]^{1/2}`$ and $`D<<[(1+2E)/|E|^3]^{1/2}`$. The sign is positive for collapsing solutions and negative for expanding ones. Note that the first expression agrees with eqn (3.12) in the limit $`E=0`$.

We now derive implicit expressions for the functions $`E_{}(D)`$ and $`E_+(D)`$, which are related to the existence of event horizons or naked singularities. Differentiating eqn (2.31) shows that when $`dV/dz=0`$ one always has

$$V=1/(S^2z)$$
(3.23)
and eqn (2.31) then gives

$$V\mathrm{\Gamma }1/(VS)=\pm \sqrt{2E+2\mathrm{\Gamma }/S},$$
(3.24)
where the positive and minus signs corresponds to the sign of $`dS/dz`$. If one also requires $`|V|=1`$ at the stationary point, eqn (3.24) gives two roots

$$S=2\mathrm{\Gamma }\pm \sqrt{4\mathrm{\Gamma }^21},|z|=8\mathrm{\Gamma }^214\mathrm{\Gamma }\sqrt{4\mathrm{\Gamma }^21},$$
(3.25)
where the plus and minus signs are distinct from the ones appearing in eqn (3.24). These roots can be real providing $`\mathrm{\Gamma }>1/2`$, corresponding to $`E>3/8`$. However, one needs to check whether both of these solutions satisfy condition (3.24).

Inserting the solutions (3.25) into eqn (2.29) gives an expression for $`D`$ in terms of $`E`$. This expression is complicated in general but it simplifies in certain regimes. For $`E>>1`$, which implies $`\mathrm{\Gamma }\sqrt{2E}>>1`$, one obtains two possible solutions:

$$S4\mathrm{\Gamma },|z|1/(16\mathrm{\Gamma }^2),D32E$$
(3.26)
and

$$S1/(4\mathrm{\Gamma }),|z|16\mathrm{\Gamma }^2,D\frac{1}{E}\left(\frac{13}{32}\mathrm{ln}\sqrt{2}\right)0.06E^1.$$
(3.27)
The first has $`\mathrm{\Gamma }>1/S`$ and therefore requires $`dS/dz>0`$ from eqn (3.24), which leads to a consistent solution. However, the second has $`\mathrm{\Gamma }<1/S`$ and requires $`dS/dz<0`$, which does not. In the limit $`E0`$, which implies $`\mathrm{\Gamma }1`$, one obtains

$$S(2\pm \sqrt{3})(1\pm 2E/\sqrt{3}),|z|(74\sqrt{3})(14E/\sqrt{3}).$$
(3.28)
The upper sign gives $`\mathrm{\Gamma }>1/S`$ and therefore requires $`dS/dz>0`$, which leads to a consistent solution as $`E`$ tends to $`0`$ from either above or below. Eqn (2.29) then gives

$$D\left(\frac{26+15\sqrt{3}}{3}\right)+E\left(\frac{109+63\sqrt{3}}{6}\right)17+36E.$$
(3.29)
Note that the constant part of this expression gives the same limiting value of $`D`$ as implied by eqn (3.11). The lower sign in eqn (3.28) gives $`\mathrm{\Gamma }<1/S`$ and requires $`dS/dz<0`$, which does not lead to a consistent solution. Finally we note that eqns (3.25) and (2.29) lead to a unique value of $`E`$ and $`D`$ for which the symmetric solution has $`V_{min}=1`$; this necessarily corresponds to a point on the lower boundary in Figure 4.

These limiting behaviours allow one to infer the rough form of the functions $`E_{}(D)`$ and $`E_+(D)`$, as indicated in Figure (4). Here we have used the fact that $`dE_+/dD`$ is positive as $`E0`$, as follows from eqn (3.29). The form of the functions in the $`E<0`$ regime can be inferred from the fact that $`E_{}(D)`$ must reach the value $`E_{}`$ mentioned in Section 3.2 when $`D=0`$ (although this value has not been calculated explicitly). Also $`E_+(D)`$ and $`E_+(D)`$ must reach the line $`E=E_{sym}(D)`$ at the same value of $`E`$ and this must clearly exceed $`3/8`$ from eqn (3.25). We note that, for sufficiently large values of $`D`$, there may be both an event horizon and a naked singularity.

The form of $`\mu t^2`$ in these solutions is shown in Figure (3c), although this gives only some of the solutions shown in Figures (3a) and (3b). It can be understood as a composite of the curves shown in Figure (1c) for $`E<0`$ and Figure (2c) for $`E>0`$. From eqns (2.33) and (3.19) to (3.21), the asymptotic form of the density profile density is given by

$$4\pi \mu r^2[\frac{1+2E}{D^3(2E)^2},\frac{2}{9D^2}]$$
(3.30)
for $`D>>[(1+2E)/E^3]^{1/2}`$ and $`D<<[(1+2E)/|E|^3]^{1/2}`$ \[cf. eqns (3.8) and (3.13).\] The $`E>0`$ solutions are everywhere underdense relative to the $`E=0`$ solutions, going to $`0`$ as $`(\mathrm{ln}|z|)^1`$ at the origin, whereas the $`E<0`$ solutions are everywhere overdense and have a second density singularity. There is a uniform core region only in the $`E=0`$ case, although CC show that this also applies for $`E<0`$ if there is pressure.

The form of $`M(z)=\mathrm{\Gamma }/S(z)`$ can be deduced immediately from Figure (3a). In the $`E>0`$ case it rises or falls monotonically, as in the $`D=0`$ case, so there is just one value of $`z`$ for which $`M=1/2`$. As $`|z|\mathrm{}`$, $`M`$ tends to a limiting value

$$M_{\mathrm{}}=\frac{\sqrt{1+2E}}{S_{\mathrm{}}(D,E)}$$
(3.31)
and the apparent horizon will be in $`z>0`$ or $`z<0`$ according to whether this is greater or less than 1/2. In the $`E<0`$ case, $`M`$ will have a minimum where $`S`$ has a maximum. This occurs at

$$|z|=\left[\frac{\pi \sqrt{1+2E}}{(2E)^{3/2}}D\right]^1$$
(3.32)
and since the minimum value is $`|E|`$, this is necessarily less than $`1/2`$. One therefore has at least two points where $`M=1/2`$, one of which is the apparent horizon for the black hole associated with the singularity at $`z_S`$. Therefore, as in the $`D=0`$ case, there will always be a black hole apparent horizon but not necessarily an event horizon. This emphasizes an important difference between the collapse singularities at $`z=1/D`$ and $`z=z_S`$: only the latter is associated with an apparent horizon, which is why only the former can be naked.

The $`D<0`$ solutions have the same form as in the $`D=0`$ case, except that the $`E<0`$ ones have a second singularity at the value $`z_S`$ given by eqn (3.15) with $`D<0`$. In this case, as $`z`$ goes from $`1/D`$ to $`z_S`$, $`S`$ first increases to some maximum value and then decreases, while $`V`$ monotonically increases from $`\mathrm{}`$ to $`+\mathrm{}`$. As in the $`E=0`$ case, such models are probably physically unrealistic since the density diverges due to shell-crossing. They are therefore not shown explicitly. The form of the (unphysical) negative-mass solutions, which only exist for $`E>0`$, is indicated by the dotted curve in Figures (3). This is similar to the $`D=0`$ case shown in Figures (2) except that asymptotically eqn (3.19) applies rather than eqn (3.8).

## 4 Conclusion

We may briefly summarize the results of our analysis as follows. (1) There are two families of spherically symmetric self-similar dust models: asymptotically flat Friedmann solutions and what we have termed asymptotically quasi-static solutions. These all represent inhomogeneous cosmological models in which the energy function $`E`$ is constant. They either expand from a Big Bang or collapse to a Big Crunch but the singularity is only at $`t=0`$ for the asymptotically Friedmann family. (2) Some of the asymptotically Friedmann models represent overdensities in a Friedmann background which recollapse to a second singularity and contain a black hole which grows as fast as the Universe. The black hole always has an apparent horizon but not necessarily an event horizon. Other asymptotically Friedmann models represent underdensities in a Friedmann background which grow as fast as the Universe. (3) The asymptotically quasi-static models can be interpreted as representing either inhomogenenous cosmological solutions (with one or two singularities) if one allows both signs of $`z`$ or self-similar collapse from an initially isothermal distibution if one allows just one sign of $`z`$, with a uniform density core in one regime and a central black hole or naked singularity in the other. (4) We have emphasized the relationship between the $`z>0`$ and $`z<0`$ solutions. Any particular asymptotically Friedmann solution is confined to one sign of $`z`$ but any asymptotically quasi-static solution necessarily spans both signs.

In an accompanying paper (Carr & Coley 1999a), it is shown that the spherically symmetric self-similar solutions with pressure share many of the qualitative features of the dust ones, especially in the supersonic regime. In particular, all of the properties (1) to (4) above still pertain. However, it should be emphasized that new types of solution arise when there is pressure. For example, there is an exact static solution and an exact Kantowski-Sachs solution, as well as families of solutions asymptotic to these. There are also asymptotically Minkowski solutions for $`\alpha >1/5`$, some of which asymptote to a finite value of $`z`$. The inclusion of pressure obviously introduces qualitatively new features in the subsonic regime, in particular the possible presence of a sonic point.

In claiming that our classification is “complete”, it should be emphasized that our considerations have been confined to similarity solutions of the simplest kind (i.e. homothetic solutions in which the similarity variable is $`zr/t`$). However, it should be noted that this is not the only type of similarity. For example, Carter & Henriksen (1989) have generalized the concept to include what they term “kinematic” self-similarity. In this context the similarity variable is of the form $`z=r/t^a`$ for $`a1`$ and the solution may contain some dimensional constant. Ponce de Leon (1993) has also introduced the closely related notion of “partial homothety”. It is not yet clear how easily the analysis of this paper can be extended to these cases. Finally it should also be emphasized that we have only been studying solutions which are homothetic everywhere. The sort of models considered by Tomita (1997), in which one patches a self-similar transition region between other non-self-similar regions, is clearly not covered here.

Acknowledgments

The author thanks Alan Coley for useful discussions and is grateful to the Yukawa Institute for Theoretical Physics at Kyoto University and the Department of Mathematics and Statistics at Dalhousie University for hospitality received during this work.

References

1. H. Bondi, 1947, MNRAS 107, 410.
2. W. B. Bonnor, 1956, J. Astrophys. 39, 143.
3. A. H. Cahill and M. E. Taub, 1971, Comm. Math. Phys. 21, 1.
4. B. J. Carr and A. A. Coley, 1999a, Phys. Rev. D., in press.
5. B. J. Carr and A. A. Coley, 1999b, Class. Quant. Grav., 16, R31.
6. B. J. Carr and S. W. Hawking, 1974, MNRAS 168, 399.
7. B. J. Carr and A. Yahil, 1990, Ap. J. 360, 330.
8. B. J. Carr and A. Whinnett, 1999, preprint.
9. B. Carter and R. N. Henriksen, 1989, Ann. Phys. Supp. 14, 47.
10. D. Christodoulou, 1984, Commun. Math. Phys. 93, 171.
11. I. H. Dwivedi and P. S. Joshi, 1997, Class. Quant. Grav. 47, 5357.
12. C. C. Dyer, 1979, MNRAS 189, 189.
13. D. M. Eardley and L. Smarr, 1979, Phys. Rev. D. 19, 2239.
14. V. T. Gurovich, 1967, Sov. Phys.-Dokl. 11, 569.
15. S. W. Hawking and G. F. R. Ellis, 1973, The Large-Scale Structure of Spacetime, Cambridge University Press.
16. R. N. Henriksen and K. Patel, 1991, Phys. Red. D. 42, 1068.
17. P. S. Joshi and I. H. Dwivedi, 1993, Phys. Rev. D. 47, 5357.
18. P. S. Joshi and T. P. Singh, 1995, Phys. Rev. D. 51, 6778.
19. R. Kantowski and R. Sachs, 1966, J. Math. Phys. 7, 443.
20. A. Krasinzki, 1997, Physics in an Inhomogeneous Universe, Cambridge University Press.
21. K. Lake, 1992, Phys. Rev. Lett. 68, 3129.
22. R. Larson, 1969, MNRAS 145, 271.
23. S. D. Maharaj, 1988, J. Math. Phys. 29, 1443.
24. M. V. Penston, 1969, MNRAS 144, 449.
25. A. Ori and T. Piran, 1990, Phys. Rev. D. 42, 1068.
26. J. Ponce de Leon, 1993, Gen. Rel. Grav. 25, 865.
27. A. M. Sintes, 1996, PhD thesis (University of Balearic Islands).
28. R. C. Tolman, 1934, Proc. Nat. Acad. Sci. USA 20, 169.
29. K. Tomita, 1995, Ap. J. 451, 1.
30. K. Tomita, 1997a, Phys. Rev. D. 56, 3341.
31. K. Tomita, 1997b, Gen. Rel. Grav. 13, 625.
32. Z. C. Wu, 1981, Gen. Rel. Grav. 13, 625.

Figures

FIGURE (1). This shows the form of (a) the scale factor $`S(z)`$, (b) the velocity function $`V(z)`$ and (c) the density function $`\mu t^2(z)`$ for the asymptotically Friedmann dust models, the arrows indicating the direction of increasing time. These are described by a single parameter $`E`$ where $`E=0`$ in the exact Friedmann case: the $`z>0`$ solutions are overdense and collapse to black holes for $`E<0`$ (with an event horizon for $`E>E_{}`$ since $`V_{min}<1`$) but they are underdense and expand forever for $`E>0`$. The $`z<0`$ solutions are just the time reverse of these. The dotted curve corresponds to a solution with negative mass; it is probably unphysical but relates to the Kantowski-Sachs solution which arises when there is pressure.

FIGURE (2). This shows the form of (a) the scale factor $`S(z)`$, (b) the velocity function $`V(z)`$ and (c) the density function $`\mu t^2(z)`$ for dust models with $`E=0`$. Two different values of $`D`$ are shown in (a) and (b) but only one in (c). These solutions necessarily span both positive and negative values of $`z`$. For $`D>0`$ they represent monotonically expanding (solid) or collapsing (broken) solutions and the latter contain a naked singularity ($`V_{min}<1`$) if $`D`$ exceeds some value $`D_+`$ (as assumed here). The $`D<0`$ models (dotted) undergo shell-crossing before encountering the singularity and are probably unphysical since $`V`$ and $`\mu `$ go negative.

FIGURE (3). This shows the form of (a) the scale factor $`S(z)`$, (b) the velocity function $`V(z)`$ and (c) the density function $`\mu t^2(z)`$ for the asymptotically quasi-static dust solutions. These are described by two parameters ($`D`$ and $`E`$) but we assume that $`D`$ is fixed and not all the solutions in (a) and (b) are shown in (c). For $`E>0`$ the solutions resemble those in Figure (2), with both monotonically expanding (solid) and collapsing (broken) solutions. The collapse singularity is naked ($`V_{min}<1`$) if $`E`$ is less than a value $`E_+(D)`$. For $`E<0`$ there are also solutions which recollapse to a black hole (solid) or emerge from a white hole (broken), as in the asymptotically Friedmann case. As $`E`$ decreases, the last solution is the symmetrical one for which $`E=E_{sym}`$ and $`z_S=1/D`$, so that the solid and broken curves coincide. The curves labelled $`E=0^+`$ and $`E=0^{}`$ show the qualitative transition as $`E`$ passes through $`0`$. The dotted curves correspond to an (unphysical) negative mass solution.

FIGURE (4). This shows the permitted regime for the parameters $`E`$ and $`D`$. The curve labelled $`E_{sym}`$ indicates the symmetric solution and all physical solutions must lie above this. The upper broken line gives the transition between different asymptotic forms for $`S_{\mathrm{}}`$ and $`V_{\mathrm{}}`$. The collapsing solutions have a naked singularity in the vertically shaded region below the line labelled $`E_+`$ and the black hole solutions have an event horizon and a particle horizon in the horizontally shaded region above the line labelled $`E_{}`$. These lines intersect on the lower boundary.


format: text

text:
# Size and area of square lattice polygons

## 1 Introduction

A self-avoiding polygon (SAP) can be defined as a walk on a lattice which returns to the origin and has no other self-intersections. The history and significance of this problem is nicely discussed in . Generally SAPs are considered distinct up to translations, so if there are $`p_n`$ SAPs of perimeter length $`n`$ there are $`2np_n`$ walks (the factor of two arising since the walk can go in two directions). In addition to enumerations by perimeter, one can also enumerate polygons by the enclosed area (or number of unit cells), or both perimeter and area. Of particular interest are the first few area-weighted moments of the perimeter generating function. Also of great interest is the mean-square radius of gyration, which measure the typical size of a SAP.

This paper builds on the work of Enting  who used transfer matrix techniques to enumerate square lattice polygons by perimeter to 38 steps. This enumeration was later extended by Enting and Guttmann to 46 steps  and then to 56 steps . This latter work also included calculations of moments of the caliper size distribution. Hiley and Sykes  obtained the number of square lattice polygons by both area and perimeter up to perimeter 18. Enting and Guttmann extended the calculation to perimeter 42 . The radius of gyration was calculated for SAPs up to 28 steps by Privman and Rudnick , using a technique based on direct counting of compact site animals on the dual lattice. Recently, Jensen and Guttmann devised an improved algorithm for the enumeration of SAPs and extended the calculation to 90 steps . The work reported here is based on generalisations of this improved algorithm. This has enabled us to extend the calculation of the radius of gyration and the first two area-weighted moments to 82 step SAPs. The generalisation of the transfer matrix technique to area-weighted moments is similar to the one used by Conway  in his calculation of series for percolation problems and lattice animals. The generalisation to the radius of gyration has to our knowledge no counterpart in the published literature, and represents a major advance in the design of efficient counting algorithms. Previous calculations of the radius of gyration were based on direct counting algorithms. The transfer matrix algorithm used in this paper is exponentially faster and thus enables us to significantly extend the series (see  for further details).

The size exponent, $`\nu `$, for self-avoiding polygons is believed to be identical to that of self-avoiding walks. This has been argued theoretically from the connection between the energy-energy and spin-spin correlation functions of the $`n`$-vector model in the limit $`n0`$, and SAPs and SAWs, respectively . Alternatively it has also been obtained from real space renormalization group arguments . The exponent describing the growth of the mean area of polygons of perimeter $`n`$ is expected to be $`2\nu `$ . Intuitively this is not surprising since it just means that the average area of a polygon is proportional to the square of the radius of gyration. So one is merely finding that for this problem the typical area and typical length scale match one another nicely. These expectations have been confirmed reasonably accurately by numerical work .

The functions we consider in this paper are: (i) the polygon generating function, $`𝒫(u)=p_nu^n`$; (ii) $`k^{th}`$ area-weighted moments of polygons of perimeter $`n`$, $`a^k_n`$; and (iii) the mean-square radius of gyration of polygons of perimeter $`n`$, $`R^2_n`$. These quantities are expected to behave as

$`p_n`$ $`=`$ $`B\mu ^nn^{\alpha 3}[1+o(1)],`$
$`a^k_n`$ $`=`$ $`E^{(k)}n^{2k\nu }[1+o(1)],`$ (1)
$`R^2_n`$ $`=`$ $`Dn^{2\nu }[1+o(1)],`$
where $`\mu =u_c^1`$ is the reciprocal of the critical point of the generating function, and $`\alpha =1/2`$ and $`\nu =3/4`$ are known exactly , though non-rigorously. It is also known  that the amplitude combination $`E^{(1)}/D`$ is universal, and that

$$BD=\frac{5}{32\pi ^2}\sigma a_0,$$
(2)
where $`a_0`$ is the area per site and $`\sigma `$ is an integer such that $`p_n`$ is non-zero only if $`n`$ is divisible by $`\sigma `$. For the square lattice $`a_0=1`$ and $`\sigma =2`$. These predictions have been confirmed numerically .

In the next section we describe the generalisation of the finite lattice method required in order to calculate the radius of gyration and area-weighted moments of self-avoiding polygons. The results of the analysis of the series are presented in Section 3.

## 2 Enumeration of self-avoiding polygons

The method used to enumerate self-avoiding polygons in this work is based on the method devised by Enting  for enumerations by perimeter and uses the enhancements of Jensen and Guttmann . In the following we first very briefly outline the original method and then show how to generalize it in order to calculate area-weighted moments and the radius of gyration. Details of the algorithm can be found in the papers cited above.

The first terms in the series for the perimeter generating function can be calculated using transfer matrix techniques to count the number of polygons spanning (in both directions) rectangles $`W+1`$ edges wide and $`L+1`$ edges long. The transfer matrix technique involves drawing a boundary through the rectangle intersecting a set of $`W+2`$ edges. For each configuration of occupied or empty edges along the boundary we maintain a (perimeter) generating function for partially completed polygons. Polygons in a given rectangle are enumerated by moving the boundary so as to add one site at a time. Due to the obvious symmetry of the lattice one need only consider rectangles with $`LW`$. Any polygon spanning such a rectangle has a perimeter of length at least $`2(W+L)`$. By adding the contributions from all rectangles of width $`WW_{\mathrm{max}}`$ (where the choice of $`W_{\mathrm{max}}`$ depends on available computational resources) and length $`WL2W_{\mathrm{max}}W+1`$, with contributions from rectangles with $`L>W`$ counted twice, the number of polygons per vertex of an infinite lattice is obtained correctly up to perimeter $`n_{\mathrm{max}}=4W_{\mathrm{max}}+2`$. The number of configurations required as $`W_{\mathrm{max}}`$ is increased grows exponentially as $`\lambda ^{W_{\mathrm{max}}}`$, where $`\lambda 2`$ for the improved algorithm . In addition to the dominant exponential growth in memory requirements there is a prefactor, which is proportional to the number of terms $`n_{\mathrm{max}}`$.

### 2.1 Area-weighted moments

Area-weighted moments can easily be calculated from the perimeter and area generating function

$$𝒞(u,v)=\underset{n,m}{}c_{n,m}u^nv^m,$$
(3)
where $`c_{n,m}`$ is the number of polygons with perimeter $`n`$ and area $`m`$. From this we get the area-weighted generating functions,

$$𝒫^{(k)}(u)=(v\frac{}{v})^k𝒞(u,v)|_{v=1}=\underset{n}{}\underset{m}{}m^kc_{n,m}u^n=\underset{n}{}p_n^{(k)}u^n,$$
(4)
and we define the average moments of area for a polygon with perimeter $`n`$

$$a^k_n=p_n^{(k)}/p_n^{(0)}=\underset{m}{}m^kc_{n,m}/p_n.$$
(5)
In order to calculate the moments of area through this approach we need to calculate a full two-parameter generating function, which generally will require a lot of computer memory. If we are only interested in the first few moments there is a much more efficient approach . We simply replace the variable $`v`$ by $`1+z`$ thus obtaining the function

$$F(u,z)=\underset{n,m}{}c_{n,m}u^n(1+z)^m=\underset{n,m}{}\underset{k=0}{\overset{m}{}}(\begin{array}{cc}m& \\ k& \end{array})c_{n,m}u^nz^m.$$
(6)
Let, $`F_i(u)`$, be the coefficient of $`z^i`$ in $`F(u,z)`$. Then we see that

$`F_0(u)`$ $`=`$ $`{\displaystyle \underset{n,m}{}}c_{n,m}u^n=𝒫(u),`$
$`F_1(u)`$ $`=`$ $`{\displaystyle \underset{n,m}{}}mc_{n,m}u^n=𝒫^{(1)}(u),`$
$`F_2(u)`$ $`=`$ $`{\displaystyle \underset{n,m}{}}m(m1)/2c_{n,m}u^n=[𝒫^{(2)}(u)𝒫^{(1)}(u)]/2,`$
and so on. Thus if we are only interested in the first and second moments of area we can truncate the series $`F(u,z)`$ at second order in $`z`$ and find the relevant moments as $`𝒫^{(1)}(u)=F_1(u)`$ and $`𝒫^{(2)}(u)=2F_2(u)+F_1(u)`$. The growth in memory requirements is still dominated by the exponential growth in the number of configurations. However, we have managed to turn the calculation of these moments from a problem with a prefactor cubic in $`W_{\mathrm{max}}`$ (the area is proportional to $`W_{\mathrm{max}}^2`$) into a problem with a prefactor linear in $`W_{\mathrm{max}}`$.

### 2.2 Radius of gyration

In the following we show how the definition of the radius of gyration can be expressed in a form suitable for a transfer matrix calculation. Note that we define the radius of gyration according to the vertices of the SAP and that the number of vertices equals the perimeter length. The radius of gyration of $`n`$ points at positions $`𝐫_i`$ is

$$n^2R_n^2=\underset{i>j}{}(𝐫_i𝐫_j)^2=(n1)\underset{i}{}(x_i^2+y_i^2)2\underset{i>j}{}(x_ix_j+y_iy_j).$$
(7)
This last expression is suitable for a transfer matrix calculation. As usual  we actually calculate the generating function, $`_g^2(u)=_np_nR^2_nn^2u^n`$. In order to do this we have to maintain five partial generating functions for each possible boundary configuration $`\sigma `$, namely

* $`P(u)`$, the number of (partially completed) polygons according to perimeter.
* $`R^2(u)`$, the sum over polygons of the squared components of the distance vectors.
* $`X(u)`$, the sum of the $`x`$-component of the distance vectors.
* $`Y(u)`$, the sum of the $`y`$-component of the distance vectors.
* $`XY(u)`$, the sum of the ‘cross’ product of the components of the distance vectors, e.g., $`_{i>j}(x_ix_j+y_iy_j)`$.

As the boundary line is moved to a new position each boundary configuration $`\sigma `$ might be generated from several configurations $`\sigma ^{}`$ in the previous boundary position. The partial generation functions are updated as follows

$`P(u,\sigma )`$ $`=`$ $`{\displaystyle \underset{\sigma ^{}}{}}u^{n(\sigma ^{})}P(u,\sigma ^{}),`$
$`R^2(u,\sigma )`$ $`=`$ $`{\displaystyle \underset{\sigma ^{}}{}}u^{n(\sigma ^{})}[R^2(u,\sigma ^{})+\delta (x^2+y^2)P(u,\sigma ^{})],`$
$`X(u,\sigma )`$ $`=`$ $`{\displaystyle \underset{\sigma ^{}}{}}u^{n(\sigma ^{})}[X(u,\sigma )+\delta xP(u,\sigma ^{})],`$ (8)
$`Y(u,\sigma )`$ $`=`$ $`{\displaystyle \underset{\sigma ^{}}{}}u^{n(\sigma ^{})}[Y(u,\sigma )+\delta yP(u,\sigma ^{})],`$
$`XY(u,\sigma )`$ $`=`$ $`{\displaystyle \underset{\sigma ^{}}{}}u^{n(\sigma ^{})}[XY(u,\sigma ^{})+\delta xX(u,\sigma ^{})+\delta yY(u,\sigma ^{})]`$
where $`n(\sigma ^{})`$ is the number of occupied edges added to the polygon and $`\delta =\mathrm{min}(n(\sigma ^{}),1)`$.

### 2.3 Further particulars

Finally a few remarks of a more technical nature. The number of contributing configurations becomes very sparse in the total set of possible states along the boundary line and as is standard in such cases one uses a hash-addressing scheme. Since the integer coefficients occurring in the series expansions become very large, the calculation was performed using modular arithmetic. Up to 8 primes were needed to represent the coefficients correctly. Further details and references are given in . The series for the radius of gyration and area-moments were calculated for SAPs with perimeter length up to 82. The maximum memory required for any given width did not exceed 2Gb. The calculations were performed on an 8 node Alpha Server 8400 with a total of 8Gb memory. The total CPU time required was about three days per prime. Obviously the calculation for each width and prime are totally independent and several calculations were done simultaneously.

In Table 1 we have listed the series for the radius of gyration and first and second area-weighted moments. The series for the radius of gyration of course agree with the terms up to length 28 computed previously , while the terms up to length 40 for the first area moment agree with the series in . The number of polygons of length $`56`$ can be found in  while those up to length 90 were reported in .

## 3 Analysis of the series

The series listed in Table 1 have coefficients which grow exponentially, with sub-dominant term given by a critical exponent. The generic behaviour is $`G(u)=_ng_nu^n(1u/u_c)^\xi ,`$ and hence the coefficients of the generating function $`g_n\mu ^nn^{\xi 1}`$, where $`\mu =1/u_c`$. To obtain the singularity structure of the generating functions we used the numerical method of differential approximants . In particular, we used this method to estimate the critical exponents (we already have very accurate estimates for $`u_c`$ from ). Since all odd terms in the series are zero and the first non-zero term is $`g_4`$ we actually analysed the function $`F(u)=_ng_{2n+4}u^n`$. Combining the relationship given above between the coefficients in a series and the critical behaviour of the corresponding generating function with the expected behaviour (1) of the mean-square radius of gyration and moments of area yields the following prediction for their generating functions:

$`_g^2(u)`$ $`=`$ $`{\displaystyle \underset{n}{}}p_{2n+4}R^2_{2n+4}n^2u^n={\displaystyle \underset{n}{}}r_nu^nR(u)(1u\mu ^2)^{(\alpha +2\nu )},`$ (9)
$`𝒫^{(k)}(u)`$ $`=`$ $`{\displaystyle \underset{n}{}}p_{2n+4}a^k_{2n+4}u^n={\displaystyle \underset{n}{}}a_n^{(k)}u^na^{(k)}(u)(1u\mu ^2)^{2(\alpha +2k\nu )}.`$ (10)
Thus we expect these series to have a critical point, $`u_c=1/\mu ^2=0.14368062927(1)`$, known to a very high degree of accuracy from the analysis in , and as stated previously the exponent $`\alpha =1/2`$, while it is expected that $`\nu =3/4`$.

Estimates of the critical point and critical exponents were obtained by averaging values obtained from second order $`[L/N;M;K]`$ inhomogeneous differential approximants. For each order $`L`$ of the inhomogeneous polynomial we averaged over those approximants to the series which used at least the first 80% - 90% of the terms of the series. We used only approximants where the difference between $`N`$, $`M`$, and $`K`$ didn’t exceed 2. Some approximants were excluded from the averages because the estimates were obviously spurious. The error quoted for these estimates reflects the spread (basically one standard deviation) among the approximants. Note that these error bounds should not be viewed as a measure of the true error as they cannot include possible systematic sources of error. In Table 2 we have listed the results of our analysis. It is evident that the estimates for $`u_c`$ and the critical exponents are in agreement with the expected behaviour. There are only some minor discrepancies in the fourth digit between the conjectured exponents and the estimates. This discrepancy is readily resolved by looking at the evidence in figure 1, where we have plotted the estimates for the critical point and exponent of $`_g^2`$. Each point in these figures represent an estimate obtained from a specific second order differential approximant with the various points obtained by varying the order of the polynomials in the approximants. It is clear that the estimates have not yet settled down to their asymptotic values and that they do converge towards the expected values as the number of terms used by the approximants is increased.

Now that the exact values of the exponents has been confirmed we turn our attention to the “fine structure” of the asymptotic form of the coefficients. In particular we are interested in obtaining accurate estimates for the amplitudes $`B`$, $`D`$ and $`E^{(1)}`$. We do this by fitting the coefficients to the assumed form (1).

The asymptotic form of the coefficients $`p_n`$ of the polygon generating function has been studied in detail previously . As argued in  there is no sign of non-analytic corrections-to-scaling exponents to the polygon generating function and one therefore finds that

$$p_n=\mu ^nn^{5/2}\underset{k=0}{}a_k/n^k.$$
(11)
This form was confirmed with great accuracy in . Estimates for the leading amplitude $`B=a_0`$ can thus be obtained by fitting $`p_n`$ to the form given in equation (11). In order to check the behaviour of such estimates we did the fitting using from 2 to 10 terms in the expansion. The results for the leading amplitude are displayed in figure 2. We notice that all fits appear to converge to the same value as $`n\mathrm{}`$, and that, as more and more correction terms are added to the fits the estimates exhibits less curvature and that the slope become smaller (although the fits using 10 terms are a little inconclusive). This is very strong evidence that (11) indeed is the correct asymptotic form of $`p_n`$. We estimate that $`B=0.5623012(1)`$.

The asymptotic form of the coefficients $`r_n`$ in the generating function for the radius of gyration has not been studied previously. When fitting to a form similar to equation (11), assuming that here are only analytic corrections-to-scaling, we find that the amplitudes of higher order terms are very large and that the leading amplitude converge rather slowly. This indicates that this asymptotic form is incorrect. We find that the coefficients fit better if we assume a leading non-analytic correction-to-scaling exponent $`\mathrm{\Delta }=3/2`$. This result confirms the prediction of Nienhuis . Note, that since the polygon generating function exponent $`2\alpha =3/2`$ a correction-to-scaling exponent $`\mathrm{\Delta }=3/2`$ is perfectly consistent with the asymptotic form (11). Because $`2\alpha +\mathrm{\Delta }`$ is an integer the non-analytic correction term becomes part of the analytic background term . We thus propose the following asymptotic form:

$$r_n=\mu ^nn[BD+\underset{k=0}{}a_k/n^{k/2}].$$
(12)
Alternative we could fit to the form

$$r_n/p_n=n^{7/2}[D+n^{5/2}\underset{k=0}{}a_k/n^{k/2}].$$
(13)
In figure 3 we show the leading amplitudes resulting from such fits while using from 1 to 10 terms in these expansions. Also shown in these figures (solid lines) are the predicted exact value of $`BD`$, given in equation 2, and the prediction for $`D`$ using the estimate for $`B`$ obtained above. As can be seen the leading amplitudes clearly converge towards their expected values and from these plots we can conclude that the prediction for $`BD`$ has been confirmed to at least 6 digit accuracy. Assuming that equation (2) is exact and using the very accurate estimate for $`B`$ we find that $`D=0.05630944(1)`$.

Fitting the coefficients for the area-weighted moments to asymptotic forms similar to equations (12) and (13) above (only the leading exponent was changed accordingly) leads to the estimates $`E^{(1)}=0.141520(1)`$ and $`E^{(2)}=0.0212505(4)`$.

As stated above the analysis of the polygon generating function is fully consistent with the prediction $`\mathrm{\Delta }=3/2`$. However, all one can conclude from the analysis is that, if non-analytic correction-to-scaling terms are present, the exponents have to be “half-integer”, so that the correction terms become part of the analytic background. The detailed analysis of the asymptotic form of the coefficients in the generating functions for the radius of gyration and area-weighted moments provide the firmest evidence to date for the actual existence of a leading non-analytic correction to scaling exponent $`\mathrm{\Delta }=3/2`$, thus confirming the theoretical predictions made by Nienhuis .

## 4 Conclusion

We have presented an improved algorithm for the calculation of the radius of gyration and area-weighted moments of self-avoiding polygons on the square lattice. This algorithm has enabled us to calculate these series for polygons up to perimeter length 82. Our extended series enables us to give very precise estimate of the critical exponents, which are consistent with the exact values $`\alpha =1/2`$ and $`\nu =3/4`$. We also obtain a very precise estimate for the amplitude $`B=0.5623012(1)`$. Analysis of the coefficients of the radius of gyration series yielded results fully compatible with the prediction $`BD=5/16\pi ^2`$. This allows us to obtain the very accurate estimate $`D=0.05630944(1)`$. From the first area-weighted moment we obtained the estimate $`E^{(1)}=0.141520(2)`$, which allows us to give a much improved estimate for the universal amplitude ratio $`E^{(1)}/D=2.51326(3)`$. We also find firm evidence for the existence of a non-analytic correction-to-scaling term with exponent $`\mathrm{\Delta }=3/2`$.

## E-mail or WWW retrieval of series

The series for the various generating functions studied in this paper can be obtained via e-mail by sending a request to I.Jensen@ms.unimelb.edu.au or via the world wide web on the URL http://www.ms.unimelb.edu.au/~iwan/ by following the instructions.

## 5 Acknowledgements

Thanks to Tony Guttmann for many valuable comments on the manuscript and the series analysis. Financial support from the Australian Research Council is gratefully acknowledged.