format: text

text:
# 1 Introduction

## 1 Introduction

Spatially homogeneous (SH) cosmologies, that is, cosmological solutions of the Einstein field equations that admit a local group of isometries acting on spacelike hypersurfaces, are of considerable importance in theoretical cosmology and have been much studied since the 1960s. These models can be used to analyze aspects of the physical Universe which pertain to or which may be affected by anisotropy in the rate of expansion, for example, the cosmic microwave background radiation, nucleosynthesis in the early Universe, and the question of the isotropization of the universe itself (see, for example, ).

Spatially homogeneous cosmologies also play an important role in attempts to understand the structure and properties of the space of all cosmological solutions of the Einstein field equations, since they are part of a symmetry-based hierarchy of cosmological models of increasing complexity, starting with the familiar Friedmann-Lemâitre models:

* Friedmann-Lemâitre cosmologies
* non-tilted SH cosmologies
* tilted SH cosmologies
* $`G_2`$ cosmologies
* $`G_1`$ cosmologies
* generic cosmologies

The terminology used in this hierarchy has the following meaning. A SH cosmology is said to be tilted if the fluid velocity vector is not orthogonal to the group orbits, otherwise the model is said to be non-tilted . A $`G_2`$ cosmology admits a local two-parameter Abelian group of isometries with spacelike orbits, permitting one degree of freedom as regards spatial inhomogeneity, while a $`G_1`$ cosmology admits one spacelike Killing vector field.

An important mathematical link between the various classes in the hierarchy is provided by the idea of representing the evolution using a state space. The physical state of a cosmological model at an instant of time is represented by a point in the state space, which is finite dimensional for classes i)-iii) and infinite dimensional otherwise. The Einstein field equations are formulated as first order evolution equations, and the evolution of a cosmological model is represented by an orbit (i.e., a solution curve) of the evolution equations in the state space. The state space of a particular class in the hierarchy is contained in the state spaces of the more general classes, which implies that the particular models are represented as special cases of the more general models. This structure opens the possibility that the evolution of a model in one class may be approximated, over some time interval, by a model in a more special class.

The models in each level of the hierarchy can be classified according to generality in various ways. For our purposes the most important of these is the algebraic classification of the isometry group of the SH models, the so-called Bianchi classification222This classification applies to SH models which admit an isometry group $`G_3`$ acting simply transitively on the hypersurfaces of homogeneity, and includes all SH models except for the Kantowski-Sachs models, which admit an isometry group $`G_4`$ acting on the hypersurfaces of homogeneity, but with no simply transitive subgroup $`G_3`$. (see, for example, , page 112). There is also a classification of the $`G_2`$ cosmologies, determined by the action of the isometry group, that is relevant for this paper. We refer to  (Table 12.4 on page 268) for details of how this classification relates to the Bianchi classification.

In this paper we assume that the matter content of the universe is a perfect fluid with equation of state $`p=(\gamma 1)\mu `$, where $`\gamma `$ is constant, the cases $`\gamma =1`$ (dust) and $`\gamma =\frac{4}{3}`$ (radiation) being of primary interest. Considerable work has been done in analyzing SH models, subject to this assumption, and a detailed, although still incomplete, description of the non-tilted models has been obtained. We refer to  (Chapters 6 & 7)and , for details. Much less is known about the tilted models, but it is evident that as one moves through stages i)-iii) in the hierarchy there is an increase in dynamical complexity and new features emerge at each stage. It is plausible to assume that as one moves to levels iv)-vi), which contain inhomogeneous models, this trend will continue. Clearly, understanding the dynamics at one level of complexity is a prerequisite for understanding the dynamics at a higher level. It is within this framework, which has a long term goal of obtaining qualitative information about the evolution of spatially inhomogeneous models, that the analysis of SH models assumes renewed importance.

For the class of tilted SH models, the Einstein field equations have been written as an autonomous DE in a number of different ways (see for example, ,  and ). Nevertheless, due to the complexity of the equations, a detailed analysis of the dynamics has not been given except in the case of a subclass of models of Bianchi type V -. Our goal in this paper is to give a qualitative analysis of the dynamics of the tilted SH cosmologies of Bianchi type II near the initial singularity and at late times, using the methods of the theory of dynamical systems. The Bianchi II cosmologies, while very special within the whole Bianchi class, nevertheless play a central role since the Bianchi II state space is part of the boundary of the state space for all higher Bianchi types (i.e. all types except for I and V). We thus expect that an analysis of the dynamics of the tilted Bianchi II class will give insight into the dynamics of the more general tilted Bianchi classes, while providing a lower bound for their dynamical complexity.

The paper is organized as follows. In Section 2 we present the evolution and constraint equations that arise from the Einstein field equations, and in Section 3 we describe the stability properties of the equilibrium points of the evolution equations, which provide the basis for the discussion of the dynamics in the singular asymptotic regime in Section 4 and in the late time asymptotic regime in Section 5. In Section 6 we discuss the implications of the results. Appendix A gives details of the derivation of the evolution equations.

Parts of the paper, in particular Sections 2 and 3, and the Appendices, are inevitably of a rather mathematical nature. Sections 4 and 5 are less technical, and in these sections we give the physical interpretations of the results. In the Introduction and in the Discussion we discuss the longer term goals of the research and its potential significance in a broader context.

The background material needed for this paper can be found in . In Appendix A it is assumed that the reader is familiar with the orthonormal frame formalism of Ellis & MacCallum  (see , Chapter 1). In addition, familiarity with some basic concepts and results from the theory of dynamical systems is assumed in Sections 2 and 3 (see , Chapter 4). We use geometrized units with $`c=1`$, $`8\pi G=1`$, and the sign conventions of .

## 2 Properties of the evolution equations

In order to write the Einstein field equations in a form amenable to dynamical systems analysis, we use the orthonormal frame formalism of Ellis & MacCallum . In this formalism the commutation functions of the orthonormal frame are used as the gravitational field variables, which has the advantage of leading directly to first order evolution equations for the gravitational field. The first step is to choose the orthonormal frame to be invariant under the group of isometries, which implies that the commutation functions depend only on a preferred time variable $`t`$. The second step is to choose the timelike frame vector $`𝐞_0`$ to be equal to the unit normal $`𝐧`$ of the group orbits, which is thus tangent to an irrotational congruence of geodesics. The third and key step is to make the commutation functions dimensionless by dividing them by the rate of expansion of the normal congruence, which leads to variables that remain bounded throughout the evolution of the models.

As described in Appendix A, the class of tilted Bianchi II cosmologies can be described by the following set of expansion-normalized variables:

$$𝐱=(\mathrm{\Sigma }_+,\mathrm{\Sigma }_{},\mathrm{\Sigma }_1,\mathrm{\Sigma }_3,N_1,v_3),$$
(2.1)
subject to one constraint of the form

$$g(𝐱)=h(v_3)\mathrm{\Omega }\mathrm{\Sigma }_3N_1=0,$$
(2.2)
where $`h(v_3)`$ is given by (A.19), i.e.

$$h(v_3)=\frac{\sqrt{3}\gamma v_3}{G},G=1+(\gamma 1)v_3^2.$$

Here $`\mathrm{\Sigma }_+,\mathrm{\Sigma }_{},\mathrm{\Sigma }_1,\mathrm{\Sigma }_3`$ are shear variables, $`N_1`$ is a spatial curvature variable and $`v_3`$ is a tilt variable. These variables determine the density parameter $`\mathrm{\Omega }`$ and the shear parameter $`\mathrm{\Sigma }`$ according to (A.17), (A.38) and (A.39). As mentioned at the end of Appendix A, we regard the off-diagonal shear variables $`\mathrm{\Sigma }_3`$ and $`\mathrm{\Sigma }_1`$ as representing the two tilt degrees of freedom.

The evolution equations for the variables (2.1), as derived in Appendix A, are given below:

$$\begin{array}{cc}\hfill \mathrm{\Sigma }_+^{}& =(2q)\mathrm{\Sigma }_+3\mathrm{\Sigma }_3^2+\frac{1}{3}N_1^2+\frac{1}{2\sqrt{3}}\mathrm{\Sigma }_3N_1v_3\hfill \\ \hfill \mathrm{\Sigma }_{}^{}& =(2q)\mathrm{\Sigma }_{}+2\sqrt{3}\mathrm{\Sigma }_1^2\sqrt{3}\mathrm{\Sigma }_3^2\frac{1}{2}\mathrm{\Sigma }_3N_1v_3\hfill \\ \hfill \mathrm{\Sigma }_1^{}& =(2q+2\sqrt{3}\mathrm{\Sigma }_{})\mathrm{\Sigma }_1\hfill \\ \hfill \mathrm{\Sigma }_3^{}& =(2q3\mathrm{\Sigma }_+\sqrt{3}\mathrm{\Sigma }_{})\mathrm{\Sigma }_3\hfill \\ \hfill N_1^{}& =(q4\mathrm{\Sigma }_+)N_1\hfill \\ \hfill v_3^{}& =\frac{v_3(1v_3^2)}{1(\gamma 1)v_3^2}(3\gamma 4\mathrm{\Sigma }_++\sqrt{3}\mathrm{\Sigma }_{}),\hfill \end{array}$$
(2.3)
where

$$q=2\left(1\frac{1}{12}N_1^2\right)\frac{1}{2}G^1\mathrm{\Omega }\left[3(2\gamma )(1v_3^2)+2\gamma v_3^2\right].$$

The auxiliary equation for $`\mathrm{\Omega }^{}`$ is

$$\mathrm{\Omega }^{}=G^1[2Gq(3\gamma 2)(2\gamma )v_3^2\gamma (\mathrm{\Sigma }_+\sqrt{3}\mathrm{\Sigma }_{})v_3^2]\mathrm{\Omega }.$$
(2.4)
The state space is the subset of $`^6`$ defined by the constraint (2.2) and the inequality $`\mathrm{\Omega }0`$, which by (A.17), (A.38) and (A.39), is equivalent to

$$\mathrm{\Sigma }_+^2+\mathrm{\Sigma }_{}^2+\mathrm{\Sigma }_1^2+\mathrm{\Sigma }_3^2+\frac{1}{12}N_1^21.$$
(2.5)
This restriction, and the fact that $`v_3^2<1`$, implies that the state space is bounded.

The evolution equations (2.3) are invariant under the transformations

$$(\mathrm{\Sigma }_+,\mathrm{\Sigma }_{},\mathrm{\Sigma }_1,\mathrm{\Sigma }_3,N_1,v_3)(\mathrm{\Sigma }_+,\mathrm{\Sigma }_{},\pm \mathrm{\Sigma }_1,\pm \mathrm{\Sigma }_3,\pm N_1,\pm v_3),$$

provided that the product $`v_3\mathrm{\Sigma }_3N_1`$ does not change sign. In addition, it follows that $`N_1,\mathrm{\Sigma }_1,\mathrm{\Sigma }_3`$ and $`v_3`$ cannot change sign along an orbit. Thus, without loss of generality, we can assume that

$$N_10,\mathrm{\Sigma }_10,\mathrm{\Sigma }_30,\text{and}v_30.$$

Taking these restrictions into account, the state space $`𝒟`$ of the tilted perfect fluid Bianchi II cosmologies is defined by the inequalities

$$N_1>0,\mathrm{\Omega }>0,0<v_3<1,\mathrm{\Sigma }_3>0,\mathrm{\Sigma }_10.$$
(2.6)
The boundary $`𝒟`$ is obtained by successively replacing the strict inequalities in (2.6) by equalities.

The evolution equations (2.3) are an autonomous DE in $`^6`$ of the form

$$𝐱^{}=𝐟(𝐱),$$

where the function $`𝐟:^6^6`$ on the right side is a rational function (note the function $`G`$ in the denominator in $`q`$, and the form of the $`v_3^{}`$ equation). For values of $`\gamma `$ that satisfy $`0<\gamma <2`$, the two functions in the denominator are strictly positive on the physical state space $`𝒟`$ and on its boundary $`𝒟`$. The DE (2.3) is thus smooth, indeed analytic, on the set $`𝒟𝒟`$. Since this set is compact and invariant, the solutions of the DE (2.3) can be extended for all $`\tau `$.

The constraint (2.2) entails a consistency requirement, namely that the equation $`g=0`$ should define an invariant set of the evolution equations (2.3). A straightforward calculation using (2.2)-(2.4) shows that

$$g^{}=(2q2\mathrm{\Sigma }_++\sqrt{3}\mathrm{\Sigma }_{})g.$$
(2.7)
It thus follows (see , Proposition 4.1 on page 29) that $`g=0`$ does indeed define an invariant set.

In analyzing a class of Bianchi cosmologies, one typically finds that the orbits in the boundary of the state space play a significant role in determining the dynamics, since orbits in the interior can shadow333i.e. approximate closely. See , page 104. orbits in the boundary. For the tilted Bianchi II models, the boundary of the state space is the union of five disjoint invariant sets, as shown in Table 2.1. The invariant sets i)-iv) describe familiar solutions but give multiple representations of them. For example, the orbits with $`\mathrm{\Sigma }_1=0`$ in the invariant set i) describe non-tilted Bianchi II cosmologies relative to a Fermi-propagated frame, while the orbits with $`\mathrm{\Sigma }_1>0`$ describe the same models, but relative to a rotating frame (see (A.37)). Likewise, the orbits with $`\mathrm{\Sigma }_1=\mathrm{\Sigma }_3=0`$ in the invariant set ii) describe (nontilted) Bianchi I cosmologies relative to a Fermi-propagated frame, while the orbits with $`\mathrm{\Sigma }_1>0`$ and/or $`\mathrm{\Sigma }_3>0`$ describe the same models, but relative to a rotating frame.

|  | Name | Restrictions | Dimension |
| --- | --- | --- | --- |
| i) | non-tilted non-vacuum Bianchi II | $`v_3=0=\mathrm{\Sigma }_3`$, $`N_1>0`$, $`\mathrm{\Omega }>0`$ | 4 |
| ii) | non-tilted non-vacuum Bianchi I | $`v_3=0=N_1`$, $`\mathrm{\Omega }>0`$ | 4 |
| iii) | Taub (vacuum Bianchi II) | $`\mathrm{\Sigma }_3=0=\mathrm{\Omega }`$, $`N_1>0`$, $`v_3<1`$ | 4 |
| iv) | Kasner (vacuum Bianchi I) | $`N_1=0=\mathrm{\Omega }`$, $`v_3<1`$ | 4 |
| v) | extreme tilt set $`(\gamma <2)`$ | $`v_3=1`$ | 4 |

Table I: The invariant sets that comprise the boundary of the state space (2.6).

The invariant set iii) contains the usual representation of the Taub vacuum solutions ( pages 137-8 and page 196) when $`\mathrm{\Sigma }_1=0=v_3`$, but it also contains multiple representations of these solutions relative to a non-Fermi-propagated frame. Likewise, the invariant set iv) contains the usual representation of the well-known Kasner vacuum solutions when $`\mathrm{\Sigma }_1=\mathrm{\Sigma }_3=v_3=0`$, but it also contains multiple representations of these solutions. The orbits in the remaining part of the boundary, the extreme tilt set v), do not correspond to spatially homogeneous cosmological solutions of the Einstein field equations since equation (A.21) breaks down when $`v_31`$ (i.e. $`v_bv^b1`$).

In addition to the above invariant sets, the DE (2.3) has one invariant set that is not a part of the boundary, namely the invariant set defined by $`\mathrm{\Sigma }_1=0`$. In view of the interpretation of $`\mathrm{\Sigma }_1`$, this invariant set corresponds to tilted Bianchi II models with one tilt degree of freedom. We shall see that the dynamics of these models is significantly simpler than the dynamics of the full class.

Technical point:

It is necessary for the subsequent analysis to determine at which points the constraint surface (2.2) is singular (i.e. the points $`𝐱^6`$ which satisfy $`g(𝐱)=0`$ and $`g(𝐱)=\mathrm{𝟎})`$, since at these points one cannot use the implicit function theorem to eliminate, locally, one of the variables. The gradient of $`g`$ is

$$g=(2h\mathrm{\Sigma }_+,2h\mathrm{\Sigma }_{},2h\mathrm{\Sigma }_1,2h\mathrm{\Sigma }_3N_1,\frac{1}{6}hN_1\mathrm{\Sigma }_3,h^{}\mathrm{\Omega }),$$
(2.8)
where the variables are listed in the order (2.1).

It follows that the surface is singular at and only at points $`𝐱`$ given by $`v_3=N_1=\mathrm{\Sigma }_3=\mathrm{\Omega }=0`$. Thus the surface is non-singular at all points of the state space (2.6), and is singular only on part of the boundary, i.e. a two-dimensional subset of the Kasner set. Indeed, if $`\mathrm{\Omega }>0`$, the constraint can be written in the form

$$h(v_3)=\frac{\mathrm{\Sigma }_3N_1}{\mathrm{\Omega }},$$

and since $`h^{}(v_3)>0`$ for $`0v_31`$ and $`0<\gamma <2`$, it follows that $`v_3`$ is determined uniquely in terms of the other variables.

It follows from (2.8) that the vectors

$$𝐞_A=\frac{}{x_A}\frac{1}{h^{}(v_3)\mathrm{\Omega }}\frac{g}{x_A}\frac{}{v_3},$$
(2.9)
where $`A=1,2,\mathrm{},5`$ and $`(x_A)=(\mathrm{\Sigma }_+,\mathrm{\Sigma }_{},\mathrm{\Sigma }_1,\mathrm{\Sigma }_3,N_1)`$, satisfy $`𝐞_Ag=0`$ and are linearly independent. These vectors thus span the tangent space to the constraint surface $`g=0`$ in $`^6`$. Equivalently, an orbit lies in the constraint surface $`g=0`$ if and only if its tangent vector is a linear combination of the $`𝐞_A`$.

## 3 Equilibrium points

In this Section we consider the local stability of the equilibrium points of the DE (2.3) and the constraint (2.2), i.e. points $`𝐱=𝐚`$ that satisfy

$$𝐟(𝐚)=\mathrm{𝟎},g(𝐚)=0.$$
(3.1)
The local stability is determined by linearizing the DE (2.3) at $`𝐱=𝐚`$, which gives

$$𝐱^{}=D𝐟(𝐚)𝐱,$$

and finding the eigenvalues of the derivative matrix $`D𝐟(𝐚)`$. The analysis is complicated by the constraint, which requires that we consider only eigenvectors that are tangent to the constraint surface, i.e. that are orthogonal to the gradient vector $`g(𝐚)`$. We shall refer to eigenvalues and eigenvectors that satisfy this condition as physical. We note that if all the physical eigenvalues have negative (positive) real parts then the equilibrium point is a local sink (source), i.e. it attracts (repels) all orbits in a neighbourhood. In addition to isolated equilibrium points, we will also encounter arcs of equilibrium points, for which one eigenvalue is necessarily zero (see for example , Section 4.3.4). In this case the criterion for a local sink (source) is that all eigenvalues other than the zero one have negative (positive) real parts.

We now list the equilibrium points, obtained by systematically solving equations (3.1). Each equilibrium point, apart from those with extreme tilt, corresponds to a self-similar solution of the Einstein field equations444The situation is analogous to the non-tilted case (, Section 5.2.3)., which we also give.

Non-vacuum equilibrium points $`(\mathrm{\Omega }>0)`$:

* Flat FL point, $`F`$

$$\mathrm{\Sigma }_+=\mathrm{\Sigma }_{}=\mathrm{\Sigma }_1=\mathrm{\Sigma }_3=N_1=v_3=0,$$

$$\mathrm{\Omega }=1,\mathrm{\Sigma }=0,q=\frac{1}{2}(3\gamma 2),0<\gamma 2.$$

Self-similar solution: the flat FL solution.
* Non-tilted point, PII

$$\mathrm{\Sigma }_+=\frac{1}{8}(3\gamma 2),N_1=\frac{3}{4}\sqrt{(3\gamma 2)(2\gamma )},\mathrm{\Sigma }_{}=\mathrm{\Sigma }_1=\mathrm{\Sigma }_3=v_3=0,$$

$$\mathrm{\Omega }=\frac{3}{16}(6\gamma ),\mathrm{\Sigma }=\frac{1}{8}(3\gamma 2),q=\frac{1}{2}(3\gamma 2),\frac{2}{3}<\gamma <2.$$

Self-similar solution: the Collins-Stewart solutions(, pages 131 and 189).
* Tilted point, PIItilt

$$\mathrm{\Sigma }_+=\frac{1}{8}(3\gamma 2),\mathrm{\Sigma }_{}=\frac{\sqrt{3}}{8}(107\gamma ),\mathrm{\Sigma }_3=\frac{1}{4}\alpha \sqrt{(11\gamma 10)(7\gamma 10)},\mathrm{\Sigma }_1=0,$$

$$N_1=\alpha \sqrt{3(5\gamma 4)(3\gamma 4)},v_3=\sqrt{\frac{(3\gamma 4)(7\gamma 10)}{(11\gamma 10)(5\gamma 4)}},$$

$$\alpha ^2=\frac{3(2\gamma )}{17\gamma 18},\frac{10}{7}<\gamma <2.$$

$$\mathrm{\Omega }=\frac{1}{4}\alpha ^2(21\gamma ^224\gamma +4),\mathrm{\Sigma }^2=1\alpha ^2(\gamma 1)(9\gamma 5),q=\frac{1}{2}(3\gamma 2).$$

Self-similar solution: first given in .
* Line of tilted points, $``$IItilt

$$\mathrm{\Sigma }_+=\frac{1}{3},\mathrm{\Sigma }_{}=\frac{1}{3\sqrt{3}},\mathrm{\Sigma }_3=\frac{2}{3}\sqrt{\frac{1}{57}(4b+1)(83b)},\mathrm{\Sigma }_1=\frac{2}{3\sqrt{3}}b,$$

$$N_1=2\sqrt{\frac{1}{57}(2b+1)(178b)},v_3=\sqrt{\frac{3(4b+1)(2b+1)}{(178b)(83b)}}$$

with

$$0<b<1,\gamma =\frac{14}{9}.$$

$$\mathrm{\Omega }=\frac{2}{171}(16b^245b+59),\mathrm{\Sigma }^2=\frac{4}{171}(2b+1)(92b),q=\frac{4}{3}.$$

Self-similar solutions: not given previously.
* Extreme tilted point, PIIextreme

$$\mathrm{\Sigma }_+=\frac{1}{3},\mathrm{\Sigma }_{}=\frac{1}{3\sqrt{3}},\mathrm{\Sigma }_3=\frac{10}{3\sqrt{57}},\mathrm{\Sigma }_1=\frac{2}{3\sqrt{3}},N_1=\frac{6}{\sqrt{19}},v_3=1$$

$$\mathrm{\Omega }=\frac{20}{57},\mathrm{\Sigma }^2=\frac{28}{57},q=\frac{4}{3},0<\gamma <2.$$
* Jacobs disc, $`𝒥`$

$$\mathrm{\Sigma }_1=\mathrm{\Sigma }_3=N_1=v_3=0,\mathrm{\Sigma }_+^2+\mathrm{\Sigma }_{}^2<1$$

$$\mathrm{\Omega }>0,\mathrm{\Sigma }<1,q=2,\gamma =2.$$

Self-similar solutions: the Jacobs stiff fluid solutions (, page 1109).

Vacuum equilibrium points $`(\mathrm{\Omega }=0)`$:

* Kasner circle, $`𝒦`$

$$\mathrm{\Sigma }_+^2+\mathrm{\Sigma }_{}^2=1,\mathrm{\Sigma }_3=\mathrm{\Sigma }_1=N_1=v_3=0$$

$$\mathrm{\Omega }=0,\mathrm{\Sigma }=1,q=2,0<\gamma 2.$$

Self-similar solutions: the Kasner vacuum solutions (, pages 132 and 188)
* Kasner circle with extreme tilt, $`𝒦_{\mathrm{𝑒𝑥𝑡𝑟𝑒𝑚𝑒}}`$

$$\mathrm{\Sigma }_+^2+\mathrm{\Sigma }_{}^2=1,\mathrm{\Sigma }_3=\mathrm{\Sigma }_1=N_1=0,v_3=1$$

$$\mathrm{\Omega }=0,\mathrm{\Sigma }=1,q=2,0<\gamma <2.$$
* Kasner lines with tilt, $`𝒦_{\mathrm{𝑡𝑖𝑙𝑡}}^\pm `$

$$\mathrm{\Sigma }_+=\sqrt{3}\mathrm{\Sigma }_{}+3\gamma 4,\mathrm{\Sigma }_{}=\frac{\sqrt{3}}{4}[3\gamma +4\pm \sqrt{(3\gamma 2)(2\gamma )}],\mathrm{\Sigma }_3=\mathrm{\Sigma }_1=0$$

$$N_1=0,0<v_3<1$$

$$\mathrm{\Omega }=0,\mathrm{\Sigma }=1,q=2,\frac{2}{3}\gamma 2.$$

Self-similar solutions: the Kasner vacuum solutions referred to a non-Fermi-propagated frame.

If $`\gamma `$ satisfies $`\frac{2}{3}<\gamma <2`$, there are two Kasner lines $`𝒦_{\mathrm{tilt}}^\pm `$, which join $`𝒦`$ to $`𝒦_{\mathrm{extreme}}`$, while if $`\gamma =\frac{2}{3}`$ or $`\gamma =2`$, the two lines coincide, with $`(\mathrm{\Sigma }_+,\mathrm{\Sigma }_{})=(\frac{1}{2},\frac{\sqrt{3}}{2})`$ and $`(\mathrm{\Sigma }_+,\mathrm{\Sigma }_{})=(\frac{1}{2},\frac{\sqrt{3}}{2})`$, respectively.

We need to know whether any of the equilibrium points are local sinks or sources. It turns out that for each value of $`\gamma `$ in the interval $`0<\gamma <2`$ one of the equilibria, with $`\mathrm{\Omega }>0`$ is a local sink555With the exception of $`F`$, many of the eigenvalues are complicated expressions in $`\gamma `$; their explicit form is unimportant for our purposes., as indicated in Table 3.1. The equilibria are related to one another by a series of bifurcations that occur as $`\gamma `$ varies. By inspection of the coordinates of these equilibrium points, we observe the following transitions:

$$F\stackrel{\gamma =\frac{2}{3}}{}\mathrm{PII}\stackrel{\gamma =\frac{10}{7}}{}\mathrm{PII}_{\mathrm{tilt}}\underset{b=0}{\overset{\gamma =\frac{14}{9}}{}}\mathrm{II}_{\mathrm{tilt}}\underset{b=1}{\overset{\gamma =\frac{14}{9}}{}}\mathrm{PII}_{\mathrm{extreme}}$$
(3.2)
| Range of $`\gamma `$ | Local sink |
| --- | --- |
| $`0<\gamma \frac{2}{3}`$ | flat FL point $`F`$ |
| $`\frac{2}{3}<\gamma \frac{10}{7}`$ | non-tilted point PII |
| $`\frac{10}{7}<\gamma <\frac{14}{9}`$ | tilted point PIItilt |
| $`\gamma =\frac{14}{9}`$ | line of tilted points $``$IItilt |
| $`\frac{14}{9}<\gamma <2`$ | extreme tilted point PIIextreme |

Table II: Local sinks in the tilted Bianchi II state space.

We can describe the mechanisms for these bifurcations, without giving full details, as follows. The linearization of the evolution equation for $`N_1`$ at $`F`$ is

$$N_1^{}=\frac{1}{2}(3\gamma 2)N_1,$$

showing that the spatial curvature variable $`N_1`$ destabilizes $`F`$ at $`\gamma =\frac{2}{3}`$. The associated eigenvector is $`𝐞_5`$ in (2.9). The linearizations of the evolution equations for $`\mathrm{\Sigma }_3`$ and $`v_3`$ at PII are

$$\mathrm{\Sigma }_3^{}=\frac{3}{8}(7\gamma 10)\mathrm{\Sigma }_3,v_3^{}=\frac{3}{8}(7\gamma 10)v_3,$$

showing that $`v_3`$ and $`\mathrm{\Sigma }_3`$ destabilize PII at $`\gamma =\frac{10}{7}`$. The associated physical eigenspace is actually one-dimensional and the eigenvector is $`𝐞_4`$. Finally, the linearization of $`\mathrm{\Sigma }_1^{}`$ at PIItilt is

$$\mathrm{\Sigma }_1^{}=\frac{3}{4}(9\gamma 14)\mathrm{\Sigma }_1,$$

showing that $`\mathrm{\Sigma }_1`$ destabilizes PIItilt at $`\gamma =\frac{14}{9}`$. Stability is transferred from PIItilt to PIextreme through the line of equilibrium points $``$IItilt, which exists only for $`\gamma =\frac{14}{9}`$. We shall refer to the bifurcation at $`\gamma =\frac{2}{3}`$ as the spatial curvature bifurcation, and, in view of the fact that $`\mathrm{\Sigma }_3`$ and $`\mathrm{\Sigma }_1`$ represent the tilt degrees freedom, we shall refer to the bifurcations at $`\gamma =\frac{10}{7}`$ and $`\gamma =\frac{14}{9}`$ as the first and second tilt bifurcations.

As regards local sources, it turns out that unless $`\gamma =2`$, none of the equilibrium points or equilibrium sets is a local source. This result follows from a careful analysis of the eigenvalues associated with the equilibrium points and sets. If $`\gamma =2`$, it turns out that a subset of the Jacobs disc shown as the shaded region in Figure I, is a local source.

Figure I: The shaded region, a subset of the Jacobs disc, is a local source in the case $`\gamma =2`$.

It is of interest to consider the invariant subset defined by $`\mathrm{\Sigma }_1=0`$, which describes the evolution of models with one tilt degree of freedom. For these models an analysis of the eigenvalues shows that there is an arc of the Kasner circle $`𝒦`$, defined by $`|\mathrm{\Sigma }_+|<\frac{1}{2}`$, $`\mathrm{\Sigma }_{}>0`$, that is a local source.

## 4 The late time asymptotic regime

We have seen that for each value of $`\gamma `$ in the range $`0<\gamma <2`$, excluding $`\gamma =\frac{14}{9}`$, there is a unique equilibrium point that is a local sink of the evolution equations, while if $`\gamma =\frac{14}{9}`$, there is an arc of equilibrium points that is a local sink. These local sinks are listed in Table 3.1, and the bifurcations that occur as $`\gamma `$ increases are displayed in equation (3.2). By definition of local sink, any orbit that enters a sufficiently small neighbourhood of the sink approaches the sink as $`\tau +\mathrm{}`$. The monotone functions in Appendix B, and numerical simulations, provide strong evidence that for a given value of $`\gamma `$ the local sink is the future attractor of the evolution equation, i.e. all orbits, except possibly a set of measure zero, approach the local sink as $`\tau +\mathrm{}`$.

The main conclusion that can be drawn from this asymptotic result is that the dynamical significance of the tilt and the shear at late times increases as the equation of state parameter increases from 0 to 2, as follows:

* in the range $`0<\gamma <\frac{2}{3}`$ the models isotropize, and since the deceleration parameter is asymptotically negative (i.e. $`\underset{\tau +\mathrm{}}{lim}q=\frac{1}{2}(3\gamma 2)<0`$), the models are inflationary,
* at $`\gamma =\frac{2}{3}`$ the spatial curvature destabilizes the flat FL equilibrium point, and for $`\gamma >\frac{2}{3}`$ the models no longer isotropize,
* at $`\gamma =\frac{10}{7}`$, the tilt destabilizes the Collins-Stewart solution, and for $`\gamma >\frac{10}{7}`$, the models are asymptotically tilted at late times, and
* if $`\gamma >\frac{14}{9}`$, the tilt is asymptotically extreme $`(v1)`$ at late times.

The $`\gamma `$-dependent limits of the dimensionless shear scalar $`\mathrm{\Sigma }`$ and the tilt variable $`v`$, as defined by (A.23), can be obtained from the list of equilibrium points in Section 3. In the physically important cases of dust $`(\gamma =1)`$ and radiation $`\left(\gamma =\frac{4}{3}\right)`$, which satisfy $`\frac{2}{3}<\gamma <\frac{10}{7}`$, the models are asymptotic to the Collins-Stewart solution i.e. they do not isotropize $`(\mathrm{\Sigma }0)`$, but the tilt becomes dynamically negligible $`(v0)`$, at late times.

## 5 The singular asymptotic regime

As mentioned in Section 3, there is no equilibrium point or equilibrium set that is a local source, except in the special case $`\gamma =2`$, or unless one restricts consideration to models with only one tilt degree of freedom (the invariant set $`\mathrm{\Sigma }_1=0`$). The implication of this fact is that a typical orbit is not past asymptotic to an equilibrium point. The situation is analogous to the case of non-tilted SH cosmologies of Bianchi types VIII and IX (see , Section 6.4), for which it has been shown that there exist infinite heteroclinic sequences based on the circle $`𝒦`$ of Kasner equilibrium points, i.e. infinite sequences of equilibrium points on $`𝒦`$, joined by special heteroclinic orbits directed into the past (see  and , Section 6.4.2, for a geometrical description of these heteroclinic sequences). These heteroclinic sequences determine the dynamics in the singular regime $`(\tau \mathrm{})`$ in the sense that a typical orbit shadows666In the language of dynamical systems, we are describing the $`\alpha `$-limit set of a typical orbit. For recent progress in proving the existence of the $`\alpha `$-limit set,we refer to  & \[17. (i.e. is approximated by) a heteroclinic sequence as $`\tau \mathrm{}`$. In physical terms, the dynamics of a typical cosmological model is approximated by a sequence of Kasner vacuum models as the singularity is approached into the past, the so-called Mixmaster oscillatory behaviour.

The present situation is more complicated due to the existence of two Kasner circles, the standard Kasner circle $`𝒦`$, and the Kasner circle $`𝒦_{\mathrm{extreme}}`$, with extreme tilt (i.e. $`v_3=1`$). The heteroclinic sequences contain orbits that join two points on $`𝒦`$, orbits that join two points on $`𝒦_{\mathrm{extreme}}`$, and orbits that join a point on $`𝒦`$ to a point on $`𝒦_{\mathrm{extreme}}`$, and vice versa. There are three families of orbits that join two points on the same Kasner circle. These orbits satisfy $`\mathrm{\Omega }=0`$, i.e.

$$\mathrm{\Sigma }_+^2+\mathrm{\Sigma }_{}^2+\mathrm{\Sigma }_1^2+\mathrm{\Sigma }_3^2+\frac{1}{12}N_1^2=1,$$

have one of $`\mathrm{\Sigma }_1,\mathrm{\Sigma }_3`$ and $`N_1`$ non-zero, and have

$$v_3=0\text{or}v_3=1.$$

The three families are given by

* $`N_1>0`$, $`\mathrm{\Sigma }_{}=C_1(\mathrm{\Sigma }_+2)`$,
* $`\mathrm{\Sigma }_3>0`$, $`\mathrm{\Sigma }_+\sqrt{3}\mathrm{\Sigma }_{}=C_2`$,
* $`\mathrm{\Sigma }_1>0`$, $`\mathrm{\Sigma }_+=C_3`$,

where $`C_1,C_2`$ and $`C_3`$ are constants. The projections of these orbits in the $`\mathrm{\Sigma }_+\mathrm{\Sigma }_{}`$-plane are shown in Figure II, with the arrows showing evolution into the past. We note that the orbits i) describe the Taub vacuum solutions of Bianchi type II (see , page 137), while the orbits ii) and iii) describe the Kasner vacuum solutions relative to a rotating frame (i.e. the angular velocity $`R_\alpha `$ of the spatial frame is non-zero; see Appendix A).

Figure II: Vacuum orbits joining points on the Kasner circle $`𝒦`$ or the extreme Kasner circle $`𝒦_{\mathrm{extreme}}`$, given by $`\mathrm{\Sigma }_+^2+\mathrm{\Sigma }_{}^2=1`$. The arrows show evolution into the past.

The orbits joining an equilibrium point on $`𝒦`$ to an equilibrium point on $`𝒦_{\mathrm{extreme}}`$ are given by

* $`\mathrm{\Sigma }_+,\mathrm{\Sigma }_{}=\text{constant}`$, $`\mathrm{\Sigma }_+^2+\mathrm{\Sigma }_{}^2=1`$, $`\mathrm{\Sigma }_1=\mathrm{\Sigma }_3=N_1=0`$,

with $`v_3^{}<0`$ (i.e. $`v_3`$ varying between 0 and 1 for evolution into the past), or $`v_3^{}>0`$. The direction of flow along these orbits is determined by the Kasner lines with tilt $`𝒦^\pm `$, which depend on the value of $`\gamma `$, and is shown in Figure III.

Figure III:777These invariant sets of orbits joining $`𝒦`$ and $`𝒦_{\mathrm{extreme}}`$ first appeared in the analysis of the tilted SH models of Bianchi type V (see  and ). Vacuum orbits joining points on the Kasner circle $`𝒦`$ and the extreme Kasner circle $`𝒦_{\mathrm{extreme}}`$. The arrows show evolution into the past.

If $`0<\gamma \frac{2}{3}`$, the flow into the past along these orbits is from $`𝒦`$ to $`𝒦_{\mathrm{extreme}}`$, while if $`\gamma =2`$, the flow into the past is from $`𝒦_{\mathrm{extreme}}`$ to $`𝒦`$. On the other hand, if $`\frac{2}{3}<\gamma <2`$, one family of orbits links $`𝒦`$ to $`𝒦_{\mathrm{extreme}}`$ while the second family does the reverse.

We can now describe the infinite heteroclinic sequences. In the case $`\frac{2}{3}<\gamma <2`$, the orbits that join successive Kasner points belong to the eight families i)-iii) with $`v_3=0`$ or $`v_3=1`$ and iv) with $`v_3^{}>0`$ or $`v_3^{}<0`$. In the case $`0<\gamma \frac{2}{3}`$, since $`\underset{\tau \mathrm{}}{lim}v_3=1`$, only orbits in the three families i)-iii) with $`v_3=1`$ are permitted.

It is not our intention to describe the detailed structure of the heteroclinic sequences, as has been done in the case of the non-tilted models of Bianchi types VIII and IX (see the references at the beginning of this section). We simply wish to point out that the evolution of tilted cosmologies of Bianchi type II in the singular asymptotic regime is governed by infinite heteroclinic sequences based on the Kasner vacuum models.

The case $`\gamma =2`$ is exceptional in that there is a local source, namely a subset of the Jacobs disc (see Figure I). We thus conjecture that a typical orbit with $`\gamma =2`$ is past asymptotic to an equilibrium point in this set. The models described by the invariant set $`\mathrm{\Sigma }_1=0`$ are also exceptional, since there is a subset of the Kasner circle $`𝒦`$ that is a local source (see the end of Section 3). We thus conjecture that a typical orbit in this set is past asymptotic to an equilibrium point in $`𝒦`$.

We conclude this section by comparing our description of the oscillatory behaviour in the singular regime using dynamical systems methods in the expansion-normalized state-space with the descriptions provided by the Hamiltonian approach (see  and , pages 63-4) and by the so-called BKL approach (, pages 533-8).

For the class of non-tilted SH models of Bianchi types VIII and IX this oscillatory behaviour is described as a succession of

* bounces off Bianchi type II potential walls, in the Hamiltonian approach,
* changes of the Kasner exponents as described by the BKL map (see , page 236), in the BKL approach, and
* vacuum Bianchi II orbits (Taub orbits) linking Kasner equilibrium points in the expansion-normalized state-space.

The presence of tilt leads to a new dynamical phenomenon which is due to the occurrence of off-diagonal shear degrees of freedom, and which is described in the above three approaches as follows:

* bounces off “centrifugal” potential walls in the Hamiltonian approach (see , pages 132-4),
* rotation of the Kasner axes, thereby permuting the Kasner exponents, in the BKL approach (see , pages 640-7), and
* non-singular Kasner orbits that link physically equivalent Kasner equilibrium points (see Figure II, ii) and iii)) in the expansion-normalized state-space.

In the expansion-normalized state space, the presence of tilt leads to an additional new phenomenon, namely transitions between states with zero tilt and states with extreme tilt, described by the non-singular Kasner orbits that link equilibrium points on $`𝒦`$ and $`𝒦_{\mathrm{extreme}}`$ (see Figure III). These orbits, that describe changes in the dynamical significance of the tilt during the oscillatory regime, do not appear to have an analogue in the other two approaches.

## 6 Discussion

In this paper we have introduced expansion-normalized variables to describe the evolution of tilted SH models of Bianchi type II, and have shown that with this choice of variable, the Einstein field equations reduce to an autonomous DE on a five dimensional compact subset of $`^6`$. Since the state space is compact, the orbits of the DE have a well-defined asymptotic behaviour as the dimensionless time $`\tau `$ tends to $`\pm \mathrm{}`$. This property enabled us to give a detailed description of the dynamics of the above class of models in the late time regime (see Section 4) and on approach to the initial singularity (see Section 5).

The description is incomplete in one respect, namely, that we have not discussed the kinematical properties of the cosmological fluid. For a non-tilted SH cosmology the only non-zero kinematical quantities of the fluid congruence are the rate of expansion and the rate of shear, since the acceleration and vorticity are necessarily zero . On the other hand, in a tilted SH cosmology, all four kinematical quantities of the fluid congruence are non-zero in general . It is essential to keep in mind, however, that in our analysis, the dimensionless shear variables (2.1) and the Hubble scalar $`H`$ describe the kinematics of the timelike congruence that is normal to the group orbits. The kinematical quantities of the fluid congruence can be expressed in terms of the variables (2.1) by adapting the formulas in  (see sections 1 and 2). The desired relations, valid for all Bianchi types, are given in Appendix C, and can easily be specialized to models of Bianchi type II.

We can draw the following conclusions from these relations:

* The cosmological fluid is expanding for all time. This fact follows from (C.1) and (C.2), and the fact that the Hubble scalar of the normal congruence is assumed to be positive.
* The vorticity of the fluid is zero888This result can be inferred directly from  (see Theorem 3.1)., but the acceleration is non-zero, provided that the pressure is non-zero $`(\gamma 1)`$.
* If the dimensionless shear of the normal congruence is small, then the dimensionless shear and acceleration of the fluid congruence are small.

It is known that the evolution of non-tilted SH cosmologies of Bianchi type II has a very simple description in the expansion-normalized state space: the orbit of a model joins two equilibrium points, or is itself an equilibrium point.999See , Section 6.3.2. We refer to  (page 53) and , for alternative descriptions. Since equilibrium points correspond to self-similar solutions of the Einstein field equations, one says that the models are asymptotically self-similar at the initial singularity $`(\tau \mathrm{})`$ and at late times $`(\tau +\mathrm{})`$. The initial asymptotic state is typically a Kasner solution and the late time asymptotic state is the Collins-Stewart perfect fluid solution. We have seen that a Bianchi II model can in general accommodate two tilt degrees of freedom, described by the two off-diagonal shear variables $`\mathrm{\Sigma }_1`$ and $`\mathrm{\Sigma }_3`$. Our analysis shows that models with only one tilt degree of freedom (i.e. $`\mathrm{\Sigma }_1=0`$) are still asymptotically self-similar. The initial asymptotic state is again a Kasner state, implying that the tilt is not dynamically significant near the initial singularity. The late time asymptotic state101010This asymptotic behaviour has previously been discussed by . depends on the equation of state parameter $`\gamma `$, and is described by the sequence of bifurcations

$$F\stackrel{\gamma =\frac{2}{3}}{}PII\stackrel{\gamma =\frac{10}{7}}{}PII_{\mathrm{tilt}},$$

a subset of (3.2). If the second tilt degree of freedom is activated (i.e. $`\mathrm{\Sigma }_10`$), we have seen that evolution is Mixmaster-like near the initial singularity, and hence the models are no longer asymptotically self-similar in the singular regime. The late time asymptotic state is described by the sequence of bifurcations (3.2), showing that if $`\gamma `$ satisfies $`\frac{14}{9}<\gamma <2`$ the tilt becomes extreme. For models containing dust $`(\gamma =1)`$ or radiation $`\left(\gamma =\frac{4}{3}\right)`$, however, the models are asymptotically self-similar, with the Collins-Stewart solution as the late time asymptote.

The above results lead to some interesting comparisons concerning the ways in which various anisotropies, in particular, spatial curvature, magnetic fields and tilt, affect the dynamics of SH cosmologies.

It is known that the anisotropic spatial curvature in a SH cosmology, as described by the variables $`N_\alpha `$, $`\alpha =1,2,3`$, affects the stability of the circle of Kasner equilibrium points (see , page 132). If all of the $`N_\alpha `$ are non-zero, corresponding to models of Bianchi types VIII and IX, all Kasner points become saddles, leading to the existence of infinite heteroclinic sequences into the past, and hence an oscillatory singular regime. It has also been found that magnetic fields affect the Kasner equilibrium points in an analogous way, although the mechanism is somewhat different in that a magnetic field with more than one degree of freedom generates off-diagonal shear components. We refer to - for details. The present paper shows that in SH models of Bianchi type II, the tilt degrees of freedom affect the Kasner circle in the same way by generating off-diagonal shear components, leading to an oscillatory singular regime. Since Bianchi II orbits are contained in the boundary of the state space of models of more general Bianchi type (all except types I and V), we expect that this mechanism will operate quite generally for tilted Bianchi models.

There are also analogies in the way in which anisotropic spatial curvature, magnetic fields and tilt affect the dynamics of SH cosmologies at late times. If the equation of state parameter satisfies $`0<\gamma <\frac{2}{3}`$, then all orbits are asymptotic in the future to the flat FL equilibrium point $`F`$, implying that the models isotropize. Since in this case the limit of the deceleration parameter $`q`$ is negative $`(\underset{\tau \mathrm{}}{lim}q=\frac{1}{2}(3\gamma 2)<0)`$, the models are inflationary. For each of the three types of anisotropy, bifurcations occur as $`\gamma `$ increases, the result of which is that these anisotropies influence the dynamics to an increasing degree at late times. First, the spatial curvature destabilizes the flat FL equilibrium point $`F`$ at the value $`\gamma =\frac{2}{3}`$, by the creation of new equilibrium points with non-zero spatial curvature. At larger values of $`\gamma `$, a magnetic field in models of Bianchi types II and VI0 destabilizes these new equilibrium points, with the result that the magnetic field becomes dynamically significant (see -). In this paper we have shown that the tilt in a Bianchi II model acts in a similar way, becoming dynamically significant if $`\gamma >\frac{10}{7}`$.

There is, however, an important difference between the effect of tilt in a SH model of Bianchi type II and the effect of spatial curvature and magnetic fields. We have mentioned that the spatial curvature destabilizes the flat FL point $`F`$ at $`\gamma =\frac{2}{3}`$. It has also been shown that a magnetic field leads to a second destabilization of $`F`$ at $`\gamma =\frac{4}{3}`$ (see ). On the other hand, the analysis in Section 3 shows that in a Bianchi II model the tilt does not destabilize $`F`$. In other words, the flat FL equilibrium point $`F`$ is stable with respect to the off-diagonal shear degrees of freedom that arise in the presence of tilt. In the language of dynamical systems the tilt increases the dimension of the stable manifold of the equilibrium point $`F`$ making it more likely that orbits will pass close to $`F`$. The physical significance of this result is that in Bianchi II models tilt increases the probability of intermediate isotropization with $`\mathrm{\Omega }1`$. Whether this result is true in general for tilted SH models requires further investigation.

## Appendix A

In this Appendix we derive the evolution equations (2.3)-(2.4) for the tilted SH models of Bianchi type II. We begin by giving the evolution equations for a general SH model, in terms of expansion-normalized variables defined relative to the timelike congruence that is normal to the group orbits. These equations follow directly from the standard orthonormal frame equations for SH models as given in  (see equation (1.90)-(1.100)), as we now describe.

We introduce a group-invariant orthonormal frame $`\{𝐞_0,𝐞_\alpha \}`$, where $`𝐞_0=𝐧`$ is the unit normal to the group orbits. The non-zero commutation functions are (, page 39)

$$\{H,\sigma _{\alpha \beta },n_{\alpha \beta },a_\alpha ,\mathrm{\Omega }_\alpha \},$$

and the energy-momentum tensor

$$T_{ab}=\mu n_an_b+2q_{(a}n_{b)}+p(g_{ab}+n_an_b)+\pi _{ab}$$
(A.1)
is described by the source terms

$$\{\mu ,p,q_\alpha ,\pi _{\alpha \beta }\},$$

relative to the chosen orthonormal frame.

The expansion-normalized commutation functions and source terms are defined by

$`\mathrm{\Sigma }_{\alpha \beta }`$ $`={\displaystyle \frac{\sigma _{\alpha \beta }}{H}},`$ $`N_{\alpha \beta }`$ $`={\displaystyle \frac{n_{\alpha \beta }}{H}},`$ $`A_\alpha `$ $`={\displaystyle \frac{a_\alpha }{H}},`$ $`R_\alpha `$ $`={\displaystyle \frac{\mathrm{\Omega }_\alpha }{H}},`$ (A.2)
$`\mathrm{\Omega }`$ $`={\displaystyle \frac{\mu }{3H^2}},`$ $`P`$ $`={\displaystyle \frac{p}{3H^2}},`$ $`Q_\alpha `$ $`={\displaystyle \frac{q_\alpha }{H^2}},`$ $`\mathrm{\Pi }_{\alpha \beta }`$ $`={\displaystyle \frac{\pi _{\alpha \beta }}{H^2}}.`$ (A.3)
We also require the shear parameter $`\mathrm{\Sigma }`$, defined by

$$\mathrm{\Sigma }^2=\frac{\sigma ^2}{3H^2},$$
(A.4)
( equation (5.17)) and the Hubble-normalized spatial curvature variables, defined by

$$𝒮_{\alpha \beta }=\frac{{}_{}{}^{3}S_{\alpha \beta }^{}}{H^2},K=\frac{{}_{}{}^{3}R}{6H^2}.$$
(A.5)
It follows from (A.2) and equations (1.28), (1.94) and (1.95) in  that

$$\mathrm{\Sigma }^2=\frac{1}{6}\mathrm{\Sigma }_{\alpha \beta }\mathrm{\Sigma }^{\alpha \beta },$$
(A.6)
and

$$𝒮_{\alpha \beta }=B_{\alpha \beta }\frac{1}{3}B_\mu ^\mu \delta _{\alpha \beta }2\epsilon ^{\mu \nu }{}_{(\alpha }{}^{}N_{\beta )\mu }^{}A_\nu ,K=\frac{1}{12}B_\mu ^\mu +A_\mu A^\mu .$$
(A.7)
where

$$B_{\alpha \beta }=2N_\alpha ^\mu N_{\mu \beta }N_\mu ^\mu N_{\alpha \beta }.$$

We introduce the usual dimensionless time variable $`\tau `$ according to

$$\frac{dt}{d\tau }=\frac{1}{H},$$
(A.8)
where $`t`$ is clock time along the normal congruence (see , page 113). In making the transition to expansion-normalized variables, the deceleration parameter $`q`$ plays an essential role (see , pages 112-3). This dimensionless scalar determines the evolution of the Hubble scalar $`H`$ according to

$$H^{}=(1+q)H,$$
(A.9)
(, equation (5.11)). Raychaudhuri’s equation (, equation (1.90), in conjunction with (A.3), (A.4), (A.8) and (A.9), leads to the following algebraic expression for $`q`$:

$$q=2\mathrm{\Sigma }^2+\frac{1}{2}(\mathrm{\Omega }+3P).$$
(A.10)
The Einstein field equations (, equations (1.91)-(1.93)), Jacobi identities (, equations (1.96)-(1.98)) and contracted Bianchi identities (, equations (1.99) and (1.100)), in conjunction with (A.2), (A.3), (A.5), (A.8) and (A.9), now yield the following sets of equations.

Gravitational evolution equations:

$`\mathrm{\Sigma }_{\alpha \beta }^{}`$ $`=(2q)\mathrm{\Sigma }_{\alpha \beta }+2\epsilon _{(\alpha }^{\mu \nu }\mathrm{\Sigma }_{\beta )\mu }R_\nu 𝒮_{\alpha \beta }+\mathrm{\Pi }_{\alpha \beta },`$ (A.11)
$`N_{\alpha \beta }^{}`$ $`=qN_{\alpha \beta }+2\mathrm{\Sigma }_{(\alpha }^\mu N_{\beta )\mu }+2\epsilon _{(\alpha }^{\mu \nu }N_{\beta )\mu }R_\nu ,`$ (A.12)
$$A_\alpha ^{}=qA_\alpha \mathrm{\Sigma }_\alpha ^\beta A_\beta +\epsilon _\alpha ^{\mu \nu }A_\mu R_\nu ,$$
(A.13)
Source evolution equations:

$$\mathrm{\Omega }^{}=(2q1)\mathrm{\Omega }3P\frac{1}{3}\mathrm{\Sigma }_\alpha ^\beta \mathrm{\Pi }_\beta ^\alpha +\frac{2}{3}A_\alpha Q^\alpha $$
(A.14)
$$Q_\alpha ^{}=2(q1)Q_\alpha \mathrm{\Sigma }_\alpha ^\beta Q_\beta \epsilon _\alpha ^{\mu \nu }R_\mu Q_\nu +3A^\beta \mathrm{\Pi }_{\alpha \beta }+\epsilon _\alpha ^{\mu \nu }N_\mu ^\beta \mathrm{\Pi }_{\beta \nu }.$$
(A.15)
Algebraic equations:

$$N_\alpha ^\beta A_\beta =0,$$
(A.16)
$$\mathrm{\Omega }=1\mathrm{\Sigma }^2K,$$
(A.17)
$$Q_\alpha =3\mathrm{\Sigma }_\alpha ^\mu A_\mu \epsilon _\alpha ^{\mu \nu }\mathrm{\Sigma }_\mu ^\beta N_{\beta \nu }.$$
(A.18)
It should be noted that the source evolution equations (i.e. the contracted Bianchi identities) are a consequence of the gravitational evolution equations and the algebraic equations, and hence contain no additional information. It is, however, convenient to use them as auxiliary equations.

It should also be noted that equations (A.11)-(A.13), with (A.10), (A.17) and (A.18) do not form a fully determined system of evolution equations. First, there is no evolution equation for the variables $`R_\alpha `$ that represent the angular velocity of the spatial frame $`\{𝐞_\alpha \}`$. One can in fact use the freedom in the choice of spatial frame, i.e. an arbitrary time-dependent rotation, to introduce a non-rotating spatial frame $`(R_\alpha =0)`$. This choice is not usually the most convenient one, however. In addition, there is neither an evolution equation nor an algebraic equation for the isotropic pressure $`P`$ and the anisotropic stress matrix $`\mathrm{\Pi }_{\alpha \beta }`$. Thus in order to obtain a fully determined system one has to specify the source and fix the spatial frame.

We now consider the case of a tilted perfect fluid, with stress energy tensor

$$T_{ab}=\stackrel{~}{\mu }u_au_b+\stackrel{~}{p}(g_{ab}+u_au_b),$$
(A.19)
and equation of state

$$\stackrel{~}{p}=(\gamma 1)\stackrel{~}{\mu }.$$
(A.20)
The 4-velocity $`𝐮`$ can be written in the form

$$u^a=\frac{1}{\sqrt{1v_bv^b}}(n^a+v^a),$$
(A.21)
where the spacelike vector $`𝐯`$ is orthogonal to the normal vector $`𝐧`$, and satisfies $`0v_bv^b<1`$. The vector $`𝐯`$ is called the tilt vector of the fluid. It has components $`(0,v^\alpha )`$ relative to the orthonormal frame $`\{𝐧,𝐞_\alpha \}`$.

We can now express the source terms $`P,Q_\alpha `$ and $`\mathrm{\Pi }_{\alpha \beta }`$ in terms of the density parameter $`\mathrm{\Omega }`$ and the tilt vector $`v_\alpha `$, as follows. First, equations (A.1) and (A.19)-(A.21) imply that

$$\mu =\frac{G}{1v^2}\stackrel{~}{\mu },$$
(A.22)
where

$$G=1+(\gamma 1)v^2,v^2=v_\alpha v^\alpha <1.$$
(A.23)
By using (A.22), in conjunction with (A.3), we now obtain

$`P`$ $`=\frac{1}{3}G^1\left[3(\gamma 1)(1v^2)+\gamma v^2\right]\mathrm{\Omega },`$ (A.24)
$`Q_\alpha `$ $`=3\gamma G^1\mathrm{\Omega }v_\alpha ,`$ (A.25)
$`\mathrm{\Pi }_{\alpha \beta }`$ $`=3\gamma G^1\mathrm{\Omega }\left(v_\alpha v_\beta \frac{1}{3}v^2\delta _{\alpha \beta }\right).`$ (A.26)
The algebraic equation (A.18), with (A.25), assumes the form

$$3\gamma G^1\mathrm{\Omega }v_\alpha =3\mathrm{\Sigma }_\alpha ^\mu A_\mu \epsilon _\alpha ^{\mu \nu }\mathrm{\Sigma }_\mu ^\beta N_{\beta \mu }.$$
(A.27)
Equations (A.17) and (A.27) express the source terms $`\mathrm{\Omega }`$ and $`v_\alpha `$ algebraically in terms of the gravitational field variable $`\mathrm{\Sigma }_{\alpha \beta },N_{\alpha \beta }`$ and $`A_\alpha `$.

It is, however, convenient to use the source evolution equations (A.14) and (A.15) to obtain evolution equations for $`\mathrm{\Omega }`$ and $`v_\alpha `$. On substituting (A.24)-(A.26) in (A.14) and (A.15) and rearranging we obtain

$$\mathrm{\Omega }^{}=G^1\left[2Gq(3\gamma 2)(2\gamma )v^2\gamma \mathrm{\Sigma }_{\alpha \beta }v^\alpha v^\beta +2\gamma A_\alpha v^\alpha \right]\mathrm{\Omega },$$
(A.28)
and111111For details of the derivation of this equation we refer to , equations (65) and (94).

$$\begin{array}{cc}\hfill v_\alpha ^{}& =\frac{v^\alpha }{[1(\gamma 1)v^2]}[(3\gamma 4)(1v^2)+(2\gamma )\mathrm{\Sigma }_{\beta \gamma }v^\beta v^\gamma \hfill \\ & +\{(2\gamma )(\gamma 1)(1v^2)\}A_\beta v^\beta ]\mathrm{\Sigma }_\alpha ^\beta v_\beta \hfill \\ & +\epsilon _\alpha ^{\mu \nu }(R_\mu +N_\mu ^\delta v_\delta )v_\nu v^2A_\alpha .\hfill \end{array}$$
(A.29)
Finally, we substitute (A.24) in (A.10) to express $`q`$ in terms of the gravitational field variables:

$$q=2\mathrm{\Sigma }^2+\frac{1}{2}G^1\mathrm{\Omega }\left[(3\gamma 2)(1v^2)+2\gamma v^2\right].$$
(A.30)
We note that (A.17) can be used to write $`q`$ in the form

$$q=2(1K)\frac{1}{2}G^1\mathrm{\Omega }[3(2\gamma )(1v^2)+2\gamma v^2]$$
(A.31)
For the class of Bianchi II cosmologies we have $`A_\alpha =0`$, and in addition we can choose the spatial frame $`\{𝐞_\alpha \}`$ to be an eigenframe of the matrix $`N_{\alpha \beta }`$, with

$$N_{11}0,N_{\alpha \beta }=0\text{otherwise}.$$
(A.32)
These restrictions, in conjunction with the constraint (A.27), imply

$$v_1=0,$$

assuming $`\mathrm{\Omega }>0`$ and $`\gamma >0`$. In addition, equation (A.12) gives $`R_2=\mathrm{\Sigma }_{13}`$, and $`R_3=\mathrm{\Sigma }_{12}`$. We are free to perform a rotation in the 23-plane to get121212This choice is the one made in , page 223.

$$v_2=0,v_30.$$

The constraint (A.27) now yields

$$\mathrm{\Sigma }_{13}=0=R_2,\mathrm{\Sigma }_{12}0,$$

$$3\gamma \mathrm{\Omega }v_3=G\mathrm{\Sigma }_{12}N_{11}.$$
(A.33)
Using these results the $`\mathrm{\Sigma }_{13}^{}`$ equation implies

$$R_1=\mathrm{\Sigma }_{23},$$

and at this stage $`R_\alpha `$ is uniquely determined in terms of $`\mathrm{\Sigma }_{\alpha \beta }`$.

We now relabel the variables as follows:

$$\mathrm{\Sigma }_+=\frac{1}{2}(\mathrm{\Sigma }_{22}+\mathrm{\Sigma }_{33}),\mathrm{\Sigma }_{}=\frac{1}{2\sqrt{3}}(\mathrm{\Sigma }_{22}\mathrm{\Sigma }_{33}),$$

$$\mathrm{\Sigma }_1=\frac{1}{\sqrt{3}}\mathrm{\Sigma }_{23},\mathrm{\Sigma }_3=\frac{1}{\sqrt{3}}\mathrm{\Sigma }_{12},$$

$$N_1=N_{11}.$$

The set of independent expansion-normalized variables is

$$(\mathrm{\Sigma }_+,\mathrm{\Sigma }_{},\mathrm{\Sigma }_1,\mathrm{\Sigma }_3,N_1,v_3),$$

subject to one constraint (A.33), which we now write in the form

$$h(v_3)\mathrm{\Omega }=\mathrm{\Sigma }_3N_1,$$
(A.34)
where

$$h(v_3)=\frac{\sqrt{3}\gamma v_3}{G},$$
(A.35)
and $`G`$ is given by

$$G=1+(\gamma 1)v_3^2,$$
(A.36)
as follows from (A.23) and the restrictions on $`v_\alpha `$.

It should be noted that the off-diagonal shear components $`\mathrm{\Sigma }_1`$ and $`\mathrm{\Sigma }_3`$ determine the angular velocity of the spatial frame according to

$$(R_\alpha )=\sqrt{3}(\mathrm{\Sigma }_1,0,\mathrm{\Sigma }_3).$$
(A.37)
The shear parameter has the simple form

$$\mathrm{\Sigma }^2=\mathrm{\Sigma }_+^2+\mathrm{\Sigma }_{}^2+\mathrm{\Sigma }_1^2+\mathrm{\Sigma }_3^2,$$
(A.38)
and the curvature parameter $`K`$ is given by

$$K=\frac{1}{12}N_1^2,$$
(A.39)
as follows from (A.7) and (A.32). The deceleration parameter $`q`$, as given by (A.31), now assumes the form

$$q=2\left(1\frac{1}{12}N_1^2\right)\frac{1}{2}G^1\mathrm{\Omega }\left[3(2\gamma )(1v_3^2)+2\gamma v_3^2\right].$$
(A.40)
The evolution equations (2.3)-(2.4) are now obtained by specializing equations (A.11), (A.12), (A.28) and (A.29) using the restrictions obtained above.

We conclude this Appendix with an important remark concerning the interpretation of the expansion-normalized variables. In Bianchi II models there are only two tilt degrees of freedom instead of the customary three (see , Table 9.5, page 211). The reason is that the matrix $`N_{\alpha \beta }`$ has two zero eigenvalues which, in conjunction with the constraint (A.27), implies that

$$N_\alpha ^\beta v_\beta =0,$$

i.e. the tilt vector lies in the two-dimensional null eigenspace of $`N_\alpha ^\beta `$. Our choice of frame, which leads to $`v_1=v_2=0`$, obscures the fact that there are two tilt degrees of freedom. The constraint (A.34) shows that the shear variable $`\mathrm{\Sigma }_3`$ represents one of the tilt degrees of freedom. It transpires that the shear variable $`\mathrm{\Sigma }_1`$ represents the second tilt degree of freedom, since if $`\mathrm{\Sigma }_1=0`$, the tilt vector is an eigenvector of $`\mathrm{\Sigma }_{\alpha \beta }`$, i.e.

$$v_{[\alpha }\mathrm{\Sigma }_{\beta ]}^\mu v_\mu =0,$$

which fixes its direction uniquely.

## Appendix B

In this Appendix we give some functions that are monotone along the orbits of the evolution equations, depending on the value of $`\gamma `$.

* Consider the function

$$Z=\frac{\alpha ^3\mathrm{\Sigma }_1^4\mathrm{\Sigma }_3^2}{\mathrm{\Omega }^3},$$
(B.1)
where

$$\alpha =\frac{Gv_3^2}{(1v_3^2)^{\frac{3}{2}(2\gamma )}},$$

and $`G`$ is given by equation (A.36). It follows131313The calculation is simpler if one uses the constraint (2.2) to express $`\mathrm{\Omega }`$ in terms of the other variables. from the evolution equations (2.3)-(2.4) that

$$Z^{}=3(9\gamma 14)Z.$$
(B.2)
It thus follows that if $`\gamma `$ satisfies $`0<\gamma <2`$, $`\gamma \frac{14}{9}`$, then $`Z`$ is a monotone function in the invariant set defined by $`\mathrm{\Omega }>0`$, $`v_3>0`$, $`\mathrm{\Sigma }_1>0`$ and $`\mathrm{\Sigma }_3>0`$. Furthermore, if $`0<\gamma <\frac{14}{9}`$, then

$$\underset{\tau +\mathrm{}}{lim}v_3\mathrm{\Sigma }_1\mathrm{\Sigma }_3=0,\underset{\tau \mathrm{}}{lim}(1v_3^2)\mathrm{\Omega }=0,$$
(B.3)
and if $`\frac{14}{9}<\gamma <2`$, then

$$\underset{\tau +\mathrm{}}{lim}(1v_3^2)\mathrm{\Omega }=0,\underset{\tau \mathrm{}}{lim}v_3\mathrm{\Sigma }_1\mathrm{\Sigma }_3=0.$$
(B.4)
* For some values of $`\gamma `$ the tilt variable $`v_3`$ is monotone. If $`\gamma `$ satisfies $`0<\gamma <\frac{2}{3}`$, the factor $`3\gamma 4\mathrm{\Sigma }_++\sqrt{3}\mathrm{\Sigma }_{}`$ in the DE for $`v_3`$ is negative (see (2.3)), and hence $`v_3`$ is monotone decreasing into the future. It thus follows that if $`0<\gamma <\frac{2}{3}`$,

$$\underset{\tau +\mathrm{}}{lim}v_3=0,\underset{\tau \mathrm{}}{lim}v_3=1.$$
(B.5)
If $`\gamma =2`$, the DE for $`v_3`$ simplifies to

$$v_3^{}=(2\mathrm{\Sigma }_++\sqrt{3}\mathrm{\Sigma }_{})v_3.$$
(B.6)
Since $`2\mathrm{\Sigma }_++\sqrt{3}\mathrm{\Sigma }_{}>0`$ unless $`(\mathrm{\Sigma }_+,\mathrm{\Sigma }_{})=(\frac{1}{2},\frac{\sqrt{3}}{2})`$, it follows that $`v_3`$ is monotone increasing into the future. Thus, if $`\gamma =2`$,

$$\underset{\tau \mathrm{}}{lim}v_3=0.$$
(B.7)
The DE (B.6) also implies that if $`\gamma =2`$, the extreme tilt set given by $`v_3=1`$ is not an invariant set, with the result that orbits can pass from the region $`v_3<1`$ to the region $`v_3>1`$ in state space.
* It follows from (2.3) and (2.4) that

$$(\beta \mathrm{\Omega })^{}=[2q(3\gamma 2)](\beta \mathrm{\Omega }),$$
(B.8)
where

$$\beta =G^1(1v_3^2)^{\frac{1}{2}(2\gamma )}.$$

The expression (A.31) for $`q`$ can be rearranged to give

$$2q(3\gamma 2)=3(2\gamma )\mathrm{\Sigma }^2+(23\gamma )K+\gamma (43\gamma )G^1v^2\mathrm{\Omega },$$

where $`K`$ is given by (A.39). Thus, if $`\gamma `$ satisfies $`0<\gamma <\frac{2}{3}`$, $`\beta \mathrm{\Omega }`$ is monotone increasing, while if $`\gamma =2`$, $`\beta \mathrm{\Omega }`$ is monotone decreasing. It follows that if $`0<\gamma <\frac{2}{3}`$, then

$$\underset{\tau +\mathrm{}}{lim}\mathrm{\Omega }=1,\underset{\tau \mathrm{}}{lim}\mathrm{\Omega }=0.$$
(B.9)
We note that the limits (B.3), (B.4), (B.5), (B.7) and (B.9) provide support for the claims made in Sections 3 and 4 concerning the asymptotic behaviour.

## Appendix C

In this appendix, we give the equations relating the kinematical quantities of the fluid congruence in terms of the dimensionless commutation functions (A.2) and (A.3) associated with the normal congruence. The fluid kinematical quantities, namely, the acceleration $`\dot{u}_a`$, the vorticity vector $`\omega _a`$ and the rate of shear tensor $`\sigma _{ab}^{\mathrm{fluid}}`$, are normalized using the Hubble scalar $`H_{\mathrm{fluid}}`$ of the fluid congruence:

$$\dot{U}_a=\frac{\dot{u}_a}{H_{\mathrm{fluid}}},W_a=\frac{\omega _a}{H_{\mathrm{fluid}}},\mathrm{\Sigma }_{ab}^{\mathrm{fl}}=\frac{\sigma _{ab}^{\mathrm{fluid}}}{H_{\mathrm{fluid}}}.$$

We now give the desired relations, which are obtained by adapting the results of  (see Sections 1 and 2)141414We thank W.C. Lim for doing these calculations.. Note the relation between our tilt vector $`v_a`$ and their vector $`\stackrel{~}{c}_a`$:

$$v_a=(\mathrm{tanh}\beta )\stackrel{~}{c}_a,$$

where

$$\mathrm{tanh}\beta =v.$$

Hubble scalar:

$$H_{\mathrm{fluid}}=BH,$$
(C.1)
where

$$B=\mathrm{cosh}\beta \frac{\left[1\frac{1}{3}(v^2+\mathrm{\Sigma }_{\alpha \beta }v^\alpha v^\beta +2A_\alpha v^\alpha )\right]}{1(\gamma 1)v^2}.$$
(C.2)
Acceleration:

$$\dot{U}_\alpha =3(\gamma 1)\mathrm{cosh}\beta v_\alpha ,\dot{U}_0=v^\alpha \dot{U}_\alpha .$$

Vorticity:

$$W_\alpha =\frac{1}{2B}(N_\alpha ^\beta v_\beta +\epsilon _\alpha ^{\mu \nu }v_\mu A_\nu +\mathrm{cosh}^2\beta N^{\mu \nu }v_\mu v_\nu v_\alpha ),W_0=v^\alpha W_\alpha .$$

Shear:

$`\mathrm{\Sigma }_{\alpha \beta }^f\mathrm{}+\delta _{\alpha \beta }`$ $`={\displaystyle \frac{\mathrm{cosh}\beta }{B}}(\mathrm{\Sigma }_{\alpha \beta }+\delta _{\alpha \beta })\mathrm{cosh}^2\beta (43\gamma )v_\alpha v_\beta `$
$`+{\displaystyle \frac{\mathrm{cosh}\beta }{B}}[N_{(\alpha }^\nu \epsilon _{\beta )\mu \nu }v^\mu +A_{(\alpha }v_{\beta )}A_\mu v^\mu \delta _{\alpha \beta }],`$
$`\mathrm{\Sigma }_{0\alpha }^f\mathrm{}`$ $`=\mathrm{\Sigma }_{\alpha \beta }^f\mathrm{}v^\beta ,\mathrm{\Sigma }_{00}^f\mathrm{}=\mathrm{\Sigma }_{\alpha \beta }^f\mathrm{}v^\alpha v^\beta .`$
## References

* MacCallum, M. A. H. (1979) General Relativity eds. S. W. Hawking and W. Israel (Cambridge University Press).
* King, A. R. and Ellis, G. F. R. (1973) Commun. Math. Phys. 31 209-42.
* MacCallum, M. A. H. (1973) Cargése Lectures in Physics Volume 6 ed. E Schatzman (Gordon & Breach).
* Wainwright, J. and Ellis, G. F. R. (1997) Dynamical Systems in Cosmology (Cambridge University Press).
* Wainwright, J. (2000) GRG, to appear.
* Rosquist, K. and Jantzen, R. T. (1988) Phys. Rep. 166 89-124.
* Bogoyavlensky, O. I. (1985) Methods in the Qualitative Theory of Dynamical Systems in Astrophysics and Gas Dynamics Springer-Verlag.
* Peresetsky, A. A. (1985) Topics in Modern Mathematics Petroviskii Seminar No. 5, ed. O. A. Oleinik (Consultants Bureau).
* Collins, C. B. and Ellis, G. F. R. (1979) Phys. Rep. 56 65-105.
* Hewitt, C. G. and Wainwright, J. (1992) Phys. Rev. D 46 4242-52.
* Harnett, D. (1996) M.Math Thesis, University of Waterloo.
* Ellis, G. F. R. and MacCallum, M. A. H. (1969) Commun. Math. Phys. 12 108-41.
* Hewitt, C. G. (1991) Class. Quantum Grav. 8 L109-14.
* Hsu, L. and Wainwright, J. (1986) Class Quantum Grav. 3 1105-24.
* Ma, P. K-H. and Wainwright, J. (1992) Relativity Today ed. Z Perjés (Nova Science Publishers); reprinted in Deterministic Chaos in General Relativity, edited by Hobill D, Burd A, Coley A, Plenum Press, 1994.
* Rendall, A. D. (1997) Class. Quantum Grav. 14 2341-2356.
* Ringström, H. (2000) Class. Quantum Grav. 17 713-731.
* Misner, C.W. (1969) Phys. Rev. Lett. 22, 1071-74.
* Misner, C.W. (1970) in Relativity, Proceedings of the Relativity Conference in the Midwest, eds. M. Carmeli, S.I. Fickler and L. Witten, Plenum Press.
* Belinskii, V. A., Khalatnikov, I. M. and Lifschitz, E. M. (1970) Adv. Phys. 19, 525-73.
* Jantzen, R.T. (1987) in Gamow Cosmology, eds. R. Ruffini and F. Melchiorri, Proc. Int. School of Physics E. Fermi, Course LXXXVI, North Holland.
* Belinskii, V. A., Khalatnikov, I. M. and Lifschitz, E. M. (1982) Adv. Phys. 31. 639-67.
* Uggla, C. (1989) Class. Quantum Grav. 6, 383-96.
* Leblanc, V. G., Kerr, D. and Wainwright, J. (1995) Class Quantum Grav. 12, 513-41.
* Leblanc, V. G. (1997) Class. Quantum Grav. 14, 2281-2301.
* Leblanc, V. G. (1998) Class. Quantum Grav. 15, 1607-26.
* Weaver, M. (2000) Class. Quantum Grav. 17, 421-434.
* van Elst, H. and Uggla, C. (1997) Class. Quantum Grav. 14, 2673-2695.


text:
# Coherent oscillation in a linear quantum system coupled to a thermal bath.

\[

## Abstract

We consider the time development of the density matrix for a system coupled to a thermal bath, in models that go beyond the standard two-level systems through addition of an energy excitation degree of freedom and through the possibility of the replacement of the spin algebra by a more complex algebra. We find conditions under which increasing the coupling to the bath above a certain level decreases the rate of entropy production, and in which the limiting behavior is a dissipationless sinusoidal oscillation that could be interpreted as the synchronization of many modes of the uncoupled system.

\]

Introduction. We explore some extensions of the problem of a two level system in a thermal bath  and find some new qualitative behavior that can be provided by these extensions. Two specific examples are the symmetrical double well in an environment with temperature on the order of the (single-well) level spacing, and the propagation of oscillating neutrinos in hot, dense media. In both cases we find behavior in which a weak coupling of the system to the bath induces some rate of entropy gain, but where an increase in the coupling beyond a certain value reduces the rate of entropy production, reducing it to zero in a limiting case. The limiting behavior can correspond to the “frozen” conditions of ref., and many other papers, in which the implicit continuous measurement of a quantity by the system-bath interaction can be said to have frozen that quantity in time. However our main focus in the present note is on circumstances that give instead a kind of dissipationless collective oscillation, which gives the appearance of a synchronization of modes that had different frequencies in the absence of the coupling to the bath. We stress, however, that there appears to be little direct connection of these results to the large literature of synchronized behavior in non-linear systems  or of stochastic resonances in forced systems .

The system itself has two coordinates, which we characterize as “horizontal” and “vertical”. The horizontal component is the “spin” coordinate of a two state system, or in cases that may be of interest in neutrino physics, a “flavor” coordinate for a 3 or 4-state system. We index these states with Greek letters. The vertical variable corresponds to a set of energy levels, $`E_i`$. The system basis states are denoted $`|E_i,\alpha `$. The dynamics of the system uncoupled to the bath is governed by the Hamiltonian,

$$H_0^{\mathrm{sys}}=\underset{i,\alpha ,\beta }{}\{E_i\delta _{\alpha ,\beta }+\lambda _{\alpha ,\beta }(E_i)\}|E_i,\alpha E_i,\beta |.$$
(1)
We take $`\lambda _{\alpha ,\beta }(E_i)E_k`$ for all values $`i,k`$. We consider a density matrix for the system $`\rho ^{\mathrm{sys}}(t)=\mathrm{Tr}_{\mathrm{bath}}[\rho ^{\mathrm{s}+\mathrm{b}}(t)]`$, where $`\rho ^{\mathrm{s}+\mathrm{b}}(t)`$ is the complete density matrix of the model. At an initial time t=0 we choose a form diagonal in the vertical indices,

$$\rho ^{\mathrm{sys}}(t=0)=\underset{i,\alpha ,\beta }{}\rho _{\alpha ,\beta }(E_i,t=0)|E_i,\alpha E_i,\beta |.$$
(2)
For these initial conditions, one of the results of a master equation derivation is that for times $`t(\mathrm{\Delta }E)^1`$, where $`\mathrm{\Delta }E`$ is a measure of the (single-well) energy spacings, the operator $`\rho ^{sys}`$ remains so nearly diagonal in the vertical (energy) indices that we can continue to describe the system by a vector in the energy space, $`\rho (E_i,t)_{\alpha ,\beta }`$, replacing ($`t=0`$) by ($`t`$) in (2). Henceforth we suppress the indices $`(\alpha ,\beta )`$ in the horizontal space, in which $`\rho (E_i,t)`$ remains an operator. For the bath-system coupling, $`H_{sb}`$, we take the dependence on the horizontal (spin-flavor variables) to be in a factorized form, $`H_{sb}=\zeta V`$, where $`\zeta `$ operates on the horizontal indices and $`V`$ depends only on the vertical coordinates and the coordinates of the bath. This coupling is also taken to be weak, in a way that allows the problem to be solved, although it can at the same time be strong in another sense, for the realization of the limits described above. At the end of this note we show that for this form of coupling the Bloch equations, generalized to include the vertical structure, are of the form,

$`{\displaystyle \frac{d}{dt}}\rho (E_i,t)`$ $`=`$ $`i[\lambda (E_i),\rho (E_i,t)]+{\displaystyle \underset{E_j}{}}\zeta \rho (E_j,t)\zeta \mathrm{\Gamma }(E_j,E_i)`$ (3)
$``$ $`{\displaystyle \frac{1}{2}}(\zeta ^2\rho (E_i,t)+\rho (E_i,t)\zeta ^2){\displaystyle }_{E_j}\mathrm{\Gamma }(E_i,E_j).`$ (4)
All of the elements in (4) are matrices in the horizontal space except for the functions $`\mathrm{\Gamma }`$. These functions obey the relation $`\mathrm{\Gamma }(E_j,E_i)=\mathrm{exp}[(E_jE_i)/T]\mathrm{\Gamma }(E_i,E_j)`$ in consequence of the thermal equilibrium of the bath. For the neutrino application, and with the right identifications, (4) is the “Quantum Kinetic Equation” derived by McKellar and Thomson . Also Loreti and Balantekin  have used a special case of (4) to discuss the energy conserving case of propagation of neutrinos through a medium with a (spatially) fluctuating electron density, in which case the collision term becomes a pure double commutator in the horizontal (flavor) space. For the double well application, we sketch at the end of this note a derivation in the context needed for our present results.

A. Symmetrical double well

The horizontal space is 2x2, with $`\sigma _3=1`$ for the particle to be on the left side and $`\sigma _3=1`$ for the right side. The barrier height is large compared to $`\mathrm{\Delta }E`$ where $`\mathrm{\Delta }E`$ is a measure of the energy spacing for the single well. The vertical states are the single well energy excitations. We choose $`\lambda _{\alpha ,\beta }(E_i)=[\sigma _1]_{\alpha ,\beta }g(E_i)`$. The energy splittings $`g(E_i)`$, which are essentially tunneling rates, have strong $`E_i`$ dependence. For the “spin” dependence of the coupling to the bath we take, $`\zeta =1+b\sigma _3`$. In the case $`b=0`$ the transition matrix elements of the operator $`V`$ for $`E_iE_j`$ in the left-hand well, considered by itself, are the same as those for the right-hand well, as would be the case in the dipole approximation of the interaction with a radiation field. For the case $`b=1`$ only the amplitude on the left-hand side of the well interacts with the bath.

The coefficient $`\mathrm{\Gamma }(E_i,E_j)`$ is the (single well) rate for a state $`E_i`$ to make a transition to state $`E_j`$. When $`T\mathrm{\Delta }E`$, where $`\mathrm{\Delta }E`$ is the scale of the vertical splittings, the bath induces transitions among the vertical states, and in this case our models differ substantively from models that do not have the vertical structure.

B. Neutrinos

The question is the time dependence of the flavor density matrix for neutrinos that oscillate and are also scattered as they traverse a dense medium, like that of the supernova core, or the early universe at a temperature greater than 1 MeV. Here the $`\lambda _{\alpha ,\beta }`$ parameters are the combination of the mass$`{}_{}{}^{2}/(2E)`$ matrix and the “index of refraction” terms, the latter energy independent and of order $`G_W`$. The functions, $`\mathrm{\Gamma }(E_i,E_j)`$, of order $`G_W^2`$, are the rates of transition for a neutrino to go from state $`E_i`$ to state $`E_j`$, as calculated with the matrix $`\zeta `$ set equal to unity. The choice $`b=0`$ would apply to $`\nu _\mu ,\nu _\tau `$ mixing in a flavor-blind medium, that is, a medium in which $`\nu _\mu `$ and $`\nu _\tau `$ have the same interactions with the environment. The choice $`b=1`$ corresponds to “active-sterile” scenarios in which one of the neutrinos does not interact with the environment at all.

Results for $`b=0`$. As an example we use a double square well with infinite barriers to the left and right of $`x=\pm 8a`$, and with a central barrier of width $`2a`$ and height $`U_0`$. We choose $`U_0`$ and the particle mass such that twenty states are bound with energies less than $`U_0`$. The low-lying states are thus very nearly the symmetrical and anti-symmetrical combinations of separate states in an infinite well of width $`7a`$. The energy splittings, $`g(E_i)`$ are easy to calculate in the limit of small tunneling. For the bath we take bosons with $`\omega =vk`$, equal energy spacing $`2\pi v/L`$ where $`L\mathrm{}`$, and with creation and annihilation operators, $`a_n^{},a_n`$. We take $`H_0=_na_n^{}a_n\omega _n+H_0^{(\mathrm{sys})}`$, and in the system-bath interaction, $`H_{sb}=V\zeta `$, we take the vertical operator $`V`$ to have off-diagonal matrix elements given in dipole expansion as ***There are diagonal matrix elements of the $`kx`$ operator as well, since the wells are not centered at $`x=0`$. In principle, diagonal matrix elements contribute to (4) in linear order, as a modification of the oscillation matrices $`\lambda `$, but in the present case these contributions cancel when $`k`$ is added to $`k`$.,

$$V_{j,k}=\frac{q}{L}\underset{n}{}[(a_n^{}+a_n)\delta _{j,k}+ik_nx_{j,k}(a_n^{}a_n)].$$
(5)
where $`q`$ is a coupling strength. We calculate the off-diagonal dipole matrix elements, $`x_{j,k}`$ using the infinite well wave-functions and define $`c_{j,k}=qk_nx_{j,k}`$, where the later formalism will justify using (near) energy conservation to express $`k_n`$ as $`\pm (E_jE_k)/v`$. The coefficients in the bath interaction part of (4) are then,

$$\mathrm{\Gamma }(E_j,E_k)=2\pi c_{j,k}c_{k,j}[\frac{\theta (E_jE_k)}{e^{(E_jE_k)/T}1}+\frac{\theta (E_kE_j)}{1e^{(E_jE_k)/T}}].$$
(6)
Incorporating these results and solving (4) for several values of $`q`$, we show in Fig.1 the probability that the particle will be on the left side at time $`t`$, with the initial condition of being on the left and in a thermal distribution at the bath temperature, which we take as $`T=5`$ in units of the single-well ground state energy.

For very strong coupling, the probability becomes simply periodic, as one mode would be in the absence of dissipation. This limit can be solved analytically, with the result that the oscillation frequency is the thermal average of the oscillation frequencies for the separate modes . In fig 2 we show the entropy, as defined by $`S=_j\mathrm{Tr}[\rho (E_j)\mathrm{log}(\rho (E_j))]`$, where the trace is over the $`2\times 2`$ space. We note that as the coupling is increased from small to moderate values the rate of entropy increase rises, but that as the coupling becomes very large the rate of entropy increase goes to zero.

We can generalize (4) to apply to a system in which the wells contain N fermions with no mutual interactions. For example, after making the appropriate modifications to the master equation for Fermi statistics, and using the same values of parameters as before, but taking $`N=5`$, all initially on the left and in thermal equilibrium, we again look at the strong coupling limit, getting results very similar to those shown if fig.1 for the density matrix of a single particle. It is surprising that one can prepare a state in which, if we looked only at number density, a set of Fermi particles could appear to be oscillating, all in phase, between the left and the right, as would an assemblage of bosons all in an “average energy” state.

For the case of neutrino mixing in a noisy environment, where the vertical space is a continuum, the generalized Bloch equations, (4) are integro-differential equations. But we find that when they are discretized with, say, around 30 energy levels, they are easily solved, with the results insensitive to the mesh size. The results are qualitatively the same as those discussed above.

There is an amusing generalization of the above considerations in the extension to a double well with a time dependent bias, $`\delta H_0^{\mathrm{sys}}=ϵ(t)\sigma _3`$. Consider a case in which we begin with $`ϵ>0`$ and $`ϵ|g(E_i)|`$ for some set of low-lying states (but $`ϵ\mathrm{\Delta }E`$). A group of particles initially localized on the left, distributed among the low-lying states, and not interacting with the medium, will stay largely on the left if the bias is maintained at the initial value. Now let the bias be slowly lowered until it is equal and opposite to the original bias. The particles will efficiently be carried to the right hand well. This is the analogue of the MSW transition in neutrino physics. It is a rather diffuse transition, since the different energy levels go through their avoided crossings at different times. With a moderate coupling to the bath, the system instead settles down rather quickly to half on the right and half on the left, with no interesting persistent order. However, if we increase the bath coupling sufficiently we regain the shift from left to right. The transition is now sharper than it was in the case without the bath. Detailed analysis of this type of behavior for the neutrino case is given in .

Results for $`b0`$. We briefly summarize some examples. In the case plotted in fig. 1, for the strongest coupling, we find little change for $`b=0.001`$. But $`b=0.005`$ starts to poison the oscillation, damping it by a factor of 20% in 5 periods. As we move to much larger $`b`$, the motion becomes more frozen, so that when, e.g., $`b=0.5`$, 98% of the probability has remained on the left hand side over a time span of five average periods of the $`b=0`$ system. A popular interpretation of frozen behavior is that the interactions with the bath are continuously measuring the value of $`\sigma _3`$ and projecting the wave function onto the initial state. Although there is no measurement implied by our equations, the results are consistent with this terminology, since in the case with $`b=0`$ the system-bath interaction does not distinguish the right-hand well from the left-hand well and it is for this case that we get the limiting behavior of undamped simple oscillations with a thermal averaged frequency.

Derivation of the master equation. We focus on the aspects that give rise to the placing of the horizontal matrices in (4) and the detailed form of the RHS of (6). The methods are closely related to those given in  and . We introduce an interaction picture based on the above division of the Hamiltonian. The system-bath coupling in this picture is denoted by $`H_I(t)`$, with vertical matrix elements,

$$[H_I(t)]_{j,k}=e^{i\lambda (E_j)t}\zeta e^{i\lambda (E_k)t}V_{j,k}^Ie^{i(E_jE_k)t}.$$
(7)
The leading terms in the system+bath density matrix, $`\rho _I^{\mathrm{s}+\mathrm{b}}`$, in this interaction picture are generated by the integral equation,

$`\rho _I^{\mathrm{s}+\mathrm{b}}(t)`$ $`=`$ $`\rho ^{\mathrm{s}+\mathrm{b}}(t=0)`$ (8)
$``$ $`{\displaystyle _0^t}dt_1{\displaystyle _0^{t_1}}dt_2[[\rho _I^{\mathrm{s}+\mathrm{b}}(t_2),H_I(t_2)],H_I(t_1)].`$ (9)
In the perturbation expansion of (9) we find a piece with an additional power of the quantities $`c_{j,k}^2[Et,E\lambda ^1]`$ in each higher order. This provides our definition of leading terms and the explanation for how a tiny $`c_{j,k}^2`$ can lead to appreciable effects over long periods of time. These terms are generated by the terms in which each successive new pair of $`H_I`$’s in the iteration solution puts the state back to the same vertical (and bath) state. These observations provide the basis for having omitted, in writing (9), all odd terms in $`H_I`$, and they also dictate the form of the density matrix to be used in solving the equation,

$$\rho _I^{\mathrm{s}+\mathrm{b}}=\rho ^{\mathrm{bath}}\underset{j}{}\rho _I(E_j,t)|jj|.$$
(10)
where,

$$\rho ^{\mathrm{bath}}=[\underset{i}{}(1e^{\omega _i/T})]\mathrm{exp}[(\underset{n}{}a_n^{}a_n\omega _n)/T],$$
(11)
and where,

$$\rho _I(E_j,t)=e^{i\lambda (E_j)t}\rho (E_j,t)e^{i\lambda (E_j)t}.$$
(12)
We substitute (10) into (9), take the bath trace, and use leading order expressions for each term in the double commutator, obtaining, for example,

$`\mathrm{Tr}`$ $`{}_{\mathrm{bath}}{}^{}[j|H_I(t_1)\rho _I^{\mathrm{s}+\mathrm{b}}(t_2)H_I(t_2)|j]=\delta (t_1t_2)`$ (13)
$`\times `$ $`{\displaystyle \underset{k}{}}\mathrm{\Gamma }(E_k,E_j)e^{i\lambda (E_j)t_1}\zeta e^{i\lambda (E_k)t_1}`$ (14)
$`\times `$ $`\rho _I(E_k,t_1)e^{i\lambda (E_k)t_1}\zeta e^{i\lambda (E_j)t_1}.`$ (15)
Similarly, we have

$`\mathrm{Tr}`$ $`{}_{\mathrm{bath}}{}^{}[j|H_I(t_1)H_I(t_2)\rho _I^{\mathrm{s}+\mathrm{b}}(t_2)|j]=\delta (t_1t_2)`$ (16)
$`\times `$ $`e^{i\lambda (E_j)t_1}\zeta ^2e^{i\lambda (E_j)t_1}{\displaystyle \underset{k}{}}\mathrm{\Gamma }(E_j,E_k)\rho _I(E_j,t_1).`$ (17)
The key to separation of the leading terms in the above is performing the sum over the modes of the scalar field first. After making the transition to the continuum, we make replacements of the form,

$`{\displaystyle _{\omega _{min}}^{\mathrm{}}}`$ $`d\omega f(\omega )e^{i(\pm E_jE_k\pm \omega )(t_1t_2)}`$ (18)
$``$ $`2\pi f[(E_kE_j)]\delta (t_1t_2)\theta (E_kE_j).`$ (19)
where neither the introduced infrared cut-off, $`\omega _{min}`$ nor the residual terms contribute to leading order, as defined above.

Using (15), (17) and their counterparts for the two other orderings in (6) , doing the $`t_2`$ integral (with the delta function symmetrically smeared, as can be justified by more accurate integrations) and differentiating with respect to $`t`$ gives,

$`{\displaystyle \frac{d}{dt}}`$ $`\rho _I(E_j,t)={\displaystyle \underset{E_k}{}}e^{i\lambda (E_j)t}\zeta e^{i\lambda (E_k)t}\rho _I(E_k,t)`$ (20)
$`\times `$ $`e^{i\lambda (E_k)t}\zeta e^{i\lambda (E_j)t}\mathrm{\Gamma }(E_k,E_j)`$ (21)
$``$ $`{\displaystyle \frac{1}{2}}(e^{i\lambda (E_j)t}\zeta ^2e^{i\lambda (E_j)t}\rho _I(E_j,t)`$ (22)
$`+`$ $`\rho _I(E_j,t)e^{i\lambda (E_j)t}\zeta ^2e^{i\lambda (E_j)t}){\displaystyle }_{E_k}\mathrm{\Gamma }(E_j,E_k).`$ (23)
Using (12) we now regain (4).

To summarize, we have developed the elements of the master equation for a bath-coupled system that is a generalization in two directions of the system of the ground states of a double-well. On the basis of this equation we have found new instances in which strong bath-system interactions can lead to more, rather than less, orderly behavior. We have not begun to address the question of identifying a realistic quantum well system in which one could make an experimental application of the conclusions. We should emphasize, however, that our limit of “large coupling” to the bath is the limit in which $`\mathrm{\Gamma }\lambda `$; it is achievable, in principle, for any strength of the system-bath coupling (whether the mechanism be photons, phonons or collisions with surrounding particles), simply by widening the barrier between the wells.

We thank Wick Haxton for hospitality at the National Institute for Nuclear Theory, where conversations among the three us began. The work of R.F.S. was supported, in part, by NSF grant PHY-9900544, N.F.B. by the Commonwealth of Australia, and R.R.V. by the Australian Research Council. R.F.S. would like to thank B. McKellar for kind hospitality during a visit to the University of Melbourne. R.F.S. thanks Horia Metiu for an enlightening conversation.