Bi-galileon theory I: motivation and formulation

We introduce bi-galileon theory, the generalisation of the single galileon model introduced by Nicolis et al. The theory contains two coupled scalar fields and is described by a Lagrangian that is invariant under Galilean shifts in those fields. This paper is the first of two, and focuses on the motivation and formulation of the theory. We show that the boundary effective theory of the cascading cosmology model corresponds to a bi-galileon theory in the decoupling limit, and argue that this is to be expected for co-dimension 2 braneworld models exhibiting infra-red modification of gravity. We then generalise this, by constructing the most general bi-galileon Lagrangian. By coupling one of the galileons to the energy-momentum tensor, we pitch this as a modified gravity theory in which the modifications to General Relativity are encoded in the dynamics of the two galileons. We initiate a study of phenomenology by looking at maximally symmetric vacua and their stability, developing elegant geometric techniques that trivially explain why some of the vacua have to be unstable in certain cases (eg DGP). A detailed study of phenomenology appears in our companion paper.


I. INTRODUCTION
When Urbain Le Verrier first noticed that the slow precession of Mercury's orbit could not be explained using perturbative Newtonian gravity and the known planets, he postulated the existence of a hitherto unseen dark planet, which he called Vulcan [1]. Of course, thanks to Einstein, we now know that what was really required was a new theory of gravity, specifically, the General Theory of Relativity [2]. 150 years later, flat rotations curves of galaxies [3] and accelerated cosmic expansion [4] have led most of us to follow Le Verrier's approach and to postulate the existence of dark matter [5] and dark energy [6] respectively. For dark matter, bullet cluster observations suggest that Le Verrier's approach is probably the right one [7]. The Le Verrier approach to dark energy is less compelling. Whereas extensions of the Standard Model offer some very good dark matter candidates, they have been nothing short of hopeless in offering a satisfactory candidate for dark energy. It is certainly worth taking a lesson from history and exploring the possibility that a new theory of gravity may be required to explain the observed cosmic acceleration.
In order to play a role in the accelerated expansion, we need gravity to deviate from General Relativity (GR) at large scales, and in particular at the scale of the cosmic horizon H −1 0 ∼ 10 26 m. Such deviations should not be felt all the way down to arbitrarily short distances because General Relativity is very successful in describing classical gravity up to the scale of the Solar System, at about 10 15 m. Modifications of General Relativity are most easily understood at the level of linearised theory. In GR, the gravitational force is mediated by a massless spin 2 particle carrying two polarisation degrees of freedom. In modified gravity we expect additional degrees of freedom to be present. For example, in Pauli-Fierz massive gravity the graviton carries 5, rather than 2, polarisation degrees of freedom [8]. Clearly if they are to play a role in accelerated expansion without spoiling the phenomenological successes of GR at shorter distances, the additional degrees of freedom should be suppressed up to Solar System scales, but not on cosmological scales. This should occur without introducing any new pathological modes. In massive gravity, the longitudinal graviton mode is problematic since it gives rise to the so called vDVZ discontinuity in the graviton propagator, even at scales where the graviton mass is negligible [9]. However, this mode can be "screened" on Solar System scales thanks to non-linear interactions and the Vainshtein mechanism [10], but only at the expense of introducing a new pathology in the form of the Boulware-Deser ghost [11,12].
Extra dimensions and the braneworld paradigm [13] have opened up some intriguing possibilities for modifying gravity at large distances (see, for example, [14][15][16][17]). The DGP model [14] is probably the most celebrated example of this. The model admits two distinct sectors: the normal branch and the selfaccelerating branch. The latter generated plenty of interest since it gave rise to cosmic acceleration without the need for dark energy [18]. However, fluctuations about the self-accelerating vacuum were found to suffer from ghost-like instabilities, ruling out this sector of the theory [19,20]. Although the normal branch is less interesting phenomenologically, it is fundamentally more healthy and is the closest thing we have to a consistent non-linear completion of massive gravity. Here the graviton is a resonance of finite width H 0 , as opposed to a massive field. At short distances, r ≪ H −1 0 the brane dynamics does not feel the width of the resonance and the theory resembles 4D GR. At large distances, r > ∼ H −1 0 , however, the theory becomes five dimensional as the resonance effectively decays into continuum Kaluza-Klein modes. The Vainshtein mechanism works well on the normal branch of DGP, screening the longitudinal graviton without introducing any new pathological modes, in contrast to massive gravity [12,21].
Much of the interesting dynamics is contained in the properties of the longitudinal scalar, π. One can identify this mode with the brane bending mode. As this mode only excites the extrinsic curvature of the brane, and not the brane or bulk geometry, it becomes strongly coupled well below the 4D Planck scale [20]. It is precisely this strong coupling that allows the Vainshtein mechanism to take effect as non-linearities become important at larger than expected distances, helping to shut down the scalar on Solar System scales [12,23]. We can isolate the dynamics of the scalar by taking the so-called decoupling limit around flat space, in which the 4D and 5D Planck scales go to infinity, but the strong coupling scale is held fixed. In this limit the tensor mode, h µν , decouples from the scalar, π, to all orders [20,22]. The π-Lagrangian now contains a plethora of information about both branches of solution, which can be studied much more easily than in the full theory.
Now an important feature of the π-Lagrangian is that it exhibits galileon invariance [24]. That is, the action is invariant under π → π + b µ x µ + c, where b µ and c are constant. This essentially follows from the fact that the brane bending contribution to the extrinsic curvature is given by K µν ≃ ∂ µ ∂ ν π. Equivalently, Poincare transformations in the bulk should not affect the physics of the brane bending mode. If the brane position is given by y = π(x), then after an infinitesimal Poincare transformation, x a → x a + ǫ a b x b + υ a the brane position is given by y ≃ π(x) + ǫ y µ x µ + υ y . Since the physics of the brane bending mode should remain unaffected it is clear that it should exhibit galileon invariance.
Motivated by the DGP model, Nicolis et al [24] developed the idea of galileon invariance and wrote down the most general Lagrangian for a single galileon degree of freedom (see [25] for some early work along similar lines). Remarkably the theory contains only five free parameters, enabling the authors, and others [26], to complete a thorough analysis. Although the Galilean symmetry is broken when we include coupling to gravity [27], covariant completions of the general galileon action can be identified with the position modulus of a probe DBI brane [28]. In any event, it seems clear that any co-dimension one brane model with large distance deviations from GR, will indeed admit a galileon description around flat space, in some appropriate limit. This follows from the discussion of the brane bending mode in the previous paragraph, and from [23], where it is argued that any viable infra-red modification of gravity leads to strong coupling, and is therefore expected to admit a non-trivial decoupling limit.
Recently there has been plenty of interest in branes of higher co-dimension and the role they might play in modifying gravity at cosmological scales [29][30][31][32][33][34][35][36][37][38]. An example of this is the cascading DGP model [29], in which one has DGP branes within DGP branes [29][30][31]. In the 6D version one has a co-dimension 2 DGP 3-brane intersecting a co-dimension one DGP 4-brane [29]. Higher co-dimension brane models seem to have extra features, such as degravitation [30,39], or perhaps even self-tuning [29,[32][33][34], in which the geometry of the 3-brane is insensitive to a large 3-brane vacuum energy. Now to specify the position of a co-dimension-N brane we need to specify the values of N "bulk" components: y 1 = π 1 (x), y 2 = π 2 (x), . . . , y N = π N (x). For a viable model giving rise to large distance deviations from GR, there should be strong coupling and therefore an appropriate limit in which the brane bending modes decouple from the tensor mode around flat space. Then the physics of each brane bending mode should be unaffected by Poincare transformations in the bulk, and we conclude that we are left with an N -galileon theory.
In this paper we will focus on the case of two coupled galileons, which we will call bi-galileon theory. Of course, this has particular relevance to co-dimension 2 brane models. Indeed, in section II, we will show explicitly that the 5D cascading cosmology model is one example of a bi-galileon theory. This model is derived from the original 6D cascading gravity model [29,30] by computing the boundary effective field theory in 5 dimensions, and taking a suitable decoupling limit. Although a detailed study of phenomenology will appear elsewhere [41], this paper is devoted to the motivation and formulation of the most general bi-galileon theory. In section III we will pitch our theory as a modification of gravity by coupling one of the galileon fields to the energy-momentum tensor. We will not couple the graviton directly to the scalars, which amounts to ignoring their backreaction on to the geometry, as was done for the case of the single galileon [24]. In section IV we will construct the general bi-galileon theory. We will also state the extension to any number of coupled galileon fields. We initiate a study of phenomenology in section V by looking at maximally symmetric vacua and their stability. Section VI contains a half-time analysis.

II. EXAMPLE: CASCADING COSMOLOGY
There have been a number of attempts to extend the DGP model to higher dimensions. Although the earliest of these ran into problems with ghosts and instabilities [40], the recently developed cascading models [29][30][31] have been more successful. These models exhibit degravitation and/or self-tuning, and seem to be free from ghosts, at least if the 3-brane tension is large enough. In this section we will compute the 4D boundary effective theory for the cascading cosmology model [31]. By taking the appropriate limit we will decouple the scalar sector and show that it corresponds to a bi-galileon theory.
The cascading cosmology model [31] is closely related to the 6D cascading DGP model [29]. The latter contains a DGP 3-brane sitting on a DGP 4-brane in 6 dimensions. The action is given by where the full 5D spacetime is made up of two copies of M, identified across a common boundary, ∂M, given by the 3-brane. Here g ab dx a dx b and q µν dx µ dx ν are the bulk and brane metric respectively, with corresponding Ricci scalars R 5 and R 4 . The extrinsic curvature of the brane is given by K µν = 1 2 L n q µν where L n is the Lie derivative with respect to the inward pointing normal, n a , to ∂M. The bulk scalar π(x) is the 5D longitudinal mode and the parameter α = 27 4m 2 6 . If we choose coordinates x a = (y, x µ ), so that the brane is at y = 0 and the bulk extends into y > 0, then it is convenient to write the bulk metric in ADM coordinates where N is the lapse function and N µ is the shift vector. Then where D µ is the covariant derivative with respect to the brane metric q µν .
To calculate the boundary effective action, we essentially follow steps completely analogous to the ones taken in [20]. In order to integrate out the bulk fields we must choose a gauge. We therefore set and add a bulk gauge fixing term where Classically this imposes the gauge F a = − 3 2 ∂ a π. To quadratic order we have where and we have re-defined N µ = h µy . The bulk equations of motion have solutions where ∆ = √ − 4 . Now we integrate out the bulk and fix the residual gauge freedom left on the brane by adding a brane gauge fixing term, This gauge fixing eliminates the mixing between N µ and f µ and we are left with We can diagonalise this action by making the following field re-definitions This gives Note that the π field is a ghost. This can be avoided by increasing the vacuum energy on the brane [29]. In addition, to quadratic order, the matter coupling is given by We now consider higher order terms. Schematically it is clear from Eq. 5 that these will take the form where a i , . . . , d i are non-negative integers satisfying a i + b i + c i + d i ≥ 3. Note that we must have a 2 ≥ 1 since terms of the form Eq. 21 stem from the boundary π term given in Eq. 5. After canonically normalising the fields, it is easy enough to check that the largest contribution comes from a term like Eq. 20, with This term is exactly like the cubic term computed in [20], and comes from the following term in the action (5) Since we want a 1 = 0, we can switch off e −3π/2 , and write this term in ADM form, Of course, we also want b 1 = 0, so we need to switch off h µν . This means q µν = η µν , and so R 4 = 0 and . Given that N ≈ 1 + 1 2 h yy , the cubic order interaction term is given by In the limit where m 5 → 0, the strongly interacting mode can be parametrized in the bulk as where 5 Ψ = 0. This satisfies the bulk equations of motion, as well as the gauge choice F a = − 3 2 ∂ a π. Indeed, the strongly interacting mode corresponds to a trivial bulk geometry, as expected [20]. Furthermore, we see from Eq. 17 that the boundary values for the fields go like N µ ∼ 1 m5 ∂ µ (π + φ), and h yy ∼ − 2 m5 ∆(π + φ) so we take Ψ = e −∆y π+φ m5 . Eq. 25 can now be written as Now if we canonically normalise, so that π = 2 3π M4 and φ = 2 3φ M4 , then we see that the cubic interaction term goes like where This represents the scale at which the 4D scalars becomes strongly coupled, and corresponds to a distance Λ −1 5 ∼ 1000km. We can now take the formal limit 2 Note that this requires that m 5 , m 6 → 0. It is easy to see that the leading order matter coupling remains. In contrast, the large interaction term Eq. 28 is the only interaction term to survive this limit. As a result, the tensorh µν and vectorÑ µ decouple from one another and from the scalars π and φ to all orders. In terms of the canonically normalised fields, the action is given by µ . Focussing on the scalar sector, the relevant piece of the action is given by The resulting equations of motion are explicitly Galilean invariant, That is, they are invariant underπ for constant b µ ,b µ , c andc [24]. Thus we confirm what we expected. Cascading cosmology admits a bigalileon description in the full decoupling limit. The limit corresponds to taking the Planck masses to infinity but holding the 4D and 5D strong coupling scales fixed. The latter correspond to 1000 km and 1000000 km respectively, whilst the Planck lengths in the original model are all submillimeter. In taking the decoupling limit we have sent all the next largest interactions to zero, so our description is only valid at scales where these additional interactions may be ignored, up to the horizon scale. These interactions may have an important role to play in understanding various aspect of the full theory (eg. Vainshtein effects), but since our goal was merely to demonstrate that a bi-galileon description exists in principle at certain scales, we do not concern ourselves with that here. The π, φ-Lagrangian associated with Eq. 32 is the analogue of the π-Lagrangian in DGP gravity [20,22]. It should capture local physics on the 3-brane, corresponding to normalisable fluctuations in the 5D bulk. Note that this means it cannot capture the self tuning behaviour of the full 6D theory. This is because the self tuning solutions correspond to non-local, non-normalisable fluctuations about the vacuum of the effective 5D theory described by the cascading cosmology model (5). We could, of course, consider normalisable fluctuations about a non-trivial self tuning solution. We expect that this would also give a bi-galileon theory in the appropriate decoupling limit. For sufficiently large 3-brane tension it would actually correspond to a ghost-free bi-galileon theory [29], unlike the example given here. However, a detailed analysis is outside of the scope of this paper.

III. BI-GALILEON MODIFICATION OF GRAVITY
The cascading cosmology model described in the previous section is an example of a modification of General Relativity by two additional galileon fields, π and φ. In a generic theory of modified gravity, the amplitude for the exchange of one graviton between two conserved sources, T µν and T ′ µν , is given by where D µναβ is the graviton propagator. Now, for General Relativity, the amplitude is given by where T = T µ µ . We will be interested in the case where gravity gets modified by two additional scalars so that locally we have Such a theory can be described by the following action whereg µν = η µν +h µν . The fluctuationh µν is identified with the graviton in General Relativity since it satisfies the same linearised field equations. Therefore, for a given source and boundary conditions, the solutionh µν coincides with the linearised GR solution. This statement is true to all orders in the limit, where the scalars decouple fromh µν . The physical graviton differs from its GR counterpart by a scalar fluctuation h µν =h µν + 2(π 1 + π 2 )η µν , since the physical metric seen by matter is given by g µν = η µν + h µν . Let us consider the decoupling limit (40), with the additional assumption that the strength of all scalar interactions are held fixed. This basically means we are neglecting the backreaction of the scalars on the geometry, so we can consider them as fields on Minkowski space. Note that we retain scalar self-interactions because they are essential for bi-galileon corrections to GR to be extrapolated from O(10 −5 ) in the solar system to O(1) at cosmological scales. Significant modifications to gravity at Hubble distances would inevitably have a role to play in the dark energy problem. The action is given by where Eh µν = − 1 2 4 h µν − 1 2h η µν + . . . is the linearised Einstein tensor, and the Lagrangian L π1,π2 only depends on the scalars π 1 and π 2 . As we can see from Eq. 31, the action for the cascading cosmology certainly takes this form in the decoupling limit. In Eq. 31 the scalar action is Galilean invariant. To generalise this, we are interested in bi-galileon modifications of gravity for which L π1,π2 is the Lagrangian for a bi-galileon theory. This means that if π i → π i + (b i ) µ x µ + c i then L π1,π2 → L π1,π2 + total derivative.
Given the form of the physical graviton (41), it is convenient to define new galileon fields π = π 1 + π 2 and ξ = π 1 − π 2 , so that only one of the scalars couples to matter. The action now takes the form where L π,ξ corresponds to a bi-galileon theory. The physical graviton is now Although only π couples to matter directly, ξ couples to matter indirectly, through its mixing with π. It is possible for the two scalars to be mixed at quadratic order, as well as at higher order. Given that we have taken all our fields to propagate on Minkowski space, one might be forgiven for thinking that we cannot say anything about cosmology, or perturbations about cosmological solutions, including de Sitter vacua. Fortunately, this is not the case. We can think of any spatial flat cosmological spacetime as a local perturbation around Minkowski space [24]. In this context, local means local in both space and time, at sub-horizon distances and over a sub-Hubble time. If we take here to be x = 0 and now to be t = 0, then for | x| ≪ H −1 and |t| ≪ H −1 , we have [24] where the Hubble scale H and its time derivativeḢ may be evaluated at t = 0. This is a perturbation about Minkowski space written in Newtonian gauge where the Newtonian potentials are For a given cosmological fluid, the corresponding GR solutions have Hubble parameter H GR . Sinceh µν agrees with the linearised GR solution, we haveh tt = −2Φ GR ,h ij = 2Φ GR δ ij , where Now in our modified theory the physical Hubble parameter (associated with h µν ) differs from the corresponding GR value, H = H GR . Due to Eq. 44, we have a non-trivial scalar The other scalar, ξ, follows from the equations of motion. Note that a Galilean transformation π → π + b µ x µ + c merely corresponds to a coordinate transformation Of particular interest are maximally symmetric vacua, for which H and H GR are constant. Indeed, if the vacuum energy is σ, we have H 2 GR = σ/3M 2 pl . It is easy to check that one scalar now takes the particularly simple formπ = − 1 4 k π x µ x µ where k π = H 2 − H 2 GR . In principle the other scalar can be anything, but we will typically assume that it takes the same form on maximally symmetric solutions, so thatξ = − 1 4 k ξ x µ x µ . Ultimately, we would like to promote our galileon theory to a fully covariant theory of gravity, coupled to a pair of scalar fields. Whilst this will inevitably break the Galilean invariance [27], we should expect to recover the galileon description in the decoupling limit M pl → ∞. One does have to be slightly careful not to introduce any higher derivatives by accident when covariantising the theory, but this can be done using additional non-minimal couplings [27]. As a step towards this, let us first transform the action (43) to "Jordan" frame by rewriting everything in terms of the physical metric (44). This gives whereL π,ξ = L π,ξ − 3M 2 pl π∂ 2 π. The covariant completion of this theory will take the form where L matter [g; Ψ n ] is the matter Lagrangian, describing matter fields, Ψ n , minimally coupled to gravity. The scalar contribution, L scalar [g; π, ξ], is constructed as follows: takeL π,ξ and let η µν → g µν , ∂ µ → ∇ µ , then add non-minimal couplings to curvature to compensate for the higher derivative terms in the equations of motion. Actually, this completion is useful in that it enables us to establish the validity of our galileon description. Recall that this neglects the backreaction of the scalars onto the geometry. From Eq. 51, we see that the scalars contribute an energy-momentum tensor The gravitational field equations arising from the covariant action (51) are given by Although the galileon description includes non-linear effects at the level of the scalar field equations, it neglects them at the level of the gravity equation. Expanding Eq. 53 about Minkowski space, g µν = η µν +h µν , we find that where ". . ." denotes terms that are O(h) suppressed relative to the terms on the LHS of Eq. 54. The gravity equation in the galileon description corresponds to neglecting everything on the RHS of Eq. 54. When is this justified? For a small graviton fluctuation, we can consistently neglect those the terms that are O(h) suppressed, denoted by ". . .". However, even then, the backreaction of the scalars on to the geometry can be significant. In other words, it is not clear that the energy momentum of the scalars, evaluated on Minkowski, can be consistently neglected. To justify neglecting this backreaction we conservatively require that T µν scalar [η, ; π, ξ] ≪ M 2 pl Eh µν . Provided this condition holds, we will see that our bi-galileon theory allows for a rich and interesting phenomenology, to be studied in detail in our companion paper [41]. As we have just discussed, we can describe cosmological solutions, such as de Sitter space, as local perturbations about Minkowski. We will also be interested in fluctuations about these solutions, and in particular about non-trivial maximally symmetric vacua. These just correspond to fluctuations in the galileon fields about their background values. As we mentioned earlier, the non-linearity of the galileon action will be essential to the success of escaping the local gravity tests and producing interesting cosmological effects simultaneously.

IV. CONSTRUCTING THE BI-GALILEON THEORY
Let us now construct the most general bi-galileon theory. Our task is to find the most general Lagrangian describing two scalars, L π,ξ , that is invariant under the Galilean transformation This amounts to requiring that the corresponding equations of motion are built exclusively out of second derivative terms ∂ µ ∂ ν π, ∂ α ∂ β ξ contracted with the background Minkowski metric . Single, or zero derivative contributions to the equations of motion break the Galilean invariance, whereas higher derivatives introduce extra degrees of freedom and ghosts. The Galilean symmetry greatly constrains the structure of the Lagrangian. For example, under an infinitesimal constant shift π → π + c, the Lagrangian L π,ξ varies as However, as this shift happens to be a Galilean transformation, we know that any change in the Lagrangian must be a total derivative. It follows that ∂L π,ξ /∂π = ∂ µ J µ , for some J µ . Now the π equation of motion goes like which is clearly a total derivative ∂ µ [J µ − ∂L π,ξ /∂(∂ µ π) + . . .]. Indeed, the π equation of motion is a total derivative at each order in π and ξ. We can apply a similar argument to the ξ equation of motion.
We can now formulate the problem as follows: at (m, n)-th order (m-th order in π and n-th order in ξ), what is the most general total derivative built out of second derivative terms? One obvious total derivative is given by It is easy to see all such terms with m + n > 4 vanish in four dimensions. We will now write out all the non-vanishing E m,n explicitly in 4D. For notational convenience, we suppress the Lorentz indices and separate Lorentz scalars by '·' when it is confusing. For example, we have ∂∂π∂∂ξ · ∂∂π∂∂ξ ≡ For m ≥ n, the non-vanishing terms are For m < n, we take E m,n = E n,m | π↔ξ . Note that in each case E m,n = ( π) m ( ξ) n + . . ..
We will now prove that Eq. 58 is the only total derivative of order (m, n) that is second order in spacetime derivatives. The proof works along the same lines as the corresponding single galileon one given in [24]. First, assume there is another suitable total derivative E ′ m,n . Now consider E ′ m,n as a Lagrangian density. As E ′ m,n is a total derivative, the resulting π and ξ equations of motion must be trivial Without loss of generality, a general term in E ′ m,n goes like 3 where i 0 + . . . + i k = m, j 0 + . . . + j k = n. Let us now consider its contribution to the resulting equations of motion (60). A term in the π equations of motion can be obtained by shifting ∂∂ (along with the Lorentz indices) from one ∂∂π to another ∂∂π or another ∂∂ξ and then dropping the π without derivatives. Therefore Eq. 61 must contribute a term proportional to in the equation of motion. As the overall equation of motion is trivial, this term must be cancelled by the variation of another term in E ′ m,n . Aside from Eq. 61, there is only one other possible term that can produce a term like Eq. 62 in the π equation of motion. It is given by It follows that this must also be in E ′ m,n . One can easily imagine a similar argument applies to ξ. Now, if we apply these arguments iteratively, we would conclude that a term proportional to ( π) m ( ξ) n must be in E ′ m,n . Then we immediately run into a contradiction, because linearly combining E m,n and E ′ m,n should produce another total derivative which is free of ( π) m ( ξ) n . This proves that the only suitable total derivative at order (m, n) is E m,n .
Returning to the full theory, the equations of motion for π and ξ in general correspond to linear combinations of the E m,n : Since these equations are derived from the same Lagrangian, L π,ξ , we anticipate an integrability condition relating the coefficients a m,n and b m,n . To construct the Lagrangian, we note that equations of motion with the correct form can be derived from the following where L m,n = (α m,n π + β m,n ξ)E m,n .
As expected there exists an integrability condition relating the coefficients in the equation of motion. One can easily check that The origin of this condition lies in the fact that πE m−1,n − ξE m,n−1 is a total derivative for n, m ≥ 1. This means we have included a number of total derivatives in the Lagrangian that do not contribute to the dynamics. Indeed, at order m + n + 1 ≥ 2 the total number of free parameters is not 2(m + n + 1) as one might naively expect, but m + n + 2 owing to the constraints given by Eq. 71. We now have the most general bi-galileon theory, given by Eqs. 65, 66 and 58. Referring back to the discussion in the previous section, we note that only π couples directly to matter, through an interaction d 4 x πT . It follows that the equations of motion are given by where the a m,n and the b m,n are related according to Eq. 71. We conclude this section by stating, without proof, the generalisation to N galileons. The analogue of E m,n is Therefore, the general N -galileon Lagrangian is given by where L m1,...,mN = (α 1 m1,...,mN π 1 + . . . + α N m1,...,mN π N )E m1,...,mN . (76)

V. MAXIMALLY SYMMETRIC VACUA AND THEIR STABILITY
Although a detailed study of bi-galileon phenomenology will be left to our companion paper [41], we will now initiate a study of maximally symmetric vacua. We will be particularly interested in establishing stability conditions associated with ghost-like excitations of the galileon fields about these vacua.
We will allow for a non-zero vacuum energy so that the background stress energy tensor is given bȳ T µ ν = −σδ µ ν . As discussed in section III, a maximally symmetric solution has a background π field, where k π = H 2 − H 2 GR is the difference between the Hubble parameters calculated in the modified gravity theory and in GR. For simplicity we assume that the background ξ field takes the same form for some constant k ξ . If our theory happens to be derived from a brane model, we expect this choice to correspond to choosing maximal symmetry in the bulk and on the branes. We now plug this into our field equations Eqs. 72 and 73. This yields where we recall that a m,n and b m,n are given by Eqs. 69 and 70. Generically these represent two independent equations in two unknowns. Being quartic polynomials, their solutions will typically contain sixteen distinct points in the real (k π , k ξ ) plane. Of course, there may be degenerate configurations with a continuum of solutions, or other configurations with no real solutions whatsoever. We will study these aspects in more detail in [41]. For the moment, let us assume that a real solution exists. Can we really trust our solution? We can think of two things in particular that could spoil the background: gravitational back reaction and/or dynamical instabilities. The former refers to backreaction of the scalars onto the geometry whereas the latter refers to ghost-like excitations. To succinctly express suitable conditions on the background, it is convenient to consider what happens when we plug our background fields directly into the action. It turns out that where we define the action polynomial, Now, it is easy to check, using Eqs. 69 and 70, that Of course this comes as no surprise since the equations of motion are designed to extremise the action. For the maximal symmetry case in question this simply corresponds to extremising L. Now let us return to the issue of backreaction of the scalars on to the geometry. This was discussed at the end of section III. Schematically, we expect Recall that to justify neglecting the backreaction, we require T µν scalar ≪ M 2 pl Eh µν ∼ M 2 pl H 2 . Assuming that we can indeed neglect the backreaction as described, we then need to ask if our background solution is stable against excitation of ghosts? To study the stability of a solution we need to examine its fluctuations, so we shift the fields so that they represent fluctuations on the appropriate background π →π + π , ξ →ξ + ξ .
We claim that the background solution has a ghost-like instability if the Hessian of L(k π , k ξ ) has a negative eigenvalue on the solution. Let us now prove this claim.
The field re-definition (85) takes us from the trivial to a non-trivial background, and gives rise to the following shift in the Lagrangian It is clear that the perturbed Lagrangian L ′ π,ξ has the same symmetries as L π,ξ , so we may once again write where L ′ m,n = (α ′ m,n π + β ′ m,n ξ)E m,n .
Note that, due to the background equations of motion, linear terms in L ′ π,ξ vanish; that is there is no contribution from m = n = 0. The equations of motion are given by with the relations a ′ m,n = (m + 1)(α ′ m,n + β ′ m+1,n−1 ) and There is a linear map between the parameters in our original theory a m,n , b m,n and the perturbed theory a ′ m,n , b ′ m,n . By direct substitution we find that for 1 m + n 4, we have where Here the combination For N < 0, we extend the definition of "factorial" using the Gamma function, N ! = Γ(N + 1), and recall that 1/Γ(N + 1) = 0 when N is a negative integer. It immediately follows that we only get contributions from i + j 4.
To study the stability of the background, we consider the leading order piece of the perturbed Lagrangian The background has a ghost-like excitations if has a negative eigenvalue. However, using Eqs. 92, 93 and the linear map Eq. 95, we see that Differentiating Eqs. 79 and 80 then using Eq. 83, one is able to show by direct comparison that Recall that the background solution contains a ghost if K has a negative eigenvalue, which is equivalent to Hess(L) also having a negative eigenvalue. Thus we have proven our claim 4 . We can understand this result geometrically. The function L(k π , k ξ ) is a smooth polynomial function of two variables, whose extrema correspond to maximally symmetric vacua in our theory. Unstable vacua happen to be those for which the Hessian of L has a negative eigenvalue on the solution. This just means that the corresponding extrema are saddles and maxima. Or put another way, stable vacua are simply minima of the action polynomial, L. Note that a vanishing eigenvalue will generically lead to strong coupling at all scales as the corresponding scalar degrees of freedom will have a vanishing kinetic term.
This geometrical insight enables us to make some generic statements about the spectrum of vacua in each case, using Morse theory. For example, if our galileon theory contains more than one vacuum solution, they cannot all be stable. This is because Morse theory tells us that you cannot have all extrema being minima if there is more than one extremum. This logic explains why the self-accelerating branch of the DGP model must have a ghost-like instability [19] if the normal branch is stable.
The action polynomial is an extremely useful construct. One can use it to build interesting theories from the bottom up, ensuring stability conditions from the outset. This will be the approach we take in our detailed study of phenomenology [41]. Given a suitable action polynomial, we can reconstruct the corresponding action for the relevant theory by computing the following coefficients: where m, n = 0, 1, 2, . . ., and we recall that we define α −1,n = β m,−1 = 0. Note that since πE m−1,n − ξE m,n−1 is a total derivative for n, m ≥ 1, we are free to set, say, β 1,n = β 2,n = β 3,n = . . . = 0, without loss of generality. However, it is clear from Eq. 101 that this choice subsequently fixes α m,n uniquely. The action polynomial can also be used to quickly calculate the equations of motion for fluctuations on a solution, given by Eqs. 90 and 91. Given the action polynomial, we can compute the coefficients in the perturbation equations directly, using the following formulae, valid for 1 m + n 4 As it is much easier to work at the level of the action polynomial, these formulae will prove invaluable to anyone wishing to study the phenomenology of maximally symmetric solutions and their fluctuations.

VI. DISCUSSION
We are now at the half way stage in our study of multiple galileon theory. We will continue with a detailed analysis of phenomenology in our companion paper [41], but for now let us take stock of where we are. We have generalised the original galileon theory [24] to include any number of coupled galileon fields. We have focused our attention on the case of two galileon fields, one of which is coupled to the trace of the energy-momentum tensor. We call this version of the theory, bi-galileon theory. Although the two galileons are coupled, perhaps even at quadratic order, we neglect any direct coupling between them and the graviton. This amounts to neglecting the backreaction of the scalars to background geometry.
Through the coupling to the energy-momentum tensor, the galileon fields lead to a long range modification of the gravitational force. We have argued that co-dimension two braneworld models with an infra-red modification of gravity will typically admit a decoupling limit corresponding to a bi-galileon theory. Indeed, we have shown explicitly how the cascading cosmology model [31] does exactly that. Of course, strictly speaking the cascading cosmology model is a 5D model with co-dimension one branes [31], although it is expected to capture many of the salient features of the full 6D cascading gravity model [29,30] which does contain co-dimension two branes. Our interest in galileon modifications of gravity need not rely on higher dimensional models of gravity. We can, in principle, reconstruct our bi-galileon theory along the lines of [27], and develop a fully covariant model of two scalar fields coupled to gravity and to matter in four dimensions. While the Galilean invariance will be broken by this procedure, we expect it to be restored in the limit M pl → ∞, so we can still learn plenty about the model simply by looking at its bi-galileon limit.
Our construction of the most general bi-galileon theory is very similar to the single galileon case, with a straightforward extension to N galileons. As was the case for the single galileon [24], there are only a finite number of terms you can write down in the equations of motion. This is because requiring Galilean symmetry, as well as no higher than second derivatives, forces us to build our equations of motion from scalar products of ∂ µ ∂ ν π and ∂ µ ∂ ν ξ. These can be thought of as matrices from which we can only build a finite number of Cayley invariants.
As a first step towards a study of the phenomenology, we have considered maximally symmetric vacua and their stability. We have found a very elegant way of describing the physics in terms of a particular function, the action polynomial, which is closely related to the on-shell background Lagrangian. Stable vacua correspond to the minima of this polynomial. Computing its Hessian on a particular solution offers us a quick and easy stability test which we will make use of in the future [41]. We can also apply Morse theory to understand the generic spread of stable and unstable vacua. For example, it becomes almost trivial to see that if there is more than one vacuum, they cannot all be stable.
Having motivated and constructed our general theory we conclude by referring the reader to our forthcoming companion paper [41] for a detailed discussion of phenomenology.
Note added: Whilst this paper was in its very final stages of preparation, [42] appeared, in which the galileon model is generalised to mixed combinations of arbitrary p-forms. It is clear that the results of section IV of this paper can be derived from the formulae given in [42], for the case of multiple scalars. Both this work and our own might be considered to be the natural generalisations of [43].