Cosmology of the selfaccelerating third order Galileon

In this paper we start from the original formulation of the galileon model with the original choice for couplings to gravity. Within this framework we find that there is still a subset of possible Lagrangians that give selfaccelerating solutions with stable spherically symmetric solutions. This is a certain constrained subset of the third order galileon which has not been explored before. We develop and explore the background cosmological evolution of this model drawing intuition from other even more restricted galileon models. The numerical results confirm the presence of selfacceleration, but also reveals a possible instability with respect to galileon perturbations.


Introduction
Supernova data [1,2] in conjunction with data from the CMB [3] suggests that the Universe is expanding at an accelerated rate. Though alternatives to accelerated expansion have been proposed in terms of inhomogeneous models (see for instance [4]), it is now most commonly asked in cosmology not whether, but how or why the Universe's expansion is accelerated.
The simplest way to get accelerated expansion is by adding a cosmological constant. This solution is now part of the standard cosmological model and it fits well with modern precision data of the cosmic microwave background, supernova data and baryon acoustic oscillations [3].
However simple and elegant the cosmological constant solution may seem, questions still remain. If this really is the right explanation for the cosmic acceleration, what is the physical origin of the constant? Why has it got the value that it has? A possible physical source is vacuum energy, but attempts to predict its value from known microscopic physics notoriously yield results more than a hundred orders of magnitude larger than the observed values [5]. Though there are proposed resolutions to this, appealing both to our lack of knowledge of fundamental physics and to anthropic arguments [5], one might alternatively look for a dynamical solution important on cosmological scales.
The alternatives to the cosmological constant come in roughly two categories. The first proposes some sort of generalised dark energy [6]. The other possibility is that of some kind of modification of gravity on large scales [7]. In most cases models of dark energy and modified gravity are dual in the sense that a conformal transformation can take a theory of one type into a theory of the other kind. The main difference lies in what the formulation is like in the frame where one assumes the background metric. If the theory in this case is minimally coupled to gravity, the theory is considered a dark energy model, whereas if it is nonminimally coupled to gravity in this frame it is considered a modified gravity model.
Recently the Galileon model [21] is a modified gravity that has attracted a good deal of attention, though it has also been formulated as a pure dark energy candidate. Motivated by the DGP and its Vainshtein mechanism, the authors of [21] constructed the Galileon model, a model of a scalar field π that has a symmetry of the equations of motion under the galilean shift: and found all five possible Lagrangians that obey this symmetry in four dimensions. The five different terms are denoted by the order of the scalar field π they contain, and there is one term unique up to total derivatives of each of the orders one to five, and none of higher than fifth order. They considered a theory built from linear combinations of these terms and found that it could yield selfaccelerating solutions in a cosmological background. By considering spherically symmetric solutions and their stability, bounds were put on the coefficients. After this promising start considerations of the cosmological evolutions of galileon theories or subsets of such theories were considered in [22], which dealt with only the behaviour of the third order galileon term, which turned out to give a slight degravitation of the cosmological constant, but gave no viable selfaccelerating solutions, and for slightly different couplings to gravity than those given in [21] in [23,24,25] which all explored the cosmology of theories of up to fourth and in the case of [24] fifth order galileon models.
Several other theoretical aspects of the galileon theory have also been considered, such as the null-energy-condition (NEC) violation [26], covariantisation [27], the Vainshtein mechanism and its effect on spherically symmetric solutions [28] and black hole accretion [29].
From theoretical considerations of the symmetries, many generalisations have been built. In [30] expressions for all the galileon invariant Lagrangians in arbitrary dimensions were derived. Generalisations to arbitrary p-form galileons were obtained in [31] and the multi-galileons, which are multiplets of scalar galileons, contained in this formalism has since been considered in [32] and [33] which mainly focused on the two-or bi-galileon case, in [34] which focused on theoretical aspects and the evasion of Derrick's theorem [35,36] and [37] which considered instabilities in spherically symmetric solutions. The paper [38] showed that the N-galileons are the goldstone bosons of the spontaneous symmetry breaking of the N + D dimensional diffeomorphism group in forming a codimension N brane of D dimensions. This property of the multigalileons was also mentioned in [32] where it was noted that the bi-galileons were the goldstones from cascading gravity from 6 dimensions.
The fact that the galileon is a natural artifact of a five dimensional theory was further reviewed in [39] which found in this a unification of the DBI [40,41,42,43,44] and galileon theories. The two theories turn out to be different nonlinear realisations of the 5D diffeomorphisms on a 4D brane. A generalisation of this theory was formulated and explored in [45].
The connection between massive gravity and higher dimensional theories are well known, and in the case of the galileon the connection to massive spin 2 fields was shown in [46].
Many other generalisations of the galileon model have also been explored. In [47] a generalisation where π → f (π) was considered in a cosmological context. Some of these generalisations break the galilean invariance. For instance a generalisation to the coefficients of the different galileon Lagrangian terms being functions of the galileon field was considered for up to third order galileons in [48], [49], [50] and [19]. Up to third order galileons have also been generalised to theories with coefficients not only dependent on the galileon field itself, but also on its kinetic term X = −g µν ∂ µ π∂ ν π/2 in [51] which considered perturbation equations in such a generalised scenario, in [52] where cosmological evolution was reviewed and in [53] which considered inflation.
Though the original formulation of the galileon was proposed as a dark energy alternative [21], the model has recently been studied as an alternative to inflation. In particular its NEC-violating properties have attracted attention in this context as this can give rise to a galileon genesis from flat space [54] yielding a theory with no Big Bang. The NECviolation was also found in the generalised theory of [53] in the inflationary scenario. In [55] a covariant version of the galileon is used to drive inflation, motivated by the fact that the constant shift symmetry is already only mildly broken in the traditional slow-roll inflation, and hence the galileon is a natural and useful place to look for generalised inflation.
We want to take the model back to its basic formulation in [21], where we feel that a simple self-accelerating solution has yet to be properly explored. In this paper we find from the theoretical constraints in [21] and [26] that a third order galileon theory with a first order term, a so-called tadpole, seems to be allowed, stable and yielding selfaccelerating solutions in a certain regime of the parameter space. We explore this simplest of the galileon selfaccelerating solution and its background evolution, as we feel that a thorough treatment of this theory seems the natural place to go first from the results of the original papers [21] and [26]. A treatment of the perturbations for this model will be postponed for a future publication.

The galileon Lagrangian
Before we start the treatment of our particular galileon model, we would like to recapitulate the formulation of the galileon model given in [21].
Looking at the action: where M P L is the four-dimensional reduced Planck mass,R is the Einstein frame Ricci scalar, T µν is the stress-energy tensor for the matter fields, 1 2ĥ µν is the coupling between gravity and the matter fields compensating for the matter not feeling the same metric as the Einstein frame one and π is the galileon. The equations of motion of the galileon Lagrangian is invariant under the galilean transformation In 4D there are 5 galilean invariant Lagrangians, one for each order of the galileon up to 5. Higher order terms involve only total derivative terms. Defining Π ≡ ∂ µ ∂ ν π, [. . .] to be the trace of the operator . . . and · to denote Lorentz-invariant contraction of indices, the five galilean invariant Lagrangians are: In general, then, the full galilean field Lagrangian takes the form For our selfaccelerating solution we want a stable deSitter attractor for a vanishing matter stress-energy tensor. Considering in this case the effect of the galileon on the cosmological evolution yields a configuration for the π in terms of the wanted asymptotic value H 0 of the Hubble parameter in this case. Plugging this into the equation of motion for the galileon leads to the following condition on the coefficients of the Lagrangian 1 Now, we want the internal galilean symmetries and external Lorentz symmetries to be unbroken, so the perturbations should obey these symmetries, hence if our equations of motion for the galilean terms are given by then the perturbations to the equations of motion must be given by where d i are the coefficients for the perturbations. To simplify the calculations we assume the coefficients d i to incorporate the zeroth order deSitter solution plus the perturbations. We can then find invertible relations between the d i s and the c i s for i = 2 . . . 5.
c 1 is then simply given by the equation (2.10). The inverse matrix becomes We should now check what constraints must be put on these coefficients for the solution to admit a stable deSitter solution with spherical symmetric solutions about compact sources. These configurations should be ghost free. We also want the fluctuations to be subluminal. Combining all this demands the conditions for the d i s are given by We should study the cosmology of such a model where the null-energy condition is viable, and see whether it can reproduce results for ΛCDM cosmology in terms of background equations and perturbations.
Actually the constraints above already ensure subluminality and no breaking of the NEC. In the coming subsection we check that this is actually true in the case considered in the article [26].

The case of L 1 = 0
In the article [26] they consider a deSitter solution with c 1 = 0 and check for the demand of no violation of the NEC conditions means that: When c 1 = 0, then c 5 is completely determined from the deSitter criterion (2.10). This means that we can rewrite our transformation matrices so that they only transform between indices 2, 3, 4. The results are By using the transformation matrices above to turn the condition on cs given in equation (2.14) into a condition on the ds, we realise that the conditions given in equation (2.13) are sufficient to prohibit the breaking of the NEC.

A possible Lagrangian for the galilean modified gravity with selfaccelerated solutions
Since the demands on the coefficients already ensures a deSitter solution, we should look at the early evolution and the evolution of perturbations in this theory. Since the assumptions also imply stability of the deSitter solution we should only check the stability in the more general non-deSitter regime. However, constraints from spherical sources compels us to choose d 4 = d 5 + c 4 = c 5 = 0. 2 Hence we are left only with c 2 and c 3 . The constraint in equation (2.13) becomes and requires c 1 to be larger than zero for the constraints on the ds to hold. Hence we get a Lagrangian: where H 0 is the value of the Hubble rate in the deSitter attractor solution. In the later numerical computations we set this value to the present day Hubble rate. We choose to work in a Jordan frame where matter is minimally coupled to the galileon, thus considering the galileon a modification to gravity as was eventually also done in the formulation made by [21]. Our final action becomes where S matter is the action of other fields like photons, neutrinos, baryons and dark matter. Because of the coupling between the galileon and gravity the resulting Einstein equations and equation of motion for the galileon field are no longer galilean invariant in a non-flat background. However, since the Jordan frame has been chosen a coupling between gravity and the galileon must be supplied. This particular choice makes the action reduce to the ordinary Einstein-Hilbert action in the limit of no galileon fields, whilst keeping the galilean symmetry in flat space. Hence this choice for a coupling to gravity is natural. It is also the same non-minimal coupling as chosen in [21] so all their stability and selfacceleration conditions hold for this model.
From this and reviewing the result of [22] we can find the equations of motion. For π it is: Using the result from [22] combined with the derivation in [56] the Einstein equations become This should enable us to study background cosmology and eventually also perturbation theory in this modified version of gravity.

Equations of motion in a homogeneous and isotropic Universe
To simplify calculations and the considerations of the equations of motion derived above we start by assuming isotropy and homogeneity. We also assume the Universe to be spatially flat. This means that we are working in a flat FRW metric background. This derivation is very similar to the one presented in section 4 of [22].
In the flat FRW metric background, the equation of motion for π in equation (3.5) now becomes: where we see that the term proportional to c 3 is equivalent to the term found in [22]. The 00 component of the Einstein equation gives the equivalent of the modified Friedmann equation The combination − 1 6 00 + 3 a 2 ii yields the Raychaudhuri equation which has a slight discrepancy of e −2π HπM 2 PL with respect to the one found in [22] which seems to be due to a difference in calculations of covariant derivatives.
As in [22] we can find an effective energy density ρ π and pressure P π for the π field.
Again the part of P π resulting from L 3 is not exactly the same as found in [22]. The same discrepancy, presumably caused by the same calculations difference, is still present.

Numerical simulation of the background equations
Before we start the numerical simulations of the background equations (3.7), (3.8) and (3.9) it is useful to rewrite the derivatives with respect to t in terms of a new variable x = ln a. Since d dt = H d dx we get for the equation of motion of the galileon: where ′ denotes a derivative with respect to x. The Friedmann and Raychaudhuri equations come out as and We now take advantage of the nice property of the Jordan frame that matter is not coupled to the galileon field. This means that the equations that determine the development of the matter or radiation from their equations of state are not altered, so ρ m = Ω m a −3 = Ω m exp (−3x) and ρ γ = Ω γ a −4 = Ω γ exp (−4x), ρ = ρ m + ρ γ and ρ + 3P = ρ m + 2ρ γ .
We use the above equations to make a set of three coupled differential equations to solve numerically in H, π and π ′ . First we use the equation: thereby transforming the set of second order equations (4.1-4.3) in π into a set of first order equations in π, π ′ and H. Then we need an equation for π ′′ and one for H ′ . However, we need at least one of these not to depend on the other. We enter the result from the first Friedmann equation into the Raychaudhuri equation, to simplify slightly. The result is: Then using (4.5) to solve for the second derivative of the galileon field we get: Plugging this into equation (4.1) we get that For reference it is nice to have the formulas for the galileon density and pressure given in terms of the ′ derivatives: −2M 2 PL He −2π H ′ π ′ + Hπ ′′ + 4M 2 PL H 2 π ′ 2 e −2π (4.10)

Results in different particular regimes
Before moving on to our full model it is interesting to review and study the behaviour in different regimes of special values for the coefficients. Some of these may be interesting for their simplicity. Others may be comparable to theories considered elsewhere. Together they provide useful intuition for our interpretation of the evolution of the full model.

All coefficients zero
If we set c 1 = c 2 = c 3 = 0, the galileon field will still be present and dynamical, and will have a special effect on the evolution of the Universe. We can see this directly from writing out the equations in this case: The galileon equation of motion gives: Adding the Friedmann equation and Einstein equations together and writing ρ tot = ρ + ρ π , P tot = P + P π , we get:ä which leads to a demand for an effective total equation of state at all times of This means that the Universe appears radiation dominated at all times. Numerical results show this, as can be seen from the evolution of H shown in figure 1.  The case of c 1 being the only non-zero constant greatly resembles the case considered above in section 5.1 with no coefficients. However we now have a linear potential making the theory start off in such a way that the evolution starts not like a radiation dominated one, but one defined by the constant c 1 if the galileon starts out very small. We see this from the equation of state of the galileon in this case: As the galileon grows the term on the right will grow leading eventually to accelerated expansion. If the galileon field on the other hand has negative values, the evolution will move towards the radiation like one as the absolute value of the field becomes larger. Adding together the Friedmann and Raychaudhuri equations we get that We note that as c 1 = 0 the theory can only approach a radiation dominated like epoch asymptotically as ρ tot gets very large.
In the tadpole only case the equations for the numerical simulations become: and H ′ is simply given by equation ( From the two plots we see that in the beginning the Universe behaves as if radiation dominated, or very nearly so. This is because early on ρ tot is very large so that we approach the asymptotic value for the equation of state parameter w eff = 1/3. As the evolution goes on it deviates from this and starts showing a smaller equation of state parameter and a smaller deceleration. This shift happens earlier, giving smaller values for the equation of state parameter and slower deceleration the larger the value of c 1 . For sufficiently large values of c 1 a period of acceleration can even be reached. However, after experiencing a minimum value of the equation of state parameter or deceleration, the equation of state parameter starts growing again eventually asymptotically towards an evolution with 0 < w eff < 1/3. We see that by making c 1 large enough we can get a value for H today equal to the one we find from current observations 3 . However, in these cases the accelerated expansion has already taken place and passed, the Universe now being in a decelerated state. We also see that most of the evolution history will change considerably in this case. Considering  The case of only c 2 different from zero is in fact a type of Brans-Dicke theory [57]. We can see this from starting with the generic Brans-Dicke action: To get this into a form similar to our own we realise that we must set: (5.8) and inserting this we see that for the second term to correspond to our c 2 term we must have: ω(φ) = c 2 8φ (5.9) I.e. a theory with only the c 2 term present is equivalent to a Brans-Dicke theory with the Brans-Dicke parameter ω(φ) ∝ 1 φ .   The background equations we need to simulate now have:

Plugging this into equation (4.1) we get that
The results for the evolution of the Hubble parameter H and the effective equation of state parameter for different values of c 2 is shown in figures 8 and 9 respectively.   We see that c 2 in fact acts as a demarcation scale between the non-galileon background theory and the theory with a galileon with all the coefficients equal to zero. The larger the value of c 2 , the later the galileon will have any impact on the evolution. However, as we raise the initial values of the galileon, the larger is the c 2 we need to compensate for this. Hence quite large initial values for the galileon will make the galileon dominate and give a Universe at early times with much smaller values of the Hubble parameter before transitioning into the radiation like phase. For instance in the case of c 2 = 50M 2 PL the initial values of π and π ′ were set to 0.5e − 10 and 0.5e − 50 respectively in the case shown in figures 8 and 9. Changing the values to 0.5e−10 and 0.5e−10 respectively had very little effect, however a change to 0.5e − 5 and 0.5e − 5 for the initial values respectively changed the evolution of the background quite considerably, resulting in a much lower value for the Hubble constant, reaching a value comparable to that in the plotted case for c 2 = 0.01M 2 PL only very recently (redshift of only a small fraction), and from there transitioning into radiation like behaviour.
Since the theory only shows transitions between a pure radiation like model and the theory without galileons, this theory can not display selfacceleration. Hence for it to give accelerated expansion a cosmological constant or some other form of dark energy would need to be added.
That the impact of the galileon decreases for increasing values of c 2 can also be seen by comparing plots of the fractional densities and densities for two different values of c 2 . This is shown for c 2 = 10M 2 PL and 10 8 M 2 PL in figures 10, 11, 12 and 13 respectively.

Non-zero value of c 3 and added cosmological constant
The theory considered in [22] corresponds to a model of non-zero c 3 , c 1 = c 2 = 0 with the addition of a cosmological constant. That is, a theory of a third order only galileon in a Since the galileon density starts out negative its absolute value |ρ π | is shown in the plot. and δ given in equation (4.7) become: and the equation for H ′ becomes: We can obtain a model equal to the one considered in [22] by considering this regime with c 1 = c 2 = 0, c 3 = 0 in a background with a cosmological constant. However, since we had a slight discrepancy in the equations as compared to their results, there might be some differences. We must also pay attention to the units, as we are working with dimensionless π fields, where as they have dimensionful π fields.
To review the results and compare to those of [22], we have plotted the evolution of the Hubble parameter H as compared to that of ΛCDM, the evolution of the energy densities ρ i and the evolution of the fractional densities Ω i in figures 14, 15 and in figure 16 for a particular value of c 3 , c 3 = 100M 2 PL /H 2 0 . As in [22] we have defined (5.14) To make the results even more comparable to those in [22] we make the Ω Λ value such that Ω m today is close to the value measured today. We also extend the simulation to z = −0.8.
There is also a cosmological constant present in this model set so as to ensure Ω m today of about the same size as measured. This corresponds quite closely to the case considered in [22]. The plot also shows the evolution of the Hubble parameter in the ΛCDM case. The results resemble closely the ones found in [22]. Similarly to [22] we also found that a correct evolution history was first obtained for relatively large values of c 3 , the models resembling ΛCDM more the larger the value of  Here the added cosmological constant has been set so as to ensure Ω m today close to its measured value, as in [22]. of c 3 . In these plots the cosmological parameter has, however, only been set to its ΛCDM value, as normalising it in the very small c 3 cases is quite hard due to the large impact of the galileon field.
The tendency is that for small enough values of c 3 the galileon field dominates the evolution too early for the Universe ever to experience a period of accelerated expansion. This should not come as a surprise after considering the case of all the coefficients set to zero in the previous subsection. In this case we have seen that the Universe will behave as if radiation dominated at all times no matter what the matter content of the Universe may be. The fact that the deviation from ΛCDM is larger with smaller c 3 is hence quite understandable, since c 3 is a demarcation scale, under which the physics still behaves as in a non-galileon Universe. The role of c 3 is hence quite similar to that of c 2 as described in section 5.3 and we would expect the model with c 2 only with a cosmological constant added to have results somewhat comparable to the results found here.
To check whether the behaviour is altered due to our slightly different results for the Raychaudhuri equation (3.9), we have performed simulations using their equation set as well as our own. They yield only slightly different results, i.e. the general behaviour is not changed and only very subtle changes in the evolution can be observed by reviewing the data very carefully.
We have also considered the effect of altering the initial values of the galileon field. Actually, as long as the initial galileon field is small and positive π start < 0.1, and the initial galileon field prime derivative is quite small and negative 0 > π ′ start > −1e − 8, the evolution histories found for each value of c 3 is quite stable, that is with the initial values within these bounds, their exact values have little impact on the ensuing evolution history of the Universe. Relatively small initial values of the galileon field and its derivative seem  like a reasonable starting point, so assuming this model to have the described behaviour seems fairly acceptable.

Numerical results for the full model background equations
Moving on to the full model, with c 1 , c 2 and c 3 obeying the constraints of equation ( which is supposed to show selfacceleration, we obtain an evolution for the Hubble parameter H as shown in figure 21. We see that the evolution of the Hubble parameter in this case is quite comparable to that of the ΛCDM model.
However, the evolution in this theory is not exactly as that of a ΛCDM model. This can be seen for instance from the evolution of the densities and fractional densities of the different components as shown in figures 22 and 23, or in the evolution of the equation of state parameter shown in figure 24. 4 This interesting behaviour may lead to possible ways of discerning this model from the ΛCDM model or at least to possibilities for constraining the values of c 2 and c 3 from distance measurements, such as supernova data or large scale structure surveys showing the evolution of the equation of state parameters over time.

Robustness of results
Not surprisingly, the model does not stably yield near ΛCDM results under all initial values of the galileon field π and its derivative. If these have values that are too large in the beginning, the early Universe evolution will be dominated by the galileon and hence not give results compatible with current knowledge from baryon nucleosynthesis etc.
However, the evolution of the background equations also exhibit a dependence on the exact values of the coefficients c 2 and c 3 . If c 2 is too small, and this seems to not only depend on the relative value with respect to c 3 , the selfacceleration starts too late in the evolution of the Universe. If c 2 is too large the selfacceleration starts too early. The fine  tuning needed for c 2 increases as c 3 gets much smaller or larger than a few M 2 PL /H 2 0 . We can summarise by saying that for the cosmological evolution to be close to our own, the values of the parameters need to be given by order a few times appropriate factors of H 0 and M PL . This makes it tempting to say that although the parameter space is constrained, the allowed values constitute a natural parameter regime.
After considering the galileon model with each term separately in sections 5.2 to 5.4 this is quite understandable. Since c 2 and c 3 are demarcation scales, keeping the theory away from the pure radiation dominated like Universe, at least one of these must be not too small to yield sensible results. Since our stability conditions dictate that c 2 > 3H 2 0 c 3 , we get that c 2 must be large enough. Since c 1 is the only term that can in fact yield selfacceleration, this term must be present, and to get the right value of the Hubble parameter today, and the right evolution, c 1 must be set by the demarcation scales. Having the demarcation scales both too large also causes trouble, since c 1 then gets very large and starts the accelerated expansion too early. . We see that after the radiation and matter dominated epochs, we get a slight increase in the effective equation of state parameter before it dips down to negative values. After this it starts increasing again. We can interpret the plot as follows. At early times the evolution follows that of the standard ΛCDM model, starting with radiation domination and then going into matter domination. As the galileon starts to dominate the evolution it has a negative density and is dominated by the coupling to gravity hence driving the theory towards a radiation dominated like epoch like in the all coefficients zero case. After only a short while of this however, the c 1 term takes over leading to a period of accelerated expansion. As time continues to go by we cross the demarcation scale given by the c 2 and c 3 entering an epoch where selfacceleration will stop and eventually turn to deceleration again.
For the initial values of the galileon and its derivative we find that they must be tuned to be quite close to zero. In the case of c 3 = 3M 2 PL /H 2 0 , c 2 = 18M 2 PL , the evolution remains stable as long as the initial value for π ′ < 1e − 20 with the initial value of π < 1e − 9. As we change only the initial value of π to around 1e − 8, the evolution shortens its matter dominated epoch and gets a more abrupt change to the selfaccelerated phase. Outside this regime the evolution becomes dominated by the galileon from a very early stage, and it becomes totally unstable. This is similar to the results found for perturbations in [50], where the initial values for the galileon had to be tuned to zero in the beginning in order for the galileon fluctuations not to take over and spoil the homogeneity and isotropy of the background.

Discussion
In this paper we have explored a certain version of the galileon model formulated in [21]. Using their derivation of stability and selfaccelerated expansion, we have found a model that both admits selfacceleration and gives stable spherically symmetric solutions, that is a third order galileon with certain constraints on the three parameters. This model has not been thoroughly studied in the literature, and we feel that it is important to provide such a study of one of the simplest proposals for a selfaccelerating galileon.
We have explored the background cosmology of this model, using numerical simulations of each term separately to build intuition. In this way we found that it is the tadpole term that leads to the selfacceleration of the solution. The second and third order terms defines a sort of demarcation scale above which the selfacceleration slows down, creating an evolution history not totally equal to that of ΛCDM. We also found that the galileon density becomes negative before it starts dominating. As it passes through zero and on to positive values the evolution of the equation of state goes through a phase where it grows slightly. This is a feature of the chosen coupling to gravity, which we have seen yields radiation dominated like solutions if present with no other galileon dynamics.
We also found that the parameter space for a solution comparable to ΛCDM is limited, but not totally constrained. The bounds on the parameters to give a ΛCDM like evolution history are, however, consistent with solar system and galaxy cluster gravitational bounds found in [21] and they lie in a somewhat natural range. We presume that further studies on the perturbations of this model will set even more stringent bounds on these parameters.
We also found that the initial values of the galileon field and its derivative must be tuned to values quite close to zero. With larger values of these initial field, the model becomes highly unstable as the galileon begins to dominate quite early, leading to cosmological histories very different from the one favoured by present day observations. This is similar to results found in [50] where the galileon perturbations had to be tuned to exactly zero in the beginning to avoid a total destruction of the homogeneity and isotropy of the background space. Further studies on the perturbations of this model will be interesting in finding out whether the same kind of anisotropic stress will be found here.
The growth of structure in the galileon Universe will also be of interest as it might provide nice signatures in both baryonic acoustic oscillations and CMB, and these can only be found through studies of the perturbation equations of this model, which we postpone for a future publication.