New Relativistic Theory for Modified Newtonian Dynamics

Alternative theories of gravity to general relativity (GR) have received immense interest in the past 20 years or so [1, 2]. The driving force behind this interest is not so much that gravity has not been tested in a large region of parameter space [3], but, more importantly, the cosmological systems residing in some parts of that region exhibit behavior from which dark matter (DM) and dark energy (DE), collectively called the dark sector, are inferred.

While most investigations deal with DE, the hypothesis that the DM phenomenon is due to gravitational degrees of freedom (d.o.f.) has received less attention [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. Earliest evidence for the existence of DM [15, 16, 17] was later supported by observations of the motion of stars within galaxies [18, 19]. Milgrom proposed [20, 21, 22] that this could, instead, result from modifying the inertia or dynamics of baryons or the gravitational law at accelerations smaller than $a_{0}\sim 1.2\times 10^{-10}{\rm m/s^{2}}$. The latter is further explored in [23] where if gradients of the potential $\Phi$ are smaller than $a_{0}$, nonrelativistic gravity is effectively governed by

Here, $G_{N}$ is the Newtonian gravitational constant, and $\rho$ the matter density. These models are referred to as modified Newtonian dynamics (MOND).

Much work has gone into deducing astrophysical consequences of MOND, its consistency with data [24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46], and alternative DM based explanations of this law [47, 48, 49, 50] It is inherently nonrelativistic and, thus, difficult to test in cosmological settings (but see [51]) as systems such as the cosmic microwave background (CMB) require a relativistic treatment. CMB physics involves only linearly perturbing a Friedmann-Lemaître-Robertson-Walker (FLRW) background, making it a particularly useful system, devoid of nonlinear modeling systematics, for testing relativistic MOND (RMOND). Relativistic theories that yield MOND behavior have been proposed [23, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67] making clear predictions regarding gravitational lensing and cosmology. In cases where the CMB and matter power spectra (MPS) have been computed, no theory has been shown to fit all of the cosmological data while preserving MOND phenomenology in galaxies [68, 69, 70, 71, 72, 73, 74, 75, 76], (though see [77]).

We present the first RMOND theory which reproduces galactic and lensing phenomenology similar to the Bekenstein-Sanders Tensor-Vector-Scalar (TeVeS) theory [53, 54] and, unlike TeVeS, successfully reproduces the key cosmological observables: CMB and MPS. We describe its construction, discuss its cosmology and show that it is devoid of ghost instabilities. We discuss open questions and possibilities toward its more fundamental grounding.

RMOND theories have always been constructed on phenomenological grounds rather than based on fundamental principles. Quite likely the reason is that the MOND law is empirical, and even the observation that it is scale invariant [78, 79] has not yet led to a definitive conclusion as to how this invariance could lead to a MOND gravitational theory. RMOND theories should obey the principle of general covariance and the Einstein equivalence principle. These are, however, do not provide any guidance as to how RMOND should look like. Indeed, many theories obeying these have nothing to do with MOND, and many RMOND theories obeying these same principles are in conflict with observations. Principle-based MOND theories include [80, 81, 82], however, these are nonrelativistic. Still, the phenomenological approach, that we also follow, can provide valuable guidance toward a more fundamental theory.

What are the necessary phenomenological facts that any successful MOND theory should lead to? It must (i) return to GR (hence, Newtonian gravity) when $\vec{\nabla}\Phi\gg a_{0}$ in quasistatic situations while (ii) reproducing the MOND law (1) when $\vec{\nabla}\Phi\ll a_{0}$. It should also (iii) be in harmony with cosmological observations including the CMB and MPS, (iv) reproduce the observed gravitational lensing of isolated objects without DM halos, and (v) propagate tensor mode gravitational waves (GWs) at the speed of light.

We consider each requirement in turn. Clearly, (i) means that when $|\vec{\nabla}\Phi|\gg a_{0}$, the standard Poisson equation $\vec{\nabla}^{2}\Phi=4\pi G_{N}\rho$ holds while (ii) means that when $|\vec{\nabla}\Phi|\ll a_{0}$ the MOND equation (1) holds. While in many cases [60, 61, 56] the transition between (i) and (ii) depends only on $|\vec{\nabla}\Phi|$, in TeVeS it is facilitated by a scalar d.o.f. $\varphi$. We follow the latter and assume that the physics encapsulated by (i) and (ii) fits within the TeVeS framework.

A template nonrelativistic action then, is

where $\Phi=\hat{\Phi}+\varphi$ is the potential that couples universally to matter, $\hat{G}$ is a constant and ${\cal Y}=|\vec{\nabla}\varphi|^{2}$. The field $\varphi$ obeys    $\vec{\nabla}\cdot[(d{\cal J}/d{\cal Y})\vec{\nabla}\varphi]=4\pi\hat{G}\rho$ while $\hat{\Phi}$ obeys the Poisson equation $\vec{\nabla}^{2}\hat{\Phi}=4\pi\hat{G}\rho$. Emergence of MOND is then ensured if ${\cal J}\rightarrow\frac{2\lambda_{s}}{3(1+\lambda_{s})a_{0}}{\cal Y}^{3/2}$ as $\vec{\nabla}\varphi\rightarrow 0$. It is in this limit that $a_{0}$ appears.

For a point source of mass $M$, the MOND-to-Newton transition occurs at $r_{M}\sim\sqrt{(G_{N}M/a_{0})}$. A MOND force $\sim\sqrt{G_{N}Ma_{0}}/r$ lends its way trivially to a Newtonian force $G_{N}M/r^{2}$ as $r\ll r_{M}$ but in the inner Solar System this is not sufficient. Corrections to $r^{-2}$ due to $\varphi$ will compete with the post-Newtonian force $\sim(G_{N}M)^{2}/r^{3}$, and these are constrained at Mercury’s orbit to less than $\sim 10^{-4}$ [83, 84]. Suppressing these may happen either through screening or tracking. In the former, $\varphi$ is screened at large $\vec{\nabla}\varphi$ so that $\Phi\approx\hat{\Phi}$ while in the latter $\varphi\rightarrow\hat{\Phi}/\lambda_{s}$, so that $G_{N}=(1+1/\lambda_{s})\hat{G}$. We model both with $\lambda_{s}$ since screening is equivalent to $\lambda_{s}\rightarrow\infty$. In terms of ${\cal J}$, tracking happens if ${\cal J}\rightarrow\lambda_{s}{\cal Y}$, while screening occurs if ${\cal J}$ has terms ${\cal Y}^{p}$ with $p\geq 3/2$ (this may be in conflict with Mercury’s orbit even as $p\rightarrow\infty$) or via higher-derivative terms absent from (2).

Consider requirement (iii), that is, successful cosmology. In (2) we have a new d.o.f. $\varphi(\vec{x})$ and we expect that the same will appear in cosmology, albeit with a time dependence, i.e. ${\bar{\phi}}(t)$. Consider a flat FLRW metric so that $g_{00}=-N^{2}$ and $g_{ij}=a^{2}\gamma_{ij}$ where $N(t)$ is the lapse function and $a(t)$ the scale factor. What should the expectation for a cosmological evolution of ${\bar{\phi}}(t)$ be? The MOND law for galaxies is silent regarding this matter. There is, however, another empirical law which concerns cosmology: the existence of sizable amounts of energy density scaling precisely as $a^{-3}$. Within the DM paradigm such a law is a natural consequence of particles obeying the collisionless Boltzmann equation. The validity of this law has been tested [85, 86] and during the time between radiation-matter equality and recombination it is valid within an accuracy of $\sim 10^{-3}$. Do scalar field models leading to energy density scaling as ${\bar{\rho}}\sim a^{-3}$ exist?

The answer is yes: shift symmetric $k$ essence. It has been shown [87] that a scalar field with Lagrangian $\sim{\cal K}(\bar{{\cal X}})$ where $\bar{{\cal X}}=\dot{{\bar{\phi}}}^{2}/N^{2}$, leads to dust (i.e. ${\bar{\rho}}\sim a^{-3}$) plus cosmological constant (CC) solutions provided ${\cal K}(\bar{{\cal X}})$ has a minimum at $\bar{{\cal X}}={\cal X}_{0}\neq 0$. Such a model is the low energy limit of ghost condensation [88, 89] although the latter also contains higher derivative terms $\sim(\square\phi)^{2}$ in its action. The FLRW action is

where $\bar{{\cal Q}}=\dot{{\bar{\phi}}}/N$ and $H=\dot{a}/a$. Interestingly, (2) and (3) are shift symmetric in $\varphi$ and ${\bar{\phi}}$ respectively.

We propose that the MOND analog on FLRW is given by (3) with

where $\Lambda$ is the CC, ${\cal K}_{2}$ and ${\cal Q}_{0}$ parameters and $(\ldots)$ denote higher powers in this expansion. Expanding in ${\cal Q}-{\cal Q}_{0}$ rather than ${\cal X}-{\cal X}_{0}$ is the most general expansion leading to dust solutions and includes the ${\cal K}(\bar{{\cal X}})$ case. The CC in this model remains a freely specifiable parameter, just as in the $\Lambda$-cold dark matter ($\Lambda$CDM) model. Following [88, 89], we call this the (gravitational) Higgs phase.

Requirement (iv), that is, correct gravitational lensing without DM, requires a relativistic theory. A minimal theory for RMOND is a scalar-tensor theory[23] with the scalar providing for a conformal factor between two metrics. However, since null geodesics are unaltered by conformal transformations, such theories cannot produce enough lensing from baryons in the MOND regime. Sanders solved the lensing problem by changing the conformal into a disformal transformation [53] using a unit-timelike vector field, incorporated by Bekenstein [54] into TeVeS. The unit-timelike vector has component $A^{0}\sim\sqrt{-g^{00}}$ and this ensures that the two metric potentials are equal (as in GR), so that solutions which mimic DM also produce the correct light deflection.

Meanwhile the anisotropic scaling of the MOND law $\sim|\vec{\nabla}\varphi|^{3}$ compared with a well-behaved cosmology implying terms like $\dot{{\bar{\phi}}}^{2}$ and $\dot{{\bar{\phi}}}^{4}$, heuristically implies (gravitational) Lorentz violation. A good way of introducing such an ingredient is via a unit-timelike vector field $A_{\mu}$, much like the spirit of the Einstein-Æther theory [90, 91], and TeVeS [53, 54].

The advanced Laser Interferometer Gravitational Observatory (LIGO) and Virgo interferometers [92] observed GWs from a binary neutron star merger. Combined with electromagnetic observations [93, 94], this strongly constrains the GW tensor mode speed to be effectively equal to that of light. By analyzing the tensor mode speed, TeVeS has been shown [95, 96, 97, 98] to be incompatible with the LIGO-Virgo observations for any choice of parameters. The necessary d.o.f. $\phi$ and $A_{\mu}$ are also ingredients of TeVeS, only there, a second metric was introduced as a combination of $g_{\mu\nu}$, $\phi$ and $A_{\mu}$. In [99], $\phi$ and $A_{\mu}$ were combined into a timelike (but not unit) vector $B_{\mu}$, and it was shown that TeVeS may be equivalently formulated with a single metric $g_{\mu\nu}$ minimally coupled to matter, and $B_{\mu}$ with a noncanonical and rather complicated kinetic term. A general class of theories based on the pair $\{g_{\mu\nu},B_{\mu}\}$ was uncovered [98] where the tensor mode speed equals the speed of light in all situations, satisfying requirement (v).

A subset of the general class [98] depends on a scalar $\phi$ and unit-timelike vector $A^{\mu}$ such that [100]

where $F_{\mu\nu}=2\nabla_{[\mu}{A}_{\nu]}$, ${J}_{\mu}={A}^{\alpha}\nabla_{\alpha}{A}_{\mu}$, and the Lagrange multiplier $\lambda$ imposes the unit-timelike constraint on ${A}_{\mu}$. In addition ${\cal F}({\cal Y},{\cal Q})$ is a free function of ${\cal Q}={A}^{\mu}\nabla_{\mu}\phi$ and ${\cal Y}=q^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\phi$ where $q_{\mu\nu}=g_{\mu\nu}+{A}_{\mu}{A}_{\nu}$ is the three-metric orthogonal to ${A}^{\mu}$. Notice that (5) is shift symmetric under $\phi\rightarrow\phi+\phi_{0}$.

On FLRW $\phi={\bar{\phi}}(t)$ while ${A}_{0}=-N$ and ${A}_{i}=0$, hence ${\cal Y}=0$ and ${\cal Q}=\bar{{\cal Q}}$. We define ${\cal K}(\bar{{\cal Q}})=-\frac{1}{2}{\cal F}(0,\bar{{\cal Q}})$ so that (5) turns precisely into (3), which we have argued that it satisfies requirement (iii).

In the weak-field quasistatic limit, we set $g_{00}=-1-2\Psi$ and $g_{ij}=(1-2\Phi)\gamma_{ij}$ and assume that ${A}^{\mu}$ aligns with the time direction so that ${A}^{0}=1-\Psi$ and ${A}^{i}=0$. The scalar is expanded as $\phi={\bar{\phi}}+\varphi$ with $\dot{\varphi}\ll|\vec{\nabla}\varphi|$ and $\dot{{\bar{\phi}}}$ may be set to its (late Universe) FLRW minimum ${\cal Q}_{0}$. Hence, ${\cal Q}=(1-\Psi){\cal Q}_{0}$. Then (5) leads to $\Psi=\Phi$ which can be subbed back to get

where ${\cal J}({\cal Y})={\cal F}({\cal Y},{\cal Q}_{0})/(2-K_{B})$. Compared with (2) a new term appears which looks like a “mass term” for $\Phi$, with $\mu=\sqrt{\frac{2{\cal K}_{2}}{2-K_{B}}}{\cal Q}_{0}$. The solution for $\Phi$ will be as obtained from (2) only for $r\lesssim r_{C}$ where $r_{C}\sim\left(r_{M}\mu^{-2}\right)^{1/3}$, and oscillatory for $r\gtrsim r_{C}$. We require $\mu^{-1}\gtrsim 1{\rm Mpc}$ so that MOND behavior according to (2) may still be attained in galaxies. Thus, the quasistatic limit has at least three parameters: $\lambda_{s}$, $a_{0}$ and $\mu$.

While matter couples only to $\Phi$, gravity comes with two potentials $\Phi$ and $\varphi$ whose action is not diagonal but contains the mixing term $J^{\mu}\nabla_{\mu}\phi\rightarrow\vec{\nabla}\Psi\cdot\vec{\nabla}\varphi$. Without the latter, $\varphi$ decouples and no modification of gravity arises in this situation, apart from $\mu^{2}\Phi$ which is akin to ghost condensation [88, 89]. Diagonalizing by setting $\Phi=\hat{\Phi}+\varphi$ and identifying $\tilde{G}=(1-\frac{K_{B}}{2})\hat{G}$ turns (6) into (2) (plus the $\mu^{2}\Phi^{2}$ term). Since, $\Psi=\Phi$, (6) leads to the right lensing whenever the solution for $\Phi$ mimics DM. This satisfies requirements (i), (ii) and (iv).

The theory just presented was constructed to lead to a FLRW universe resembling $\Lambda$CDM. Given a general ${\cal K}({\cal Q})$, we define the energy density as $8\pi\tilde{G}{\bar{\rho}}={\cal Q}\frac{d{\cal K}}{d{\cal Q}}-{\cal K}$ and pressure as $8\pi\tilde{G}{\bar{P}}={\cal K}$ so that the usual FLRW equations are satisfied. The field equation for ${\bar{\phi}}$ may be integrated once to give $\frac{d{\cal K}}{d{\cal Q}}=\frac{I_{0}}{a^{3}}$ for initial condition $I_{0}$. When ${\cal K}$ obeys the expansion (4), then ${\cal Q}={\cal Q}_{0}+I_{0}/a^{3}+\ldots$, so that ${\bar{\rho}}={\bar{\rho}}_{0}/a^{3}+\ldots$, where $8\pi\tilde{G}{\bar{\rho}}_{0}={\cal Q}_{0}I_{0}$. The pressure is ${\bar{P}}=w_{0}{\bar{\rho}}_{0}/a^{6}+\ldots$ where $w_{0}=\frac{8\pi\tilde{G}{\bar{\rho}}_{0}}{4{\cal Q}_{0}^{2}{\cal K}_{2}}$ is the equation of state at $a=1$, that is, $w=w_{0}/a^{3}+\ldots$ so that ${\bar{P}}=w{\bar{\rho}}$. A time-varying $w$ implies an adiabatic sound speed $c_{\rm ad}^{2}=d{\bar{P}}/d{\bar{\rho}}=\frac{d{\cal K}/d{\cal Q}}{{\cal Q}\;d^{2}{\cal K}/d{\cal Q}^{2}}$ and if ${\cal K}$ obeys (4) then $c_{\rm ad}^{2}=2w_{0}/a^{3}+\ldots$. Clearly, $w\geq 0$ and $c_{\rm ad}^{2}\geq 0$, where the zero point is reached as $a\rightarrow\infty$. As the solution depends on the initial condition $I_{0}$, the density ${\bar{\rho}}$ is not (classically) predicted.

For a proper cosmological matter era in the Higgs phase we need $w_{0}$ to be sufficiently small. Observations [85, 86] give $w\lesssim 0.02$ at $a\sim 10^{-4}$, hence, $w_{0}\lesssim 2\times 10^{-14}$. Meanwhile, $\mu^{-1}\gtrsim{\rm Mpc}$ in order not to spoil the MOND behavior, leading to $w_{0}>\frac{3H_{0}^{2}{\rm Mpc}^{2}\Omega_{0}}{2(2-K_{B})}\gtrsim 10^{-8}$. Unless the effect of the $\mu$ term in (6) is alleviated in some future theory, the Higgs phase cannot be extended too long in the past, and higher terms in (4) must be taken into consideration. Within the present setup, one can arrange this with a function ${\cal K}({\cal Q})$ which suppresses $w$ and $c_{\rm ad}^{2}$ during most of the cosmic evolution. Examples are ${\cal K}=2{\cal K}_{2}{\cal Z}_{0}^{2}\left[\cosh({\cal Z})-1\right]$ (“Cosh function”) and ${\cal K}=2{\cal K}_{2}{\cal Z}_{0}^{2}\left[e^{{\cal Z}^{2}}-1\right]$ (“Exp function”) where ${\cal Z}=({\cal Q}-{\cal Q}_{0})/{\cal Z}_{0}$.

The tight coupling of baryons to photons in the early Universe leads to Silk damping and wipes out all small-scale structure in baryons, preventing the formation of galaxies in the late Universe. Within GR, cold DM sustains the gravitational potentials during the tight coupling period, driving the formation of galaxies and affecting the relative peak heights of the CMB as further corroborated by e.g. the Planck satellite [101]. Checking whether this theory fits the CMB and MPS spectra requires studying linear fluctuations on FLRW.

We consider scalar modes in the Newtonian gauge so that $g_{00}=-(1+2\Psi)$, $g_{0i}=0$ and $g_{ij}=a^{2}(1-2\Phi)\gamma_{ij}$ and perturb the scalar as $\phi={\bar{\phi}}+\varphi$ and the vector as ${A}_{\mu}=\{-1-\Psi,\,\vec{\nabla}_{i}\alpha\}$. The perturbed Einstein, vector and scalar equations, then depend on the new scalar modes $\varphi$ and $\alpha$ and their derivatives. The shear equation remains as in GR, as do the usual perturbed Boltzmann equations for baryon, photons and neutrinos, since they couple only to $g_{\mu\nu}$.

Setting $\chi\equiv\varphi+\dot{{\bar{\phi}}}\alpha$, $\gamma\equiv\dot{\varphi}-\dot{{\bar{\phi}}}\Psi$, $E\equiv\dot{\alpha}+\Psi$ and defining the density contrast $\delta$ and momentum divergence $\theta$ via

the Einstein equations take the same form as in GR, i.e. $\delta G^{0}_{\;\;0}=8\pi G\sum_{I}{\bar{\rho}}_{I}\delta_{I}$ and $\delta G^{0}_{\;\;j}=-8\pi G\sum_{I}({\bar{\rho}}_{I}+{\bar{P}}_{I})\vec{\nabla}_{j}\theta_{I}$ where the index $I$ runs over all matter species including the new variables $\delta$ and $\theta$. These obey standard fluid equations

but with nonstandard pressure contrast:

Hence, the resulting system is not equivalent to a dark fluid: the nonstandard pressure, thus defined, does not close under the fluid variables but, rather, depends on the vector field perturbations $\alpha$ and $E$. The latter evolves with

Cosmologically, the necessary additional free parameters to $\Lambda$CDM are $\lambda_{s}$ (influencing the effective cosmological gravitational strength), $K_{B}$, ${\cal K}_{2}$ (or equivalently $w_{0}$) and ${\cal Q}_{0}$. These fix $\mu$ appearing in the quasistatic regime. More elaborate functions ${\cal K}({\cal Q})$ introduce further parameters, e.g. ${\cal Z}_{0}$ in the case of the “Cosh” or “Exp” functions above. Note that $a_{0}$ does not appear in the linear cosmological regime but will play a role once nonlinear terms from ${\cal F}({\cal Y},{\cal Q})$ kick in.

In Figs. 1 and 2 we show the CMB and MPS in the case of a “Cosh”, an “Exp” and a “Higgs-like” function ${\cal K}({\cal Q})=\frac{{\cal K}_{2}}{4{\cal Q}_{0}^{2}}\left({\cal Q}^{2}-{\cal Q}_{0}^{2}\right)^{2}$, computed numerically by evolving the FLRW background and linearized equations using our own Boltzmann code [102], which is in excellent agreement with other codes, see [103] for a comparison. We have used adiabatic initial conditions [104] and a standard initial power spectrum $P_{0}=A_{s}k^{n_{s}}$ with amplitude $A_{s}$ and spectral index $n_{s}$. The MPS has an additional bias parameter $b$. We used RECFAST version 1.5 for modeling recombination and have boosted $k$ sampling, time sampling and $\ell$ sampling accuracy for ensuring robust results. The detailed cosmology and the dependence of the spectra on the parameters will be investigated elsewhere [104]. For a wide range of parameters, this relativistic MOND theory is consistent with the CMB measurements from Planck. This happens because $c_{\rm ad}^{2}$ and $w$ are small enough so that $\Pi\rightarrow 0$ and we get dustlike evolution as $\dot{\delta}=3\dot{\Phi}-\frac{k^{2}}{a^{2}}\theta$ and $\dot{\theta}=\Psi$, while the vector field decouples.

Now, we consider stability of the theory on Minkowski spacetime. We expand $g^{\mu\nu}=\eta^{\mu\nu}-h^{\mu\nu}$, split ${A}_{\mu}=(-1+\frac{1}{2}h^{00},{A}_{i})$ and let $\nabla_{\mu}\phi=({\cal Q}_{0}+\dot{\varphi},\vec{\nabla}\varphi)$ with $h^{\mu\nu}$, ${A}_{i}$ and $\varphi$ being small perturbations. Expanding (5) to second order gives

where we have used the desired late Universe limit for which $\partial^{2}\bar{{\cal F}}/\partial{\cal Q}^{2}\rightarrow-2d^{2}{\cal K}/d{\cal Q}^{2}=-4{\cal K}_{2}$ and $\partial{\cal F}/\partial{\cal Q}=\bar{{\cal F}}=0$. We set $\partial{\cal F}/\partial{\cal Y}=(2-K_{B})\lambda_{s}$ as a free parameter which is zero in the MOND limit but nonzero in the GR limit when reached by tracking. Inspecting (13), the tensor mode action is as in GR as expected.

For vector modes we choose the gauge $h^{00}=0$, $h^{0i}=W^{i}$ and $h^{ij}=0$ while ${A}_{i}=\beta_{i}$ and $\varphi=0$ where $W^{i}$ and $\beta^{i}$ are transverse. Setting all modes $\propto e^{i(-\omega t+\vec{k}\cdot\vec{x})}$, the dispersion relation for $\beta^{i}$ is $\omega^{2}=k^{2}+{\cal M}^{2}$ where their mass is ${\cal M}^{2}=\frac{(2-K_{B})(1+\lambda_{s}){\cal Q}_{0}^{2}}{K_{B}}$, hence, they are healthy if $0<K_{B}<2$ and $\lambda_{s}>-1$. They decouple from $T_{\mu\nu}$ and are not expected to be generated to leading order by compact objects.

Considering scalar modes in the Newtonian gauge we set $h^{00}=-2\Psi$, $h^{0i}=0$ and $h^{ij}=-2\Phi\gamma^{ij}$ while ${A}_{i}=\vec{\nabla}_{i}\alpha$ and find the dispersion relations $\omega^{2}=0$ and $\omega^{2}=\frac{(2-K_{B})}{{\cal K}_{2}K_{B}}(1+\frac{1}{2}K_{B}\lambda_{s})k^{2}+{\cal M}^{2}$. Thus, we require that ${\cal K}_{2}>0$ in addition to the vector stability conditions. Only two normal modes exist implying the presence of constraints. These are revealed through a Hamiltonian analysis which also shows that these conditions lead to a positive Hamiltonian [106, 107] for the $\omega\neq 0$ modes. The $\omega=0$ case leads to a constant mode with zero Hamiltonian but, also, to a mode varying linearly with $t$. The Hamiltonian for the latter is positive for momenta larger than $\sim\mu$ and otherwise negative, also requiring that $\lambda_{s}>0$. Such instabilities are likely akin to Jeans-type instabilities and do not cause quantum vacuum instability at low momenta [108].

MOND has enjoyed success in fitting galactic rotation curves [24, 25, 27, 29, 44] and reproducing the baryonic Tully-Fisher relation [31, 33, 45]. The radial acceleration relation (RAR) [41] finds a comfortable interpretation within MOND. Studies of MOND with galaxy clusters [30, 109, 110, 38, 111, 112] report that either $a_{0}$ is larger in clusters and/or an additional dark component is necessary even when the MOND prescription is used. These studies, however, use the classic modified-inertia MOND while the theory presented here has additional features warranting its separate testing with clusters. We note that a RAR for clusters was reported [112], similar to the galaxy one albeit with $a_{0}$ a factor of 10 higher. MOND has been tested with dwarf spheroidal galaxies where discrepancies for some [26] were later dismissed with improved data [28, 113, 114, 115, 116]. There, good agreement was reported, except for Draco and Carina where the fits are quite poor [26, 114, 116, 117]. It is argued [113] that those two might be systems not in equilibrium. The global stability of M33 has been tested [118] with positive results while wide-binary data do not yet yield a decisive test [119].

We have shown how the cosmological regime of this theory reproduces the CMB and MPS power spectra on linear scales and that MOND-like behavior emerges in the quasistatic approximation. The latter is expected to hold for virialized objects, however, how such objects emerge from the underlying density field, i.e. how the two regimes connect, is an open problem. This will happen at a scale which is expected to depend on $a_{0}$, $\mu$ and ${\cal Q}_{0}$ and quite likely the nonlinear $\sim\nabla(\nabla\phi)^{2}/a_{0}$ term coming from ${\cal F}$ will play a role. It is reasonable to expect that on mildly nonlinear scales, the quasistatic regime is not yet reached.

We remark that $A_{\mu}$ also contains a pure vector mode perturbation which is expected to behave similarly as in the Einstein-Æther theory [90, 91]. This may lead to imprints on the $B$-mode CMB polarization signal [120].

Setting $\tilde{M}_{p}^{2}=1/(8\pi\tilde{G})$ and canonically normalizing as $\tilde{\phi}=\sqrt{2{\cal K}_{2}}\tilde{M}_{p}\phi$ in (4), the FLRW action (3) becomes

where $\Lambda_{c}^{2}=\tilde{M}_{p}\sqrt{2{\cal K}_{2}}{\cal Q}_{0}$. Considering the MOND limit in (5) gives $\tilde{M}_{p}^{2}{\cal F}/2\rightarrow|\vec{\nabla}\tilde{\phi}|^{3}/\Lambda_{0}^{2}$ where $\Lambda_{0}^{2}=12\left[{\cal K}_{2}(1+1/\lambda_{s})/(2-K_{B})\right]^{3/2}M_{p}a_{0}$. This scale is indicative of the energy scale above which quantum corrections may be important and below which we can trust the classical theory. Since $a_{0}\sim H_{0}/6$ then $\Lambda_{0}\gtrsim{\rm meV}\sim(0.1{\rm mm})^{-1}$. Newton’s $r^{-2}$ law has been tested down to $\sim 52\mu{\rm m}$ [121] and the curves in Figs.1 and 2 have $\Lambda_{0}^{-1}\lesssim 100{\rm nm}$.

Absence of ghosts to quadratic order signifies a healthy theory that could arise as a limit of a more fundamental theory. We do not have such a theory at present but we discuss a case that may bring us closer. The vector in (5) does not seem to obey gauge invariance but in the quadratic action (13) it does so through mixing with diffeomorphisms of $h_{\mu\nu}$. This is not an accident. Let us normalize via $\hat{A}_{\mu}=M_{{\rm GGC}}{A}_{\mu}$ for some scale $M_{{\rm GGC}}$ and insert the term $-\frac{1}{4}\frac{\tilde{M}_{p}^{4}}{M_{{\rm GGC}}^{4}}\lambda^{2}$. Varying with $\lambda$ and using the constraint to eliminate $\lambda$ from the action, perform a Stückelberg transformation $\hat{A}_{\mu}\rightarrow\hat{A}_{\mu}+\nabla_{\mu}\xi/M_{{\rm GGC}}$ and define the covariant derivative acting on “angular field” $\xi$ as ${\cal D}_{\mu}\xi=\nabla_{\mu}\xi/M_{{\rm GGC}}+\hat{A}_{\mu}$. The action turns to $S=S_{EH}+\int d^{4}x\sqrt{-g}\bigg{\{}-\frac{1}{4g_{{\rm GGC}}^{2}}\hat{F}_{\mu\nu}\hat{F}^{\mu\nu}+\frac{1}{4}({\cal D}^{\mu}\xi{\cal D}_{\mu}\xi+M_{{\rm GGC}}^{2})^{2}\bigg{\}}$ plus $\phi$-dependent terms, where $\hat{F}_{\mu\nu}=\nabla_{\mu}\hat{A}_{\nu}-\nabla_{\nu}\hat{A}_{\mu}$, $g_{{\rm GGC}}^{2}=\frac{M_{{\rm GGC}}^{2}}{K_{B}\tilde{M}_{p}^{2}}$. The resulting action is that of the gauged ghost condensate (GGC) [122] or bumblebee field [123, 124] which has been proposed as a healthy gauge-invariant theory of spontaneous Lorentz violation. The Einstein-Æther theory, part of (5), is the (healthy) decoupling limit of GGC by taking $M_{{\rm GGC}}\rightarrow\infty$ if $0<K_{B}<2$ (in our notation) [122]. It is argued [122] that $M_{{\rm GGC}}$ can be as high as $10^{12}{\rm GeV}$.

Given that $\phi$ is shift symmetric it is natural to charge it under this symmetry similar to $\xi$ letting ${\cal D}_{\mu}\phi=\nabla_{\mu}\phi/M_{{\rm GGC}}+\hat{A}_{\mu}$. Interestingly, we may identify ${\cal Q}-{\cal Q}_{0}\rightarrow{\cal D}^{\mu}\xi{\cal D}_{\mu}\phi$ while the term ${J}^{\mu}\nabla_{\mu}\phi\rightarrow\hat{F}^{\mu\nu}{\cal D}_{\mu}\xi{\cal D}_{\nu}\phi$, both multiplied by appropriate constants. The terms involving ${\cal Y}$ may be constructed using $\left(g^{\mu\nu}+{\cal D}^{\mu}\xi{\cal D}^{\nu}\xi/M_{{\rm GGC}}^{4}\right){\cal D}_{\mu}\phi{\cal D}_{\nu}\phi$. Although extending our work as such does not explain the MOND term ${\cal Y}^{3/2}$, it may provide promising directions for further improvements.