The phase space analysis of modified gravity (MOG)

We investigate the cosmological consequences of a scalar-vector-tensor theory of gravity known as MOG. In MOG, in addition to metric tensor, there are two scalar fields $G(x)$ and $\mu(x)$, and one vector field $\phi_{\alpha}(x)$. 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 $\Lambda$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, strongly phantom crossing happens. Also we find that this theory in its original form ($\Lambda\neq 0$), possesses a true sequence of cosmological epochs. Albeit we show that, surprisingly, there are two radiation dominated epochs $f_5$ and $f_6$, two matter dominated phases $f_3$ and $f_4$, and two late time accelerated eras $f_{12}$ and $f_{7}$. 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)\sim t^{0.46}$ and $t^{0.52}$, respectively, which are slower than the standard case, i.e. $a(t)\sim t^{2/3}$. 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 consequently the weak equivalence principle, in principle, can be violated. Naturally, the free parameters of MOG are chosen such that 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 somehow 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 is recently 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 ina e-mail: sara.jamali@stu.um.ac.ir b e-mail: mroshan@um.ac.ir vestigated 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 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 that how is the cosmological behavior of MOG considering the variety of fields that have been incorporated. On the other hand there are a few papers considering cosmology of MOG. For example in [10] the Noether symmetries of the cosmic pint-like Lagrangian of MOG has 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] for relevant works.
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 provide a fast and reliable procedure to numerically solve the field equations. Note that field equations of MOG are drastically complicated than Einstein's general relativity (GR), see equations (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 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 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 [13].
This paper is organized as follows. In section 2 we briefly introduce MOG and the modified Friedmann equations. In section 3.1 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 section 3.2 we bring back the cosmological constant and analyze the system. In section 4 we study the phase space of the system at infinity. Finally, conclusions are drawn in sec 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 Λ as 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 anti-symmetric tensor reminiscent of the Maxwell's 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. Albeit 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 Where G µν is the Einstein tensor and J α is a "fifth force" matter current defined as 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's 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 the ordinary matter, and 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 T µν = (ρ + p)u µ u ν + pg µν Finally, bearing in mind that χ 2 = 2/G and ψ = ln µ, we find the following Friedmann equationṡ Where 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 mentioning that scalar fields µ and G have negative contributions to the total energy density. More specifically the kinetic termsμ 2 µ 2 anḋ G 2 G 2 appear with negative sing 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 equations (7)-(11) as followṡ Note that using the field equation of the vector field we have replaced φ 0 , the only non-zero 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. equations (12)- (15), are drastically complicated than the standard Friedmann equations. However, as we shall see, despite of 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 [14]. 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 [14] for more details. In this case one may straightforwardly set the matter current as J = κρ m . Where κ is an another positive coupling constant. It is noteworthy that in a non-homogeneous spacetime because of the coupling between matter and the vector field, these relations are not true and one may expect significant departures from Λ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 [15] and [16]. Therefore, regarding the energy exchange between matter and the vector field in a non-homogeneous background, it seems interesting to check the possibility that if MOG can explain the accelerated expansion without invoking the cosmological constant and just by taking into account the matter inhomogeneities. This issue can be a matter of study for future works. Therefore, in what follows 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 equations (12) and (13) as where G N is the Newtonian gravitational constant and H =ȧ/a is the Hubble parameter and a dot denotes derivative with respect to cosmic time t. Where ρ 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 = pDE ρDE . It is also useful to write the effective equation of state parameter ω eff , that 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 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 prime denotes 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, type and the stability of each fixed point, can be found using the eigenvalues of the stability matrix [17]. Now let us define the following dimensionless variables Substituting these dynamical variables into equation (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 [18] and references therein. For more recent works we refer the reader to [19]. There is also a non-minimally coupled scalar field χ in MOG. However, in the following we explicitly show that there is no transitions from anti-gravity to gravity in the context of MOG. More specifically, if the evolution starts from a 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 form 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 equation (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 that the Λ term dominates and soĠ = 0 andG > 0, see also equation (14). However one may naturally expect that Λ becomes important only at the late times. Therefore it is very unlikely to 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 β = GN 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 set x = 0 in the equations of motion. In this case, the critical points (y, r, m, z) of equations (23)-(26) are listed ) ω eff = 5 18 Surprisingly the fixed points are numbers and there is no free parameter to be constrained. This point is also clear from equations (23)-(27), where with the aid 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 coordinate of the fixed points [13]. 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 a chance to make the model more consistent with the cosmological observations. However this is not the case in MOG.
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 , 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, sign of the eigenvalues of p 1 are negative and reminders are zero. However, we can not 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.
p 3 : G-Matter dominated (GMD) era: 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 [20]. However, because of the appearance of φ 0 (t) in the denominator of (31), it is not trivial to claim that such screening effect occur in MOG. Albeit screening effects are not just for scalar fields, and can happen in theories with vector fields [21]. 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 matter dominated epoch has been replaced by a stable dark energy dominated era [22].
p 4 : Gµ-Matter dominated (GµMD) era: 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 [23] and [24]. We have a 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: 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 matter dominated 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 radiation dominated era is larger than in the matter era (p 3 ).
p 6 : Gµ-Radiation dominated (GµRD) era: This point corresponds to an unstable Gµ-Radiation epoch for which the eigenvalues are (− 9 5 , 1, − 9 10 , 9 10 ). 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 . Where 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 energy dominated epoch, the strong phantom crossing is inconsistent with the observation, see [25] for more details. Also 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 . We see that for the later case system falls in the stable GMD point which is grossly inconsistent with the current cosmological observations.
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 at the deep radiation dominated universe, we see that the dynamics reaches the stable p 3 point and stay there forever. On the other hand for a slightly different initial conditions the fixed point p 6 is realized. This fact explicitly shows that the dynamics starts from or close to an unstable point, i.e. p 5 . Of course it is not needed 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 set Ω 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 ends at strongly phantom attractor p 7 . Note that the density parameters are related to our dynamical system variables as Furthermore, in Fig. 2 we have shown 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 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 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 a unaccelerated expansion. It is also somehow inconsistent with the expectation that MOG should lead to stronger force in he weak filed 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. Albeit, 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 case, it seems that MOG does not possesses a true consequences of 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. 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 can not be considered as a dark energy model in the sense that its extra fields can not 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 introduced in section 3.1 and find the relevant fixed points (y, r, m, z, x). Setting to zero the right hand side of equations (23)-(27) we find the following critical points: ) In what follows 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. 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 confront again the question that does this point correspond to a suitable matter dominated era? If not, then this fact put a serious doubt on the cosmological validity of this model even if it pass the local experiments and address the dark matter problem in 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 is 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 too 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 impacts 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 than the standard case. Therefore the growth rate of matter perturbations will be different from that of ΛCDM model. So observational data of galaxy clustering may help to distinguish the consequences of MOG. To summarize, existence of a nonstandard 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 [26].
f 5 : Radiation-dominated era: This point is an 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 Solid lines reaches the de-sitter like attractor f12 through the f3 while dashed lines ends at strongly phantom attractor f7 passing f4. A common feature in the dynamics of these different initial conditions is that the scalar fields remain constant during the early stages of the universe.
fields starts and eventually they dominate the evolution. Note that f 6 also have 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 initial condition in a way that starting point of the evolution of the universe be 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 can not 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 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 center manifold theorem does not reveal the character of these point. 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 to 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.
f 12 : de Sitter-like attractor: This point corresponds to the stable dark energy dominated universe where the dynamics is dominated with the cosmological constant and the scalar field G. Eigenvalues are (− 16 3 , − 14 3 , − 13 3 , − 13 3 , − 13 3 ) and the cosmic scale factor grows exponentially. Also this point corresponds to the exact solution G(t) ∝ e −4t/3 and µ=constant. We recall that in the absence of the cosmological constant there is no stable dark energy dominated era with ω eff = −1.
f 13,14 : GµΛ era: These points correspond to an epoch at which the scalar fields together with the cosmological constant dominates the evolution and the ordinary matter and radiation energy densities are zero. The cosmic scale factor grows exponentially and for other functions we have G(t) ∝ e −2t/7 and µ(t) ∝ e −22t/7 . Also the relevant eigenvalues are (−6.03, − 30 7 , − 23 7 , − 23 7 , 2.74) which clearly show that these points are unstable. Obviously these points can not be considered as early time unstable fixed points (i.e. points which can show an inflationary period). Because 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 that for a very early time unstable fixed point all the eigenvalues should be positive. It is also clear in The long-dashed line shows that universe starting from radiation dominated f6 and passing through unstable points f7,14 and f8, f10, toward unstable era f13 and ends at a late time attractor f12. The green solid line also starts from f6, passing through q7 and q3, ends at strongly phantom f7. Note that there are unstable fixed points at infinity, q3,4 and q7, 9. toward an unstable matter dominated epoch, f 3 or f 4 and finish in a stable dark energy dominated point, f 12 . Note that strongly phantom attractor, f 7 is not acceptable as a standard late time solution. It is clear form Fig. 5 that for specific initial conditions phase trajectory of the system passes through unaccelerated eras f 8 and f 10 . However this trajectory is not acceptable because it does not realize the matter dominated points. Also we have plotted the behavior of ω DE and ω eff in Fig. 4. It is important to mention that although ω eff remains always greater than −1, ω DE becomes smaller than −1 about present time. In other words, MOG can provide a slightly phantom solutions at present time.
Let us summarize the main results of this section. We explicitly showed that MOG possesses true cosmological sequence of the cosmological epochs. In fact there is a standard radiation dominated point f 5 . Also there are two unstable matter dominated phases which are not standard in the sense that the cosmic scale factor grows slowly than the standard case. On the other hand there is a stable late time solutions f 12 .

Phase space analysis at infinity
It is important to stress that, form equation (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 words, we have to check the fixed points lying in the infinity of the dynamical system. As we shall show, such a study will ensure us that there is nothing special, such as a regular bounce, at the moment when the expansion rate of the universe become zero. We know that even if we find some fixed points in the infinity, it does not mean that a true cosmological trajectory will realize that points. However, finding these points and analyzing their stability will provide a better understand-ing for the general behavior of the system. This section Also if x r = m r = z r = 0 (namely for q 1 and q 2 ), one may divide (35) by √ 1 + R 2 and take the limit of R → ∞. In this case we find (y r + r r ) ≤ 0 (37) 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−10 . The corresponding eigenvalues of the stability matrix constructed from equations (34) for q 3,4 and q 7−10 are ( 1 10 , − 1 20 , − 1 20 , 0, 0) and ( 1207 570 , − 1207 1140 , − 1207 1140 , 0, 0), respectively. Note that due to symmetry of equations, the eigenvalues of stability matrix for q 7 to q 10 are the same. This means that accepted fixed points are unstable critical point.
One may provide a simple interpretation for these point. In fact they are "middle" time, and not an early or late time, point where the scalar fields G and µ dominates the evolution and the expansion rate becomes zero for a moment. In other words if q 3,4 or q 7−10 was a stable point, then the universe could enter an 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 it's viability as a dark energy model using the so-called dynamical system method. In fact, we checked the possibility that if 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. equations (23)-(27), and found the corresponding fixed points in two different cases: Λ = 0 and Λ = 0.
In section 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, 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 can not recover standard cosmological epochs in MOG.
In section 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 equation (1). Also note that there are two coupling constants ω 0 and κ in this theory. Surprisingly these parameters do not appear in our phase space analysis. In other words, our fixed points are numbers and there is no free parameter to be constrained. In this sense MOG behaves like Λ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 matter dominated epochs of MOG, namely f 3 and f 4 are now unstable as expected. Interestingly, there are also two early time radiation dominated and 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 matter dominated epochs f 3 or f 4 and end at 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 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 existence of the extra fields, in large scales, behave like an effective repulsive "force" and cancel out the attractive nature of the gravitation, and consequently cosmic acceleration vanishes.
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 considerations. We leave this issue as a 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. A stability of p 10 using Center Manifold Theory

Acknowledgments
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 [17]. In other words, since zero eigenvalues reveals no information about the qualitative and stability behavior of the system, one should find a way to get information about 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.