Cosmology of the Galileon extension of Bekenstein's theory of relativistic Modified Newtonian Dynamics

TL;DR

This paper investigates a relativistic modified gravity model called g-TeVeS, which aims to replicate MOND phenomena without dark matter, but finds it struggles to match standard cosmological behavior in certain parameter regimes.

Contribution

It analyzes the cosmological behavior of the g-TeVeS model and identifies parameter space restrictions affecting its ability to reproduce standard cosmology.

Findings

g-TeVeS cannot fully replicate Friedmannian cosmology in its estimated viable parameter space.

Closer agreement with standard cosmology is possible in other regions of the parameter space.

The model employs a Galileon-induced Vainshtein mechanism to suppress modifications in strong gravity regimes.

Abstract

A generalization of Bekenstein's Tensor-Vector-Scalar (TeVeS) model of modified gravity has recently been proposed as an alternative to dark matter. This model -- which we will refer to as g-TeVeS -- utilizes a Galileon-induced Vainshtein mechanism to suppress modifications to General Relativity in strong gravity regimes and so avoids the need to introduce the baroque kinetic terms that typically exist in relativistic models of Modified Newtonian Dynamics (MOND). We explore the behavior of this model in spacetimes with exact Friedmann-Robertson-Walker (FRW) symmetry. The ability of the theory to recover MOND phenomenology places restrictions on the theory's parameter space and it is found that within an estimate of this area of parameter space the theory cannot successfully approximate the Friedmannian cosmological behavior of cold dark matter. It is found that much closer agreement may be recovered in other regions of the theory's parameter space and the reasons for this are discussed.

Figures

Tables

Table 1. Table 1: Best-fit values for various models, quoted to two significant figures: Model A: Restricted parameters with ρ c = 0 subscript 𝜌 𝑐 0 \rho_{c}=0 ; Model B: Restricted parameters with ρ c ≠ 0 subscript 𝜌 𝑐 0 \rho_{c}\neq 0 ; Model C: Unrestricted parameters with ρ c = 0 subscript 𝜌 𝑐 0 \rho_{c}=0 ; Model D: Unrestricted parameters with ρ c ≠ 0 subscript 𝜌 𝑐 0 \rho_{c}\neq 0 .

	$𝒮$	${ϵ}_c$	$1/{\tilde{a}}_0$	$8k_{\mathrm{gal}}$	$\Lambda$	${\rho }_c/{\rho }_b$	${\varphi }_i$	$y_i$	$G_C/G-1$
$A$	$7.0\times {10}^{-2}$	$4.7\times {10}^{-2}$	$5.4$	$1.6\times {10}^{-24}$	$2.3$	$0.0$	$6.0\times {10}^{-2}$	$1.2\times {10}^{-8}$	$6.0\times {10}^{-2}$
$B$	$3.9\times {10}^{-5}$	$1.6\times {10}^{-13}$	$5.4$	$1.0\times {10}^{-25}$	$1.6$	$5.2$	$1.4$	$5.6\times {10}^{-9}$	$4.8\times {10}^{-3}$
$C$	$1.7\times {10}^{-3}$	$2.7\times {10}^1$	$1.7\times {10}^{-4}$	$-7.2\times {10}^{-39}$	$1.3\times {10}^{-1}$	$0.0$	$8.9\times {10}^{-3}$	$-3.1\times {10}^{-5}$	$-1.9\times {10}^{-2}$
$D$	$9.9\times {10}^{-7}$	$2.3\times {10}^{13}$	$8.0\times {10}^{-1}$	$-6.7\times {10}^{-35}$	$2.1$	$5.4$	$2.4\times {10}^{-5}$	$-1.1\times {10}^{-8}$	$-2.1\times {10}^{-4}$

Equations

$$\vec{\nabla }\cdot \left[{\mu }_m(x)\vec{\nabla }\Phi \right]=4\pi G_N{\rho }_b$$

$$S$$

$$S_{\mathrm{grav}}$$

$$S_A$$

$$S_{\varphi }$$

$$+\frac{2k_{\mathrm{gal}}}{3}{\widetilde{ϵ}}^{\alpha \beta \gamma \delta }{\widetilde{ϵ}}^{\mu \nu \rho \sigma }{\tilde{\nabla }}_{\alpha }\varphi {\tilde{\nabla }}_{\mu }\varphi {\tilde{\nabla }}_{\nu }{\tilde{\nabla }}_{\beta }\varphi {\tilde{R}}_{\gamma \delta \rho \sigma }+m_{\varphi }^2{\varphi }^2)$$

$$K^{\alpha \beta \mu \nu }$$

$$+c_3{\tilde{g}}^{\alpha \nu }{\tilde{g}}^{\beta \mu }+c_4{A}^{\alpha }{A}^{\beta }{\tilde{g}}^{\mu \nu })$$

$$g_{\mu \nu }\equiv e^{-2\varphi }{\tilde{g}}_{\mu \nu }-2\sinh (2\varphi )A_{\mu }{A}_{\nu }$$

$$S_{\varphi }^{(\mathrm{TeVeS})}=-\frac{1}{8\pi G}\int d^{4}x\sqrt{-\tilde{g}}f(X)$$

$$S_{\Lambda }=-\frac{1}{8\pi G}\int d^{4}x\sqrt{-g}\Lambda$$

$$\mathbf{g}=-\left(\vec{\nabla }{\Phi }_N+\vec{\nabla }\varphi \right)$$

$$\frac{d\varphi }{dr}={\left(\sqrt{\frac{8k_{\mathrm{gal}}}{r^2}+\frac{r^2}{G_{N}M{a}_0}+{\left(\frac{{ϵ}_{\mathrm{c}}{r}^2}{2G_{N}M}\right)}^2}+\frac{{ϵ}_{\mathrm{c}}{r}^2}{2G_{N}M}\right)}^{-1}$$

$$G$$

$$r_V$$

$$r_M$$

$$a_0\sim 1.2\times 10^{-10}m/s^2\sim c{H}_0/6$$

$$r_{SOL}\sim 10^9{G}_N{M}_{\text{astrosun}}$$

$$k_{\mathrm{gal}}\gtrsim \frac{10^{36}(G_N{M}_{\text{astrosun}}{)}^3}{8a_0}$$

$$k_{\mathrm{gal}}<\frac{10^3{G}_N{M}_{\text{astrosun}}}{8a_0^3}$$

$${\left(2\times 10^{-9}\frac{1}{a_0}\right)}^4\lesssim k_{\mathrm{gal}}<{\left(4\times 10^{-6}\frac{1}{a_0}\right)}^4$$

$${ϵ}_{\mathrm{c}}<0.1$$

$${\tilde{g}}_{\mu \nu }d{x}^{\mu }\otimes d{x}^{\nu }=-{\tilde{N}}^2(t)d{t}^2+b^2(t){\gamma }_{ij}d{x}^i\otimes d{x}^j$$

$${\gamma }_{ij}d{x}^i\otimes d{x}^j=\frac{1}{1-\kappa r^2}d{r}^2+r^{2}d{S}_{(2)}^2$$

$$A_{\mu }=\mp \tilde{N}{\delta }_{\mu }^0.$$

$$g_{\mu \nu }d{x}^{\mu }\otimes d{x}^{\nu }=-N^2(t)d{t}^2+a^2(t){\gamma }_{ij}d{x}^i\otimes d{x}^j$$

$$N$$

$$a$$

$$s=\ln b$$

$$\tilde{h}$$

$$y$$

$$S$$

$$\displaystyle+8k_{\rm gal}\frac{b}{\tilde{N}^{3}}y^{3}\tilde{h}\left(\kappa+\frac{b^{2}\tilde{h}^{2}}{\tilde{N}^{2}}\right)+\pi_{s}\left(\frac{ds}{dt}-\tilde{h}\right)+\pi_{\phi}\left(\frac{d\phi}{dt}-y\right)\bigg{\}}+S_{m}$$

$$\delta S_m$$

Full text

We explore the behavior of this model in spacetimes with exact Friedmann-Robertson-Walker (FRW) symmetry. The ability of the theory to recover MOND phenomenology places restrictions on the theory’s parameter space and it is found that within an estimate of this area of parameter space the theory cannot successfully approximate the Friedmannian cosmological behavior of cold dark matter. It is found that much closer agreement may be recovered in other regions of the theory’s parameter space and the reasons for this are discussed.

.1 Introduction

Milgrom’s observation Milgrom (1983) that a wide variety of the astrophysical phenomena usually attributed to the effects of dark matter can instead be accounted for by a modification to the dynamics of visible matter has provided an intriguing hint that something may be missing in our understanding of gravity and/or inertia. In its original formulation, Milgrom’s Modified Newtonian Dynamics (MOND) was non-relativistic in the same sense that Newton’s theory of gravity is. If Newtonian gravity is a limiting form of General Relativity, what is MOND a limiting form of? One formulation of MOND is as a modified Poisson equation:

[TABLE]

where $\Phi$ is the gravitational potential felt by non-relativistic test particles, $x=|\vec{\nabla}\Phi|/a_{0}$, $\rho_{b}$ is the density of baryonic matter, $G_{N}$ is the locally measured value of Newton’s gravitational constant, $a_{0}$ is a constant with the dimensions of acceleration and $\mu_{m}(x)$ is a function subject to the limiting forms $\mu_{m}\rightarrow 1$ as $x\gg 1$ and $\mu_{m}\rightarrow x$ as $x\ll 1$ but is otherwise unspecified; an explicit form such as $\mu_{m}(x)=x/(1+x)$ has usually been chosen but such forms lack theoretical motivation 111There do, however, exist examples where non-relativistic MOND phenomenology is recovered with a specific counterpart of $\mu_{m}(x)$, whose form is fixed by heuristic theoretical considerations Milgrom (1999); Klinkhamer and Kopp (2011); Pazy and Argaman (2012)..

A significant amount of research has gone into studying the consequences of (1) in astrophysical systems Sanders and McGaugh (2002); Bekenstein and Magueijo (2006); Gentile et al. (2007); Wu et al. (2009); Dai et al. (2010); Zhao and Famaey (2010); Magueijo and Mozaffari (2012); Famaey and McGaugh (2012); Hees et al. (2016); Margalit and Shaviv (2016); Pereira et al. (2016); Ko (2016); however, the lack of a fully-relativistic formulation of the theory makes it difficult to know the realm of the equation’s validity. It was also attempted to derive the MOND formula from fundamental theory Verlinde (2016). A number of relativistic theories that recover MOND-like phenomenology have been proposed Bekenstein and Milgrom (1984); Bekenstein (2004); Navarro and Van Acoleyen (2006); Sanders (2007); Zlosnik et al. (2007); Milgrom (2009); Blanchet and Marsat (2011); Deffayet et al. (2011a); Woodard (2015); Khoury (2015); Kim et al. (2016). All of these examples possess a similar ambiguity to (1) in that all possess a function in the Lagrangian that must be chosen by hand. This makes it difficult to know exactly what a given theory predicts. Each of the examples involve the introduction of new degrees of freedom into physics and it can be that the line between their interpretation as an additional ‘dark force’ in nature or simply a type of dark matter becomes blurred. Indeed, alternatively there have been a number of attempts to produce results for $\Phi$ similar to that in solutions of (1) by instead evoking a dark matter with exotic dynamics and coupling to matter Blanchet and Le Tiec (2008); Berezhiani and Khoury (2016, 2015).

Recently a model that produces MOND-like phenomenology whilst avoiding the use of an unspecified function has been proposed by Babichev, Deffayet, and Esposito-Farese Babichev et al. (2011). This theory uses Bekenstein’s TeVeS theory as basis Bekenstein (2004); Skordis (2009) but extends it with the addition of a Galileon-type term Deffayet et al. (2011b) and the removal of the free function. The Galileon term leads to the Vainshtein screening Vainshtein (1972) of the force generated by the scalar field around the high-curvature environment of the solar system. We will therefore refer to the current model as g-TeVeS for Galileon extended Tensor-Vector-Scalar theory. We also note that applying the use of screening mechanisms to MOND has also been investigate recently in Burrage et al. (2016) by utilizing the symmetron mechanism Hinterbichler and Khoury (2010).

The action for this theory is as follows:

[TABLE]

[TABLE]

where $\tilde{R}$ is the Ricci-scalar corresponding to the connection $\tilde{\nabla}_{\mu}$ compatible with the metric $\tilde{g}_{\mu\nu}$ (i.e. $\tilde{\nabla}_{\mu}\tilde{g}_{\alpha\beta}=0$), $\tilde{g}$ is the determinant of $\tilde{g}_{\mu\nu}$ and $\tilde{\epsilon}^{\alpha\beta\gamma\delta}$ is the completely antisymmetric tensor built using it, and $\chi$ represent the matter fields in the universe,

[TABLE]

where $A^{\mu}\equiv\tilde{g}^{\mu\nu}A_{\nu}$, and matter is taken to couple to the metric $g_{\alpha\beta}$:

[TABLE]

By comparison, the action of Bekenstein’s TeVeS theory differs only in the form of $S_{\phi}$ Bekenstein (2004). In TeVeS there is no Galileon term and, in its ‘diagonal frame’ formulation Skordis (2009), the equivalent action to $S_{\phi}$ takes the form:

[TABLE]

where $X\equiv\tilde{g}^{\mu\nu}\tilde{\nabla}_{\mu}\phi\tilde{\nabla}_{\nu}\phi$ and- as in the case of the function $\mu_{m}$ in (1)- the function $f$ is required to take certain limiting forms but it is otherwise unspecified; that there are no free functions present in (5) is what represents a simplification over TeVeS.

We now discuss the individual terms in (2) in more detail. The first term, $S_{{\rm grav}}$ is simply the Einstein-Hilbert term for the metric $\tilde{g}_{\mu\nu}$ along with a cosmological constant term. One may additionally consider a cosmological term with respect to the metric $g_{\mu\nu}$:

[TABLE]

The action $S_{A}$ is the Einstein-Aether action Jacobson and Mattingly (2001). This is the most general action for a fixed-norm timelike vector field coupled to $\tilde{g}_{\mu\nu}$ which produces field equations second-order in time. The first term in $S_{\phi}$, proportional to the constant $\epsilon_{\rm c}$ is a canonical kinetic term; it is included as a step towards providing a well-defined Cauchy problem for the theory in the limit of weak fields Bruneton (2007). The second term, proportional to the constant $1/\tilde{a}_{0}$ is the term that provides an analog to the limit where $\mu_{m}(x)\sim x$ in (1); it is this limit which is responsible for a dark matter-like effect in astrophysical systems. The third term is is an example of the ‘Paul’ Galileon Lagrangian of the Fab Four scalar field models Charmousis et al. (2012). Its role is to suppress the additional force on test bodies in high curvature regions that would otherwise exist due to the $1/\tilde{a}_{0}$ term. These effects will be discussed in more detail in Section .2. Finally, the term proportional $m_{\phi}^{2}$ is a mass term for the scalar field.

.2 The parameter space

We now briefly discuss the restrictions on the parameters of the g-TeVeS theory that were derived by the authors of Babichev et al. (2011). It is found in the theory that on galactic and sub-galactic scales, the non-relativistic acceleration $\mathbf{g}$ felt by a test body obeys the relation

[TABLE]

where $\vec{\nabla}$ is the spatial gradient operator and $\Phi_{N}$ is the Newtonian potential due to baryonic matter; thus, spatial gradients of the field $\phi$ provide an extra force acting on test bodies and the model is constructed to make that force resemble the phenomenology implied by MOND. It is found in Babichev et al. (2011) that in spherical symmetry, outside of a gravitating object with mass $M$, $d\phi/dr$ takes the form:

[TABLE]

Importantly, it has been assumed that the following equalities hold:

[TABLE]

i.e. that the parameters $\{G,\tilde{a}_{0}\}$ that appear in the action (2) are to be identified respectively with the locally measured value of Newton’s constant $G_{N}$ and the scale $a_{0}$ appearing in MOND phenomenology. This is a non-trivial assumption. It is known in the case of TeVeS that $G_{N}$ and $a_{0}$ will generally depend on a number of additional quantities such as the ‘background’ value of the scalar field $\phi$ (in the case of $a_{0}$), certain constants in the tensor $K^{\alpha\beta\mu\nu}$ should they be non-zero, and constants in the scalar field action Bekenstein (2004); Bekenstein and Sagi (2008); Skordis (2009). In this work we also assume that the equalities (12) hold though a more comprehensive analysis involving an extended version of the Parameterized Post-Newtonian Vanshteinian (PPNV) formalism Avilez-Lopez et al. (2015) may show otherwise (as is the case of TeVeS).

The authors of Babichev et al. (2011) identify two distinct length scales $r_{V}$ and $r_{M}$ $\gg r_{V}$ which mark transitions between forms of $d\phi/dr$:

[TABLE]

The parameters should be such that for regions in the vicinity of gravitating sources where there is little evidence for dark matter, then $|\vec{\nabla}\phi|/|\vec{\nabla}\Phi_{N}|\ll 1$. For example, when considering the gravity due to the sun, if scales $r$ within the solar system are less than $r_{V}$ then the additional force due to $\phi$ is proportional to $r$ and actually decreases the closer one gets to the sun. This is the Vainshtein mechanism created by the presence of the Galileon term.

It is required that as one moves to larger distances from a gravitating source then one reaches scales $r>r_{V}$ (e.g. when considering the motion of stars towards the outer regions of galaxies) then the additional acceleration due to $\phi$ is approximately $-\sqrt{G_{N}Ma_{0}}/r$ where

[TABLE]

and $H_{0}$ is the measured value of the Hubble constant today. However, this acceleration cannot persist on much smaller astrophysical scales such as those of the Solar System as its presence would result in effects on planetary orbits that are not observed Sanders (2006). If, however, the parameter $k_{\rm gal}$ is non-zero and typical planetary distances from the sun are significantly smaller than the radius $r_{V}$ then the additional acceleration provided by $\phi$ will be of magnitude $r/\sqrt{8k_{\rm gal}}$. Thus, if $k_{\rm gal}$ is sufficiently big then the additional acceleration experienced by objects in the solar system will be currently undetectable.

We assign a rough measure $r_{SOL}$ of the maximum distance from the sun probed by experiment to be of the order of lengthscales associated with Pluto’s orbit:

[TABLE]

and we require that $r_{V}$ for the sun is assumed greater than $r_{SOL}$ so that the anomalous acceleration is of the form $r/\sqrt{8k_{\rm gal}}$ for $r<r_{SOL}$. This gives the condition:

[TABLE]

However, as the ratio $r_{V}/r_{M}$ grows as $k_{\rm gal}^{1/4}$, there is an upper bound on how big $k_{\rm gal}$ can be. It is also observed that systems with baryonic mass as little as $10^{3}M_{\astrosun}$ display mass discrepancy effects so the $r_{V}$ should be smaller than $r_{M}$ for these systems:

[TABLE]

Hence the constant $k_{\rm gal}$ is restricted from both sides:

[TABLE]

The effects of an acceleration dominated by the $1/r$ term have been observed at distances from the gravitating source corresponding to at least $10r_{M}$. In the g-TeVeS model, the MOND regime is exited for $r>r_{M}/\epsilon_{\rm c}$ and so we have:

[TABLE]

In summary then, the value $\tilde{a}_{0}$ is essentially fixed by accounting for the motion of stars within galaxies which display evidence for dark matter; the value $k_{\rm gal}$ is restricted to be small enough so that the Galileon term does not dominate in the rather low baryonic mass systems which still display evidence for dark matter but that nonetheless is big enough to dominate in the Solar System; the value $\epsilon_{\rm c}$ is restricted to be small enough so that the MOND term dominates for a sufficiently great span of scales where there is evidence for dark matter.

Additionally there are the four constants $\{c_{i}\}$ in (4). In the canonical Einstein-Aether theory (the theory which is a limit of (2) as $\phi\rightarrow\mathrm{cst.}=0$), it is known that in the quasistatic weak-field regime that the parameters $c_{1}$ and $c_{4}$ contribute to a modification of the relation between $G_{N}$ and $G$; these and other parameters are involved in deviations of the theory’s parameterized post-Newtonian (PPN) parameters from those of General Relativity Foster and Jacobson (2006). It likely must be assumed then that the $\{c_{i}\}$ in general take values such that they have a negligible effect on PPN parameters, including the ratio $G_{N}/G$. Once again, this would involve an extended version of the PPNV formalism Avilez-Lopez et al. (2015) and is left for a future investigation.

A comprehensive comparison between g-TeVeS and the standard $\Lambda$-cold dark matter cosmological model ($\Lambda\mathrm{CDM}$) would require development of its cosmological perturbation theory and the computation of its predictions for the growth of cosmic structure. In the present paper we take the first steps in this direction and examine how close the expansion history of the universe in g-TeVeS can be to that of $\Lambda\mathrm{CDM}$. We will see whether there exist parameters in the restricted parameter space discussed in this section which lead to similarity to $\Lambda\mathrm{CDM}$ and we will also explore the remainder of the parameter space, it is important to do so because- as argued- it is not yet clear what restrictions the astrophysical limit of the theory places on parameters.

.3 FRW Symmetry

We now derive the field equations of the g-TeVeS model assuming FRW symmetry. In this symmetry the Einstein metric $\tilde{g}_{\mu\nu}$ can be cast into the following form:

[TABLE]

where one can further adopt the coordinates

[TABLE]

where $\kappa$ is the spatial curvature scalar of co-moving spatial hypersurfaces and $dS^{2}_{(2)}$ the line element of a $2$-sphere.

The field $A_{\mu}$, is unit-timelike with respect to the metric $\tilde{g}_{\mu\nu}$ due to the Lagrangian constraint in the action (4). In FRW symmetry this is satisfied by the following ansatz, which we will adopt:

[TABLE]

Without loss of generality we will choose the minus sign option, so that the vector $A^{\mu}\equiv\tilde{g}^{\mu\nu}A_{\nu}$ ‘points’ in the future direction. The scalar field $\phi$ is assumed to depend only on the coordinate $t$. The matter frame metric $g_{\mu\nu}$ then takes the form

[TABLE]

where

[TABLE]

We define new variables in order to perform a Hamiltonian formulation of the action. Rather than the scale factor $b$ we will work with a variable $s$ defined by

[TABLE]

assuming that $b$ does not pass through zero.

We furthermore introduce the variables:

[TABLE]

and enforce the constraints via Lagrangian multiplier terms, $\pi_{s}$ and $\pi_{\phi}$. Using some useful results provided in Appendix .10, the FRW-reduced action can be shown to be:

[TABLE]

where following Skordis (2008) we define $K_{F}\equiv 1+(c_{1}+3c_{2}+c_{3})/2$. A way to explicitly include matter is by considering the definition of the matter action in terms of its variation which leads to the stress-energy tensor; in particular we have that

[TABLE]

We take the matter content to be described a perfect fluid, where for each matter species $T_{\mu\nu}=\rho u_{\mu}u_{\nu}+P(g_{\mu\nu}+u_{\mu}u_{\nu})$ where $u_{\mu}\equiv g_{\mu\nu}u^{\nu}$ and $g_{\mu\nu}u^{\mu}u^{\nu}=-1$. After integration over the trivial spatial integral in the action yields the FRW-reduced form

[TABLE]

where we have assumed that $u^{\mu}$ aligns with the direction of cosmic time basis vector $(\partial_{t})^{\mu}$ i.e. $u^{\mu}=(\frac{1}{N},0,0,0)=(\frac{e^{-\phi}}{\tilde{N}},0,0,0)$ and $u_{\mu}=(-e^{\phi}\tilde{N},0,0,0)$.

By requiring stationarity of (30) under small variations of fields, we can proceed to derive the equations of motion. Varying with respect to $s,\phi,\pi_{s},\pi_{\phi},\tilde{N}$ give Hamilton’s equations and the Friedmann equation/Hamiltonian constraint:

[TABLE]

where henceforth $\rho$ and $P$ refer to the total matter density and pressure, including a potential matter frame cosmological constant $\Lambda$. Finally, varying with respect to $\tilde{h}$ and $y$ gives the constraint equations

[TABLE]

For later use, we introduce the following variables:

[TABLE]

The quantity $H$ is the rate of change of $\ln a$ with respect to the matter frame proper time $\tau$ where $d\tau=Ndt$; this is the matter frame Hubble parameter. The quantity $\tilde{H}$ is the Einstein frame Hubble parameter, expressing the rate of change of $\ln b$ with respect to Einstein frame proper time $\tilde{\tau}$ where $d\tilde{\tau}=\tilde{N}dt$ and $\tilde{Y}$ is the rate of change of $\phi$ with respect to $\tilde{\tau}$. Additionally, we find the following relation between the two Hubble parameters:

[TABLE]

Hence $\tilde{H}^{2}=e^{2\phi}H^{2}+2e^{\phi}H\tilde{Y}+\tilde{Y}^{2}$ so that (37) leads to the matter frame Friedmann equation

[TABLE]

In finding solutions to the equations of motion, for simplicity we will from now on set $\tilde{\Lambda}=0$ and $\kappa=0$ and allow the matter frame cosmological constant $\Lambda$ to be non-zero. Additionally we set $m_{\phi}=0$ as was done in Babichev et al. (2011).

.4 Energy Densities

We now look a bit more closely at notions of the energy density of the scalar field $\phi$. This will later aid comparison with the dark sector of the $\Lambda\mathrm{CDM}$ model. It’s helpful to write the matter frame Friedmann equation (42) in a more familiar form:

[TABLE]

where

[TABLE]

We use the subscript $r$ to denote a radiation component and $\rho_{dust}\equiv\rho_{b}+\rho_{c}$ where $b$ denotes the baryonic contribution and $c$ denotes the contribution due to non-baryonic dust (i.e. a cold dark matter component should one exist). In writing (42) in the form (43) we have entirely separated contributions to $H^{2}$ that depend on $\phi$ and those who do not. Furthermore we can develop a notion of fractional energy density of $\Omega_{I}$ of a cosmological component by writing (43) as

[TABLE]

where the index $I$ is over the various types of matter as well as $\phi$. Following earlier TeVeS literature Bourliot et al. (2007), the Einstein frame Friedmann equation may be written as

[TABLE]

where $\tilde{\varrho}_{r}=\rho_{r}$ etc. for matter species and

[TABLE]

This equation may be written as

[TABLE]

We can rewrite (46) using (41) to yield

[TABLE]

where

[TABLE]

This equation may be written as

[TABLE]

By inspection,

[TABLE]

The measures of energy density due to the scalar field $\rho_{\phi}$ and $\tilde{\varrho}_{\phi}$ may differ considerably. The presence of the Galileon term implies that the signs of these quantities are not positive-definite and $\rho_{\phi}$ is not positive-definite even in its absence.

.5 Numerical Solution Strategy

Aside from a number of specific situations we will look at, it will be necessary to numerically integrate the equations of motion- this will involve the evolution of the set $\{s,\phi,\pi_{s},\pi_{\phi},\tilde{h},y\}$, i.e. a total of six variables. We also have the three constraints (37), (38) and (39), hence, there are three independent variables to be evolved. That is, the space of initial data is three-dimensional. Finally, we have the four evolution equations (28), (29), (33) and (34) for $s,\phi,\pi_{s},\pi_{\phi}$ respectively, but of course not all three are independent. Additionally, a specific form of $\tilde{N}$ must be chosen. For integration, first we choose an initial time $t_{0}$ and set initial conditions for all variables as follows:

Freely specify values $\{s(t_{0}),\phi_{i}\equiv\phi(t_{0}),y_{i}\equiv y(t_{0})\}$ as initial conditions and choose a form of the function $\tilde{N}$. 2. 2.

Use the constraint (37) to find $\tilde{h}(t_{0})$. Of the three roots to this cubic polynomial in $\tilde{h}$, there is one positive root which reduces to the result that would appear in the absence of the Galileon term when we take the limit $k_{\rm gal}\rightarrow 0$ and we choose this root. If, for example, the contribution of the Galileon term to the right hand side of (37) is small and positive, there will be a positive root $\tilde{h}_{2}\gg\tilde{h}_{1}$ representing another solution to the field equations but it will deviate significantly from $\Lambda\mathrm{CDM}$ and will be neglected. 3. 3.

Use the constraint (38) to determine $\pi_{s}(t_{0})$. 4. 4.

Use the constraint (39) to determine $\pi_{\phi}(t_{0})$.

At this point, all initial conditions are consistently set. We proceed to integrate the four equations (28), (29), (33) and (34) as follows. For each time $t_{n}$, repeat the following steps:

Consider a nearby moment $t_{n+1}>t_{n}$. If $t_{n+1}<t_{\mathrm{end}}$ then evolve to a nearby moment $t=t_{n+1}$ to determine $\{s(t_{n+1}),\phi(t_{n+1}),\pi_{s}(t_{n+1}),\pi_{\phi}(t_{n+1})\}$ using (28), (29), (33) and (34). 2. 2.

Use the constraints (38) and (39) to determine $\tilde{h}(t_{n+1})$ and $y(t_{n+1})$. This involves solving a pair of coupled equations cubic in $\tilde{h}$ and $y$. Care must be taken that a different branch of $\{\tilde{h},y\}$ solutions is not found here (much like how it was appropriate to find $\tilde{h}_{1}$ rather than $\tilde{h}_{2}$ as a root in the initial value equation for $\tilde{h}$). 3. 3.

Check that the constraint (37) is satisfied to reasonable numerical precision. This should guarantee that the previous step has avoided skipping onto a different branch of solutions for $\{\tilde{h},y\}$.

.6 Regimes and Approximate Solutions

We now consider the behavior of the scalar field and metric in some specific situations. For ease of calculation we’ll use a spacetime gauge where $\tilde{N}=1$ and so we identify $t=\tilde{\tau}$ where $\tilde{\tau}$ measures Einstein frame cosmic proper time. The evolution equation for $\pi_{\phi}$ (34) can be integrated to yield:

[TABLE]

where again $\rho$ and $P$ refer to the total matter densities and pressures and $\pi_{\phi(0)}$ is a constant. From the $\pi_{\phi}$ constraint (39) and defining $\sigma_{Y}={\rm Sign}[\tilde{Y}]$, we then have:

[TABLE]

where

[TABLE]

Solutions to (55) are the cosmological analog of the solution (11), in other words, we are dealing with a temporal Vainshtein mechanism.

We now identify two main regimes of evolution. When

[TABLE]

we shall refer to this as the Canonical Regime as here the canonical kinetic term is dominating the evolution of the field $\phi$. When

[TABLE]

we will refer to this as the Nonlinear Regime. The transition between the linear and nonlinear regimes occurs when $\tilde{Y}_{lin}\sim\tilde{Y}_{nonlin}$ i.e. when $|{\cal C}{\cal A}|\sim 1$ and the Nonlinear Regime corresponds to $|{\cal C}{\cal A}|\gg 1$.

Within the Nonlinear Regime limit, there exist two further distinct sub-regimes, depending on the dominant terms inside ${\cal A}$. The first of these is when

[TABLE]

In this limit the MOND term contribution to ${\cal A}$ dominates and hence we will refer this to the MOND Regime.

Alternatively, when

[TABLE]

then ${\cal A}(\tilde{\tau})$ is dominated by the Galileon term and so we refer to this as the Galileon Regime.

The g-TeVeS parameter restrictions of Section .2 enforce that ${\cal A}(\tilde{\tau})$ is positive-definite during the Galileon era. Therefore, if these restrictions are adopted then the solution (59) only possesses real solutions if ${\cal C}$ is negative. By inspection $J$ is negative-definite for matter sources with equation of state $w\geq-1$ and so the function ${\cal C}$ may only be non-positive with a suitably chosen $\pi_{\phi(0)}$.

Consider two universes with identical ${\cal C}(t)$ and with parameters enabling similar evolution of the metric but with one with a sufficiently large ${\cal A}(t)$ to push it into the nonlinear regime. The ratio of $\tilde{Y}$ in the universe where it is in the nonlinear regime to that of a universe where it is in the linear regime is:

[TABLE]

As $|{\cal A}_{{\rm nonlin}}{\cal C}|\gg 1$, then the ratio (62) is below unity so that $\tilde{Y}_{{\rm nonlin}}$ is suppressed compared to $\tilde{Y}_{lin}$. This is a realization of a temporal Vainshtein mechanism. In particular as $b\rightarrow 0$, then $\tilde{H}\rightarrow\infty$ so that $\tilde{Y}\rightarrow 0$ and the theory reduces to GR in the early universe.

We now discuss some approximate solutions to the theory; these will aid the interpretation of numerical results.

.6.1 Canonical Regime with matter fluid with constant $w_{mat}$

First we consider the limit where the scalar field is inside the Canonical Regime and there is a single matter source with constant equation of state $w_{mat}\equiv P_{mat}/\rho_{mat}$. Consider the following ansatz:

[TABLE]

Then we have $\tilde{Y}=n\tilde{H}$ and it is found that the evolution equations for $\pi_{\phi}$ and $\pi_{s}$ subject to the Hamiltonian constraint are consistent if

[TABLE]

which yields an Einstein frame Friedmann equation:

[TABLE]

It follows then that

[TABLE]

The interpretation of this is that the scalar field tracks the matter component so that $\tilde{\varrho}_{\phi}/\rho_{mat}$ is a constant. This situation is identical to ordinary TeVeS Skordis et al. (2006); Bourliot et al. (2007).

Recall now that $H=e^{-\phi}(\tilde{H}-\tilde{Y})=e^{-\phi}(1-n)\tilde{H}$ and so $\tilde{H}^{2}=e^{2\phi}H^{2}/(1-n)^{2}$, and also $b=e^{\phi}a$, and so

[TABLE]

where we have defined an effective Newton’s constant $G_{eff}$:

[TABLE]

Note that $G_{eff}$ is not positive-definite and will generally have a time dependence via $\phi$. In the dust-dominated regime ($w_{mat}=0$) and cosmological constant-dominated regime ($w_{mat}=-1$), the effective Newton’s constant in the Canonical regime is negative for $0<\epsilon_{c}/K_{F}<1/3$ and this should be excluded.

We see from (68) that in the Canonical Regime, the ratio $\epsilon_{\rm c}/K_{F}$ controls two important effects. The first is a constant rescaling of the effective Newton’s constant in the matter frame Friedmann equation (i.e. the explicit dependence of $G_{eff}$ on $n$), the second is a time-dependent rescaling of the effective Newton’s constant via the factor $e^{-4\phi}$. From (63) we see that $d\phi/d\ln b$ is proportional to $n$ and so a larger value of $K_{F}/\epsilon_{\rm c}$ will lead to a more dramatic growth in $\phi$.

.6.2 General Regime with negligible gravitation of scalar field kinetic terms

Again we consider a dominant gravitating matter source with constant $w_{mat}$. Now assume that there is a situation where $\tilde{Y}$ has a negligible impact on the Einstein-frame Friedmann equation (37) and evolution equation (33) for $s$ i.e. that

[TABLE]

where $\rho_{w}\equiv\rho_{(w)0}/a^{3(1+w_{mat})}$. Moreover (34) gives

[TABLE]

The above equations suggest a linear relation between $\pi_{\phi}$ and $\pi_{s}$ so that we may take the ansatz $\pi_{\phi}=\beta\pi_{s}+C_{\pi}$ – where $C_{\pi}$ is a constant– and comparing (69) and (71) we have that

[TABLE]

while from (38) and (39) and having in mind that $|K_{F}\tilde{H}|\gg 4|k_{\rm gal}\tilde{Y}^{3}\tilde{H}^{2}|$

[TABLE]

Given this, in the Canonical Regime for $C_{\pi}=0$ we have \tilde{Y}\mathrel{\vbox{ \offinterlineskip\halign{\hfil#\cr\propto\cr\kern 2.0pt\cr\sim\cr\kern-2.0pt\cr}}}\tilde{H} which is consistent with the results of Section .6.1.

In the non-linear regime we find the following behaviour:

MOND regime

In the MOND Regime we find two distinct behaviors which depend on the size of the constant $C_{\pi}$. If $3|\beta K_{F}\tilde{H}|\gg|C_{\pi}|b^{-3}$ and $\tilde{H}>0$ we have $\tilde{Y}\simeq-{\rm Sign}(\beta)\sqrt{3\tilde{a}_{0}|\beta|K_{F}}\tilde{H}^{1/2}$. In the radiation era this corresponds to $\tilde{Y}\sim 1/a$ and $\phi\sim C_{1}+C_{2}a$ while in the matter era it corresponds to $\tilde{Y}\sim a^{-3/4}$ and $\phi\sim C_{1}+C_{2}a^{3/4}$ 222The constants $C_{1}$ and $C_{2}$ are of course not the same in all of the considered cases.. In both cases, the evolution of $\phi$ is negligible with respect to the Friedman equation where $\phi$ can be taken to be approximately constant.

If $3|\beta K_{F}\tilde{H}|\ll|C_{\pi}|b^{-3}$ then $\tilde{Y}\sim a^{-3/2}$ irrespective of the equation of state of matter. However, the evolution of $\phi$ in the radiation era is $\phi\sim C_{1}+C_{2}\sqrt{a}$ while in the matter era $\phi\sim C_{1}+C_{2}\ln a$. The later, is akin to the tracker solutions found in section .6.1.

It is easy to check that both types of behaviors are consistent and may also be derived using (56) and (57) into (59) under the a posteriori justified assumption $\phi\sim C_{1}$.

Galileon regime

In the Galileon Regime we find once again two distinct behaviors which depend on the size of the constant $C_{\pi}$. If for $3|\beta K_{F}\tilde{H}|\gg|C_{\pi}|b^{-3}$ we have $\tilde{Y}^{2}\tilde{H}^{2}\simeq-\beta K_{F}/(4k_{\rm gal})$. Clearly, this case can only be satisfied for $k_{\rm gal}<0$. Furthermore, during the radiation era $\tilde{Y}\sim a^{2}$ and $\phi\sim C_{1}+C_{2}a^{4}$ while in the matter era $\tilde{Y}\sim a^{3/2}$ and $\phi\sim C_{1}+C_{2}a^{3}$.

If $3|\beta K_{F}\tilde{H}|\ll|C_{\pi}|b^{-3}$ we have that $\tilde{Y}\sim a^{3/2}$ and $\phi\sim C_{1}+C_{2}a^{7/2}$ while in the matter era $\tilde{Y}\sim a^{3/4}$ and $\phi\sim C_{1}+C_{2}a^{9/4}$.

In all cases, the constants $C_{2}$ multiplying the time-dependent part of $\phi$ are tiny and are in addition multiplied by tiny values of the scale factor, so that $\phi$ may in fact be taken to be constant during the galileon era.

.7 Comparison To $\Lambda\mathrm{CDM}$

A given set of parameters and initial conditions of the model will produce an expansion history $H(a)$. The behavior of this function will reveal the extent to which the model can imitate the effect of dark matter on the cosmological background expansion. A simple measure of the extent to which this agrees with the $\Lambda\mathrm{CDM}$ ‘concordance model’ expansion history $H_{\Lambda\mathrm{CDM}}(a)$ is via the quantity

[TABLE]

where $a_{i}$ is the physical scale factor at the beginning of integration of the equations of motion (taken to be deep in the radiation era), $a_{n}$ is its value today and $H_{\Lambda CDM}(a)$ is determined by General Relativity’s Friedmann equation

[TABLE]

where spatial curvature has been assumed to be zero. For the fiducial $\Lambda\mathrm{CDM}$ model to which we compare to, we use the following cosmological parameters:

[TABLE]

The $\Lambda\mathrm{CDM}$ Hubble parameter $H_{0}$ today is taken to be $67\mathrm{km}\mathrm{s}^{-1}\mathrm{Mpc}^{-1}$. We use units where this value is equal to unity and hence, in particular, we have $H_{0}/a_{0}\simeq 5.44$. We take the scale factor today $a_{n}$ to also be equal to unity and in looking at g-TeVeS models we take $\rho_{b}$ and $\rho_{r}$ in general to be assumed known and identical to their corresponding values in the $\Lambda\mathrm{CDM}$ case.

Exploration of how ${\cal S}$ varies across the parameter space of the k-mouflage model will help us understand what features of the model are responsible for producing a similar or disimilar cosmology to that of $\Lambda\mathrm{CDM}$. The parameters of the g-TeVeS model cosmological background evolution equations are as follows:

$\epsilon_{\rm c}$: the canonical kinetic term coefficient 2. 2.

$\tilde{a}_{0}$: the MOND kinetic term coefficient 3. 3.

$k_{\rm gal}$: the Galileon kinetic term coefficient 4. 4.

$K_{F}$: Einstein-aether field constant. 5. 5.

$\phi_{i}(a_{i})$: the initial value of the scalar field at the start of integration 6. 6.

$y_{i}(a_{i})$: the coordinate time derivative $d\phi/dt$ at the start of integration.

Of the matter content, we allow $\Lambda$ and- in some cases- $\rho_{c}$ to vary and for simplicity we assume that spatial curvature $\kappa$ is negligible. Consider the matter frame Friedmann equation at some early time deep in the radiation era:

[TABLE]

where the dots denote assumed subdominant contributions due to other species of matter and time derivatives of $\phi$. If the Friedmann equation takes this form across a span of time where $\tilde{Y}\ll H$ deep into the radiation era then the expansion of the universe will be approximately that of General Relativity with a value of Newton’s constant $G_{C}=Ge^{-4\phi_{i}}/K_{F}$. This implies though that $G_{C}$ will not necessarily be the same as the locally measured value of Newton’s constant $G_{N}$ which itself is assumed in Babichev et al. (2011) to be equal to $G$. We will continue this assumption in the optimization procedure. Constraints on the primordial helium abundance from Big Bang Nucleosynthesis (BBN) Carroll and Lim (2004) constrain the difference between these two quantities to be in the following interval:

[TABLE]

This should be regarded as a conservative constraint; it is likely that using restrictions from BBN in conjunction with data from the cosmic microwave background (CMB) will be able to place much more stringent restrictions on $G_{C}/G_{N}$ Avilez and Skordis (2014). If the above assumptions hold then the above bounds constrain the deviation of $e^{-4\phi_{i}}/K_{F}$ from unity and this is implemented in the procedure to try and find minima of the quantity (74). Instead of varying the pair $\{\phi_{i},K_{F}\}$ we equivalently vary $\{\phi_{i},G_{C}\}$.

The final quantity that we will consider is the (assumed dust-like) dark matter density $\rho_{c}$ (or equivalently the constant ratio $\rho_{c}/\rho_{b}$). For a modified gravity alternative to dark matter, it may be hoped that the new gravitational degrees of freedom of the theory may producde a dark matter effect by themselves; failing that, it is instructive to see how much dark matter may be needed in addition (not necessarily the amount present in $\Lambda\mathrm{CDM}$).

We consider the following four models:

• •

Model A: Minimization of ${\cal S}$ given the restrictions on the g-TeVeS theory’s parameters described in Section .2 and the restriction that $\rho_{c}=0$.

• •

Model B: Minimization of ${\cal S}$ given the restrictions on the g-TeVeS theory’s parameters described in .2 and an unrestricted $\rho_{c}$

• •

Model C: Minimization of ${\cal S}$ without any restrictions on the g-TeVeS theory’s parameters and the restriction that $\rho_{c}=0$

• •

Model D: Minimization of ${\cal S}$ without any restrictions on the g-TeVeS theory’s parameters and an unrestricted $\rho_{c}$.

To minimize ${\cal S}$ we employ a Markov Chain Monte-Carlo (MCMC) procedure where random initial values for $V=\{\epsilon_{\rm c},\tilde{a}_{0},k_{\rm gal},G_{C},\Lambda,\rho_{c}/\rho_{b},\phi_{i},y_{i}\}$ are chosen subject to the restrictions placed upon the particular model then new values are generated and adopted if they lead to a decrease in ${\cal S}$.

The best-fit parameters for each model are contained in Table 1 and for each of the best-fits given, the ratio of their predicted $H(a)$ vs physical scale factor $a$ is plotted in Figure 1. We additionally show the corresponding plot for the General Relativistic case without Dark Matter and with $\Omega_{\Lambda}=1-\Omega_{b}-\Omega_{r}$ and this corresponds to a value ${\cal S}=0.11$. We stress that we make no comparison with data and that these “best-fits” are with respect to closeness to our fiducial $\Lambda\mathrm{CDM}$ model. We now discuss the features of the best-fit for each model:

Model A

The best-fit value of ${\cal S}$ for this model is $0.07$, and so in the sense of (74) the Hubble parameter is somewhat closer to $\Lambda\mathrm{CDM}$ than the General Relativistic model without dark matter. However, as can be seen from Figure 1 , the best-fit model has a Hubble parameter substantially smaller than that of $\Lambda\mathrm{CDM}$ for much of the cosmological matter era, meaning that the scalar field $\phi$ is unable to provide an additional dust-like dark matter component during this era.

It may be verified that for the entirety of the evolution of $\phi$, the field $\tilde{Y}$ is in the Nonlinear Regime with $|\tilde{Y}|\gg|{\cal C}{\cal A}|$ (see Section .6), being in the Galileon Regime until around $a\sim 10^{-5}$ before transitioning to the MOND regime; the feature where the $H/H_{\Lambda CDM}$ momentarily is greater than unity shortly before $a\sim 1$ corresponds to $\tilde{Y}$ changing sign but quickly returning to the MOND regime. It was found in Section .6.1 that if the theory is in the Canonical Regime that the effective Newton’s constant $G_{eff}$ is negative for $0<\epsilon_{\rm c}/K_{F}<1/3$; given that $\epsilon_{\rm c}$ is restricted to be less than $0.1$ for Model A, it is possible that this tends to an avoidance of the Canonical Regime. We will see in the case of Model C that the theory being in the Canonical Regime for much of the matter era is vital to producing a dark matter effect.

Model B

This model retains the parameter restrictions of Model A but allows for a non-zero $\rho_{c}$. The best-fit value of ${\cal S}$ for this model is $3.9\times 10^{-5}$, considerably closer to $\Lambda\mathrm{CDM}$ than that of Model A; inspection of Figure 1 shows the corresponding physical Hubble parameter $H(a)$ to be extremely similar to that to the $\Lambda\mathrm{CDM}$ Hubble parameter. From Table 1 we see that a value of $\rho_{c}=5.2\rho_{b}$ is preferred- a contribution of dark matter comparable to that in $\Lambda CDM$.

It may be verified that the behavior of $|\tilde{Y}|$ for this fit is roughly similar to that of Model A and hence the theory remains in the Nonlinear Regime throughout cosmological evolution but interestingly Model B is the one model of the four considered here which appears to favour a value of $\phi_{i}$ of order unity.

Model C

For this model the astrophysical restrictions are not used on the parameter space and there is no additional dark matter. A value of ${\cal S}=1.7\times 10^{-3}$ is found and hence the closeness to the $\Lambda\mathrm{CDM}$ model compared to Model A is rather greater. For this model, the evolution of $\tilde{Y}$ and $\phi$ as a function of physical scale factor $a$ are plotted in Figure 2 as well as the quantity $1/|{\cal A}|$ of (56) which tracks the dominant source of nonlinearity; this quantity increases during the Galileon Regime and is approximately constant during the MOND Regime. We see then from Figure 2 that the theory is initially in the Galileon Regime and during this time $\tilde{Y}$ grows as $1/\tilde{H}\sim a^{2}$ (using the approximate proportionality of scale factors $b$ and $a$ due to the near-constancy of $\phi$ during this era); this evolution is in accordance with the results of (73), suggesting that its assumptions such as the linear relation between scalar field momentum $\pi_{\phi}$ and gravitational momentum $\pi_{s}$ are good ones. If the solution (73) is an accurate approximation then we see that in the radiation era the scalar field in the Galileon Regime acts like an effective cosmological constant in the Einstein frame Friedmann equation (37). At later times during the radiation era, the scalar field transitions to the MOND Regime and Figure 2 shows $\tilde{Y}$ decreasing as $\tilde{H}^{1/2}\sim a^{-1}$ here, again consistent with the results of (73) and whilst these approximations hold then during this regime, terms in $\tilde{Y}$ produce an additional, dust-like contribution to (37).

As the dominant non-scalar field contribution to the background evolution becomes the baryonic dust, $\rho_{\phi}$ will no longer vary as $a^{-3}$ during the MOND Regime within the approximation (73) and the optimization procedure prefers parameters that produce a transition to the Canonical Regime at roughly the time when baryons begin to dominate over radiation; this leads to a preferred value of $\tilde{a}_{0}\sim 5.9\times 10^{3}$ and this value is substantially different to the MOND acceleration scale $a_{0}\simeq 5.44$. However, we emphasize that the relation between $\tilde{a}_{0}$ and $a_{0}$ may be such that they need not be numerically close.

From the solutions of Section .6.1, we see that in the limit of constant $\phi$, the scalar field produces an effective rescaling of Newton’s constant and thus acts like an additional gravitating source whose density varies as $a^{-3}$. Crucially though, $\phi$ is found to vary substantially during the era where the density of baryons is greater than that of radiation (see Figure 2). The effective dark matter-like contribution to $H$ in the matter era is then produced by a time-dependent rescaling of the effective Newton’s constant. From (63) we see that the rate of growth of $\phi$ is determined by the ratio $K_{F}/\epsilon_{\rm c}$; this ratio must take the right value so that the rescaling of Newton’s constant due to $\phi$ is not too big or too small. The best-fit favors a value $\epsilon_{\rm c}/K_{F}\sim 27$. Additionally, the best-fit involves a negative value of $k_{\rm gal}$; we note however that similar but slightly higher values of ${\cal S}$ can be reached with positive $k_{\rm gal}$ and so it is not clear how much cosmology favors one sign of $k_{\rm gal}$ over another. Bear in mind, however, that the astrophysical Vainshtein mechanism operates only for $k_{\rm gal}>0$ and so negative $k_{\rm gal}$ are of dubious importance.

Model D

Finally we consider a model which involves no restrictions on the g-TeVeS parameters and allows for a dark matter component. An excellent fit is achieved, with ${\cal S}=9.9\times 10^{-7}$ and, as with Model B, it can be seen from Figure 1 that the physical Hubble parameter $H(a)$ is extremely close to $H_{\Lambda\mathrm{CDM}}(a)$. The preferred values for $\rho_{c}$ and $\Omega_{\Lambda}$ are close to the values chosen for the $\Lambda\mathrm{CDM}$ cosmology used. This indicates that the modification of gravity due to the fields $\phi$ and $A_{\mu}$ has been suppressed.

It can be verified that for this best-fit, the field $\tilde{Y}$ enters the Canonical Regime deep in the radiation era and never spends time in the MOND Regime; this is in part due to the large value of $\epsilon_{\rm c}\sim 2.3\times 10^{13}$ and this value also leads to a far smaller gravitational contribution of $\tilde{Y}$, as we expect from the results of Section .6.1.

.7.1 The initial data for $\tilde{Y}$ and the existence of pathological solutions

It is found that a wide variety of $\{\epsilon_{\rm c},\tilde{a}_{0},k_{\rm gal},G_{C},\Lambda,\rho_{c}/\rho_{b},\phi_{i},y_{i}\}$ yield universes that do not reach the defined present cosmological moment $a=1$. This is illustrated in Figure (3) which shows the evolution of $|\tilde{Y}|$ for varying initial values of $|\tilde{Y}|$ and for parameters that obey the restrictions of constants of the Model A best-fit. Each of the curves that cease continuing past a particular scale factor $a<1$ represent a cosmos that cannot evolve past that scale factor; this is not due to the expansion rates $H$ or $\tilde{H}$ getting stuck at zero but rather because the equations of motion no longer possess real solutions at the following moments. Curves whose evolution is thwarted are those where $|\tilde{Y}|$- having initially been positive- reaches zero whilst the Galileon term dominates over the MOND term; curves which persist to evolve after reaching [math] are those which reach [math] when the MOND term dominates. Roughly speaking, in the limit of Galileon domination, the evolution equation for $\tilde{Y}$ may be written in the form:

[TABLE]

hence $\tilde{Y}^{2}(t)=\tilde{Y}_{i}^{2}+2\int_{t_{i}}^{t}\xi dt$. If the integral contribution is negative and decreasing as $t$ increases then at some time $t_{0}$, $\tilde{Y}^{2}$ will reach zero and then evolution will cease as $\tilde{Y}^{2}$ cannot continue on to take negative values. If, instead, the MOND term dominates then the evolution equation for $\tilde{Y}$ takes the form

[TABLE]

This does not encounter the same problem as that with (79). Upon reaching [math] following an initial, positive $\tilde{Y}_{i}$, $\tilde{Y}$ can pass to a negative value, for which the solution takes the form $-\tilde{Y}^{2}(t)=2\int_{t_{0}}^{t}\zeta dt$, where $t_{0}$ is the moment that $\tilde{Y}$ reaches [math], and so a continually negative $\zeta(t)$ does not cause evolution to stop.

Even for curves which persist, it can be verified that of them in Figure 3 correspond to universes that eventually begin collapsing. The values of $\{\epsilon_{\rm c},\tilde{a}_{0},k_{\rm gal},G_{C},\Lambda,\rho_{c}/\rho_{b},\phi_{i},z_{i}\}$ are such for these curves that the collapse of the universe for these universes begins after the present moment $a=1$.

As can further be seen from Figure 3, many curves at early times follow a solution $|\tilde{Y}|\sim b^{3/2}$ which is consistent with the approximate solution (73) assuming $\phi\sim\mathrm{cst.}$ (hence $a\propto b$) and $\tilde{H}\sim 1/b^{2}\sim 1/a^{2}$ i.e. that the universe is radiation dominated. For the larger values of $|\tilde{Y}_{i}|$ we see a deviation from this solution before $|\tilde{Y}|$ begins decreasing. The decreasing of $|\tilde{Y}|$ marks the onset of the MOND regime and the deviation from $|\tilde{Y}|\sim a^{3/2}$ marks a new solution where $|\tilde{Y}|\sim a^{1/2}$ and is having a significant impact on the background evolution. It is found that higher yet initial values of $|\tilde{Y}|$ to collapse of the universe before $a=1$. These different behaviors are shown in Figure (4) for a wide range of initial data $\{\phi_{i},\tilde{Y}_{i}\}$ for the Model A best-fit constant parameters. Initial data in the lower (blue) region of the plot all lead to thwarted evolution. Initial data in the upper (green) region correspond to universes where the universe expands for a period but collapses before $a=1$. The upper region is constructed from two, independent possibilities: the region above the (interrupted) upper straight line in Figure (4) denotes parameter space where no real, positive solutions to the Einstein frame Friedmann equation (37) exist; the additional region which splits the middle region where evolution persists to the present scale factor involves solutions where the universe evolves for some time with positive $\tilde{h}$ but where ultimately the physical scale factor $H$ passes through zero to a negative value at a moment earlier than the present. The accompanying, unusual, evolution of $H$ for one such universe is shown explicitly in Figure 5.

Negative initial $\tilde{Y}_{i}$ for Model A parameters are found to generally lead to similar behavior but with a substantially narrower range of parameters leading to universes which persist. For less restricted parameter ranges- particularly those with $k_{\rm gal}<0$- these restrictions on the allowable $\tilde{Y}_{i}$ don’t necessarily still apply; the reason for this is likely that when $k_{\rm gal}<0$ then in the Galileon Regime there exist solutions of (73 )for $\tilde{Y}^{2}$ even when $\xi=0$ - which is not the case for $k_{\rm gal}>0$.

The presence of thwarted evolution and (too early) collapse of universes represents a challenge to the optimization procedure as many universes become simply unsuitable for comparison with $\Lambda CDM$.

.8 Comparison to TeVeS

Having seen FRW solutions of the k-mouflage Galileon model, we now look at how results compare to those of the original TeVeS theory. This will help show what comparative features of g-TeVeS and TeVeS are responsible for dark matter-type effects in cosmology. As was discussed in Section .1, the two theories differ in their form of the scalar field action $S_{\phi}$; an equivalent formulation of TeVeS has the following form of scalar field action:

[TABLE]

where the auxiliary field $\mu$ has been has been introduced. In FRW symmetry we have:

[TABLE]

From this we have contributions to the $s$,$\tilde{N}$,$y$ equations of motion and a new equation of motion-obtained by varying with respect to $\mu$:

[TABLE]

It may then be possible to obtain $\mu(y)$ from this equation. In Bekenstein (2004), Bekenstein proposed a specific choice of $W$ which yielded MOND-like phenomenology on astrophysical scales. This form is:

[TABLE]

where $\hat{\mu}\equiv\mu/\mu_{0}$, where $\mu_{0}$ and $\xi_{B}$ are constant parameters 333The constant $\xi_{B}$ is related to Bekenstein’s constant $\ell_{B}$ and $\mu_{0}$ via $\xi_{B}=\frac{\mu_{0}^{3}}{64\pi\ell_{B}^{2}}$, according to the conventions of  Bourliot et al. (2007).. The FRW cosmological regime spans from $\hat{\mu}=2$ to $+\infty$. Applying the optimization procedure to this model- and again assuming $G=G_{N}$ 444Strictly speaking $G\neq G_{N}$ in TeVeS, but rather $G_{N}/G=2/\mu_{0}+2/(2-c_{1}+c_{4})$. However, as we are interested in comparing to g-TeVeS where the exact relation is at present unknown, we keep using $G=G_{N}$. In any case, the exact relation between the two constants will only become important when placing constraints using data, which is not what we do in the present work.- we find that ${\cal S}$ is minimized for the following parameter values:

[TABLE]

where all dimensionful parameters are expressed in units where $H_{0}=1$.

This is a comparatively good fit to the Model C best-fit unrestricted g-TeVeS model and it can be seen from Figures 6 and 7 that each fit produces similar expansion histories and $\Omega_{I}$, $\tilde{\omega}_{I}$ quantities. It is apparent from Figure 7 that TeVeS reduces to a model dominated by a canonical kinetic term at around between $b=10^{-3}$ and $10^{-2}$ (when the baryons begin to dominate the Einstein frame expansion) due to the tracking of the dominant matter component in agreement with equation (66).. It may be shown that in the Canonical Regime of TeVeS, the number $\mu_{0}/K_{F}$ corresponds to $\epsilon_{\rm c}/K_{F}$; for the TeVeS best-fit and the Model C g-TeVeS best-fit (see Section .7), these quantities take similar values (approximately $25$ and $27$ respectively). Furthermore, the value of $\phi$ at $a=1$ is similar in both cases: $-0.527$ (TeVeS) and $-0.532$ (g-TeVeS).

.9 Conclusions and Discussion

In this paper we have examined the solutions to the g-TeVeS theory in FRW symmetry. The theory has been constructed to replicate the success of MOND which is a modification to Newtonian gravity and as such involves only visibile matter and the Newtonian gravitational field. The relativistic theory, however, involves the introduction of new degrees of freedom in the gravitational sector and these can not just complicate the relation between constants appearing in actions and observed quantities but they can also modify MOND itself. A manifestation of the former is that the numbers $G$ and $a_{0}$ associated with MOND may depend on the local value for the scalar field $\phi$. We have seen that if the theory is to have a dark matter-type effect in the background cosmology then this will likely involve significant growth of $\phi$ over cosmological history. Therefore it may be inaccurate to assume that $|\phi|\ll 1$ on astrophysical scales and this would endanger the assumptions made in Babichev et al. (2011) that the bare Newton’s constant $G$ and acceleration scale $\tilde{a}_{0}$ are to be identified with the measured Newton’s constant and MOND acceleration scale.

It is found that the theory can give a reasonable match to the $\Lambda\mathrm{CDM}$ expansion history for parameters $\{\epsilon_{\rm c},k_{\rm gal},\tilde{a}_{0}\}$ that differ orders of magnitude from the suggested restrictions on them given in Babichev et al. (2011) however it remains an open question whether these values (accompanying the substantial deviation of $\phi$ from [math]) can lead to MOND-like phenomenology.

We can isolate three effects that seem relevant to the theory resembling a cold dark matter component. The first is that when the gravitational effect of $\phi$ is sub-dominant in the radiation era, if the theory is in the MOND Regime, then the scalar field’s kinetic terms approximate an additional dust component; this does not persist into the (baryonic) matter era and so the MOND Regime must be exited and the Canonical Regime reached. In this regime, the growth of $\phi$ is substantial and this can lead to an increase in Newton’s constant. If the growth is of the right amount, then this increase can resemble- by making matter heavier roughly by a factor of the $\Lambda\mathrm{CDM}$ ratio $(\rho_{b}+\rho_{c})/\rho_{b}$- an additional cold dark matter component. If approximation to $\Lambda\mathrm{CDM}$ pushes the theory towards being in the MOND regime during much of the radiation era tends to suppress $k_{\rm gal}$ relative to $\tilde{a}_{0}$; indeed similarly good fits for Model C were possible if $k_{\rm gal}=0$. The cosmological moment of transition between MOND Regime and Canonical Regime is further dictated by $\epsilon_{\rm c}$ and $\tilde{a}_{0}$ so their values must be tuned to find the right moment. The amount of rescaling of Newton’s constant chiefly depends on $\epsilon_{\rm c}/K_{F}$ and so it must be tuned in isolation.

Though it was not possible to find cosmologies similar to that of $\Lambda\mathrm{CDM}$ when parameters were restricted to obey the astrophysical constraints of Babichev et al. (2011), one may wonder whether additional terms may be added to the action (2) to modify the cosmological behaviour of the theory. We have investigated the effect of allowing for a non-zero $m_{\phi}$ and it appears that this term does not significantly improve the best-fit in the case of the restricted parameters of Model A. A reason for this is perhaps that even in the absence of MOND and Galileon terms in (2), the scalar sector of a canonical kinetic term together with a mass term is not equivalent to adding a massive scalar field because of the different gravitational effect of the scalar field in Einstein and matter frames; though the contribution of the scalar field to the Einstein frame Friedmann equation (37) may approximate that of a massive scalar field in some limits, this is not the case for the matter frame Friedmann equation (42) with ‘cross’ terms such as $H\tilde{Y}$ present. Due to the presence of the aether field $A_{\mu}$, one can build additional terms to add to the gravitational Lagrangian that project out time-like variations of $\phi$; for example, a term $A^{\mu}\tilde{\nabla}_{\mu}\phi A^{\nu}\tilde{\nabla}_{\nu}\phi|A^{\sigma}\tilde{\nabla}_{\sigma}\phi|$ will behave identically to the MOND kinetic term in the background cosmology but conceivably will have only a marginal effect on astrophysical scales. However, it is not clear how to consistently add or not add such terms to (2).

We have additionally checked that similar behavior existed in the original TeVeS theory. For Bekenstein’s original choice of free function, the theory possesses a limiting form describing a Canonical Regime and the theory’s parameters can be chosen such that a similar amount of rescaling of Newton’s constant as in best-fit g-TeVeS models can occur and it is found that parameters of the theory may be chosen that this regime exists for much of the matter era. Prior to the matter era, the theory transitions to a Nonlinear Regime which differs from the MOND Regime and does not give an additional dust-like contribution during the radiation era. This is perhaps the reason why the g-TeVeS theory gives a slightly closer fit to $\Lambda\mathrm{CDM}$ than TeVeS with Bekenstein’s function but the effect is small, as can be seen from the similarity of the results shown in Figures 6 and 7.

An important step in further investigating the cosmological consequences of g-TeVeS will be to develop the theory of cosmological perturbations as was done for the TeVeS theory Skordis et al. (2006); Skordis (2006). This will make it possible to deduce the effect of the new gravitational degrees of freedom on the cosmological microwave background and growth of large scale structure and will likely provide a powerful test of the theory.

Acknowledgements.

We thank Stéphane Ilić and Michael Kopp for many useful discussions. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) / ERC Grant Agreement n. 617656 “Theories and Models of the Dark Sector: Dark Matter, Dark Energy and Gravity.”

.10 Quantities in FRW Symmetry

Some useful quantities can be calculated:

[TABLE]

where $\gamma^{i}_{jk}$ are the Christoffel symbols corresponding to the metric $\gamma_{ij}$ and spatial derivatives acting on it. From this it follows that:

[TABLE]

For the scalar field ansatz $\phi=\phi(t)$, the only non-vanishing part of $\tilde{\nabla}_{\mu}\phi$ will be $\tilde{\nabla}_{t}\phi=d\phi/dt$. Given our previous results, it follows that:

[TABLE]

All other components of $\tilde{\nabla}_{\mu}\tilde{\nabla}_{\nu}\phi$ being zero. Additionally we have that:

[TABLE]

with all other components of $\tilde{\nabla}_{\mu}A_{\nu}$ vanishing.