The phase-space analysis of modified gravity (MOG)

We investigate the cosmological consequences of a scalar-vector-tensor theory of gravity known as modified gravity (MOG). In MOG, in addition to metric tensor, there are two scalar fields G(x) and μ(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu (x)$$\end{document}, and one vector field ϕα(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{\alpha }(x)$$\end{document}. Using the phase space analysis, we explore the cosmological consequences of a model of MOG and find some new interesting features which are absent in Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda $$\end{document}CDM model. More specifically we study the possibility that if the extra fields of this theory behave like dark energy to explain the cosmic speedup. More interestingly, with or without cosmological constant, a strongly phantom crossing occurs. Also we find that this theory in its original form (Λ≠0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda \ne 0$$\end{document}) possesses a true sequence of cosmological epochs. However, we show that, surprisingly, there are two radiation-dominated epochs, f5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_5$$\end{document} and f6\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_6$$\end{document}, two matter-dominated phases, f3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_3$$\end{document} and f4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_4$$\end{document}, and two late time accelerated eras, f12\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{12}$$\end{document} and f7\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{7}$$\end{document}. Depending on the initial conditions the universe will realize only three of these six eras. However, the matter-dominated phases are dramatically different from the standard matter-dominated epoch. In these phases the cosmic scale factor grows as a(t)∼t0.46\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a(t)\sim t^{0.46}$$\end{document} and t0.52\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t^{0.52}$$\end{document}, respectively, which are slower than the standard case, i.e. a(t)∼t2/3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a(t)\sim t^{2/3}$$\end{document}. Considering these results we discuss the cosmological viability of MOG.


Introduction
We investigate the cosmological consequences of a modified theory of gravity known as MOG in the relevant literature [1]. MOG is a relativistic theory which exploits three kinds of gravitational fields, i.e. tensor, scalar and vector fields. More specifically, in addition to the metric tensor, MOG possesses two scalar fields G(x), μ(x), and a Proca vector field φ α (x). The vector field is directly coupled to the matter fields. Therefore, this theory is not a metric theory of gravity and consea e-mail: sara.jamali@stu.um.ac.ir b e-mail: mroshan@um.ac.ir quently the weak equivalence principle, in principle, can be violated. Naturally, the free parameters of MOG are chosen such as to make the theory consistent with the experimental tests of the equivalence principle. The main motivation for introducing this theory is to solve the dark matter enigma. It is claimed that MOG can explain the flat rotation curve of the spiral galaxies without adding any dark matter halo [2,3]. Also this theory explains the matter discrepancy in the galaxy clusters [4]. It is worthy to mention that it is not the first time that some modifications in the gravitational law can address the above mentioned problems. For an explicit example we refer the reader to Modified Newtonian Dynamics (MOND) [5] and its relativistic generalizations such as Tensor-Scalar-Vector theory (TeVeS) [6]. It has recently been claimed that MOG is more successful than MOND in explaining the flat rotation curves [7]. Also the local stability of spiral galaxies in MOG has been investigated in [8]. The gravitational Jeans instability for molecular clouds has been studied in [9].
Our purpose in this paper is to study the cosmological behavior of a MOG model. It is important mentioning that like f (R) gravity, MOG may refer to a large class of models corresponding to different energy contributions for the scalar and vector fields. In other words by changing the kinetic and potential energy contributions of the fields, one may construct a new model of MOG. In this paper we restrict ourselves to a MOG model presented in [10].
The astrophysical consequences of this theory, more specifically the astrophysical issues relevant to the dark matter problem, have been widely investigated. Since MOG is still considered as an alternative theory to dark matter particles and has not been ruled out yet, it seems necessary to check its cosmological consequences. We know that adding only a single scalar field to a gravitational theory can lead to significant outcomes in the cosmological issues. For example we recall the quintessence model and Brans-Dicke theory. Therefore it is natural to ask about the cosmological behavior of MOG considering the variety of fields that have been incorporated. On the other hand there are a few papers con-sidering the cosmology of MOG. For example in [10] the Noether symmetries of the cosmic point-like Lagrangian of MOG have been studied and some exact cosmological solutions have been found. Also in [11] the perturbation growth in the context of MOG has been studied. See [12][13][14] for relevant work.
In order to check the main cosmological features of MOG, and the minimum requirements that it must possess, we use the dynamical system method (or the phase-space analysis). This method provides a fast and reliable procedure to numerically solve the field equations. Note that the field equations of MOG are drastically more complicated than Einstein's general relativity (GR); see Eqs. (12)- (15). More importantly, this method enables us to project the dynamics into a compact region and explore the most important "'events"' that can be happen. One of the necessary requirements that a cosmological model should satisfy is the existence of a true sequence of cosmological epochs. More specifically, the cosmic evolution should start with a radiation-dominated phase. After this phase there should be a proper matter-dominated phase which is long enough to allow for the structure formation and fast enough to be consistent with the observations of the age of the universe. Finally the universe should enter an accelerated epoch consistent with the relevant observations such as the Supernovae type Ia data. Fortunately, the dynamical system method is an excellent tool for checking this important requirement. This method has been applied to various alternative theories and cosmological models, for example see [15][16][17][18][19][20][21][22]. This paper is organized as follows. In Sect. 2 we briefly introduce MOG and the modified Friedmann equations. In Sect. 3 we introduce the dynamical system variables and the autonomous first order differential equations. Also we find the critical points and explore their stability and physical relevance. In this section we assume that the cosmological constant is zero. In Sect. 3.2 we bring back the cosmological constant and analyze the system. In Sect. 4 we study the phase space of the system at infinity. Finally, conclusions are drawn in Sect. 5.

Modified Friedmann equations in MOG
Let us start with an action for MOG presented in [10], where R is the Ricci scalar, is a positive constant corresponding to the cosmological constant in the Einstein-Hilbert action. It is noteworthy that, although in this paper we will denote by the cosmological constant, it can be considered as the mass term for the scalar field χ . In other words, it is not exactly the cosmological constant and one may present different interpretations for its appearance in the action. Also ω 0 denotes a positive coupling constant, S M is the matter action, and B μν = ∇ μ φ ν − ∇ ν φ μ is an antisymmetric tensor reminiscent of the Maxwell tensor in electrodynamics. The new scalar fields χ and ψ are related to G and μ introduced in [1] as χ 2 = 2/G and ψ = ln μ; see [10] for more details. It should be stressed that the scalar field G, in principle, can be negative. This means that χ can be a pure imaginary function. However, the Lagrangian density in the above action remains always real. However, χ is an auxiliary function for writing the action in a more common and compact form, and our main scalar field is G. The potential V φ is chosen as V φ ∝ e 2ψ φ β φ β . This means that μ appears as a time dependent mass for the vector field and plays a central role for addressing the dark matter problem [1].
Varying the action with respect to the fields, one can find the relevant field equations, Here G μν is the Einstein tensor and J α is a "fifth force" matter current defined as a nonzero J α means that there is a coupling between matter and the vector field φ μ . This coupling can, in principle, lead to a violation of the Einstein equivalence principle. In this paper we assume that ∇ α J α = 0. This is an extra assumption and in principle one may study different versions of MOG in which this conservation equation is violated. Also, the total energy-momentum tensor is defined as where T μν is the energy-momentum tensor for ordinary matter, and 32π ∇ μ ψ∇ ν ψ − 1 2 g μν ∇ α ψ∇ α ψ . In order to study the cosmological consequences of MOG, we assume a flat Friedmann-Robertson-Walker (FRW) metric where a(t) is the cosmic scale factor. Also we assume that the cosmic fluid can be characterized by an ideal fluid with energy density distribution ρ, the pressure p and the velocity four vector u μ . In this case the energy-momentum tensor is Finally, bearing in mind that χ 2 = 2/G and ψ = ln μ, we find the Friedmann equations: Here a dot stands for a derivative with respect to time t, and ρ includes both matter and radiation contributions, i.e. ρ = ρ m + ρ r . It is important to mention that the scalar fields μ and G have negative contributions to the total energy density. More specifically the kinetic termsμ 2 μ 2 andĠ 2 G 2 appear with negative sign in (7). As we will show, this fact leads to some phantom features in this model. It is obvious that if we change the sign of the kinetic terms in the action, then the cosmological consequences of this model, in principle, will change.
In the following we restrict ourselves to the potential V φ = − 1 2 μ 2 φ α φ α . This is the original potential of MOG presented in [1]. In this case after some algebraic manipulations, we rewrite Eqs. (7)-(11) as follows: Note that using the field equation of the vector field we have replaced φ 0 , the only nonzero component of the vector field, with 16π J (t) ω 0 μ 2 , where J = J 0 is the time component of the matter current J α . It is clear that these non-linear differential equations, i.e. Eqs. (12)- (15), are drastically more complicated than the standard Friedmann equations. However, as we shall see, despite this complexity the dynamical system approach provides a fast numerically stable integration of the equations.
As we have already mentioned the equivalence principle can be violated in this theory. Consequently the ordinary energy-momentum tensor T μν is not conserved [23]. However, fortunately in the isotropic and homogeneous FRW space-time and with the assumption that ∇ α J α = 0, T μν is conserved and one may use the standard relations between energy densities and the scale factor, i.e. ρ m ∝ a −3 and ρ r ∝ a −4 ; see [23] for more details. In this case one may straightforwardly set the matter current as J = κρ m . Here κ is another positive coupling constant. It is noteworthy that in a non-homogeneous space-time because of the coupling between matter and the vector field, these relations are not true and one may expect significant departures from the CDM model. We recall that there are several attempts in the literature to find a relationship between the cosmic-speed up and the inhomogeneities in the distribution of matter, for example see [24,25]. Therefore, regarding the energy exchange between matter and the vector field in a nonhomogeneous background, it seems interesting to check the possibility that MOG can explain the accelerated expansion without invoking the cosmological constant just by taking into account the matter inhomogeneities. This issue can be a matter of study for future work. Therefore, in the following we work in an isotropic and homogeneous background. Now let us consider MOG as a dark energy model. In order to find the equation of state parameter of dark energy, i.e. ω DE , we rewrite Eqs. (12) and (13) as where G N is the Newtonian gravitational constant and H = a/a is the Hubble parameter and a dot denotes derivative with respect to cosmic time t. Here ρ DE and p DE are defined as Now it is possible to express the equation of state parameter of dark energy by noticing that ω DE = p DE ρ DE . It is also useful to write the effective equation of state parameter ω eff , which conveniently is defined to include all components of the energy budget of the cosmos, namely

Phase-space analysis of MoG
In order to apply the phase-space analysis to MOG, we transform the field equations (12)-(15) into the autonomous form , where x is the column vector constituted by an appropriate set of new variables and f(x) is the corresponding column vector of the autonomous differential equations. Also a prime denotes the derivative with respect to ln a. The fixed points x c of the system satisfy x = 0, and in order to determine the stability of these points, we perturb the system around the fixed points as x = x c +δ, where δ is a column vector for the perturbations. Expanding the autonomous equations up to linear order in perturbations, we have δ = M δ where M is the stability matrix. Finally, the type and the stability of each fixed point, can be found using the eigenvalues of the stability matrix [26]. Now let us define the following dimensionless variables: Substituting these dynamical variables into Eq. (12), we find a constraint equation, If we assume that G > 0, then it is clear form the definition of Q that it is a positive parameter. Therefore we rewrite the constraint equation as After some algebraic manipulations, the modified Friedmann equations take the following form: It is important mentioning that although G (or equivalently χ ) can be considered as a time dependent gravitational constant, its sign is not necessarily positive. This means that, in principle, anti-gravity is possible in MOG. In fact it is well known that in non-minimally coupled scalar-tensor theories of gravity, the anti-gravity regime can exist; see [27] and the references therein. For more recent work we refer the reader to [28][29][30]. There is also a non-minimally coupled scalar field χ in MOG. However, in the following we explicitly show that there is no transition from anti-gravity to gravity in the context of MOG. More specifically, if the evolution starts from an anti-gravity regime, it will remain permanently at that phase. In other words, if G starts with a negative value, then its sign will not change during the cosmic evolution. Therefore, it has to start from a positive value and one can be sure that y and r are also positive quantities during the whole thermal history. To show this fact more precisely, let us rewrite Eq. (12) as In the early universe we can neglect the term. In this case, supposing a negative value for G, one finds thatĠ is also negative for an expanding universe. However, there may exist a minimum for G and after that it can increase and finally become positive. Note that there is a positive term in the right hand side of (28). Therefore, in principle, there is a point where the term dominates and soĠ = 0 andG > 0; see also Eq. (14). However, one may naturally expect that becomes important only at the late times. Therefore it is very unlikely that we have an anti-gravity to gravity transition at the early stages of the universe. Finally we deduce that y and r are positive quantities. More specifically, we shall show that if the evolution starts with suitable initial conditions including a positive G,Ġ can be negative or positive during the cosmic evolution but G remains positive.
Using the introduced dynamical variables, the equation of state parameter of dark energy, ω DE and the effective equation of state parameter, ω eff , take the following simple form: Although ω eff can be written just in terms of phase space coordinates, ω DE contains a new variable β = G N G . It is clear form the definition of ω DE that, in principle, the denominator can become zero. Consequently this parameter can become infinite. More specifically, this is the case for some initial conditions considered in this paper; for example, see Fig. 2.

MOG without cosmological constant ( = 0)
As we have already mentioned, one of the main purposes of the current paper is to check if MOG can be considered as a dark energy model. To investigate if extra fields of MOG can play the role of dark energy, we set the cosmological constant to zero in the autonomous differential equations. It is equivalent to setting x = 0 in the equations of motion. In this case, the critical points (y, r, m, z) of Eqs. (23)-(26) are listed thus: Surprisingly the fixed points are numbers and there is no free parameter to be constrained. This point is also clear from Eqs. (23)- (27), where because of the special choice of the dynamical variables, free parameters do not appear in the autonomous differential equations. In fact, in principle, the free parameters appear in the coordinates of the fixed points [15][16][17][18][19][20][21][22]. Since each fixed point corresponds to an exact solution for the fields of the theory, existence of the free parameters in the fixed points provides an opportunity to make the model more consistent with the cosmological observations. However, this is not the case in MOG.
In this sense MOG behaves like the CDM model where the fixed points are ( R , m , ) = (1, 0, 0), (0, 1, 0) and (0, 0, 1), where the are cosmic density parameters. Therefore, as we shall see MOG leads to clear cosmological consequences as in the case of the CDM model.
p 1,2 : G-μ dominated curves p 1 (m > 0) and p 2 (m < 0) correspond to two distinct curves in the phase space. Every point on these curves is a fixed point. In this case there is a constraint on z as 1.046 < z < 22.95 to keep m real. These curves cover a wide range of ω eff , from non-phantom, to slightly phantom and strongly phantom as −14.3 ≤ ω eff ≤ 0.3. Eigenvalues of the stability matrix for p 1 and p 2 are , respectively. It is easy to show that for p 1 one of the eigenvalues is positive for z < 3. Therefore, every point on the curve p 1 in this interval is unstable. On the other hand, in the case of p 2 , for points in the interval 1.072 < z < 14.93, there is at least one positive eigenvalue. Thus, these points are unstable. Note that, for 3 < z < 22.95, the signs of the eigenvalues of p 1 are negative and the remainders are zero. However, we cannot simply decide that these points are stable. In fact, because of the existence of a zero eigenvalue, our first order perturbation analysis does not work and one has to use other methods, such as the center manifold theory, in order to reliably determine the stability of such a point. In the case of fixed points p 10 we have used the center manifold theorem and we write the results in the Appendix A.
This point corresponds to an expanding epoch in which the radiation density is zero and μ is constant. In other words, only matter and the scalar field G dominate the evolution. Using the corresponding ω eff one can easily show that a(t) ∝ t 6/13 , G(t) ∝ t −8/13 and φ 0 ∝ t −18/13 . In this phase the vector field mass (μ) is constant. It is noteworthy that the vector field's equation in MOG can be written as Therefore the vector field mass directly depends on the energy density. This situation is reminiscent of chameleon scalar fields where the scalar field's mass depends on the environment's density. It is well known that this property leads to a screening effect for hiding the scalar fields effect in the local experiments [31,32]. However, because of the appearance of φ 0 (t) in the denominator of (31), it is not trivial to claim that such a screening effect occur in MOG. However, screening effects are not just for scalar fields, and they can happen in theories with vector fields [33]. It is clear that this solution ( p 3 ) is completely different from the standard matter-dominated phase for which a(t) ∝ t 2/3 . Also, for p 3 eigenvalues are (− 13 6 , − 13 6 , −1, − 1 3 ). Therefore, surprisingly p 3 is an attractive/stable critical point. A simple interpretation is that, in the context of MOG and in the absence of , the universe can enter a permanent matter-dominated era. This is grossly inconsistent with the cosmological observations that imply that a matter-dominated epoch has been replaced by a stable dark energy-dominated era [34].
This critical point corresponds to an exact solution a(t) ∝ t 12/23 , G(t) ∝ t −10/23 , μ(t) ∝ t −18/23 and φ 0 is constant. Interestingly this is exactly a solution that has been obtained in [10] using the Noether symmetry approach. Also, p 4 is an unstable critical point since the eigenvalues are (−2.28, −2.04, −1, 0.26). From the stability point of view, it seems to be a true matter-dominated era. We recall that the radiation and matter eras are expected to be unstable critical points. However, p 4 is significantly different from the standard matter era for which the cosmic scale factor grows as a(t) ∝ t 2/3 . This situation is reminiscent of metric f (R) cosmology where there is a "φMDE" regime for which the cosmic scale factor does not follow the standard behavior [35,36]. We have the same fixed point ( f 3 with a slightly different ω eff ) even in the presence of the term. We shall discuss more about these fixed points in the next section.
p 5 : Radiation-dominated era The case p 5 corresponds to a standard radiation-dominated epoch whose ω eff = 1 3 . Therefore the cosmic scale factor grows as a(t) ∝ t 1/2 . In this era the vector field mass remains constant. Note that in GμMD era this mass increases with time. Also the eigenvalues are (−1, −1, 1, 0), which establishes an unstable radiation era. Therefore, although there is no standard matterdominated phase in MOG without ,there is a standard radiation era. It is also interesting that, unlike in the standard CDM model, the expansion rate in the radiationdominated era is larger than in the matter era ( p 3 ).
In this case, a(t) ∝ t 10/19 , G(t) ∝ t 2/19 and μ(t) ∝ t −10/19 . Therefore unlike the matter-dominated phases, G(t) increases with time in this era. The behavior of this radiation dominated era is slightly different from the standard case. Note that for this fixed point we have ω eff = 4 15 . p 7 : Strongly phantom attractor Eigenvalues for this point are (−18, −17, −9, −9). Therefore, p 7 shows a stable dark energy-dominated era. p 7 corresponds to a strongly phantom behavior with ω eff = −19 3 . It is easy to verify that the cosmic scale factor vary as a(t) ∝ (t rip − t) −1/8 . Here t is smaller than the constant t rip . In fact if t = t rip , the universe ends up with a finite-time. Although this is a stable dark energydominated epoch, the strong phantom crossing is inconsistent with the observation; see [37] for more details. Also the other fields grow as G(t) ∝ (t rip − t) −1/4 and μ(t) ∝ (t rip − t) 5/4 . It is evident that phantom crossing can occur in MOG; see Fig. 4. It is important to mention that the universe may not enter this phase. In fact, the phase-space trajectory of the system can end at p 3 before reaching the phantom attractor p 6 . In Fig. 1 we have explicitly shown this fact. For some different initial conditions in the deep radiation-dominated universe, we see that the dynamics reaches the stable p 3 point and stays there forever. On the other hand for slightly different initial conditions the fixed point p 6 is realized. This fact explicitly shows that the dynamics starts from an unstable point or close to it, i.e. p 5 . Of course there is no need to set the initial condition very close to that of solid lines in Fig. 1 in order to find p 6 . For example if one sets M to 3.7 × 10 −4 , p 6 is still realized by the system. More specifically, one may find initial conditions which lead to late time solutions p 3 or p 7 . We have shown the system's evolution for such a set of initial conditions in Fig. 1. Solid lines end at the strongly phantom attractor p 7 . Note that the density parameters are related to our dynamical system variables as Furthermore, in Fig. 2 we show the evolution of ω eff and ω DE . p 8,9 : Unstable unaccelerated era These critical points are on the lines p 1,2 . In the case of p 8,9 the effective equation of state parameter is − 1 3 and the eigenvalues are (−2(1 ∓ √ 10), 1, 0, 0), respectively. In this case, there is no acceleration and the cosmic scale factor grows uniformly with cosmic time. Also the other functions vary as G(t) ∝ t 2 and μ(t) ∝ t ± √ 10 . It is worthy to mention that such a behavior for the cosmic scale factor is impossible in the context of CDM model. In . We see that for the latter case the system falls in the stable GMD point which is grossly inconsistent with the current cosmological observations fact, in the standard model the accelerationä(t) can vanishes only at a single moment. However, in MOG the effect of the extra fields can effectively appear as "repulsive" force which can eliminate the attractive nature of the gravity and provide an unaccelerated expansion. It is also somewhat inconsistent with the expectation that MOG should lead to a stronger force in the weak field limit. We know that modified theories which try to address the flatness problem of the rotation curves of the spiral galaxies have to strengthen the gravitational force. However, it is necessary to stress that the evolution does not necessarily enter these epochs, i.e. P 8,9 for a cosmologically viable trajectory. p 10,11 : Unstable de Sitter-like era These points are also on the lines p 1,2 . The effective equation of state parameter for these points is ω eff = −1. Therefore p 10,11 correspond to an epoch in which the cosmic scale factor grows exponentially. The scalar fields vary with time as G(t) ∝ e 3t and μ(t) ∝ e ±t √ 39/2 . Also the eigenvalues are (−3 ∓ √ 78, −1, 0, 0). One may certainly conclude that p 11 is an unstable fixed point. On the other hand, the stability of non-hyperbolic point p 10 can be shown using center manifold theory. In the Appendix A we use this theory to specify the stability character of p 10 . Also in Fig. 5 we showed that p 10 is also an unstable critical point. Therefore, there is no late time stable de Sitter phase in the cosmic evolution of MOG, when is zero. Now let us summarize the general cosmological behavior of MOG in the absence of the cosmological constant (or equivalently when the scalar field G is not massive). In this  Fig. 2 The evolution of ω eff and ω DE for the initial conditions presented in Fig. 1 for the solid curves. Although the effective equation of state parameter vary smoothly with time and shows the strongly phantom behavior, the dark energy's equation of state parameter experiences singularities at r + y = 1 β case, it seems that MOG does not possess true consequences for the cosmological phases. The fixed points are: unstable radiation-dominated ( p 5 ) or unstable GμRD ( p 6 ), unstable GμMD point ( p 4 ) followed by the late time strongly phantom attractor p 6 . There is an interesting feature in the dynamics of this model. In fact there is a matter-dominated attractor GMD, i.e. the point p 3 . In other words, regarding the initial conditions, the universe can enter this matter-dominated phase and stay there forever. Although one can choose initial conditions for which the evolution does not realize p 3 , the late time attractor p 7 also is not physically accepted. In Fig.  1, we have plotted the evolution of the density parameters for two set of initial condition. In the absence of , MOG cannot be considered as a dark energy model in the sense that its extra fields cannot play the role of dark energy. Also, as we have already mentioned, the standard matter dominated era is replaced with the GμMD epoch at which the scale factor grows as a(t) ∝ t 12/23 . Clearly this behavior is far away from the standard case. We will discuss more about this important point in the next section.

MOG with cosmological constant ( = 0)
In this section we explore the original version of MOG, i.e. = 0. In this case, x = 0 and we use the same dynamical variables as introduced in Sect. 3.1 and find the relevant fixed points (y, r, m, z, x). Setting to zero the right hand side of Eqs. (23)- (27) we find the following critical points: In the following we shall study the stability of the above mentioned fixed points and discuss their physical interpretation.
We emphasize again that we are checking the possibility that if MOG can be a cosmologically viable theory. The cases f 1,2 indicate two different curves in the phase space. Every point on these curves is a fixed point and the relevant eigenvalues are the same as lines p 1,2 in the previous section. On the other hand, f 3 is similar to the GMD phase (where the matter and the scalar field G dominate the evolution i.e. p 3 ) introduced in the previous section. However, it is interesting that unlike the p 3 point, f 3 is an unstable fixed point. In other words, existence of the cosmological constant change the character of this fixed point. Eigenvalues of the stability matrix are (− 13 6 , − 13 6 , 13 6 , −1, − 1 3 ) and the cosmic scale factor grows as a(t) ∝ t 6/13 and for the other fields we have G(t) ∝ t −8/13 , μ is constant and φ 0 ∝ t −18/13 . We again are confronted with the question: does this point correspond to a suitable matter-dominated era? If not, then this fact puts serious doubt on the cosmological validity of this model even if it passes the local experiments and addresses the dark matter problem on a galactic scale. . Therefore f 4 is an unstable critical point. Furthermore, the scale factor grows as a(t) ∝ t 12/23 and for other functions we have μ(t) ∝ t −18/23 , φ 0 is constant and G(t) ∝ t −10/23 . It is interesting that there are two possible matter-dominated phases in MOG, i.e. f 3 and f 4 . The expansion rate in f 4 is slightly faster than f 3 but still slower than the standard matter-dominated case. We emphasize that this point is the main result of the current paper. The behavior of MOG in the matter-dominated era is crucial because this theory is an alternative theory for dark matter particles. Therefore, it should possesses an appropriate matter era in which structure formation happens without any need to cold dark matter particles.
As we mentioned before, this fact may put a serious doubt on the viability of this model. In fact a slower expansion rate, in principle, changes the duration of the matter-dominated phase. Consequently, there may be some impact on the cosmic microwave background (CMB) observations, for example on the angular size of the sound horizon. On the other hand, the strength of the gravitational force is different in MOG from the standard case. Therefore the growth rate of matter perturbations will be different from that of the CDM model. So observational data of galaxy clustering may help to distinguish the consequences of MOG. To summarize, existence of a non-standard expansion rate in the matter dominated epoch put a serious constraint on the viability of MOG but does not necessarily rule out MOG. One need to explore the cosmic structure formation in the context of MOG and carefully check the impacts of this theory on the CMB observations, in order to make a reliable decision about the viability of the theory [38].

f 5 : Radiation-dominated era
This point is a standard unstable radiation-dominated era for which ω eff = 1 3 . Also the relevant eigenvalues are: (2, −1, −1, 1, 0). We recall that even in the absence of there was an unstable radiation-dominated era p 5 . These points, namely p 5 and f 5 , are the only eras at which the scalar fields are constant. It is worthy to mention that at these epochs the vector field varies rapidly with time, as φ 0 ∝ t −2 , since it is directly coupled to the radiation density. However, after this phase the "running" of the scalar fields starts and eventually they dominate the evolution. Note that f 6 shows the same behavior as p 6 , i.e. an unstable GμRD era which is slightly different from the standard radiation era since ω eff = 4 15 . One can choose the initial condition in a way that the starting point of the evolution of the universe is either f 5 or f 6 . f 7 : Strongly phantom attractor In this case eigenvalues are (−18, −17, −9, −9, −8). Therefore, f 7 is a stable critical point. This point is exactly similar to p 7 presented in the previous section. Thus all our analysis for p 7 are also true for f 7 . We just mention that f 7 is different from a standard phantom crossing era and it cannot be considered as an acceptable late time solution, because of the strongly phantom crossing behavior. The initial conditions will specify the final fate of the system, i.e. f 12 or f 7 . It is interesting that MOG provides two kinds of late time accelerated expansions: de Sitter or a strongly phantom solution.
Furthermore, f 8,9 are the same as the unaccelerated solutions p 8,9 . In this case the eigenvalues are (−2 ∓ 2 √ 10, 1, 1, 0, 0), respectively, which show that f 8,9 are unstable points. For these solutions the scale factor uniformly increases with time. Our analysis in the previous section for p 8,9 is also true for these solutions. Also f 10,11 with eigenvalues (−3∓ √ 78, −1, 0, 0, 0) are similar to p 10,11 . It is clear that f 11 is unstable and f 10 is a non-hyperbolic fixed point. We tried to check the stability of f 10 using the center manifold theory as we did for p 10 . However, unlike p 10 , we realized that even the center manifold theorem does not reveal the character of these points. Therefore one needs to use more advanced methods such as the normal form theory. On the other hand we are considering a five-dimensional manifold and these methods become very complicated and are out of the scope of this paper. Thus we have to rely on our phase-space trajectories to decide about the stability of this point. We see in Fig. 5 that f 10 appears as an unstable critical point. This is also consistent with the result that we have already obtained for p 10 . Note that f 8,9 and f 10,11 are special cases of the lines f 1,2 with z = 2 and z = 3, respectively. , which clearly show that these points are unstable. Obviously these points cannot be considered as early time unstable fixed points (i.e. points which can show an inflationary period). The contribution of the cosmological constant is comparable to the other components while we know that does not play an important role in the early universe. Also one may expect for a very early time unstable fixed point all the eigenvalues should be positive. It is also clear in

Phase-space analysis at infinity
It is important to stress that from Eq. (12) one may straightforwardly conclude that H , in principle, can become zero. In such a moment, our dynamical variables become infinite. Therefore it seems necessary to check this possibility. In other this case we find (y r + r r ) ≤ 0.
Therefore our fixed points should satisfy the relevant constraint. Now it is straightforward to verify that the only allowed fixed points at infinity are q 3,4 and q 7 -q 10 . The corresponding eigenvalues of the stability matrix constructed from Eq. (34) for q 3,4 and q 7 -q 10 are ( 1 10 , − 1 20 , − 1 20 , 0, 0) and ( 1207 570 , − 1207 1140 , − 1207 1140 , 0, 0), respectively. Note that due to the symmetry of the equations, the eigenvalues of stability matrix for q 7 to q 10 are the same. This means that accepted fixed points are unstable critical points.
One may provide a simple interpretation for these points. In fact they are "middle" time, and not an early or late time point, where the scalar fields G and μ dominate the evolution and the expansion rate becomes zero for a moment. In other words if q 3,4 or q 7 -q 10 were stable points, then the universe could enter a static phase and stay there forever. In Fig. 5 we have shown two different trajectories which realize this point. More specifically the blue dotted curve starts from the radiation dominated point f 5 and after passing the infinity fixed points q 9 and q 4 , falls into the late time point f 7 . The green solid line shows another path which starts from GμRD point, f 6 , passes the unstable infinity fixed points q 7 and q 3 , and eventually reaches the strongly phantom stable point f 7 .

Conclusion
In this paper we have considered the cosmological behavior and consequences of a scalar-vector-tensor theory of gravity, known as MOG in the literature. Although this theory is known as an alternative theory for dark matter particles, we have investigated its viability as a dark energy model using the so-called dynamical system method. In fact, we checked the possibility that the extra fields of MOG can play the role of dark energy. We first derived the autonomous equations of the relevant dynamical system, i.e. Eqs. (23)- (27), and we found the corresponding fixed points in two different cases: = 0 and = 0. In Sect. 3.1, we showed that in the absence of , there is not a standard late time epoch. More specifically, the evolution starts from the unstable radiation-dominated epoch p 5 , then reaches the unstable matter-dominated epoch p 4 and eventually ends at the stable late time strongly phantom attractor p 7 , which is not physically acceptable. It is worth mentioning that depending on initial conditions, the universe could enter the stable GMD era ( p 3 ). Obviously such an initial condition leads to a wrong cosmological behavior. Although p 5 is a standard radiation-dominated era, the matter-dominated epoch p 4 and the late time attractor p 7 are not standard. In fact, in p 4 the scale factor grows as a(t) ∝ t 12 23 , that is, slower than the standard case in which a(t) ∝ t 2 3 . On the other hand, for the late time attractor p 7 , ω eff is − 19 3 , which is not a standard late time solution. It shows that in the absence of , one cannot recover standard cosmological epochs in MOG.
In Sect. 3.2, we investigate the original form of the theory, i.e. with nonzero . It is noteworthy that in this theory the cosmological constant can be considered as the mass of scalar field G (or equivalently χ ); see Eq. (1). Also note that there are two coupling constants ω 0 and κ in this theory. Surprisingly these parameters do not appear in our phasespace analysis. In other words, our fixed points are numbers and there is no free parameter to be constrained. In this sense MOG behaves like the CDM model where the fixed points do not include any free parameter.
The first effect of the nonzero cosmological constant is to change the character of the stable matter-dominated phase p 3 (or equivalently f 3 ). In other words, two matterdominated epochs of MOG, namely f 3 and f 4 , are now unstable as expected. Interestingly, there are also two early time radiation-dominated f 5 and f 6 and also two late time attractors f 4 and f 7 . In fact, depending on the initial conditions, universe could start from standard radiation era f 5 or GμRD epoch, f 6 , which is slightly different from the standard radiation-dominated era, continues toward matterdominated epochs f 3 or f 4 and end at the strongly phantom, non-physical, attractor f 7 or the standard stable de Sitter epoch f 4 . Also, the phantom crossing behavior in MOG is clearly seen, as shown in Fig. 4. As we already mentioned, a true cosmological path could start from a standard radiation-dominated epoch f 5 , continue toward one of the unstable matter-dominated GMD ( f 3 ) or GμMD ( f 4 ), and finally reach the stable late time solution f 4 as shown in Fig. 3.
Furthermore, there is a new feature in the cosmology of MOG that is absent in the standard model. In fact there are unstable eras p 8,9 (or equivalently f 8,9 ) in which the universe expands uniformly. In these eras the extra fields, on large scales, behave like an effective repulsive "force" and cancel out the attractive nature of the gravitation, and consequently cosmic acceleration vanishes (Figs. 1, 4).
As a final remark, we emphasize that there is no standard matter dominated phase in the model of MOG studied in this paper. In fact with or without cosmological constant the scale factor grows as a(t) ∼ t 0.5 instead of the conventional t 2/3 . In principle, this may lead to inconsistencies with CMB and large scale structure formation observations. Therefore, it seems necessary to investigate this issue with more careful consideration. We leave this issue for future study. It is important to mention that even if the above mentioned observations rule out the original form of MOG, one may find some special models of MOG (with different self-interaction potentials) which can pass the observations. Acknowledgments This work is supported by Ferdowsi University of Mashhad under Grant No. 39640(04/11/1394). We thank John Moffat, Viktor Toth, Martin Green and especially Fatimah Shojai for insightful comments. Sara Jamali thanks Luca Amendola and Sohrab Rahvar for valuable discussions.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .
Appendix A: Stability of p 10 using center manifold theory The stability of a hyperbolic point x 0 of the non-linear systemẋ = f (x) can be determined by the behavior of the linear systemẋ = Ax, where A = Df(x 0 ) and D denotes differentiation with respect to x. On the other hand, the stability of a non-hyperbolic point can be determined by recognizing the behavior of the center manifold near that point [26]. In other words, since zero eigenvalues reveal no information as regards the qualitative and stability behavior of the system, one should find a way to get information as regards this part of the system.
Considering the phase space of MOG, we find the nonhyperbolic fixed point p 9 . For the sake of simplicity, we change the coordinate in a way that p 9 lies at the origin. More specifically, the simple change of coordinates as z → z + 3 and m → m + 13 8 would transfer p 9 to the origin. In this case the transformed autonomous equations are  Here F = F(z, m, h(z, m)) and K = K (z, m, h(z, m)). Solving this matrix equation, one finds the expansion coefficient, a, b, . . . and the behavior of center manifold using the following equation: X = CX + F(X, h(X)).
As already mentioned, X stands for the zero components. It is noteworthy that in the expansion of y and r with respect to z and m, we look for the first nonzero coefficient, which will determine the behavior of the center manifold. Using this method for p 10 , one finds m ≈ 1.63228m 2 z 2 + 1.17611m 2 z + 0.406759m 2 − 5.08273mz 2 − 3.39935mz, z ≈ −0.342063m 2 z 2 − 0.357055m 2 z +0.138952m 2 + 1.49451mz 2 + 2.03504mz + m.
Using these equations, in Fig. 6, we have plotted a twodimensional phase space in the z-m plane. One may straightforwardly conclude that the origin, p 10 , represents an unstable critical point.