# Aether scalar tensor theory: Hamiltonian Formalism

Marianthi Bataki<EMAIL_ADDRESS>Department of Physics, University
of Cyprus 1, Panepistimiou Street, 2109, Aglantzia, Cyprus  CEICO, Institute
of Physics of the Czech Academy of Sciences, Na Slovance 1999/2, 182 00,
Prague, Czechia  Constantinos Skordis<EMAIL_ADDRESS>CEICO, Institute of
Physics of the Czech Academy of Sciences, Na Slovance 1999/2, 182 00, Prague,
Czechia  Tom Zlosnik<EMAIL_ADDRESS>Institute of Theoretical
Physics and Astrophysics, University of Gdańsk, 80-308 Gdańsk, Poland

###### Abstract

The Aether Scalar Tensor (AeST) theory is an extension of General Relativity
(GR), proposed for addressing galactic and cosmological observations without
dark matter. By casting the AeST theory into a $3+1$ form, we determine its
full non-perturbative Hamiltonian formulation and analyse the resulting
constraints. We find the presence of four first class and four second class
constraints and show that the theory propagates six degrees of freedom on
general backgrounds. We determine the case of small perturbations around
Minkowski space as an example stemming from our analysis.

##  I Introduction

Einstein’s theory of General Relativity (GR) has been extraordinarily
successful and remains the paradigmatic theory of gravity. Data from a wide
range of astrophysical systems, from precision tests of gravity in the solar
system Will (2014) to recent measurements of gravitational waves Abbott  _et
al._ (2017) are in accord with its predictions. Although there are expected to
be corrections to the theory that become important at extremely high energy
scales/short length scales (for instance in the early universe Starobinsky
(1980)), a question mark hangs over whether corrections or additional
structure in the gravitational sector may also become important on
cosmological and certain astrophysical scales. The reason for this is the
considerable evidence that matter in many astrophysical and cosmological
systems moves as if under the gravitational pull of a greater amount of matter
than seems to exist in those systems. A notable early example of this was the
observation of the motion of stars orbiting in the outskirts of spiral
galaxies Rubin  _et al._ (1980), the rotational velocity of the stars being
too large to be due to the gravity of visible matter interior to their orbits.
Similar ‘mass discrepancies’ have since been consistently seen in other
galaxies Lelli  _et al._ (2019), clusters of galaxies Hahn  _et al._ (2007),
and on the largest cosmological scales Aghanim  _et al._ (2020).

The mainstream explanation for these effects is that there exists an
additional, dark matter in the universe, which at present outweighs known
matter by a factor of over five Aghanim  _et al._ (2020). Despite there being
many proposed candidates for what this matter may be, see Jungman  _et al._
(1996); Bertone  _et al._ (2005); Marsh (2016); Green and Kavanagh (2021) for
reviews, evidence for it remains comprised entirely of its gravitational
effect on known matter. As such, there remains the possibility that what is
being observed may not be the effect of the presence of additional (dark)
matter, but that of additional gravitational degrees of freedom leading to a
modification of the way known matter affects the gravitational field. For this
to manifest, an extension of GR must be at play.

Modified Newtonian Dynamics (MOND) is a non-relativistic framework proposed by
Milgrom Milgrom (1983), as a way of addressing galactic observations without
dark matter. In one formulation, Newton’s second law of motion is changed
below an acceleration scale $a_{0}\sim 1.2\times 10^{-10}m/s^{2}$ while non-
relativistic gravity is governed by Poisson’s equation. In another
formulation, Newton’s second law is kept but the gravitational equation
determining the non-relativistic potential $\Phi$ from the matter density is
generalized and departs from Poisson’s equation at low potential gradients
determined by $a_{0}$. Several instances of the second formulation exist,
starting from a single potential formulation of Bekenstein and Milgrom
Bekenstein and Milgrom (1984), to other later formulations involving
additional potentials Milgrom (2010, 2023).

The MOND proposal leads to a wide variety of predictions that can be compared
against astrophysical data Famaey and McGaugh (2012); Banik and Zhao (2022)
though has the inherent restriction that in the absence of a relativistic
completion, it has not been clear what the realm of validity of MOND is. This
has prompted the construction of a variety of extensions of GR Bekenstein and
Milgrom (1984); Bekenstein (1988); Sanders (1997); Bekenstein (2004); Navarro
and Van Acoleyen (2006); Zlosnik  _et al._ (2007); Sanders (2005); Milgrom
(2009); Babichev  _et al._ (2011); Deffayet  _et al._ (2011); Blanchet and
Marsat (2011); Sanders (2011); Mendoza  _et al._ (2012); Woodard (2015);
Khoury (2015); Blanchet and Heisenberg (2015); Hossenfelder (2017); Burrage
_et al._ (2019); Milgrom (2019); D’Ambrosio  _et al._ (2020); Käding (2023),
which lead to MOND behaviour at the quasistatic weak-field limit111 Extended
dark matter models Blanchet (2007); Blanchet and Le Tiec (2009); Berezhiani
and Khoury (2016, 2015); Kaplinghat  _et al._ (2016); Kamada  _et al._ (2017)
have also been proposed to accommodate for MOND behaviour, keeping GR as the
description of gravity so that only two massless tensor-mode gravitational
waves are present. These models automatically have a $\Lambda$CDM cosmology in
accordance with observations.. However, none of these extensions have been
shown to fit the observations of the Cosmic Microwave Background radiation as
most recently reported by the Planck Surveyor satellite Aghanim  _et al._
(2020). Moreover, they usually do not lead to a gravitational wave tensor mode
speed equalling the speed of light, although exceptions can be found, as
required by data Abbott  _et al._ (2017); Savchenko  _et al._ (2017).

Starting from a general class of theories based on a metric, a unit-timelike
vector field and a scalar field Skordis and Złośnik (2019), the AeST theory
Skordis and Zlosnik (2021) was constructed to have a massless spin-2 graviton
which propagates at the speed of light while tending to MOND in the weak-field
quasistatic regime relevant to galaxies. It was further required to have a
Friedman-Lemaitre-Robertson-Walker (FLRW) behaviour extremely close to that of
$\Lambda$CDM by having features akin to shift-symmetric k-essense Scherrer
(2004) and ghost condensate theory Arkani-Hamed  _et al._ (2004, 2007). The
closeness to $\Lambda$CDM persists also when linear fluctuations around FLRW
are included, which enable the theory to provide a similarly good match to
precision cosmological data (for example the linear matter power spectrum and
CMB temperature and polarization anisotropy power spectra) to the $\Lambda$CDM
model. It was further shown that linear fluctuations on Minkowski space
propagate two massless tensor modes, two massive vector modes and one massive
scalar mode, all of which are healthy provided certain constraints on the
theory parameters are satisfied Skordis and Zlosnik (2022). A sixth mode was
shown to have a linear $t$-dependence and to have positive Hamiltonian for
momenta larger than a mass scale which observationally is $\lesssim
10^{-30}eV$ and negative otherwise. This behaviour is akin to a Jeans
instability and do not cause quantum vacuum instability at low momenta
Gumrukcuoglu  _et al._ (2016). Further studies of the AeST theory have been
performed in Bernardo and Chen (2023); Kashfi and Roshan (2022); Mistele
(2022); Mistele  _et al._ (2023); Tian  _et al._ (2023); Llinares (2023);
Verwayen  _et al._ (2023); Mistele (2023).

Despite these promising features of the theory, it is crucial that it can
match the success of GR in all cases where it has been tested, while fitting
observations in the regimes where GR fails without the addition of dark
matter. This will require finding solutions to the theory in systems that
might not be describable by linear perturbations propagating on highly
symmetric backgrounds. Towards these ends, an important first step will be to
cast the equations as first order evolution equations in time - this will
enable both analytical and numerical solutions for more complicated to be more
easily found.

In this paper we develop the canonical/Hamiltonian formulation of the theory
and this will allow us to put the theory’s equations of motion in the form of
Hamilton’s first-order equations of motion. The completion of the canonical
analysis also enables clarification of other issues, such as, the number of
degrees of freedom that the theory possesses and whether the theory is an
example of an irregular system, that is, a theory where the canonical
structure varies throughout phase space. A manifestation of the latter can be
that perturbations around some backgrounds describe different number of
degrees of freedom than perturbations around other backgrounds, see Alexandrov
_et al._ (2013, 2021) for examples.

The outline of the paper is as follows: In Section II we introduce the theory;
in Section III we introduce the so-called 3+1 formalism and its associated
decomposition of fields; in Section IV we cast the theory in Hamiltonian form
and perform a full constraint analysis in section V; in Section VI we restrict
the full non-perturbative Hamiltonian to the case of small perturbations
around a Minkowski space solution and in doing so demonstrate the recovery of
the results previously found in Skordis and Zlosnik (2022). Finally, in
Section VII we present our conclusions.

##  II The theory

The theory depends on a metric $g_{\mu\nu}$ universally coupled to matter so
that the Einstein equivalence principle is obeyed, a scalar field $\phi$ and a
unit time-like vector field $\hat{A}^{\mu}$ 222We depart from previous
expositions of this theory Skordis and Zlosnik (2021, 2022) and denote the
4-dimensional vector field by $\hat{A}^{\mu}$. We reserve the symbol $A^{\mu}$
(without the ‘hat’) for the projected vector field on the 3-dimensional
hypersurface introduced in section III, see equation (10), as this will
feature more prominently than $\hat{A}^{\mu}$ in the present work.. The action
is

$\displaystyle S=$ $\displaystyle\int
d^{4}x\frac{\sqrt{-g}}{16\pi\tilde{G}}\bigg{\\{}R-2\Lambda-\frac{K_{{B}}}{2}\hat{F}^{\mu\nu}\hat{F}_{\mu\nu}+(2-K_{{B}})\left(2J^{\mu}\nabla_{\mu}\phi-{\cal
Y}\right)-{\cal F}({\cal Y},{\cal
Q})-\lambda(\hat{A}^{\mu}\hat{A}_{\mu}+1)\bigg{\\}}+S_{m}[g]$ (1)

where $g$ is the metric determinant, $\nabla_{\mu}$ the covariant derivative
compatible with $g_{\mu\nu}$, $R$ is the Ricci scalar, $\Lambda$ is the
cosmological constant, $\tilde{G}$ is the bare gravitational strength,
$K_{{B}}$ is a constant and $\lambda$ is a Lagrange multiplier imposing the
unit time-like constraint on $A_{\mu}$. In addition, we have defined the
tensors $J^{\mu}=\hat{A}^{\nu}\nabla_{\nu}\hat{A}^{\mu}$ and
$\hat{F}_{\mu\nu}=2\nabla_{[\mu}\hat{A}_{\nu]}$ while the matter action
$S_{m}$ is assumed not to depend explicitly on $\phi$ or $\hat{A}^{\mu}$.
Defining

$\displaystyle\hat{q}_{\mu\nu}\equiv$ $\displaystyle
g_{\mu\nu}+\hat{A}_{\mu}\hat{A}_{\nu}$ (2)

the function ${\cal F}({\cal Y},{\cal Q})$ depends on the scalars

$\displaystyle{\cal Q}\equiv$ $\displaystyle\hat{A}^{\mu}\nabla_{\mu}\phi$ (3)

and

$\displaystyle{\cal Y}\equiv$
$\displaystyle\hat{q}^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\phi$ (4)

The function ${\cal F}$ is subject to conditions so that the cosmology of the
theory is compatible with $\Lambda$CDM on FRLW spacetimes and a MOND limit
emerges in quasistatic situations.

On a flat FLRW background the metric takes the form
$ds^{2}=-dt^{2}+a^{2}\gamma_{ij}dx^{i}dx^{j}$ where $a(t)$ is the scale factor
and $\gamma_{ij}$ is a flat spatial metric. The vector field reduces to
$\hat{A}^{\mu}=(1,0,0,0)$ while $\phi\rightarrow{\bar{\phi}}(t)$ leading to
${\cal Q}\rightarrow\bar{{\cal Q}}=\dot{{\bar{\phi}}}$ and ${\cal
Y}\rightarrow 0$, so that we may define ${\cal K}(\bar{{\cal
Q}})\equiv-\frac{1}{2}{\cal F}(0,\bar{{\cal Q}})$. We require that ${\cal
K}(\bar{{\cal Q}})$ has a minimum at ${\cal Q}_{0}$ (a constant) so that we
may expand it as ${\cal K}={\cal K}_{2}\left(\bar{{\cal Q}}-{\cal
Q}_{0}\right)^{2}+\ldots$, where the $(\ldots)$ denote higher terms in this
Taylor expansion. This condition leads to ${\bar{\phi}}$ contributing energy
density scaling as dust $\sim a^{-3}$ akin to Scherrer (2004); Arkani-Hamed
_et al._ (2004), plus small corrections which tend to zero when
$a\rightarrow\infty$. In principle, ${\cal K}$ could be offsetted from zero at
the minimum ${\cal Q}_{0}$, i.e. ${\cal K}({\cal Q}_{0})={\cal K}_{0}$,
however, such an offset can always be absorbed into the cosmological constant
$\Lambda$ and thus we choose ${\cal K}_{0}=0$ by convention, implying the same
on the parent function ${\cal F}$.

In the quasistatic weak-field limit we may set the scalar time derivative to
be at the minimum ${\cal Q}_{0}$, as is expected to be the case in the late
universe. This means that we may expand $\phi={\cal Q}_{0}t+\varphi$.
Moreover, in this limit ${\cal F}\rightarrow(2-K_{{B}}){\cal J}({\cal Y})$,
with ${\cal J}$ defined appropriately as ${\cal J}({\cal
Y})\equiv\frac{1}{2-K_{{B}}}{\cal F}({\cal Y},{\cal Q}_{0})$. It turns out
that MOND behaviour emerges if ${\cal
J}\rightarrow\frac{2\lambda_{s}}{3(1+\lambda_{s})a_{0}}|{\cal Y}|^{3/2}$ where
$a_{0}$ is Milgrom’s constant and $\lambda_{s}$ is a constant which is related
to the Newtonian/GR limit. Specifically, there are two ways that GR can be
restored: (i) screening and (ii) tracking. In the former, the scalar is
screened at large gradients $\vec{\nabla}\varphi$, where
$\vec{\nabla}\leftrightarrow\vec{\nabla}_{i}$ is the spatial gradient on a
flat background $\gamma_{ij}$, and in the latter, $\lambda_{s}\varphi$ becomes
proportional to the Newtonian potential, leading to an effective Newtonian
constant

$G_{{N}}=\frac{1+\frac{1}{\lambda_{s}}}{1-\frac{K_{{B}}}{2}}\tilde{G}.$ (5)

Screening may be achieved either through terms in ${\cal J}\sim{\cal Y}^{p}$
with $p>3/2$ or through galileon-type terms which must be added to (1). Either
way, for our purposes in this article, we may model screening as
$\lambda_{s}\rightarrow\infty$.

##  III 3+1 formalism

###  III.1 ADM decomposition

####  III.1.1 Decomposition of the metric and its derivatives

As a necessary first step towards constructing the Hamiltonian formalism for
the theory, we must make a distinction between space and time. Specifically,
we assume that for the region of spacetime of interest, there exists a global
time coordinate $t(x^{\mu})$. Given this, we may define a ‘flow of time’
vector field $t^{\mu}$ which satisfies $t^{\mu}\nabla_{\mu}t=1$. We use the
notation $\dot{f}\equiv\partial_{t}f$ for some field $f$. Furthermore, we may
define a vector field $n^{\mu}$ which is normal to surfaces of constant $t$;
as such, this field is timelike and may be defined so that it has unit-norm
i.e. $g_{\mu\nu}n^{\mu}n^{\nu}=-1$. We may expand this time field as
$t^{\mu}=Nn^{\mu}+N^{\mu}$, where we have introduced the lapse function $N$
and shift vector $N^{\mu}$, which are given respectively by
$N=-t^{\mu}n_{\mu}$ and $N_{\mu}=q_{\mu\nu}t^{\nu}$. We coordinatize surfaces
of constant $t$ by spatial coordinates $x^{i}$, where $i,j,k$ will be used
throughout to denote spatial coordinate indices. The full spacetime metric may
be decomposed as:

$g_{\mu\nu}=-n_{\mu}n_{\nu}+q_{\mu\nu},$ (6)

where $q_{\mu\nu}$ is the metric on the spatial hypersurface (and therefore,
for example $q_{\mu\nu}n^{\mu}=0$). Note the difference between
$\hat{q}_{\mu\nu}$ defined in (2) and $q_{\mu\nu}$ defined in (6).

It is useful to define a spatial derivative
$\hat{\partial}_{\mu}=q_{\mu}^{\;\;\nu}\partial_{\nu}$ and covariant
derivative $D_{\mu}$ compatible wiht $q_{\mu\nu}$, i.e.
$D_{\alpha}q_{\mu\nu}=0$. Specifically:

$D_{\mu}q_{\alpha\beta}=\hat{\partial}_{\mu}q_{\alpha\beta}-\gamma^{\nu}_{\mu\alpha}g_{\nu\beta}-\gamma^{\nu}_{\mu\beta}g_{\alpha\nu}=0$
(7)

where $\gamma^{\mu}_{\alpha\beta}$ are Levi-Civita symbols associated with the
metric $q_{\mu\nu}$ and derivative $\hat{\partial}_{\mu}$. Given $q_{\mu\nu}$
and $n^{\mu}$, a further useful quantity is extrinsic curvature tensor
$K_{\mu\nu}$, defined as

$\displaystyle K_{\mu\nu}\equiv\frac{1}{2}{\cal
L}_{n}q_{\mu\nu}=q_{\mu}^{\;\;\alpha}\nabla_{\alpha}n_{\nu}$ (8)

In component form, we need

$\displaystyle K_{ij}$
$\displaystyle=\frac{1}{2N}\left(\dot{q}_{ij}-D_{i}N_{j}-D_{j}N_{i}\right)$
(9)

while the components of the metric, $n^{\mu}$ and the Christoffel connection
are displayed in appendix A.

We will adhere to the convention that spatial indices are always lowered and
raised with the spatial metric $q_{ij}$, i.e. $K^{i}_{\;\;j}=q^{ik}K_{kj}$.

####  III.1.2 Decomposition of the vector field

For the vector field $\hat{A}_{\mu}$ we consider a similar decomposition

$\hat{A}_{\mu}=\chi n_{\mu}+A_{\mu}$ (10)

where $A_{\mu}\equiv q_{\mu}^{\phantom{\mu}\nu}\hat{A}_{\nu}$ and

$\chi=-n^{\mu}\hat{A}_{\mu}$ (11)

In component form we find

$\displaystyle\hat{A}_{0}=$ $\displaystyle-N\chi+N^{i}A_{i}\;,$
$\displaystyle\hat{A}_{i}=$ $\displaystyle A_{i}$ (12)
$\displaystyle\hat{A}^{0}=$ $\displaystyle\frac{\chi}{N}\;,$
$\displaystyle\hat{A}^{i}=$ $\displaystyle A^{i}-\frac{\chi}{N}N^{i}$ (13)

leading to $N^{i}A_{i}=A_{0}$, while

$\chi=\sqrt{1+|\vec{A}|^{2}}$ (14)

where $|\vec{A}|^{2}\equiv A^{i}A_{i}=q^{ij}A_{i}A_{j}$, and have taken the
positive sign of the square root by convention.

###  III.2 The $3+1$ action

We now present the necessary steps in writing the action (1) in $3+1$ form.
One of the steps involves solving for $\dot{\phi}$ in terms of the canonical
momenta $\delta S/\delta\dot{\phi}$ and this will involve having to invert
potentially very complicated combinations of funtions $\partial{\cal
F}/\partial{\cal Y}$ and $\partial{\cal F}/\partial{\cal Q}$. Instead we can
move this structure into elsewhere in the theory by introducing auxiliary
fields $\mu$ and $\nu$ such that we set

$\displaystyle{\cal F}({\cal Y},{\cal Q})=-\nu{\cal Q}^{2}+\mu{\cal Y}+{\cal
U}(\nu,\mu)$ (15)

For the scalar field $\phi$, we then find the scalars ${\cal Q}$ and ${\cal
Y}$ as

$\displaystyle{\cal Q}$ $\displaystyle=\chi\sigma+A^{i}D_{i}\phi$ (16)
$\displaystyle{\cal Y}$ $\displaystyle=|\vec{A}|^{2}\sigma^{2}+2\chi\sigma
A^{i}D_{i}\phi+\left(q^{ij}+A^{i}A^{j}\right)D_{i}\phi D_{j}\phi$ (17)

where we have defined

$\displaystyle\sigma$
$\displaystyle=\frac{1}{N}\left(\dot{\phi}-N^{i}D_{i}\phi\right)$ (18)

Consider now the vector-dependent terms in (1) involving $\hat{F}_{\mu\nu}$
and $J^{\mu}$. These depend on the derivatives of $\hat{A}^{\mu}$ which are
displayed in appendix A. Using those relations and letting

$F_{ij}\equiv 2D_{[i}A_{j]}=\hat{F}_{ij}$ (19)

and

$\displaystyle F_{i}\equiv\frac{1}{N}\hat{F}_{0i}=$
$\displaystyle\frac{1}{N}\left[\dot{A}_{i}+D_{i}\left(N\chi-N^{j}A_{j}\right)\right],$
(20)

we define the “magnetic” aspect of $A^{i}$ as

$\displaystyle B^{k}=\frac{1}{2}\epsilon^{kij}F_{ij},$ (21)

with inverse $F_{ij}=\epsilon_{ijk}B^{k}$ and the “electric” aspect of $A^{i}$
as

$\displaystyle E_{i}=F_{i}+\frac{1}{N}\epsilon_{ijk}N^{j}B^{k}.$ (22)

With the above relations and again using (134) and (135) from appendix A we
find

$\displaystyle J^{0}=$ $\displaystyle\frac{1}{N}\vec{A}\cdot\vec{E}$ (23)
$\displaystyle J^{i}=$ $\displaystyle\chi
E^{i}-\frac{\vec{A}\cdot\vec{E}}{N}N^{i}-\epsilon^{ijk}A_{j}B_{k}$ (24)

so that the 3+1 form of (1) as

$\displaystyle S=$ $\displaystyle\int
d^{4}x\frac{N\sqrt{q}}{16\pi\tilde{G}}\bigg{\\{}{\cal
R}+|K^{2}|-|K|^{2}-2\Lambda+K_{{B}}\left(|\vec{E}|^{2}-|\vec{B}|^{2}\right)+2(2-K_{{B}})\sigma\vec{A}\cdot\vec{E}+2(2-K_{{B}})\left(\chi\vec{E}-\vec{A}\times\vec{B}\right)\cdot\vec{D}\phi$
$\displaystyle+\nu\left(\chi\sigma+\vec{A}\cdot\vec{D}\phi\right)^{2}-\left(2-K_{{B}}+\mu\right)\left[|\vec{A}|^{2}\sigma^{2}+2\chi\sigma\vec{A}\cdot\vec{D}\phi+|\vec{D}\phi|^{2}+\left(\vec{A}\cdot\vec{D}\phi\right)^{2}\right]-{\cal
U}(\nu,\mu)\bigg{\\}}+S_{m}[g]$ (25)

where ${\cal R}$ is the Ricci scalar corresponding to the spatial metric
$q_{ij}$, and $S$ is a functional of $(q_{ij},A_{i},\phi,\mu,\nu,N,N^{i})$,

##  IV Hamiltonian formulation

The first step in passing to the Hamiltonian formulation is to determine the
canonical momenta which are

$\displaystyle\Pi^{ij}\equiv$ $\displaystyle\frac{\delta
S}{\delta\dot{q}_{ij}}=\frac{\sqrt{q}}{16\pi\tilde{G}}\left(K^{ij}-Kq^{ij}\right)$
(26) $\displaystyle\Pi^{i}\equiv$ $\displaystyle\frac{\delta
S}{\delta\dot{A}_{i}}=\frac{\sqrt{q}}{8\pi\tilde{G}}\bigg{[}K_{{B}}E^{i}+(2-K_{{B}})\left(\sigma
A^{i}+\chi D^{i}\phi\right)\bigg{]}$ (27) $\displaystyle\Pi\equiv$
$\displaystyle\frac{\delta
S}{\delta\dot{\phi}}=\frac{\sqrt{q}}{8\pi\tilde{G}}\bigg{[}(2-K_{{B}})\vec{A}\cdot\vec{E}+\nu\sigma$
$\displaystyle-\left(2-K_{{B}}+\mu-\nu\right)\vec{A}\cdot\left(\sigma\vec{A}+\chi\vec{D}\phi\right)\bigg{]}$
(28)

Letting $\hat{\Pi}\equiv\Pi^{ij}q_{ij}$ the inverse relations are

$\displaystyle K^{ij}=$
$\displaystyle\frac{16\pi\tilde{G}}{\sqrt{q}}\left(\Pi^{ij}-\frac{1}{2}\hat{\Pi}q^{ij}\right)$
(29) $\displaystyle\Xi\sigma=$
$\displaystyle\frac{8\pi\tilde{G}}{\sqrt{q}}\left[\Pi-\frac{2-K_{{B}}}{K_{{B}}}\vec{A}\cdot\vec{\Pi}\right]$
$\displaystyle+\left(2\frac{2-K_{{B}}}{K_{{B}}}+\mu-\nu\right)\chi\vec{A}\cdot\vec{D}\phi$
(30) $\displaystyle K_{{B}}E^{i}=$
$\displaystyle\frac{8\pi\tilde{G}}{\sqrt{q}}\Pi^{i}-(2-K_{{B}})\left(\sigma
A^{i}+\chi D^{i}\phi\right)$ (31)

where

$\Xi=\chi^{2}\nu-\left(2\frac{2-K_{{B}}}{K_{{B}}}+\mu\right)|\vec{A}|^{2}$
(32)

while the canonical momenta for $\mu$ and $\nu$ are identically zero:

$\displaystyle\Pi^{(\mu)}\equiv$ $\displaystyle\frac{\delta
S}{\delta\dot{\mu}}\approx 0$ (33) $\displaystyle\Pi^{(\nu)}\equiv$
$\displaystyle\frac{\delta S}{\delta\dot{\nu}}\approx 0$ (34)

Using (29), (30) and (31) to remove $K_{ij}$, $E_{i}$ and $\sigma$ from (25)
leads to the Hamiltonian form of the action,

$\displaystyle S$ $\displaystyle=\int
d^{4}x\bigg{\\{}\Pi^{ij}\dot{q}_{ij}+\Pi^{i}\dot{A}_{i}+\Pi\dot{\phi}+\Pi^{(\mu)}\dot{\mu}+\Pi^{(\nu)}\dot{\nu}$
$\displaystyle-N{\cal H}-N^{i}{\cal
H}_{i}-\lambda^{(\mu)}\Pi^{(\mu)}-\lambda^{(\nu)}\Pi^{(\nu)}\bigg{\\}},$ (35)

where we have added Lagrange multipliers $\lambda^{(\mu)}$ and
$\lambda^{(\nu)}$ imposing the constraints (33) and (34) and where

$\displaystyle{\cal H}_{i}=$ $\displaystyle-2D_{j}\Pi^{j}_{\;\;i}+\Pi
D_{i}\phi-\vec{D}\cdot\vec{\Pi}\;A_{i}-\epsilon_{ijk}\Pi^{j}B^{k}$
$\displaystyle+\Pi^{(\mu)}D_{i}\mu+\Pi^{(\nu)}D_{i}\nu$ (36)

is the diffeomorphism constraint and

$\displaystyle{\cal H}=$
$\displaystyle\frac{8\pi\tilde{G}}{\sqrt{q}}\left[2\Pi^{ij}\Pi_{ij}-\hat{\Pi}^{2}+\frac{1}{2K_{{B}}}|\vec{\Pi}|^{2}+\frac{C_{1}^{2}}{2\Xi}\right]+\chi\left[\frac{C_{1}C_{2}}{\Xi}\vec{A}\cdot\vec{D}\phi+\vec{D}\cdot\vec{\Pi}-\frac{2-K_{{B}}}{K_{{B}}}\vec{\Pi}\cdot\vec{D}\phi\right]+\frac{\sqrt{q}}{16\pi\tilde{G}}\bigg{\\{}-{\cal
R}+2\Lambda+K_{{B}}|\vec{B}|^{2}$
$\displaystyle+\left[\frac{C_{2}^{2}\chi^{2}}{\Xi}+2-K_{{B}}+\mu-\nu\right]\left[\vec{A}\cdot\vec{D}\phi\right]^{2}+2(2-K_{{B}})\vec{A}\times\vec{B}\cdot\vec{D}\phi+\left[2-K_{{B}}+\mu+\frac{(2-K_{{B}})^{2}}{K_{{B}}}\chi^{2}\right]|\vec{D}\phi|^{2}+{\cal
U}\bigg{\\}},$ (37)

the Hamiltonian constraint. We have defined

$\displaystyle C_{1}\equiv\Pi-\frac{2-K_{{B}}}{K_{{B}}}\vec{A}\cdot\vec{\Pi},$
(38) $\displaystyle C_{2}\equiv 2\frac{2-K_{{B}}}{K_{{B}}}+\mu-\nu.$ (39)

Note that to arrive at the action (35) we have defined the coefficients
multiplying $(\Pi^{(\mu)},\Pi^{(\nu)})$ to be
$(\lambda^{(\mu)}+N^{i}D_{i}\mu,\lambda^{(\nu)}+N^{i}D_{i}\nu)$ which we are
free to do at the beginning of the constraint analysis. This leads to the 2nd
line in (36).

The combination

${\cal H}_{\rm pri}=N{\cal H}+N^{i}{\cal
H}_{i}+\lambda^{(\mu)}\Pi^{(\mu)}+\lambda^{(\nu)}\Pi^{(\nu)}$ (40)

is the _primary_ Hamiltonian density and it is a sum of constraints on the
phase space which is coordinatized by
$\\{q_{ij},\Pi^{ij},A_{i},\Pi^{i},\varphi,\Pi,\mu,\Pi^{(\mu)},\nu,\Pi^{(\nu)}\\}$.
These constraints are obtained by varying (35) with
$\\{N,N^{i},\lambda^{(\mu)},\lambda^{(\nu)}\\}$ and are given respectively by:

$\displaystyle{\cal H}$ $\displaystyle\approx 0$ (41a) $\displaystyle{\cal
H}_{i}$ $\displaystyle\approx 0$ (41b) $\displaystyle\Pi^{(\mu)}$
$\displaystyle\approx 0$ (41c) $\displaystyle\Pi^{(\nu)}$
$\displaystyle\approx 0$ (41d)

where $\approx$ denotes that these phase space functions need only vanish on
the submanifold of phase space that they define (‘weakly vanish’). The next
step is to check whether these constraints are preserved by the time evolution
generated by ${\cal H}_{\rm pri}$.

##  V The propagation of constraints

###  V.1 Poisson brackets

For a phase space coordinatized by fields $Q_{I}(t,{\vec{x}})$ and
$P_{I}(t,{\vec{x}})$ (where here indices $I,J,\dots$ label the different
fields), it is useful to introduce the Poisson bracket defined for quantities
$A(Q_{I},P_{J})$, $B(Q_{I},P_{J})$. If $\tau_{i\ldots j}^{k\ldots l}(x)$ is a
general tensor field and ${\cal F}[\tau_{i\ldots j}^{k\ldots l}(x)]$ a
functional of $\tau_{i\ldots j}^{k\ldots l}(x)$, then the functional
derivative of ${\cal F}[\tau_{i\ldots j}^{k\ldots l}(x)]$ wrt $\tau_{i\ldots
j}^{k\ldots l}(x)$, in three dimensions is defined as

$\displaystyle\frac{\delta{\cal F}[\tau_{i\ldots j}^{k\ldots
l}({\vec{x}})]}{\delta\tau_{a\ldots b}^{c\ldots d}({\vec{y}})}$
$\displaystyle=\lim_{\epsilon\rightarrow 0}\frac{1}{\epsilon}\bigg{\\{}{\cal
F}\big{[}\tau_{i\ldots j}^{k\ldots l}({\vec{x}})$
$\displaystyle+\epsilon\delta^{(3)}({\vec{y}}-{\vec{x}})\delta^{k}_{S\\{c}\dots\delta^{l}_{d\\}}\delta^{S\\{a}_{i}\dots\delta^{b\\}}_{j}\big{]}$
$\displaystyle-{\cal F}\left[\tau_{i\ldots j}^{k\ldots
l}({\vec{x}})\right]\bigg{\\}}$ (42)

with $\delta^{(3)}({\vec{y}}-{\vec{x}})$ being the three-dimensional Dirac
delta-function and $S\\{ab..cd\\}$ applies the symmetries of the tensor field
$\tau$ e.g. if $\tau^{ijkl}=\tau^{jikl}$, $\tau^{ijkl}=-\tau^{ijlk}$ then
$S\\{ijkl\\}=(ij)[kl]$. With this definition, the Poisson bracket of
$A(Q_{I},P_{J})$ and $B(Q_{I},P_{J})$ is defined as

$\displaystyle\left\\{A,B\right\\}\equiv\sum_{I}\int d^{3}x\left(\frac{\delta
A}{\delta Q_{I}}\frac{\delta B}{\delta P_{I}}-\frac{\delta A}{\delta
P_{I}}\frac{\delta B}{\delta Q_{I}}\right)$ (43)

For time evolution according to any general Hamiltonian $H$ (corresponding to
general Hamiltonian density $\mathscr{H}$)

$H=\int d^{3}x\mathscr{H}$ (44)

we have Hamilton’s equations for quantities $f(Q_{I},P^{I})$ on phase space

$\displaystyle\dot{f}$ $\displaystyle=\left\\{f,H\right\\}.$ (45)

In our case, we test the time evolution of the constraints according to
$H_{\rm pri}$, that is, letting ${\cal C}_{I}$ being any of the constraints in
the set $\\{{\cal H},{\cal H}_{i},\Pi^{(\mu)},\Pi^{(\nu)}\\}$, we require that
$\dot{{\cal C}}_{I}\approx 0$ with $H\rightarrow H_{\rm pri}$ in (44) and
(45).

It is generally more straightforward to evaluate the Poisson bracket of
_smeared_ constraints i.e. we define, for some ‘test’ function $N({\vec{x}})$:

$\displaystyle{{\cal H}}[{N}]\equiv\int d^{3}yN({\vec{y}})\;{\cal
H}[Q_{I}({\vec{y}}),P^{I}({\vec{y}})].$ (46)

With the above definition, (40) is rewritten as

$H_{\rm pri}={{\cal H}}[{N}]+{{\cal
H}_{i}}[{N^{i}}]+{\Pi^{(\mu)}}[{\lambda^{(\mu)}}]+{\Pi^{(\nu)}}[{\lambda^{(\nu)}}].$
(47)

To proceed with the evaluation of Poisson brackets, it is useful to have at
hand the following two results. If ${\cal H}$ depends algebraically on $Q_{I}$
then $\frac{\delta{\cal H}[N]}{\delta
Q_{I}({\vec{x}})}=N({\vec{x}})\frac{\partial{\cal H}}{\partial
Q_{I}}({\vec{x}})$. Furthermore, given arbitrary tensor fields
$\sigma^{i\ldots j}_{k\ldots l}$ and $\tau_{i\ldots j}^{k\ldots l}$ then

$\displaystyle\frac{\delta}{\delta\tau_{b_{1}\ldots b_{s}}^{a_{1}\ldots
a_{r}}}\int d^{3}y\,\sigma^{d_{1}\ldots d_{s}}_{c_{1}\ldots c_{r}}\,{\cal
L}_{\vec{N}}\tau_{d_{1}\ldots d_{s}}^{c_{1}\ldots c_{r}}$
$\displaystyle=-{\cal L}_{\vec{N}}\sigma^{d_{1}\ldots d_{s}}_{c_{1}\ldots
c_{r}}\delta_{S\\{a_{1}}^{c_{1}}\ldots\delta_{a_{r}\\}}^{c_{r}}\delta^{S\\{b_{1}}_{d_{1}}\ldots\delta^{b_{s}\\}}_{d_{s}}$
(48)

where ${\cal L}_{\vec{N}}$ is the Lie derivative according to the vector field
$\vec{N}$.

To satisfy (45) there are ten Poisson brackets to be evaluated among the
smeared constraints ${\cal C}_{I}$. Of these, three are trivial, namely,
$\left\\{\Pi^{(\mu)},\Pi^{(\mu)}\right\\}=\left\\{\Pi^{(\nu)},\Pi^{(\nu)}\right\\}=\left\\{\Pi^{(\mu)},\Pi^{(\nu)}\right\\}=0$
vanish strongly. Let us now evaluate the remaining seven.

###  V.2 Poisson brackets involving a diffeomorphism constraint

Starting from (36) it can be shown that

$\displaystyle{\vec{{\cal H}}}[{\vec{N}}]$ $\displaystyle\overset{b}{=}\int
d^{3}x\big{(}\Pi^{ij}{\cal L}_{\vec{N}}q_{ij}+\Pi{\cal
L}_{\vec{N}}\phi+\Pi^{i}{\cal L}_{\vec{N}}A_{i}$
$\displaystyle+\Pi^{(\mu)}{\cal L}_{\vec{N}}\mu+\Pi^{(\nu)}{\cal
L}_{\vec{N}}\nu\big{)}$ (49)

where $\overset{b}{=}$ means equal to up to a boundary term. As in Romano
(1993), given a general _smeared tensor density_ ${\tau_{i\ldots j}^{k\ldots
l}}[{\sigma^{i\ldots j}_{k\ldots l}}]$ one finds

$\displaystyle\left\\{{\vec{{\cal H}}}[{\vec{N}}],{\tau_{i\ldots j}^{k\ldots
l}}[{\sigma^{i\ldots j}_{k\ldots l}}]\right\\}=$ $\displaystyle{\tau_{i\ldots
j}^{k\ldots l}}[{\,{\cal L}_{\vec{N}}\,\sigma^{i\ldots j}_{k\ldots l}\,}]$
(50)

and choosing $\sigma\rightarrow N^{i}$ and $\tau\rightarrow{\cal H}_{i}$ we
find the Poisson bracket between two diffeomorphism constraints as

$\displaystyle\left\\{{\vec{{\cal H}}}[{\vec{N}}],{\vec{{\cal
H}}}[{\vec{N}^{\prime}}]\right\\}=$ $\displaystyle{\vec{{\cal H}}}[{{\cal
L}_{\vec{N}}\vec{N}^{\prime}}]$ (51)

Similarly choosing $\sigma\rightarrow N$ and $\tau\rightarrow{\cal H}$ we find
the Poisson bracket between the diffeomorphism and the Hamiltonian constraint
as

$\displaystyle\left\\{{\vec{{\cal H}}}[{\vec{N}}],{{\cal H}}[{N}]\right\\}=$
$\displaystyle{{\cal H}}[{{\cal L}_{\vec{N}}N}]$ (52)

These two Poisson brackets are identical to the ones for GR, however, the
phase-space is now enlarged due to the additional fields $A_{i}$, $\phi$,
$\mu$, and $\nu$.

Lastly, setting $\sigma\rightarrow\lambda^{(A)}$ and
$\tau\rightarrow\Pi^{(A)}$ (with $A$ denoting either $\mu$ or $\nu$) trivially
gives $\left\\{{\vec{{\cal
H}}}[{\vec{N}}],{\Pi^{(\mu)}}[{\lambda^{(\mu)}}]\right\\}={\Pi^{(\mu)}}[{{\cal
L}_{\vec{N}}\lambda^{(\mu)}}]$ and $\left\\{{\vec{{\cal
H}}}[{\vec{N}}],{\Pi^{(\nu)}}[{\lambda^{(\nu)}}]\right\\}={\Pi^{(\nu)}}[{{\cal
L}_{\vec{N}}\lambda^{(\nu)}}]$ which vanish weakly.

###  V.3 Poisson brackets of two Hamiltonian constraints

Now consider the Poisson bracket of two Hamiltonian constraints, that is
$\left\\{{{\cal H}}[{N}],{{\cal H}}[{N^{\prime}}]\right\\}$. We let

$\displaystyle\frac{\delta{{\cal H}}[{N}]}{\delta Q_{I}}=$ $\displaystyle
N{\cal A}_{(Q_{I})}+{\cal A}_{(Q_{I})}^{i}D_{i}N$
$\displaystyle+\frac{\sqrt{q}}{16\pi\tilde{G}}\delta_{(Q_{I},q_{ij})}\left[q^{ij}\vec{D}^{2}N-D^{i}D^{j}N\right]$
(53) $\displaystyle\frac{\delta{{\cal H}}[{N}]}{\delta P_{I}}=$ $\displaystyle
N{\cal B}_{(P_{I})}+{\cal B}_{(P_{I})}^{i}D_{i}N$ (54)

where $Q_{I}=\\{q_{ij},A_{i},\phi\\}$, $P_{I}=\\{\Pi^{ij},\Pi^{i},\Pi\\}$,
$\delta_{(Q_{I},q_{ij})}=1$ if $Q_{I}\rightarrow q_{ij}$ and zero otherwise,
and the exact form of ${\cal A}_{(Q_{I})}$, ${\cal A}_{(Q_{I})}^{i}$, ${\cal
B}_{(P_{I})}$ and ${\cal B}_{(P_{I})}^{i}$ is displayed in appendix B.

Using equation (43) directly and making use of (53) and (54), the cross terms
proportional to $NN^{\prime}$ vanish due to the antisymmetry of the Poisson
bracket. In addition, from appendix B we have that ${\cal
A}_{(q_{ij})}^{k}={\cal B}^{i}_{(\Pi^{jk})}={\cal B}^{i}_{(\Pi)}=0$ while
${\cal A}_{(A_{j})}^{i}=-{\cal A}_{(A_{i})}^{j}$ and ${\cal
B}^{i}_{(\Pi^{j})}=-\chi\delta^{i}_{\;\;j}$, so that after some integrations
by parts we find

$\displaystyle\left\\{{{\cal H}}[{N}],{{\cal H}}[{N^{\prime}}]\right\\}$
$\displaystyle=\int d^{3}x\bigg{\\{}\chi\left[D_{j}{\cal
A}_{(A_{k})}^{j}-{\cal A}_{(A_{k})}\right]$ $\displaystyle\
+\frac{\sqrt{q}}{16\pi\tilde{G}}\left[q^{ij}D^{k}{\cal
B}_{(\Pi^{ij})}-q^{jk}D^{i}{\cal B}_{(\Pi^{ij})}\right]$ $\displaystyle\
-{\cal A}_{(A_{i})}^{k}\left[D_{i}\chi+{\cal B}_{(\Pi^{i})}\right]$
$\displaystyle\ -{\cal B}_{(\Pi)}{\cal
A}_{(\phi)}^{k}\bigg{\\}}\left(ND_{k}N^{\prime}-N^{\prime}D_{k}N\right)$ (55)

Finally, plugging in all the expressions from appendix B leads to

$\displaystyle\left\\{{{\cal H}}[{N}],{{\cal H}}[{N^{\prime}}]\right\\}$
$\displaystyle=\int
d^{3}x\bigg{[}-2D_{i}\Pi^{ik}+\Pi\vec{D}^{k}\phi-A^{k}\vec{D}\cdot\vec{\Pi}$
$\displaystyle+\Pi^{i}D^{k}A_{i}-\Pi^{i}D_{i}A^{k}\bigg{]}\left(ND_{k}N^{\prime}-N^{\prime}D_{k}N\right)$
(56)

so that after integrations by parts we get

$\displaystyle\left\\{{{\cal H}}[{N}],{{\cal
H}}[{N^{\prime}}]\right\\}\approx{{\cal
H}_{i}}[{ND^{i}N^{\prime}-N^{\prime}D^{i}N}]$ (57)

Once again, this Poisson bracket is identical to the ones for GR albeit for an
enlarged phase-space.

###  V.4 Secondary constraints

Now consider the Poisson bracket of ${{\cal H}}[{N}]$ with
${\Pi^{(\mu)}}[{\lambda^{(\mu)}}]$ and ${\Pi^{(\nu)}}[{\lambda^{(\nu)}}]$
leading to

$\displaystyle\left\\{{{\cal
H}}[{N}],{\Pi^{(\mu)}}[{\lambda^{(\mu)}}]\right\\}$ $\displaystyle=\int
d^{3}xN\lambda^{(\mu)}\frac{\partial{{\cal H}}}{\partial\mu}$ (58)

and

$\displaystyle\left\\{{{\cal
H}}[{N}],{\Pi^{(\nu)}}[{\lambda^{(\nu)}}]\right\\}$ $\displaystyle=\int
d^{3}xN\lambda^{(\nu)}\frac{\partial{{\cal H}}}{\partial\nu}$ (59)

respectively. Hence, for general $N$, $\lambda^{(\mu)}$ and $\lambda^{(\nu)}$,
the requirement these Poisson brackets vanish weakly implies the additional
secondary constraints on phase space,

$\displaystyle{\cal S}^{(\mu)}\equiv\frac{\partial{{\cal
H}}}{\partial\mu}=\frac{\sqrt{q}}{16\pi\tilde{G}}\left({\cal
Y}+\frac{\partial{\cal U}}{\partial\mu}\right)$ $\displaystyle\approx 0$ (60)

and

$\displaystyle{\cal S}^{(\nu)}\equiv\frac{\partial{{\cal
H}}}{\partial\nu}=\frac{\sqrt{q}}{16\pi\tilde{G}}\left(-{\cal
Q}^{2}+\frac{\partial{\cal U}}{\partial\nu}\right)$ $\displaystyle\approx 0$
(61)

where remember that $\Xi=\Xi(\mu,\nu)$ from (32). Equations (60) and (61)
correspond to the Euler-Lagrange equations for the auxiliary fields $\mu$ and
$\nu$. In other words, given a prescribed ${\cal U}(\mu,\nu)$ one can use (60)
and (61) to determine $\mu({\cal Q},{\cal Y})$ and $\nu({\cal Q},{\cal Y})$,
where ${\cal Q}$ and ${\cal Y}$ are to be evaluated in phase space. Using (30)
followed by (17) and (16) to collect terms together, one finds that

$\displaystyle{\cal Q}$
$\displaystyle=\frac{1}{\Xi}\left[\left(2\frac{2-K_{{B}}}{K_{{B}}}+\mu\right)\vec{A}\cdot\vec{D}\phi+\frac{8\pi\tilde{G}}{\sqrt{q}}\chi
C_{1}\right]$ (62)

and

$\displaystyle{\cal Y}$
$\displaystyle=|\vec{D}\phi|^{2}+\left(\vec{A}\cdot\vec{D}\phi\right)^{2}+\frac{|\vec{A}|^{2}}{\Xi^{2}}\left[\frac{8\pi\tilde{G}}{\sqrt{q}}C_{1}+C_{2}\chi\vec{A}\cdot\vec{D}\phi\right]^{2}$
$\displaystyle+\frac{2\chi}{\Xi}\left(\frac{8\pi\tilde{G}}{\sqrt{q}}C_{1}+C_{2}\chi\vec{A}\cdot\vec{D}\phi\right)\vec{A}\cdot\vec{D}\phi,$
(63)

respectively. This procedure then reconstructs ${\cal F}({\cal Q},{\cal Y})$
through (15) in terms of phase space variables with the help of (62) and (63).

| Primary | Secondary
---|---|---
| ${\cal H}_{i}$ | ${\cal H}$ | $\Pi^{(\mu)}$ | $\Pi^{(\nu)}$ | ${\cal S}^{(\mu)}$ | ${\cal S}^{(\nu)}$
${\cal H}_{i}$ | ${\cal H}_{i}$ | ${\cal H}$ | $\Pi^{(\mu)}$ | $\Pi^{(\nu)}$ | ${\cal S}^{(\mu)}$ | ${\cal S}^{(\nu)}$
${\cal H}$ |  | ${\cal H}_{i}$ | ${\cal S}^{(\mu)}$ | ${\cal S}^{(\nu)}$ | $U^{(\mu)}$ | $U^{(\nu)}$
$\Pi^{(\mu)}$ |  |  | $0$ | $0$ | $C^{(\mu)(\mu)}$ | $C^{(\nu)(\mu)}$
$\Pi^{(\nu)}$ |  |  |  | $0$ | $C^{(\mu)(\nu)}$ | $C^{(\nu)(\nu)}$
${\cal S}^{(\nu)}$ |  |  |  |  | $E^{(\mu)(\mu)}$ | $E^{(\mu)(\nu)}$
${\cal S}^{(\mu)}$ |  |  |  |  |  | $E^{(\nu)(\nu)}$
| First Class | Second Class
Table 1: Table of constraints and their Poisson Brackets, showing whether they
vanish strongly, weakly, or not at all. The classification into
primary/secondary and first/second class is marked. Note also that the
combination $H_{\rm FC}$ defined through (74) is first class, with
$\lambda^{(A)}$ being functions of all the phase space variables as determined
through (73), even though some individual parts of $H_{\rm FC}$ are second
class.

Having found the secondary constraints (60) and (61), the analysis is not
necessarily finished. We now define the secondary Hamiltonian through

$H_{\rm sec}=H_{\rm pri}+\int d^{3}x\left[u^{(\mu)}{\cal
S}^{(\mu)}+u^{(\nu)}{\cal S}^{(\nu)}\right]$ (64)

where $u^{(\mu)}$ and $u^{(\nu)}$ are Lagrange multipliers enforcing the
secondary constraints ${\cal S}^{(\mu)}$ and ${\cal S}^{(\nu)}$ respectively.
We then check that all constraints (primary and secondary) are preserved in
time by taking their Poisson bracket with $H_{\rm sec}$. Since the Poisson
bracket of all constraints with $H_{\rm pri}$ is weakly vanishing by default,
it is sufficient to consider the Poisson brackets of ${\cal S}^{(\mu)}$ and
${\cal S}^{(\nu)}$ with any constraint in the set $\\{{\cal H}_{i},{\cal
H},\Pi^{(\mu)},\Pi^{(\nu)},{\cal S}^{(\mu)},{\cal S}^{(\nu)}\\}$, implying
eleven brackets in total.

We set the index $A=\\{\mu,\nu\\}$ and consider collectively the vector
${{\cal S}^{(A)}}[{u^{(A)}}]$. From equation (50 we have:

$\displaystyle\left\\{{{\cal S}^{(A)}}[{u^{(A)}}],{\vec{{\cal
H}}}[{\vec{N}}]\right\\}$ $\displaystyle=-{{\cal S}^{(A)}}[{{\cal
L}_{\vec{N}}u^{(A)}}]$ (65)

which therefore vanish weakly. Hence, $\vec{{\cal H}}$ remains a first class
constraint. The remaining nine brackets do not vanish, meaning that all other
constraints are second class. Consider first

$\displaystyle\partial_{t}{\Pi^{(A)}}[{\lambda^{(A)}}]\approx\sum_{B}\left\\{{\Pi^{(A)}}[{\lambda^{(A)}}],{{\cal
S}^{(B)}}[{u^{(B)}}]\right\\}\approx 0.$ (66)

The involved Poisson brackets are evaluated as

$\displaystyle\left\\{{\Pi^{(A)}}[{\lambda^{(A)}}],{{\cal
S}^{(B)}}[{u^{(B)}}]\right\\}=-{C^{AB}}[{\lambda^{(A)}u^{(B)}}]$ (67)

where

$\displaystyle C^{AB}\equiv\frac{\partial{\cal S}^{(B)}}{\partial
A}=\begin{pmatrix}\frac{\partial{\cal
S}^{(\mu)}}{\partial\mu}&\frac{\partial{\cal S}^{(\mu)}}{\partial\nu}\\\
\frac{\partial{\cal S}^{(\nu)}}{\partial\mu}&\frac{\partial{\cal
S}^{(\nu)}}{\partial\nu}\end{pmatrix}$ (68)

Then (67) gives two homogeneous equations for two unknowns, the Lagrange
multipliers $u^{(\mu)}$ and $u^{(\nu)}$, and therefore implies that they must
both vanish.

Consider now

$\displaystyle\partial_{t}{{\cal H}}[{N}]\approx\sum_{B}\left\\{{{\cal
H}}[{N}],{{\cal S}^{(B)}}[{u^{(B)}}]\right\\}\approx 0.$ (69)

The involved Poisson brackets are evaluated as

${U^{B}}[{\xi^{(B)}}]\equiv\left\\{{{\cal H}}[{N}],{{\cal
S}^{(B)}}[{\xi^{(B)}}]\right\\}=\left\\{{{\cal H}}[{N}],\frac{\partial{{\cal
H}}[{\xi^{(B)}}]}{\partial B}\right\\}$ (70)

for arbitrary functions $\xi^{(A)}$ (not necessarily the Lagrange multipliers
$u^{(A)}$), which then allows us to use the ${\cal A}$ and ${\cal B}$
coefficients defined in (53) and (54) and displayed in appendix B. Explicitly,
after some integrations by parts and further computations we find

$\displaystyle U^{A}$ $\displaystyle=N\bigg{(}-{\cal
B}_{(\Pi^{ij})}\frac{\partial{\cal A}_{(q_{ij})}}{\partial A}+{\cal
A}_{(A_{i})}\frac{\partial{\cal B}_{(\Pi^{i})}}{\partial A}-{\cal
B}_{(\Pi^{i})}\frac{\partial{\cal A}_{(A_{i})}}{\partial A}$
$\displaystyle+{\cal A}_{(\phi)}\frac{\partial{\cal B}_{(\Pi)}}{\partial
A}+D_{i}{\cal B}_{(\Pi)}\frac{\partial{\cal A}^{i}_{(\phi)}}{\partial
A}\bigg{)}.$ (71)

It is crucial to note that $U^{(A)}$ can vanish neither strongly nor weakly.
This can be easily seen by focussing on the ${\cal
B}_{(\Pi^{ij})}\frac{\partial{\cal A}_{(q_{ij})}}{\partial A}$ term, which is
the only term containing $\Pi^{ij}$ and thus cannot be cancelled by other
terms. Back to (70), we set $\xi^{(B)}\rightarrow u^{(B)}\approx 0$ which
results in $\partial_{t}{{\cal
H}}[{N}]={U^{(\mu)}}[{u^{(\mu)}}]+{U^{(\nu)}}[{u^{(\nu)}}]\approx 0$ without
requiring any additional constraints or equations involving the Lagrange
multipliers.

Lastly, consider $\partial_{t}{{\cal S}^{(A)}}[{\xi^{(A)}}]$. Using our
definitions from (67) and (70), we have

$\displaystyle\partial_{t}{{\cal S}^{(A)}}[{\xi^{(A)}}]=$
$\displaystyle-{U^{A}}[{\xi^{(A)}}]-\sum_{B}{C^{AB}}[{\xi^{(A)}\lambda^{(B)}}]$
$\displaystyle+\sum_{B}\left\\{{{\cal S}^{(A)}}[{\xi^{(A)}}],{{\cal
S}^{(B)}}[{u^{(B)}}]\right\\}$ (72)

Now, the last Poisson bracket depends linearly on $u^{(B)}$ (and its
derivative) and the same for $\xi^{(A)}$, hence, it vanishes when
$u^{(B)}\approx 0$. Require then that $\partial_{t}{{\cal
S}^{(A)}}[{\xi^{(A)}}]\approx 0$ leads to the linear system of equations

$\displaystyle U^{A}+\sum_{B}C^{AB}\lambda^{(B)}\approx 0$ (73)

If the matrix $C^{AB}$ is invertible then the two equations above _determine_
the Lagrange multipliers $\lambda^{A}$ which then become functions of all the
phase-space variables. Then the Hamiltonian analysis is complete and no
further constraints in phase-space are required, with the conclusion that the
theory possesses three first class primary constraints ${\cal H}_{i}$, three
second class primary constraints ${\cal H}$ and $\Pi^{(A)}$, and two second
class secondary constraints ${\cal S}^{(A)}$. We list all the constraints and
their Poisson brackets in table 1. On the other hand, if the matrix has a
vanishing determinant, there exists a left null eigenvector which when applied
to the consistency relation may produce further constraints, called _tertiary_
constraints. Given the forms of (62),(63) and equations of motion (60),(61) it
is to be expected that $C^{AB}$ is indeed invertible in general situations.

###  V.5 Hamiltonian evolution

Having found all the constraints and solved for the Lagrange multipliers, we
may now form the _first class_ Hamiltonian $H_{\rm FC}$. We find this as the
secondary Hamiltonian with the Lagrange multiplies subbed-in. Given that
$u^{(A)}\approx 0$, we have that $H_{\rm sec}\approx H_{\rm pri}$, hence,

$H_{\rm FC}=\int d^{3}x\left[N{\cal H}+N^{i}{\cal
H}_{i}+\lambda^{(\mu)}\Pi^{(\mu)}+\lambda^{(\nu)}\Pi^{(\nu)}\right]$ (74)

with $\lambda^{(\mu)}$ and $\lambda^{(\nu)}$ being functions of the phase-
space variables $\\{q_{ij},A_{i},\phi,\Pi^{ij},\Pi^{i},\Pi,\mu,\nu\\}$, that
is, $\Pi^{(\mu)}$ and $\Pi^{(\nu)}$ are absent from $\lambda^{(A)}$.
Hamilton’s equations are then found using (45) with $H\rightarrow H_{\rm FC}$.
We find

$\displaystyle\dot{q}_{ij}\approx$ $\displaystyle N{\cal
B}_{(\Pi^{ij})}+2D_{(i}N_{j)}$ (75a) $\displaystyle\dot{\Pi}^{ij}\approx$
$\displaystyle-N{\cal
A}_{(q_{ij})}-2\Pi^{k(j}D_{k}N^{i)}+D_{k}\left(N^{k}\Pi^{ij}\right)$
$\displaystyle+\frac{\sqrt{q}}{16\pi\tilde{G}}\left(D^{i}D^{j}N-q^{ij}\vec{D}^{2}N\right)$
(75b) $\displaystyle\dot{A}_{i}\approx$ $\displaystyle N{\cal
B}_{(\Pi^{i})}-\chi D_{i}N+N^{j}D_{j}A_{i}+A_{j}D_{i}N^{j}$ (75c)
$\displaystyle\dot{\Pi}^{i}\approx$ $\displaystyle-N{\cal A}_{(A_{i})}-{\cal
A}^{j}_{(A_{i})}D_{j}N+D_{k}\left(N^{k}\Pi^{i}\right)-\Pi^{j}D_{j}N^{i}$ (75d)
$\displaystyle\dot{\phi}\approx$ $\displaystyle N{\cal
B}_{(\Pi)}+N^{i}D_{i}\phi$ (75e) $\displaystyle\dot{\Pi}\approx$
$\displaystyle D_{i}\left(\Pi N^{i}-N{\cal A}^{i}_{(\phi)}\right)$ (75f)
$\displaystyle\dot{\mu}\approx$ $\displaystyle\lambda^{(\mu)}+N^{i}D_{i}\mu$
(75g) $\displaystyle\dot{\nu}\approx$
$\displaystyle\lambda^{(\nu)}+N^{i}D_{i}\nu$ (75h)

along with the constraints (41), (60) and (61), and where to reiterate,
$\lambda^{(A)}$ are functions of all the phase-space variables obtained after
solving (73).

##  VI An example: Small perturbations around Minkowski space

The Hamiltonian density is rather complicated in general so for illustration
it is useful to have a look at a simple example. We consider a background
Minkowski space and the linear evolution of small perturbations to this
background.

###  VI.1 Background Minkowski solution

Before discussing the choice of function ${\cal U}(\nu,\mu)$, let us first
determine what conditions a Minkowski background imposes on the fields, in the
general sense. Note that a cosmological constant $\Lambda$ can always be
absorbed into the definition of ${\cal U}(\nu,\mu)$, hence, we set $\Lambda=0$
without any loss of generality. We use an overbar to denote the values that
fields take in the background and we fix the background gauge to

$\displaystyle\bar{N}=1\qquad\bar{N}^{i}=0.$ (76)

Moreover, for Minkowski space we have $\bar{q}_{ij}=\delta_{ij}$ and then
(75a) trivially gives $\bar{\Pi}^{ij}=0$ (and so $\bar{{\cal
B}}_{(\Pi^{ij})}=0$). Now, the background solution should not violate the
isometries of Minkowski space, including rotational invariance, hence, all
three dimensional vector fields should vanish:

$\displaystyle\bar{A}_{i}=\bar{\Pi}^{i}=D_{i}{\bar{\phi}}=0$ (77)

so that (75c) and (75d) are trivially satisfied (and $\bar{{\cal
B}}_{(\Pi^{i})}=\bar{{\cal A}}_{(A_{i})}=\bar{{\cal A}}^{i}_{(A_{j})}=0$,
while $\bar{{\cal B}}^{i}_{(\Pi^{j})}=-\delta^{i}_{\;\;j}$). The conditions
(77) in turn imply that $\bar{\chi}=1$ and ${\bar{\phi}}={\bar{\phi}}(t)$
leading in addition to $\bar{{\cal Y}}=0$. Hence, from (32) we find
$\bar{\Xi}={\bar{\nu}}$ and from (38) we get $\bar{C}_{1}=\bar{\Pi}$.

Meanwhile, (77) leads to $\bar{{\cal A}}^{i}_{(\phi)}=0$ (and hence
$\bar{{\cal A}}_{(\phi)}=0$) so that from (71) we find $\bar{U}^{A}=0$ and
hence, (73) leads to $\lambda^{(A)}=0$. This last condition then implies
through (75g) and (75h) that ${\bar{\mu}}$ and ${\bar{\nu}}$ are functions of
$\vec{x}$ only. Moreover, (75f) implies that $\bar{\Pi}$ is also a function of
$\vec{x}$ only.

Now from (75b) we have that $\bar{{\cal A}}_{(q_{ij})}=0$ and this imposes
that $\bar{{\cal U}}=(8\pi\tilde{G})^{2}\bar{\Pi}^{2}/{\bar{\nu}}$. Finally,
we have that $\bar{{\cal B}}_{(\Pi)}=8\pi\tilde{G}\bar{\Pi}/{\bar{\nu}}$ is
time-independent, so that (75e) may be integrated to get ${\bar{\phi}}={\cal
Q}_{0}t$, where ${\cal Q}_{0}=8\pi\tilde{G}\bar{\Pi}/{\bar{\nu}}$ is a
constant. We thus set

$\displaystyle 8\pi\tilde{G}\bar{\Pi}={\cal Q}_{0}{\bar{\nu}}$ (78)

without loss of generality, implying that $\bar{\Pi}$ and ${\bar{\nu}}$ have
the same $\vec{x}$ dependence, and that ${\cal Q}_{0}$, being a constant, is
independent of ${\bar{\nu}}$.

Now (62) gives $\bar{{\cal Q}}={\cal Q}_{0}$ so that the constraint (61) gives
$\partial{\cal U}/\partial{\bar{\nu}}={\cal Q}_{0}^{2}$. Since we also have
that that $\partial{\cal U}/\partial{\bar{\mu}}=0$ from (61), we find that

$\displaystyle\bar{{\cal U}}={\cal Q}_{0}^{2}{\bar{\nu}}$ (79)

which is consistent with our discussion in the previous paragraph. Note that
the above condition, is a condition that ${\cal U}$ must satisfy in order to
have Minkowski solutions but it does not completely determine the general form
of the function ${\cal U}$.

Finally, before discussing function choices, let us connect with the original
function ${\cal F}$. From (15) we find the relations

$\displaystyle\mu=\frac{\partial{\cal F}}{\partial{\cal Y}}$ (80)
$\displaystyle\nu=-\frac{\partial{\cal F}}{\partial{\cal Q}^{2}}$ (81)

Hence, since $\bar{{\cal Y}}=0$ and $\bar{{\cal Q}}={\cal Q}_{0}$ is a
constant, this implies that ${\cal F}$ and its derivatives will also be, at
best, constants. We immediately get that both ${\bar{\mu}}$ and ${\bar{\nu}}$
are constants and hence, so is $\bar{\Pi}$. Incidentally, (15) returns
$\bar{{\cal F}}=0$.

###  VI.2 Choice of function

We now turn to the choice of function ${\cal U}(\nu,\mu)$ in order to pave the
way for departures from Minkowski. Our goal here is to compare to Skordis and
Zlosnik (2022), hence, we restrict ourselves to cases where the function
${\cal F}({\cal Y},{\cal Q})$ of (1) takes the form

$\displaystyle{\cal F}({\cal Y},{\cal Q})$
$\displaystyle=(2-K_{{B}})\lambda_{s}{\cal Y}-2{\cal K}_{2}\left({\cal
Q}-{\cal Q}_{0}\right)^{2}$ (82)

where $\lambda_{s}$ and ${\cal K}_{2}$ are constants. The above functional
form is motivated by making sure that at the weak field quasistatic limit
Newtonian gravity is recovered and that the large-scale FLRW cosmology admits
dust solutions. Exact Minkowski is recovered when ${\cal Y}=0$ and ${\cal
Q}={\cal Q}_{0}$ in accordance with the previous subsection.

We then find that

$\displaystyle\mu=(2-K_{{B}})\lambda_{s}$ (83)

and

$\displaystyle\nu=2{\cal K}_{2}\left(1-\frac{{\cal Q}_{0}}{{\cal Q}}\right)$
(84)

which inverts to

$\displaystyle{\cal Q}=\frac{{\cal Q}_{0}}{1-\frac{\nu}{2{\cal K}_{2}}}$ (85)

so that

$\displaystyle{\cal U}=$ $\displaystyle\frac{{\cal
Q}_{0}^{2}\nu}{1-\frac{\nu}{2{\cal K}_{2}}}$ (86)

Using the above functional form into the constraint (60) leads to ${\cal
Y}=0$, while constraint (61) returns back (85). The last relation also leads
to

$\displaystyle{\bar{\nu}}=0,$ (87)

and hence,

$\displaystyle\bar{\Pi}=0.$ (88)

###  VI.3 Hamiltonian density to 2nd order

We perturb the lapse function as

$\displaystyle N=$ $\displaystyle 1+\Psi$ (89)

where $\Psi$ is a small perturbation and let the shift $N_{i}=h_{i}$ be a pure
perturbation (zero background). Additionally, we perturb our phase-space
variables as

$\displaystyle q_{ij}=$ $\displaystyle\delta_{ij}+h_{ij}$ (90)
$\displaystyle\Pi^{ij}=$ $\displaystyle\varpi^{ij}$ (91) $\displaystyle
A_{i}=$ $\displaystyle\alpha_{i}$ (92) $\displaystyle\Pi^{i}=$
$\displaystyle\varpi^{i}$ (93) $\displaystyle\phi=$ $\displaystyle{\cal
Q}_{0}t+\varphi$ (94) $\displaystyle\Pi=$ $\displaystyle\frac{{\cal
Q}_{0}}{8\pi\tilde{G}}{\bar{\nu}}+\varpi$ (95) $\displaystyle\nu=$
$\displaystyle{\bar{\nu}}+{\delta\nu}$ (96)

where $h_{ij}$, $\alpha_{i}$, $\varphi$, $\varpi^{ij}$, $\varpi^{i}$, $\varpi$
and ${\delta\nu}$ are small perturbations, and we further define $h\equiv
h_{i}^{\;\;i}$. We raise and lower indices using the background metric
$\delta_{ij}$.

The variable $\mu$ is fixed to (83) as it is a constant. We keep ${\bar{\nu}}$
in all calculations where the background otherwise vanishes and take the limit
${\bar{\nu}}\rightarrow 0$ only at the end. We first calculate the secondary
constraints, the first of which leads to ${\cal S}^{(\mu)}=0$ for the chosen
function. To expand the constraint (61), that is ${\cal S}^{(\nu)}\approx 0$,
to first order and determine ${\delta\nu}$ we need ${\cal Q}$ and ${\cal U}$
to first order, however, since the latter is also needed for the Hamiltonian
constraint to second order, we compute that to get

$\displaystyle{\cal Q}=$ $\displaystyle{\cal
Q}_{0}\left(1-\frac{1}{2}h-\frac{{\delta\nu}}{{\bar{\nu}}}+8\pi\tilde{G}\frac{\varpi}{{\cal
Q}_{0}{\bar{\nu}}}\right),$ (97) $\displaystyle{\cal U}=$ $\displaystyle{\cal
Q}_{0}^{2}\left[{\bar{\nu}}+{\delta\nu}+\frac{{(\delta\nu)^{2}}}{2{\cal
K}_{2}}\right].$ (98)

Thus, into (61) we find

$\displaystyle{\delta\nu}=\frac{1}{1+\frac{{\bar{\nu}}}{2{\cal
K}_{2}}}\left(8\pi\tilde{G}\frac{\varpi}{{\cal
Q}_{0}}-\frac{1}{2}{\bar{\nu}}h\right).$ (99)

We now compute the diffeomorphism constraint (36) which very easily leads to

$\displaystyle N^{i}{\cal H}_{i}\approx$ $\displaystyle\frac{{\cal
Q}_{0}{\bar{\nu}}}{8\pi\tilde{G}}h^{i}\vec{\nabla}_{i}\varphi-2h^{i}\vec{\nabla}_{j}\varpi^{j}_{\;\;i}$
(100)

Finally, we compute the Hamiltonian constraint (37). Since
$\vec{A}\cdot\vec{D}\phi$ is second order, we need to zeroth order
$C_{2}=(2-K_{{B}})(2+K_{{B}}\lambda_{s})/K_{{B}}$ and $C_{1}$ and $\Xi$ to
second order. These are calculated as

$\displaystyle C_{1}=$ $\displaystyle\frac{{\cal
Q}_{0}}{8\pi\tilde{G}}{\bar{\nu}}+\varpi-\frac{2-K_{{B}}}{K_{{B}}}\alpha^{i}\varpi_{i}$
(101) $\displaystyle\Xi=$
$\displaystyle{\bar{\nu}}+{\delta\nu}+\left[{\bar{\nu}}-\frac{(2-K_{{B}})(2+K_{{B}}\lambda_{s})}{K_{{B}}}\right]\alpha_{i}\alpha^{i}$
(102)

and into (37) we find

$\displaystyle N{\cal H}=$ $\displaystyle\frac{{\cal
Q}_{0}^{2}}{8\pi\tilde{G}}{\bar{\nu}}+{\cal Q}_{0}\varpi+\frac{{\cal
Q}_{0}^{2}}{8\pi\tilde{G}}{\bar{\nu}}\Psi+8\pi\tilde{G}\left[2\varpi^{ij}\varpi_{ij}-\hat{\varpi}^{2}+\frac{1}{2K_{{B}}}|\vec{\varpi}|^{2}+\frac{1}{4{\cal
K}_{2}}\frac{1}{1+\frac{{\bar{\nu}}}{2{\cal
K}_{2}}}\varpi^{2}\right]-\frac{2-K_{{B}}}{K_{{B}}}{\cal
Q}_{0}\alpha^{j}\varpi_{j}$ $\displaystyle-\frac{{\cal
Q}_{0}{\bar{\nu}}}{4{\cal K}_{2}}\frac{1}{1+\frac{{\bar{\nu}}}{2{\cal
K}_{2}}}\varpi h-\frac{{\cal
Q}_{0}^{2}{\bar{\nu}}}{16\pi\tilde{G}}\left[\frac{1}{8}\frac{1-\frac{{\bar{\nu}}}{2{\cal
K}_{2}}}{1+\frac{{\bar{\nu}}}{2{\cal
K}_{2}}}h^{2}+\frac{1}{4}h^{ij}h_{ij}+|\vec{\alpha}|^{2}\right]+{\cal
Q}_{0}\Psi\varpi-\frac{{\cal
Q}_{0}{\bar{\nu}}}{8\pi\tilde{G}}\vec{\alpha}\cdot\vec{\nabla}\phi+\Psi\vec{\nabla}_{i}\varpi^{i}-\frac{2-K_{{B}}}{K_{{B}}}\vec{\varpi}\cdot\vec{\nabla}\varphi$
$\displaystyle+\frac{1}{16\pi\tilde{G}}\bigg{[}-\frac{1}{4}|\vec{\nabla}h|^{2}+\frac{1}{2}\vec{\nabla}_{i}h\vec{\nabla}_{j}h^{ij}-\frac{1}{4}\vec{\nabla}_{k}h^{ij}\left(\vec{\nabla}_{i}h^{k}_{\;\;j}+\vec{\nabla}_{j}h^{k}_{\;\;i}-\vec{\nabla}^{k}h_{ij}\right)+\vec{\nabla}_{i}\Psi\vec{\nabla}_{j}h^{ij}-\vec{\nabla}\Psi\cdot\vec{\nabla}h\bigg{]}$
$\displaystyle+\frac{1}{16\pi\tilde{G}}\bigg{[}K_{{B}}|\vec{B}|^{2}+(2-K_{{B}})\frac{2+\lambda_{s}K_{{B}}}{K_{{B}}}|\vec{\nabla}\varphi+{\cal
Q}_{0}\vec{\alpha}|^{2}\bigg{]}.$ (103)

We may now set ${\bar{\nu}}=0$ and combine the two to get the second order
Hamiltonian as

$\displaystyle H^{(2)}=$ $\displaystyle\int
d^{3}x\bigg{\\{}8\pi\tilde{G}\left[2\varpi^{ij}\varpi_{ij}-\hat{\varpi}^{2}+\frac{1}{2K_{{B}}}|\vec{\varpi}|^{2}+\frac{1}{4{\cal
K}_{2}}\varpi^{2}\right]-\frac{2-K_{{B}}}{K_{{B}}}\vec{\varpi}\cdot\left(\vec{\nabla}\varphi+{\cal
Q}_{0}\vec{\alpha}\right)+\frac{1}{16\pi\tilde{G}}\bigg{[}-\frac{1}{4}|\vec{\nabla}h|^{2}$
$\displaystyle+\frac{1}{2}\vec{\nabla}_{i}h\vec{\nabla}_{j}h^{ij}-\frac{1}{4}\vec{\nabla}_{k}h^{ij}\left(\vec{\nabla}_{i}h^{k}_{\;\;j}+\vec{\nabla}_{j}h^{k}_{\;\;i}-\vec{\nabla}^{k}h_{ij}\right)+K_{{B}}|\vec{B}|^{2}+(2-K_{{B}})\frac{2+\lambda_{s}K_{{B}}}{K_{{B}}}|\vec{\nabla}\varphi+{\cal
Q}_{0}\vec{\alpha}|^{2}\bigg{]}$ $\displaystyle+\Psi\left[{\cal
Q}_{0}\varpi+\vec{\nabla}\cdot\vec{\varpi}+\frac{1}{16\pi\tilde{G}}\left(\vec{\nabla}^{2}h-\vec{\nabla}_{i}\vec{\nabla}_{j}h^{ij}\right)\right]-2h^{i}\vec{\nabla}_{j}\varpi^{j}_{\;\;i}\bigg{\\}}$
(104)

Having found the second order Hamiltonian, we compare with the results from
Skordis and Zlosnik (2022), which only find the scalar mode Hamiltonian. The
Lagrangian and corresponding Hamiltonian in Skordis and Zlosnik (2022) is
multiplied by $16\pi\tilde{G}$ compared with what we have in this work, and we
should do the same here for proper comparison. Doing that leads to a rescaling
of canonical momenta by $16\pi\tilde{G}$, that is, we define
$\tilde{\varpi}^{ij}\equiv 16\pi\tilde{G}\varpi^{ij}$ and similarly for
$\varpi^{i}$ and $\varpi$.

We define the traceless operator
$D_{ij}\equiv\vec{\nabla}_{i}\vec{\nabla}_{j}-\frac{1}{3}\vec{\nabla}^{2}\delta_{ij}$
and expand the perturbations in scalar modes as

$\displaystyle h_{ij}=$ $\displaystyle-2\Phi\delta_{ij}+D_{ij}\eta$ (105)
$\displaystyle\tilde{\varpi}^{ij}=$
$\displaystyle-\frac{1}{6}P_{\Phi}\delta^{ij}+\frac{3}{2\vec{\nabla}^{4}}D^{ij}P_{\eta}$
(106) $\displaystyle A_{i}=$ $\displaystyle\vec{\nabla}_{i}\alpha$ (107)
$\displaystyle\tilde{\varpi}^{i}=$
$\displaystyle-\vec{\nabla}^{i}\frac{1}{\vec{\nabla}^{2}}P_{\alpha}$ (108)
$\displaystyle h_{i}=$ $\displaystyle-\vec{\nabla}_{i}\zeta$ (109)

so that $\\{\Phi,P_{\Phi}\\}$, $\\{\eta,P_{\eta}\\}$,
$\\{\alpha,P_{\alpha}\\}$ and $\\{\varphi,P_{\varphi}\\}$ form canonical
pairs, that is, $\int
dtd^{3}x\left[\dot{h}_{ij}\varpi^{ij}+\dot{\alpha}_{i}\varpi^{i}\right]=\int
dtd^{3}x\left[\dot{\Phi}P_{\Phi}+\dot{\eta}P_{\eta}+\dot{\alpha}P_{\alpha}\right]$.
Using the above expressions into (104) we find the scalar-mode Hamiltonian as

$\displaystyle H^{(2)}=$ $\displaystyle\int
d^{3}x\bigg{\\{}-\frac{1}{24}P_{\Phi}^{2}+\frac{3}{2}\left|\frac{1}{\vec{\nabla}^{2}}P_{\eta}\right|^{2}+\frac{1}{4K_{{B}}}|\vec{\nabla}\frac{1}{\vec{\nabla}^{2}}P_{\alpha}|^{2}$
$\displaystyle-\frac{2-K_{{B}}}{K_{{B}}}P_{\alpha}\left(\varphi+{\cal
Q}_{0}\alpha\right)-2\left|\vec{\nabla}\left(\Phi+\frac{1}{6}\vec{\nabla}^{2}\eta\right)\right|^{2}$
$\displaystyle+\frac{1}{8{\cal
K}_{2}}P_{\varphi}^{2}+(2-K_{{B}})\frac{2+\lambda_{s}K_{{B}}}{K_{{B}}}|\vec{\nabla}\left(\varphi+{\cal
Q}_{0}\alpha\right)|^{2}$ $\displaystyle+\Psi\left[{\cal
Q}_{0}P_{\varphi}-P_{\alpha}-4\vec{\nabla}^{2}\left(\Phi+\frac{1}{6}\vec{\nabla}^{2}\eta\right)\right]$
$\displaystyle+2\zeta\left(\frac{1}{6}\vec{\nabla}^{2}P_{\Phi}-P_{\eta}\right)\bigg{\\}}$
(110)

Our comparison is, however, not yet finished because in Skordis and Zlosnik
(2022) one of the variables used is not $\varphi$ but the combination
$\chi\equiv\varphi+{\cal Q}_{0}\alpha$. We thus perform a canonical
transformation to new canonical pairs $\\{\chi,P_{\chi}\\}$ and
$\\{{\tilde{\alpha}},P_{\tilde{\alpha}}\\}$ defined through
${\tilde{\alpha}}=\alpha$, $P_{\chi}=P_{\varphi}$ and
$P_{\tilde{\alpha}}=P_{\alpha}-{\cal Q}_{0}P_{\varphi}$ to get

$\displaystyle H^{(2)}=$ $\displaystyle\int
d^{3}x\bigg{\\{}-\frac{1}{24}\left|P_{\Phi}\right|^{2}+\frac{3}{2k^{4}}\left|P_{\eta}\right|^{2}+\frac{1}{8{\cal
K}_{2}}\left|P_{\chi}\right|^{2}+\frac{1}{4k^{2}K_{{B}}}\left|{\cal
Q}_{0}P_{\chi}+P_{\tilde{\alpha}}\right|^{2}-\frac{2-K_{{B}}}{2K_{{B}}}[\chi\left(P_{\tilde{\alpha}}^{*}+{\cal
Q}_{0}P_{\chi}^{*}\right)+c.c.]$
$\displaystyle-2k^{2}\left|\Phi-\frac{1}{6}k^{2}\eta\right|^{2}+(2-K_{{B}})\frac{2+\lambda_{s}K_{{B}}}{K_{{B}}}k^{2}|\chi|^{2}+\Psi
C_{\Psi}^{*}+\Psi^{*}C_{\Psi}+\zeta
C_{\zeta}^{*}+\zeta^{*}C_{\zeta}\bigg{\\}}$ (111)

where

$\displaystyle C_{\Psi}\equiv$
$\displaystyle-\frac{1}{2}P_{\tilde{\alpha}}+2k^{2}\left(\Phi-\frac{1}{6}k^{2}\eta\right)\approx
0$ (112) $\displaystyle C_{\zeta}\equiv$
$\displaystyle-\frac{1}{6}k^{2}P_{\Phi}-P_{\eta}\approx 0$ (113)

are two constraints imposed by the Lagrange multipliers $\Psi$ and $\zeta$.
The resulting Hamiltonian (111) is identical to the one found in Skordis and
Zlosnik (2022), up to some symbol relabeling.

##  VII Conclusions

In this article, we presented the general Hamiltonian analysis for AeST
theory, which extends GR with the inclusion of a unit-timelike vector field
$\hat{A}^{\mu}$ and scalar field $\phi$ in addition to the metric
$g_{\mu\nu}$. To simplify the computations we further introduced two auxiliary
fields $\mu$ and $\nu$ in order to avoid inversion of the free function ${\cal
F}$ which is part of the theory. Our analysis revealed the existence of four
first class constraints and four second class constraints. The first class
constraints consist of the three (primary) constraints ${\cal H}_{i}$ defined
in (36) and the first class Hamiltonian $H_{\rm FC}$ defined in (74), which is
the linear combination of the primary Hamiltonian constraint ${\cal H}$
defined in (37) (by itself second class), ${\cal H}_{i}$ and also
$\Pi^{(\mu)}$ and $\Pi^{(\nu)}$. The four second class constraints are the two
canonical momenta $\Pi^{(\mu)}$ and $\Pi^{(\nu)}$ of the auxiliary fields
$\mu$ and $\nu$, see (33) and (34), and the two secondary constraints ${\cal
S}^{(\mu)}$ and ${\cal S}^{(\nu)}$ defined through (60) and (61). The
existence of these second class constraints arises from the presence of the
auxiliary fields $\mu$ and $\nu$. See table 1 for a summary of these
cosntraints.

We may use the constraint analysis to count the number of physical degrees of
freedom. We have six variables in the spatial metric $q_{ij}$, three in
$A_{i}$ and one for each of $\phi$, $\mu$ and $\nu$, that is $12$ in total.
Counting in the canonical momenta doubles this to $24$. We subtract the four
second class constraints and twice the number of first class constraints which
remove the gauge redundant degrees of freedom, that is, we subtract $12$
degrees of freedom because of the constraints. We finally divide by two to
find six physical degrees of freedom.

Taking linear perturbations around a Minkowski background we expanded the
Hamiltonian to quadratic order and recovered the same results found in Skordis
and Zlosnik (2022) using different methods. In the process we showed that the
number of perturbative degrees of freedom found in Skordis and Zlosnik (2022)
matches the number found here using the full non-linear theory. Our formalism
may be used to compute the quadratic Hamiltonian of AeST theory on other
backgrounds in order to determine whether those backgrounds are stable or not.
Of particular interest are the case of de Sitter space and static spherically
symmetric configurations which we leave for a future work.

###### Acknowledgements.

TZ was supported by the the Czech Science Foundation (GAČR) grant No.
20-28525S until 30.09.2022., and from 01.11.2022. by the project No.
2021/43/P/ST2/02141 co-funded by the Polish National Science Centre and the
European Union Framework Programme for Research and Innovation Horizon 2020
under the Marie Skłodowska-Curie grant agreement No. 945339. CS was partially
supported 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).

##  Appendix A Useful Results

In the ADM formulation the vector field aligned with the time direction is

$\displaystyle t^{\mu}=(1,\vec{0})$ (114)

while the normal to the hypersurface is decomposed as

$\displaystyle n^{\mu}=$ $\displaystyle\frac{1}{N}(1,-\vec{N})\;,$
$\displaystyle n_{\mu}=$ $\displaystyle(-N,\vec{0})$ (115)

and the projector to the hypersurface as

$\displaystyle q^{0}_{\phantom{0}0}=$ $\displaystyle q^{0}_{\phantom{0}i}=0$
$\displaystyle q^{i}_{\phantom{i}0}=$ $\displaystyle N^{i}$ $\displaystyle
q^{i}_{\phantom{i}j}=\delta^{i}_{\phantom{i}j}$ (116)

The metric then has components

$\displaystyle g_{00}=$ $\displaystyle-N^{2}+|\vec{N}|^{2}\;,$ $\displaystyle
g^{00}=$ $\displaystyle-\frac{1}{N^{2}}\;,$ (117) $\displaystyle g_{0i}=$
$\displaystyle N_{i}\;,$ $\displaystyle g^{0i}=$
$\displaystyle\frac{1}{N^{2}}N^{i}\;,$ (118) $\displaystyle g_{ij}=$
$\displaystyle q_{ij}$ $\displaystyle g^{ij}=$ $\displaystyle
q^{ij}-\frac{N^{i}N^{j}}{N^{2}}$ (119)

The Christoffel connection splits into:

$\displaystyle\Gamma^{0}_{00}=$
$\displaystyle\frac{1}{N}\left[\dot{N}+N^{i}D_{i}N+N^{i}N^{j}K_{ij}\right]$
(120) $\displaystyle\Gamma^{0}_{0i}=$
$\displaystyle\frac{1}{N}\left[D_{i}N+N^{j}K_{ji}\right]$ (121)
$\displaystyle\Gamma^{0}_{ij}=$ $\displaystyle\frac{1}{N}K_{ij}$ (122)
$\displaystyle\Gamma^{i}_{00}=$
$\displaystyle\dot{N}^{i}+N^{j}D_{j}N^{i}+2NK^{ij}N_{j}+ND^{i}N$
$\displaystyle-\frac{N^{i}}{N}\left[\dot{N}+N^{k}N^{l}K_{kl}+N^{j}D_{j}N\right]$
(123) $\displaystyle\Gamma^{i}_{0j}=$ $\displaystyle
NK^{i}_{\phantom{i}j}-\frac{1}{N}N^{i}N^{k}K_{kj}+D_{j}N^{i}-\frac{N^{i}}{N}D_{j}N$
(124) $\displaystyle\Gamma^{k}_{ij}=$
$\displaystyle\gamma^{k}_{ij}-\frac{N^{k}}{N}K_{ij}$ (125)

The components of the extrinsic curvature $K_{\mu\nu}$ can then be calculated
to be:

$\displaystyle K_{00}=$ $\displaystyle N^{i}N^{j}K_{ij}$ (126) $\displaystyle
K_{0i}=K_{i0}=$ $\displaystyle N^{j}K_{ij}$ (127) $\displaystyle
K_{ij}=N\Gamma^{0}_{ij}=$
$\displaystyle\frac{1}{2N}\left(\dot{q}_{ij}-D_{i}N_{j}-D_{j}N_{i}\right)$
(128)

The spacetime-covariant derivative of $\hat{A}^{\mu}$ is given by:

$\displaystyle\nabla_{\nu}\hat{A}^{\mu}$
$\displaystyle=\partial_{\nu}\hat{A}^{\mu}+\Gamma^{\mu}_{\nu\sigma}\hat{A}^{\sigma}$
(129)

and individual components are

$\displaystyle\nabla_{0}\hat{A}^{0}=$
$\displaystyle\frac{\dot{\chi}}{N}+\frac{1}{N}A^{i}\left(D_{i}N+N^{j}K_{ij}\right)$
(130) $\displaystyle\nabla_{i}\hat{A}^{0}=$
$\displaystyle\frac{D_{i}\chi}{N}+\frac{1}{N}K_{ik}A^{k}$ (131)
$\displaystyle\nabla_{0}\hat{A}^{i}=$
$\displaystyle\dot{A}^{i}+\left[NK^{i}_{\phantom{i}j}+D_{j}N^{i}-\frac{N^{i}}{N}\left(N^{k}K_{kj}+D_{j}N\right)\right]A^{j}$
$\displaystyle-\frac{\dot{\chi}}{N}N^{i}+\chi\left(K^{ij}N_{j}+D^{i}N\right)$
(132) $\displaystyle\nabla_{i}\hat{A}^{j}=$ $\displaystyle D_{i}A^{j}+\chi
K^{j}_{\phantom{i}i}-\frac{N^{j}}{N}\left(D_{i}\chi+K_{ik}A^{k}\right)$ (133)

From the above relations we then construct $J^{\mu}$ as

$\displaystyle J^{0}=$
$\displaystyle\frac{\chi}{N^{2}}\left[\dot{\chi}+A^{i}D_{i}N-N^{i}D_{i}\chi\right]$
$\displaystyle+\frac{A^{i}}{N}\left[K_{ij}A^{j}+D_{i}\chi\right]$ (134)
$\displaystyle J^{i}=$
$\displaystyle\frac{\chi}{N}\bigg{\\{}\dot{A}^{i}+\left[NK^{i}_{\phantom{i}j}+D_{j}N^{i}-\frac{N^{i}}{N}D_{j}N\right]A^{j}$
$\displaystyle-\frac{\dot{\chi}}{N}N^{i}+\chi
D^{i}N-N^{j}\left[D_{j}A^{i}-\frac{N^{i}}{N}D_{j}\chi\right]\bigg{\\}}$
$\displaystyle+A^{j}\left[D_{j}A^{i}+\chi
K^{i}_{\phantom{i}j}-\frac{N^{i}}{N}\left(D_{j}\chi+K_{jk}A^{k}\right)\right]$
(135)

##  Appendix B Coefficients for the variations of the smeared Hamiltonian
constraint

In evaluating the smeared Hamiltonian constraint it is useful to define the
coefficients ${\cal A}_{(Q_{I})}$, ${\cal A}_{(Q_{I})}^{i}$, ${\cal
B}_{(P_{I})}$ and ${\cal B}_{(P_{I})}^{i}$ where
$Q_{I}=\\{q_{ij},A_{i},\phi\\}$ and $P_{I}=\\{\Pi^{ij},\Pi^{i},\Pi\\}$. See
(53) and (54). Using (37) and the variables $\Xi$, $C_{1}$, $C_{2}$ defined
through (32), (38) and (39) respectively the ${\cal A}$-coefficients are found
to be

$\displaystyle{\cal A}_{(q_{ij})}=$
$\displaystyle\frac{8\pi\tilde{G}}{\sqrt{q}}\left\\{4\Pi^{k(i}\Pi_{k}^{\;\;j)}-2\hat{\Pi}\Pi^{ij}+\frac{1}{2K_{{B}}}\Pi^{i}\Pi^{j}-\frac{C_{1}^{2}C_{2}}{2\Xi^{2}}A^{i}A^{j}-\frac{1}{2}\left[\Pi^{kl}\left(2\Pi_{kl}-\hat{\Pi}q_{kl}\right)+\frac{1}{2K_{{B}}}|\vec{\Pi}|^{2}+\frac{C_{1}^{2}}{2\Xi}\right]q^{ij}\right\\}$
$\displaystyle+\frac{1}{2\chi}\left(\frac{2-K_{{B}}}{K_{{B}}}\vec{\Pi}\cdot\vec{D}\phi-\frac{C_{1}C_{2}}{\Xi}\vec{A}\cdot\vec{D}\phi-\vec{D}\cdot\vec{\Pi}\right)A^{i}A^{j}-\frac{\chi
C_{1}C_{2}}{\Xi}\left(\frac{C_{2}}{\Xi}\vec{A}\cdot\vec{D}\phi
A^{i}A^{j}+A^{(i}D^{j)}\phi\right)$
$\displaystyle+\frac{\sqrt{q}}{32\pi\tilde{G}}\bigg{\\{}-{\cal R}+2\Lambda-
K_{{B}}|\vec{B}|^{2}+\left(\frac{C_{2}^{2}\chi^{2}}{\Xi}+2-K_{{B}}+\mu-\nu\right)\left(\vec{A}\cdot\vec{D}\phi\right)^{2}+2(2-K_{{B}})\vec{A}\times\vec{B}\cdot\vec{D}\phi$
$\displaystyle+\left(2-K_{{B}}+\mu+\frac{(2-K_{{B}})^{2}}{K_{{B}}}\chi^{2}\right)|\vec{D}\phi|^{2}+{\cal
U}\bigg{\\}}q^{ij}+\frac{\sqrt{q}}{16\pi\tilde{G}}\bigg{\\{}R^{ij}+K_{{B}}B^{i}B^{j}-\frac{(2-K_{{B}})^{2}}{K_{{B}}}|\vec{D}\phi|^{2}A^{i}A^{j}$
$\displaystyle-2\left(\frac{C_{2}^{2}\chi^{2}}{\Xi}+2-K_{{B}}+\mu-\nu\right)\vec{A}\cdot\vec{D}\phi
A^{(i}D^{j)}\phi+2(2-K_{{B}})B_{l}\left(A_{k}\epsilon^{lk(i}D^{j)}\phi-\epsilon^{lk(i}A^{j)}D_{k}\phi\right)$
$\displaystyle-\frac{C_{2}^{2}}{\Xi^{2}}\left(2\frac{2-K_{{B}}}{K_{{B}}}+\mu\right)\left(\vec{A}\cdot\vec{D}\phi\right)^{2}A^{i}A^{j}-\left(2-K_{{B}}+\mu+\frac{(2-K_{{B}})^{2}}{K_{{B}}}\chi^{2}\right)D^{i}\phi
D^{j}\phi\bigg{\\}},$ (136) $\displaystyle{\cal A}_{(A_{i})}=$
$\displaystyle\frac{8\pi\tilde{G}}{\sqrt{q}}\,\frac{C_{1}}{\Xi}\left[\frac{C_{1}\,C_{2}}{\Xi}A^{i}-\frac{2-K_{{B}}}{K_{{B}}}\Pi^{i}\right]-\frac{1}{\chi}\left\\{\frac{2-K_{{B}}}{K_{{B}}}\vec{\Pi}\cdot\vec{D}\phi-\frac{C_{1}C_{2}}{\Xi}\vec{A}\cdot\vec{D}\phi-
D_{j}\Pi^{j}\right\\}A^{i}$ $\displaystyle-\frac{\chi
C_{2}}{\Xi}\bigg{\\{}\frac{2-K_{{B}}}{K_{{B}}}\vec{A}\cdot\vec{D}\phi\Pi^{i}-C_{1}\left(D^{i}\phi+\frac{2C_{2}}{\Xi}\vec{A}\cdot\vec{D}\phi
A^{i}\right)\bigg{\\}}$
$\displaystyle+\frac{\sqrt{q}}{8\pi\tilde{G}}\bigg{\\{}\epsilon^{ijk}\left[K_{{B}}D_{j}B_{k}+(2-K_{{B}})B_{j}D_{k}\phi\right]+\left[\frac{C_{2}^{2}\chi^{2}}{\Xi}+2-K_{{B}}+\mu-\nu\right]\vec{A}\cdot\vec{D}\phi
D^{i}\phi$
$\displaystyle+\frac{C_{2}^{2}}{\Xi^{2}}\left[2\frac{2-K_{{B}}}{K_{{B}}}+\mu\right]\left(\vec{A}\cdot\vec{D}\phi\right)^{2}A^{i}+\frac{(2-K_{{B}})^{2}}{K_{{B}}}|\vec{D}\phi|^{2}A^{i}+(2-K_{{B}})D_{j}(A^{j}D^{i}\phi-A^{i}D^{j}\phi)\bigg{\\}},$
(137) $\displaystyle{\cal A}^{i}_{(A_{j})}=$
$\displaystyle\frac{\sqrt{q}}{8\pi\tilde{G}}\left[-K_{{B}}\epsilon^{ijk}B_{k}+(2-K_{{B}})\left(\vec{A}^{i}\vec{D}^{j}\phi-\vec{A}^{j}\vec{D}^{i}\phi\right)\right],$
(138) $\displaystyle{\cal A}^{i}_{(\phi)}=$
$\displaystyle\frac{2-K_{{B}}}{K_{{B}}}\chi\Pi^{i}-\frac{\chi}{\Xi}C_{1}C_{2}A^{i}-\frac{\sqrt{q}}{8\pi\tilde{G}}\bigg{\\{}(2-K_{{B}})\epsilon^{ijk}\vec{A}_{j}\vec{B}_{k}+\left[\frac{\chi^{2}}{\Xi}C_{2}^{2}+2-K_{{B}}+\mu-\nu\right]\left(\vec{A}\cdot\vec{D}\phi\right)A^{i}$
$\displaystyle+\left(2-K_{{B}}+\mu+\frac{(2-K_{{B}})^{2}}{K_{{B}}}\chi^{2}\right)D^{i}\phi\bigg{\\}}$
(139)

and

$\displaystyle{\cal A}^{i}_{(q_{jk})}=$ $\displaystyle 0$ (140)
$\displaystyle{\cal A}_{(\phi)}=$ $\displaystyle D_{i}{\cal A}^{i}_{(\phi)}.$
(141)

The ${\cal B}$-coefficients are

$\displaystyle{\cal B}_{(\Pi^{ij})}=$
$\displaystyle\frac{16\pi\tilde{G}}{\sqrt{q}}\left(2\Pi_{ij}-\hat{\Pi}q_{ij}\right)$
(142) $\displaystyle{\cal B}^{i}_{(\Pi^{jk})}=$ $\displaystyle 0$ (143)
$\displaystyle{\cal B}_{(\Pi^{i})}=$ $\displaystyle-
D_{i}\chi-\chi\frac{2-K_{{B}}}{K_{{B}}}\bigg{\\{}D_{i}\phi+\frac{1}{\Xi}C_{2}\,\vec{A}\cdot\vec{D}\phi\;A_{i}\bigg{\\}}$
$\displaystyle+\frac{8\pi\tilde{G}}{K_{{B}}\sqrt{q}}\bigg{\\{}\Pi_{i}-\frac{2-K_{{B}}}{\Xi}C_{1}A_{i}\bigg{\\}}$
(144) $\displaystyle{\cal B}^{i}_{(\Pi^{j})}=$
$\displaystyle-\chi\delta^{i}_{\;\;j}$ (145) $\displaystyle{\cal B}_{(\Pi)}=$
$\displaystyle\frac{1}{\Xi}\bigg{\\{}\frac{8\pi\tilde{G}}{\sqrt{q}}C_{1}+\chi
C_{2}\,\vec{A}\cdot\vec{D}\phi\bigg{\\}}$ (146) $\displaystyle{\cal
B}^{i}_{(\Pi)}=$ $\displaystyle 0$ (147)

## References

  * Will (2014) C. M. Will, Living Rev. Rel. 17, 4 (2014), arXiv:1403.7377 [gr-qc] .
  * Abbott  _et al._ (2017) B. Abbott _et al._ (Virgo, LIGO Scientific), Phys. Rev. Lett. 119, 161101 (2017), arXiv:arXiv:1710.05832 [gr-qc] .
  * Starobinsky (1980) A. A. Starobinsky, Phys. Lett. B 91, 99 (1980).
  * Rubin  _et al._ (1980) V. C. Rubin, N. Thonnard, and W. K. Ford, Jr., Astrophys. J. 238, 471 (1980).
  * Lelli  _et al._ (2019) F. Lelli, S. S. McGaugh, J. M. Schombert, H. Desmond, and H. Katz, Mon. Not. Roy. Astron. Soc. 484, 3267 (2019), arXiv:1901.05966 [astro-ph.GA] .
  * Hahn  _et al._ (2007) O. Hahn, C. Porciani, C. M. Carollo, and A. Dekel, Mon. Not. Roy. Astron. Soc. 375, 489 (2007), arXiv:astro-ph/0610280 .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), arXiv:1807.06209 [astro-ph.CO] .
  * Jungman  _et al._ (1996) G. Jungman, M. Kamionkowski, and K. Griest, Phys. Rept. 267, 195 (1996), arXiv:hep-ph/9506380 .
  * Bertone  _et al._ (2005) G. Bertone, D. Hooper, and J. Silk, Phys. Rept. 405, 279 (2005), arXiv:hep-ph/0404175 .
  * Marsh (2016) D. J. E. Marsh, Phys. Rept. 643, 1 (2016), arXiv:1510.07633 [astro-ph.CO] .
  * Green and Kavanagh (2021) A. M. Green and B. J. Kavanagh, J. Phys. G 48, 043001 (2021), arXiv:2007.10722 [astro-ph.CO] .
  * Milgrom (1983) M. Milgrom, Astrophys. J. 270, 365 (1983).
  * Bekenstein and Milgrom (1984) J. Bekenstein and M. Milgrom, Astrophys. J. 286, 7 (1984).
  * Milgrom (2010) M. Milgrom, Mon. Not. Roy. Astron. Soc. 403, 886 (2010), arXiv:0911.5464 [astro-ph.CO] .
  * Milgrom (2023) M. Milgrom, (2023), arXiv:2305.19986 [gr-qc] .
  * Famaey and McGaugh (2012) B. Famaey and S. McGaugh, Living Rev. Rel. 15, 10 (2012), arXiv:arXiv:1112.3960 [astro-ph.CO] .
  * Banik and Zhao (2022) I. Banik and H. Zhao, Symmetry 14, 1331 (2022), arXiv:2110.06936 [astro-ph.CO] .
  * Bekenstein (1988) J. D. Bekenstein, Phys. Lett. B202, 497 (1988).
  * Sanders (1997) R. H. Sanders, Astrophys. J. 480, 492 (1997), astro-ph/9612099 [astro-ph] .
  * Bekenstein (2004) J. D. Bekenstein, Phys. Rev. D70, 083509 (2004), [Erratum: Phys. Rev.D71,069901(2005)], arXiv:astro-ph/0403694 [astro-ph] .
  * Navarro and Van Acoleyen (2006) I. Navarro and K. Van Acoleyen, JCAP 0609, 006 (2006), arXiv:gr-qc/0512109 [gr-qc] .
  * Zlosnik  _et al._ (2007) T. G. Zlosnik, P. G. Ferreira, and G. D. Starkman, Phys. Rev. D75, 044017 (2007), arXiv:astro-ph/0607411 [astro-ph] .
  * Sanders (2005) R. Sanders, Mon. Not. Roy. Astron. Soc. 363, 459 (2005), arXiv:astro-ph/0502222 .
  * Milgrom (2009) M. Milgrom, Phys. Rev. D80, 123536 (2009), arXiv:arXiv:0912.0790 [gr-qc] .
  * Babichev  _et al._ (2011) E. Babichev, C. Deffayet, and G. Esposito-Farese, Phys. Rev. D84, 061502 (2011), arXiv:arXiv:1106.2538 [gr-qc] .
  * Deffayet  _et al._ (2011) C. Deffayet, G. Esposito-Farese, and R. P. Woodard, Phys. Rev. D84, 124054 (2011), arXiv:arXiv:1106.4984 [gr-qc] .
  * Blanchet and Marsat (2011) L. Blanchet and S. Marsat, Phys. Rev. D84, 044056 (2011), arXiv:arXiv:1107.5264 [gr-qc] .
  * Sanders (2011) R. H. Sanders, Phys. Rev. D 84, 084024 (2011).
  * Mendoza  _et al._ (2012) S. Mendoza, T. Bernal, J. C. Hidalgo, and S. Capozziello, AIP Conf. Proc. 1458, 483 (2012), arXiv:1202.3629 [gr-qc] .
  * Woodard (2015) R. P. Woodard, Can. J. Phys. 93, 242 (2015), arXiv:arXiv:1403.6763 [astro-ph.CO] .
  * Khoury (2015) J. Khoury, Phys. Rev. D91, 024022 (2015), arXiv:arXiv:1409.0012 [hep-th] .
  * Blanchet and Heisenberg (2015) L. Blanchet and L. Heisenberg, Phys. Rev. D 91, 103518 (2015), arXiv:1504.00870 [gr-qc] .
  * Hossenfelder (2017) S. Hossenfelder, Phys. Rev. D 95, 124018 (2017), arXiv:1703.01415 [gr-qc] .
  * Burrage  _et al._ (2019) C. Burrage, E. J. Copeland, C. Käding, and P. Millington, Phys. Rev. D 99, 043539 (2019), arXiv:1811.12301 [astro-ph.CO] .
  * Milgrom (2019) M. Milgrom, Phys. Rev. D 100, 084039 (2019), arXiv:1908.01691 [gr-qc] .
  * D’Ambrosio  _et al._ (2020) F. D’Ambrosio, M. Garg, and L. Heisenberg, (2020), arXiv:2004.00888 [gr-qc] .
  * Käding (2023) C. Käding, Astronomy 2, 128 (2023), arXiv:2304.05875 [astro-ph.CO] .
  * Blanchet (2007) L. Blanchet, Class. Quant. Grav. 24, 3529 (2007), arXiv:astro-ph/0605637 .
  * Blanchet and Le Tiec (2009) L. Blanchet and A. Le Tiec, Phys. Rev. D 80, 023524 (2009), arXiv:0901.3114 [astro-ph.CO] .
  * Berezhiani and Khoury (2016) L. Berezhiani and J. Khoury, Phys. Lett. B753, 639 (2016), arXiv:arXiv:1506.07877 [astro-ph.CO] .
  * Berezhiani and Khoury (2015) L. Berezhiani and J. Khoury, Phys. Rev. D92, 103510 (2015), arXiv:arXiv:1507.01019 [astro-ph.CO] .
  * Kaplinghat  _et al._ (2016) M. Kaplinghat, S. Tulin, and H.-B. Yu, Phys. Rev. Lett. 116, 041302 (2016), arXiv:1508.03339 [astro-ph.CO] .
  * Kamada  _et al._ (2017) A. Kamada, M. Kaplinghat, A. B. Pace, and H.-B. Yu, Phys. Rev. Lett. 119, 111102 (2017), arXiv:1611.02716 [astro-ph.GA] .
  * Savchenko  _et al._ (2017) V. Savchenko _et al._ , Astrophys. J. 848, L15 (2017), arXiv:arXiv:1710.05449 [astro-ph.HE] .
  * Skordis and Złośnik (2019) C. Skordis and T. Złośnik, Phys. Rev. D100, 104013 (2019), arXiv:1905.09465 [gr-qc] .
  * Skordis and Zlosnik (2021) C. Skordis and T. Zlosnik, Phys. Rev. Lett. 127, 161302 (2021), arXiv:2007.00082 [astro-ph.CO] .
  * Scherrer (2004) R. J. Scherrer, Phys. Rev. Lett. 93, 011301 (2004), arXiv:astro-ph/0402316 [astro-ph] .
  * Arkani-Hamed  _et al._ (2004) N. Arkani-Hamed, H.-C. Cheng, M. A. Luty, and S. Mukohyama, JHEP 05, 074 (2004), hep-th/0312099 [hep-th] .
  * Arkani-Hamed  _et al._ (2007) N. Arkani-Hamed, H.-C. Cheng, M. A. Luty, S. Mukohyama, and T. Wiseman, JHEP 01, 036 (2007), arXiv:hep-ph/0507120 .
  * Skordis and Zlosnik (2022) C. Skordis and T. Zlosnik, Phys. Rev. D 106, 104041 (2022), arXiv:2109.13287 [gr-qc] .
  * Gumrukcuoglu  _et al._ (2016) A. E. Gumrukcuoglu, S. Mukohyama, and T. P. Sotiriou, Phys. Rev. D 94, 064001 (2016), arXiv:1606.00618 [hep-th] .
  * Bernardo and Chen (2023) R. C. Bernardo and C.-Y. Chen, Gen. Rel. Grav. 55, 23 (2023), arXiv:2202.08460 [gr-qc] .
  * Kashfi and Roshan (2022) T. Kashfi and M. Roshan, JCAP 10, 029 (2022), arXiv:2204.05672 [gr-qc] .
  * Mistele (2022) T. Mistele, JCAP 11, 008 (2022), arXiv:2103.16954 [gr-qc] .
  * Mistele  _et al._ (2023) T. Mistele, S. McGaugh, and S. Hossenfelder,  (2023), 10.1051/0004-6361/202346025, arXiv:2301.03499 [astro-ph.GA] .
  * Tian  _et al._ (2023) S. Tian, S. Hou, S. Cao, and Z.-H. Zhu, Phys. Rev. D 107, 044062 (2023), arXiv:2302.13304 [gr-qc] .
  * Llinares (2023) C. Llinares, (2023), arXiv:2302.12032 [astro-ph.GA] .
  * Verwayen  _et al._ (2023) P. Verwayen, C. Skordis, and C. Bœhm, (2023), arXiv:2304.05134 [astro-ph.CO] .
  * Mistele (2023) T. Mistele, (2023), arXiv:2305.07742 [gr-qc] .
  * Alexandrov  _et al._ (2013) S. Alexandrov, K. Krasnov, and S. Speziale, JHEP 06, 068 (2013), arXiv:1212.3614 [hep-th] .
  * Alexandrov  _et al._ (2021) S. Alexandrov, S. Speziale, and T. Zlosnik, Class. Quant. Grav. 38, 175011 (2021), arXiv:2104.03753 [gr-qc] .
  * Romano (1993) J. D. Romano, Gen. Rel. Grav. 25, 759 (1993), arXiv:gr-qc/9303032 [gr-qc] .