Cosmological dynamics of viable f(R) theories of gravity
CCosmological dynamics of viable f ( R ) theories of gravity Sulona Kandhai ∗ and Peter K. S. Dunsby †
1, 2 Astrophysics, Cosmology and Gravity Centre (ACGC),Department of Mathematics and Applied Mathematics,University of Cape Town, Rondebosch 7701, Cape Town, South Africa. South African Astronomical Observatory, Observatory 7925, Cape Town, South Africa
A complete analysis of the dynamics of the Hu-Sawicki modification to General Relativity ispresented. In particular, the full phase-space is given for the case in which the model parametersare taken to be n = 1 , c = 1 and several stable de Sitter equilibrium points together with anunstable “matter-like” point are identified.We find that if the cosmological parameters are chosen to take on their ΛCDM values today, thisresults in a universe which, until very low redshifts, is dominated by an equation of state parameterequal to , leading to an expansion history very different from ΛCDM. We demonstrate that thisproblem can be resolved by choosing ΛCDM initial conditions at high redshift and integrating theequations to the present day. PACS numbers: 04.50.Kd, 98.80.-k, 98.80.Cq, 12.60.-i
I. INTRODUCTION
Over the past decade, our ability to make high preci-sion measurements of the cosmic microwave backgroundand distance measurements of Type Ia supernovae haveled to the remarkable conclusion that the expansion rateof the universe is accelerating.The most widely accepted cosmological model of theuniverse, based on General Relativity –
The Concor-dance Model – attributes this acceleration to a myste-rious, exotic form of energy density called dark energy ,which dominates the energy density of the universe at thepresent time. While this model is by far the simplest andmost successful, in terms of fitting the current cosmolog-ical data, there are a number of shortcomings associatedwith this theory, such as fine tuning issues, the coinci-dence problem and the need for unknown, exotic formsof energy.Naturally, such fundamental problems begs the ques-tion of whether or not our understanding of physics iscomplete on the scales which we are applying it; for ex-ample, could the late time acceleration of the universe bea manifestation of the break down of the gravitational in-teraction on cosmological scales? It is for this reason thatviable alternative theories of gravitation are currently re-ceiving a great deal of attention.A very popular class of alternatives to the ΛCDMmodel involve modifying, the Einstein-Hilbert action byconsidering terms containing higher order curvature in-variants. This results in a gravitational theory whichproduces predictions that differ from General Relativity only late in the matter dominated era of the universe,producing an effective dark energy term which leads tothe observed late time acceleration, as well as retaining ∗ [email protected] † peter.dunsby [at] uct.ac.za the success of General Relativty on solar system scales.These f ( R ) theories of gravity, by their construction, in-troduce an additional scalar degree of freedom which hasa coupling to matter that is extremely weak [6],[31],[7].This scalar degree of freedom results in a long range fifthforce which predicts that the curvature of the space timein the neighbourhood of a localised energy density is un-coupled from the energy density. Taking this effect intoconsideration, several conditions must be satisfied in or-der for a f ( R ) theory to be compatible in both the lowand high curvature regimes [31]. It is important that theexpansion history in the early universe matches what isobtained from General Relativity, in order to be consis-tent with big bang nucleosynthesis, as well as the factthat the f ( R ) theory should generate a late time acceler-ation compatible with current observations, without theintroduction of a cosmological constant.If we take the function f ( R ) to be made up of the usuallinear Lagrangian plus some corrective terms encapsu-lated in a second function g ( R ): f ( R ) = R + g ( R ), thenthe above conditions can be expressed as the followinglimits on g ( R ): lim R →∞ g ( R ) = Constant , lim R → g ( R ) = 0.The first limit corresponds to epochs during which themodel f ( R ) should resemble a late time ΛCDM universe,containing an effective cosmological constant, for whichthe Lagrangian looks like R − f ( R ) model chosen mustbehave close to standard GR.So far, the most compelling models proposed, whichsatisfy these limits are broken power law functions, due tothe fact that such functions can be parametrised to givethe required behaviour in the relevant regimes. Severalmodels have been proposed which are designed to satisfythe above limits [31, 61] and we will focus on the first ofthese, introduced originally by Hu & Sawicki [31].The main difficulty which plagues the investigation ofsuch models results from the inherent complexity associ-ated with the inclusion of higher order invariants of thecurvature tensor into the gravitational action. It is con- a r X i v : . [ g r- q c ] O c t sequently extremely difficult if not impossible to derivethe cosmological dynamics and obtain exact solutions forsuch theories.One method, which has been found to be particularlyeffective in overcoming many of the issues encounteredwith higher order theories, is the so-called dynamical sys-tems approach to cosmology [2, 10, 12]. This techniquerelies on the ability to express the cosmological equa-tions as a set of autonomous differential equations, usinga set of suitably chosen variables. In this paper we usethese techniques to proved a detailed description of thecosmological dynamics of the Hu & Sawicki model anddetermine under what conditions a ΛCDM-like expansionhistory can be obtained.In Section II below we briefly present the field equa-tions for f ( R ) gravity in a metric formalism. SectionIII discusses the Hu-Sawicki model, which is the focusof analyses performed below. We present the dynamicalsystems approach to f ( R ) gravity, give the general formof the dynamical system for a general f ( R ) and carefullydetail the analysis and its results in Section IV, wherewe discuss the fixed points, and their stability, both ina compact and non-compact phase space. The expan-sion history for the particular Hu-Sawicki model, consid-ered in this paper is calculated in Section V to providea comparison with the Concordance model, as well as toillustrate the effects of fixing the initial conditions on theviability of the model. Finally, we review the results pre-sented and discuss conclusions and various limitations ofthe analysis in Section VI. Throughout this paper naturalunits are assumed ( (cid:126) = c = k B = 8 πG = 1). II. f ( R ) COSMOLOGY
In principle there is no a-priory reason to restrict thegravitational Lagrangian. In fact, it is quite possible thatthe addition of higher powers of R and correspondinginvariants may improve the characterisation of gravita-tional fields near regions where R → ∞ [4]. Furtherimpetus to explore results of nonlinear gravitational La-grangians is due to the fact that every theory attemptingto unify the fundamental interactions require either thatthere are non minimal couplings to the geometry or thathigher order curvature invariants appear in the action[6],[4].In higher order gravity theories, the Einstein-Hilbertaction is modified by considering the addition of higherorder, most commonly, second order , curvature invari-ants. Second order modifications result in quadratic La-grangians, which involve some of the four possible cur-vature invariants of the Ricci scalar, Ricci and Riemanntensors; R , R ab R ab , R abcd R abcd and ε iklm R ikst R stlm , here ε iklm is the completely antisymmetric tensor, of rankfour. In fact, by making use of the following identities, which are true for all 4-dimensional space times:( δ/δg ab ) (cid:90) d x √− g ( R abcd R abcd − R ab R ab + R ) = 0 , (1)( δ/δg ab ) (cid:90) d x √− gε abcd R abcd R efcd = 0 , (2)as well as the fact that the following identity is true forspacetimes with maximally symmetric spatial sections,( δ/δg ab ) (cid:90) d x √− g (3 R ab R ab − R ) = 0 , (3)it is possible to deduce that any quadratic modificationsto the Lagrangian can in general be represented by com-binations of powers of the Ricci scalar. The modifiedLagrangian is then given by: L = √− g ( f ( R ) + L m ) , (4)where L m is the standard matter Lagrangian. ClearlyGeneral Relativity is recovered when f ( R ) = R − G αβ + f (cid:48) R αβ − (cid:0) f − (cid:3) f (cid:48) (cid:1) g αβ − ∇ α ∇ β f (cid:48) = T αβ . (5)Here primes denote derivatives with respect to the Ricciscalar. We are only concerned with modifications toGR at late times, so T αβ is described completely bythe matter dominated stress-energy tensor. For thespatially flat ( k = 0) FLRW metric with signature( − + ++), where the matter is described as a perfectfluid with equation of state w = pρ , the independentfield equations in terms of the function f ( R ) are (dotsrepresent temporal derivatives): The Raychaudhury equation H + 3 H = − f (cid:48) (cid:16) p m + 2 H ˙ f (cid:48) + ( f − Rf (cid:48) )+ ˙ R f (cid:48)(cid:48)(cid:48) + ¨ Rf (cid:48)(cid:48) (cid:17) , (6)where H represents the usual expansion rate, H = ˙ aa , the Friedmann equation , H = f (cid:48) (cid:16) ρ m + ( Rf (cid:48) − f ) − H ˙ f (cid:48) (cid:17) . (7) The trace equation Rf (cid:48)(cid:48) = ρ (1 − w ) + f (cid:48) R − f − Hf (cid:48)(cid:48) ˙ R − f (cid:48)(cid:48)(cid:48) ˙ R , (8)In addition the Energy Conservation equation for stan-dard matter is given by˙ ρ = − H (1 + w ) ρ. (9)Finally an expression for the Ricci scalar in terms of H can be obtained by combining the Raychaudhuri and theFriedmann equations: R = 6 ˙ H + 12 H . (10) III. THE HU-SAWICKI MODEL
The Hu-Sawicki model (HS model) for f ( R )-gravity[31] was designed specifically to overcome several prob-lems faced by other f ( R ) functions considered as alter-natives to the Lagrangian of General Relativity. It isconstructed in such a way that it satisfies the limitsmentioned in the previous section, and therefore, givenan appropriate choice of parameters, is able to recoverstandard General Relativity at high redshifts, as well asmimic ΛCDM at low redshifts. See [31] for details onthis class of models. Below is a short summary of theparameters which enter this theory. The functional formof f ( R ) is characterised by the correction g ( R ), which isa general class of broken power laws: f ( R ) = R − m c (cid:0) Rm (cid:1) n c (cid:0) Rm (cid:1) n + 1 . (11)Here n > m is taken to be: m ≡ κ ¯ ρ M pc ) − (cid:18) Ω m h . (cid:19) , (12)where ¯ ρ is the average density of the present epoch and κ = 8 πG . c and c are dimensionless parameters,whose relationship will be shown to be associated withthe present matter density parameter, Ω m .The sign of the second derivative of the correction, g ( R ), is required to be strictly positive for high curva-tures relative to density in order to guarantee that athigh density, the solution is stable in this regime. g RR ≡ d g ( R ) dR > R (cid:29) m . (13)The above condition is important to ensure that at highredshifts, the General Relativity values of quantities,such as Ω m h are recovered. For convenience we define g R ≡ dg/dR .The key appeal of this class of models is that, explicitly,there is no cosmological constant term. However, if weconsider the limit of very high curvature relative to m :lim m /R → f ( R ) ≈ R − c c m + c c m (cid:18) m R (cid:19) n , (14)it can be seen that, for a fixed value of c c , the limitingcase of c c → c c , the curvature is frozen to a fixed value and no longerdecreases with matter density; this gives rise to a set of models which exhibit late time acceleration similar toΛCDM [31].It follows that [31] admits the following relationshipbetween the parameters and the energy densities of thecosmological constant, ˜Ω Λ , and matter, ˜Ω m : c c ≈ Λ ˜Ω m . (15)The remaining two parameters n and c c are free to choosein order to determine the expansion history and howwell it fits observations, compared to the ΛCDM model.Larger values of n allows the model to mimic ΛCDM un-til later in the expansion history, while smaller values of c /c leads to a reduction in general deviations from it.Because both the Hubble parameter and the critical den-sity depend on the correction g R , ˜Ω m only becomes theactual value of the matter density in the following limit:lim c /c → ˜Ω m = Ω m . (16)The matter density in the physical units, however, re-mains unchanged.
1. Fixing the parameter values today
Considering the ΛCDM model with flat spatial geom-etry, we have: R ≈ m (cid:32) a + 4 ˜Ω Λ ˜Ω m (cid:33) , (17)and for the derivative of the correction: g R = − n c c (cid:18) m R (cid:19) n +1 . (18)The values for today are: R ≈ m (cid:18) m − (cid:19) , g R ≈ − n c c (cid:18) m − (cid:19) − n − (19)for | g R | (cid:28)
1. [31] uses the following specific values:˜Ω m = 0 . , ˜Ω Λ = 0 . . (20)Therefore we have: R = 41 m , g R ≈ − n c c (cid:18) (cid:19) n +1 . (21)Figure 9 in [31] show the range and combinations of pa-rameter values, g R and n , which are acceptable and en-able the model to satisfy solar system tests. IV. DYNAMICAL SYSTEMS APPROACH TO f ( R ) GRAVITY
We follow the general strategy for applications of thedynamical systems approach to the cosmology of fourthorder gravity, developed in [1, 2, 10, 12]. This provides away of analysing the dynamics of any analytic function f ( R ) which is invertible for the Ricci scalar. Followingthe approach of [2, 26, 27], a compact analysis of thephase space for a general f ( R ) theory is performed, bydefining a strictly positive normalisation to pull any so-lutions at infinity into a finite volume. In order to makethis possible, an appropriately normalised time variableneeds to be defined. Ensuring that this time variableis strictly non-decreasing amounts to guaranteeing thatthe expansion-normalisation adopted is strictly positive.Thus, the sign of time is maintained, whether study-ing expanding or collapsing cosmologies. This meansthat any negative contribution to the Friedmann equa-tion must be absorbed into the normalisation. In the casewhere a quantity can be both positive and negative, eachoption must be studied in a separate sector of the phasespace. The full phase space can then be reconstructedby simply aligning the various sectors along their com-mon boundaries. A phase space which is constructed inthis way can therefore include all static, bouncing andre-collapsing models.In this paper we will only consider sectors of the phasespace where R ≥
0, simply because negative Ricci scalarvalues are not of any real physical interest. We require f (cid:48) , f (cid:48)(cid:48) >
0, and the matter density is assumed to be non-negative. The function f under consideration here, isalways positive. Therefore, in (7), the term which con-tains f attached to a negative sign must be absorbed intothe positive-definite normalisation. A. Compact phase space analysis of f ( R ) gravity Rewriting the modified Friedmann equation (7) in thefollowing way allows a convenient definition of the nor-malised variables: (cid:32) H + 32 ˙ f (cid:48) f (cid:48) (cid:33) + 32 ff (cid:48) = 3 ρ m f (cid:48) + 32 R + (cid:32)
32 ˙ ff (cid:48) (cid:33) . (22)The left hand side of the above equation is a positivedefinite quantity. Quite naturally, assigning each term in(22) a name, we obtain the following dynamical variables: x = 32 ˙ f (cid:48) f (cid:48) D , v = 32 RD , (23) y = 32 ff (cid:48) D , Ω = 3 ρ m f (cid:48) D , Q = 3 HD ,
Figure 1:
Top panel - this gives a perspective of the 3-dimensional compact phase space of the Hu-Sawicki modelfor n = 1 , c = 1. This plot demonstrates the anti-symmetrybetween the expanding and contracting regions of the phasespace, defined by the Q = 0 plane, as well as the possibil-ity that orbits can cross this plane – implying the existenceof universe models with bounce behaviour. Bottom panel -This gives the phase plane of the invariant submanifold v = 0,corresponding to universe models with zero Ricci curvature. where D represents the normalisation which compactifiesthe phase space, and takes the form D = (cid:32) H + 32 ˙ f (cid:48) f (cid:48) (cid:33) + 32 ff (cid:48) . (24)The normalised time variable is then defined as follows: ddτ ≡ D ddt . (25)Thus, (22) and (24) establish two independent constraintequations for our system:1 = Ω + x + v , (26)1 = ( Q + x ) + y . (27)The dynamical variables defined by (23) constitute thecoordinates of our compact phase space and the bound-aries of this phase space are defined by the above twoconstraints as follows: − ≤ x ≤ , ≤ Ω ≤ , − ≤ Q ≤ , (28)0 ≤ v ≤ , ≤ y ≤ . The construction of D ensures that the dynamical vari-ables are well defined when H = 0, thus we expect allstatic, expanding, collapsing and bounce solutions to beincluded. Expanding and collapsing universes are con-nected across the Q = 0 boundary. B. The General Propagation Equations
Differentiating the dynamical variables (23) with re-spect to τ and substituting the independent cosmolog-ical equations, (7) − (8), leads to a set of 5 first orderautonomous differential equations. The dimensionalityof the system can be reduced by using the constraintequations (26) and (27) to eliminate y and Ω. Below, weshow the general propagation equations for the reduced3-dimensional autonomous system, whereΓ ≡ f (cid:48) f (cid:48)(cid:48) R (29)specifies the specific f ( R ) model. In order to close thesystem, Γ must be expressed in terms of the dynami-cal variables. This implies that the above system char-acterises a general dynamical system for any modifiedgravity cosmology defined by a function f ( R ), which isinvertible in terms of the dynamical variables, in otherwords, Γ needs to be expressed as a function of ( x, Q, v ): dvdτ = − v (cid:2) ( Q + x ) (cid:0) v + 4 xQ − (cid:0) − v − x (cid:1) (1 + 3 w ) (cid:1) − Q − x + 2 x Γ ( v − (cid:3) , (30) dxdτ = 16 (cid:2) − x v Γ + (cid:0) − v − x (cid:1) (1 − w ) + 2 v + 4 (cid:0) x − (cid:1) (cid:0) − Q − xQ (cid:1) + x ( Q + x ) (cid:0)(cid:0) − v − x (cid:1) (1 + 3 w ) − v (cid:1)(cid:3) , (31) dQdτ = 16 (cid:2) − xQ + (5 + 3 w ) Qx (1 − xQ ) − Q (1 − w ) − Qx (1 + 3 w ) − vQ (1 + w ) ( Q + x ) + 2 v (1 − Γ Qx )] . (32)Clearly, v = 0 corresponds to an invariant submanifold;solutions which start their evolution in this plane willremain there forever. This submanifold corresponds touniverse models for which the Ricci scalar vanishes. Dueto the existence of an invariant submanifold, no globalattractor can exist for cosmological systems defined bythe above compact dynamical variables.It is useful to express the cosmological equations termsof the dynamical variables. By first expressing the sec-ond time derivative of f (cid:48) , given by the trace equation, interms of these variables:¨ f (cid:48) f (cid:48) = H Q [(1 − w )Ω + 2 v − y − xQ ] . (33)and then using (33) and the constraint equations, we canwrite the modified Raychaudhuri equation (6) in the fol- lowing useful form:˙ H = − H Q (cid:0) − y − x (cid:1) . (34)In this way, once the fixed points have been determined,the expansion history at these points can be found di-rectly by integrating this equation. C. Dynamical Systems analysis of the Hu-Sawickimodel
Let us now focus on the HS model as described in [31],where f ( R ) is given by: f ( R ) = R − m c (cid:0) Rm (cid:1) n c (cid:0) Rm (cid:1) n + 1 . (35)We define the following relations: m ≡ λH , r ≡ RH , h ≡ HH , (36)which enables us to write (35) in a more convenient form: f ( R ) = CH (cid:34)(cid:16) rC (cid:17) − c (cid:0) rC (cid:1) n c (cid:0) rC (cid:1) n + 1 (cid:35) . (37)Here c , c and n are constant parameters to be con-strained by observations and the dimensionless param-eter λ is related to the ratio of m defined by (12) andthe present day value of the Hubble parameter.We perform the compact dynamical systems analysisoutlined in the previous section is to (35) or (37) for thecase n = 1. It has been shown that n remains uncon-strained by current cosmological data [53] and, thereforedoes not compromise generality at this stage. For sim-plicity we also set c = 1. This has the added benefit ofleading to a Γ which can be expressed entirely in termsof the dynamical systems variables:Γ ≡ f (cid:48) f (cid:48)(cid:48) R = 12 vy ( v − y ) . (38)The expression for Γ associated with this specific model(38) is then substituted into equations (30) - (32), anda dynamical systems analysis is performed to obtain theequilibrium points of the system. The equilibrium points,as well as their corresponding exact solutions for the scalefactor are given in Table I. The Hartman-Grobman the-orem is used to determine the stability of these fixedpoints, where possible. Some points obtained are non-hyperbolic, and in this case we resort to the Center Man-ifold Theorem or other techniques to clarify their nature.In Figure 1, we present a view of the compact phasespace by showing several trajectories between the fixedpoints; in the upper panel we show the entire 3D phasespace, and in the lower panel we show the invariant submanifold corresponding to v = 0. D. Stationary points, stability and exact solutions
For the HS model, with n = c = 1 and keeping anarbitrary equation of state w , the fixed points for theentire phase space, as well as the exact solution of thescale factor at each point, are summarised in Table I.Note that although Table II includes a stability classi-fication of the fixed points identified for universes domi-nated respectively by radiation, matter and a cosmolog-ical constant, the bulk of the analysis is focused on uni-verses dominated by a dust fluid ( w = 0). Accordingly,the phase-portraits below are all generated for w = 0.The points for which the v and/or y coordinates are ex-actly zero make up a subset of fixed points which appearin the dynamical systems analysis for any arbitrary f ( R )theory of gravity since Γ = 0 and the resulting equations are independent of f ( R ). For example, we find the cor-responding fixed points for which v = 0 that appear in[2] for f ( R ) = R + αR n .
1. Stability
The stability of all the fixed points, save A ± , K ± and J , could be inferred using the Hartman-Grobman the-orem. The fixed points A ± are non-hyperbolic, so theCenter Manifold Theorem was used to determine the sta-bility of A +1 Perturbation theory was employed to findthe stability near A − .
2. Non analytic points K + and J Considering the structure of Γ, as expressed in (38)in terms of the variables, it is clear that when v = y ,this term becomes undefined. Thus, analytically, it isnot possible to identify and analyse stationary points inthe phase-space when this occurs. The matter-like points K ± , which lie on the intersection of the y = v and x = 0planes as detailed in [2] are also present here. The point J which lies on the v = y = 1 surface corresponds to asaddle point and represents a static universe.The respective stabilities of these points were inferredby inspection; by examining the time evolution of thecoordinates, when initial conditions are chosen near thefixed points. E. Exact solutions
The rate of expansion, H , and the deceleration param-eter, q , are related by the following expression:˙ H = − (1 + q ) H , (39)where q = − ¨ aa ˙ a . This relationship can be used to deter-mine the time evolution of the scale factor at each of theequilibrium points, provided the deceleration parameterat that equilibrium point is known. In order to find thevalue of q at the i th equilibrium point, we need to express q in terms of the compact variables x, y, v, Q, Ω .It follows, from the Raychaudhuri equation (34) andthe constraint equations, that q i = 1 − z i Q i . (40)For the cases in which Q (cid:54) = 0, direct integration of(39) results in an expression describing the evolution of The system is significantly less complicated using the finite vari-ables defined in Section IV F given by (44), therefore the flowof the Center Manifold at A + is analysed using this coordinatesystem. Point Coordinate (
Q, x, y, v,
Ω) Scale factor evolution, a ( t ) A ± (cid:104) ± √ , , , , (cid:105) a ( t ) = a e H ( t − t ) B ± (cid:2) ± , ± , , , (cid:3) a ( t ) = a e H ( t − t ) C ( w = − (cid:20) , (cid:113) w +1)1+3 w , − w , , − w (cid:21) a ( t ) = a D ( w = − (cid:20) , − (cid:113) w +1)1+3 w , − w , , − w (cid:21) a ( t ) = a E [0 , , , , a ( t ) = a F [0 , − , , , a ( t ) = a G ± [ ± , ∓ , , , a ( t ) = a (2 H ( t − t ) + 1) H ± ( w ≤ ) (cid:104) ∓ w − , ± w − w − , , , −
49 3 w − w − (cid:105) a ( t ) = a (2 H ( t − t ) + 1) J [0 , , , , a ( t ) = a K ± (cid:104) ± √ , , , , (cid:105) a ( t ) = a (cid:0) H ( t − t ) + 1 (cid:1) Table I: Each equilibrium point in the phase space has an expanding and collapsing version, denoted by a subscript +/-. Fixedpoints with Q = 0, correspond to static universes. The points H ± , C and D depend on the equation of state parameter, w ; H ± only lies in the phase space for − ≤ w ≤ , and while C and D lie in the phase space for all values of w ≤ −
1. For ourpurposes − ≤ w ≤
1, therefore we only consider the case when w = − C , D , E , F , G ± , H ± all lieon the invariant submanifold v = 0. Other w dependent fixed points were found, but are not included in this analysis as theyonly lie within the phase space for unphysical values of the equation of state parameter.Point Eigenvalues of Jacobian Stability w = 0 w = w = − A + (cid:104) , − √ , − √ − √ w (cid:105) Attractor Attractor Attractor A − (cid:104) , √ , √ + √ w (cid:105) Repellor Repellor Repellor B + (cid:2) − , − , − − w (cid:3) Attractor Attractor Attractor B − (cid:2) , , + w (cid:3) Repellor Repellor Repellor C (cid:20) (cid:113) w )1+3 w , − (cid:113) w )1+3 w , (cid:113) w )1+3 w (cid:21) (cid:3) (cid:3) Saddle D (cid:20) (cid:113) w )1+3 w , − (cid:113) w )1+3 w , − (cid:113) w )1+3 w (cid:21) (cid:3) (cid:3) Saddle E (cid:2) , , (cid:3) Repellor Repellor Repellor F (cid:2) − , − , − (cid:3) Attractor Attractor Attractor G + (cid:2) , , − w (cid:3) Repellor Repellor Repellor G − (cid:2) − , − , − + 2 w (cid:3) Attractor Attractor Attractor H + (cid:104) − w − , − w +79( w − , − w − w − (cid:105) Saddle Saddle Saddle H − (cid:104) w +79( w − , w − , w − w − (cid:105) Saddle Saddle Saddle J Saddle Saddle Saddle K + Unstablespiral Unstable Unstable K − Stable spiral Unstable UnstableTable II: In the table above we summarise the stability of each equilibrium point, for equation of state parameters correspondingto dust, radiation and a cosmological constant. A + and A − are non-hyperbolic and consequently we use the Center ManifoldTheorem to obtain their stability. Points C , D only exist physically in the phase space for w = −
1, and the eigenvalues forthese points are all equal to zero. The stability of these points were therefore determined by inspection – they are self-evidentlysaddle points. The points J and K ± lie on the plane y = v ; for these points Γ, and thus the system, is undefined and thereforethere are no analytic eigenvalues for these points. The classification of these points was therefore done by inspection. For thecases w = , −
1, it was only possible to determine whether or not the points were stable. the scale factor for each equilibrium point of the com-pact dynamical system. This can be done for the points A ± , B ± , G ± , H ± , and K ± .Fixed points with deceleration parameter q = − q = 1 represent universes which appear to be “radiation-like” , as the scale factor is proportional to the square Note that this “radiation-like” behaviour refers to the properties root of cosmic time.One of the short-comings of the dynamical systemsapproach to cosmology, as outlined in [11], is that itis possible for the analysis to admit fixed points whichcorrespond to solutions of the dynamical system (as de-fined, for example, by (30) - (32)) but do not satisfy thecosmological equations. In many cases, constants of in-tegration, which emerge in families of solutions to thecosmological equations, result in additional constraints,which must be satisfied by all physical points of the sys-tem. Setting the derivatives of the dynamical variablesequal to zero; x (cid:48) = F ( x ) = 0 (41)implies either F ( x ) = 0 (42)or x (cid:48) = 0 ⇒ x = constant . (43)Solutions to (41) may result from solving either of theequations (42) or (43), where the latter now represents aset of constraints imposed on the system [11]. For thisreason, it is important to verify that the solutions ob-tained satisfy the cosmological equations.The points C , D , E , F , J , K , all lie on the non-invariant Q = 0 submanifold, describing solutions for which thescale factor has no time dependence, and so representstatic universes.For the HS model with n = c = 1, we find that thenon-static, analytic fixed points belong to one of two scalefactor solutions: radiation-like or de Sitter like expan-sion. The fact that other cosmological evolutions do notappear as stationary phase states does not imply that thismodel does not allow them – it may be that the choice ofvariables places analytic limitations on what can appearas stationary solutions.The expanding versions of A and B are stable equilib-rium states, and correspond to de Sitter scale factor ex-pansion histories. Several interesting orbits exist, whichoriginate near an unstable “radiation-like” point, for ex-ample G + or H + , and evolve towards one of these accel-erated expansion points. In Figure 2 four such orbits arepresented. Extremely fine tuned initial conditions are re-quired for a trajectory to pass close to the non-analytic,unstable expanding matter point K + and evolve towardeither of the de Sitter late-time attractors. It is morenatural for trajectories to begin near one of the radiationlike points, G + or H + and asymptote towards one of thede Sitter points. This feature of this class of models willbe clarified in Section V. of the curvature fluid - at these points the curvature fluid causesthe expansion of the universe to scale with time in the same waythat ordinary radiation would: a ( t ) ∝ √ t , exhibiting “radiation-like” behaviour. Figure 2: Here two examples of solution orbits are presented.The orbit represented in the left panel begins near the radia-tion like repeller point G + , is pulled toward the saddle radia-tion like point H + and eventually evolves toward the de Sitterlike attractor B + . The orbit represented in the right panel alsobegins near G + , but evolves towards the de Sitter point, A + .The existence of these orbits show that this model naturallyadmits solutions which have a late time de Sitter expansion,which can produce expansion histories which look like a DarkEnergy fluid is dominating the universe at late times. F. Non-compact phase-space analysis of a HSmodel
Integrating the cosmological dynamical system, an ex-pansion history for the Hu-Sawicki model is obtained,which we can compare its predictions for the Hubble pa-rameter, deceleration parameter, total density and equa-tion of state parameter with those of the Concordancemodel. In order to calculate the expansion history, a noncompact phase space must be constructed, such that theHubble rate H is no longer a dynamical variable. We,therefore, make use of the following non-compact coordi-nates:˜ x = ˙ f (cid:48) f (cid:48) H , ˜ v = 16 RH , (44)˜ y = 16 ff (cid:48) H , ˜Ω m = 13 ρ m f (cid:48) H . (45)The coordinate transformation between compact andnon-compact variables is as follows:˜ x = 2 xQ , (46)˜ u = uQ , where ˜ u represents all of ˜ y, ˜ v and ˜Ω in turn.The dimensionality of the dynamical system has beenreduced as there is no dynamical variable which repre-sents the normalised volume expansion. The Friedmannequation at (7) can be reshuffled to give a constraintequation in terms of the above variables:1 = ˜Ω + ˜ v − ˜ x − ˜ y . (47)Since we are interested in integrating the system with re-spect to redshift, the differential equations correspondingto the non-compact dynamical system can be obtainedby differentiating the non-compact variables (44) withrespect to redshift, z : d ˜ xdz = 1( z + 1) (cid:104) ( − w ) ˜Ω + ˜ x + (1 + ˜ v ) ˜ x − v + 4˜ y (cid:105) , (48) d ˜ ydz = − z + 1) [˜ v ˜ x Γ − ˜ x ˜ y + 4 ˜ y − y ˜ v ] , (49) d ˜ vdz = − ˜ v ( z + 1) [(˜ x Γ + 4 − v )] , (50) d ˜Ω dz = 1( z + 1) (cid:104) ˜Ω ( − w + ˜ x + 2 ˜ v ) (cid:105) , (51)where Γ is still given by (38), due to the fact that therelationship between v and y , and ˜ v and ˜ y is preservedduring the coordinate transformation from compact tonon-compact variables.In terms of the dynamical variables, we also have anexpression for the dimensionless expansion rate of theuniverse, h , given by the Raychaudhuri equations, dhdz = hz + 1 [(2 − ˜ v )] , (52)where h = HH .The four interesting non-boundary points and their co-ordinates in the compact and non-compact phase spaces are tabulated below in Table III, along with their stabilityand their scale factor evolution, in a universe dominatedby dust. The purpose of this finite analysis is to pro-duce an expansion history for the particular HS modelunder investigation with which the ΛCDM model can becompared. It is therefore only important to consider the expanding version of the fixed points. That is to say, weneed only look at the non-compact fixed points corre-sponding to the expanding versions of the compact fixedpoints obtained earlier. V. EXPANSION HISTORY FOR THEHU-SAWICKI MODEL, WITH n = 1 , c = 1 As already pointed out, the utility of the finite dy-namical systems analysis is to provide a fast numericallystable integration of the cosmological equations linkingthe early and late time evolution (represented by fixedpoints in the phase-space of these models). This expan-sion history can then be compared to the ΛCDM modelin order to determine whether it provides a good alterna-tive to the standard model. We will see below that thisdepends critically on the choice of parameter values andthe initial conditions.
A. Initial conditions
The expansion history for a universe governed by theHS model with n = 1 and c = 1 is calculated, byperforming an integration of the differential equationsrepresenting the dynamical system (48) − (51). Themodel was parametrised based on the considerations of[31] wherein fiducial restrictions were placed on the rela-tionship between c and c : c c ≈ Λ Ω m , (53)to control the ratio of matter density to the cosmologicalconstant via the parameters c and c . As is done in [31],the following values for the densities at the present epochare used: Ω Λ = 0 . , Ω m = 0 . . (54)The initial values (at the present epoch) for the Ricciscalar, R , and the derivative of the function, f (cid:48) ( R ) ≡ f R ,can be given by R ≈ m (cid:18) m − (cid:19) ,f R ≈ − g R = 1 − n c c (cid:18) m − (cid:19) − n − , (55)for | g R | (cid:28)
1. Substituting (54) into (55), we have: R = 41 m = 41 λH . (56)0 Point Non-Compact(˜ x, ˜ y, ˜ v, ˜Ω) Compact( Q, x, y, v,
Ω) Stability ( w = 0) Scale factor solution˜ A + [0 , , ,
0] [ √ , , , ,
0] Attractor a ( t ) = a e H ( t − t ) ˜ B + [1 , , ,
0] [ , , , ,
0] Attractor a ( t ) = a e H ( t − t ) ˜ G + [ − , , ,
0] [2 , − , , ,
0] Repellor a ( t ) = a (2 H ( t − t ) + 1) ˜ H + [1 − w, , , − w ] (cid:104) w − , w − w − , , , −
49 3 w − w − (cid:105) Saddle a ( t ) = a (2 H ( t − t ) + 1) ˜ K + (cid:2) , , , (cid:3) (cid:104) √ , , , , (cid:105) Unstable spiral a ( t ) = a (cid:0) H ( t − t ) + 1 (cid:1) Table III: There exist two stable equilibrium points for which the scale factor increases exponentially into the future; ˜ A + and ˜ B + . It follows that these states can be associated with the late time accelerated expansion of the universe attributed toan effective f ( R ) Dark Energy. The scale factor evolution of the points, ˜ G + and ˜ H + depend on the square root of time, andtherefore represent “radiation” -like universes. The non-analytic matter point is included, as it also appears in the non-compactphase space on the plane ˜ y = ˜ v . In order to calculate the initial values for the dynamicalvariables (˜ x, ˜ y, ˜ v, ˜Ω), expressions for each are required interms of the quantities q , r and h (with r and h definedby (36)), for which the initial values can be obtained byusing (55) and (54). With n = 1, c = 1 we obtain: f (cid:48) ( r ) = c r ( c r + 2 λ )( c r + λ ) , r = 6(1 − q ) h , (57)from which we can obtain the deceleration parameter: q = 6 h − r h . (58)It follows that, in terms of the Ricci scalar, the Hubbleparameter, the deceleration parameter and c and C , thedynamical variables are:˜ v = 1 − q = r h , ˜Ω = Ω m ( z +1) h f (cid:48) ( r ) , ˜ y = r ( r c + λ ) h ( r c +2 λ ) , ˜ x = Ω + v − y − . (59)Ω denotes the initial value of the dynamical variable Ω,and is not to be confused with Ω m , which represents thematter density parameter of the universe at the presentepoch. Using (53) and (55), and the fact that m ≡ CH ,we find the following values for the parameters : n = 1 , c = 1 , c = 119 , λ = 0 .
24 (60)and the following initial values: h = 1 , r = 9 .
840 (61) f ( R ) theories commonly suffer singularities in the expansion his-tory, where the first derivative of the Hubble parameter diverges.The parameters chosen above enable a singularity free expansionhistory within the redshift range which we are interested. [ ? ]details theoretical constraints on the parameters of this model,as well as investigating certain observational constraints. ˜ x ,HS = − . , ˜ y ,HS = 1 . , ˜ v ,HS = 1 . , ˜Ω ,HS = 0 . . (62)The above set of coordinates correspond to the presentepoch , z = , as calculated from the parameter valuesand constraints outlined in [31].The values for the present epoch in the compact phasespace are found to be: Q ,HS = 0 . , x ,HS = − . , y ,HS = 0 . ,v ,HS = 0 . , Ω ,HS = 0 . . (63)Below, in Figure 3, an orbit in the 3D compactphase space is presented, which has its initial values at( Q , x , v ) as given by (63). B. Comparing the Hu-Sawicki Model ( n = 1 , c = 1 )with Λ CDM
Figure 4 shows the redshift evolution of the dimension-less Hubble parameter and the deceleration parameter forthe specific HS model considered above, in comparisonwith a ΛCDM model parametrised by the same valuesfor Ω m and Ω Λ as above, i.e. z = 0 or “today”.
1. Hubble parameter, h The dimensionless Hubble rate for the ΛCDM modelis givel by: h ( z ) = (cid:112) Ω m (1 + z ) + Ω Λ . (64)and is compared to the solution of the differential equa-tion for h ( z ) given by (52).1 Figure 3: The fixed points represented by A + , B + , G + , H + and K + in the above plot correspond to the compact fixedpoints given in Table I. The matter point K + lies on the plane y = v for which Γ is undefined. The crossed square indicatesthe point which corresponds to the present epoch as given by(63), today ( z = 0), for the model defined by the parametersdefined in [31]. It can be seen that this point is close to thede Sitter stationary solution A + . From (cid:2) the orbit evolvesforward in time toward A + . In its past, it passes by theunstable radiation-like stationary state, H + . Preceding thepoint H + , the orbit evolves from the unstable static universephase state E .
2. The deceleration parameter, q The deceleration parameter for the ΛCDM model isgiven by the following expression: q = 1 H (cid:18)
12 Ω m (1 + z ) − Λ (cid:19) , (65)where as the deceleration parameter calculated for thisspecific HS model is given in terms of the dynamical vari-ables as q = 1 − ˜ v. (66)
3. Total equation of state parameter w T It is interesting to consider the redshift evolution ofthe total equation of state parameter, w T : w T = p T ρ T = p eff ρ m + ρ eff . (67)Using the trace and Friedmann equations, this can bewritten in terms of the dynamical systems variables de-fined in (44) above: w T = 13 (1 − v ) (68)The right panel in Figure 4 shows the behaviour of thetotal equation of state parameter. It is interesting tonote that this total equation of state parameter asymp-totes toward the value of w = , indicating a universewhich is dominated, for most of its history, by a cur-vature fluid exhibiting “radiation-like” behaviour , andonly now , at very low redshifts does a change in the formof the dominant energy density take place. This result iscounterintuitive, owing to the fact that the entire analy-sis was performed with respect to a dust only universe.However, it is consistent with the dynamical systemsanalysis which have unstable radiation like phase statesfrom which the corresponding solution orbit begins. Thismodel produces a late time negative equation of statewhich is consistent with the dynamical systems analysisshowing an approach toward a universe having scale fac-tor which evolves exponentially with time. However, thisresult is not consistent with idea that the HS model ap-proximation holds throughout the expansion history ofthe universe, as it predicts that a radiation fluid energydensity dominates even at relatively low redshifts. Thisis a highly unsatisfactory match to observations and, ofcourse, ΛCDM predictions. This result is specific to themodel considered, where n = 1 , c = 1. Choosing the ini-tial parameter values to be exactly equal to their ΛCDMvalues today ( z = 0) is what compromises the agree-ment between this model and the ΛCDM model. To ob-tain a better match, an appropriate adjustment of theseinitial values from their corresponding ΛCDM values isrequired. In order to understand why this is so, considera plot of the correction g ( R ) generated by the model withinitial values specified at (61) for z = , i.e. the ΛCDMvalues corresponding to “today”, in Figure 5.In order for the model investigated to mimic ΛCDMbehaviour, that is, exhibit a GR+ cosmological constant nature, the initial value of the Ricci scalar should lie onthe plateau corresponding to a constant value of g ( R ).However, for the specific model and density parametersabove, the initial value of R/H at z = 0, the presentepoch, is 9 .
84 as stated in (61). This value of
R/H is not high enough to place the correction initially onthe plateau, and therefore it does not allow the model( n = 1 , c = 1) to mimic ΛCDM behaviour. In fact,as can be seen in Figure 5, the Ricci scalar determinedby this model, for initial parameter values given by (61),2 Figure 4: Above, we plot in red the Hu-Sawicki model ( n = 1 , c = 1 , c = 1 / , C = 0 . z = 0, and in blue, the ΛCDM model. In (a) the dimensionless Hubble parameter for ΛCDM given by equation (64) iscompared with the dimensionless Hubble rate obtained from the integration of the equations extracted from the non-compactdynamical systems analysis (the solution to the dynamical system (48) − (52)). At very low redshifts, the two models coincidefor a short interval, but after around z = 0 .
05 they deviate substantially. In (b), the deceleration parameter for ΛCDM given byequation (65) is compared with the deceleration parameter determined by the solution of the dynamical system at (48) − (52),given by equation (66). The deceleration parameter determined by this specific HS model supports late time acceleration aswell as a stabilised deceleration at higher redshifts. However, there is a large discrepancy in the values for q predicted by theHS model and that predicted by ΛCDM for this choice of parameters. (c) shows the redshift evolution of the total equation ofstate parameter determined by the HS model. It exhibits odd behaviour; at high redshifts it tends towards a constant value of indicating radiation domination-like behaviour, even though the entire analysis was performed assuming a dust filled universeof w = 0. At low redshifts it decreases sharply to indicate the domination of a fluid with negative pressure, which reachesa minimum at -1.37 and begins to increase towards the present epoch to about -0.76. This is consistent with the dynamicalsystems analysis which shows an early time radiation-like repeller, no matter point, and a late time stable de Sitter phase state,as a result of the fluid with equation of state w ≈ − . decreases with redshift. This indicates that integratingfrom z = 0 only drives the value of g ( R ) further awayfrom its ΛCDM plateau limit. C. Initial conditions at z = 20 To remedy this, we can choose initial values which in-creases the value of r so that it sits comfortably on theplateau of g ( R ), therefore generating an effective cosmo-logical constant. This can be easily done by fixing ourinitial redshift to be z = 20 - see Figure 6. In this case,the initial values for the Ricci scalar, R and the non-compact dynamical variables are : h = 47 . , r = 6677 . , ˜ x ,HS = 0 . , ˜ y ,HS = 0 . , ˜ v ,HS = 0 . , ˜Ω ,HS = 0 . . (69)The above initial values are extremely close to the mat-ter point K + , which is consistent with the fact that, athigher redshifts, we expect the dust fluid will dominatethe equation of state of the universe. To illustrate this point, the solution trajectory with initial values given by(69) is presented in Figure 7.At this point, it is worth noting that the value of z =20 is sufficiently high in redshift to begin the integrationof the dynamical system. As we increase redshift, thevalues of the coordinates tend closer and closer to thematter point K + , more and more quickly, to such anextent that it makes little difference whether we beginat z = 1000 or z = 20, as we have chosen. Figure 8demonstrates this fact, showing the dynamical variablesas a function of redshift. The asymptote of the variablestoward their respective matter point, K + , coordinates isobvious.To demonstrate the effects of placing the initial valueof the correction initially on its function’s plateau, thecorresponding plots of the dimensionless Hubble param-eter, the deceleration parameter and the equation of stateparameter is presented in Figure 9.It is clear that placing the correction g ( R ) initially on the constant valued plateau permits the specific HSmodel analysed to mimic the ΛCDM model relativelywell. By doing this, we have corrected the asymptote3 Figure 5: The main plot shows the behaviour of the correc-tion. To fulfil its purpose, g ( R ) is constructed to tend towardzero as r → r in order to mimic the observedcosmological constant behaviour, which is well described, sofar, by the ΛCDM model [31]. The plateau of constant g ( R )offers a simulated cosmological constant term to the gravita-tional Lagrangian. The dotted black line indicates the initialvalue of the Ricci scalar, R/H ( z = 0) . The inner plot showsthe Ricci scalar defined as r = R/H as a function of redshift.Figure 6: We show the correction, when initial conditions arefixed at z = 20. In this case the initial value of g ( R ) sits onthe constant valued plateau of g ( R ) which represents a sim-ulated cosmological constant.The dotted black line indicatesthe initial value of the Ricci scalar, R/H ( z = 20) . Theinserted panel gives the redshift evolution of the expansionnormalised Ricci scalar r = R/H , which is now increasingwith redshift as expected, leading to a past expansion historyvery close to ΛCDM. of the equation of state from to 0, we have man-aged to simulate the Hubble parameter, h , of the ΛCDM Figure 7: The above plot shows the solution trajectory of thecompact dynamical systems equations, integrated from theinitial values specified at (69) for z = 20. The crossed circlesymbol represents the starting point ( z = 20) and the crossedsquare symbol represents the coordinate values at the presentepoch ( z = 0). The fixed point K + represents the matterfixed point at ( x, y, v, Ω) = (0 , , , A + represents a de Sitter point. This trajectory clearly beginsextremely close to the matter point, is affected by the spiralnature of the neighbourhood of K + on the Q - v plane andthen evolves neatly toward the de Sitter point. This showsan expansion evolution which is similar to ΛCDM, in thatit evolves from a state which is matter dominated to a statewhich is filled with a fluid causing exponentially acceleratedlate time expansion. model very closely, and there is an improvement in theearly time behaviour of the deceleration parameter. Fromthese results it can be concluded that the HS model canproduce a very good simulation of the ΛCDM model,however, this occurs at the expense of the present ΛCDMvalues of q , h and r . In order to obtain the desired be-haviour; for the cosmology to be dominated by a dustfluid in the past, the values of q Λ CDM , h Λ CDM and r Λ CDM corresponding to the present epoch ( z = 0) cannot be used as starting values for the integration. Theseparameters must be adjusted to enable g ( R ) to initiallyassume its plateau value.4 Figure 8: The (˜ x, ˜ y, ˜ v ) coordinates of the matter points, K + are (cid:0) , , (cid:1) , indicated on the plots by black dotted line. Thedifference between these values and the initial values of the dynamical variables ˜ x, ˜ y and ˜ v is too small to resolve at a redshiftgreater than 20. We therefore assert that beginning the integration at this redshift is sufficient to determine the behaviour ofthis model in the ΛCDM regime.Figure 9: Here we provide a comparison of the Hu-Sawicki model with n = 1 , c = 1 , c = 1 / , C = 0 .
24 (red dashed line) andthe Λ
CDM model (64) (solid blue line), where initial conditions are fixed at z = 20. The left panel gives the dimensionlessHubble parameter which is almost indistinguishable from Λ CDM . The middle panel gives the deceleration parameters, givingboth late time acceleration today and q HS ⇒ q Λ CDM ≈ at high redshifts. The right panel gives a total effective equation ofstate which clearly asymptotes towards w Total = 0 at high redshifts (as expected) as well as tending toward a negative value inthe low redshift regime, thus behaving like dark energy. These results represent a dramatic improvement on what was obtainedfor initial conditions at z = 0, shown in figure 4, which was unable to produce a matter dominated evolution at high redshift.
1. Effects of increasing n Note that the parameter n plays an important rolein the derivative of g ( R ). Large values of n result in a steeper , slope of the transition between the limiting val-ues of g ( R ). Therefore increasing n − making the tran-sition from g ( R ) → g ( R ) → Constant more rapid − could place the initial value for g ( R ) at the presentepoch ( z = 0) on the constant plateau of the correctionfunction. Figure 10 shows a plot of the correction, g ( R )for n = 3. In this plot, it is clear that the initial value of r at z = 0 results in a value of g ( R ) which sits well onthe constant part of the function, enabling a model whichclosely resembles the ΛCDM model at high curvature. Figure 10: Above, is a plot of the correction g ( R ) for n =3. It shows a significantly steeper slope, enabling a morerapid transition from the GR limit to the ΛCDM limit of thefunction. This slope enables the HS model (for n > z = . However, n (cid:54) = 1 is not permitted bya dynamical systems analysis performed with the variablesdefined by (23) or (44), therefore other means must be usedto compute its resulting expansion history. The inserted panelshows the behaviour of the corresponding Ricci scalar. VI. CONCLUSIONS
In this paper we present for the first time a detailedanalysis of the cosmological dynamics of the Hu-Sawicki(HS) model. In order to achieve this, we first expressedthe fourth-order cosmological equations in terms of a setof generalised dimensionless expansion normalised dy-namical variables which, span the phase space of thisproblem. While these equations are completely generalfor any f ( R ) theory of gravity, the choice of theory isfixed by expressing the term Γ (which appears in theequations) in terms of the set of dynamical variables,thus closing the system of equations. In the case of the HS model, this can only be done for the particular choiceof model parameters: n = 1 and c . While this may ap-pear somewhat restrictive, it is easy to show using othernumerical experiments, that this choice does not lead toa loss of generality. Indeed, choosing n >
1, simply leadsto an improved fit to the standard ΛCDM expansion his-tory [53].A compact dynamical systems analysis was performedusing a positive definite normalisation, which brings anyequilibrium points at infinity to the boundary of the re-gion defined by the range of the dynamical systems vari-ables. From this analysis, twelve equilibrium points wereidentified, each having an expanding (
Q >
0) and con-tracting (
Q <
0) version. Four de Sitter like stationarystates were found A ± and B ± , two within the expandingpart of the phase space, with the other two correspond-ing to the collapsing versions. Focusing on the expand-ing versions of these points, both A + and B + are stable stationary phase states. Two unstable “radiation like”states - G + and H + exist in the expanding sector of thephase space. There are several orbits which connect these“radiation” points to the de Sitter phase states, offeringtrajectories which resemble the chronological evolution ofthe scale factor of our universe.On the surface y = v , a very interesting, non-analytic, matter-like point, K + , was identified which after carefulanalysis, turns out to be an unstable spiral point. Theredoes exist a trajectory which evolves from K + toward ade Sitter phase state at A + , the existence of which indi-cates that the Hu-Sawicki model provides a modificationto gravity able to produce expansion histories consistentwith the ΛCDM model.In order to consider the expansion history generatedby the HS model, a finite dynamical systems analysiswas performed to confirm that the non-boundary pointsobtained in the compact analysis were, in fact, finitepoints. The finite points obtained in this analysis arethe analogies of the de Sitter and “radiation like” points, A ± , B ± , G ± and H ± . It was found that for n = 1, the expansion historygenerated from the initial values and parameter valuesspecified in [31] for the present epoch resulted in a uni-verse which is dominated until very low redshifts by anequation of state parameter equal to , implying thatfixing the initial conditions to be exactly those of theΛCDM conditions today leads to an expansion historyvery different from ΛCDM and high redshifts ( z > n = 1, the slope of thecorrection g ( R ) defined by the HS model does not al-low for the initial value of g ( R ) to sit on the constantvalued plateau of the function. It is this plateau whichmimics the late time cosmological constant behaviour ofthe model. We showed that by setting initial conditionsto be the same as ΛCDM model at high redshift, leadsto a value of g ( R ) which is equal to its “plateau” valueand consequently leads to an expansion history which ismuch closer to the ΛCDM model. Moreover it produces adimensionless Hubble parameter that is nearly indistin-6guishable from ΛCDM, with deviations only appearingin the dynamics of the deceleration parameter. Most im-portantly, it gives an equation of state parameter whichbehaves like CDM for most of the expansion history, andbegins to tend toward -1 to indicate the domination ofan effective cosmological constant term, at low redshifts.Unfortunately for values of n different from 1, the dy-namical systems formalism presented here does not leadto a set of autonomous equations which can be analysedin the same way. However, we expect that the key fea-tures found in the n = 1 case to remain the same. Neverthe less finding a set of variables which is able to facili-tate the dynamical systems analysis a general HS modelin order to perform a complete analysis of the general phase space of these models remains an interesting prob-lem which has been partially addressed in a recent paper[65]. Acknowledgments
We would like to thank David Bacon for a comprehen-sive reading of the manuscript and for his useful com-ments. S.K. is grateful to the NRF and the Faculty ofScience, University of Cape Town for financial support.P. K. S. D. thanks the National Research Foundation(NRF) for financial support. [1] Abdelwahab, M et al. Class. Quantum Grav. et al. Phys. Rev. D An Introduction toDynamical Systems . Cambridge University Press, 1994.Print.[4] Barrow, J. D. and Ottewill, A. C. (1983)
J. Phys. A:Math. Gen. Beyond Einstein Gravity:A Survey of Gravitational Theories for Cosmology andAstrophysics . New York: Springer, 2011. Print.[7] Cardone, V. F. et al. JCAP02 (2012)030 arXiv:1201.3272[astro-ph][8] Calcagni, G. et al. Nucl. Phys. B (2006) 404-438[9] Carloni, S. et al. Class. Quantum Grav. et al. Gen.Rel.Grav. (2009) 1757-1776arXiv:0706.0452 [gr-qc][11] Carloni, S. et al. (2008) arXiv:0812.2211 [astro-ph][12] Carloni, S. et al. J. Phys. A: Math. Theor. An Introduction to Mod-ern Astrophysics, nd Ed.
Pearson Education, 2006.Print.[14] Carroll, S. M. et al. Phys.Rev. D (2005) 063513[15] Chevallier, M. and Polarski, D. Int.J.Mod.Phys. D (2001) 213 arXiv:0009008 [gr-qc][16] Chiba, T. (2003) Phys.Lett. B , 1-3 arXiv:0307338[astro-ph][17] Clifton, T. et al. Phys.Rept. (2012) 1-189arXiv:1106.2476 [astro-ph.CO][18] de Block, W. J. G. Adv.Astron.
Class. Quantum Grav. Dynamical Theory of Groups and Fields .New York: Gordon and Breach, 1965. Print.[21] Dodelson, S.
Modern Cosmology.
Academic Press, 2003.Print[22] Dunsby, P. K. S. et al. (2010) arXiv:1005.2205 [gr-qc][23] Faraoni, V. (2008) arXiv:0810.2602 [gr-qc][24] Fleury, P. et al. Phys.Rev. D (2013) 123526 arXiv:1302.5308 [astro-ph.CO][25] Gannon, T. Nucl.Phys. B (2003) 335-358arXiv:0305070 [hep-th][26] Goheer, N. et al. Class. Quantum Grav. et al. Class. Quantum Grav. Phys.Rev. D (1999)023502 arXiv:9811068 [gr-qc][29] Hobson, M. P. , Efstathiou, G. P., Lasenby, A. N. GeneralRelativity : An Introduction for Physicists . CambridgeUniversity Press, 2006. Print.[30] Hu, W.
Annals Phys. (2003) 203-225 arXiv:0210696[astro-ph][31] Hu, W. et al . Phys.Rev. D (2007) 064004arXiv:0705.1158 [astro-ph][32] Igari, K. Publ. RIMS, Kyoto Univ (1974) 613-629.[33] Kehagias, A. (2013) arXiv:1312.1155 [hep-th][34] Khoury, J. and Weltman, A. (2004) Phys.Rev.Lett. ,171104 arXiv:0309300 [astro-ph][35] Kofman, L. (1996) arXiv:9605155 [astro-ph][36] Langlois, D. Lect.Notes Phys. (2010) 1-57arXiv:1001.5259 [astro-ph.CO][37] Liddle, A.
An Introduction to Modern Cosmology, nd Ed.
John Wiley & Sons, 2003. Print.[38] Linder, E. V.
Phys.Rev.Lett (2003) 091301arXiv:0208512 [astro-ph][39] Martinelli, M. et al . Phys.Rev. D (2009) 123516arXiv:0906.2350 [astro-ph][40] Mizohata, S. et al. Lectures on Cauchy Problem Bombay,1965.[41] Mota, D. F. et al. (2004)
Mon.Not.Roy.Astron.Soc. ,291 arXiv:0309273 [astro-ph][42] Peebles, P. J. E
Principles of Physical Cosmology.
NewJersey: Princeton University Press, 1993. Print.[43] Peebles, P. J. E. and Ratra, B.
Rev. Mod. Phys. (2003)559-606[44] Percival, I. and Richards, D. Introduction to Dynamics .Cambridge: Cambridge University Press, 1994. Print.[45] Perko, L.
Texts in Applied Mathematics 7: DifferentialEquations and Dynamical Systems
New York: Springer,2001. Print.[46] Perlmutter, S. et al Nature (1998) 51-54[47] Planck collaboration (Ade, P.A.R. et al. ) arXiv:1303.5076 [astro-ph.CO][48] Planck collaboration (Ade, P.A.R. et al. )arXiv:1303.5062 [astro-ph.CO][49] Planck collaboration (Ade, P.A.R. et al. )arXiv:1303.5082 [astro-ph.CO][50] Planck collaboration (Ade, P.A.R. et al. )arXiv:1303.5075 [astro-ph.CO][51] Rajagopal, S. and Kumar, A. (2013) arXiv:1303.6026 [gr-qc][52] Saito, K. and Ishibashi, A. PTEP (2013) 013E04arXiv:1209.5159 [gr-qc][53] Santos, B. et al . A & A A31 (2012) arXiv:1207.2478[astro-ph][54] Scarpa, R.
AIP Conf.Proc. (2006) 253-265arXiv:0601478 [astro-ph][55] Schmidt, F. et al. Phys.Rev. D (2009) 083505arXiv:0908.2457 [astro-ph.CO][56] Schneider, P. Extragalactic Astronomy and Cosmology:An Introduction.
Berlin: Springer, 2006. Print [57] Shafieloo, A. et al. Phys. Rev. D (2011) 063519arXiv:1107.1033 [astro-ph.CO][58] Sotiriou, T. P. and V. Faraoni (2010) Rev.Mod.Phys. ,451-497 arXiv:0805.1726 [gr-qc][59] Sotiriou, T. P. (2010) Phys.Lett. B , 225-228arXiv:0805.1160 [gr-qc][60] Sotiriou, T. P. (2007) arXiv:0710.4438 [gr-qc][61] Starobinsky, A. A. JETP Lett. (2007) 157-163arXiv:0706.2041 [astro-ph][62] Starobinsky, A. A. astro-ph/9605155 UH-IFA-96-28[63] Wainwright, J. and Ellis, G. Dynamical Systems in Cos-mology . Cambridge: Cambridge University Press, 1997.Print.[64] Weinberg, S.
Gravitation and Cosmology: Principles andApplications of the General Theory of Relativity
JohnWiley & Sons, 1972. Print.[65] S. Carloni,