<figure><img src="image_1.png" /></figure>

HAL
open science

# Revisiting metric perturbations in tensor-vector-scalar theory

Martin Feix

<figure><img src="image_2.png" /></figure>

## To cite this version:

Martin Feix. Revisiting metric perturbations in tensor-vector-scalar theory. Physical Review D, 2016, 93 (10), pp.104039. 10.1103/PhysRevD.93.104039. hal-01331647

HAL Id: hal-01331647

https://hal.sorbonne-universite.fr/hal-01331647v1

Submitted on 7 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Revisiting metric perturbations in tensor-vector-scalar theory

Martin Feix*

CNRS, UMR 7095 & UPMC, Institut d'Astrophysique de Paris, 98 bis Boulevard Arago, 75014 Paris, France
and Department of Physics, Technion—Israel Institute of Technology,

Technion City, 32000 Haifa, Israel

(Received 13 March 2012; published 20 May 2016)

I revisit cosmological perturbations in Bekenstein's tensor-vector-scalar theory (TeVeS). Considering only scalar modes in the conformal Newtonian gauge, the extra degrees of freedom are expressed in a way suitable for studying modifications at the level of the metric potentials. Assuming a universe in the matter-dominated phase, I discuss the mechanism responsible for boosting structure growth and confirm the vector field as its key ingredient. Using a semianalytic approach, I further characterize the evolution of density perturbations and the potentials on sub- and superhorizon scales.

DOI: 10.1103/PhysRevD.93.104039

# I. INTRODUCTION

Despite the remarkable success of the standard cosmological model [1], open questions related to the nature of cosmic acceleration and dark matter [2] still leave room to ponder alternative possibilities. Among the various proposals, it has been suggested that the gravitational dynamics may be different from the predictions of general relativity (GR). Additionally, motivated by considerations beyond the field of cosmology, this has triggered theoretical developments aimed at finding viable models, and there exist now a large number of modified theories which typically introduce extra degrees of freedom in the gravitational sector [3].

The phenomenology of such modified gravity theories is quite rich and can involve complex and nonintuitive results which are hard to understand from the equations of motion alone. A way forward is to use appropriate parametrizations which capture the new fields and their dynamics, but allow one to view their effects from a different angle and to explore how familiar quantities and observables are affected. Here I want to follow this approach using the example of Bekenstein's tensor-vector-scalar theory (TeVeS) [4] which was originally constructed as a relativistic extension for the modified Newtonian dynamics (MOND) paradigm [5] and has been subject to numerous studies in the literature [6]. Considering perturbations around a spatially flat Friedmann-Robertson-Walker (FRW) background, I will express the additional degrees of freedom in a way suitable for studying modifications at the level of the metric potentials in the conformal Newtonian gauge. This particular ansatz will then be used to study several aspects of the cosmological evolution such as the modified growth of density perturbations and the gravitational slip which has been identified as a generic feature of modified gravity theories [7]. To allow an

analytic treatment of the cosmological background, I will adopt an Einstein-de Sitter (EdS) universe which only contains pressureless matter and provides an excellent approximation to the matter era of a realistic universe. The approach taken in this work differs in motivation from the recent attempts of parametrizing deviations from GR in a model-independent fashion [8] which are interesting in their own right.

The paper is organized as follows: Starting with a brief review of TeVeS in Sec. II, its background cosmology and the modified EdS model are presented in Sec. III. Considering perturbations in the conformal Newtonian gauge, I then introduce a parametrization of the new variables and discuss the evolution of cosmological perturbations in Sec. IV. Finally, I conclude in Sec. V. For clarity, some of the material involving lengthy expressions is presented in an appendix. If not stated otherwise, I will assume the notation of Ref. [9] throughout this work.

# II. FUNDAMENTALS OF TEVES

## A. Fields and action

In its original form, TeVeS [4,10] is a bimetric gravity theory which is based on three dynamical fields: an Einstein metric $\tilde{g}_{\mu\nu}$, a timelike vector field $A_\mu$ such that

$$ \tilde{g}^{\mu\nu} A_{\mu} A_{\nu} = -1, \qquad (1) $$

and a scalar field $\phi$. The gravity-matter coupling involves a second metric $g_{\mu\nu}$ which is obtained from

$$ g_{\mu\nu} = e^{-2\phi} \tilde{g}_{\mu\nu} - 2A_{\mu}A_{\nu} \sinh(2\phi). \qquad (2) $$

The frames delineated by the metric fields $\tilde{g}_{\mu\nu}$ and $g_{\mu\nu}$ are called the Einstein frame and matter frame, respectively. The geometric part of the action is the same as in GR,

<footnote>*feix@iap.fr</footnote>
$$S_g = \frac{1}{16\pi G} \int \tilde{g}^{\mu\nu} \tilde{R}_{\mu\nu} \sqrt{-\tilde{g}} d^4 x, \quad (3)$$

where $\tilde{R}_{\mu\nu}$ is the Ricci tensor of $\tilde{g}_{\mu\nu}$ and $\tilde{g}$ the determinant of $\tilde{g}_{\mu\nu}$. The vector field's action $S_v$ reads as follows,

$$S_v = -\frac{1}{32\pi G} \int [K_B F^{\mu\nu} F_{\mu\nu} - \lambda(A_\mu A^\mu + 1)] \sqrt{-\tilde{g}} d^4 x, \quad (4)$$

where $F_{\mu\nu} = \tilde{\nabla}_{\mu}A_{\nu} - \tilde{\nabla}_{\nu}A_{\mu}$ and indices are raised and lowered with respect to $\tilde{g}_{\mu\nu}$, i.e. $A^{\mu} = \tilde{g}^{\mu\nu}A_{\nu}$. Here the constant $K_B$ describes the coupling of the vector field to gravity, and $\lambda$ is a Lagrangian multiplier enforcing the normalization condition of $A_\mu$. The action $S_s$ of the scalar field $\phi$ involves an additional nondynamical scalar field $\mu$ and takes the form

$$S_s = -\frac{1}{16\pi G} \int [\mu h^{\mu\nu} \tilde{\nabla}_{\mu} \phi \tilde{\nabla}_{\nu} \phi + V(\mu)] \sqrt{-\tilde{g}} d^4 x, \quad (5)$$

where $h^{\mu\nu} = \tilde{g}^{\mu\nu} - A^\mu A^\nu$ and $V(\mu)$ is an initially arbitrary (potential) function. As the field $\mu$ is related to the invariant $h^{\mu\nu} \tilde{\nabla}_{\mu} \phi \tilde{\nabla}_{\nu} \phi$, however, it could in principle be eliminated from the action. Finally, the matter action is given by

$$S_m = \int \mathcal{L}_m[g, \Upsilon^B, \nabla\Upsilon^B] \sqrt{-\tilde{g}} d^4 x, \quad (6)$$

where $\Upsilon^B$ is a generic collection of matter fields. By construction, world lines are geodesics of $g_{\mu\nu}$ rather than $\tilde{g}_{\mu\nu}$. Despite its explicit bimetric construction, TeVeS may be written in pure tensor-vector form [11] and provide a particular example of general Aether-type theories [12].

## B. Choice of the potential, quasistatic systems and relation to MOND

The behavior of TeVeS in the nonrelativistic limit depends on the assumed potential $V$. Originally, Bekenstein made the choice

$$V(\mu) = \frac{3\mu_0^2}{128\pi l_B^2} [\hat{\mu}(4 + 2\hat{\mu} - 4\hat{\mu}^2 + \hat{\mu}^3) + 2\log(1 - \hat{\mu})^2], \quad (7)$$

where the constant $l_B$ corresponds to a length scale, $\hat{\mu} = \mu/\mu_0$, and $\mu_0$ is a dimensionless constant [13]. GR is then recovered in the limit $K_B \to 0$ and $l_B \to \infty$. Applying the usual approximations for weak fields and quasistatic systems, one finds $V' \equiv dV/d\mu < 0$, and therefore $0 < \mu < \mu_0$. Using that also $V(\mu) < 0$ for the given range, the metric $g_{\mu\nu}$ turns out to be identical to the metric obtained in GR if the nonrelativistic gravitational potential is replaced by

$$\begin{aligned} W &= \Xi \Phi_N + \phi, \\ \Xi &= e^{-2\phi_C} (1 + K_B/2)^{-1}, \end{aligned} \quad (8)$$

where $\phi_C$ is the cosmological value of $\phi$ at the time the system in question breaks away from the cosmological expansion, and $\Phi_N$ is the Newtonian potential generated by the matter density $\rho$ [14]. In this approximation, one further has

$$h^{\mu\nu} \tilde{\nabla}_{\mu} \phi \tilde{\nabla}_{\nu} \phi \to (\nabla \phi)^2 \equiv \|\nabla \phi\|_2^2, \quad (9)$$

and the equation of the scalar field reduces to

$$\nabla \cdot (\mu \nabla \phi) = 8\pi G\rho. \quad (10)$$

Equations (8) and (10) correspond to the MOND paradigm [4]. If $\mu \to \mu_0$, the theory reaches its Newtonian limit, and the measured gravitational constant $G_N$ is given by

$$G_N = \frac{\mu_0 + 2 - K_B}{\mu_0(1 - K_B/2)} G. \quad (11)$$

Similarly, the theory reaches its MONDian limit as $\mu \to 0$, and the acceleration constant $a_0$ can be expressed in terms of the TeVeS and potential parameters,

$$a_0 = \frac{\sqrt{6} e^{\phi_C}}{2l_B \sqrt{\pi \mu_0}} \frac{G}{G_N}. \quad (12)$$

As can be seen from above, $a_0$ depends on $\phi_C$ and may therefore, in principle, change with time [15].

Different choices of $V$ and their implications for bound structures have been studied in Refs. [16–18]. Focusing on cosmological scales, I will assume the following class of potentials in this paper [18],

$$V_n(\mu) = \frac{3\mu_0^2}{32\pi l_B^2} \left[ \frac{n+4+(n+1)\hat{\mu}}{(n+1)(n+2)} (\hat{\mu}-2)^{n+1} + \frac{(-1)^n}{2} \log(1-\hat{\mu})^2 + \sum_{m=1}^n \frac{(-1)^{n-m}}{m} (\hat{\mu}-2)^m \right], \quad (13)$$

where $n \ge 2$ [19]. Adopting different values of $n$, Fig. 1 illustrates the resulting potential shape as a function of $\hat{\mu}$. Such a generalized potential reduces to Bekenstein's toy model if $n = 2$. The derivative of $V_n(\mu)$ takes a simpler form and can be expressed as

$$V'_n(\mu) = \frac{3\mu_0}{32\pi l_B^2} \frac{\hat{\mu}^2 (\hat{\mu}-2)^n}{\hat{\mu}-1}. \quad (14)$$

Requiring that $V'$ is single-valued and $V' \ge 0$ [4], one is always free to choose between two possible potential
<figure><img src="image_4.png" /><figcaption>FIG. 1. Illustration of the generalized potential function $V_n(\mu)$ given by Eq. (13) for $n = 2$ (solid line), 3 (dotted line), and 6 (dashed line).</figcaption></figure>

branches. Here I shall use the branch ranging from the extremum at $\mu = 2\mu_0$ to infinity. Under these preliminaries, the potential in Eq. (13) gives rise to tracker solutions of the scalar field [18], with a background evolution similar to other theories involving tracker fields [20]. I will further elaborate on this behavior in Sec. III B.

### III. COSMOLOGICAL BACKGROUND

#### A. Evolution equations

Imposing the usual assumptions of an isotropic and homogeneous spacetime, both $g_{\mu\nu}$ and $\tilde{g}_{\mu\nu}$ are given by FRW metrics with scale factors $a$ and $b = ae^{\bar{\phi}}$, respectively, where $\bar{\phi}$ is the background value of the scalar field [4,9]. Adopting a spatially flat universe, the modified Friedmann equation in the matter frame reads

$$3H^2 = 8\pi G_{\text{eff}}(\bar{\rho}_{\phi} + \bar{\rho}), \qquad (15)$$

where the physical Hubble parameter is $H = \dot{a}/a^2$ and the overdot denotes the derivative with respect to conformal time. Here $\bar{\rho}$ corresponds to the FRW background density of the fluid and the scalar field density takes the form

$$\bar{\rho}_{\phi} = \frac{e^{2\bar{\phi}}}{16\pi G}(\bar{\mu}V' + V). \qquad (16)$$

The effective gravitational coupling strength is given by

$$G_{\text{eff}} = Ge^{-4\bar{\phi}} \left( 1 + \frac{d\bar{\phi}}{d \log a} \right)^{-2}, \qquad (17)$$

which is generally time-varying through its dependence on the scalar field $\bar{\phi}$. Just as in GR, the energy density $\bar{\rho}$ evolves according to

$$\dot{\bar{\rho}} = -3\frac{\dot{a}}{a}(1+w)\bar{\rho}, \qquad (18)$$

where $w$ is the equation-of-state (EoS) parameter of the fluid. In the case of multiple background fluids, i.e. $\bar{\rho} = \sum_i \bar{\rho}_i$, the relative densities $\Omega_i$ are defined as

$$\Omega_i = 8\pi G_{\text{eff}} \frac{\bar{\rho}_i}{3H^2} = \frac{\bar{\rho}_i}{\bar{\rho} + \bar{\rho}_{\phi}}. \qquad (19)$$

The evolution of the scalar field $\phi$ is governed by

$$\ddot{\bar{\phi}} = \dot{\bar{\phi}}\left(\frac{\dot{a}}{a} - \dot{\bar{\phi}}\right) - \frac{1}{U}\left[3\bar{\mu}\frac{\dot{b}}{b}\dot{\bar{\phi}} + 4\pi G a^2 e^{-4\bar{\phi}}(\bar{\rho} + 3\bar{P})\right], \qquad (20)$$

where $\bar{P}$ is the fluid's background pressure, and the function $U$ is related to the potential $V$,

$$U(\bar{\mu}) = \bar{\mu} + 2 \frac{V'}{V''}. \qquad (21)$$

In addition, the scalar field obeys the constraint equation,

$$\dot{\bar{\phi}}^2 = \frac{1}{2} a^2 e^{-2\bar{\phi}} V', \qquad (22)$$

which can be inverted to obtain $\bar{\mu}(a, \bar{\phi}, \dot{\bar{\phi}})$. For later use, I also introduce the relation

$$2\frac{\dot{a}\dot{b}}{a\bar{b}} - \frac{\ddot{b}}{b} - \bar{\mu}\dot{\bar{\phi}}^2 = 4\pi G a^2 e^{-4\bar{\phi}}(\bar{\rho} + \bar{P}), \qquad (23)$$

which follows from combining Eq. (15) with Eq. (22) and the corresponding Raychaudhuri equation [9]. For a broad class of potentials $V$, one typically finds $\exp \bar{\phi} \approx 1$ and $\bar{\rho}_{\phi} \ll 1$ throughout cosmological history [4,18,21]. Therefore, the background evolution is very similar to the standard case of GR, with only small corrections induced by the scalar field.

#### B. Tracker solutions of the scalar field

For the potentials specified by Eq. (13), it has been found that the scalar field exhibits a (stable) tracking behavior and synchronizes its energy density with the dominant component of the universe [18,21]. Tracking occurs as $V'$ tends to its zero point where $\bar{\mu} = 2\mu_0$, and the evolution of the field $\bar{\phi}$ is then approximately given by
$$ \bar{\phi} = \bar{\phi}_0 + \frac{|1 + 3w|}{2\beta\mu_0|1 - w| - |1 + 3w|} \log a, \quad (24) $$

where $\bar{\phi}_0$ is an integration constant and $\beta = \pm 1$. The sign depends on the matter fluid's EoS parameter $w$ and Eq. (20). The density $\bar{\rho}_\phi$ then exactly scales like that of the fluid, and the relative density parameter $\Omega_\phi$ turns approximately into a constant,

$$ \Omega_\phi = \frac{(1 + 3w)^2}{6\mu_0(1 - w)^2}. \quad (25) $$

The right-hand side of Eq. (24) slightly differs from the expression presented in Ref. [18]. In Appendix A, I show that Eq. (24) is indeed the correct result.

During tracking, $\bar{\mu}$ can be expressed as $\bar{\mu} = 2\mu_0(1 + \epsilon)$ with $0 < \epsilon \ll 1$. Using $V'(2\mu_0) = 0$ and expanding $V'$ to lowest order in $\epsilon$, Eq. (22) leads to

$$ \epsilon = \frac{1}{2} \left( \frac{16\pi l_B^2}{3\mu_0} \frac{e^{2\bar{\phi}}}{a^2} \dot{\bar{\phi}}^2 \right)^{1/n} . \quad (26) $$

It turns out that this is the only stage at which the constant $l_B$ enters the evolution equations. Taking the time derivative of the above yields the useful relation

$$ \dot{\bar{\phi}}\dot{\epsilon} = \frac{2}{n} \left( \dot{\bar{\phi}}^2 - \dot{\bar{\phi}}\frac{\dot{a}}{a} + \ddot{\bar{\phi}} \right) \epsilon. \quad (27) $$

Stable tracking requires $\epsilon$ to asymptotically decrease to zero, i.e. $\epsilon \to 0$. Therefore, one has the condition $\dot{\epsilon} < 0$ which may be used to infer the proper sign of the parameter $\beta$ in Eq. (24).

### C. Modified Einstein-de Sitter cosmology

In what follows, I shall assume a universe entirely made of pressureless matter with perfect tracking of the scalar field, corresponding to the EdS model in GR. Setting $\bar{P} = w = 0$ fixes $\beta = -1$, and thus the scalar field can be written as

$$ \bar{\phi} = \bar{\phi}_0 - \frac{1}{2\mu_0 + 1} \log a. \quad (28) $$

To find the proper value of $\beta$, one may either insert Eq. (24) into Eq. (20) or use the argument presented in Appendix A. Since the fluid evolves according to Eq. (18), the density takes the form $\bar{\rho} = \bar{\rho}_0 a^{-3}$, where $\bar{\rho}_0$ is the background density's value today. Exploiting Eq. (25) allows one to rewrite the modified Friedmann equation in the matter frame as

$$ H^2 = H_0^2 a^{-3+4/(2\mu_0+1)}, \quad (29) $$

where I have used

$$ H_0^2 = e^{-4\bar{\phi}_0} \frac{8\pi G\bar{\rho}_0}{3} \left(1 + \frac{1}{6\mu_0 - 1}\right) \left(1 - \frac{1}{2\mu_0 + 1}\right)^{-2}. \quad (30) $$

The deviation of the Hubble expansion from the ordinary EdS case is entirely characterized by the parameter $\mu_0$. For several reasons, $\mu_0$ should take a rather large value on the order of 100–1000 [4]. Thus, this deviation will be small, typically at the percent level in the range of practical interest.

## IV. METRIC PERTURBATIONS IN TEVES

### A. Preliminaries

#### 1. Matter-frame perturbations

Now I will turn to metric perturbations around a spatially flat FRW spacetime in TeVeS. The starting point is the set of linear perturbation equations derived in Ref. [9]. For simplicity, I shall consider only scalar modes and work within the conformal Newtonian gauge. Perturbations of the metric are then characterized by two scalar potentials $\Psi$ and $\Phi$, and the line element in the matter frame is given by

$$ ds^2 = a^2[-(1 + 2\Psi)d\tau^2 + (1 - 2\Phi)\delta_{ij}dx^i dx^j]. \quad (31) $$

Similarly, one needs to consider perturbations of the other fields: While the fluid perturbation variables are defined in the usual way; i.e., the density perturbation, for instance, is expressed in terms of the density contrast $\delta$,

$$ \rho = \bar{\rho} + \delta\rho = \bar{\rho}(1 + \delta), \quad (32) $$

and the scalar field is perturbed as

$$ \phi = \bar{\phi} + \varphi, \quad (33) $$

where $\varphi$ is the scalar field perturbation. Finally, the perturbed vector field is written as

$$ A_{\mu} = ae^{-\bar{\phi}}(\bar{A}_{\mu} + \alpha_{\mu}), \quad (34) $$

where $\bar{A}_\mu = (1, 0, 0, 0)$ and

$$ \alpha_{\mu} = (\Psi - \varphi, \nabla\alpha). \quad (35) $$

The time component of the vector field perturbation is constrained to be a combination of metric and scalar field perturbations, which is a consequence of the unit-norm condition in Eq. (1). Therefore, one needs to consider only the longitudinal perturbation component $\alpha$.
## 2. Einstein-frame perturbations

Instead of using Eq. (31), one may also express perturbations in the Einstein frame [9,21]. In this case, metric perturbations are written as

$$ \tilde{g}_{00} = -b^2 e^{-4\bar{\phi}}(1 + 2\tilde{\Psi}), \quad (36) $$

$$ \tilde{g}_{0i} = -b^2 \partial_i \tilde{\zeta}, \quad (37) $$

$$ \tilde{g}_{ij} = b^2(1 - 2\tilde{\Phi})\delta_{ij}, \quad (38) $$

where the Einstein-frame perturbations are given by

$$ \tilde{\Psi} = \Psi - \varphi, \quad (39) $$

$$ \tilde{\Phi} = \Phi - \varphi, \quad (40) $$

$$ \tilde{\zeta} = (e^{-4\bar{\phi}} - 1)\alpha. \quad (41) $$

To avoid lengthy expressions in the perturbed field equations, it is convenient to work with variables from both frames. Since the equation governing the evolution of the scalar field perturbation is of second order, it is further helpful to introduce an auxiliary field $\gamma$ which allows one to split the scalar field equation into a system of two first-order equations [9]. For the present gauge choice, the field $\gamma$ is given by

$$ \gamma = 2e^{\bar{\phi}} \frac{U}{a} (-\dot{\varphi} + \dot{\bar{\phi}} \tilde{\Psi}). \quad (42) $$

## B. Parametrizing the new degrees of freedom

To explore how the new degrees of freedom modify cosmological dynamics, it is helpful to define closure relations which may be used to eliminate perturbations of the scalar and vector field from the evolution equations. Consider a new function $B_\varphi$ such that

$$ \log B_\varphi = \tilde{\Psi} - \varphi. \quad (43) $$

The function $B_\varphi$ is generally time- and scale-dependent and describes the magnitude of scalar field perturbations relative to the metric potential $\tilde{\Psi}$. Its definition is not completely arbitrary, but motivated by the algebraic structure of the perturbation equations. The idea is to try to express the new variables in terms of metric and matter fluid perturbations, which typically yields equations that are trivially satisfied because of vanishing background terms. The form of these equations can then serve as a hint for defining suitable relations. Since such a procedure is neither fundamental nor unique, however, one should just take Eq. (43) as an educated guess.

Similarly, I define a relation for the variables $\gamma$ and $\tilde{\zeta}$,

$$ \log B_{\gamma} \propto 3 \dot{\tilde{\Phi}} + k^2 \tilde{\zeta} - \frac{ae^{-\bar{\phi}}}{2U} \gamma, \quad (44) $$

where $\mathbf{k}$ is the conformal wave vector, and $k = |\mathbf{k}|$. In what follows, the proportionality constant is set to unity for simplicity. $B_\gamma$ characterizes a combination of both the scalar and the vector field perturbations. Because of the term $k^2\tilde{\zeta}$, the relative contribution due to the vector field perturbation $\alpha$ is scale-dependent and becomes negligible as $k \to 0$. Likewise, it will dominate the expression on very small scales within the horizon. Hence $B_\gamma$ may be solely associated with vector field perturbations if one considers the theory on subhorizon scales. The functions $B_\varphi$ and $B_\gamma$ are strictly positive by construction and, thus, well defined throughout cosmological history.

Finally, let me also introduce

$$ \begin{aligned} \tilde{B}_{\varphi} &= \log B_{\varphi}, & \tilde{B}_{\varphi} \frac{\dot{a}}{a} &= \frac{\dot{B}_{\varphi}}{B_{\varphi}}, \\ \tilde{B}_{\gamma} \frac{\dot{a}}{a} &= \log B_{\gamma}, & \tilde{B}_{\gamma} \frac{\dot{a}^2}{a^2} &= \frac{\dot{B}_{\gamma}}{B_{\gamma}}, \end{aligned} \quad (45) $$

which allows one to express the modifications in a more convenient way. In terms of these new functions, the dynamics of GR is recovered in the limit

$$ \tilde{B}_{\varphi} \rightarrow \tilde{\Psi}, \quad \tilde{B}_{\gamma} \frac{\dot{a}}{a} \rightarrow 3\dot{\tilde{\Phi}}, \quad \bar{\phi} \rightarrow 0. \quad (46) $$

## C. Subhorizon scales

In the following, I shall assume the previously discussed modified EdS cosmology with perfect tracking of the scalar field. This allows one to use the corresponding background expressions presented in Sec. III C and considerably simplifies the analysis of the modified equations. Since $\bar{\mu}$ resides close to its minimum in this case, i.e. $\bar{\mu} = 2\mu_0(1 + \epsilon)$ with $\epsilon \ll 1$, one may further exploit the two first-order expressions,

$$ \frac{\bar{\mu}}{U} = 1 - \frac{2}{n}\epsilon \quad (47) $$

and

$$ \frac{2\mu_0}{U} = 1 - \frac{n+2}{n}\epsilon, \quad (48) $$

which are useful to rewrite terms involving the field $U$.

## 1. Modified potentials

Adopting an EdS universe together with the closure relations presented in Sec. IV B, one may now write metric
perturbations solely in terms of auxiliary functions and matter variables. The resulting equations are quite lengthy and can be found in Appendix B. As a first application, I consider scales much smaller than the horizon, i.e. $aH/k \ll 1$. Inserting the logarithmic approximation for $\bar{\phi}$ specified by Eq. (28), I expand the Einstein-frame potentials $\tilde{\Psi}$ and $\tilde{\Phi}$ in powers of $aH/k$. To second order, this yields

$$ \tilde{\Psi} = A_0 + A_2 \frac{a^2 H^2}{k^2} + \mathcal{O}(\epsilon), \quad (49) $$

$$ \tilde{\Phi} = B_0 + B_2 \frac{a^2 H^2}{k^2} + \mathcal{O}(\epsilon), \quad (50) $$

where

$$ A_0 = \frac{[K_B - 2(1 - e^{4\bar{\phi}})]\tilde{B}_{\bar{\phi}}}{2(1 - e^{4\bar{\phi}}) - K_B(1 - 2e^{4\bar{\phi}})}, \quad (51) $$

$$ B_0 = \frac{e^{4\bar{\phi}} K_B \tilde{B}_{\bar{\phi}}}{2(1 - e^{4\bar{\phi}}) - K_B(1 - 2e^{4\bar{\phi}})}, \quad (52) $$

and $A_2, B_2$ are complicated expressions involving $\delta, \bar{\phi}$, and the functions in Eq. (45). Compared to GR, where the EdS model gives $A_0 = B_0 = 0$ and $A_2 = B_2 = -3\delta/2$, the potentials exhibit a sophisticated dependence on $\bar{\phi}$ which is expected to have a significant impact on their evolution.

Previous studies of TeVeS cosmology found that structure can form more efficiently than in GR and identified a vector instability, which occurs for small enough values of $K_B$, as the key ingredient for enhanced growth [21–23]. Scalar field perturbations, on the other hand, were argued to play only a negligible role for structure formation. Here I seek to study this mechanism at the level of the modified potentials. To this end, I will also consider full numerical solutions which were obtained with the modified Boltzmann code of Ref. [21]. These solutions serve as a consistency check of the used assumptions, and allow one to determine the auxiliary functions $B_\phi$ and $B_\gamma$ and to characterize their effect on the potentials. The solver adopts standard adiabatic initial conditions where perturbations of the TeVeS scalar and vector fields are set to zero. Assuming three neutrinos with a mass of 2.3 eV and no dark matter, all calculations in this work use a spatially flat cosmology with density parameters $\Omega_m \approx 0.04$ and $\Omega_\nu \approx 0.15$, dimensionless Hubble constant $h = 0.74$, and fixed TeVeS parameters $l_B = 100$ Mpc and $n = 2$.

Let me now focus on $\tilde{\Psi}$ which dictates the dynamics of nonrelativistic matter. Since $|\bar{\phi}| \ll 1$ and $\mu_0 \gg 1$, Eq. (49) can be further simplified by taking the limit $\mu_0 \to \infty$ [24] and keeping only leading-order terms with respect to $\bar{\phi}$. In the process, one has to be careful when dealing with

<figure><img src="image_5.png" /><figcaption>FIG. 2. Evolution of $\tilde{\Psi}$ for $K_B = 1$ (top) and $K_B = 0.07$ (bottom). Assuming $\mu_0 = 250$ and a fixed wave number $k = 0.5$ Mpc$^{-1}$, the figure shows the full numerical (solid lines) and semianalytic (dashed lines) results described in the text. The dotted line indicates the contribution of the zeroth-order term $A_0$.</figcaption></figure>

expressions involving the constant $K_B$, especially if $K_B$ is comparable to $\bar{\phi}$:

(a) $K_B \gtrsim 1$. In this case, it turns out that

$$ A_0 \approx \tilde{B}_{\bar{\phi}} \quad (53) $$

and

$$ A_2 \approx -2(\tilde{B}_\gamma + \hat{B}_\gamma - 10\tilde{B}_\phi - 4\hat{B}_\phi). \quad (54) $$

If the sum of $\tilde{B}_\gamma$ and the two time-derivative terms is not too large, the zeroth-order contribution, $A_0$, will dominate, and $\tilde{\Psi} \approx \tilde{B}_\phi$ during matter domination. The top panel of Fig. 2 shows a comparison between the numerically computed potential and $\tilde{B}_\phi$, assuming $K_B = 1$, $\mu_0 = 250$, and a wave number $k = 0.5$ Mpc$^{-1}$. Indeed, the function $\tilde{B}_\phi$ matches the true potential quite well, which suggests a negligible contribution from vector field perturbations as long as $K_B$ is large enough. Since the scalar field perturbation typically satisfies $\varphi \ll \tilde{\Psi}$, the approximation $\tilde{\Psi} \approx \tilde{B}_\phi$ also follows independently from the definition of $B_\phi$ and is valid for any choice of $K_B$.

(b) $K_B \ll 1$. To leading order, one finds

$$ A_0 \approx \frac{1 - 2\tilde{K}^{-1}}{1 + 2\tilde{K}^{-1}} \tilde{B}_{\bar{\phi}} \quad (55) $$
and

$$A_2 \approx 4 \frac{5 - 41 \tilde{K}^{-1} + 56 \tilde{K}^{-2} - 32 \tilde{K}^{-3}}{(1 + 2 \tilde{K}^{-1})^2} \tilde{B}_{\varphi} + \tilde{K}^{-1} \frac{4(\tilde{B}_{\gamma} + \tilde{B}_{\gamma} - 5\tilde{B}_{\varphi}) - 6\delta}{1 + 2\tilde{K}^{-1}} - 2 \frac{\tilde{B}_{\gamma} + \tilde{B}_{\gamma} - 4\tilde{B}_{\varphi}}{1 + 2\tilde{K}^{-1}}, \quad (56)$$

where $\tilde{K}^{-1}K_B = 1 - e^{4\tilde{\phi}}$. There are two main differences with respect to the previous case. First, $\delta$ explicitly appears in the equation, sourcing the potential directly and not only through the functions $B_\varphi$ and $B_\gamma$. Second, the terms $A_0$ and $A_2$ now strongly depend on the ratio $\tilde{K}^{-1}$. Calculations for the background yield $\tilde{\phi} < 0$ and $\tilde{K}^{-1} > 0$, resulting in a suppression of $A_0$ relative to higher-order contributions. Hence, the second-order term, specified by $A_2$, can be expected to be important. The bottom panel of Fig. 2 shows $\tilde{\Psi}$ together with the semianalytic approximation using Eqs. (55) and (56) for $K_B = 0.07$, $\mu_0 = 250$, and $k = 0.5$ Mpc$^{-1}$. The full numerical result is again matched very well, but this time the contribution of $A_0$ is much smaller and negative, growing up to a fraction of roughly 30% near the end of matter domination. Considering the term $A_2$, only expressions related to $B_\gamma$ and $\delta$ turn out as relevant sources.

The above demonstrates that the algebraic mechanism triggered by small values of the vector coupling, $K_B \sim |4\tilde{\phi}|$, is an important element in generating enhanced growth. Over the range $0.07 < K_B < 1$, the full numerical analysis shows that the ratio $\tilde{B}_\gamma/\tilde{B}_\varphi$ changes only by a factor of less than 3–4 (and similarly for $\tilde{B}_\gamma/\tilde{B}_\varphi$), not enough to explain why the terms involving $B_\varphi$ and its time derivative are suppressed. Neglecting the new contributions and changes in the above equations does not lead to an augmentation of growth. In accordance with Ref. [22], this points toward the implementation of the TeVeS vector field as the key ingredient for boosting structure formation. Once enhanced growth occurs, the potential is found to be sourced by $B_\varphi$, $B_\gamma$, and $\delta$ at comparable levels. As the evolution equations of the potentials, the scalar and vector fields, and the density contrast are coupled to each other, however, this complex interplay is hardly surprising.

How do the potentials $\Psi$ and $\tilde{\Psi}$ evolve with cosmic time? To find an answer, I start from the observation that $\tilde{B}_\varphi \sim \tilde{B}_\varphi$ in the late matter era (and similarly for $\tilde{B}_\gamma$ and $\tilde{B}_\varphi$). A closer inspection suggests the relation

$$\tilde{B}'_{\varphi} \propto a^{p-1} \tilde{B}_{\varphi}, \quad (57)$$

<figure><img src="image_6.png" /><figcaption>FIG. 3. Evolution of $\Psi$ for $k = 0.5$ Mpc$^{-1}$ and different combinations of $K_B$ and $\mu_0$. The figure illustrates the numerical result (solid lines) together with the analytic approximation (dashed lines) valid in the late matter era.</figcaption></figure>

where the prime denotes the derivative with respect to the scale factor $a$, and the constant $p$ depends on the parameters of the theory. Remarkably, one finds

$$p^{-1} = \mu_0 K_B, \quad (58)$$

which, using the tracking solution Eq. (28) and $|\tilde{\phi}| \ll 1$, can be related to the background scalar field through

$$a^p \propto e^{-2\tilde{\phi}/K_B}. \quad (59)$$

Integrating by separation of variables then leads to

$$\log \tilde{B}_{\varphi} = c_1 a^p + c_2, \quad (60)$$

where $c_1$ is fixed by the proportionality constant and $c_2$ is a constant of integration. Since $\varphi \ll \tilde{\Psi}$, one has $\Psi \approx \tilde{\Psi} \approx \tilde{B}_\varphi$, and the potentials approximately satisfy Eq. (60). Figure 3 illustrates the evolution of $\Psi$ together with the analytic approximation, assuming $k = 0.5$ Mpc$^{-1}$ and different combinations of $K_B$ and $\mu_0$. The simple models agree well with the numerical results for $a \sim 0.01$ and later. In the case of $\mu_0 K_B = 250$, the additional growth is strongly suppressed and $\Psi$ basically approaches a constant. The result is found to be independent of scale as long as $aH/k \ll 1$ is satisfied. Given the complexity of TeVeS, the appearance of such an effective description is astonishing. Since the growth of perturbations is mainly determined by the ratio $\tilde{\phi}/K_B$, one may expect this behavior also for other
parameter choices. Exploring its origin is beyond the scope of this paper, but might be addressed in future work.

## 2. Gravitational slip

Next, I want to address how the TeVeS growth mechanism manifests itself in a difference between the matter-frame potentials $\Psi$ and $\Phi$. Introducing the gravitational slip as

$$ \xi_G = \Psi - \Phi = \tilde{\Psi} - \tilde{\Phi}, \qquad (61) $$

I will, just as before, consider the limit $\mu_0 \to \infty$ and keep only leading-order terms with respect to $\bar{\phi}$. All comparisons to full numerical results assume $\mu_0 = 250$ and $k = 0.5$ Mpc$^{-1}$:

(a) $K_B \gtrsim 1$. Collecting the relevant terms yields

$$ A_0 - B_0 \approx (e^{4\bar{\phi}} - 1)\tilde{B}_{\phi} \qquad (62) $$

and

$$ A_2 - B_2 \approx \tilde{B}_{\gamma} + 4\tilde{B}_{\phi} + \hat{B}_{\phi}. \qquad (63) $$

Unlike the case of $\tilde{\Psi}$, a factor of $e^{4\bar{\phi}} - 1$ suppresses the term $A_0 - B_0$, and hence it should become comparable to higher-order contributions. The top panel of Fig. 4 compares $\xi_G$ to its semianalytic approximation using Eqs. (62) and (63) for $K_B = 1$. Qualitatively, the approximation is in agreement with the full numerical result. As expected, the term $A_0 - B_0$ is not dominant

<figure><img src="image_7.png" /><figcaption>FIG. 4. Same as Fig. 2, but now for the slip $\xi_G$. The dotted lines indicate the contribution of the zeroth-order term $A_0 - B_0$.</figcaption></figure>

and increases to a maximum relative contribution of around 25% at the end of the matter era. The numerical results show that most of $\xi_G$ is due to the second-order expression specified by $A_2 - B_2 \approx \tilde{B}_{\gamma}$ which, according to Eq. (44), can be associated with the vector field perturbation $\tilde{\xi}$.

(b) $K_B \ll 1$. A straightforward calculation leads to

$$ A_0 - B_0 \approx -\frac{2\tilde{K}^{-1}}{1 + 2\tilde{K}^{-1}} \tilde{B}_{\phi} \qquad (64) $$

and

$$ A_2 - B_2 = 2 \frac{2 - 29\tilde{K}^{-1} + 80\tilde{K}^{-2} - 56\tilde{K}^{-3}}{(1 + 2\tilde{K}^{-1})^2} \tilde{B}_{\phi} + \tilde{K}^{-1} \frac{4(\hat{B}_{\gamma} - 4\hat{B}_{\phi}) + 3(2\tilde{B}_{\gamma} - \delta)}{1 + 2\tilde{K}^{-1}} + \frac{\tilde{B}_{\gamma} + \hat{B}_{\phi}}{1 + 2\tilde{K}^{-1}}. \qquad (65) $$

For sufficiently small values of $K_B$, the effect on $\xi_G$ is very similar to the case of $\tilde{\Psi}$, but now both the first and second-order terms increase to nearly the same degree. The bottom panel of Fig. 4 shows $\xi_G$ together with the semianalytic approximation using Eqs. (64) and (65) for $K_B = 0.07$. Again, the full numerical result is well reproduced, and the second-order term, determined by $A_2 - B_2$, is dominated by the terms involving $B_{\gamma}$, its time derivative, and $\delta$. The relative contribution of $A_0 - B_0$ is almost identical to the previous case $K_B = 1$. Apart from serving as a consistency check and validation of the used assumptions, the above results confirm that the mechanisms responsible for boosting the potential $\Psi$ and for generating a considerable slip $\xi_G$ are, indeed, the very same.

## 3. Growth of density perturbations

To conclude the section on subhorizon scales, consider the growth of density perturbations during the matter-dominated phase. Since only the gravitational part is modified, one may follow the usual derivation of the growth equation and obtains

$$ a^2 \delta'' + \frac{a}{2} \left( 3 + \frac{4}{2\mu_0 + 1} \right) \delta' + \left( \frac{k}{aH} \right)^2 \Psi = 0. \qquad (66) $$

If $p = (\mu_0 K_B)^{-1}$ is chosen small enough, the analytic expression of $\Psi$ given by Eq. (60) can be approximated by its first-order expansion. Using $aH \propto a^{-1/2}$ then yields a relation of the form
<figure><img src="image_8.png" /><figcaption>FIG. 5. Growth of density perturbations for different values of $K_B$, $\mu_0 = 250$, and $k = 0.5$ Mpc$^{-1}$. The top panel shows the evolution of $\delta$ (solid lines) together with analytic models based on approximate solutions to the modified growth equation (dashed lines; see text). The bottom panel depicts the models' relative accuracy for $K_B = 0.07$ (solid line), 0.1 (dashed line), and 1 (dotted line).</figcaption></figure>

$$ \left(\frac{k}{aH}\right)^2 \Psi \propto a^q, \quad q-1 \propto p. \qquad (67) $$

Substituting this into Eq. (66) and remembering $\mu_0 \gg 1$, the solution for $\delta$ can be written as the sum of a power law, i.e. $\delta \propto a^q$, and a linear combination of solutions to the homogeneous system,

$$ \delta_1 \propto a^{-1/2}, \quad \delta_2 = \text{const} \qquad (68) $$

The approximation is expected to hold even slightly beyond the validity of the adopted first-order expansion since $\log \Psi$ effectively grows proportionally to $\log a$ over the time range of interest in cases where $p \le 0.05$.

The top panel of Fig. 5 illustrates the evolution of $\delta$ together with these analytic models (dashed lines) for $k = 0.5$ Mpc$^{-1}$, $\mu_0 = 250$, and different choices of $K_B \ge 0.07$. Starting from $a \sim 0.02$, the numerical results nicely tend toward the analytic solutions which were matched near the end of matter domination at $a = 0.2$. Their relative deviation is indicated in the figure's bottom panel and is in accordance with Sec. IV C 1, where similar domains of validity have been found. If $p \ll 1$, the constant $q$ approaches unity and one obtains the usual GR EdS growing mode $\delta \propto a$, which is the case for $K_B = 1$. The other cases exhibit enhanced growth and are characterized by power-law solutions with $q > 1$. Again, it is remarkable

that the complex dynamics of TeVeS allows a surprisingly simple, effective description of perturbation growth for a wide range of parameters.

For $\mu_0 K_B \sim 10$ and smaller, the simple power-law approximation will break down and one must resort to Eq. (60) when solving the growth equation, leading to analytic solutions for $\delta$ which involve exponential integrals. As is seen from the topmost line in Fig. 3, $\Psi$ grows faster than a power law in these cases, and one may expect a similar behavior for $\delta$.

### D. Superhorizon scales

On scales much larger than the horizon, terms proportional to $k^2$ may be safely neglected in the perturbation equations, i.e. $k \to 0$. Considering the Hamiltonian constraint and the evolution equation for $\delta$, one obtains

$$ -3 \frac{\dot{b}}{b} \left( \dot{\bar{\Phi}} + \frac{\dot{b}}{b} \tilde{\Psi} \right) + \frac{a e^{-\bar{\Phi}}}{2} \dot{\bar{\phi}}_{\gamma} = 4\pi G a^2 e^{-4\bar{\Phi}} \bar{\rho} (\delta - 2\varphi) \qquad (69) $$

and

$$ \dot{\delta} = 3\dot{\bar{\Phi}} = 3(\dot{\bar{\Phi}} + \dot{\varphi}). \qquad (70) $$

As usual, the above implies that the combination $\delta - 3\Phi$ is conserved over time.

Using the above equations together with the closure relations from Sec. IV B, it is possible to arrive at an equation governing the evolution of $\tilde{\Phi}$. The corresponding derivation is sketched in Appendix B 2, and the final expression takes the form

$$ \tilde{\Phi}'' + \frac{13}{10a} \tilde{\Phi}' = \frac{2}{15a^2} \left( \tilde{B}_{\gamma} + 2\hat{B}_{\gamma} - \frac{3}{2} \hat{B}_{\varphi} \right), \qquad (71) $$

where it is again assumed that $\mu_0 \gg 1$, and also $\tilde{\Phi} \approx \Phi$ since $\varphi \ll \tilde{\Phi}$ is typically satisfied. There is no explicit dependence on $K_B$ which can only enter through $B_\varphi$ and $B_\gamma$. Formally, the solution to Eq. (71) can be expressed in terms of the auxiliary functions on the right-hand side which are unknown a priori. The numerical analysis, however, suggests that the right-hand terms nearly cancel each other,

$$ \tilde{B}_{\gamma} + 2\hat{B}_{\gamma} - \frac{3}{2}\hat{B}_{\varphi} \approx 0. \qquad (72) $$

Considering the resulting homogeneous equation then leads to the solutions

$$ \Phi_1 \propto a^{-3/10}, \quad \Phi_2 = \text{const}, \qquad (73) $$
<figure><img src="image_9.png" /><figcaption>FIG. 6. Evolution of $\Phi$ for $k = 5 \times 10^{-4}$ Mpc$^{-1}$, $\mu_0 = 250$, and $K_B = 0.1$. The top panel shows the numerical result (solid line) and the analytic model described in the text (dashed line). The dotted lines indicate how $\Phi$ changes within the range $0.07 < K_B < 1$. The bottom panel depicts the relative accuracy of the analytic model.</figcaption></figure>

and one is left with the usual result that $\Phi_2$ and $\delta_2$ are frozen for modes which have not yet entered the horizon. Compared to the situation in GR, where $\Phi_1 \propto a^{-5/2}$, the decaying mode evolves much more slowly and may not die out fast enough to be discarded in the late matter era.

The top panel of Fig. 6 shows the evolution of the potential $\Phi$ together with the analytic model (matched near $a = 0.05$) based on Eq. (73) for $k = 5 \times 10^{-4}$ Mpc$^{-1}$, $\mu_0 = 250$, and $K_B = 0.1$. The numerical result is basically independent of $K_B$, as indicated by the dotted lines which quantify the change of $\Phi$ over the range $0.07 < K_B < 1$. The analytic model is very accurate, with a relative error less than approximately 2–3% throughout matter domination. As already suspected, the decaying mode is important even at late times around $a \sim 0.1$. Hence, $\Phi$ and $\delta$ are no longer conserved, but characterized by a rather slow decay outside the horizon. The insensitivity to $K_B$ on these large scales can also be understood from Eq. (44) which suggests that the impact of vector field perturbations becomes negligible as $k \to 0$.

## V. CONCLUSIONS

In this work, I have revisited the evolution of cosmological perturbations in Bekenstein's TeVeS theory. Considering only scalar modes in the conformal Newtonian gauge, I have introduced two auxiliary functions, $B_\phi$ and $B_\gamma$, which allow one to express perturbations

of the TeVeS scalar and vector fields in a way suitable for studying modifications due to these new degrees of freedom at the level of the metric potentials.

Assuming a universe in the matter-dominated phase, I have examined the theory's behavior on scales well inside and outside the horizon, respectively. Deriving approximate expressions for the potentials on subhorizon scales, I have adopted a semianalytic approach to identify and describe the mechanism responsible for the superlinear growth of density perturbations. An important element of this growth mechanism is the ratio between $\bar{\phi}$, the background scalar field, and the coupling constant $K_B$, triggering significant changes in the potentials if $K_B \lesssim 4|\bar{\phi}|$. The analysis confirms the implementation of the TeVeS vector field as the key to enhanced structure formation, in agreement with the result of Ref. [22]. The evolution of $\delta$, the density contrast, and the potential $\Psi$ during the matter era has been explored semianalytically and is well approximated in terms of simple power-law solutions for a wide range of parameters. This is striking in view of the theory's complexity and does warrant further investigation. In the limit of superhorizon scales, vector field perturbations effectively decouple from the equations. Unlike the situation in GR, the potential $\Phi$ and $\delta$ are not conserved outside the horizon, but characterized by a slowly decaying mode which remains important even at late times.

The approach taken here is not restricted to the framework of TeVeS, and might be useful to study the properties of other modified gravity theories with additional degrees of freedom.

## ACKNOWLEDGMENTS

The author likes to thank Jacob D. Bekenstein, Cosimo Fedeli, David F. Mota, Stephen E. Rafter, and the anonymous referees for useful discussions and comments which significantly improved the manuscript. The author is particularly grateful to Constantinos Skordis for help with TeVeS perturbations and his modified Boltzmann code. This work was partially supported by a fellowship from the Lady Davis Foundation and by the Grant Spin(e) ANR-13-BS05-0005 of the French National Research Agency. The author acknowledges support through a fellowship from the Minerva Foundation.

## APPENDIX A: SCALAR FIELD EVOLUTION DURING TRACKING

In the following, I will assume the generalized potential defined in Eq. (13) and adopt the notation and definitions used in Ref. [18]. There it has been found that the scalar field evolves during tracking as

$$ \bar{\phi} = \bar{\phi}_0 + \phi_1 \log a, \quad (A1) $$

where
$$ \phi_1 \equiv \frac{d\bar{\phi}}{d \log a} \qquad (A2) $$

is approximately constant. Following the derivation presented in Ref. [18], one shows that

$$ \frac{\phi_1}{1 + \phi_1} = \frac{\beta}{2\mu_0} \sqrt{\left(\frac{1 + 3w}{1 - w}\right)^2}, \qquad (A3) $$

where $\beta = \pm 1$ denotes the sign of the scalar field's time derivative, i.e.

$$ \beta \equiv \operatorname{sgn} \dot{\bar{\phi}}. \qquad (A4) $$

To see that the sign in Eq. (A3) is chosen appropriately, one uses Eq. (A1) and finds that

$$ \beta = \operatorname{sgn}\left(\phi_1 \frac{\dot{a}}{a}\right) = \operatorname{sgn}\phi_1 = \operatorname{sgn}\beta, \qquad (A5) $$

where I have assumed that $|\phi_1| \ll 1$ for the last equality. Note that this is justified because of the requirement $\mu_0 \gg 1$ for viable cosmological models. When evaluating the square root in Eq. (A3), one needs to take into account that the argument's sign depends on the actual choice of $w$. Therefore, one has

$$ \sqrt{\left(\frac{1 + 3w}{1 - w}\right)^2} = \frac{|1 + 3w|}{|1 - w|}, \qquad (A6) $$

which eventually gives the result in Eq. (24). During tracking, the field $\bar{\mu}$ (see Sec. III B) evolves as $\bar{\mu} = 2\mu_0(1 + \epsilon)$, where

$$ \log \epsilon \propto - \frac{2\phi_1 + 3(1+w)}{n} \log a, \qquad (A7) $$

and thus $2\phi_1 + 3(1+w) > 0$ emerges as a condition for stable tracking. For a universe dominated by a cosmological constant $\Lambda$, one has $w = -1$ and therefore $\beta = 1$.

Since the time derivative of $\bar{\phi}$ changes its sign when passing from the matter to the $\Lambda$ era (resulting in $\bar{\rho}_{\phi}$ momentarily going to zero) [21], it follows that $\beta = -1$ during matter domination. This result is in accordance with previous work [4,22] where it has been shown that $\bar{\phi}$ decreases with time during the matter era.

## APPENDIX B: PARAMETRIZED PERTURBATION EQUATIONS

In what follows, I will adopt the linear TeVeS perturbation equations for scalar modes in the conformal Newtonian gauge [9]. Further, I will assume a spatially flat spacetime and a universe filled with pressureless matter only, corresponding to the modified EdS cosmology introduced in Sec. III C. In this case, the fluid's pressure components may be neglected ($w = C_s = \Sigma = 0$, where $C_s$ and $\Sigma$ are the fluid's sound speed and shear, respectively) and the background density evolves as $\bar{\rho} \propto a^{-3}$. As usual, the equations are expressed in Fourier space using the conformal wave vector $\mathbf{k}$ in accordance with the coordinate system specified in Sec. IVA.

### 1. General case

The goal is to express the metric potentials in terms of matter fluid variables, using the closure relations from Sec. IV B. Since there remain three fields, i.e. the two metric potentials and the vector field perturbation $\alpha$, this requires finding three linearly independent equations. For the first equation, I eliminate $\dot{\alpha}$ between the corresponding vector field and propagation equations, which leads to

$$ \tilde{\Phi} + (1 - e^{4\bar{\phi}})E - e^{4\bar{\phi}}\tilde{\Psi} - 4\dot{\bar{\phi}}\alpha + e^{4\bar{\phi}}\left(\frac{\dot{a}}{a} + 5\dot{\bar{\phi}}\right)\tilde{\zeta} = 0. \qquad (B1) $$

Differentiating the above and substituting all remaining time derivatives by a suitable combination of the perturbation equations eventually gives

$$ \begin{aligned} & e^{4\bar{\phi}}k^2\tilde{\zeta} + \left[ \frac{(1 - e^{4\bar{\phi}})^2}{K_B} \left( \frac{\bar{\mu}\dot{\bar{\phi}}^2}{1 - e^{4\bar{\phi}}} + 8\pi G a^2 e^{-4\bar{\phi}}\bar{\rho} \right) - 8 \left( 1 + 2 \frac{1 + e^{4\bar{\phi}}}{1 - e^{4\bar{\phi}}} \right) \dot{\bar{\phi}}^2 - 4\ddot{\bar{\phi}} \right] \alpha + 4\dot{\bar{\phi}} \left( -1 + \frac{1 + e^{4\bar{\phi}}}{1 - e^{4\bar{\phi}}} \right) \tilde{\Phi} \\ & + e^{4\bar{\phi}} \left[ 4 \frac{1 + e^{4\bar{\phi}}}{1 - e^{4\bar{\phi}}} \dot{\bar{\phi}} \left( \frac{\dot{a}}{a} + 5\dot{\bar{\phi}} \right) + \frac{\ddot{a}}{a} - 2 \frac{\dot{a}^2}{a^2} - 4\dot{\bar{\phi}} \frac{\dot{a}}{a} + 5\dot{\bar{\phi}}^2 + 5\ddot{\bar{\phi}} \right] \tilde{\zeta} + 4\pi G a^2 e^{-4\bar{\phi}} \bar{\rho} \left[ 1 + 3e^{4\bar{\phi}} - 2 \frac{(1 - e^{4\bar{\phi}})^2}{K_B} \right] \theta - e^{4\bar{\phi}} \frac{\dot{a}}{a} \hat{B}_{\phi} \\ & + \left[ \left( 1 + 3e^{4\bar{\phi}} - \frac{1 - e^{4\bar{\phi}}}{K_B} \right) \bar{\mu} \dot{\bar{\phi}} - e^{4\bar{\phi}} \left( 4 \frac{\dot{b}}{b} + 5\dot{\bar{\phi}} + 4 \frac{1 + e^{4\bar{\phi}}}{1 - e^{4\bar{\phi}}} \dot{\bar{\phi}} \right) \right] \tilde{\Psi} - \bar{\mu} \dot{\bar{\phi}} \left( 1 + 3e^{4\bar{\phi}} - \frac{1 - e^{4\bar{\phi}}}{K_B} \right) \tilde{B}_{\phi} - e^{4\bar{\phi}} \frac{\dot{a}}{a} \tilde{B}_{\gamma} = 0. \end{aligned} \qquad (B2) $$
where $\theta$ denotes the fluid's velocity potential, $k = |\mathbf{k}|$, and $\tilde{\zeta}$ is related to $\alpha$ through Eq. (41). To find a second equation, one may start from Eq. (44). Similar as before, I take its time derivative and use the perturbation equations to recast the resulting expression into a more convenient form. A bit of algebra then reveals

$$ -e^{-4\bar{\phi}}k^2\tilde{\Phi} + e^{-4\bar{\phi}}\left(1 - \frac{\bar{\mu}}{2U}\right)k^2\tilde{\Psi} - e^{-4\bar{\phi}}\frac{\bar{\mu}}{2U}\dot{\bar{\phi}}k^2\alpha + \left[\frac{\bar{\mu}}{U}\dot{\bar{\phi}} - 3\left(\bar{\mu}\dot{\bar{\phi}} - \frac{\dot{b}}{b}\right) + \frac{\dot{U}}{U}\right](k^2\tilde{\zeta} + 12\pi G a^2 e^{-4\bar{\phi}}\bar{\rho}\theta) \\ + \left[3\left(\frac{\dot{b}}{b} - 5\dot{\bar{\phi}}\right)\frac{\dot{b}}{b} + \frac{4\pi}{U}Ga^2e^{-4\bar{\phi}}\bar{\rho} - 6\frac{\dot{b}}{b} + 3\left(\bar{\mu}\dot{\bar{\phi}} - \frac{\dot{b}}{b}\right)\left(\frac{\bar{\mu}}{U}\dot{\bar{\phi}} - 3\left(\bar{\mu}\dot{\bar{\phi}} - \frac{\dot{b}}{b}\right) + \frac{\dot{U}}{U}\right)\right]\tilde{\Psi} - \frac{4\pi}{U}Ga^2e^{-4\bar{\phi}}\bar{\rho}\delta - 3\frac{\dot{a}\dot{b}}{a b}\hat{B}_{\varphi} \\ - \frac{\dot{a}^2}{a^2}\hat{B}_{\gamma} + \left[\frac{e^{-4\bar{\phi}}}{2U}(\bar{\mu}k^2 - 16\pi G a^2\bar{\rho}) - 3\bar{\mu}\dot{\bar{\phi}}\left(\frac{\bar{\mu}}{U}\dot{\bar{\phi}} - 3\left(\bar{\mu}\dot{\bar{\phi}} - \frac{\dot{b}}{b}\right) + \frac{\dot{U}}{U}\right)\right]\tilde{B}_{\varphi} + \left(3\bar{\mu}\dot{\bar{\phi}} - 5\frac{\dot{b}}{b} - 2\dot{\bar{\phi}} - \frac{\dot{U}}{U}\right)\frac{\dot{a}}{a}\tilde{B}_{\gamma} = 0. \quad (B3) $$

Finally, the last equation is obtained from eliminating $E$ between Eq. (B1) and the Hamiltonian constraint equation. Together with the relations presented in Sec. IV B, one ends up with

$$ -\left(2 - \frac{K_B}{1 - e^{4\bar{\phi}}}\right)k^2\tilde{\Phi} - K_B \frac{e^{4\bar{\phi}}}{1 - e^{4\bar{\phi}}}k^2\tilde{\Psi} + K_B \left(\frac{\dot{b}}{b} - 4\frac{e^{4\bar{\phi}}}{1 - e^{4\bar{\phi}}}\dot{\bar{\phi}}\right)k^2\alpha - 2e^{4\bar{\phi}}\left(\frac{\dot{b}}{b} - U\dot{\bar{\phi}}\right)(k^2\tilde{\zeta} + 12\pi G a^2 e^{-4\bar{\phi}}\bar{\rho}\theta) - 8\pi G a^2 \bar{\rho}\delta \\ - 2e^{4\bar{\phi}}\left[3(\bar{\mu} + U)\frac{\dot{b}}{b}\dot{\bar{\phi}} - 3\bar{\mu}U\dot{\bar{\phi}}^2 - 8\pi G a^2 e^{-4\bar{\phi}}\bar{\rho}\right]\tilde{\Psi} + 2e^{4\bar{\phi}}\left[3\bar{\mu}\dot{\bar{\phi}}\left(\frac{\dot{b}}{b} - U\dot{\bar{\phi}}\right) - 8\pi G a^2 e^{-4\bar{\phi}}\bar{\rho}\right]\tilde{B}_{\varphi} - 2e^{4\bar{\phi}}U\dot{\bar{\phi}}\frac{\dot{a}}{a}\tilde{B}_{\gamma} = 0. \quad (B4) $$

## 2. Superhorizon limit

For scales well outside the horizon, the perturbation equations can be simplified by considering the limit $k \to 0$. Taking the time derivative of Eq. (44),

$$ 3\ddot{\tilde{\Phi}} - \frac{\dot{a}^2}{a^2}\hat{B}_{\gamma} - \left(\frac{\dot{a}}{a} - \frac{\dot{U}}{U} - \dot{\bar{\phi}}\right)\left(3\dot{\bar{\Phi}} - \frac{\dot{a}}{a}\tilde{B}_{\gamma}\right) = \frac{ae^{-\bar{\phi}}}{2U}\dot{\gamma}, \quad (B5) $$

the first step is to substitute $\dot{\gamma}$ with the help of the perturbed scalar field equations. Next, one needs to eliminate terms involving the potential $\tilde{\Psi}$. This can be achieved by appropriately combining the time derivatives of the modified propagation equations with the remaining perturbation equations. Finally, using Eqs. (69) and Eq. (70) allows one to derive an equation governing the evolution of $\tilde{\Phi}$.

Defining the auxiliary quantities

$$ \tilde{P} \equiv 3 \frac{\dot{b}^2}{b^2} - 4\pi G a^2 e^{-4\bar{\phi}}\bar{\rho} \quad (B6) $$

and

$$ \tilde{Q} \equiv U \left[ 2 \frac{\ddot{b}}{b} + \frac{\dot{b}}{b} \left( 5 \dot{\bar{\phi}} - \frac{\dot{b}}{b} \right) \right], \quad (B7) $$

the result can be expressed as

$$ (\tilde{P} - 3\tilde{Q})\ddot{\tilde{\Phi}} + \tilde{P}\frac{\dot{a}\dot{b}}{a b}\hat{B}_{\varphi} + \left[\tilde{P}\left((3\bar{\mu} + 2)\dot{\bar{\phi}} - \frac{\dot{b}}{b}\right) - 3\tilde{Q}\right. \\ \times \left((2 + U^{-1})\frac{\dot{b}}{b} + \frac{\dot{U}}{U} + (1 + \bar{\mu}U^{-1})\dot{\bar{\phi}}\right)]\dot{\tilde{\Phi}} + \tilde{Q}\frac{\dot{a}^2}{a^2}\hat{B}_{\gamma} \\ \left. + \frac{\dot{a}}{a}\left[\tilde{Q}\left(2\frac{\dot{b}}{b} + \dot{\bar{\phi}} + \frac{\dot{U}}{U}\right) + \tilde{P}\left(\frac{\dot{b}}{b} - \bar{\mu}\dot{\bar{\phi}}\right)\right]\right]\tilde{B}_{\gamma} = 0. \quad (B8) $$

Adopting the modified EdS cosmology, this eventually leads to the expression in Eq. (71).
[1] M. Tegmark, M. R. Blanton, M. A. Strauss, F. Hoyle, D. Schlegel, R. Scoccimarro, M. S. Vogeley, D. H. Weinberg, I. Zehavi, A. Berlind, T. Budavari, A. Connolly, D. J. Eisenstein, D. Finkbeiner, J. A. Frieman et al., Astrophys. J. 606, 702 (2004); J. Benjamin, C. Heymans, E. Semboloni, L. van Waerbeke, H. Hoekstra, T. Erben, M. D. Gladders, M. Hetterscheidt, Y. Mellier, and H. K. C. Yee, Mon. Not. R. Astron. Soc. 381, 702 (2007); J. Guy, M. Sullivan, A. Conley, N. Regnault, P. Astier, C. Balland, S. Basa, R. G. Carlberg, D. Fouchez, D. Hardin et al., Astron. Astrophys. 523, A7 (2010); P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, J. G. Bartlett et al. (Planck Collaboration), arXiv: 1502.01589.

[2] G. Bertone, D. Hooper, and J. Silk, Phys. Rep. 405, 279 (2005); J. A. Frieman, M. S. Turner, and D. Huterer, Annu. Rev. Astron. Astrophys. 46, 385 (2008).

[3] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Phys. Rep. 513, 1 (2012).

[4] J. D. Bekenstein, Phys. Rev. D 70, 083509 (2004).

[5] M. Milgrom, Astrophys. J. 270, 384 (1983); 270, 365 (1983); J. Bekenstein and M. Milgrom, Astrophys. J. 286, 7 (1984).

[6] B. Famaey and S. S. McGaugh, Living Rev. Relativity 15, 10 (2012).

[7] E. Bertschinger, Astrophys. J. 648, 797 (2006); P. G. Ferreira and C. Skordis, Phys. Rev. D 81, 104020 (2010); E. Bertschinger, Phil. Trans. R. Soc. A 369, 4947 (2011).

[8] T. Baker, P. G. Ferreira, C. Skordis, and J. Zuntz, Phys. Rev. D 84, 124018 (2011); J. Zuntz, T. Baker, P. G. Ferreira, and C. Skordis, J. Cosmol. Astropart. Phys. 6 (2012) 032; T. Baker, Phys. Rev. D 85, 044020 (2012).

[9] C. Skordis, Phys. Rev. D 74, 103513 (2006).

[10] C. Skordis, Classical Quantum Gravity 26, 143001 (2009).

[11] T. G. Zlosnik, P. G. Ferreira, and G. D. Starkman, Phys. Rev. D 74, 044037 (2006).

[12] T. G. Zlosnik, P. G. Ferreira, and G. D. Starkman, Phys. Rev. D 75, 044017 (2007).

[13] Note the change of notation: Here $\mu$ is related to Bekenstein's auxiliary scalar field $\sigma_B$ through $\mu = 8\pi G\sigma_B^2$, and the coupling constant of the scalar field, i.e. Bekenstein's $k_B$, is redefined in terms of $\mu_0$ using $\mu_0 = 8\pi/k_B$.

[14] For $\phi_C \neq 0$, the metric $g_{\mu\nu}$ is not asymptotically Minkowskian. However, this is easily remedied by an appropriate rescaling of coordinates.

[15] J. D. Bekenstein and E. Sagi, Phys. Rev. D 77, 103512 (2008).

[16] H. S. Zhao and B. Famaey, Astrophys. J. 638, L9 (2006).

[17] H. Zhao, Int. J. Mod. Phys. D 16, 2055 (2007).

[18] F. Bourliot, P. G. Ferreira, D. F. Mota, and C. Skordis, Phys. Rev. D 75, 063508 (2007).

[19] This class of potentials will modify the dynamics of quasistatic systems if $n \neq 2$, but this can be avoided by resorting to a piecewise definition of $V(\mu)$.

[20] B. Ratra and P. J. E. Peebles, Phys. Rev. D 37, 3406 (1988); P. G. Ferreira and M. Joyce, Phys. Rev. Lett. 79, 4740 (1997); C. Skordis and A. Albrecht, Phys. Rev. D 66, 043523 (2002).

[21] C. Skordis, D. F. Mota, P. G. Ferreira, and C. Bøhm, Phys. Rev. Lett. 96, 011301 (2006).

[22] S. Dodelson and M. Liguori, Phys. Rev. Lett. 97, 231301 (2006).

[23] P. Zhang, M. Liguori, R. Bean, and S. Dodelson, Phys. Rev. Lett. 99, 141302 (2007).

[24] The value of $\mu_0$ still matters for the growth of vector perturbations and enters through the auxiliary functions and $\bar{\phi}$.