On the oscillations and future asymptotics of locally rotationally symmetric Bianchi type III cosmologies with a massive scalar field

We analyse spatially homogenous cosmological models of locally rotationally symmetric Bianchi type III with a massive scalar field as matter model. Our main result concerns the future asymptotics of these spacetimes and gives the dominant time behaviour of the metric for all solutions for late times. This metric is forever expanding in all directions, however in one spatial direction only at a logarithmic rate, while at a power-law rate in the other two. Although the energy density goes to zero, it is matter dominated in the sense that the metric components differ qualitatively from the corresponding vacuum future asymptotics. Our results rely on a conjecture for which we give strong analytical and numerical support. For this we apply methods from the theory of averaging in nonlinear dynamical systems. This allows us to control the oscillations entering the system through the scalar field by the Klein-Gordon equation in a perturbative approach.


Introduction
Spatially homogenous cosmologies have been investigated in vacuum [19] and with various matter models such as perfect fluids [5,23], collisionless kinetic gases (Vlasov matter) [1,16,17] (recently [2,6,11]) magnetic fields or scalar fields [5,17] (recently [8,10,25]), in particular with a focus on their asymptotic dynamics and stability. A fruitful approach has been to formulate the respective systems in terms of expansion normalised (or Hubble normalised) dimensionless quantities. Typically the Raychaudhuri equation decouples, which is the evolution equation for the Hubble scalar. The latter describes the overall rate of spatial expansion. The asymptotics of the remaining reduced system is then typically given by equilibrium points and often can be determined by a dynamical systems analysis; cf [5,14,23].
With the present work we contribute to spatially homogenous scalar field cosmology. In the literature a lot of attention has been given to scalar fields with an exponential potential. One of the reasons is that in this case the above approach is applicable since the Raychaudhuri equation decouples due to the symmetry that the exponential is its own derivative; cf [5,Sec IV,A]. Scalar fields with harmonic potentials, such as massive scalar fields which satisfy the Klein-Gordon equation, have also been investigated [5,Sec III,E]. Here however the Raychaudhuri equation fails to decouple, which is why the above approach is not applicable and the resulting system proves to be difficult to analyse in a dynamical systems approach. On the other hand, other approaches, such as a local stability analysis, are also difficult to apply due to oscillations entering the system via the Klein-Gordon equation. Nevertheless there have been successful enquiries going this route, such as the recent works on various nonlinear stability problems for the Einstein-Klein-Gordon system concerning Minkowski spacetime [10,12,26] and the Milne model [8,25]. While those results concern vacuum dominated future asymptotics, in the present work we consider a class of models where the future is matter dominated. By this we mean that the dominant behaviour of the metric for large times differs qualitatively from that of the corresponding vacuum spacetime.
In this paper we introduce a new approach to the field and apply techniques from the theory of averaging in non-linear dynamical systems [21]. The core idea is to construct a time average of the original system, in the solution of which the oscillations are smoothed out. The oscillations are thus viewed as perturbations and it turns out that the Hubble scalar plays the role of a time dependent perturbation parameter which controls the magnitude of error between the full and the averaged solutions. Most importantly, the Raychaudhuri equation of the averaged system decouples, and the remaining reduced system can be analysed with the well known dynamical systems techniques mentioned above. Consequently we reduce the analysis of the original system to that of a corresponding averaged system.
In general the averaged solution represents only an approximation to the full solution, and error estimates in terms of the perturbation parameter are typically only valid for limited time. However in cases where either system is attracted to an equilibrium point the validity of such error estimates can be extended to all times; cf [21,Chap 5]. It turns out that we are looking at such a case, and since we can show that our perturbation parameter, the Hubble scalar, goes to zero the full and the averaged solutions converge in the limit when time goes to infinity.
The spatially homogenous cosmologies we are considering are those of locally rotationally symmetric (LRS) Bianchi type III, which in different terminology are also referred to as LRS Bianchi type VI −1 ; cf [16,20,23]. As matter model we consider a massive scalar field. Hence we are dealing with the Einstein-Klein-Gordon system within the LRS Bianchi III symmetry class. Our main results consider the future asymptotics of these systems. We express these in terms of the dominant behaviour of physical and geometric quantities for large times, such as the shear parameter and the energy density (Theorem 1) as well as the metric components (Theorem 2). The premises of the main theorems contain the assumption that a conjecture holds (Conjecture 1). For the latter we give strong analytical and numerical support. A closely related statement is proven in [7]. Apart from supporting the assumption and conjecture, the numerics also convincingly demonstrates agreement with all our analytical results.
Finally we mention that one of our motivations to consider LRS Bianchi III Einstein-Klein-Gordon cosmologies was the result of [15] which concerns Einstein-Vlasov cosmologies with massive particles within the same symmetry class. The finding was a forever expanding matter dominated future asymptotics with one scale factor increasing linear with time, and the other one increasing slower than any power law. We were interested if the case of a massive scalar field would yield the same future asymptotics. Our results gives an affirmative answer to this question.
We begin with a brief introduction to periodic averaging in Section 2. In Section 3 we introduce the LRS Bianchi III Einstein-Klein-Gordon system and formulate an averaging conjecture for it; cf Conjecture 1. The conjecture gives error estimates comparing the full and a corresponding averaged solution. We then provide analytical support for this conjecture in Section 4. Under the assumption that the conjecture holds, we then derive our main results in Section 5: the future asymptotics of LRS Bianchi III Einstein-Klein-Gordon cosmologies; cf theorems 1 and 2. Section 6 is then concerned with numerical support for our results and assumptions. Finally, we conclude with a discussion of our results and an outlook on open problems and further applicability of the averaging techniques introduced here to the field of mathematical cosmology in Section 7.
In particular for Section 4 we assume some familiarity with standard methods of dynamical systems theory; cf [14]. However, because it has not been used to a wider extend in the literature of mathematical cosmology, we put a brief introduction to centre manifold analysis in Appendix B, following [4,Chap 1]. We apply these tools in Subsection 5.3. Appendix A contains two technical definitions which we refer to in Section 2.

Periodic averaging and the Van der Pol equation
We restrict the following outline to a brief summary of those aspects of the theory of averaging which are relevant to the analysis of the present work, and follow Sanders, Verhulst and Murdock [21], in particular its chapters 2 and 5.
In Subsection 2.1 we outline the basic idea of periodic averaging. We specify what we mean by a problem in standard form and for such give two theorems which estimate the error between the full and the corresponding averaged solution. The theory is then applied to the Van der Pol equation in Subsection 2.2. This is not only an illustrative example, but is also of direct relevance to our problem since the Klein-Gordon equation, though coupled to the Einstein equations, can be considered as a specific example of a Van der Pol equation.
2.1. Periodic averaging. The theory of averaging studies initial value problems of the general formẋ with x, f (x, t, ) ∈ R n , where plays the role of a, usually small, perturbation parameter. Typically one would then perform a Taylor expansion of f in around = 0. For the simplest form of averaging, periodic averaging, the zeroth order term usually vanishes, and one is typically looking at problems of the standard forṁ with f 1 and f [2] T -periodic in t. The exponents correspond to the respective perturbative order, and the square bracket marks the remainder of the series; cf [21,p 13,Notation 1.5.2].
To first order, the theory is then concerned with the question to what degree solutions of (1) can be approximated by the solutions of an associated averaged systemż The basic result is given by the following theorem: Lemma 1 ([21, p 31, Thm 2.8.1]). Let f 1 be Lipschitz continuous, let f [2] be continuous, and let 0 , D, L be as in Definition A.1. Then there exists a constant c > 0 such that ||x(t, ) − z(t, )|| < c for 0 ≤ ≤ 0 and 0 ≤ t ≤ L/ , and where || . || denotes the norm ||u|| : In other words, the error one makes by approximating the full system (1) by the averaged system (2) will be of order on timescales of order −1 . When the solutions of the full or averaged system are attracted by an asymptotically stable critical point, the domain of approximation might be extendable to all times; cf [21,Chap 5]. For instance: Lemma 2 (Eckhaus/Sanchez-Palencia [21, p 101, Thm 5.5.1]). Consider the initial value problemẋ with a, x ∈ D ⊂ R n and f 1 T -periodic in t. Suppose f 1 is a KBM-vector field ( Definition A.2) producing the averaged equatioṅ where z = 0 is an asymptotically stable critical point in the linear approximation, f 1 is continuously differentiable with respect to z in D and has a domain of attraction D o ⊂ D. Then for any compact K ⊂ D o and for all a ∈ K

The Van der Pol equation.
A demonstrative example of first order periodic averaging, which is relevant to the analysis of the present paper, is given by the class of Van der Pol equationsφ with g sufficiently smooth. (4) describes a harmonic oscillator with generally nonlinear feedback and damping. An amplitude-phase (variation of constants) transformation yields a system in standard form (1) where the arguments of g are understood to be substituted by the transformation. Consequently, the averaged system (2) is given by 1 2π 0 cos(s − ϕ) g r sin(s − ϕ), r cos(s − ϕ) ds (7) and f 1 ϕ (r) defined analogously. Note that f 1 is independent of ϕ. Specialising now to the specific Van der Pol equation with g(φ,φ) = (1 − φ 2 )φ as an illustrative example we obtain the averaged system ṙ ϕ = By Lemma 1 we know that the error between [r, ϕ] T and [r, ϕ] T will be of order on timescales of order −1 . Since f 1 is independent of ϕ we restrict to the decoupled equationṙ, 2 which has the two equilibrium points r = 0 and r = 2. The equilibrium 1 We follow here the notation of [21] and put bars above the variables to indicate that they belong to the averaged system. 2 In full rigorosity one has to justify this, which can be done by using a coordinate transformation φ = r sin(θ),φ = r cos(θ) instead, which yields a single averaged equation dr/dθ = f 1 1 (r), with f 1 1 (r) of the same form as here; cf [ In this section we lay out the basic setup for our analysis. We introduce the LRS Bianchi III Einstein-Klein-Gordon system in Subsection 3.1. In Subsection 3.2 we then bring this system into a form closely resembling the standard form of periodic averaging problems discussed in the previous section. We call this the quasi-standard form and formulate an averaging conjecture for it in Subsection 3.3 which seeks to widen the applicability of the lemmas of the previous section to our case. At the end of the section we also derive the monotonicity and future asymptotics for one of our dynamical quantities, the Hubble scalar; cf Lemma 3. Finally we state the implications of the latter in conjunction with our averaging conjecture; cf Proposition 1.
3.1. The basic system. An LRS Bianchi III metric can be written in the form where g H 2 denotes the 2-metric of negative constant curvature on hyperbolic 2space; cf [3,App B.4]. We adopt the frame choice and variables of [18,Sec 2] and [15,Sec 3] and the resulting dynamical systems formulation of the Einstein equations for the metric (9) and an energy-momentum tensor of the form [T a b ] = diag(−ρ, p, p, p). For a Klein-Gordon field of mass 1 and the metric (9) we have where the field φ is subject to the Klein-Gordon equation 2 g φ = φ; cf [17, Sec 3.1]. Note that due to spatial homogeneity, the Einstein equations force φ to be independent of the spatial coordinates. Some of the equations are decoupled. For our analysis it suffices to restrict to the reduced coupled part of the LRS Bianchi III Einstein-Klein-Gordon system which is given bẏ with the deceleration parameter H := (ȧ a +2˙b b )/3 denotes the Hubble scalar, ie a measure of the overall isotropic rate of spatial expansion. The corresponding evolution equation (11) is also referred to as the Raychaudhuri equation. HΣ + := (ȧ a −˙b b )/3 denotes the only independent component of the shear tensor, ie a measure of anisotropy in the rate of spatial expansion. Finally, Ω := ρ/(3H 2 ) defines a rescaled energy density which is nonnegative by definition. Since we are interested in non-vacuum solutions we consider Ω to be positive.
The variables are subject to the Hamiltonian constraint Consequently, Σ + and Ω are bounded and take values in the range (−1, 1) and (0, 1) respectively. The Klein-Gordon field φ is unbounded, and so is H a priori. However we will restrict our interest to the case of an (initially) expanding universe, ie to H > 0.
3.2. The system in quasi-standard form. The Klein-Gordon equation (13) has the form of a Van der Pol equation (4) with g(φ,φ) = −3φ and where H(t) plays the role of . The system (11)-(13) therefore describes a harmonic oscillator with nonlinear damping, and where the time dependence of the latter is governed by the coupling of the Einstein equations with the Klein-Gordon equation via H.
Following Subsection 2.2 we apply an amplitude-phase transformation (5) to formulate the Klein-Gordon equation (13) as two first order equations. Furthermore, we transform The result is the system (11), (12) together with the first order Klein-Gordon equationsΩ and where now q is understood to be expressed by (16).
where f 1 (x, t) is given by the square brackets in (12), (17) and (18), and f [2] (x, t) by the square bracket in (11). Note that f 1 , f [2] are independent of H. We see that (19) is resembling the standard form (1) with H(t) playing the role of the perturbation parameter . The crucial difference is however, that H is time dependent and itself subject to an evolution equation which is part of the system. Because of this we say that (19) has quasi-standard form.
3.3. Averaging conjecture. Because (19) has only quasi-standard form and not standard form, we cannot directly apply lemmas 1 or 2 to obtain rigorous error estimates. However we formulate the following conjecture: Conjecture 1. Consider some arbitrary non-vacuum (Ω > 0) initial value satisfying the Hamiltonian constraint (15) and H > 0. Let [H(t), x(t)] T denote the respective solution to (19) and let z(t) denote the solution to the corresponding aver- containing the Σ + and Ω components of the corresponding 3-vectors x(t), z(t). Then there exists a t * such that In Section 4 we provide analytical support for this conjecture and Section 6 also contains numerical support. Under the premise that it holds, we then derive the future asymptotics of Bianchi III Einstein-Klein-Gordon cosmologies in Section 5. The following observations are key: For an initial value with positive H, H is strictly monotonically decreasing with t, and lim t→∞ H(t) = 0.
Proof. The sign of the Raychaudhuri equation (11) is determined by the sign of the factor in square brackets. In Subsection 3.1 we discussed that Σ + ∈ (−1, 1) and Ω ∈ (0, 1). With this we see from (16) that −(1 + q) < 0. From (11) we also see thatḢ = 0 for H = 0. Together this gives the statement of the lemma.
Hence, assuming that Conjecture 1 holds, Proposition 1 states that the error between X(t) and Z(t) goes to zero as t → ∞. Consequently the future asymptotics of the averaged system coincides with that of the full system with respect to Σ + and Ω in this limit. As we show in Section 5, modulo a calculation of the rate of approximation this suffices to derive the asymptotic form of the metric. In Section 6 we find agreement of these results with numerics, and also give convincing numerical support for Conjecture 1. Before all that however, we give strong analytical support for Conjecture 1 in the following section.

Analytical support for conjecture 1
In this section we provide analytical support for Conjecture 1. In Subsection 4.1 we argue that, because of Lemma 3, for large t we can truncate terms of second order in H and work with the resulting first order system in good approximation. We formulate Assumption 1 under which this holds. We then perform a qualitative analysis of the corresponding averaged system and identify its attractors in Subsection 4.2; cf Lemma 4. In Subsection 4.3 we then use this together with lemmas 1 and 2 to give error estimates between the first order and the averaged system. Finally, in Subsecton 4.4 we close our line of arguments heuristically by stating that an infinite series of such first order approximations should yield Conjecture 1 in the continuum limit.
4.1. The first order approximation. As stated in Subsection 3.1 we consider H to be positive initially. Hence, from Lemma 3 we know that H(t) becomes arbitrarily small as t becomes arbitrarily large. This suggests that for sufficiently large t, (19) may be viewed as a perturbative series in H around H = 0 in analogy to (1). Hence for sufficiently large t = t * it seems natural to truncate the second order term to yield the first order approximation 3 where H * denotes the value H(t * ) = H(t * ) which is constant due to the first equation.
This step is heuristic since we did not provide proof that we can truncate in good approximation, and we also did not specify what exactly we mean by the latter. We thus assume the following: Assumption 1. The error between x(t) and y(t) resulting from truncating (19) to Remark 1. Section 6 includes numerical support for the validity of Assumption 1.

4.2.
Qualitative analysis of the averaged system. The second equation of (21) has standard form (1). The associated averaged equation (2) is given by In (21) [H, y] T represents the same quantities as [H, x] T in (19). We rename the vector to indicate the correspondence of the solutions to the respective systems. together withφ = 0. Since the latter equation is decoupled, we focus on the reduced system (22). By the Hamiltonian constraint (15) and the positivity of Ω, the state-space is given by the bounded region Figure 1. The system (22) admits the four equilibrium points T = (−1, 0), Q = (1, 0), D = (1/2, 0) and F = (0, 1). To formulate the following lemma we define the LRS Bianchi I state-space Lemma 4. D is the future attractor for the flow of (22) in closure(X ) \ X I . F is the future attractor for the flow in X I . There is a separatrix from F to D. T (Q) is the past attractor for solutions in closure(X ) with initial data left (right) of the separatrix.
Proof. A check of the eigenvalues of the linearisation of (22) at these points reveals that T, Q and F are hyperbolic. In particular T and Q are local sources while F is a local saddle with its stable manifold tangential to Σ + direction. D is degenerate with a stable manifold tangent to Σ + direction and a centre manifold tangent to the vector [−1/2, 1] T . We prove in Lemma 5 that the flow on that arm of the centre manifold which intersects the state-space is stable. Hence D acts as a local sink for the flow in the closure of the state-space. The information gathered suffices to draw the qualitative flow diagram of Figure 1 from which the statements of the lemma are apparent.

4.3.
Error estimates for the first order approximation. For some initial value, let y(t) denote the solution of the first order approximation (21) and let z(t) denote the solution of the corresponding averaged system given by (22) together withφ = 0. Then from Lemma 1 we know that y(t) − z(t) = O(H * ) on time scales of O(H −1 * ). Furthermore, by Lemma 4 we have a case of averaging with attraction and, in analogy to the discussion of Subsection 2.2, we can extend the validity of this error estimate for all times for the Σ + and Ω components. Hence, denoting the respective 2-vectors by capital letters we have 4.4. Closing the argument. The steps taken so far are as follows: (a) We considered the LRS Bianchi III Einstein-Klein-Gordon system in quasistandard form (19) with solution [H(t), x(t)] T , and for which we showed in We thus argued that for sufficiently large t = t * we should be able to truncate and work with the first order system (21) with solution [H(t), y(t)] T in good approximation. (c) Since the previous argument is heuristic, we made Assumption 1 under which it is valid. That is, we assumed that the error from truncation between x(t) and Hence, we have an O(H * ) truncation error estimate between x(t) and y(t) for t > t * from Assumption 1 in step (c). Furthermore we have an O(H * ) averaging error estimate between Y(t) and Z(t) for t > t * from step (f). In combination we have where X(t) denotes the 2-vector containing the Σ + and Ω components of x(t).
Now recall that H * denotes the value H(t) at the time of truncation t = t * in (c). However because of Lemma 3 in (a), we can choose t * arbitrarily large and H * arbitrarily small. In the continuum limit of infinitely many such truncation times this should yield (20) and thus Conjecture 1. We support this statement with the following line of arguments: Let t * be an arbitrary truncation time in (c) which yields the first order approximation (21) with parameter H * = H(t * ). Consider now the infinite sequence (t k ) ∞ k=0 , with t k := t * + k and > 0. Clearly we have t k+1 > t k for all k. Hence because of Lemma 3 we can consider each of them as separate truncation times, yielding the infinite sequence of first order approximations with parameters (H k = H(t k )) ∞ k=0 . For each we can then make Assumption 1 and perform steps (c)-(f) to yield the infinite sequence of error estimates , which is the statement of Conjecture 1.
Remark 2. For the scope of the present work, we refrain from giving rigorous proofs of Assumption 1 and of the continuum limit above. Because of this, the content of this section does not qualify as proof of Conjecture 1 but rather as strong analytical support of it. Section 6 contains further numerical support. For our analysis on the LRS Bianchi III future asymptotics in Section 5 below we work under the premise that the conjecture holds. Table 1. The exact solutions associated with the equilibrium points of the reduced averaged system (22). a(t), b(t) denote the scale factors of the metric. c 1 , c 2 ∈ R.

The future asymptotics of LRS Bianchi III Einstein-Klein-Gordon cosmologies
In the previous section we provided analytical support for Conjecture 1. Under the premise that it holds, in this section we derive the future asymptotics of LRS Bianchi III Einstein-Klein-Gordon cosmologies. In Subsection 5.1 we give the exact solutions associated with the equilibrium points of the averaged system. In particular we show that D can be associated with the Bianchi III form of flat spacetime. We then interpret the latter in the context of LRS Bianchi III matter solutions in Subsection 5.2, and argue that in order to determine the asymptotic form of the metric we also have to find the rate of approximation of solutions to D. We achieve this by means of a centre manifold analysis in Subsection 5.3; cf Lemma 6. After that, we present our main results in Subsection 5.4 which give the dominant future asymptotic behaviour of LRS Bianchi III Einstein-Klein-Gordon solutions; cf theorems 1 and 2. Finally, in Subsection 5.5 we compare this behaviour to the corresponding future asymptotics in the vacuum case. We point out the quantitative difference and thus speak of the Klein-Gordon future asymptotics to be matter dominated. At the end of this section we also take this opportunity to shed some light on the question when the rate of approximation to the fixed point has to be taken into account and when not. 5.1. Exact solutions associated with the equilibrium points. Each equilibrium point of the averaged system (22) can be associated with an exact solution of the Einstein equations. This is true since for equilibrium solutions the values of Σ + and Ω are by definition constant and hence the evolution equations become particularly simple and can be solved exactly. The results are summarised in Table 1. In the following we lay out the calculation at the example of D.
Starting with the Raychaudhuri equation (11) and evaluating it at D = (1/2, 0) givesḢ From the evolution equations for the spatial metric components (cf for instance [17, (2.29)]) and the variable definitions of H and Σ + given in Subsection 3.1 we havė for the evolution equations for the scale factors of the LRS Bianchi III metric (9). Evaluating these at D yieldsȧ = 0 andḃ = 3bH(t)/2, and using (25) we get where here and in the following c 1 , c 2 ∈ R + . This is the Bianchi III form of flat spacetime; cf [23, p 193, (9.7)].
An analogous calculation yields the scale factors for the solutions associated with T and Q. For T we get a(t) = c 1 t and b(t) = c 2 , which is the Taub Kasner solution; cf [23, Sec 6.2.2 and p 193, (9.6)]. It is spatially flat and of Bianchi type I. For Q we get a(t) = c 1 t −1/3 and b(t) = c 2 t 2/3 , which is the non-flat LRS Kasner Bianchi I solution; cf [23, Sec 6.2.2 and Sec 9.1.1 (2)]. Both are vacuum solutions in the sense that Ω = 0.
From the interpretation of Σ + as the shear parameter given in Subsection 3.1 we can readily infer that F corresponds to an isotropically expanding solution since it has vanishing shear. But we also calculate the asymptotic rate of expansion for large t. The Raychaudhuri equation (11) evaluated at F readṡ Recalling that in the averaged systemφ = 0 (cf Subsection 4.2) we choose ϕ(t) = ϕ = 0 and solve (29) with the initial condition H(0) → ∞ to obtain 4 Evaluating (26)  In Lemma 4we showed that D is the future attractor for solutions in X of the reduced averaged system (22), ie for LRS Bianchi III non-vacuum solutions. Assuming that Conjecture 1 holds this implies that lim t→∞ X(t) = D as well, ie that the Σ + and Ω coordinates of non-vacuum solutions of the original system (19) are asymptotic to the same values.
In Subsection 5.1 we associated with D the Bianchi III form of flat spacetime which is characterised by the scale factors (27). We can however not infer from this, that the scale factors for solutions in X approach (27) in the limit. This is because D fixes the limit form of (26) and thus only the time derivative of the scale factors. Hence the latter will coincide with that of the Bianchi III form of flat spacetime in the limit, but not necessarily the scale factors themselves. We can however obtain the asymptotic form of the scale factors by taking into account the rate of approximation of solutions in X to D.
As stated in the proof of Lemma 4, D is degenerate with a stable manifold tangent to Σ + direction and a centre manifold tangent to the vector [−1/2, 1] T . We can thus get the stability of D for solutions in X as well as the approximation rate for these solutions from a centre manifold analysis (cf Appendix B) as follows.

5.3.
Centre manifold analysis of D. We start by relating (22) to a system of the form (B1)-(B2). Firstly, let us go back for a moment to step (a) of the line of arguments of Subsection 4.4, and hence to the original equation in quasi-standard form (19). We transform to a rescaled time τ according to For convenience we denote the time transformed quantities by the same letters as before, but view them now as functions of τ instead of t. Hence we write the resulting system as and denote its solution by [H(τ ), x(τ )] T . Applying now the analogous steps to (b)-(d) of Subsection 4.4, we obtain the reduced averaged system where again all quantities are functions of τ .  (19) and (31) we also say that the latter is the guiding system of the former.
(34) is of the form (B1)-(B2). Hence, from Lemma B.1 we know that there exists a local centre manifold y = h(x) tangent to the origin. From the eigenvector to the eigenvalue 0 we know that it is tangent to the x-direction. As described at the end of Appendix B, we can approximate the centre manifold to arbitrary degree. For our purposes a lowest non-vanishing order approximation is sufficient. Choosing φ(x) = − 1 6 x 2 and applying the function M defined by (B7), we have The flow on the centre manifold is governed by (B3) which in our case reduces to ∂ τ u = f (u, h(u)). Plugging in (35) and (36) yields and by Lemma B.2 the solution to this equation governs the local stability of the origin of (34) and the local flow around it.
Lemma 5. D is locally asymptotically stable for the flow of (22) in X .
Proof. The zero solution of (37) is attracting solutions with positive u. From Lemma B.2(b) we then get that the zero solution of (34) is locally asymptotically stable for the flow with x ≥ 0. The claim of the lemma then follows from the statement of Remark 4, noting that x ≥ 0 corresponds to Ω ≥ 0; cf (33).
Remark 5. Lemma 5 fills the blank in the proof of Lemma 4. Remark 6. Since there is no attracting structure other than D for the flow in X , the statement of Lemma 5 actually holds globally for that flow, not only locally; cf the proof of Lemma 4.
Remark 7. D is a saddle-node and its local flow is topologically equivalent to the diagram in [14, p 149, Fig 3].

The future asymptotics of LRS Bianchi III non-vacuum solutions.
Let X(τ ) = [Σ + (τ ), Ω(τ )] T denote the respective components of the solution to the guiding system (31) of (19). Let Z(τ ) = [Σ + (τ ), Ω(τ )] T denote the solution to the reduced part of the corresponding averaged guiding system (32). Lemma 6 gives the behaviour of Z(τ ) for large τ . Assuming that Conjecture 1 holds, this gives also an estimate To quantify this estimate we need to determine the dominant behaviour of H for large times.
Lemma 7. Assuming that Conjecture 1 holds, H(t) ≈ 2 3t for large t. Proof. We know from Lemma 6 (or Lemma 4) that lim t→∞ Z(t) = D. Provided that Conjecture 1 holds we then know from Proposition 1 that also lim t→∞ X(t) = D. In Subsection 5.1 we established that the exact solution associated with D is the Bianchi III form of flat spacetime, and we calculated its Hubble scalar in (25). The Hubble scalar of a solution which approaches D will thus approach the Hubble scalar of the Bianchi III form of flat spacetime. Since X approaches D we thus can infer that the corresponding Hubble scalar has the form H(t) = 2 3t (1 + o (1)). The lowest order approximation concludes the proof. Alternatively we can follow from Lemma 6 and Conjecture 1 that for large τ .
From this result we immediately get the asymptotic form of the time transformation (30).
Combining now lemmas 7 and 8 with (39) we have We are now ready to formulate our main results. Theorem 1. Assuming that Conjecture 1 holds, the shear parameter Σ + and the rescaled energy density Ω of LRS Bianchi III Einstein-Klein-Gordon solutions (cf Subsection 3.1) behave like and Ω(t) ≈ 3 2 ln t for large t.
Proof. This is an immediate consequence of (40) since t −1 falls off faster than (ln t) −1 .
Remark 8. The behaviour of Σ + , Ω found in Theorem 1 is identical to what was found in [15, p 1290] in the LRS Bianchi III Einstein-Vlasov case.
What is left to do is to translate the above result to an asymptotic behaviour of the metric in terms of metric time. Proof. Transforming t → τ according to (30) the evolution equations for the scale factors read ∂ τ a(τ ) = a(τ )(1 − 2Σ + (τ )) and ∂ τ b(τ ) = b(τ )(1 + Σ + (τ )). Plugging in the result of Theorem 1 we obtain for large τ , with positive constants c 1 , c 2 , c 2 . In the last approximation for b(τ ) we acknowledged that for large τ the τ dependence is dominated by the exponential such that the square-root can be absorbed into the integration constant in good approximation. Transforming to metric time using Lemma 8 concludes the proof.  (27) and Table 1. Comparing this to our result (41) of Theorem 2 we see that the two solutions agree with respect to the scale factor b(t), which is expanding at a linear rate with t in both cases. However they disagree in the functional form of the scale factor a(t), which is constant in the vacuum case, but forever expanding at a logarithmic rate with t in the Klein-Gordon case. Because of this qualitative difference we say that the Klein-Gordon future is matter dominated.
We take this opportunity to emphasize a circumstance, which albeit having been implicitly taken into account in the literature, to our knowledge has not been pointed out explicitly so far. We discussed in Subsection 5.1 that one can identify the equilibrium points of the reduced dynamical system with a corresponding exact solution to the Einstein equations. Furthermore, in Subsection 5.2 we pointed out that if a solution to the reduced dynamical system converges to an equilibrium point, then this does not necessarily imply that the associated metric converges to the exact solution associated with that point, but that the convergence rate to the point may have to be taken into account as well.
The question now is when the convergence rates are important and when not. We illustrate this at the example at hand. In our case, both the LRS Bianchi III vacuum and Klein-Gordon solutions converge to the equilibrium point D, as can be seen from Figure 1 and our analysis above. D is associated with the Bianchi III form of flat spacetime. However only the vacuum solutions are asymptotic to the latter. For the Klein-Gordon solutions we had to take into account the convergence rates and arrived at the qualitatively different result (41) in Theorem 2.
To illustrate why the convergence rates do not make a difference in the vacuum case, we go over it in the framework of our setup. From Figure 1 and the proof of Lemma 4 we know that generic LRS Bianchi III vacuum solutions lie on the stable manifold of D, and that the eigenvalue associated with that attraction is −3/2. By [14, Sec 2.9, Thm 1] we then have |Σ + − 1/2| ≤ κ exp(−3τ /2) with κ > 0. Transforming again (26) to rescaled time, we have a = a(1 − 2Σ + ). Plugging into this the convergence estimate we get a = ∓2κa exp(−3τ /2), where the sign depends on whether we consider solutions with Σ + ≷ 1/2. Integration yields An analogous calculation can be done for b.
What we have shown is, that even if one takes into account the rate of approximation for the vacuum solutions, this does not alter the result from that associated with the equilibrium point solution. The reason is that the exponential convergence rate along the stable manifold resulted in the exponent in the equation for a(t) to decay as t −1 and thus to vanish for large t.
The general feature appears to be that the exponential convergence rate to hyperbolic sinks, (or on the stable manifolds of hyperbolic saddles or degenerate equilibrium points), results in the convergence rates to not alter the result from the exact solution associated with the fixed point. In cases of degenerate equilibrium points, the generally slower rate of approximation however can make a qualitative difference. In the vast majority of cases in the literature on spatially homogenous cosmology so far, the relevant equilibrium points have been hyperbolic. In these cases the natural association of the equilibrium point solution with the future asymptotic metric works. With degenerate equilibrium points one however has to be more careful as shows the case of [15] as well as the case of the present work.

Numerical support for the results
In this section we test the agreement of our analytical results and assumptions against numerics. We give the numerical setup and equations in Subsection 6.1 and elaborate on the chosen time variable and simulation range. In subsections 6.2-6.4 we then present the simulation outcome in the form of flow plots ( Figure 3) and function plots in normal ( Figure 4) and logarithmic scales ( Figure 5). The plots convincingly demonstrate agreement with our analytical results and assumptions. Finally, in Subsection 6.5 we briefly investigate and comment on the possible past asymptotics at the hand of Figure 6. 6.1. Simulation setup. We simulate two systems: Firstly we solve the guiding system (31) of (19) and simultaneously also integrate the differential form of (30) to capture the transformation between t and τ . Hence we can recover the solution of the original system (19), ie the quantities as functions  Table 2. The six initial data sets for our simulations of the systems (42)-(44) and (45)-(47). We essentially vary (Σ + (0), Ω(0)), ie the location of the initial data point projected into X , except for iv where we also vary H(0).
of t. For a better overview we write out the full system explicitly: with q given by (16). Note that as we work in the guiding system all quantities including t are seen as functions of τ .
As initial data we pick the six sets given in Table 2 and refer to the corresponding solutions by the Roman numerals in the first column. These data satisfy the Hamiltonian constraint (15) and thus, projected into the Σ + , Ω plane, lie in X . We integrate the averaged system (45)-(47) in the interval τ ∈ [−40, 40] and plot the corresponding solutions in orange in the figures of this section. The solutions to the full system (42)-(44) are plotted in blue, and since we are primarily interested in the future asymptotics, we integrate these from τ = 0 to τ ≈ 10 for all figures except Figure 6 where we also integrate them to τ = −40. The cutoff at τ ≈ 10 has the the following reason: As will become clear from the discussion below, the solutions to the full system oscillate with exponentially increasing frequency; cf also Lemma 8. 6 Hence the numerical grid has to be refined to smaller and smaller step sizes. The cutoff occurs before roundoff errors become an issue, which is at τ ≈ 10 for all solutions. Despite this limitation, our numerical results give a clear demonstration of the validity of our analytical results.  Figure 2 shows a flow plot of the solutions in the H, Σ + , Ω domain. The solutions to the full system oscillate around the solutions to the averaged system with a decaying envelope, and the plot demonstrates how all solutions converge to the point D located at (Σ + , Ω) = (1/2, 0) in the H = 0 plane. One can also see that, figuratively speaking, this convergence occurs in three steps as follows: Firstly the solutions approach the H = 0 plane, then the centre manifold of D and finally they approach D itself along its centre manifold. On the one hand this demonstrates that the truncation from the full system to the first order approximation performed in Subsection 4.1, as a consequence of which H was approximated by a constant, was legitimate, and hence directly supports Assumption 1. On the other hand, it directly supports Lemma 5. Figure 3 shows the projections of the trajectories into H = 0.  Figure 3b) and thus support Lemma 5. That this approach is slow in comparison to the speeds of the trajectories away from the local centre manifold can be seen as follows: While the endpoints of the trajectories to the full solutions are reached at τ ≈ 10 the trajectories of the averaged solutions do not advance much further until they terminate at τ = 40, and this advance even becomes smaller and smaller the closer the initial data is to D.
6.3. Function plots. Figure 4 restricts to solution i and shows H, Σ + , Ω and t plotted against τ together with the respective averaged quantities, as well as separate error plots between the full and averaged solutions. The function plots of figures 4a, 4b and 4d demonstrate nicely the oscillations of the full quantities on top of the averaged ones. The corresponding error plots demonstrate the fast convergence of the full to the averaged quantities. In particular, we see that the decay of the errors of Σ + , Ω indeed follow an envelope proportional to H(τ ) as conjectured in Conjecture 1. 7 The error of H on the other hand appears to decay with H(τ ) 2 , which we would await intuitively from the fact that the Raychaudhuri equation (left equation of (42)) is of second order in H(τ ) and vanishes to first order.
Comparing Figure 4a to figures 4b and 4d we have again confirmation of the faster decay of H to 0 in comparison to the convergence rate of (Σ + , Ω) to D. In other words, we have again confirmation that the truncation to the first order approximation is legitimate; cf Subsection 4.1. 7 There is a small offset of the mean errors of Σ + , Ω from zero which converges away slower than O(H(τ )). In our numerical experiments we found that the magnitude of this error depends on the chosen initial phase ϕ(0). We leave a deeper analysis of this phenomenon as an open problem. The function plot in Figure 4c shows the exponential growth of t with τ , supporting Lemma 8. It is apparent from this plot that in a simulation in t instead of τ one would have to evolve to very large times in order to enter the asymptotic regime, which we already pointed out in Subsection 6.1. This is demonstrated even better in the following subsection.
6.4. Future asymptotics. In Figure 5 we show plots of the same quantities as in Figure 4, however in logarithmic (figures 5a, 5c) and double logarithmic scalings (figures 5b, 5d), together with the analytically determined asymptotic behaviours (dashed). Due to the logarithmic scaling, the asymptotic behaviours are determined by straight lines. This allows for a direct comparison of the numerics to the statements of Theorem 1, and for an indirect comparison to the statements of Theorem 2. Note also the higher plot range up to τ = 40. The full solutions are however still plotted up to τ ≈ 10, where we terminated the integration. However this range is sufficient, especially in combination with the continuation of the averaged solution, to convincingly demonstrate the asymptotics.
From the plots it is apparent how all quantities indeed approach the respective analytically determined asymptotic behaviours for large τ . It is also apparent that this happens faster for H (and thus also for t) than for Σ + and Ω, which we already stressed in the preceding subsections.
From (c) we stress again the very large values of t required to enter the asymptotic regime. In fact, of all our solutions, solution i plotted here turned out to be the one which reaches the asymptotic regime the fastest, which is why we chose it for figures 4 and 5. 6.5. Past asymptotics. Though our focus here lies on the future asymptotics, it is also interesting to investigate what happens if we integrate our initial data backwards in time. In Figure 6 we thus show again the flow plot of Figure 3(a), however this time the full solutions are also integrated to τ = −40 into the past. As we would expect since H(τ ) is strictly monotonically decreasing, the oscillations around the averaged solutions become larger and larger into the past, and quickly reach the dimensions of the state-space. When this happens, each of our solutions approaches the Bianchi I boundary. Hence one may conjecture that asymptotically into the past LRS Bianchi III Einstein-Klein-Gordon solutions will behave like LRS Bianchi I solutions, and potentially will converge to the LRS Kasner solutions associated with T, Q; cf Subsection 5.2. However we have to be cautious with this intuition, since as stressed already beforehand, X does not represent the full reduced state-space. A deeper investigation is necessary in order to make a concrete statement about the past asymptotics. We leave this as an open question.

Discussion
We conclude with a brief summary of our main results (Subsection 7.1), some remarks on the novelty and utility of our averaging approach (Subsection 7.2) and an outlook on open problems (Subsection 7.3). 7.1. Summary of the main results. Our main results give the behaviour of LRS Bianchi III Einstein-Klein-Gordon solutions for large times in terms of the shear variable and energy density (Theorem 1) and in terms of the metric components (Theorem 2). The theorems have Conjecture 1 as premise, for which we gave strong numerical and analytical support. The numerics also convincingly demonstrates agreement with the conclusions of the theorems.
Our results are global in the sense that all solutions show this behaviour future asymptotically. The asymptotic metric we found is forever expanding in all spatial directions. The scale factor of the 2-hyperboloidal part of the spatial geometry expands ∝ t while the radial part only expands ∝ ln t. This result matches with that of [15] where massive Vlasov matter was considered.
Due to the anisotropic expansion the shear variable is non-vanishing in the limit for time to infinity. The energy density on the other hand converges to vacuum. Despite that, we speak of the future asymptotics to be matter dominated, since it is qualitatively different from the future asymptotics in the vacuum case, where the first scale factor goes to a constant; cf Subsection 5.5. 7.2. Utility of the averaging approach. Apart from our main results, another aspect of our work which is novel in the context of mathematical cosmology is the approach we took. The application of averaging methods allowed us to control the oscillations stemming from the Klein-Gordon equation and at the same time resulted in a dynamical system with decoupled Raychaudhuri equation such that then standard methods of dynamical systems in mathematical cosmology could be applied.
Due to the decay of the perturbation parameter to zero, ie of the Hubble scalar, and since the averaged system is attracted by an equilibrium point, the oscillations in the full solution decay and vanish in the limit for time to infinity. Consequently, the full and the averaged solutions converge in that limit. Our results for the limit thus concern the solutions to the full system and are not a mere averaging approximation.
7.3. Outlook. We expect these techniques to also be applicable to other classes of spatially homogenous Einstein-Klein-Gordon cosmologies. Furthermore they also should be a useful tool in cases where oscillations or other perturbations play a role. For instance, we found striking similarities between the system discussed here, and the evolution equations encountered in [9,13,24]. These papers discuss the future asymptotics of different classes of spatially homogenous perfect fluid cosmologies, and encountered oscillatory behaviour and a phenomenon called selfsimilarity breaking; cf [22]. Averaging methods might be suited to shed new light on this phenomenon.
Though our results are robust given the strong agreement with numerics and the strong analytical support, it is desirable to rigorously prove Conjecture 1. Section 4 suggests one route to do so, which isvia providing rigorous proofs to Assumption 1 and the limiting process outlined at the end of that section. Finally, we point out the proof of a statement closely related to that of Conjecture 1 in the closely related work [7]. (If the vector field f (x, t) contains parameters we assume that the parameters and the initial conditions are independent of , and that the limit is uniform in the parameters.)

Appendix B. Centre manifold analysis
Here we closely follow the discussion in [4,Chap 1]. Given an autonomous dynamical system ∂ τ x = f (x), x ∈ R n a set S ⊂ R n is a local invariant manifold for the system if for initial data x(0) = x 0 ∈ S the corresponding solution x(τ ) is in S for |τ | < T where T > 0. If we can always choose T = ∞, then we say that S is an invariant manifold. Consider now the system ∂ τ x = Ax + f (x, y), x ∈ R n , f ∈ C 2 (B1) ∂ τ y = By + g(x, y), y ∈ R m , g ∈ C 2 (B2) with constant matrices A, B such that all eigenvalues of A have zero real parts while all eigenvalues of B have negative real parts. Further f (0, 0) = 0, f (0, 0) = 0, g(0, 0) = 0, g (0, 0) = 0, where f , g denote the Jaccobians. Note that the system exhibits an equilibrium point at the origin.
If y = h(x) is a (local) invariant manifold of (B1)-(B2) and h is smooth, then it is called a (local) centre manifold if h(0) = 0, h (0) = 0. Where it is clear from the context, we also speak of a centre manifold and really mean a local centre manifold.
In other words, solutions of (B1)-(B2) with initial value sufficiently close to the origin, for large times approach solutions of (B3), ie on the centre manifold, at an asymptotic rate.
To solve (B3) we need to know the centre manifold y = h(x) or at least an approximation of it for small |x|. Substituting y = h(x(τ )) into (B2) yields h (x) Ax + f (x, h(x)) = Bh(x) + g(x, h(x)).

(B6)
For the centre manifold we have h(0) = 0, h (0) = 0. Consequently, to obtain the centre manifold one has to solve (B6) to these conditions. In general this is impossible analytically. However, with the following theorem we can always approximate the centre manifold up to an arbitrarily high accuracy. Before we formulate the theorem we define the following map M on functions φ : R n → R m which are C 1 in a neighbourhood of the origin: