# Relativistic Khronon Theory in agreement with Modified Newtonian Dynamics and Large-Scale Cosmology

Luc Blanchet, Constantinos Skordis

► To cite this version:

Luc Blanchet, Constantinos Skordis. Relativistic Khronon Theory in agreement with Modified Newtonian Dynamics and Large-Scale Cosmology. Journal of Cosmology and Astroparticle Physics, 2024, 11, pp.040. 10.1088/1475-7516/2024/11/040. [hal-04576609]

HAL Id: hal-04576609

https://hal.science/hal-04576609v1

Submitted on 17 Jul 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

PAPER • OPEN ACCESS

# Relativistic Khronon theory in agreement with modified Newtonian dynamics and large-scale cosmology

To cite this article: Luc Blanchet and Constantinos Skordis JCAP11(2024)040

View the [article online](#) for updates and enhancements.

You may also like

- Anti-Newtonian dynamics and self-induced Bloch oscillations of correlated particles
  Stefano Longhi

- Khronometric Theories of Modified Newtonian Dynamics
  Éanna. É. Flanagan

- Gas Expulsion in MOND: The Possible Origin of Diffuse Globular Clusters and Ultra-faint Dwarf Galaxies
  Xufen Wu and Pavel Kroupa
# Relativistic Khronon theory in agreement with modified Newtonian dynamics and large-scale cosmology

Luc Blanchet ID a and Constantinos Skordis ID b,c

<sup>a</sup>GReCO, Institut d'Astrophysique de Paris, UMR 7095, CNRS, Sorbonne Université,
98<sup>bis</sup> boulevard Arago, 75014 Paris, France

<sup>b</sup>CEICO, Institute of Physics of the Czech Academy of Sciences,
Na Slovance 1999/2, 182 21, Prague, Czech Republic

<sup>c</sup>Department of Physics, University of Oxford,
Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, U.K.

E-mail: luc.blanchet@iap.fr, skordis@fzu.cz

**ABSTRACT:** We propose an extension of General Relativity (GR) based on a space-time foliation by three-dimensional space-like hypersurfaces labeled by the Khronon scalar field $\tau$. We show that this theory (i) leads to modified Newtonian dynamics (MOND) at galactic scales for stationary systems; (ii) recovers GR plus a cosmological constant in the strong field regime; (iii) is in agreement with the standard cosmological model and the observed cosmic microwave background anisotropies at linear cosmological scales, where the theory reduces to a subset of the generalized dark matter (GDM) model. We compute the second order action on a Minkowski background and show that it contains the usual tensor modes of GR and a scalar degree of freedom with dispersion relation $\omega = 0$. We find that the deconstrained Hamiltonian is bounded from below for wavenumbers larger than $\sim 10^{-31}$ eV and unbounded for smaller wavenumbers.

**KEYWORDS:** gravity, modified gravity

ARXIV ePRINT: 2404.06584

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

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Action and field equations</td><td>3</td></tr><tr><td>2.1 General definitions</td><td>3</td></tr><tr><td>2.2 Covariant action and field equations</td><td>4</td></tr><tr><td>3 Nonrelativistic limit</td><td>6</td></tr><tr><td>3.1 Post Newtonian expansions for the fields</td><td>6</td></tr><tr><td>3.2 PN expansion of the field equations</td><td>8</td></tr><tr><td>4 Cosmology</td><td>12</td></tr><tr><td>4.1 Friedmann-Lemaître-Robertson-Walker (FLRW) cosmology</td><td>12</td></tr><tr><td>4.2 Linear perturbations on a FLRW background</td><td>15</td></tr><tr><td>4.3 Khronon cosmology and the MOND limit</td><td>19</td></tr><tr><td>5 Linear stability on Minkowski space</td><td>23</td></tr><tr><td>5.1 The normal modes</td><td>23</td></tr><tr><td>5.2 The Hamiltonian formulation</td><td>25</td></tr><tr><td>6 Discussions</td><td>27</td></tr><tr><td>6.1 The Hořava khronometric theory</td><td>27</td></tr><tr><td>6.2 Degrees of freedom and further discussion about the dispersion relation</td><td>29</td></tr><tr><td>7 Conclusion</td><td>30</td></tr><tr><td>A Formulation in adapted coordinates</td><td>31</td></tr></table>

## 1 Introduction

More than forty years ago Milgrom [1–3] proposed an empirical formula, dubbed MOND for “MOdified Newtonian Dynamics”, able to fit with astonishing precision the rotation curves of galaxies and the Tully-Fisher relation between the baryonic mass and the asymptotic rotation velocity of galaxies. The MOND formula permits to resolve a number of challenges faced by the standard cosmological model Λ-CDM (Cold Dark Matter plus a cosmological constant Λ) at galactic scales [4]. It suggests a modification of fundamental physics, notably the gravitational sector described by general relativity (GR), in a regime of low accelerations, below a certain critical acceleration scale $a_0$, empirically measured at the value $a_0 \simeq 1.210^{-10} \text{ m/s}^2$.

The MOND formula has prompted the construction of several extensions of GR [5, 6], which purport to fit cosmological observations without the need of dark matter. A well-studied example is the Tensor-Vector-Scalar (TeVES) theory of Bekenstein and Sanders [7–9], which extends GR with a time-like vector field (akin to Einstein-Aether theory [10, 11]) and a scalar field; see [12] for a review. Einstein-Aether theories, originally motivated by the
phenomenology of local Lorentz invariance violation [10, 11], were extended to lead to the MOND formula at the regime of galaxies [13, 14]. Further GR extensions with a MOND limit were proposed in [15–28] and the cosmological and astrophysical predictions of several of these theories have been extensively studied [29–38]. The recent Aether-Scalar-Tensor (AeST) theory proposed by Skordis and Zlosnik [39] based on [40], has been the first example of a GR extension reproducing the MOND formula in galaxies and simultaneously being in agreement with the standard cosmological model Λ-CDM, i.e. predicting the full spectrum of anisotropies of the cosmic microwave background (CMB).<fnref n="1" />

A simple extension of GR leading to MOND behaviour in galaxies was proposed by Blanchet and Marsat [18], hereafter BM theory, as well as an independent variant by Sanders [19], based on only two dynamical fields, the metric $g_{\mu\nu}$ and a scalar field $\tau$ called the Khronon, and labeling a family of space-like three-dimensional hypersurfaces. All ordinary matter fields in this theory are universally coupled to the metric. However, in the local freely falling frame associated with the matter fields, the presence of the Khronon field induces an effective violation of the local Lorentz invariance (LLI). The BM theory was inspired by the Hořava-Lifshitz approach [43–45] for a possible violation of LLI and a completion of GR in a (power-counting) renormalizable theory at high energy. However, its motivation is quite different from that of [43–45], as it is not concerned with the problem of quantum gravity, and postulates the violation of LLI at low energy, in the MOND regime for weak gravitational fields. This idea has been generalized by Bonetti and Barausse [46] and the BM theory has recently been further discussed by Flanagan [47].

The main drawback of the BM theory, that concerns us here, is that it does not explain the effects commonly attributed to dark matter, which occur at large scales in cosmological perturbations. We therefore extend the BM theory by adding, in a quite natural way, a kinetic term for the Khronon scalar field. Following [48–50] we know that a large class of kinetic $k$-essence terms are able to mimic dark matter in perturbations around an homogeneous and isotropic Friedman-Lemaître-Robertson-Walker (FLRW) background. Hence, we shall prove that with this crucial addition the model, for an appropriate choice of the Khronon kinetic term, is in agreement with the standard cosmological model at large scales (to first order in cosmological perturbation), and in particular retrieves the full observed spectrum of anisotropies of the CMB, while still providing the relevant MOND limit.

More precisely, we show that the Khronon equations can be recast into the framework of the generalized dark matter (GDM) model [51, 52]. The Khronon behaves like a GDM fluid with time-dependent equation of state in the background, and with zero viscosity and non-zero, $k$-dependent sound speed in linear perturbations. Using the GDM framework we investigate the allowed constraints on cosmology and the CMB (namely the cosmological dark matter which should be close enough to a “pure dust” model) which are consistent with recovering the MOND phenomenology at galactic scales. In particular we find that a kinetic term for the Khronon based on the Dirac-Born-Infeld (DBI) form of the action reproduces in a natural way all observations.

<footnote n="1">¹A different class of MOND theories attempt to modify the properties of dark matter itself, rather than the gravitational law. An example is the dipolar dark matter model [41, 42] which is also in agreement with the standard cosmological model Λ-CDM in first order cosmological perturbations.</footnote>
Finally, we also show that the theory recovers GR plus a cosmological constant in the strong field, high acceleration regime, and in particular has the same parametrized post-Newtonian (PPN) parameters as GR for tests in the Solar System.

The paper is organized as follows. In section 2 we describe the covariant action of the theory and derive the field equations. We discuss in appendix A an equivalent formulation of the theory in the coordinate system adapted to the Khronon, that is, when the time coordinate is identified with the Khronon field itself. In section 3 we construct the nonrelativistic slow motion limit of the theory and discuss how Newtonian behaviour at high-acceleration (and GR in the strong-field high-acceleration limit) is reached and how MOND behaviour at the low-acceleration limit is obtained. In section 4 we show how some specific choices for the new term in the action that we propose in the present work leads to cosmology in agreement with observations of large scale structure and the CMB radiation. In section 5 we expand the action to second order in perturbations on a Minkowski background, determine the normal modes and compute the deconstrained Hamiltonian (in particular we find that the theory has the same tensorial modes as GR). Section 6 ends with a short discussion of the Hořava khronometric theory and presents our conclusions.

We use a (–+++) metric signature convention and curvature convention of MTW [53]. We denote general coordinates by $x^\mu$ which are taken to have dimensions of length, such that, for a time coordinate $t$ we have $x_0 = ct$ with $c$ being the speed of light. Partial $\partial_\mu$ and covariant derivatives $\nabla_\mu$ are with respect to $x^\mu$ and have dimensions of inverse length. We denote by $\gamma_{ij}$ the Euclidean metric in general coordinates (reserving $\delta_{ij}$ for the Kronecker symbol in cartesian coordinates), and by $\vec{\nabla}_i$ the covariant derivative associated with $\gamma_{ij}$. We use the vector calculus notation $\vec{\nabla}_i \leftrightarrow \vec{\nabla}$ and $\vec{\nabla}u \cdot \vec{\nabla}v \equiv \vec{\nabla}_i u \vec{\nabla}^i v = \gamma^{ij} \vec{\nabla}_i u \vec{\nabla}_j v$ for any scalars $u$ and $v$. In cosmology we denote by $\gamma_{ij}^\kappa$ a spatially homogeneous and isotropic metric of spatial curvature $\kappa$, hence $\kappa = 0$ corresponds to a spatially flat Universe, so that $\gamma_{ij}^0 = \gamma_{ij}$, and $\kappa > 0$ ($\kappa < 0$) to a positively (negatively) spatially curved Universe. We have $\gamma_{ij}^\kappa = \gamma_{ij}(1 + \frac{\kappa}{4}|\vec{x}|^2)^{-2}$ in a spherical coordinate system.

## 2 Action and field equations

### 2.1 General definitions

The dynamical degrees of freedom of the theory comprise the metric $g_{\mu\nu}$, the Khronon scalar field $\tau$ which is assumed to have units of time, and the ordinary matter fields $\Psi$. The Khronon defines a spatial foliation through the orthonormal unit-timelike vector field $n^\mu$ defined via

$$n_{\mu} = -\frac{c}{Q} \nabla_{\mu} \tau, \qquad (2.1)$$

with the useful scalar quantity

$$Q \equiv c \sqrt{-g^{\mu\nu} \nabla_{\mu} \tau \nabla_{\nu} \tau}. \qquad (2.2)$$

Both $n_\mu$ and $Q$ are dimensionless. By construction one can see that $n^\mu n_\mu = -1$ and that $n^\mu$ is future pointing, $n^0 > 0$. The vector $n^\mu$ serves to define the (covariant) spacelike
acceleration of the congruence as<fnref n="2" />

$$A_{\mu} = c^2 n^{\nu} \nabla_{\nu} n_{\mu} = -c^2 q_{\mu}^{\nu} \nabla_{\nu} \ln \mathcal{Q}, \qquad (2.3)$$

where the factor of $c^2$ ensures that it has units of acceleration, and where the projector onto the family of spatial hypersurfaces is denoted as

$$q_{\mu}^{\nu} \equiv \delta_{\mu}^{\nu} + n_{\mu} n^{\nu}, \qquad (2.4)$$

so that $n^{\mu} A_{\mu} = 0$. Note that the spatial foliation (2.1)-(2.3) is invariant by arbitrary reparametrization of the Khronon field: $\tau \to \tau' = F(\tau)$. However the kinetic term for the Khronon, eq. (2.7) below, will break the reparametrization invariance of the model.

Both in the BM theory [18] and in the Sanders variant [19], the acceleration (2.3) is used to lead to MOND behaviour in the nonrelativistic and low acceleration regime. This is achieved through a function $\mathcal{J}(\mathcal{Y})$ of the covariant acceleration (2.3) squared denoted

$$\mathcal{Y} \equiv \frac{A_{\mu} A^{\mu}}{c^4}, \qquad (2.5)$$

that appears in the action functional of the theory. The function $\mathcal{J}$ has the dimension of an inverse squared length and is chosen such that when $\mathcal{Y}$ is large (the high acceleration limit) GR is recovered, while when $\mathcal{Y}$ is small (the low acceleration limit) and in addition we assume slow-motion quasi stationary systems, MOND behaviour is reproduced. We discuss the details of how this happens, as well as choices of $\mathcal{J}(\mathcal{Y})$ which have this property, in section 3.

## 2.2 Covariant action and field equations

The covariant action of the theory is

$$S = \frac{c^3}{16\pi G} \int d^4 x \sqrt{-g} [R - 2\mathcal{J}(\mathcal{Y}) + 2\mathcal{K}(\mathcal{Q})] + S_m[\Psi, g], \qquad (2.6)$$

where $R$ is the four-dimensional Ricci scalar, $g$ is the metric determinant and $G$ the gravitational coupling strength. The matter fields $\Psi$ are universally coupled to the metric $g_{\mu\nu}$ and as such the Einstein equivalence principle holds (we regard the Khronon field as part of the gravitational sector). Moreover, in the high acceleration limit GR is recovered and as such, the strong equivalence principle (SEP) holds there, however, it does not hold at the low acceleration limit which includes MOND and the $\mu$-domination regime.

The function $\mathcal{J}(\mathcal{Y})$ appearing in (2.6), where $\mathcal{Y}$ is defined by (2.5), has already been postulated in [18] and [19] and is part of the original BM theory. Our new addition is the function $\mathcal{K}(\mathcal{Q})$ which depends on the scalar field $\mathcal{Q}$ defined in (2.1) and serves as a kinetic term for the Khronon field $\tau$. Following [39], our main assumption about this function is that it admits a Taylor expansion around the value $\mathcal{Q} = 1$. For definiteness we specifically adopt in this section

$$\mathcal{K}(\mathcal{Q}) = \mu^2 (\mathcal{Q} - 1)^2, \qquad (2.7)$$

<footnote n="2">²In [18, 19] the acceleration is denoted $a^{\mu}$; however, here we reserve $a$ for the scale factor in cosmology.</footnote>
and leave more general functions having such a Taylor expansion investigated in section 4.3. Here the constant $\mu$ has dimensions of inverse length, and such a term in the action plays a vital role in the success of AeST theory in fitting cosmological observations [39]. A similar term appears in the case of ghost condensation, which is an extension of GR associated with the breaking of time diffeomorphisms [49] and in the case of shift-symmetric $k$-essence [48]. In both of these last two cases the relevant term appearing in the function $\mathcal{K}$ is, rather, $\mathcal{K}(X) = \mu^2(X+1)^2$ with $X = c^2\nabla_\mu\phi\nabla^\mu\phi$ for a generic scalar $\phi$ (using our metric signature conventions), $\phi = \tau$ and $X = -Q^2$ in our case.

The term $\mathcal{K}(Q) = \mathcal{K}(\sqrt{-X})$ in the action (2.6) belongs to the general class of Horndeski theories [54], as a particular case of the $k$-essence term $L_2^H \equiv G_2(\phi, X)$. The acceleration squared (2.5) is part of the more general class of DHOST theories, as we have

$$ \mathcal{Y} = \frac{L_4^{(2)}}{X^2} - \frac{L_5^{(3)}}{X^3}, \qquad (2.8) $$

where $L_4^{(2)} \equiv \phi^\nu\phi_{\mu\nu}\phi^\rho\phi_\rho^\mu$ and $L_5^{(3)} \equiv (\phi^\mu\phi^\nu\phi_{\mu\nu})^2$ with standard notations [55, 56]. However the action (2.6) contains a function of $\mathcal{Y}$, involving higher powers of $\mathcal{Y}$, see (3.23) below, and the theory lies outside the class of DHOST theories.

By varying the action with respect to the metric we obtain the following generalized Einstein field equations

$$ G^{\mu\nu} + (\mathcal{J} - \mathcal{K}) g^{\mu\nu} - \frac{2}{c^4} \mathcal{J} \mathcal{Y} A^\mu A^\nu + \left[ \frac{2}{c^2} \nabla_\rho (\mathcal{J} \mathcal{Y} A^\rho) - \mathcal{Q} \mathcal{K}_\mathcal{Q} \right] n^\mu n^\nu = \frac{8\pi G}{c^4} T^{\mu\nu}, \quad (2.9) $$

where $G^{\mu\nu} = R^{\mu\nu} - \frac{1}{2}R g^{\mu\nu}$ is the Einstein tensor, $T^{\mu\nu} = \frac{2}{\sqrt{-g}} \delta S_m / \delta g_{\mu\nu}$ is the matter stress-energy tensor, and where $\mathcal{J} \equiv d\mathcal{J}/d\mathcal{Y}$ and $\mathcal{K}_\mathcal{Q} \equiv d\mathcal{K}/d\mathcal{Q}$ denote the derivatives with respect to the arguments $\mathcal{Y}$ and $\mathcal{Q}$. Taking the space-time trace of eq. (2.9) we obtain

$$ -R + 4\mathcal{J} - 2\mathcal{Y}\mathcal{J}\mathcal{Y} - 4\mathcal{K} + \mathcal{Q}\mathcal{K}_\mathcal{Q} - \frac{2}{c^2}\nabla_\mu(\mathcal{J}\mathcal{Y}A^\mu) = \frac{8\pi G}{c^4}T. \qquad (2.10) $$

We define the stress-energy tensor $\mathcal{T}^{\mu\nu}$ of the Khronon field to be such that (by definition)

$$ G^{\mu\nu} = \frac{8\pi G}{c^4} [T^{\mu\nu} + \mathcal{T}^{\mu\nu}], \qquad (2.11) $$

and we read directly from eq. (2.9)

$$ \mathcal{T}^{\mu\nu} = \frac{c^4}{8\pi G} \left\{ -(\mathcal{J} - \mathcal{K}) g^{\mu\nu} + \frac{2}{c^4} \mathcal{J} \mathcal{Y} A^\mu A^\nu + \left[ -\frac{2}{c^2} \nabla_\rho (\mathcal{J} \mathcal{Y} A^\rho) + \mathcal{Q} \mathcal{K}_\mathcal{Q} \right] n^\mu n^\nu \right\}. \quad (2.12) $$

From the contracted Bianchi identity $\nabla_\nu G^{\mu\nu} \equiv 0$ and the equations of motion satisfied by the matter fields $\Psi$, which imply the covariant conservation of the matter tensor $T^{\mu\nu}$, we infer that the stress-energy tensor of the Khronon must also be covariantly conserved,

$$ \nabla_\nu \mathcal{T}^{\mu\nu} = 0. \qquad (2.13) $$

Finally we vary with respect to the Khronon field and obtain the covariant conservation law $\nabla_\mu S^\mu = 0$ where the current vector $S^\mu$ reads as<fnref n="3" />

$$ S^\mu = \frac{2}{\mathcal{Q}} \mathcal{J} \mathcal{Y} A^\mu n^\nu \nabla_\nu \ln \mathcal{Q} + \left[ \frac{2}{\mathcal{Q}} \nabla_\nu (\mathcal{J} \mathcal{Y} A^\nu) - c^2 \mathcal{K}_\mathcal{Q} \right] n^\mu. \qquad (2.14) $$

<footnote n="3">³The Khronon equation $\nabla_\mu S^\mu = 0$ reduces to eq. (2.8) in the BM model [18] when $\mathcal{K} = 0$.</footnote>
As a check of the consistency of the formalism, we can derive the conservation equation of the Khronon stress-energy tensor (2.13) without reference to the Bianchi identity, i.e. directly from the explicit expression (2.12) and the Khronon equation $\nabla_{\mu}S^{\mu} = 0$.

The formulation of the theory in the adapted coordinate system for which $t = \tau$ (unitary gauge), is relagated to appendix A.

# 3 Nonrelativistic limit

## 3.1 Post Newtonian expansions for the fields

We investigate the nonrelativistic (or slow motion) limit of the model, in the case of an isolated matter system (e.g. the Milky Way galaxy), and determine under which conditions we can recover the MOND phenomenology in the weak acceleration regime. In this approximation, we introduce two a priori different scalar potentials $\phi$ and $\psi$, as well as the “gravitomagnetic” vector potential $\zeta$, and make the usual post-Newtonian (PN) ansatz on the metric components generated by the isolated system, using the ADM form (A.1),

$$N = 1 + \frac{\phi}{c^2} + \mathcal{O}(c^{-4}), \qquad (3.1a)$$

$$N_i = \frac{4}{c^3} \zeta_i + \mathcal{O}(c^{-5}), \qquad (3.1b)$$

$$q_{ij} = \gamma_{ij} \left( 1 - \frac{2\psi}{c^2} \right) + \mathcal{O}(c^{-4}), \qquad (3.1c)$$

where $\gamma_{ij}$ is a Euclidean metric in general coordinates, and with usual notation for the PN remainder terms $\mathcal{O}(c^{-n})$. The PN ansatz (3.1) is standard in GR and follows from the leading PN order of the stress-energy tensor $T^{\mu\nu}$ for ordinary matter fields, namely

$$T^{00} = \rho_m c^2 + \mathcal{O}(c^0), \qquad (3.2a)$$

$$T^{0i} = \rho_m v_m^i c + \mathcal{O}(c^{-1}), \qquad (3.2b)$$

$$T^{ij} = \rho_m v_m^i v_m^j + p_m \gamma^{ij} + \mathcal{O}(c^{-2}), \qquad (3.2c)$$

where $\rho_m$, $p_m$ and $v_m^i$ are the mass density, pressure and coordinate velocity of the matter field in the Newtonian approximation. Here $\rho_m$ is the coordinate density, to which $\rho$ [see eqs. (A.5)] reduces in the Newtonian limit, and obeys the ordinary continuity equation

$$\dot{\rho}_m + \nabla_i (\rho_m v_m^i) = 0. \qquad (3.2d)$$

The dot denotes the time derivative and $\nabla_i$ is the covariant derivative associated with $\gamma_{ij}$.

In order to justify the PN ansatz (3.1) for the present theory we must also take into account the Khronon field equation. In particular, we must ensure that the stress-energy tensor $T^{\mu\nu}$ of the Khronon field given by eq. (2.12), admits the same leading PN behaviour as for the ordinary matter in (3.2). As shown by Flanagan [47] this will be satisfied if we
assume that the expansion of the Khronon field is of the type

$$ \tau = t + \frac{\sigma(\mathbf{x}, t)}{c^2} + \mathcal{O}(c^{-4}), \qquad (3.3) $$

where $\sigma$ is the leading order perturbation. We shall check this point in eq. (3.14) below.

Inserting (3.3) into (2.1) we obtain the PN expansion of the scalar field $Q$ and the foliation's unit vector $n_\mu$ as

$$ Q = 1 - \frac{\Xi}{c^2} + \mathcal{O}(c^{-4}), \qquad (3.4a) $$

$$ n_0 = -1 - \frac{1}{c^2}(\Xi + \dot{\sigma}) + \mathcal{O}(c^{-4}), \qquad (3.4b) $$

$$ n_i = -\frac{1}{c}\left(1 + \frac{\Xi}{c^2}\right) \vec{\nabla}_i \sigma + \mathcal{O}(c^{-5}), \qquad (3.4c) $$

where we have defined (following the notation in [47])

$$ \Xi \equiv \phi - \dot{\sigma} + \frac{1}{2} |\vec{\nabla} \sigma|^2. \qquad (3.5) $$

We obtain in turn the acceleration vector $A_\mu$ defined by (2.3) as

$$ A_0 = \frac{1}{c} \vec{\nabla} \sigma \cdot \vec{\nabla} \Xi + \mathcal{O}(c^{-3}), \qquad (3.6a) $$

$$ A_i = \vec{\nabla}_i \Xi + \mathcal{O}(c^{-2}). \qquad (3.6b) $$

To fully determine the slow motion limit, we need to know the PN order of the functions $\mathcal{J}$ and $\mathcal{K}$ in the action. We find that

$$ \mathcal{Y} = \frac{1}{c^4} |\vec{\nabla} \Xi|^2 + \mathcal{O}(c^{-5}), \qquad (3.7) $$

and as we show in (3.23) below, in the deep MOND limit $\mathcal{J} \sim \Lambda - \mathcal{Y} + c^2 a_0^{-1} \mathcal{Y}^{3/2} + \mathcal{O}(\mathcal{Y}^2)$, hence, the PN order of this function is $\mathcal{J} = \mathcal{O}(c^{-4})$. This is consistent with the fact that the cosmological constant scales with the MOND acceleration constant $a_0$ like $\Lambda \sim a_0^2/c^4$ and is thus formally $\mathcal{O}(c^{-4})$.

The PN order of the function $\mathcal{K}$ is similarly determined by appealing to the definition (2.7) with the constant $\mu$ having a non-zero finite limit in the PN expansion, i.e. $\mu = \mathcal{O}(c^0)$. In particular, with eq. (3.4a) this implies the following leading PN orders

$$ \mathcal{K} = \frac{\mu^2 \Xi^2}{c^4} + \mathcal{O}(c^{-6}), \qquad (3.8a) $$

$$ \mathcal{K}_Q = -\frac{2\mu^2\Xi}{c^2} + \mathcal{O}(c^{-4}), \qquad (3.8b) $$

$$ \mathcal{Q}\mathcal{K}_Q = -\frac{2\mu^2\Xi}{c^2} + \mathcal{O}(c^{-4}). \qquad (3.8c) $$

The above PN orders hold for the choice we have made for the function $\mathcal{K}$ in the action, i.e. the $(Q-1)^2$ case. But obviously, if we include in $\mathcal{K}$ higher powers of $Q-1$ as we do in section 4.3 this will not change the leading PN limit (3.8).
## 3.2 PN expansion of the field equations

### 3.2.1 The *ij* components of the Einstein equation

Considering first the *ij* components of the field equation (2.9), using the facts that $\mathcal{J}$ and $\mathcal{K}$ are small quantities $\mathcal{O}(c^{-4})$ together with $n_i = \mathcal{O}(c^{-1})$, we obtain<fnref n="4" />

$$ \Delta (\phi - \psi) \gamma_{ij} - \vec{\nabla}_i \vec{\nabla}_j (\phi - \psi) = \mathcal{O}(c^{-2}), \quad (3.9) $$

where $\Delta \equiv \gamma^{ij}\vec{\nabla}_i\vec{\nabla}_j$ is the Laplace operator. The previous equation yields, in the case of regular isolated matter systems,

$$ \phi = \psi + \mathcal{O}(c^{-2}). \quad (3.10) $$

The equality of the two potentials $\phi$ and $\psi$ in the nonrelativistic limit is very important for the viability of the theory as an alternative to dark matter, as it implies that the light deflection and the gravitational lensing is given by the same formula as in GR. That is, for any nonrelativistic baryonic distribution of matter for which the forces $\vec{\nabla}\phi$ are in harmony with observations as if dark matter was present, the same matter distribution will also lead to the correct gravitational lensing signal, again as if dark matter was present. This fact is actually true for the general class of modified Einstein-Aether theories (in particular the AeST theory) and was first noticed in [13].

### 3.2.2 The 00 Einstein equation

Next we investigate the equation satisfied by the potential $\phi$. This follows from the 00 component of the Einstein field equations (2.9). In this case we find that the Khronon kinetic term $\mathcal{K}$ in the action contributes to the nonrelativistic limit through the term $\propto c^2\mathcal{Q}\mathcal{K}_Q$. Using eq. (3.8c) and $\psi = \phi$ we obtain

$$ \Delta\phi + \vec{\nabla} \cdot (\mathcal{J}_y \vec{\nabla}\Xi) + \mu^2\Xi = 4\pi G\rho_m + \mathcal{O}(c^{-2}). \quad (3.11) $$

### 3.2.3 The 0i Einstein equation

The equation for the gravitomagnetic potential $\zeta_i$ is found using the 0i Einstein equation (2.9), and to leading PN order it gives

$$ \Delta\zeta^i - \vec{\nabla}^i\vec{\nabla}_j\zeta^j - \vec{\nabla}^i\dot{\phi} - [\vec{\nabla} \cdot (\mathcal{J}_y \vec{\nabla}\Xi) + \mu^2\Xi]\vec{\nabla}^i\sigma = 4\pi G\rho_m v_m^i + \mathcal{O}(c^{-2}), \quad (3.12) $$

where once again we have used the equality of the potentials $\phi$ and $\psi$ through (3.10).

$$ G^{00} = \frac{2}{c^2} \Delta\psi + \mathcal{O}(c^{-4}), $$

$$ G^{0i} = \frac{2}{c^3} [\Delta\zeta^i - \vec{\nabla}^i\vec{\nabla}_j\zeta^j + \vec{\nabla}^i\dot{\psi}] + \mathcal{O}(c^{-5}), $$

$$ G^{ij} = \frac{1}{c^2} [\Delta(\phi - \psi)\gamma^{ij} - \vec{\nabla}^i\vec{\nabla}^j(\phi - \psi)] + \mathcal{O}(c^{-4}). $$

<footnote n="4">⁴We note for reference the leading PN behaviour of the components of the Einstein tensor:</footnote>
### 3.2.4 Consistency of the PN field equations

In analogy with the matter equations (3.2) and following [47], we define the Khronon “mass density”

$$ \rho_{\tau} = -\frac{1}{4\pi G} \left[ \vec{\nabla} \cdot (\mathcal{J}_{y} \vec{\nabla} \Xi) + \mu^2 \Xi \right], \qquad (3.13a) $$

and the Khronon “velocity” field

$$ v_{\tau}^{i} = -\vec{\nabla}^{i} \sigma. \qquad (3.13b) $$

Such definitions are motivated by the fact that the components of the Khronon stress-energy tensor (2.12) in leading PN order take the same form as for ordinary matter, i.e.

$$ \mathcal{T}^{00} = \rho_{\tau} c^2 + \mathcal{O}(c^0), \qquad (3.14a) $$

$$ \mathcal{T}^{0i} = \rho_{\tau} v_{\tau}^{i} c + \mathcal{O}(c^{-1}), \qquad (3.14b) $$

$$ \mathcal{T}^{ij} = \mathcal{O}(c^0). \qquad (3.14c) $$

The above PN scaling of the Khronon tensor justifies the PN expansion we have adopted for the metric (3.1) and Khronon field (3.3).

With the previous definitions, we may rewrite the Einstein field equations (3.11) and (3.12) in the useful form

$$ \Delta\phi = 4\pi G(\rho_m + \rho_\tau) + \mathcal{O}(c^{-2}), \qquad (3.15a) $$

$$ \Delta\zeta^i - \vec{\nabla}^i \vec{\nabla}_j \zeta^j - \vec{\nabla}^i \dot{\phi} = 4\pi G (\rho_m v_m^i + \rho_\tau v_\tau^i) + \mathcal{O}(c^{-2}), \qquad (3.15b) $$

which shows that the Einstein field equations, together with the continuity equation (3.2d) for matter, imply the continuity equation for the Khronon quantities (3.13), i.e.

$$ \dot{\rho}_{\tau} + \vec{\nabla}_{i} (\rho_{\tau} v_{\tau}^{i}) = \mathcal{O}(c^{-2}). \qquad (3.16) $$

The latter equation is nothing but the slow motion approximation of the Khronon field equation $\nabla_{\mu}S^{\mu} = 0$ since from eq. (2.14) we have

$$ S^0 = -4\pi G \rho_{\tau} + \mathcal{O}(c^{-2}), \qquad (3.17a) $$

$$ S^i = -4\pi G \frac{\rho_{\tau} v_{\tau}^i}{c} + \mathcal{O}(c^{-3}). \qquad (3.17b) $$

and the Kronon equation reduces to $\dot{S}^0 + c \partial_i S^i = \mathcal{O}(c^{-2})$. This shows the consistency of the Khronon equation with the Einstein field equations in the slow motion approximation.

### 3.2.5 Check in adapted coordinates

By performing a coordinate transformation $\{x^{\mu}\} \rightarrow \{\bar{x}^{\mu}\}$ such that $\bar{t} \equiv \bar{x}^0/c = t + \sigma/c^2$ and $\bar{x}^i = x^i$, one obtains coordinates adapted to the foliation — the so-called unitary gauge as used in appendix A. Indeed the Khronon field becomes

$$ \bar{\tau} = \bar{t} + \mathcal{O}(c^{-4}), \qquad (3.18) $$
while the ADM form of the metric in these coordinates reads

$$ \bar{N} = 1 + \frac{\Xi}{c^2} + \mathcal{O}(c^{-4}), \qquad (3.19a) $$

$$ \bar{N}_i = \frac{1}{c} \vec{\nabla}_i \sigma + \frac{1}{c^3} \left[ 4\zeta_i + 2\phi \vec{\nabla}_i \sigma - \sigma \vec{\nabla}_i \dot{\sigma} \right] + \mathcal{O}(c^{-5}), \qquad (3.19b) $$

$$ \bar{q}_{ij} = \gamma_{ij} \left( 1 - \frac{2\psi}{c^2} \right) - \frac{1}{c^2} \vec{\nabla}_i \sigma \vec{\nabla}_j \sigma + \mathcal{O}(c^{-4}), \qquad (3.19c) $$

where we recall that $\Xi = \phi - \dot{\sigma} + \frac{1}{2}|\vec{\nabla}\sigma|^2$. With adapted coordinates (to leading PN order) we can directly insert (3.19) into the variational equations for the lapse, shift and spatial metric in appendix A. In particular we find that the momentum constraint equation (A.7) is identically satisfied to leading order ($\propto 1/c$) while the Hamiltonian equation (A.6) and evolution equation (A.8) are equivalent to the equality of potentials (3.10) together with the equation (3.11). More work should also permit to confirm eq. (3.12) by working out the momentum constraint equation (A.7) to next-to-leading order $\propto 1/c^3$.

To summarize, the MOND dynamics is recovered in the unitary gauge as well, as it should be, since physics does not depend on the gauge choice. This was the approach taken in ref. [18]. However, as was realized in ref. [47] and we have verified in this section, using the unitary gauge is incompatible with the standard PN expansion in (3.1), and rather, the expansion of the shift vector acquires the term $\vec{\nabla}_i\sigma/c$ which dominates over the usual $\zeta_i/c^3$ assumed in eq. (3.1b) which is the standard case. Doing that, one does recover the full PN limit, including the MOND dynamics, in a consistent way, however, the resulting expansion is more complicated to compute because of the $\vec{\nabla}_i\sigma/c$ term in $\bar{N}_i$. Ref. [18] used the unitary gauge but without the $c^{-1}$ term in the $\bar{N}_i$ expansion, concluding that the MOND dynamics is recovered. This would be equivalent to setting $\sigma = 0$ throughout section 3.2, and this is consistent only in stationary systems. Outside stationarity, $\sigma$ must necessarily be included to avoid inconsistencies with solving the PN equations. This is seen through the Khronon equation (3.16), which when $\sigma = 0$ implies that $\dot{\phi} = 0$ and from (3.15a) this latter condition is consistent with the PN expansion only if $\dot{\rho}_m = 0$, that is, for stationary systems.<fnref n="5" />

### 3.2.6 Stationary configurations

We now consider the case where the system is, in addition, stationary. This means that $\dot{\phi} = \dot{\sigma} = 0$ for the fields and $\dot{\rho}_m = \dot{v}_m^i = 0$ for the matter, thus, the continuity equation for matter reduces to the constraint $\vec{\nabla}_i(\rho_m v_m^i) = 0$. Furthermore, the Khronon equation (3.16) reduces to $\vec{\nabla} \cdot (\rho_\tau \vec{\nabla}\sigma) = 0$ which we can solve by choosing $\sigma = 0$. This is the unitary gauge, where the Khronon field has no perturbation and it is given by the time coordinate exactly.

<footnote n="5">⁵In GR, one could start the expansion of $\bar{N}_i$ at order $c^{-1}$, but one then finds that the $c^{-1}$ term can be removed by a gauge transformation and so the standard PN expansion is justified there. This is also correct in extensions of GR, provided one includes perturbations for all additional fields. However, imposing gauge conditions on the additional fields (e.g. the unitary gauge), removes the freedom to gauge away the $c^{-1}$ term in $\bar{N}_i$. Put it differently, as is relevant to our case, the same gauge transformation which removes the $c^{-1}$ term from $\bar{N}_i$, necessarily introduces a perturbation to the Khronon.</footnote>
The choice of the unitary gauge then leads to $\Xi = \phi$ and the acceleration vector (3.6) becomes

$$A_0 = \mathcal{O}(c^{-3}), \qquad (3.20a)$$

$$A_i = \vec{\nabla}_i \phi + \mathcal{O}(c^{-2}). \qquad (3.20b)$$

Thus the acceleration of the congruence of unit vector field $n_\mu$ reduces to the physical Newtonian acceleration, and from (3.11) the equation for the Newtonian potential $\phi$ becomes the modified Poisson (or modified Helmholtz) equation

$$\vec{\nabla} \cdot [(1 + \mathcal{J}_Y) \vec{\nabla} \phi] + \mu^2 \phi = 4\pi G \rho_m + \mathcal{O}(c^{-2}). \qquad (3.21)$$

### 3.2.7 Choice of function $\mathcal{J}(\mathcal{Y})$

The equation (3.21) is still too general as we have not yet specified the function $\mathcal{J}(\mathcal{Y})$. Indeed, one can get a wide variety of solutions through such a choice, many of which may not necessarily lead to either Newtonian nor MONDian dynamics and may not fit observations altogether. Our goal is to specify the conditions that $\mathcal{J}$ needs to satisfy in order that both Newtonian and MONDian regimes emerge under appropriate conditions. The fact that we shall recover the MOND regime only for stationary systems is still satisfying because most tests of MOND (rotation curves of galaxies, Tully-Fisher relation, etc.) are performed in stationary situations.

From the first term in (3.21) we readily identify the MOND interpolating function $f(y)$, where $y \equiv |\nabla\phi|/a_0$, as$^6$

$$f(y) = 1 + \mathcal{J}_Y(\mathcal{Y}), \qquad (3.22)$$

and from (3.20) we can just take $\mathcal{Y} = A_\mu A^\mu / c^4$ to be $\mathcal{Y} \simeq a_0^2 y^2 / c^4$ in the nonrelativistic limit. Admissible functions $\mathcal{J}(\mathcal{Y})$ are those for which the MOND function is linear, $f(y) \simeq y$, for low accelerations $y \ll 1$, while it tends to a constant in the strong field regime when $y \gg 1$. The latter constant can be set to be 1 so that $G$ represents the measured Newton's gravitational constant.

Integrating (3.22) we obtain the required behaviour of the function $\mathcal{J}$ in order for MOND to emerge in the limit $y \ll 1$:

$$\mathcal{J} = \Lambda - \mathcal{Y} + \frac{2c^2}{3a_0} \mathcal{Y}^{3/2} + \mathcal{O}\left(y^4/a_0^2\right), \qquad (3.23)$$

where $\Lambda$ is the cosmological constant. Note that the low acceleration regime $y \ll 1$ is also the regime of linear cosmological perturbations, since the acceleration $A_\mu$ is a first order quantity in cosmological perturbations. Thus, $\mathcal{Y}$ is second order in cosmology and therefore, $\Lambda$ is the cosmological constant measured in cosmological observations, e.g. the CMB.

With the above condition on the function $\mathcal{J}$,$^7$ the equation (3.21) reduces to the usual MOND equation, i.e., having the Bekenstein-Milgrom form [5], apart from the presence of

<footnote>$^6$We denote the MOND function as $f(y)$ instead of the usual $\mu(y)$ to avoid conflict with the constant scale $\mu$ in eq. (3.21). We may also denote $\chi_s(y) = f(y) - 1 = \mathcal{J}_Y(\mathcal{Y})$ the “gravitational susceptibility” of the dark matter (in the sense of [57]).</footnote>

<footnote>$^7$Note that $\mathcal{J}$ is always of order $\mathcal{O}(c^{-4})$, and $\mathcal{J}_Y$ of order $\mathcal{O}(c^0)$, irrespective of the expansions (3.23) or (3.26). The constant $a_0$ serves as a second expansion parameter, which determines corrections to MOND or to Newton, and is independent of the PN parameter $c^{-1}$.</footnote>
the “mass” term $\mu^2\phi$. This term is absent from the BM theory, and is the result of the function $K(Q)$, in the same way as in AeST theory [39] (see also [58]). As such, the $\mu^2\phi$ term is connected to having viable cosmological solutions and provides a link between the stationary slow motion approximation and cosmology, as we discuss further in section 4.

We expect the phenomenology of spherically symmetric solutions to the Khronon theory to be very similar to that in the case of the AeST theory, studied in [58]. The effect of the $\mu^2\phi$ term in eq. (3.21) is then to induce an oscillatory behaviour for the potential $\phi$ on scales larger than the critical scale $r_C$ given by

$$r_C \sim \left(\frac{r_M}{\mu^2}\right)^{1/3} \quad \text{with} \quad r_M \sim \sqrt{\frac{GM}{a_0}}, \qquad (3.24)$$

and can be neglected when the distance $r$ with respect to the center of the source (a typical galaxy) is much smaller than $r_C$. Here $r_M$ denotes the MOND transition radius ($M$ being the mass of the source) at which the physical acceleration $|\nabla\phi|$ equals $a_0$. The result (3.24) follows from eq. (3.21) in the MOND regime where $1 + \mathcal{J}_y \simeq |\nabla\phi|/a_0$, by investigating the scale at which the $\mu^2\phi$ term will dominate over the first term (“$\mu$-domination” regime [58]). Since the goal is to describe both the Newtonian and MOND regimes of the galaxy we must ensure that $r_C \gg r_M$. With $a_0 \simeq 1.2 \ 10^{-10} \text{m/s}^2$ the MOND critical acceleration scale and $M \sim 10^{11} M_\odot$ for the mass of the galaxy, we have $r_M \sim 10 \text{ kpc}$ and we generously choose

$$\mu^{-1} \gg 100 \text{ kpc}. \qquad (3.25)$$

For instance we could take $\mu^{-1} \gtrsim 1 \text{ Mpc}$, or even larger.

Finally, in order to recover the Newtonian law in the high acceleration regime $y \gg 1$, we impose that $\mathcal{J}$ behaves as

$$\mathcal{J} = \Lambda_{\infty} + \mathcal{O}\left(\frac{a_0}{y}\right), \qquad (3.26)$$

where the constant $\Lambda_\infty$ is a priori different from $\Lambda$ in (3.23), and can be seen to be the effective cosmological constant in the strong-field (GR) regime. With (3.26) imposed the model will in fact recover the full GR in this regime. In particular we conclude that the theory has the same parametrized post-Newtonian (PPN) limit as GR and is therefore viable in the Solar System and for binary pulsar tests; see eq. (6.6) where $\beta = \lambda = 0$, and with $\alpha = 0$ in the high acceleration regime. Notice that for the choice we have made previously for the constant $\mu^{-1}$, say $\mu^{-1} \gtrsim 1 \text{ Mpc}$, the “mass” term $\mu^2\phi$ is negligible at the scale of galaxies, and is a fortiori negligible at the much smaller scale of the Solar System, and has no incidence on the PPN parameters.

# 4 Cosmology

## 4.1 Friedmann-Lemaître-Robertson-Walker (FLRW) cosmology

### 4.1.1 General FLRW background

We consider first a FLRW background cosmology in the synchronous coordinate system. For simplicity, we choose $c = 1$ from this section onwards. The metric is

$$ds^2 = -dt^2 + a^2 \gamma_{ij} dx^i dx^j, \qquad (4.1)$$
where $t$ is the cosmic time and $\gamma_{ij}^{\kappa}$ is a spatially homogeneous and isotropic metric of spatial curvature $\kappa$. The symmetries of the FLRW spacetime also require that the Khronon field is a function of time only, that is, $\tau = \bar{\tau}(t)$, so that the background value of $Q$ is

$$ \bar{Q} = \dot{\bar{\tau}}, \qquad (4.2) $$

with the dot denoting derivatives with respect to cosmic time $t$, and the overbar always denoting the background FLRW value. Furthermore, the acceleration vanishes in the background, $\bar{A}_{\mu} = 0$.

Defining the Hubble parameter as $H \equiv \dot{a}/a$, the Einstein equations (2.9) lead to the Friedmann equations

$$ 3H^2 + \frac{3\kappa}{a^2} - \Lambda = 8\pi G \sum_{I \neq \mathcal{K}} \bar{\rho}_I + \bar{Q}\bar{\mathcal{K}}_Q - \bar{\mathcal{K}}, \qquad (4.3a) $$

$$ -\left(2\dot{H} + 3H^2 + \frac{\kappa}{a^2}\right) + \Lambda = 8\pi G \sum_{I \neq \mathcal{K}} \bar{P}_I + \bar{\mathcal{K}}, \qquad (4.3b) $$

where $\bar{\mathcal{K}} \equiv \mathcal{K}(\bar{Q})$ and $\bar{\mathcal{K}}_Q \equiv \mathcal{K}_Q(\bar{Q})$, and as indicated the index $I \neq \mathcal{K}$ runs over all types of ordinary matter species, i.e., excluding the contribution of the Khronon field. The cosmological constant $\Lambda$ is the one appearing in eq. (3.23). The Khronon contribution makes it possible to defining the Khronon energy density and pressure respectively as

$$ \bar{\rho}_{\mathcal{K}} = \frac{1}{8\pi G} (\bar{Q}\bar{\mathcal{K}}_Q - \bar{\mathcal{K}}), \qquad (4.4a) $$

$$ \bar{P}_{\mathcal{K}} = \frac{\bar{\mathcal{K}}}{8\pi G}. \qquad (4.4b) $$

Thus, on a FLRW background, the Khronon can be casted as a perfect fluid with time-dependent equation of state $w(t)$ given by

$$ w \equiv \frac{\bar{P}_{\mathcal{K}}}{\bar{\rho}_{\mathcal{K}}} = \frac{\bar{\mathcal{K}}}{\bar{Q}\bar{\mathcal{K}}_Q - \bar{\mathcal{K}}}. \qquad (4.5) $$

As the equation of state $w$ is time-dependent, this also prompts the definition of the adiabatic speed of sound $c_{\text{ad}}^2$ defined as

$$ c_{\text{ad}}^2 \equiv \frac{d\bar{P}_{\mathcal{K}}}{d\bar{\rho}_{\mathcal{K}}} = \frac{\bar{\mathcal{K}}_Q}{\bar{Q}\bar{\mathcal{K}}_Q\bar{Q}}, \qquad (4.6) $$

which obeys the well-known equation $\dot{w} = 3(w - c_{\text{ad}}^2)(1 + w)H$, see e.g. [52].

The above casting of the Khronon field into a perfect fluid, is consistent with the Bianchi identities and energy conservation, $\dot{\bar{\rho}}_{\mathcal{K}} + 3H(1+w)\bar{\rho}_{\mathcal{K}} = 0$, which can also be derived from the Khronon equation $\nabla_{\mu}\mathcal{S}^{\mu} = 0$, where (2.14) implies $\bar{\mathcal{S}}^0 = -c^2\bar{\mathcal{K}}_Q$ and $\bar{\mathcal{S}}^i = 0$. Specifically, these equations may be integrated once to give

$$ \bar{\mathcal{K}}_Q = \frac{I_0}{a^3}, \qquad (4.7) $$

where $I_0$ is a constant set by initial conditions. Note that the non-integrated version of the Khronon equation with (2.14), leads in the FLRW background to

$$ \dot{\bar{Q}} + 3Hc_{\text{ad}}^2\bar{Q} = 0, \qquad (4.8) $$

a relation which can be sometimes useful in reducing formulas.
### 4.1.2 FLRW background in adapted coordinates

Rather than the synchronous coordinate system one can also use the adapted coordinates of appendix A where $t = \tau$. In this case, one cannot set the lapse to 1 as in the synchronous case, and the metric must take the form

$$ds^2 = -\bar{N}^2 d\tau^2 + a^2 \gamma_{ij} dx^i dx^j, \qquad (4.9)$$

where $\bar{N} = \bar{N}(\tau) = 1/\bar{Q}(\tau)$ is at this point unspecified. This leads to the Friedmann equations

$$3 \frac{\mathcal{H}^2}{\bar{N}^2} + \frac{3\kappa}{a^2} - \Lambda + \bar{K} + \bar{N} \bar{K}_N = 8\pi G \sum_{I \neq K} \bar{\rho}_I, \qquad (4.10a)$$

$$-\frac{1}{\bar{N}^2} [2\mathcal{H}' + 3\mathcal{H}^2(1 - 2c_{\text{ad}}^2)] - \frac{\kappa}{a^2} + \Lambda - \bar{K} = 8\pi G \sum_{I \neq K} \bar{P}_I, \qquad (4.10b)$$

where we pose $\mathcal{H} = \frac{a'}{a} = \bar{N}H$, and the dependence of $\bar{K}$ on $\bar{N}$ is via the condition $\bar{Q} = 1/\bar{N}$. The Khronon equation (2.14) in the background is equivalent to $\bar{N}' = 3\mathcal{H}c_{\text{ad}}^2\bar{N}$, and yields

$$\bar{N}^2 \bar{K}_N = -\frac{I_0}{a^3}, \qquad (4.11)$$

which becomes an algebraic relation determining $\bar{N}(a)$. This in turn, is used in (4.10) turning it into a modified Friedmann equation and removing also the fluid interpretation of the Khronon (recall that $I$ runs over all types of matter species but excludes the Khronon).

### 4.1.3 Approximate dust solutions

Specifying a function $\mathcal{K}(Q)$, completely determines the form of $\bar{\rho}_K$, $\bar{P}_K$, $w$ and $c_{\text{ad}}^2$. In addition, eq. (4.7) may then be inverted to find a solution for $\bar{Q}(a)$, and equivalently (4.11) may be inverted to find $\bar{N}(a)$. This then, when inserted into (4.4), determines the exact dependence of the Khronon energy density and pressure in terms of the scale factor $a$. Thus the form of $\mathcal{K}$ fully determines the type of cosmological evolution that ensues, and, indeed, depending on this form, one can have cosmologies which are not necessarily compatible with observations in the absence of a dark matter component. Examples are, the stiff-fluid case where $\mathcal{K} \propto Q^2$ ($w = 1$) and the radiation-type case $\mathcal{K} \propto Q^4$ ($w = 1/3$).

Our interest is to determine the conditions on $\mathcal{K}(Q)$ that may lead to approximate dust solutions for the Khronon, in accordance with observations. Let us first say that exact dust solutions, i.e. $w = 0$, are impossible as they would imply that $\mathcal{K} = 0$. A simple case emerges following [48, 49]. In the so-called shift-symmetric $k$-essence model, when the action for a scalar field $\tau$ takes the form $(X+1)^2$ with $X = (\nabla\tau)^2$ the usual canonical kinetic term, there are approximate dust solutions in cosmology. Since $X = -Q^2$ for the Khronon, it is clear that approximate dust solutions exist if $\mathcal{K}(Q)$ is expandable as a Taylor series around $Q = 1$. This was the motivation for the specific choice made in eq. (2.7).

We first focus on this specific case, that we refer to as the “quadratic” function, which we recall here:

$$\mathcal{K}(Q) = \mu^2(Q - 1)^2. \qquad (4.12)$$
Then $K_Q = 2\mu^2(Q-1)$ and eq. (4.7) can be solved for $\overline{Q}$ to get

$$ \overline{Q} = 1 + \frac{I_0}{2\mu^2 a^3}, \qquad (4.13) $$

where $I_0$ is a constant. Hence we find that the energy density of the Khronon is

$$ \bar{\rho}_{\mathcal{K}} = \frac{\mu^2}{8\pi G} (\bar{Q}^2 - 1) = \frac{I_0}{8\pi G a^3} \left(1 + \frac{\tilde{w}_0}{a^3}\right), \qquad (4.14) $$

where we posed $\tilde{w}_0 = \frac{I_0}{4\mu^2}$. The constant $I_0$ is directly related to the Khronon energy density today $\bar{\rho}_{\mathcal{K},0} = \frac{I_0}{8\pi G}(1 + \tilde{w}_0)$ (i.e. when $a=1$); in particular it is positive meaning also that $\overline{Q} \ge 1$. We can rewrite

$$ \bar{\rho}_{\mathcal{K}} = \frac{\bar{\rho}_{\mathcal{K},0}}{a^3} \frac{1 + \frac{\tilde{w}_0}{a^3}}{1 + \tilde{w}_0} = \frac{\bar{\rho}_{\mathcal{K},0}}{a^3} \left(1 - w_0 + \frac{w_0}{a^3}\right), \qquad (4.15) $$

where we have also introduced the value of the Khronon equation of state today $w_0 = \frac{\tilde{w}_0}{1+\tilde{w}_0}$, see (4.17). If the Khronon is to give full account of the total dark matter in this model the Khronon energy density today should be

$$ \bar{\rho}_{\mathcal{K},0} = \frac{3H_0^2}{8\pi G} \Omega_{\mathcal{K},0}, \qquad (4.16) $$

where $H_0$ is the Hubble-Lemaître parameter and $\Omega_{\mathcal{K},0} \sim 0.26$ is the measured fraction of dark matter. From (4.15), the Khronon energy density is approximately that of dust, i.e. up to some additional correction scaling as $a^{-6}$, while the equation of state evolves as

$$ w = \frac{\overline{Q} - 1}{\overline{Q} + 1} = \frac{\tilde{w}_0}{\tilde{w}_0 + a^3}. \qquad (4.17) $$

In the deep past, $a \to 0$ and $w \to 1$ so that the Khronon with the quadratic functional form of $K$ behaves as a stiff fluid. The adiabatic speed of sound is likewise given as

$$ c_{\text{ad}}^2 = 1 - \frac{1}{\overline{Q}} = \frac{2\tilde{w}_0}{2\tilde{w}_0 + a^3}. \qquad (4.18) $$

The equation of state $w$ for the quadratic potential will be shown in figure 1 together with the results from other potentials investigated in section 4.3.

## 4.2 Linear perturbations on a FLRW background

We now consider linear perturbations on a FLRW background and show that one can recast the Khronon equations into the generalized dark matter (GDM) model [51, 52]. This is a model defined only on the FLRW background and linearized perturbation level (and thus without a unique non-linear completion), and determined by three parametric functions: the equation of state $w(t)$, the sound speed $c_s^2(t, k)$ and viscosity $c_{\text{vis}}^2(t, k)$; the last two parameters appear only at the linearized perturbed level and may depend on the wavenumber $k$.

Keeping the synchronous coordinate system for the background FLRW dynamics, and considering only scalar modes, the metric is

$$ ds^2 = -(1+2\phi)dt^2 - 2\nabla_i\zeta dt dx^i + a^2[(1-2\psi)\gamma_{ij}^\kappa + \hat{D}_{ij}\nu]dx^i dx^j, \qquad (4.19) $$
where $\hat{D}_{ij} \equiv \vec{\nabla}_i \vec{\nabla}_j - \frac{1}{3}\gamma_{ij}^\kappa \vec{\nabla}^2$ is a traceless covariant derivative operator ($\vec{\nabla}_i \gamma_{jk}^\kappa \equiv 0$). We also perturb the Khronon field to first order as

$$ \tau = \bar{\tau} + \sigma, \qquad (4.20) $$

leading to $\bar{Q} = \dot{\bar{\tau}}$ as in the previous section, and $Q = \bar{Q}(1 - \phi) + \dot{\sigma}$, while the normal $n_\mu$ to the spacetime foliation and acceleration $A_\mu$ components are found as

$$ n_0 = -1 - \phi, \qquad n_i = -\frac{1}{\bar{Q}} \vec{\nabla}_i \sigma, \qquad (4.21a) $$

$$ A_0 = 0, \qquad A_i = \vec{\nabla}_i \Upsilon, \qquad (4.21b) $$

respectively, where we have posed

$$ \Upsilon \equiv \phi - \partial_t \left( \frac{\sigma}{\bar{Q}} \right) = \phi - \frac{1}{\bar{Q}} \left( \dot{\sigma} + 3Hc_{\text{ad}}^2 \sigma \right). \qquad (4.22) $$

From the perturbation variables for both the metric and Khronon one can define in a usual way three independent gauge invariant variables as

$$ \Psi \equiv \psi + \frac{1}{6} \vec{\nabla}^2 \nu + H \left( \zeta + \frac{a^2}{2} \dot{\nu} \right), \qquad (4.23a) $$

$$ \Phi \equiv \phi - \partial_t \left( \zeta + \frac{a^2}{2} \dot{\nu} \right), \qquad (4.23b) $$

together with $\Upsilon$ for the Khronon, which is directly gauge invariant since the acceleration (4.21b) vanishes in the background.

With these, we determine the perturbed Khronon equation (2.14) as

$$ \vec{\nabla}^2 (\dot{\Upsilon} + H\Upsilon) = 4\pi G a^2 (1+w) \bar{\rho}_K \left\{ \partial_t \left[ \frac{1}{c_{\text{ad}}^2} \left( \phi - \frac{\dot{\sigma}}{\bar{Q}} \right) + 3\psi \right] - \frac{1}{a^2} \vec{\nabla}^2 \left( \zeta - \frac{\sigma}{\bar{Q}} \right) \right\}. \qquad (4.24) $$

This equation is gauge invariant, as it can be entirely expressed in terms of the gauge invariant variables $\Psi$, $\Phi$ and $\Upsilon$. In the case of an ordinary matter species “I”, we find the stress-energy tensor components as

$$ T_{I0}^0 = -\bar{\rho}_I (1 + \delta_I), \qquad (4.25a) $$

$$ T_{Ii}^0 = -(\bar{\rho}_I + \bar{P}_I) \vec{\nabla}_i \theta_I, \qquad (4.25b) $$

$$ T_{I0}^i = -\frac{1}{a^2} (\bar{\rho}_I + \bar{P}_I) \vec{\nabla}_i (\zeta - \theta_I), \qquad (4.25c) $$

$$ T_{Ij}^i = \bar{\rho}_I (w_I + \Pi_I) \delta_{ij} + (\bar{\rho}_I + \bar{P}_I) \hat{D}_j^i \Sigma_I, \qquad (4.25d) $$

where the density contrast is defined as $\delta_I \equiv \delta \rho_I / \bar{\rho}_I$ and where $\theta_I$ is the velocity divergence, $\Pi_I$ the pressure contrast and $\Sigma_I$ the anisotropic stress. Applying these definitions to the
stress-energy tensor (2.12) of the Khronon field,<fnref n="8" /> we find

$$ \delta_{\kappa} = \frac{1+w}{c_{\text{ad}}^2} \left( \frac{\dot{\sigma}}{\dot{Q}} - \phi \right) + \frac{\vec{\nabla}^2 \Upsilon}{4\pi G a^2 \bar{\rho}_{\kappa}}, \quad (4.26a) $$

$$ \theta_{\kappa} = \frac{\sigma}{\dot{Q}}, \quad (4.26b) $$

$$ \Pi_{\kappa} = (1+w) \left( \frac{\dot{\sigma}}{\dot{Q}} - \phi \right), \quad (4.26c) $$

$$ \Sigma_{\kappa} = 0. \quad (4.26d) $$

Combining the above results for $\delta_{\kappa}$ and $\theta_{\kappa}$ with eq. (4.22) we note that the perturbation variable $\Upsilon$ for the Khronon can be obtained by solving the equation

$$ \vec{\nabla}^2 \Upsilon = 4\pi G \bar{\rho}_{\kappa} a^2 \left\{ \frac{1+w}{c_{\text{ad}}^2} \Upsilon + \Delta_{\kappa} \right\}, \quad (4.27) $$

where we have introduced the (necessarily gauge-invariant) co-moving density contrast of the Khronon as

$$ \Delta_{\kappa} \equiv \delta_{\kappa} + 3H(1+w)\theta_{\kappa}. \quad (4.28) $$

From eqs. (4.26), the Khronon fluid is no longer an adiabatic fluid in first order perturbation, and we have $\Pi_{\kappa} = c_{\text{ad}}^2 \delta_{\kappa} + \Pi_{\text{nad}}$ where the non-adiabatic pressure is

$$ \Pi_{\text{nad}} = - \frac{c_{\text{ad}}^2}{4\pi G a^2 \bar{\rho}_{\kappa}} \vec{\nabla}^2 \Upsilon. \quad (4.29) $$

Furthermore, after further calculation, one can show that

$$ \Pi_{\text{nad}} = (c_s^2 - c_{\text{ad}}^2) \Delta_{\kappa}, \quad (4.30) $$

where $c_s^2(t, \vec{x})$ is the speed of sound (generally a function of both time and space) and $\Delta_{\kappa}$ is the co-moving density contrast (4.28). The relation (4.30) is non-trivial and does not in general hold for any stress-energy tensor of the type (4.25); when it does (as in our case), the stress-energy tensor takes the form of the GDM model. Specifically, the speed of sound is found to be (in Fourier space)<fnref n="9" />

$$ c_s^2(t, k) = c_{\text{ad}}^2 \left[ 1 + \frac{c_{\text{ad}}^2 k^2}{4\pi G a^2 \bar{\rho}_{\kappa} (1+w)} \right]^{-1}, \quad (4.31) $$

and so as $k \to 0$, we have $c_s^2 \to c_{\text{ad}}^2$, that is, the Khronon is adiabatic in the limit $k \to 0$.

Turning to the Einstein equations (2.9), perturbed to linear order about the FLRW background, it is convenient to define the variable $\eta \equiv \psi + \frac{1}{6}\vec{\nabla}^2\nu$. With the previous

<footnote n="8">⁸Although we can call it $T_{\kappa}^{\mu\nu} \equiv T^{\mu\nu}$ to be consistent with (4.25)-(4.26), at linear order in perturbation the Khronon stress-energy tensor also includes the contribution from the function $\mathcal{J}$ through the term $\propto \nabla_{\rho} (\mathcal{J}_{\nu} A^{\rho})$ in (2.12), where we can use $\mathcal{J}_{\kappa} = -1$ with this approximation, see eq. (3.23).</footnote>

<footnote n="9">⁹From eq. (4.27) we have $\Upsilon = -\frac{c_s^2}{1+w}\Delta_{\kappa}$ in Fourier space.</footnote>
definitions, they are given, see e.g. eqs. (2.13) in [52], by the two constraint equations

$$-3H(\dot{\psi} + H\phi) + \frac{1}{a^2}[(\vec{\nabla}^2 + 3\kappa)\eta + H\vec{\nabla}^2\zeta] = 4\pi G \sum_I \bar{\rho}_I \delta_I, \quad (4.32a)$$

$$\dot{\eta} + H\phi + \frac{\kappa}{2}\left(\dot{\nu} + \frac{2}{a^2}\zeta\right) = 4\pi G \sum_I (\bar{\rho}_I + \bar{P}_I)\theta_I, \quad (4.32b)$$

that are found through the components $G_0^0$ and $G_i^0$ of the Einstein tensor, and the two propagation equations found by splitting $G^i_j$ into a trace part

$$\ddot{\psi} + 3H\dot{\psi} + H\dot{\phi} + [2\dot{H} + 3H^2]\phi + \frac{1}{3a^2}\vec{\nabla}^2[\phi - \eta - \dot{\zeta} - H\zeta] - \frac{\kappa}{a^2}\eta = 4\pi G \sum_I \bar{\rho}_I \Pi_I, \quad (4.33a)$$

and the traceless part

$$\ddot{\nu} + 3H\dot{\nu} + \frac{2}{a^2}(\eta - \phi + \dot{\zeta} + H\zeta) = 16\pi G \sum_I (\bar{\rho}_I + \bar{P}_I)\Sigma_I. \quad (4.33b)$$

In eqs. (4.32)–(4.33) the sums over index “I” include the Khronon field $I = K$ through the definitions (4.26) together with the background values in section 4.1.3. Furthermore, switching to the fluid variables for the Khronon, turns (4.24) into the fluid continuity equation for the density contrast

$$\dot{\delta}_K = 3H(w - c_s^2)\delta_K - 9H^2(1+w)(c_s^2 - c_{\text{ad}}^2)\theta_K - \frac{1+w}{a^2}\vec{\nabla}^2(\zeta - \theta_K) + 3(1+w)\dot{\psi}, \quad (4.34)$$

while a fluid Euler equation is found by combining the expression (4.26a) with the time derivative of (4.26b) as

$$\dot{\theta}_K = c_s^2 \left( 3H\theta_K + \frac{\delta_K}{1+w} \right) + \phi = 3Hc_{\text{ad}}^2\theta_K + \frac{\Pi_K}{1+w} + \phi. \quad (4.35)$$

The above equations for $\delta_K$ and $\theta_K$ are the standard shearless fluid equations. Thus we conclude, up to linear order around FLRW, that the Khronon behaves as a GDM fluid [51, 52] with zero viscosity $c_{\text{vis}}^2 = 0$, non-zero sound speed $c_s^2$ in perturbations given by (4.31), and time-dependent equation of state $w = w(a)$. We note, however, that once second order or higher perturbations are included, the correspondence with a perfect fluid will be lost. Moreover, the Khronon field does not contain a vector mode perturbation (being a scalar), thus, it does not lead to a purely vectorial velocity component even within the linearized fluid approximation, in sharp contrast with standard CDM.

For any function $K(Q)$ in the action for which $w \to 0$ and $c_{\text{ad}}^2 \to 0$, in the sense that they are arbitrarily small throughout the history of the Universe for $a \gtrsim 10^{-5}$ (i.e. for which the background solution to the Khronon energy density is approximately that of dust), the linear perturbations of the Khronon also behave as linear perturbations of dust. Thus, the approximate dust solutions discussed in 4.1.3 can be extended to the linearized regime, so that the Khronon model is expected to fit large scale structure and CMB data equally well as the $\Lambda$-CDM concordance model. Specifically, we can conservatively use constraints on GDM taken from [59, 60] and apply them to the Khronon model, however, due to the $k$-dependent sound speed which doesn't fall under the models studied in [59, 60], a proper comparison of the Khronon model with data is left for a future investigation.
## 4.3 Khronon cosmology and the MOND limit

### 4.3.1 Tension between the cosmology and MOND in the case of the quadratic function

As discussed above, with any choice for the function $K(Q)$ the Khronon field behaves like GDM on a FLRW background plus linear perturbations. Furthermore, the functional choice $K = \mu^2(Q - 1)^2$, the “quadratic” function, leads to approximate dust solutions in the late universe which can make the Khronon in agreement with cosmological observations, notably the CMB and large scale structure, provided these dust solutions can be extended as far back as the radiation-matter equality. However, the equation of state (4.17) for the quadratic function is not exactly zero but tends to $w \to 1$ in the Early Universe. The turning point may be found from (4.17) and is around $\tilde{w}_0 \sim a^3$.<fnref n="10" /> The constraints on GDM found in [59, 60] require that $w \lesssim 0.0164$ around $a \sim 10^{-4.5}$ at the 99% level (see the yellow shaded region in figure 1) with smaller values required until $a \sim 10^{-2}$ after which the constraints become weaker. It is thus sufficient to set

$$ \tilde{w}_0 \approx w_0 < 5.3 \times 10^{-16} \quad (\text{for cosmology}), \quad (4.36) $$

which represents the equation of state of the Khronon today. Now, from our discussion below (4.13), we have $I_0 \approx 8\pi G \bar{\rho}_{K,0}$ and $\tilde{w}_0 = I_0/(4\mu^2)$ so that combining the two (and restoring the relevant factors of c) we find

$$ \tilde{w}_0 = \frac{2\pi G \bar{\rho}_{K,0}}{\mu^2 c^2} = \frac{3H_0^2 \Omega_{K,0}}{4\mu^2 c^2}. \quad (4.37) $$

Given that from observations $H_0 \sim 70 \text{ km/s/Mpc} \sim 2.33 \times 10^{-4} \text{ c/Mpc}$ and $\Omega_{K,0} \sim 0.26$ we find that $\tilde{w}_0 \sim 1.06 \times 10^{-8}(\mu^{-1}/\text{Mpc})^2$ and thus cosmology places the bound

$$ \mu^{-1} \lesssim 0.22 \text{ kpc}. \quad (4.38) $$

Such a bound is incompatible with having a MOND limit in galaxies, as even for our own galaxy, we should have MOND behaviour out to tens of kpc. Indeed, reversing the problem and requiring

$$ \mu^{-1} \gtrsim 100 \text{ kpc}, \quad (4.39) $$

a very optimistic bound coming from matching the MOND phenomenology,<fnref n="11" /> leads to the constraint

$$ \tilde{w}_0 \gtrsim 1.06 \times 10^{-10} \quad (\text{for galaxies}). \quad (4.40) $$

Clearly then, these two inequalities on $\tilde{w}_0$ given by (4.36) and (4.40) are in conflict. Hence, the exact functional dependence for $K$ chosen in (4.12) cannot be in simultaneous harmony with observations of galaxies and with cosmology. Below, we investigate how different choices of the function $K(Q)$ may save the situation.

<footnote n="10">¹⁰From (4.17) the inflection point of the curve $w(a)$ occurs at $\tilde{w}_0 = 2a^3$.</footnote>

<footnote n="11">¹¹The requirement may in fact be even larger than (4.39) when matching to real data [see eq. (3.25)].</footnote>
### 4.3.2 Case of functions with higher (cubic and quartic) powers

From the discussion in [48], what we generally need in order to have dust solutions in the late Universe is that $\mathcal{K}$ be expandable as a Taylor series around $Q = 1$, that is, admissible functions should be expandable as

$$ \mathcal{K} = 2\mu^2 \sum_{n=1}^{+\infty} \frac{\mathcal{K}_{n+1}}{n+1} (Q-1)^{n+1}, \qquad (4.41) $$

where $\mathcal{K}_{n+1}$ are dimensionless constants and $\mathcal{K}_2 = 1$ to match our convention in (4.12).

To study the cosmological behaviour for a single power indexed by $n$, let us consider the functional form $\mathcal{K} = \frac{2\mu^2}{n+1}\mathcal{K}_{n+1}(Q-1)^{n+1}$ for integer $n \ge 2$. Then from (4.7) we find that

$$ \bar{Q} = 1 + \left( \frac{I_0}{2\mu^2 \mathcal{K}_{n+1} a^3} \right)^{1/n}, \qquad (4.42) $$

so that the energy density and pressure are found from (4.4) to be

$$ \bar{\rho}_{\mathcal{K}} = \frac{I_0}{8\pi G a^3} \left[ 1 + \frac{n}{n+1} \left( \frac{I_0}{2\mu^2 \mathcal{K}_{n+1} a^3} \right)^{1/n} \right], \qquad (4.43) $$

$$ \bar{P}_{\mathcal{K}} = \frac{I_0}{8\pi G(n+1)a^3} \left( \frac{I_0}{2\mu^2 \mathcal{K}_{n+1} a^3} \right)^{1/n}, \qquad (4.44) $$

and the equation of state is

$$ w = \left( \frac{I_0}{2\mu^2 \mathcal{K}_{n+1} a^3} \right)^{1/n} \left[ n + 1 + n \left( \frac{I_0}{2\mu^2 \mathcal{K}_{n+1} a^3} \right)^{1/n} \right]^{-1}. \qquad (4.45) $$

Clearly then, in the early universe, as $a \to \infty$, the Khronon once again behaves approximately as dust, however, the next-to-leading order behaviour is different. Importantly, for large $a$, the equation of state $w$ scales as $a^{-3/n}$ and, hence, for $n > 1$, this powerlaw is shallower than for the quadratic case, allowing us a better chance to reconcile the phenomenology of galaxies (MOND) with cosmology. Moreover, in the PN limit the relevant quantity to check (see section 3.2) is

$$ c^2 \mathcal{K}_Q = \frac{2\mu^2 \mathcal{K}_{n+1} (-\Xi)^n}{c^{2n-2}}, \qquad (4.46) $$

and specifically for $n = 2$ we get that $c^2 \mathcal{K}_Q \sim \Xi^2/c^2$ — which represents a 1PN correction negligible at the level of eq. (3.11). Thus for any $n \ge 2$ the function $\mathcal{K}$ does not contribute to the lowest PN limit and so we have more freedom to choose the coefficients $\mathcal{K}_{n+1}$ in order to pass the cosmological constraints (4.36). In the extreme case where we disallow the $n = 1$ term altogether, no conflict occurs between the MOND limit in galaxies and cosmology.

Realistically, though, it is unlikely that $n = 1$ can be ignored as, unless fine-tuning is present, typical functional forms will always have the $n = 1$ term (i.e., the quadratic function) in the expansion (4.41). In that case, we expect the Universe to go through different phases where $w$ scales initially as $a^{-3}$ (the $n = 1$ case) and then starts to gradually turn to different scalings until it settles to a final value in the early Universe. With this, we can expect to
<figure><img src="image_4.png" /></figure>

**Figure 1.** The Khronon equation of state parameter $w(a)$ versus scale factor $a$ for several models. Displayed are two quadratic models $K = \mu^2(Q-1)^2$ (solid red curves), one with $\mu^{-1} \sim 223$ kpc which is consistent with a MOND limit but which cannot fit the CMB, and the other with $\mu^{-1} \sim 1$ kpc which can fit the CMB but is inconsistent with a MOND limit. Introducing a cubic term with $K_3 \sim 10^5$ (green dot) or a quartic term with $K_4 \sim 10^6$ (black dash-dot) makes $\mu^{-1} \sim 223$ kpc compatible with the CMB. The two DBI-inspired models (4.47) with $\mu^{-1} \sim 223$ kpc and $\lambda_D \sim 30$, and $\mu^{-1} \sim 22.3$ Mpc and $\lambda_D \sim 1$ (blue dashed) are compatible with both a MOND limit and the CMB. The yellow shaded region is allowed by the data: it corresponds to the 99% credible regions of the binned “var-w” model of [59, 60], using only Planck Power Spectrum data, and with the equation of state piecewise parametrized into eight redshift bins.

reconcile galaxies with cosmology, and two specific examples are shown in figure 1: one case combining $n=1$ and $n=2$ with $K_3 \sim 10^5$, and one case combining $n=1$ and $n=3$ with $K_4 \sim 10^6$ (both cases assuming $\mu^{-1} \sim 223$ kpc). Unfortunately, we find that the required numerical values for $K_3$ or $K_4$ are unnaturally large. More generally, it seems that if the expansion (4.41) terminates at a finite order, large numbers for $K_{n+1}$ are inevitable, which warrants our next investigation below.

### 4.3.3 Dirac-Born-Infeld (DBI) type of functions

We would like to explore a function which has a well-defined Taylor expansion at $Q=1$ in order to reproduce the dust solutions we have explored in the previous sub-sections, see eq. (4.41), but which tends to a different behaviour when $Q \gg 1$. A particular case is the
Dirac-Born-Infeld (DBI) inspired function [61–63] of the form

$$K(Q) = \frac{2\mu^2}{\lambda_D} \left[ 1 - \sqrt{1 - \lambda_D(Q-1)^2} \right], \qquad (4.47)$$

where $\mu$ is the usual parameter as above and $\lambda_D$ is a new dimensionless parameter. At $Q=1$ we have $K=0$ and more specifically, for small $Q-1$ it admits the expansion $K \sim \mu^2(Q-1)^2 + \cdots$ which matches the quadratic function (4.12) and thus reproduces the required dust phenomenology in the late Universe; furthermore, as we have seen in (4.46), the PN limit is still safe. More generally, using (4.7), we find that

$$\bar{Q} = 1 + \frac{1}{\sqrt{\lambda_D + \left(\frac{2\mu^2 a^3}{I_0}\right)^2}}. \qquad (4.48)$$

Hence, as $a \to \infty$ we have $\bar{Q} \to 1$ (late Universe) while as $a \to 0$ then $\bar{Q} \to 1 + \lambda_D^{-1/2}$ (early Universe). Thus $\bar{Q}$ is bounded through $1 \le \bar{Q} < 1 + \lambda_D^{-1/2}$.

The energy density and pressure are found from (4.4) to be

$$\bar{\rho}_{\kappa} = \frac{1}{8\pi G} \left( \frac{I_0}{a^3} + \frac{2\mu^2}{\lambda_D} \left[ -1 + \sqrt{1 + \lambda_D \left( \frac{I_0}{2\mu^2 a^3} \right)^2} \right] \right), \qquad (4.49a)$$

$$\bar{P}_{\kappa} = \frac{1}{8\pi G} \frac{2\mu^2}{\lambda_D} \left[ 1 - \frac{1}{\sqrt{1 + \lambda_D \left( \frac{I_0}{2\mu^2 a^3} \right)^2}} \right], \qquad (4.49b)$$

respectively. By expanding (4.49a) we find

$$\bar{\rho}_{\kappa} = \frac{I_0}{8\pi G a^3} + \mathcal{O}\left(\frac{1}{a^6}\right) \quad (\text{when } a \to \infty, \text{ late universe}), \qquad (4.50a)$$

$$\bar{\rho}_{\kappa} = \frac{I_0}{8\pi G a^3} \left(1 + \lambda_D^{-1/2}\right) + \mathcal{O}\left(a^0\right) \quad (\text{when } a \to 0, \text{ early universe}). \qquad (4.50b)$$

The constant $I_0$ is determined by the Khronon energy density today (for given values of the constants $\mu$ and $\lambda_D$), which is in turn related to the measured fraction of dark matter:

$$\bar{\rho}_{\kappa,0} = \frac{1}{8\pi G} \left( I_0 + \frac{2\mu^2}{\lambda_D} \left[ -1 + \sqrt{1 + \lambda_D \left( \frac{I_0}{2\mu^2} \right)^2} \right] \right) = \frac{3H_0^2}{8\pi G} \Omega_{\kappa,0}. \qquad (4.51)$$

Defining the function of the scale factor

$$q(a) \equiv \frac{1}{\sqrt{\lambda_D + \left(\frac{2\mu^2 a^3}{I_0}\right)^2}}, \qquad (4.52)$$

i.e. $q = \bar{Q} - 1$, such that $0 \le q < \lambda_D^{-1/2}$, we find from eqs. (4.49) that the exact evolution of the equation of state and adiabatic sound speed as functions of the scale factor are

$$w = \frac{\lambda_D q^2 + (1+q)\sqrt{1 - \lambda_D q^2} - 1}{2 + (1 + \lambda_D)q}, \qquad (4.53a)$$

$$c_{\text{ad}}^2 = \frac{q(1 - \lambda q^2)}{1 + q}. \qquad (4.53b)$$
Alternatively it is interesting to pose $\mathcal{Z}(a) \equiv \frac{I_0}{2\mu^2 a^3} = q(1 - \lambda_D q^2)^{-1/2}$ and the equation of state and adiabatic sound speed become

$$w = \frac{\mathcal{Z}}{1 + \lambda_D \mathcal{Z}^2 + (1 + \mathcal{Z})\sqrt{1 + \lambda_D \mathcal{Z}^2}}, \quad (4.54a)$$

$$c_{\text{ad}}^2 = \frac{\mathcal{Z}}{(1 + \lambda_D \mathcal{Z}^2) \left[ \sqrt{1 + \lambda_D \mathcal{Z}^2 + \mathcal{Z}} \right]}. \quad (4.54b)$$

The great improvement of the DBI model with respect to the previous “$K_2 + K_{n+1}$” models, in that this model behaves approximatively as dust not only today but also in the early universe, hence it will easily pass the constraints from the CMB. In the late Universe, when $a$ is large, we may expand the above relations to get

$$w = \frac{I_0}{4\mu^2 a^3} + \mathcal{O}\left(\frac{1}{a^6}\right), \quad c_{\text{ad}}^2 = \frac{I_0}{2\mu^2 a^3} + \mathcal{O}\left(\frac{1}{a^6}\right), \quad (4.55a)$$

while in the early Universe, when $a$ is small, they expand as

$$w = \frac{2\mu^2 a^3}{I_0 \sqrt{\lambda_D} (1 + \sqrt{\lambda_D})} + \mathcal{O}(a^6), \quad c_{\text{ad}}^2 = \frac{4\mu^4 a^6}{I_0^2 \lambda_D (1 + \sqrt{\lambda_D})} + \mathcal{O}(a^{12}), \quad (4.55b)$$

respectively. The turning point between the two regimes depends on the parameters $\mu$ and $\lambda_D$. For fixed $\mu$, increasing $\lambda_D$ moves the turning point to later times and vice-versa. Increasing $\mu$, generally moves the turning point to earlier times. This behaviour may be crudely seen by equating the early and late expansions of $w(a)$ above as a first approximation; the precise relations can be found numerically, leading to monotonic dependence of $a_{\text{turn}}$ as a function of the parameters.

The equation of state for the DBI like function is plotted in figure 1 (blue dashed curves) for the two values $\mu^{-1} \sim 22.3$ Mpc and $\mu^{-1} \sim 223$ kpc that are compatible with the MOND phenomenology, i.e. for which the constraint (4.39) is satisfied. As we see from figure 1 the DBI function can be made compatible with the CMB as well, with reasonable values of the parameter $\lambda_D$. For instance the DBI function with $\mu^{-1} \sim 22.3$ Mpc and $\lambda_D \sim 1$ fulfills well the purpose of matching both MOND and the correct cosmology and CMB anisotropies. We note, however, that the constraints from [59, 60] do not directly apply, as they correspond to purely time-dependent sound speed and viscosity, while in our case the sound speed is also $k$-dependent. In this sense, our choice of using the “var-w” model is the most conservative with regards placing rough bounds on the present models. A proper comparison of the Khronon model with data is left for a future investigation.

# 5 Linear stability on Minkowski space

## 5.1 The normal modes

We now turn to the question of the linear stability of small fluctuations on a Minkowski background. We thus expand the metric as

$$g^{\mu\nu} = \eta^{\mu\nu} - h^{\mu\nu}, \quad (5.1)$$
and the Khronon field around a background value $\bar{\tau} = t$, as

$$ \tau = t + \sigma, \qquad (5.2) $$

so that $Q = 1 - \phi + \dot{\sigma}$. Hence, only the quadratic part of the function $\mathcal{K}$ in the limit $Q \to 1$ will contribute to this order in the action, that is, one can consider only $\mathcal{K}(Q) = \mu^2(\phi - \dot{\sigma})^2$. Meanwhile $A_i = \vec{\nabla}_i(\phi - \dot{\sigma})$ so that $\mathcal{Y} = |\vec{\nabla}(\phi - \dot{\sigma})|^2$. Hence, to lowest order the MOND term does not contribute. In order to be general enough, let us set $\mathcal{J} = -\alpha\mathcal{Y} + \cdots$ where $\alpha$ is a constant which in the MOND limit equals 1 and in the Newtonian limit equals 0 (in section 6.1 we discuss the more general case related to the khronometric Hořava gravity). Given this, we find that the action (2.6) when expanded to second order in fluctuations reads

$$ S = \frac{1}{16\pi G} \int d^4 x \left\{ -\frac{1}{2} \bar{\nabla}_{\mu} h \bar{\nabla}_{\nu} h^{\mu\nu} + \frac{1}{4} \bar{\nabla}_{\rho} h \bar{\nabla}^{\rho} h + \frac{1}{2} \bar{\nabla}_{\mu} h^{\mu\rho} \bar{\nabla}_{\nu} h^{\nu}_{\rho} - \frac{1}{4} \bar{\nabla}^{\rho} h^{\mu\nu} \bar{\nabla}_{\rho} h_{\mu\nu} \right. \\ \left. + 2\alpha |\vec{\nabla}(\dot{\sigma} - \phi)|^2 + 2\mu^2 (\dot{\sigma} - \phi)^2 + 8\pi G T_{\mu\nu} h^{\mu\nu} \right\}. \qquad (5.3) $$

Notice that since the only new terms compared to GR are the ones containing the term $\dot{\sigma} - \phi$, that is, only a scalar mode, both tensor modes and vector modes behave exactly like in GR. The former propagate at the speed of light while the latter are non-dynamical.

We thus focus on the scalar mode action. We take the metric from (4.19) and set $a = 1$ and $\kappa = 0$ for a Minkowski background as well as setting the matter source to zero (i.e. $T_{\mu\nu} = 0$) leading to the scalar mode action as<fnref n="12" />

$$ S = \frac{1}{16\pi G} \int d^4 x \left\{ 6 \left( \frac{1}{6} \vec{\nabla}^2 \dot{\nu} - \dot{\psi} \right) \left( \frac{1}{6} \vec{\nabla}^2 \dot{\nu} + \dot{\psi} \right) + 4 \left( \frac{1}{6} \vec{\nabla}^2 \dot{\nu} + \dot{\psi} \right) \vec{\nabla}^2 \zeta + 2 |\vec{\nabla} \psi|^2 - \frac{2}{3} \psi \vec{\nabla}^4 \nu \right. \\ \left. + 4 \left( \vec{\nabla}^2 \psi + \frac{1}{6} \vec{\nabla}^4 \nu \right) \phi + \frac{1}{18} |\vec{\nabla} (\vec{\nabla}^2 \nu)|^2 + 2\alpha |\vec{\nabla} (\dot{\sigma} - \phi)|^2 + 2\mu^2 (\dot{\sigma} - \phi)^2 \right\}. \qquad (5.4) $$

Choosing the Newtonian gauge for which $\zeta = \nu = 0$ (and $\psi = \Psi$, $\phi = \Phi$) and switching to Fourier space by defining the Fourier variables via

$$ \mathcal{O}(\vec{x}) = \sqrt{8\pi G} \int \frac{d^3 k}{(2\pi)^3} e^{i\vec{k} \cdot \vec{x}} \mathcal{O}_{\vec{k}}, \qquad (5.5) $$

for any operator $\mathcal{O}$, leads to

$$ S = \int \frac{dt \, d^3 k}{(2\pi)^3} \left[ -3|\dot{\psi}_{\vec{k}}|^2 + k^2|\psi_{\vec{k}}|^2 - k^2 (\phi_{\vec{k}}^* \psi_{\vec{k}} + \psi_{\vec{k}}^* \phi_{\vec{k}}) + (\alpha k^2 + \mu^2) |\dot{\sigma} - \phi|^2 \right]. \qquad (5.6) $$

To find the normal modes we let all variables be proportional to $e^{i\omega t}$ and, defining the vector in field space through $\mathcal{W} \equiv \{\psi, \phi, \sigma\}$, the above action is written as

$$ S = \int \frac{dt \, d^3 k}{(2\pi)^3} \mathcal{W}^\dagger \mathcal{M} \mathcal{W}, \qquad (5.7) $$

$$ S = \frac{1}{8\pi G} \int d^4 x \left\{ -3\dot{\Psi}^2 + |\vec{\nabla}\Psi|^2 + 2\Phi\vec{\nabla}^2\Psi + \alpha|\vec{\nabla}\Upsilon|^2 + \mu^2\Upsilon^2 \right\}. $$

<footnote n="12">¹²In terms of the gauge-invariant variables (4.22)–(4.23) (with $a = 1$) we have</footnote>
where $\mathcal{M}$ is the matrix

$$ \mathcal{M} = \begin{pmatrix} -3\omega^2 + k^2 & -k^2 & 0 \\ -k^2 & \alpha k^2 + \mu^2 & i(\alpha k^2 + \mu^2)\omega \\ 0 & -i(\alpha k^2 + \mu^2)\omega & (\alpha k^2 + \mu^2)\omega^2 \end{pmatrix}. \qquad (5.8) $$

The normal modes are determined by setting $\det \mathcal{M} = 0$, leading to the condition $\omega^2 = 0$. Hence, just like GR, there are no propagating wavelike modes, that is modes which evolve as $e^{i\omega t}$ with non-zero $\omega$. However, there are two non-propagating modes with $\omega = 0$, one of which is dynamical, that is, of the form $A_{\vec{k}} + B_{\vec{k}}t$ (see [64] for the existence of such modes in generic massive gravity theories). The constant mode $A_{\vec{k}}$ will not lead to any instability and the linearly unstable mode $B_{\vec{k}}t$ only leads to a mild instability (not exponential). By considering the Hamiltonian description below, we show that this is nothing but a Jeans instability which is only exhibited at very low momenta and is thus harmless.

## 5.2 The Hamiltonian formulation

We start again from the scalar mode action (5.4) before fixing the gauge, and switch to Fourier space to get

$$ S = \int \frac{dt \, d^3 k}{(2\pi)^3} \left\{ -3|\dot{\psi}_{\vec{k}}|^2 + \frac{1}{12} k^4 |\dot{\nu}_{\vec{k}}|^2 + k^2 \left[ \left( \frac{1}{6} k^2 \dot{\nu}_{\vec{k}} - \dot{\psi}_{\vec{k}} \right) \zeta_{\vec{k}}^* + \text{c.c.} \right] \right. \\ \left. + k^2 |\psi_{\vec{k}} - \frac{1}{6} k^2 \nu_{\vec{k}}|^2 - k^2 \left[ \left( \psi_{\vec{k}} - \frac{1}{6} k^2 \nu_{\vec{k}} \right) \phi_{\vec{k}}^* + \text{c.c.} \right] + (\alpha k^2 + \mu^2) |\dot{\sigma}_{\vec{k}} - \phi_{\vec{k}}|^2 \right\}, \quad (5.9) $$

where c.c. means complex conjugate of the term within the square bracket. Notice that (5.9) does not contain any time derivatives of the fields $\zeta_{\vec{k}}$ and $\phi_{\vec{k}}$, hence, we expect those to lead to constraints. The canonical momenta are found as

$$ P_{\vec{k}}^{(\psi)} = -2 \left( 3 \dot{\psi}_{\vec{k}} + k^2 \zeta_{\vec{k}} \right), \qquad (5.10) $$

$$ P_{\vec{k}}^{(\nu)} = \frac{1}{6} k^4 \left( \dot{\nu}_{\vec{k}} + 2 \zeta_{\vec{k}} \right), \qquad (5.11) $$

$$ P_{\vec{k}}^{(\sigma)} = 2(\alpha k^2 + \mu^2)(\dot{\sigma}_{\vec{k}} - \phi_{\vec{k}}), \qquad (5.12) $$

and, after a Legendre transformation, the Hamiltonian is found to be

$$ H = \int \frac{d^3 k}{(2\pi)^3} \left\{ -\frac{1}{12} |P_{\vec{k}}^{(\psi)}|^2 + \frac{3}{k^4} |P_{\vec{k}}^{(\nu)}|^2 + \frac{1}{4(\alpha k^2 + \mu^2)} |P_{\vec{k}}^{(\sigma)}|^2 \right. \\ \left. - k^2 |\psi_{\vec{k}} - \frac{1}{6} k^2 \nu_{\vec{k}}|^2 + C_{\zeta} \zeta_{\vec{k}}^* + C_{\phi} \phi_{\vec{k}}^* + \text{c.c.} \right\}, \qquad (5.13) $$

where $C_{\phi}$ and $C_{\zeta}$ are two constraints

$$ C_{\phi} \equiv \frac{1}{2} P_{\vec{k}}^{(\sigma)} + k^2 \psi_{\vec{k}} - \frac{k^4}{6} \nu_{\vec{k}} \approx 0, \qquad (5.14a) $$

$$ C_{\zeta} \equiv -P_{\vec{k}}^{(\nu)} - \frac{k^2}{6} P_{\vec{k}}^{(\psi)} \approx 0, \qquad (5.14b) $$

imposed by the non-dynamical fields $\phi$ and $\zeta$, respectively. As usual we use the symbol $\approx$ to denote weakly vanishing constraints (those that vanish only on-shell).
We require that the constraints are preserved by time evolution according to the Hamiltonian $H = \int \frac{d^3k}{(2\pi)^3} \mathcal{H}$, with $\mathcal{H}$ being the Hamiltonian density in Fourier space. We define the Poisson brackets on phase space as

$$ \{f,g\} = (2\pi)^3 \int d^3k \left[ \sum_A \left( \frac{\delta f}{\delta X^A} \frac{\delta g}{\delta P_{XA}^*} - \frac{\delta g}{\delta X^A} \frac{\delta f}{\delta P_{XA}^*} \right) \right], \quad (5.15) $$

where $X^A \equiv \{\psi_{\vec{k}}, \nu_{\vec{k}}, \sigma_{\vec{k}}\}$ and $P_{XA} \equiv \{P_{\vec{k}}^{(\psi)}, P_{\vec{k}}^{(\nu)}, P_{\vec{k}}^{(\sigma)}\}$. The Poisson brackets define the time evolution of a variable $f$ via $\dot{f} = \{f, \mathcal{H}\}$, and applied to the constraints (5.14) give

$$ \dot{C}_{\phi} = C_{\zeta}, \quad (5.16a) $$

$$ \dot{C}_{\zeta} = 0. \quad (5.16b) $$

In other words, the constraints are preserved by time evolution on-shell. Therefore as one might expect, the stability of the primary constraints in the absence of gauge fixing does not create new constraints. Having ensured the stability of constraints in the Hamiltonian, we can now simplify the system by employing gauge fixing.

In the Hamiltonian formulation, primary first-class constraints generate gauge transformations. The infinitesimal change of a phase space quantity $f$ under this gauge transformation generated by the constraint $C_a$ (with $a = \zeta, \phi$) is given by:

$$ \Delta f = \{f, C_a^*[\epsilon_a]\}, \quad (5.17) $$

where we have introduced the smearing $C_a^*[\epsilon_a]$ of a constraint $C_a^*$ with test function $\epsilon_a$:

$$ C_a^*[\epsilon_a] = \int \frac{d^3k}{(2\pi)^3} \epsilon_{a,\vec{k}} C_{a,\vec{k}}^*. \quad (5.18) $$

Consider the following gauge transformations generated by the constraints $C_\zeta$ and $C_\phi$,

$$ \Delta \nu = \{\nu, C_{\zeta}^*[\epsilon_{\zeta}]\} = -\epsilon_{\zeta}, \quad (5.19a) $$

$$ \Delta P_{\nu} = \{P_{\nu}, C_{\phi}^{*}[\epsilon_{\phi}]\} = \frac{1}{6} k^4 \epsilon_{\phi}. \quad (5.19b) $$

Thus, we may set $\nu$ and $P_\nu$ to zero by a gauge transformation by choosing $\epsilon_\zeta = \nu$ and $\epsilon_\phi = -\frac{6}{k^4}P_\nu$. We then check what constraints are placed on the Lagrange multipliers $\zeta, \phi$ by this gauge fixing. We invoke two new gauge fixing constraints:

$$ G_{\nu} \equiv \nu \approx 0, \quad (5.20a) $$

$$ G_{P_{\nu}} \equiv P_{\nu} \approx 0, \quad (5.20b) $$

and find

$$ \{G_{\nu}, \mathcal{H}\} = \frac{6}{k^4} G_{P_{\nu}} - 2\zeta_{\vec{k}}, \quad (5.21a) $$

$$ \{G_{P_{\nu}}, \mathcal{H}\} = \frac{1}{3} k^4 (\phi_{\vec{k}} - \psi_{\vec{k}}) + \frac{1}{18} k^6 G_{\nu}. \quad (5.21b) $$

Therefore the following gauge restrictions are placed on the Lagrange multipliers: $\zeta = 0$ and $\phi = \psi$. We recognize these conditions, respectively, as a restriction to the Newtonian
gauge and the content of the Einstein field equation here dictating equality between metric potentials in this gauge. We may adopt these conditions alongside the constraints $G_ν$, $G_{Pν}$ in the Hamiltonian (5.13) and the primary constraints, yielding in addition

$$P_{\vec{k}}^{(\sigma)} \approx -2k^2 \psi_{\vec{k}}, \qquad (5.22a)$$

$$P_{\vec{k}}^{(\psi)} \approx 0, \qquad (5.22b)$$

so that the deconstrained Hamiltonian density is

$$\mathcal{H}^{(\text{dec})} = \int \frac{d^3 k}{(2\pi)^3} \frac{(1-\alpha)k^2 - \mu^2}{\alpha k^2 + \mu^2} k^2 |\psi_{\vec{k}}|^2. \qquad (5.23)$$

In the Newtonian case $\alpha = 0$, hence, the Hamiltonian is positive for large momenta, that is for $k \ge \mu$. Thus, the linear instability, discussed in the previous section and due to the $\mu^2(Q-1)^2$ term in (2.6), is a Jeans-type instability occuring only when $k < \mu$. We discuss more extensibly below.

# 6 Discussions

## 6.1 The Hořava khronometric theory

The model we have proposed is reminiscent of, and indeed was inspired [18, 19] by the Hořava-Lifshitz (HL) gravity [43–45]. In eq. (3.23) we found that in order to account for the MONDian low acceleration regime when $Y \to 0$, the function $\mathcal{J}(Y)$ takes essentially the form $\mathcal{J} = -Y + \mathcal{O}(Y^{3/2})$, so that the gravitational part of the Lagrangian in adapted coordinates, see eq. (A.3), is

$$\mathcal{L}_g = \frac{c^4}{16\pi G} \sqrt{q} N \left[ \mathcal{R} + K_{ij} K^{ij} - K^2 + 2Y + \mathcal{O}(Y^{3/2}) \right], \qquad (6.1)$$

where we recall that $Y = A_i A^i / c^4$ in adapted coordinates, and we skip other terms like $\mathcal{K}(1/N)$ for this discussion. The theory (6.1) is only invariant under the subgroup of diffeomorphisms leaving invariant the preferred time foliation. In particular the term $2Y$ implies a breaking of the local Lorentz invariance (LLI).

On the other hand the HL gravity postulates a more general Lagrangian involving three arbitrary constants $\lambda$, $\beta$ and $\alpha$ labeling LLI-violating terms such that

$$\mathcal{L}_g^{\text{HL}} = \frac{c^4}{16\pi G} \sqrt{q} N \left[ \mathcal{R} + (1 - \beta) K_{ij} K^{ij} - (1 + \lambda) K^2 + 2\alpha Y + \dots \right]. \qquad (6.2)$$

The term $2\alpha Y$ was added in [44, 45] to provide stability of the scalar degree of freedom, within the so-called “non-projectable” version of the theory for which the lapse $N$ depends not only on time as in the original “projectable” version [43], but also on space coordinates. In particular the constant $\alpha$ must satisfy [44]

$$0 < \alpha < 1, \qquad (6.3)$$

while a non-zero $\beta$ leads to a speed of the tensor gravitational wave (GW) different from $c$. In addition to the explicit terms in (6.2), the ellipsis contain terms which are of high order
(fourth and sixth) in spatial derivatives but with no time derivatives, and which ensure the power-counting renormalizability of the theory at high energy; these terms are suppressed below some high energy scale. Introducing the Khronon (which is the Stückelberg field associated with the broken diffeomorphism invariance), the covariant version of (6.2) reads

$$ \mathcal{L}_g^{\text{HL}} = \frac{c^4}{16\pi G} \sqrt{-g} \left[ R - \frac{1}{3}(\beta + 3\lambda)\Theta^2 - \beta\sigma_{\mu\nu}\sigma^{\mu\nu} + 2\alpha\mathcal{Y} + \dots \right], \quad (6.4) $$

where $\Theta = K = \nabla_\mu n^\mu$ is the expansion and $\sigma_{\mu\nu} = K_{\mu\nu} - \frac{1}{3}\Theta q_{\mu\nu}$ is the shear of the spacetime congruence $n_\mu$, with $q_{\mu\nu} \equiv g_{\mu\nu} + n_\mu n_\nu$ and the extrinsic curvature $K_{\mu\nu} = \nabla_\mu n_\nu + c^{-2}n_\mu A_\nu$ (symmetric in $\mu\nu$ in the case of a hypersurface-orthogonal congruence). The theory (6.4) is also recognized as a sub-class of Einstein-Aether theories [10, 11] for which the Aether $n_\mu$ is hypersurface-orthogonal.

Taking the model (6.2)–(6.4) at face, the effective Newton’s constant (as measured in a Cavendish type experiment) in the low energy limit, i.e. setting the dots to zero, reads

$$ G_N = \frac{G}{1 - \alpha}, \quad (6.5) $$

while the PPN parameters of the theory are the same as for GR, with the exception of the preferred-frame parameters $\alpha_1$ and $\alpha_2$ given by [45, 65] (see [46] for further discussions)

$$ \alpha_1 = -\frac{8(\alpha - \beta)}{1 - \beta}, \quad (6.6a) $$

$$ \alpha_2 = \frac{(\alpha - \beta)[2\alpha(1 + \beta + 2\lambda) - \beta(3 + \beta + 3\lambda) - \lambda]}{(1 - \alpha)(1 - \beta)(\beta + \lambda)}. \quad (6.6b) $$

Furthermore the cosmological equations (on a FLRW background) are

$$ 3 \left(1 + \frac{\beta}{2} + \frac{3\lambda}{2}\right) H^2 + \frac{3\kappa}{a^2} = 8\pi G \sum_I \bar{\rho}_I, \quad (6.7a) $$

$$ -\left(1 + \frac{\beta}{2} + \frac{3\lambda}{2}\right) \left(2\dot{H} + 3H^2\right) - \frac{\kappa}{a^2} = 8\pi G \sum_I \bar{P}_I. \quad (6.7b) $$

The term $2\alpha\mathcal{Y} \propto A_i A^i$ in the action corresponds to a linear perturbation of the FLRW background and has no effect on (6.7).

The present model (6.1) looks to be a particular case of the HL gravity. However there is the crucial difference that in the HL model the violation of local Lorentz invariance is motivated by the completion of GR at high energy, while in (6.1) the violation of Lorentz invariance is supposed to be active only in the low energy, weak acceleration sector of the theory. In particular the value $\alpha = 1$ which is required in eq. (6.1) seems to be incompatible with the measured Newton’s constant (6.5) and the PPN parameters (6.6). This is because the term $2\mathcal{Y}$ is required to cancel the Newtonian gravity in the low acceleration regime and to replace it by the MOND gravity, thanks to the term $\sim \mathcal{Y}^{3/2}$ in (6.1). However, there is no contradiction, as the measured value for $G_N$ happens in the high acceleration regime, where $\alpha = 0$, following the function $\mathcal{J}(\mathcal{Y})$ given by eq. (3.26). Since we have also $\beta = \lambda = 0$ we recover GR and in particular the same PPN limit as GR.
## 6.2 Degrees of freedom and further discussion about the dispersion relation

At all instances the theory has three dynamical degrees of freedom when expanded on a Minkowski background to quadratic order. However, only two are propagating which are the usual massless tensor modes as in GR. The third degree of freedom is non-propagating (it has $\omega = 0$ and so there is no associated wave), but still dynamical as its solution is of the form $A_{\vec{k}} + B_{\vec{k}}t$. Such modes have been previously identified to exist in generic Lorentz-violating massive gravity theories [64], and also exist in the case of AeST [39, 66], and for the ghost condensate theory [49] in the limit when the higher-derivative interactions are ignored. In our case, there are no higher derivative interaction terms in the action, which are also quadratic in the fields (and so contribute to the usual dispersion relation corresponding to a linear wave equation). However, the higher order terms (for instance, the MOND term), will effectively lead to a non-linear propagation equation for the third degree of freedom, that is, the non-propagating $\omega = 0$ mode will become propagating. We expect the full (non-perturbative) Hamiltonian analysis to reveal the same three dynamical degrees of freedom, similarly to [67] but less complicated and can be used to study the perturbative Hamiltonian dynamics on more general backgrounds.

Let us discuss further the Hamiltonian and the linear instability found above. This is due to the “condensate” term $\mu^2(Q-1)^2$ which equals $\mu^2\phi^2$ in the unitary gauge. In the Newtonian case where $\alpha = 0$, we have that the Hamiltonian is positive for $k > \mu$ but in the MOND case, where $\alpha = 1$ it would seem that the Hamiltonian is always negative even as $k \to \infty$. However, eq. (5.23) is of the same form as the AeST Hamiltonian for the $\omega = 0$ mode analysed in [66]. There it was shown that including the MOND (non-linear) term, leads once more to a Hamiltonian bounded from below for $k \ge \mu$, and still unbounded from below for $k < \mu$. Thus, the same type of Jeans instability persists also at the MOND limit. Taking $\mu^{-1} \sim 22.3$ Mpc as in figure 1 the deconstrained Hamiltonian is bounded from below for wavenumbers larger than $\sim 3 \times 10^{-31}$ eV.

There can be further terms which can influence the conclusions regarding the Hamiltonian. A specific case is the HL extension discussed above, with parameters $\beta$ and $\lambda$ (we set $\beta = 0$ to have GW travelling at the speed of light as required by observations). Including the perturbations coming from (6.4), we find that the third mode in this case is propagating and has dispersion relation

$$ \omega^2 = \frac{\lambda [(1 - \alpha)k^2 - \mu^2]}{(2 + 3\lambda)(\alpha k^2 + \mu^2)} k^2, \qquad (6.8) $$

which in the limit $\mu \to 0$ is consistent with [45]. Notice that the instability structure is unchanged from the $\beta = \lambda = 0$ case of section 5.2, that is, we have that $\omega^2 < 0$ for $\sqrt{1 - \alpha}k < \mu$.<fnref n="13" /> However, in this case, the instability turns from linear to exponential. The corresponding deconstrained Hamiltonian is

$$ \mathcal{H}^{(\text{dec})} = \int \frac{d^3 k}{(2\pi)^3} \left\{ \frac{\lambda}{4(2+3\lambda)} |P_{\vec{k}}^{(\psi)}|^2 + \frac{(1-\alpha)k^2 - \mu^2}{\alpha k^2 + \mu^2} k^2 |\psi_{\vec{k}}|^2 \right\}, \qquad (6.9) $$

<footnote n="13">¹³This seems to depend on the sign of $\lambda$ and $2+3\lambda$. However, $2+3\lambda > 0$ to have a well-defined Friedman equation (6.7) and $\lambda > 0$ to have $\omega^2 > 0$ in the $\sqrt{1-\alpha}k > \mu$ case.</footnote>
which is well-defined in the limit $\lambda \to 0$ and corresponds to (5.23). Assuming that both $\alpha$ and $\lambda$ are small, the maximum instability rate $\Gamma_{\max} = -i\omega_{\max}$ is $\Gamma_{\max} = \mu\sqrt{\lambda/8}$ which we require to be smaller than $H_0$. Of particular interest is the subcase $\lambda = 2\alpha$ for which $\omega^2 = k^2$ (ignoring small $O(\lambda)$ corrections) in the $k \gg \mu$ limit (corresponding to an equivalent case discussed in [45]), that is, the Khronon perturbations are also luminal at high momenta. In this case, the PPN parameters (6.6) are $\alpha_1 = -8\alpha$ while $\alpha_2 = \frac{4\alpha^2}{1-\alpha} \ll \alpha_1$ which is negligible. The observational bounds on $\alpha_1$ impose $\alpha \lesssim 10^{-6}$, and hence, $\mu \lesssim 1000 H_0$, corresponding to $\mu^{-1} \gtrsim 10$ Mpc. This bound is consistent with having a MOND limit and can be weaker for smaller $\alpha_1$.

# 7 Conclusion

We have proposed (extending previous works [18, 19]) an extension of GR to account for the phenomenology of dark matter at the scale of galaxies as summarized by the MOND formula [1–3]. Although this theory cannot be considered as fundamental (there is an arbitrary function in the action which is not explained by fundamental physics), it owns a number of attractive features:

1. It is based on only two dynamical fields: the metric and the scalar Khronon field (plus ordinary matter fields);

2. It recovers MOND at the scale of galaxies (in the weak acceleration regime), with however the restriction to systems being stationary;

3. In the strong acceleration regime the theory agrees with GR and in particular has the same PPN limit as GR for tests in the Solar System and in binary pulsars;

4. The theory has no propagating GW with helicity 0 or helicity 1, so GW are the same as in GR;

5. Last but not least, it can be arbitrarily close to the $\Lambda$-CDM cosmological model at the level of linear cosmological perturbations, where it retrieves the full observed spectrum of CMB anisotropies (for a wide range of parameters $\lambda_D$ and $\mu$).

We note the importance of checking the validity of the MOND limit in galaxies when (3.23) is taken as a low energy effective field theory expansion which comes with a new scale. A naive estimate is that the MOND term in (3.23), that is, $M_p^2|\vec{\nabla}\phi|^3/a_0 \to |\vec{\nabla}\phi_c|^3/\Lambda_0^2$ where $\phi_c \sim \phi/M_p$ is the canonically normalized version of $\phi$, leads to a strong-coupling scale of $\Lambda_0 \sim \sqrt{a_0 M_p}$. Given that $a_0 \sim H_0/6$ then $\Lambda_0 \sim \text{meV}$ corresponding to $\Lambda_0^{-1} \sim 10^{-4}$ m, which is clearly much smaller than the scale of galaxies. However, since there are additional dynamical degrees of freedom in the theory, this naive estimate is bound to change. We leave a more thorough investigation of this estimate to a future study.

We also highlight that to next order in cosmological perturbations, the MOND terms will become relevant and the correspondence with $\Lambda$-CDM will break down. In that case, it is necessary to use N-body simulations to determine the non-linear cosmological large scale structure, and it would be interesting to check where (and how) quasistatic MONDian sources might emerge in such a setup. This is left for a future investigation.
Acknowledgments

We thank Gilles Esposito-Farèse, Eanna Flanagan, Elias Kiritsis, Ignacy Sawicki and Leonardo Trombetta for interesting discussions. We have received support from the Barrande mobility programme (project number 8J21FR028). C.S. acknowledges support by the European Structural and Investment Funds and the Czech Ministry of Education, Youth and Sports (MSMT) (Project CoGraDS-CZ.02.1.01/0.0/0.0/15003/0000437) and by the Royal Society Wolfson Visiting Fellowship “Testing the properties of dark matter with new statistical tools and cosmological data”.

A Formulation in adapted coordinates

It is always possible to adopt a coordinate system for which the coordinate time $t \equiv x^0/c$ is equal to the Khronon field: $t = \tau(x)$. Hence in this coordinate system $Q = 1/N$ where $N = (-g^{00})^{-1/2}$ is the lapse. The unit vector is $n_\mu = (-N, 0)$. Introducing the shift $N_i = g_{0i}$ and the spatial metric $q_{ij} = g_{ij}$, we have the usual 3+1 form for the metric,

$$ds^2 = -c^2 N^2 dt^2 + q_{ij} (dx^i + cN^i dt) (dx^j + cN^j dt). \qquad (A.1)$$

The acceleration takes the 3-dimensional expression (with $A^0 = 0$)

$$A^i = c^2 D^i \ln N, \qquad (A.2a)$$

$$\mathcal{Y} = \frac{A_i A^i}{c^4} = D_i \ln N D^i \ln N, \qquad (A.2b)$$

where $D^i = q^{ij} D_j$ (acting here on a scalar) denotes the covariant derivative compatible with the spatial metric, i.e. $D_i q_{jk} = 0$. The covariant action (2.6) becomes, after discarding a total divergence (with $\sqrt{-g} = N\sqrt{q}$)

$$S = \frac{c^4}{16\pi G} \int dt \, d^3\mathbf{x} \sqrt{q} N \left[ \mathcal{R} + K_{ij} K^{ij} - K^2 - 2\mathcal{J}(\mathcal{Y}) + 2\mathcal{K}(\mathcal{Q}) \right] + S_m[\Psi, N, N_i, q_{ij}], \quad (A.3)$$

where $\mathcal{R}$ is the 3-dimensional scalar curvature of the 3-metric $q_{ij}$ and $K_{ij}$ is the extrinsic curvature

$$K_{ij} = \frac{1}{2N} \left( \frac{1}{c} \partial_t q_{ij} - D_i N_j - D_j N_i \right). \qquad (A.4)$$

In such adapted coordinates the Khronon field has disappeared and the independent dynamical degrees of freedom are geometrical: $N$, $N_i$ and $q_{ij}$ (and the matter fields $\Psi$). See section 6.1 for a generalization of the action inspired by Hořava gravity.

For the variation of the minimally and universally coupled matter action we adopt a notation similar to that in [18]

$$\rho \equiv -\frac{1}{c^2\sqrt{q}} \frac{\delta S_m}{\delta N}, \qquad J^i \equiv \frac{1}{c\sqrt{q}} \frac{\delta S_m}{\delta N_i}, \qquad T^{ij} \equiv \frac{2}{N\sqrt{q}} \frac{\delta S_m}{\delta q_{ij}}, \qquad (A.5)$$

where $\rho$, $J^i$ and $T^{ij}$ reduce to the usual notions of mass density, current density and spatial stresses in the case of a perfect fluid. The variation with respect to the lapse $N$ gives

$$\mathcal{R} + K^2 - K_{ij} K^{ij} - 2\mathcal{J} + 4\mathcal{Y} \mathcal{J}_y + 2\mathcal{K} - 2\mathcal{Q} \mathcal{K}_Q + \frac{4}{c^2} D_i (\mathcal{J}_y A^i) = \frac{16\pi G}{c^2} \rho. \qquad (A.6)$$
Next, the variation with respect to $N_i$ yields the momentum constraint equation as in GR,

$$D_j(K^{ij} - q^{ij}K) = -\frac{8\pi G}{c^3} J^i. \qquad (A.7)$$

Varying with respect to the spatial metric $q_{ij}$ gives

$$\frac{1}{cN} D_t (K^{ij} - q^{ij} K) + \frac{2}{N} D_k [N^{(i} (K^{j)k} - q^{j)k} K)] + 2K^{ik} K_k^j - KK^{ij} - \frac{1}{2} q^{ij} (K^{kl} K_{kl} + K^2) \\ + G^{ij} - \frac{1}{N} (D^i D^j N - q^{ij} D_k D^k N) - \frac{2}{c^4} \mathcal{J} Y A^i A^j + (\mathcal{J} - K) q^{ij} = \frac{8\pi G}{c^4} T^{ij}, \qquad (A.8)$$

where $G^{ij} = R^{ij} - \frac{1}{2}q^{ij}R$ is the 3-dimensional Einstein tensor, and we use the convenient notation $D_t \equiv \partial_t - cN^k D_k$. The trace of eq. (A.8), combined with the constraint equation (A.7), gives

$$-\frac{4}{cN} D_t K - 3 K_{ij} K^{ij} - K^2 - \mathcal{R} + \frac{4}{N} D_i D^i N - 4 \mathcal{Y} \mathcal{J}_Y + 6 (\mathcal{J} - \mathcal{K}) = \frac{16\pi G}{c^4} \left( T_i^i + \frac{2c}{N} N_i J^i \right). \qquad (A.9)$$

Finally, by eliminating $\mathcal{R}$ between (A.6) and (A.9) we obtain the following Poisson-like equation,

$$\frac{1}{c^2} D_i [(1 + \mathcal{J}_Y) A^i] + \mathcal{J} + \mathcal{Y} - \mathcal{K} - \frac{1}{2} \mathcal{Q} \mathcal{K}_\mathcal{Q} - \frac{D_t K}{cN} - K^{ij} K_{ij} = \frac{4\pi G}{c^2} \left( \rho + \frac{2}{cN} N_i J^i + \frac{1}{c^2} T_i^i \right), \qquad (A.10)$$

which is the generalization of eq. (3.13) in the BM theory.

## References

[1] M. Milgrom, *A modification of the Newtonian dynamics as a possible alternative to the hidden mass hypothesis*, Astrophys. J. 270 (1983) 365 [INSPIRE].

[2] M. Milgrom, *A modification of the Newtonian dynamics: Implications for galaxies*, Astrophys. J. 270 (1983) 371 [INSPIRE].

[3] M. Milgrom, *A modification of the Newtonian dynamics: implications for galaxy systems*, Astrophys. J. 270 (1983) 384 [INSPIRE].

[4] B. Famaey and S. McGaugh, *Modified Newtonian Dynamics (MOND): Observational Phenomenology and Relativistic Extensions*, Living Rev. Rel. 15 (2012) 10 [arXiv:1112.3960] [INSPIRE].

[5] J. Bekenstein and M. Milgrom, *Does the missing mass problem signal the breakdown of Newtonian gravity?*, Astrophys. J. 286 (1984) 7 [INSPIRE].

[6] J.D. Bekenstein, *Phase Coupling Gravitation: Symmetries and Gauge Fields*, Phys. Lett. B 202 (1988) 497 [INSPIRE].

[7] R.H. Sanders, *A stratified framework for scalar-tensor theories of modified dynamics*, Astrophys. J. 480 (1997) 492 [astro-ph/9612099] [INSPIRE].

[8] J.D. Bekenstein, *Relativistic gravitation theory for the modified Newtonian dynamics paradigm*, Phys. Rev. D 70 (2004) 083509 [Erratum ibid. 71 (2005) 069901] [astro-ph/0403694] [INSPIRE].

[9] R.H. Sanders, *A tensor-vector-scalar framework for modified dynamics and cosmic dark matter*, Mon. Not. Roy. Astron. Soc. 363 (2005) 459 [astro-ph/0502222] [INSPIRE].
[10] T. Jacobson and D. Mattingly, *Gravity with a dynamical preferred frame*, Phys. Rev. D **64** (2001) 024028 [gr-qc/0007031] [INSPIRE].

[11] T. Jacobson and D. Mattingly, *Einstein-aether waves*, Phys. Rev. D **70** (2004) 024003 [gr-qc/0402005] [INSPIRE].

[12] C. Skordis, *The Tensor-Vector-Scalar theory and its cosmology*, Class. Quant. Grav. **26** (2009) 143001 [arXiv:0903.3602] [INSPIRE].

[13] T.G. Złośnik, P.G. Ferreira and G.D. Starkman, *Modifying gravity with the æther: An alternative to dark matter*, Phys. Rev. D **75** (2007) 044017 [astro-ph/0607411] [INSPIRE].

[14] A. Halle, H.S. Zhao and B. Li, *A nonuniform dark energy fluid: Perturbation equations*, Astrophys. J. Suppl. **177** (2008) 1 [arXiv:0711.0958].

[15] M. Milgrom, *Bimetric MOND gravity*, Phys. Rev. D **80** (2009) 123536 [arXiv:0912.0790] [INSPIRE].

[16] E. Babichev, C. Deffayet and G. Esposito-Farèse, *Improving relativistic MOND with Galileon k-mouflage*, Phys. Rev. D **84** (2011) 061502 [arXiv:1106.2538] [INSPIRE].

[17] C. Deffayet, G. Esposito-Farèse and R.P. Woodard, *Nonlocal metric formulations of MOND with sufficient lensing*, Phys. Rev. D **84** (2011) 124054 [arXiv:1106.4984] [INSPIRE].

[18] L. Blanchet and S. Marsat, *Modified gravity approach based on a preferred time foliation*, Phys. Rev. D **84** (2011) 044056 [arXiv:1107.5264] [INSPIRE].

[19] R.H. Sanders, *Hiding Lorentz Invariance Violation with MOND*, Phys. Rev. D **84** (2011) 084024 [arXiv:1105.3910] [INSPIRE].

[20] S. Mendoza, T. Bernal, J.C. Hidalgo and S. Capozziello, *MOND as the weak-field limit of an extended metric theory of gravity*, AIP Conf. Proc. **1458** (2012) 483 [arXiv:1202.3629] [INSPIRE].

[21] R.P. Woodard, *Nonlocal metric realizations of MOND*, Can. J. Phys. **93** (2015) 242 [arXiv:1403.6763] [INSPIRE].

[22] J. Khoury, *Alternative to particle dark matter*, Phys. Rev. D **91** (2015) 024022 [arXiv:1409.0012] [INSPIRE].

[23] L. Blanchet and L. Heisenberg, *Dark Matter via Massive (bi-)Gravity*, Phys. Rev. D **91** (2015) 103518 [arXiv:1504.00870] [INSPIRE].

[24] S. Hossenfelder, *Covariant version of Verlinde’s emergent gravity*, Phys. Rev. D **95** (2017) 124018 [arXiv:1703.01415] [INSPIRE].

[25] C. Burrage, E.J. Copeland, C. Käding and P. Millington, *Symmetron scalar fields: Modified gravity, dark matter, or both?*, Phys. Rev. D **99** (2019) 043539 [arXiv:1811.12301] [INSPIRE].

[26] M. Milgrom, *Noncovariance at low accelerations as a route to MOND*, Phys. Rev. D **100** (2019) 084039 [arXiv:1908.01691] [INSPIRE].

[27] F. D'Ambrosio, M. Garg and L. Heisenberg, *Non-linear extension of non-metricity scalar for MOND*, Phys. Lett. B **811** (2020) 135970 [arXiv:2004.00888] [INSPIRE].

[28] C. Käding, *Lensing with Generalized Symmetrons*, Astronomy **2** (2023) 128 [arXiv:2304.05875] [INSPIRE].

[29] C. Skordis, D.F. Mota, P.G. Ferreira and C. Boehm, *Large Scale Structure in Bekenstein's theory of relativistic Modified Newtonian Dynamics*, Phys. Rev. Lett. **96** (2006) 011301 [astro-ph/0505519] [INSPIRE].
[30] F. Bourliot, P.G. Ferreira, D.F. Mota and C. Skordis, *The cosmological behavior of Bekenstein's modified theory of gravity*, Phys. Rev. D 75 (2007) 063508 [astro-ph/0611255] [INSPIRE].

[31] S. Dodelson and M. Liguori, *Can Cosmic Structure form without Dark Matter?*, Phys. Rev. Lett. 97 (2006) 231301 [astro-ph/0608602] [INSPIRE].

[32] B. Li, J.D. Barrow, D.F. Mota and H.S. Zhao, *Testing Alternative Theories of Dark Matter with the CMB*, Phys. Rev. D 78 (2008) 064021 [arXiv:0805.4400] [INSPIRE].

[33] C. Skordis, *Generalizing tensor-vector-scalar cosmology*, Phys. Rev. D 77 (2008) 123502 [arXiv:0801.1985] [INSPIRE].

[34] J. Zuntz et al., *Vector field models of modified gravity and the dark sector*, Phys. Rev. D 81 (2010) 104015 [arXiv:1002.0849] [INSPIRE].

[35] X.-D. Xu, B. Wang and P. Zhang, *Testing the tensor-vector-scalar Theory with the latest cosmological observations*, Phys. Rev. D 92 (2015) 083505 [arXiv:1412.4073] [INSPIRE].

[36] D.-C. Dai and D. Stojkovic, *Comment on 'Covariant version of Verlinde's emergent gravity', Phys. Rev. D 96 (2017) 108501 [arXiv:1706.07854] [INSPIRE].

[37] T.G. Złośnik and C. Skordis, *Cosmology of the Galileon extension of Bekenstein's theory of relativistic modified Newtonian dynamics*, Phys. Rev. D 95 (2017) 124023 [arXiv:1702.00683] [INSPIRE].

[38] L. Tan and R.P. Woodard, *Structure Formation in Nonlocal MOND*, JCAP 05 (2018) 037 [arXiv:1804.01669] [INSPIRE].

[39] C. Skordis and T. Złośnik, *New Relativistic Theory for Modified Newtonian Dynamics*, Phys. Rev. Lett. 127 (2021) 161302 [arXiv:2007.00082] [INSPIRE].

[40] C. Skordis and T. Złośnik, *Gravitational alternatives to dark matter with tensor mode speed equaling the speed of light*, Phys. Rev. D 100 (2019) 104013 [arXiv:1905.09465] [INSPIRE].

[41] L. Blanchet and A. Le Tiec, *Model of Dark Matter and Dark Energy Based on Gravitational Polarization*, Phys. Rev. D 78 (2008) 024031 [arXiv:0804.3518] [INSPIRE].

[42] L. Blanchet and A. Le Tiec, *Dipolar Dark Matter and Dark Energy*, Phys. Rev. D 80 (2009) 023524 [arXiv:0901.3114] [INSPIRE].

[43] P. Hořava, *Quantum Gravity at a Lifshitz Point*, Phys. Rev. D 79 (2009) 084008 [arXiv:0901.3775] [INSPIRE].

[44] D. Blas, O. Pujolas and S. Sibiryakov, *Consistent Extension of Hořava Gravity*, Phys. Rev. Lett. 104 (2010) 181302 [arXiv:0909.3525] [INSPIRE].

[45] D. Blas, O. Pujolas and S. Sibiryakov, *Models of non-relativistic quantum gravity: The good, the bad and the healthy*, JHEP 04 (2011) 018 [arXiv:1007.3503] [INSPIRE].

[46] M. Bonetti and E. Barausse, *Post-Newtonian constraints on Lorentz-violating gravity theories with a MOND phenomenology*, Phys. Rev. D 91 (2015) 084053 [Erratum ibid. 93 (2016) 029901] [arXiv:1502.05554] [INSPIRE].

[47] É.É. Flanagan, *Khronometric Theories of Modified Newtonian Dynamics*, Astrophys. J. 958 (2023) 107 [arXiv:2302.14846] [INSPIRE].

[48] R.J. Scherrer, *Purely kinetic k-essence as unified dark matter*, Phys. Rev. Lett. 93 (2004) 011301 [astro-ph/0402316] [INSPIRE].

[49] N. Arkani-Hamed, H.-C. Cheng, M.A. Luty and S. Mukohyama, *Ghost condensation and a consistent infrared modification of gravity*, JHEP 05 (2004) 074 [hep-th/0312099] [INSPIRE].
[50] N. Arkani-Hamed et al., *Dynamics of gravity in a Higgs phase*, JHEP 01 (2007) 036
[hep-ph/0507120] [INSPIRE].

[51] W. Hu, *Structure formation with generalized dark matter*, Astrophys. J. 506 (1998) 485
[astro-ph/9801234] [INSPIRE].

[52] M. Kopp, C. Skordis and D.B. Thomas, *Extensive investigation of the generalized dark matter model*, Phys. Rev. D 94 (2016) 043512 [arXiv:1605.00649] [INSPIRE].

[53] C. Misner, K. Thorne and J. Wheeler, *Gravitation*, Freeman, San Francisco, U.S.A. (1973).

[54] G.W. Horndeski, *Second-order scalar-tensor field equations in a four-dimensional space*, Int. J. Theor. Phys. 10 (1974) 363 [INSPIRE].

[55] D. Langlois, *Dark energy and modified gravity in degenerate higher-order scalar-tensor (DHOST) theories: A review*, Int. J. Mod. Phys. D 28 (2019) 1942006 [arXiv:1811.06271] [INSPIRE].

[56] A. De Felice, S. Mukohyama and K. Takahashi, *Avoidance of Strong Coupling in General Relativity Solutions with a Timelike Scalar Profile in a Class of Ghost-Free Scalar-Tensor Theories*, Phys. Rev. Lett. 129 (2022) 031103 [arXiv:2204.02032] [INSPIRE].

[57] L. Blanchet, *Gravitational polarization and the phenomenology of MOND*, Class. Quant. Grav. 24 (2007) 3529 [astro-ph/0605637] [INSPIRE].

[58] P. Verwayen, C. Skordis and C. Bœhm, *Aether Scalar Tensor (AeST) theory: quasistatic spherical solutions and their phenomenology*, Mon. Not. Roy. Astron. Soc. 531 (2024) 272 [arXiv:2304.05134] [INSPIRE].

[59] M. Kopp, C. Skordis, D.B. Thomas and S. Ilić, *Dark Matter Equation of State through Cosmic History*, Phys. Rev. Lett. 120 (2018) 221102 [arXiv:1802.09541] [INSPIRE].

[60] S. Ilić, M. Kopp, C. Skordis and D.B. Thomas, *Dark matter properties through cosmic history*, Phys. Rev. D 104 (2021) 043520 [arXiv:2004.09572] [INSPIRE].

[61] M. Born and L. Infeld, *Foundations of the new field theory*, Nature 132 (1933) 1004.1 [INSPIRE].

[62] M. Born and L. Infeld, *Foundations of the new field theory*, Proc. Roy. Soc. Lond. A 144 (1934) 425 [INSPIRE].

[63] P.A.M. Dirac, *An extensible model of the electron*, Proc. Roy. Soc. Lond. A 268 (1962) 57
[INSPIRE].

[64] S.L. Dubovsky, *Phases of massive gravity*, JHEP 10 (2004) 076 [hep-th/0409124] [INSPIRE].

[65] D. Blas and H. Sanctuary, *Gravitational Radiation in Hořava Gravity*, Phys. Rev. D 84 (2011) 064004 [arXiv:1105.5149] [INSPIRE].

[66] C. Skordis and T. Złośnik, *Aether scalar tensor theory: Linear stability on Minkowski space*, Phys. Rev. D 106 (2022) 104041 [arXiv:2109.13287] [INSPIRE].

[67] M. Bataki, C. Skordis and T. Złośnik, *Aether scalar-tensor theory: Hamiltonian formalism*, Phys. Rev. D 110 (2024) 044015 [arXiv:2307.15126] [INSPIRE].