# Higher-dimensional minimal theory of mass-varying massive gravity and its
cosmological consequences

Ahmad Khoirul Falah<EMAIL_ADDRESS>Theoretical Physics
Laboratory, Theoretical High Energy Physics Research Division, Faculty of
Mathematics and Natural Sciences, Institut Teknologi Bandung, Jl. Ganesha no.
10 Bandung, Indonesia, 40132  Andy Octavian Latief<EMAIL_ADDRESS>Physics
of Magnetism and Photonics Research Division, Faculty of Mathematics and
Natural Sciences, Institut Teknologi Bandung, Jl. Ganesha no. 10 Bandung,
Indonesia, 40132  Husin Alatas<EMAIL_ADDRESS>Theoretical Physics
Division, Department of Physics, IPB University (Bogor Agricultural
University), Jl. Meranti, Kampus IPB Darmaga, Bogor 16680, Indonesia  Bobby
Eka Gunara<EMAIL_ADDRESS>(Corresponding author) Theoretical Physics
Laboratory, Theoretical High Energy Physics Research Division, Faculty of
Mathematics and Natural Sciences, Institut Teknologi Bandung, Jl. Ganesha no.
10 Bandung, Indonesia, 40132

###### Abstract

In this paper we construct higher-dimensional minimal theory of mass-varying
massive gravity (MTMVMG) where the masslike scalar potential is coupled to a
vielbein potential, unlike in the previous literature where it is coupled to
metric, such that the number of graviton degrees of freedom in the theory is
the same as in general relativity. We then study the cosmological aspects of
this theory and show that it has eight critical points: five in the massless
sector and three in the massive sector. In contrast to the standard theory of
mass-varying massive gravity where the graviton mass asymptotically approaches
zero at late times, hence making the contribution of massive gravity to the
late-time cosmic expansion minimum, the MTMVMG can provide good descriptions
both in the massless and massive sectors. Especially, there are at least two
interesting possible scenarios for the late-time cosmology in the theory: the
dark energy is either due to the constant graviton mass which comes from the
scalar field that becomes frozen after the reheating era, or due to the
quintessence paradigm where the scalar field is dynamic. Therefore, if the
accelerating expansion of the universe in the massless sector can be explained
by standard quintessence paradigm, in the massive sector it has to be
explained by the nontrivial interplay between quintessence and massive
gravity.

##  I Introduction

Several alternative theories modifying Einstein’s general relativity have been
proposed in the past decades as parts of an effort to solve long-standing
cosmological problems such as dark energy, dark matter, and inflation. One of
these theories, based on the assumption that graviton might have nonzero mass
and hence later coined as the theory of massive gravity, was originally
visioned by M. Fierz and W. Pauli in 1939 [1]. However, the theory did not
have a continuous transition to general relativity in the limit of zero
graviton mass, an issue known as the van Dam–Veltman–Zakharov discontinuity
[2, 3]. This problem was remedied by Vainshtein’s nonlinear mechanism in 1972
[4], but the nonlinear terms then gave rise to another problem called the
Boulware-Deser ghost [5]. Building on several previous attempts [6, 7], this
problem was then finally resolved in 2010 by C. de Rham, G. Gabadadze, and A.
Tolley (dRGT) [8, 9], resulting in a Lorentz-invariant, ghost-free nonlinear
theory of massive gravity [10, 11, 12, 13, 14, 15, 16, 17, 18]. (See also
Refs. [19, 20] for reviews.)

Unfortunately, the dRGT theory also had some serious challenges: there was no
stable homogeneous and isotropic cosmological solutions [21, 22, 23, 24],
together with other pathologies such as Higuchi bound [25, 26] and positivity
bound [27, 28]. To overcome this, the minimal theory of massive gravity (MTMG)
was then proposed by A. De Felice and S. Mukohyama in 2016 [29] by imposing
some constraints which suppress the five degrees of freedom in the original
dRGT theory such that there are only two degrees of freedom, both of them are
tensor modes, as in the case of general relativity, but now the theory is not
Lorentz invariant. It has the same Friedmann–Lemaître–Robertson–Walker (FLRW)
equations as in the dRGT theory, but now the FLRW background is stable [30].
There are two branches of solutions. The first is the self-accelerating
branch, which is phenomenologically the same with the $\Lambda$CDM cosmology,
except that the accelerating expansion of the universe is now caused by the
graviton mass term, not necessarily by the cosmological constant. The second
is the normal branch, which is phenomenologically different from general
relativity in the scalar and tensor sectors, leading to nontrivial dynamics
which could be tested against the predictions of general relativity [30, 31].

Another attempt to modify the dRGT theory was done by Q.-G. Huang, Y.-S. Piao,
and S.-Y. Zhou in 2012 by coupling the graviton potentials to a scalar field
$\psi$, enabling the graviton to have a varying mass [32]. This theory of
mass-varying massive gravity (MVMG) is again Boulware-Deser ghost-free and the
Lorentz invariance is satisfied. Varying the graviton mass can lead to
interesting cosmological behaviors both in the inflationary and late-time era.
Specifically, the graviton mass will asymptotically approach zero at late
times due to the dynamics of the theory, hence there is no need to fine tune
the graviton mass to a very small number to be in line with the cosmological
bounds for the graviton mass in the present time [33]. However, this may be a
disadvantage, since it means that the contribution of massive gravity to
explain the cosmic expansion at late times is minimum [34].

Another disadvantage of the MVMG theory is that it has many graviton degrees
of freedom, which may lead to instabilities. Therefore, to suppress these
degrees of freedom, we will follow in this paper the method of Refs. [29, 30].
First, we define the precursor theory by writing the MVMG action using the
vielbein formulation in the Arnowitt–Deser–Misner (ADM) formalism, but here we
generalize the theory to the case of higher dimensions. We also adopt the
vielbein potential in Ref. [35] and couple it to the masslike scalar potential
$W(\psi)$, unlike in the previous literature where it is coupled to metric.
The purpose of this is such that the number of graviton degrees of freedom in
the theory is the same as in general relativity. We then perform the Legendre
transformation to the precursor action to obtain the precursor Hamiltonian.
After imposing the nontrivial constraints to the theory, we will obtain the
minimal theory of mass-varying massive gravity (MTMVMG), where the number of
graviton degrees of freedom in this theory becomes ${D(D-3)}/2$, as in the
$D$-dimensional general relativity, where $D$ is the number of spacetime
dimensions. For $D=4$, the number of graviton degrees of freedom in the MTMVMG
is two, in contrast to the MVMG theory where there are five graviton degrees
of freedom. Therefore, in the light of minimally modified gravity (MMG), a
modified theory of gravity with two local gravitational degrees of freedom,
discussed in Refs. [36, 37], where the type-I MMG is for the theory in which
there exists an Einstein frame and the type-II MMG is for the one in which
there is no Einstein frame [38, 39], the MTMVMG can be viewed as an extended
type-II MMG theory.

To study the cosmological aspects of this theory, we can perform the Legendre
transformation again to get the expression for the MTMVMG action, which then
can be used to obtain the Friedmann-Lemaître equations. We take both the
scalar potential and the graviton mass couplings to have exponential forms,
and find that there are eight critical points in the theory: five in the
massless sector and three in the massive sector. Therefore, the MTMVMG theory
can have both massless and massive sectors even in the late-time era, in
contrast to the ordinary MVMG theory where, as mentioned previously, the
dynamics of the theory will lead the graviton mass to asymptotically approach
zero at late times. This makes the MTMVMG a richer theory, which then can give
us good descriptions of both the inflationary and late-time era. Especially,
there are at least two interesting possible scenarios for the late-time
cosmology: the dark energy is either due to the constant graviton mass which
comes from the scalar field $\psi_{\infty}$ that becomes frozen after the
reheating era, or due to the quintessence paradigm where the scalar field
$\psi$ is dynamic. Therefore, if the accelerating expansion of the universe in
the massless sector can be explained by standard quintessence paradigm, in the
massive sector it has to be explained by the nontrivial interplay between
quintessence and massive gravity.

This paper is organized as the following: In Sec. II we derive the precursor
action by writing the MVMG theory using the vielbein formulation. In Sec. III
we then construct the minimal theory by imposing $D$-constraints. In Sec. IV
we derive the Friedmann-Lemaître equations for our model. We then perform in
Sec. V the dynamical analysis around the critical points and check in Sec. VI
their local and global existences. In Sec. VII we discuss the cosmological
implication on the inflationary expansion and late-time acceleration. We then
conclude the paper and write several remarks in Sec. VIII. Detailed
calculations are presented in the Appendix.

##  II Precursor Theory

In this section we will construct the MVMG action by replacing the graviton
mass with the masslike scalar potential $W(\psi)$. However, we will adopt the
vielbein potential from Ref. [35] and couple it to $W(\psi)$, in contrast to
Ref. [32] where $W(\psi)$ is coupled to metric. We will then define the action
for the precursor theory, which will be needed later to construct the MTMVMG
action.

Let us first consider two $D$-dimensional Lorentzian manifolds
$(\mathcal{M},g)$ and $(\mathcal{M}_{0},g_{0})$ parametrized by the coordinate
systems $x^{\mu}$ and $y^{a}$ with $\mu,a=0,\ldots,D-1$, respectively. Using
the ADM formalism, the metrics $g$ and $g_{0}$ can be written as

$\displaystyle ds^{2}_{(d)}$ $\displaystyle=$
$\displaystyle-N^{2}dt^{2}+\gamma_{ij}(dx^{i}+N^{i}dt)(dx^{j}+N^{j}dt),$ (1)
$\displaystyle ds^{2}_{(b)}$ $\displaystyle=$ $\displaystyle-
N_{0}^{2}d\tau^{2}+\gamma_{0ij}(dy^{i}+N_{0}^{i}d\tau)(dy^{j}+N_{0}^{j}d\tau),$
(2)

where $\gamma_{ij}$ and $\gamma_{0ij}$ are the components of the induced
spatial metrices, $N$ and $N_{0}$ are the lapse functions, $N^{i}$ and
$N_{0}^{i}$ are the shift vectors, and $i,j=1,\ldots,D-1$. The subscripts
$(d)$ and $(b)$ denote the dynamic and the background, respectively. Let us
then introduce the vielbeins $\tensor{E}{{}^{A}_{\mu}}$ and
$\tensor{E}{{}_{0}^{A}{}_{a}}$ such that the metrics (1) and (2) can be
written as
$g_{\mu\nu}=\eta_{AB}\tensor{E}{{}^{A}_{\mu}}\tensor{E}{{}^{B}_{\nu}}$ and
$g_{0ab}=\eta_{AB}\tensor{E}{{}_{0}^{A}{}_{a}}\tensor{E}{{}_{0}^{B}{}_{b}}$,
where $A,B=0,\ldots,D-1$ and $\eta_{AB}$ is the flat Minkowski metric.
Therefore, the vielbeins $\tensor{E}{{}^{A}_{\mu}}$ and
$\tensor{E}{{}_{0}^{A}{}_{a}}$ have the form

$\tensor{E}{{}^{A}_{\mu}}=\begin{pmatrix}N&0\\\
N^{k}\tensor{e}{{}^{I}_{k}}&\tensor{e}{{}^{I}_{i}}\end{pmatrix},\qquad\tensor{E}{{}_{0}^{A}{}_{a}}=\begin{pmatrix}N_{0}&0\\\
N_{0}^{k}\tensor{e}{{}_{0}^{I}{}_{k}}&\tensor{e}{{}_{0}^{I}{}_{i}}\end{pmatrix},$
(3)

whose duals are given by

$\tensor{E}{{}^{\mu}_{A}}=\begin{pmatrix}\frac{1}{N}&0\\\
-\frac{N^{i}}{N}&\tensor{e}{{}^{i}_{I}}\end{pmatrix},\qquad\tensor{E}{{}_{0}^{a}{}_{A}}=\begin{pmatrix}\frac{1}{N_{0}}&0\\\
-\frac{N_{0}^{i}}{N_{0}}&\tensor{e}{{}_{0}^{i}{}_{I}}\end{pmatrix},$ (4)

such that they satisfy

$\displaystyle\tensor{E}{{}^{A}_{\mu}}\tensor{E}{{}^{\mu}_{B}}=\tensor{\delta}{{}^{A}_{B}},$
$\displaystyle\tensor{E}{{}_{0}^{A}{}_{a}}\tensor{E}{{}_{0}^{a}{}_{B}}=\tensor{\delta}{{}^{A}_{B}},$
$\displaystyle\tensor{E}{{}^{\mu}_{A}}\tensor{E}{{}^{A}_{\nu}}=\tensor{\delta}{{}^{\mu}_{\nu}},$
$\displaystyle\tensor{E}{{}_{0}^{a}{}_{A}}\tensor{E}{{}_{0}^{A}{}_{b}}=\tensor{\delta}{{}^{a}_{b}}.$

Here $\tensor{e}{{}^{I}_{i}}$ and $\tensor{e}{{}_{0}^{I}{}_{i}}$ are the
spatial vielbeins, with $i,I=1,\ldots,D-1$, and $\tensor{e}{{}^{i}_{I}}$ and
$\tensor{e}{{}_{0}^{i}{}_{I}}$ are their duals, respectively.

Now, we consider a smooth embedding
$\phi:\mathcal{M}\longrightarrow\mathcal{M}_{0}$ such that we could have
pulled back quantities

$\displaystyle\tensor{\tilde{E}}{{}^{A}_{\mu}}(x)$ $\displaystyle=$
$\displaystyle\frac{\partial\phi^{a}}{\partial
x^{\mu}}\tensor{E}{{}_{0}^{A}{}_{a}}(\phi(x)),$ (5) $\displaystyle
f_{\mu\nu}(x)$ $\displaystyle=$ $\displaystyle\frac{\partial\phi^{a}}{\partial
x^{\mu}}\frac{\partial\phi^{b}}{\partial x^{\nu}}g_{0ab}(\phi(x)),$ (6)

so that $\mathcal{M}$ obeys diffeomorphism rules through a Stückelberg field
$\phi^{a}$. In the unitary gauge where
$\phi^{a}=\tensor{\delta}{{}^{a}_{\mu}}x^{\mu}$, we simply recover
$\tensor{\tilde{E}}{{}^{A}_{\mu}}=\tensor{E}{{}_{0}^{A}{}_{\mu}}$ and
$f_{\mu\nu}=g_{0\mu\nu}$. Note that the original formulation of the dRGT
theory has a Minkowski background, $g_{0ab}=\eta_{ab}$. Using the vielbein
formulation in Ref. [35], the ghost-free potential related to the graviton
mass has the form

$\sum_{n=0}^{D}\frac{c_{n}}{n!(D-n)!}\hat{\epsilon}_{A_{1}A_{2}\cdots
A_{D}}\tensor{\bm{E}}{{}^{A_{1}}}\wedge\cdots\wedge\tensor{\bm{E}}{{}^{A_{n}}}\wedge\tensor{\bm{\tilde{E}}}{{}^{A_{n+1}}}\wedge\cdots\wedge\tensor{\bm{\tilde{E}}}{{}^{A_{D}}}$
(7)

with one-forms $\bm{E}^{A}=\tensor{E}{{}^{A}_{\mu}}dx^{\mu}$ and
$\bm{\tilde{E}}^{A}=\tensor{\tilde{E}}{{}^{A}_{\mu}}dx^{\mu}$ the dRGT mass
term for arbitrary background and dimension. The quantity $\hat{\epsilon}$
denotes the Levi-Civita symbol in the flat spacetime. Furthermore, inspired by
[32], the graviton mass is replaced by a function of a scalar field $W(\psi)$,
where $\psi(x)$ is well defined on $\mathcal{M}$. Coupling $W(\psi)$ to the
potential in Eq. (7) and adding this coupling term to the
Einstein–Klein–Gordon action, we obtain the ghost-free MVMG action as the
following,

$\displaystyle S_{\text{MVMG}}$ $\displaystyle=$
$\displaystyle\int_{\mathcal{M}}d^{D}x\,\det{(E)}\,\left(\frac{M_{\text{Pl}}^{D-2}}{2}R(E)-\frac{1}{2}\partial^{\mu}\psi\partial_{\mu}\psi-V(\psi)\right)$
(8)
$\displaystyle-\,\frac{1}{4}\int_{\mathcal{M}}{W(\psi)}\left(\sum_{n=0}^{D}\frac{c_{n}}{n!(D-n)!}\hat{\epsilon}_{A_{1}A_{2}\cdots
A_{D}}\tensor{\bm{E}}{{}^{A_{1}}}\wedge\cdots\wedge\tensor{\bm{E}}{{}^{A_{n}}}\wedge\tensor{\bm{\tilde{E}}}{{}^{A_{n+1}}}\wedge\cdots\wedge\tensor{\bm{\tilde{E}}}{{}^{A_{D}}}\right),$

where $V(\psi)$ is the scalar potential function.

The action for the precursor theory can be obtained by simply substituting the
vielbeins in (3) and (4) to the action (8),

$\displaystyle S_{\text{pre}}$ $\displaystyle=$ $\displaystyle\int
d^{D}x\,N\det{(e)}\Bigg{[}\frac{M_{\text{Pl}}^{D-2}}{2}\left({}^{(D-1)}R(e)+K^{ij}K_{ij}-K^{2}\right)+\frac{1}{2N^{2}}\dot{\psi}^{2}-\frac{1}{2}\partial^{i}\psi\partial_{i}\psi+\frac{N^{i}N^{j}}{2N^{2}}\partial_{i}\psi\partial_{j}\psi-\frac{N^{i}}{N^{2}}\dot{\psi}\partial_{i}\psi-V(\psi)$
(9)
$\displaystyle-\,W(\psi)\left(|\det{(X)}|\frac{M}{N}\sum_{n=0}^{D-1}c_{n}\mathcal{S}_{n}(Y)+\sum_{n=0}^{D-1}c_{D-n}\mathcal{S}_{n}(X)\right)\Bigg{]},$

where ${}^{(D-1)}R(e)$ is the spatial Ricci scalar, while $K^{ij}$ and $K$ are
the second fundamental form and the mean curvature, respectively. For any
function $\psi$, we define $\dot{\psi}\equiv{\partial\psi}/{\partial t}$ and
$\partial_{i}\psi\equiv{\partial\psi}/{\partial x^{i}}$. Here, the quantities
$\\{M,M^{i},\tensor{\tilde{e}}{{}^{I}_{i}}\\}$ are the pull back of
$\\{N_{0},N_{0}^{i},\tensor{e}{{}_{0}^{I}{}_{i}}\\}$ via Stückelberg field
$\phi$ as a background image on $\mathcal{M}$ given by Eq. (5) and (6),
respectively. The elements $\mathcal{S}_{n}$ are the $n$th order symmetric
polynomials which depend on either
$\tensor{Y}{{}^{I}_{J}}\equiv\tensor{e}{{}^{I}_{k}}\tensor{\tilde{e}}{{}^{k}_{J}}$
or
$\tensor{X}{{}^{I}_{J}}\equiv\tensor{\tilde{e}}{{}^{I}_{k}}\tensor{e}{{}^{k}_{J}}$
(see Appendix A for more discussions). Note that the action (9) violates the
local Lorentz symmetry because we have used the ADM vielbeins in Eqs. (3) and
(4).

##  III Minimal Theory of MVMG

The discussion in this section is divided into two parts. First, we construct
the Hamiltonian for the MTMVMG and identify some constraints which restricts
the graviton degrees of freedom. Second, we discuss the MTMVMG action which is
constructed from the precursor action discussed in the previous section but
with some additional constraints.

###  III.1 MTMVMG Hamiltonian and some constraints

Let us first consider the spatial vielbeins $\tensor{e}{{}^{I}_{i}}$ and the
scalar field $\psi$ as the canonical variables which correspond to the
conjugate momenta defined as

$\displaystyle\tensor{\pi}{{}^{i}_{I}}$ $\displaystyle\equiv$
$\displaystyle\frac{\delta
S_{\text{pre}}}{\delta\tensor{\dot{e}}{{}^{I}_{i}}}$ (10) $\displaystyle=$
$\displaystyle\det{(e)}\,M_{\text{Pl}}^{D-2}(K^{ij}-K\gamma^{ij})\delta_{IJ}\tensor{e}{{}^{J}_{j}},$
$\displaystyle\pi$ $\displaystyle\equiv$ $\displaystyle\frac{\delta
S_{\text{pre}}}{\delta\dot{\psi}}$ (11) $\displaystyle=$
$\displaystyle\det{(e)}\,\left(\frac{1}{N}\dot{\psi}-\frac{N^{i}}{N}\partial_{i}\psi\right).$

We can switch from Lagrangian to Hamiltonian by performing the Legendre
transformation in order to see some constraints of the theory. As it is well-
known in the vielbein language that the lapse function $N$ and the shift
vector $N^{i}$ appear as Lagrange multipliers enforcing the diffeomorphism
constraints [35], namely $\mathcal{R}_{0}\approx 0$ and
$\mathcal{R}_{i}\approx 0$. These constraints are called the primary
constraints of the first kind, which then enable us to construct another set
of constraints called the secondary constraints of the first kind [40]. Note
that there are only $D-2$ independent secondary constraints, since two of them
can be obtained from the others. We denote them as
$\tilde{\mathcal{C}}_{\tau}$ $(\tau=1,\ldots,D-2)$, together with their
Lagrange multipliers $\lambda^{\tau}$.

Additionally, as studied in Refs. [29, 30], from Eq. (10) we also have another
set of primary constraints of the second kind $\mathcal{P}^{[MN]}$ that lead
to the secondary constraints of the second kind $\mathcal{Z}^{[MN]}$ in the
phase space, together with their Lagrange multipliers $\alpha_{MN}$ and
$\beta_{MN}$, with $M,N=1,\ldots,(D-1)(D-2)/2$. These secondary constraints
are necessary since the primary constraints should be preserved with respect
to the time evolution. Therefore, the precursor Hamiltonian can be written
down as

$\displaystyle H_{\text{pre}}$ $\displaystyle=$ $\displaystyle\int
d^{D-1}x\,\Big{(}-N\mathcal{R}_{0}-N^{i}\mathcal{R}_{i}$ (12)
$\displaystyle+\,W(\psi)N\mathcal{H}_{0}+W(\psi)M\mathcal{H}_{1}+\tilde{\lambda}^{\tau}\tilde{\mathcal{C}}_{\tau}$
$\displaystyle+\,\alpha_{MN}\mathcal{P}^{[MN]}+\beta_{MN}\mathcal{Z}^{[MN]}\Big{)},$

where

$\displaystyle\mathcal{R}_{0}$ $\displaystyle\equiv$
$\displaystyle\det{(e)}\,\frac{M_{\text{Pl}}^{D-2}}{2}{}^{(D-1)}R(e)$ (13)
$\displaystyle-\,\frac{1}{2\det{(e)}\,M_{\text{Pl}}^{D-2}}\left[\tensor{\pi}{{}^{i}_{I}}\tensor{\pi}{{}^{I}_{i}}-\frac{1}{D-2}(\tensor{\pi}{{}^{i}_{I}}\tensor{e}{{}^{I}_{i}})^{2}\right],$
$\displaystyle-\,\frac{1}{2\det{(e)}}\pi^{2}-\frac{\det{(e)}}{2}\partial_{i}\psi\partial^{i}\psi$
$\displaystyle-\,\det{(e)}\,V(\psi),$ $\displaystyle\mathcal{R}_{i}$
$\displaystyle\equiv$
$\displaystyle\nabla_{j}(\tensor{\pi}{{}^{j}_{I}}\tensor{e}{{}^{I}_{i}})-\pi\partial_{i}\psi,$
(14) $\displaystyle\mathcal{H}_{0}$ $\displaystyle\equiv$
$\displaystyle\det{(e)}\,\sum_{n=0}^{D-1}c_{D-n}\mathcal{S}_{n}(X),$ (15)
$\displaystyle\mathcal{H}_{1}$ $\displaystyle\equiv$
$\displaystyle\det{(e)}\,|\det{(X)}|\sum_{n=0}^{D-1}c_{n}\mathcal{S}_{n}(Y),$
(16)

and

$\displaystyle\mathcal{P}^{[MN]}$ $\displaystyle\equiv$
$\displaystyle(\tensor{e}{{}^{M}_{j}}\delta^{KN}-\tensor{e}{{}^{N}_{j}}\delta^{KM})\tensor{\pi}{{}^{j}_{K}},$
(17) $\displaystyle\mathcal{Z}^{[MN]}$ $\displaystyle\equiv$
$\displaystyle(\tensor{e}{{}^{M}_{j}}\delta^{KN}-\tensor{e}{{}^{N}_{j}}\delta^{KM})\tensor{\tilde{e}}{{}^{j}_{K}}.$
(18)

Since the constraints above remove some graviton degrees of freedom, we can
construct a theory in which the spatial graviton degrees of freedom coincide
with the standard general relativity. Inspired by Ref. [29], we can impose the
$D$-constraints in unitary gauge given by

$\displaystyle\mathcal{C}_{0}$ $\displaystyle\equiv$
$\displaystyle\\{\mathcal{R}_{0},H_{1}\\}_{\text{PB}}-W(\psi)\frac{\partial\mathcal{H}_{0}}{\partial
t}\approx 0,$ (19) $\displaystyle\mathcal{C}_{i}$ $\displaystyle\equiv$
$\displaystyle\\{\mathcal{R}_{i},H_{1}\\}_{\text{PB}}\approx 0,$ (20)

where $\\{\cdots\\}_{\text{PB}}$ denotes the Poisson bracket and

$H_{1}\equiv\int d^{D-1}x\,W(\psi)M\mathcal{H}_{1}.$ (21)

Note that these constraints consist of two new constraints and $D-2$
independent constraints $\tilde{\mathcal{C}}_{\tau}$ which already exist in
the precursor theory. Moreover, the constraints in Eqs. (19) and (20) imply
that the theory admits the Lorentz symmetry violation. The Hamiltonian of the
MTMVMG theory then reads

$\displaystyle H_{\text{MTMVMG}}$ $\displaystyle=$ $\displaystyle\int
d^{D-1}x\,\Big{[}-N\mathcal{R}_{0}-N^{i}\mathcal{R}_{i}$ (22)
$\displaystyle+\,W(\psi)(N\mathcal{H}_{0}+M\mathcal{H}_{1})+\lambda\mathcal{C}_{0}+\lambda^{i}\mathcal{C}_{i}$
$\displaystyle+\,\alpha_{MN}\mathcal{P}^{[MN]}+\beta_{MN}\mathcal{Z}^{[MN]}\Big{]}.$

Thus, in total we have $D^{2}-D+2$ constraints, which means that the number of
spatial graviton degrees of freedom are $D(D-3)/2$, as in the $D$-dimensional
general relativity.

###  III.2 MTMVMG action

To construct the MTMVMG action, one has to employ the Legendre transformation
on the Hamiltonian functional (22). It will be shortly discussed in this
subsection, but its detailed derivations will be presented in the Appendix B.

As discussed in the previous subsection, we should have the nontrivial
$D$-constraints (19) and (20) in the MTMVMG theory. In order to have a
consistent theory, we have to modify the conjugate momenta (10) and (11) to

$\displaystyle\frac{\tensor{\pi}{{}^{i}_{I}}}{\det{(e)}}$
$\displaystyle\equiv$ $\displaystyle
M_{\text{Pl}}^{D-2}(K^{ij}\delta_{IJ}\tensor{e}{{}^{J}_{j}}-K\tensor{e}{{}^{i}_{I}})$
(23)
$\displaystyle-\,\lambda\frac{W(\psi)}{2}\frac{M}{N}\Theta^{ij}\delta_{IJ}\tensor{e}{{}^{J}_{j}},$
$\displaystyle\frac{\pi}{\det{(e)}}$ $\displaystyle\equiv$
$\displaystyle\frac{\dot{\psi}}{N}-\frac{N^{i}}{N}\partial_{i}\psi-\lambda\frac{dW}{d\psi}\frac{M}{N}\Phi,$
(24)

with

$\displaystyle\Theta^{ij}$ $\displaystyle\equiv$
$\displaystyle-|\det{(X)}|\delta^{IK}(\tensor{e}{{}^{i}_{K}}\tensor{\tilde{e}}{{}^{j}_{J}}+\tensor{e}{{}^{j}_{K}}\tensor{\tilde{e}}{{}^{i}_{J}})$
(25)
$\displaystyle\times\,\sum_{n=1}^{D-1}\sum_{m=1}^{n}(-1)^{m}c_{n}\tensor{\left(Y^{m-1}\right)}{{}^{J}_{I}}\mathcal{S}_{n-m}(Y),$
$\displaystyle\Phi$ $\displaystyle\equiv$
$\displaystyle|\det{(X)}|\sum_{n=1}^{D-1}c_{n}\mathcal{S}_{n}(Y).$ (26)

The notation $\tensor{(M^{m})}{{}^{I}_{J}}$ means

$\tensor{(M^{m})}{{}^{I}_{J}}\equiv\tensor{M}{{}^{I}_{K_{1}}}\tensor{M}{{}^{K_{1}}_{K_{2}}}\cdots\tensor{M}{{}^{K_{m-1}}_{J}}.$
(27)

The modifications (23) and (24) imply that the MTMVMG theory modifies both the
kinetic part and the mass term. As we will see later, this also provides a
class of solutions which coincides in the dRGT theory in the FLRW background.

For the sake of convenience, let us first introduce the following tensors

$\tensor{\mathcal{K}}{{}^{i}_{j}}\equiv\tensor{\tilde{e}}{{}^{i}_{I}}\tensor{e}{{}^{I}_{j}},\qquad\tensor{\mathcal{\bar{K}}}{{}^{i}_{j}}\equiv\tensor{e}{{}^{i}_{I}}\tensor{\tilde{e}}{{}^{I}_{j}}.$
(28)

satisfying
$\tensor{\mathcal{K}}{{}^{i}_{k}}\tensor{\mathcal{\bar{K}}}{{}^{k}_{j}}=\tensor{\delta}{{}^{i}_{j}}$,
which correspond to the spatial metrics by

$\tensor{\mathcal{K}}{{}^{i}_{k}}\tensor{\mathcal{K}}{{}^{k}_{j}}=\tilde{\gamma}^{il}\gamma_{lj},\qquad\tensor{\mathcal{\bar{K}}}{{}^{i}_{k}}\tensor{\mathcal{\bar{K}}}{{}^{k}_{j}}=\gamma^{il}\tilde{\gamma}_{lj},$
(29)

where
$\tilde{\gamma}_{ij}=\delta_{IJ}\tensor{\tilde{e}}{{}^{I}_{i}}\tensor{\tilde{e}}{{}^{J}_{j}}$
is the spatial metric on $\mathcal{M}$. Performing the Legendre
transformation, we obtain the MTMVMG action,

$S_{\text{MTMVMG}}=S_{\text{pre}}+S_{\lambda},$ (30)

where

$\displaystyle S_{\lambda}$ $\displaystyle=$
$\displaystyle\frac{2}{M_{\text{Pl}}^{D-2}}\int
d^{D}x\,N\sqrt{\gamma}\left(\lambda\frac{W(\psi)}{4}\frac{M}{N}\right)^{2}\left(\gamma_{ik}\gamma_{jl}-\frac{1}{D-2}\gamma_{ij}\gamma_{kl}\right)\Theta^{ij}\Theta^{kl}$
(31) $\displaystyle+\,\frac{1}{2}\int
d^{D}x\,N\sqrt{\gamma}\left(\lambda\frac{dW}{d\psi}\frac{M}{N}\right)^{2}\Phi^{2}-\int
d^{D}x\,\sqrt{\gamma}\left[\lambda\mathcal{\bar{C}}_{0}+\lambda^{i}\tensor{\mathcal{C}}{{}_{i}}\right],$

and

$\displaystyle\mathcal{\bar{C}}_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{2}W(\psi)M\left(\gamma_{ik}\gamma_{jl}-\frac{1}{D-2}\gamma_{ij}\gamma_{kl}\right)\Theta^{kl}(K^{ij}-K\gamma^{ij})$
(32)
$\displaystyle+\,W(\psi)|\det{(\mathcal{\bar{K}})}|\sum_{n=1}^{D-1}\sum_{m=1}^{n}(-1)^{m}c_{n}\tensor{\left(\mathcal{K}^{m-1}\right)}{{}^{k}_{l}}\tensor{\tilde{\zeta}}{{}^{l}_{k}}\mathcal{S}_{n-m}(\mathcal{K})$
$\displaystyle-\,M|\det{(\mathcal{\bar{K}})}|\frac{dW}{d\psi}\left(\frac{\dot{\psi}}{N}-\frac{N^{i}}{N}\partial_{i}\psi\right)\sum_{n=1}^{D-1}c_{n}\mathcal{S}_{n}(\mathcal{K}),$
$\displaystyle\tensor{\mathcal{C}}{{}_{i}}$ $\displaystyle=$
$\displaystyle-W(\psi)\nabla_{k}M|\det{(\mathcal{\bar{K}})}|\sum_{n=1}^{D-1}\sum_{m=1}^{n}(-1)^{m}c_{n}\tensor{\left(\mathcal{K}^{m}\right)}{{}^{k}_{i}}\mathcal{S}_{n-m}(\mathcal{K})$
(33)
$\displaystyle-\,M|\det{(\mathcal{\bar{K}})}|\partial_{i}\psi\frac{dW}{d\psi}\sum_{n=1}^{D-1}c_{n}\mathcal{S}_{n}(\mathcal{K}),$

with

$\displaystyle\Theta^{ij}$ $\displaystyle=$
$\displaystyle-2|\det{(\mathcal{\bar{K}})}|\gamma^{il}\sum_{n=1}^{D-1}\sum_{m=1}^{n}(-1)^{m}c_{n}\tensor{\left(\mathcal{K}^{m}\right)}{{}^{j}_{l}}\mathcal{S}_{n-m}(\mathcal{K})$
(34) $\displaystyle\Phi$ $\displaystyle=$
$\displaystyle|\det{(\mathcal{\bar{K}})}|\sum_{n=1}^{D-1}c_{n}\mathcal{S}_{n}(\mathcal{K}).$
(35)

It is worth mentioning that the additional term in (23) implies that we have
to set $\alpha_{MN}=\beta_{MN}=0$ [30, 41]. This is so because the tensors
$\mathcal{P}^{[MN]}$ and $Y^{[MN]}$ are antisymmetric, while the tensor
$\Theta^{ij}$ is symmetric.

We could extend the action (30) by adding the matter field,

$S_{\text{MTMVMG-M}}=S_{\text{pre}}+S_{\lambda}+S_{\text{matter}}.$ (36)

Here, we consider the matter field part $S_{\text{matter}}$ to be the perfect
fluid whose energy-momentum tensor has the form

$T_{\mu\nu}=\frac{2}{\sqrt{-g}}\frac{\delta S_{\text{matter}}}{\delta
g^{\mu\nu}}=\rho_{m}U_{\mu}U_{\nu}+P_{m}\left(g_{\mu\nu}+U_{\mu}U_{\nu}\right),$
(37)

where $U^{\mu}$ and $\rho_{m}$ are the unit velocity of the fluid and the
energy density, respectively. The pressure $P_{m}$ is given by the state
equation of matter fields,

$P_{m}=w_{m}\rho_{m},$ (38)

with $w_{m}$ is a real constant [42]. In the standard higher-dimensional
cosmology, we particularly have $w_{m}=\frac{1}{D-1}$ (radiation), $w_{m}=0$
(dust), and $w_{m}=-1$ (vacuum) [43].

##  IV Friedmann-Lemaître Equations

In this section we consider a cosmological model in the MTMVMG theory. Our
starting point is to write down the metric ansatz for the dynamic and the
background manifolds which are spatially flat,

$\displaystyle ds_{(d)}^{2}$ $\displaystyle=$ $\displaystyle
g_{\mu\nu}dx^{\mu}dx^{\nu}=-N^{2}(t)dt^{2}+a^{2}(t)\delta_{ij}dx^{i}dx^{j},$
$\displaystyle ds_{(b)}^{2}$ $\displaystyle=$ $\displaystyle
f_{\mu\nu}dx^{\mu}dx^{\nu}=-M^{2}(t)dt^{2}+\tilde{a}^{2}(t)\delta_{ij}dx^{i}dx^{j}.$

In the case at hand, the action (36) simplifies to

$\displaystyle S_{\text{MTMVMG-M}}$ $\displaystyle=$ $\displaystyle-\int
d^{D}x\,a^{D-1}\Bigg{\\{}(D-1)M_{\text{Pl}}^{D-2}NH^{2}-\frac{1}{2N}\dot{\psi}^{2}+NV(\psi)$
(39)
$\displaystyle-\,W(\psi)\bigg{(}M\sum_{n=0}^{D-1}c_{n}A_{n}u^{D-n-1}+N\sum_{n=0}^{D-1}c_{D-n}A_{n}u^{n}\bigg{)}$
$\displaystyle+\,\frac{\lambda^{2}M^{2}}{N}\Bigg{[}\frac{D-1}{D-2}\frac{2}{M_{\text{Pl}}^{D-2}}\left(\frac{W(\psi)}{4}\sum_{n=1}^{D-1}c_{n}B_{n}u^{D-n-1}\right)^{2}$
$\displaystyle+\,\frac{1}{2}\left(\frac{dW}{d\psi}\sum_{n=1}^{D-1}c_{n}A_{n}u^{D-n-1}\right)^{2}\Bigg{]}-\lambda\Bigg{[}\frac{\dot{\psi}M}{N}\frac{dW}{d\psi}\sum_{n=1}^{D-1}c_{n}A_{n}u^{D-n-1}$
$\displaystyle+\,(D-1)W(\psi)M\left(H-uH_{f}\right)\sum_{n=1}^{D-1}c_{n}B_{n}u^{D-n-1}\Bigg{]}\Bigg{\\}}+S_{\text{matter}}.$

Here, we have introduced the quantities

$H\equiv\frac{\dot{a}}{Na},\qquad
H_{f}\equiv\frac{\dot{\tilde{a}}}{M\tilde{a}},\qquad
u\equiv\frac{\tilde{a}}{a},$ (40)

and parameters

$\displaystyle A_{n}$ $\displaystyle\equiv$
$\displaystyle\frac{(D-1)!}{(D-n-1)!n!},$ (41) $\displaystyle B_{n}$
$\displaystyle\equiv$
$\displaystyle\sum_{m=1}^{n}(-1)^{m}\frac{(D-1)!}{(D-n+m-1)!(n-m)!}.$ (42)

Then, the variation with respect to the lapse function $N(t)$ gives us the
first Friedmann-Lemaître equation,

$\frac{1}{2}(D-1)(D-2)H^{2}=\frac{1}{M_{\text{Pl}}^{D-2}}(\rho_{m}+\rho_{\text{MG}}+\rho_{\lambda}),$
(43)

where

$\displaystyle\rho_{\text{MG}}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2N^{2}}\dot{\psi}^{2}+V(\psi)+W(\psi)\sum_{n=0}^{D-1}c_{n+1}A_{n}u^{D-n-1},$
(44) $\displaystyle\rho_{\lambda}$ $\displaystyle\equiv$
$\displaystyle\frac{\lambda\dot{\psi}M}{N^{2}}\frac{dW}{d\psi}\sum_{n=1}^{D-1}c_{n}A_{n}u^{D-n-1}-\frac{D-1}{D-2}\frac{\lambda^{2}M^{2}W^{2}(\psi)}{2M_{\text{Pl}}^{D-2}N^{2}}\left(\sum_{n=1}^{D-1}c_{n}B_{n}u^{D-n-1}\right)^{2}$
(45) $\displaystyle+\,\frac{(D-1)\lambda
HMW(\psi)}{N}\sum_{n=1}^{D-1}c_{n}B_{n}u^{D-n-1}+\frac{\lambda^{2}M^{2}}{2N^{2}}\left(\frac{dW}{d\psi}\right)^{2}\left(\sum_{n=1}^{D-1}c_{n}A_{n}u^{D-n-1}\right)^{2}.$

The variation with respect to the scale factor $a(t)$ gives us the second
Friedmann-Lemaître equation,

$(D-2)\frac{\dot{H}}{N}+\frac{1}{2}(D-1)(D-2)H^{2}=-\frac{1}{M_{\text{Pl}}^{D-2}}(P_{m}+P_{\text{MG}}+P_{\lambda}),$
(46)

where

$\displaystyle P_{\text{MG}}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2N^{2}}\dot{\psi}^{2}-V(\psi)-\frac{W(\psi)}{N}\sum_{n=1}^{D-1}\left(Mc_{n}+Nc_{n+1}\right)nA_{n}u^{D-n-1},$
(47) $\displaystyle P_{\lambda}$ $\displaystyle\equiv$
$\displaystyle\frac{\lambda^{2}M^{2}W^{2}(\psi)}{2M_{\text{Pl}}^{D-2}N^{2}}\left(\sum_{n=1}^{D-1}\frac{D-2n-1}{D-2}c_{n}B_{n}u^{D-n-1}\right)\left(\sum_{n=1}^{D-1}c_{n}B_{n}u^{D-n-1}\right)$
(48)
$\displaystyle-\,\frac{\lambda^{2}M^{2}}{2N^{2}}\left(\frac{dW}{d\psi}\right)^{2}\left(\sum_{n=1}^{D-1}\frac{D-2n-1}{D-1}c_{n}A_{n}u^{D-n-1}\right)\left(\sum_{n=1}^{D-1}c_{n}A_{n}u^{D-n-1}\right)$
$\displaystyle-\,\frac{\lambda
MW(\psi)H_{f}}{N}\sum_{n=1}^{D-1}\Bigg{(}(D-n-1)\frac{M}{N}+(n-1)u\Bigg{)}c_{n}B_{n}u^{D-n-1}$
$\displaystyle-\,\left(\frac{\dot{\lambda}MN+\lambda(\dot{M}N-M\dot{N})}{N^{3}}\right)W(\psi)\sum_{n=1}^{D-1}c_{n}B_{n}u^{D-n-1}$
$\displaystyle+\,\frac{\lambda
M\dot{\psi}}{N^{2}}\frac{dW}{d\psi}\sum_{n=1}^{D-1}c_{n}\left(\frac{nA_{n}}{D-1}+B_{n}\right)u^{D-n-1}.$

From the variation with respect to $\psi(t)$, we obtain the equation of
motions,

$\displaystyle\frac{1}{N^{2}}\ddot{\psi}+\left(\frac{(D-1)H}{N}-\frac{\dot{N}}{N^{3}}\right)\dot{\psi}+\frac{dV}{d\psi}+\frac{1}{N}\frac{dW}{d\psi}\Bigg{\\{}\sum_{n=0}^{D-1}\left(Mc_{n}+Nc_{n+1}\right)A_{n}u^{D-n-1}$
$\displaystyle\quad+\,\frac{\lambda
M}{N}\left(NH\sum_{n=1}^{D-1}nc_{n}A_{n}u^{D-n-1}+MH_{f}\sum_{n=1}^{D-1}(D-n-1)c_{n}A_{n}u^{D-n-1}\right)$
$\displaystyle\quad+\,\left(\frac{\dot{\lambda}M}{N}+\frac{\lambda(\dot{M}N-M\dot{N})}{N^{2}}\right)\sum_{n=1}^{D-1}c_{n}A_{n}u^{D-n-1}-\lambda
M\Bigg{[}\frac{\lambda
M}{N}\frac{d^{2}W}{d\psi^{2}}\left(\sum_{n=1}^{D-1}c_{n}A_{n}u^{D-n-1}\right)^{2}$
$\displaystyle\quad-\,\frac{D-1}{D-2}\frac{\lambda
MW(\psi)}{4M_{\text{Pl}}^{D-2}N}\left(\sum_{n=1}^{D-1}c_{n}B_{n}u^{D-n-1}\right)^{2}+(D-1)\left(H-uH_{f}\right)\sum_{n=1}^{D-1}c_{n}B_{n}u^{D-n-1}\Bigg{]}\Bigg{\\}}=0.$
(49)

Performing the variation with respect to $\lambda$ will give us

$\displaystyle
W(\psi)\left(uH_{f}-H\right)\left(\sum_{n=1}^{D-1}c_{n}B_{n}u^{D-n-1}\right)+\frac{\lambda}{D-2}\frac{W^{2}(\psi)M}{M_{\text{Pl}}^{D-2}N}\left(\sum_{n=1}^{D-1}c_{n}B_{n}u^{D-n-1}\right)^{2}$
$\displaystyle\quad-\,\frac{\dot{\psi}}{D-1}\frac{dW}{d\psi}\left(\sum_{n=1}^{D-1}c_{n}A_{n}u^{D-n-1}\right)-\frac{\lambda}{D-1}\frac{M}{N}\left(\frac{dW}{d\psi}\right)^{2}\left(\sum_{n=1}^{D-1}c_{n}A_{n}u^{D-n-1}\right)^{2}=0,$
(50)

which relates the Hubble rate of the background spacetime with the Hubble rate
of the dynamical spacetime. In the simple case where $\psi$ is trivial and
$D=4$, Eq. (50) corresponds to the branches of solutions discussed in [30]. In
the case of $D=5,6$ with trivial $\psi$, one has to solve the qubic and the
quartic polynomials, respectively, while for $D\geq 7$ the solutions of the
polynomials are still unknown. For nontrivial $\psi$, it is still unknown
whether such branches exist. These aspects will be considered elsewhere.

##  V Dynamical System Analysis

Let us consider a special case where the couplings $W(\psi)$ and $V(\psi)$
have the form

$\displaystyle W(\psi)$ $\displaystyle=$ $\displaystyle
W_{0}\exp\left(-\frac{\lambda_{W}\psi}{\sqrt{M_{\text{Pl}}^{D-2}}}\right),$
(51) $\displaystyle V(\psi)$ $\displaystyle=$ $\displaystyle
V_{0}\exp\left(-\frac{\lambda_{V}\psi}{\sqrt{M_{\text{Pl}}^{D-2}}}\right),$
(52)

which has been considered in four dimensional cases [44, 45] where the
constants $V_{0},W_{0}>0$ and $\lambda_{V},\lambda_{W}\geq 0$. Moreover, the
form of the scalar potential $V(\psi)$ may provide inflationary expansion of
the early universe model and has been well studied in the context of dynamical
systems in Ref. [46]. Note that for $\lambda_{W}=0$ we have the MTMVMG theory
with ordinary constant graviton mass [29, 30].

For the rest of the paper we simply take a branch of solutions of Eq. (50)
where $\lambda=0$. Setting the lapse functions $N(t)=M(t)=1$, we introduce the
autonomous variables,

$\displaystyle x_{\rho}$ $\displaystyle\equiv$
$\displaystyle\left(\frac{2\rho_{m}}{M_{\text{Pl}}^{D-2}\lambda_{D}^{2}H^{2}}\right)^{1/2},$
(53) $\displaystyle x_{\psi}$ $\displaystyle\equiv$
$\displaystyle\frac{\dot{\psi}}{\sqrt{M_{\text{Pl}}^{D-2}}\lambda_{D}H},$ (54)
$\displaystyle x_{V}$ $\displaystyle\equiv$
$\displaystyle\left(\frac{2V(\psi)}{M_{\text{Pl}}^{D-2}\lambda_{D}^{2}H^{2}}\right)^{1/2},$
(55) $\displaystyle x_{W}$ $\displaystyle\equiv$
$\displaystyle\left(\frac{2W(\psi)}{M_{\text{Pl}}^{D-2}\lambda_{D}^{2}H^{2}}\right)^{1/2},$
(56)

with $\lambda^{2}_{D}\equiv(D-1)(D-2)$ such that the equations of motion (46),
(49), and (50) can be written down in terms of the autonomous variables,

$\displaystyle\frac{2}{D-1}x_{\psi}^{\prime}$ $\displaystyle=$
$\displaystyle(1-w_{m})x_{\psi}^{3}-(1+w_{m})x_{\psi}x_{V}^{2}$ (57)
$\displaystyle-\,[f_{2}(u)+w_{m}f_{1}(u)]x_{\psi}x_{W}^{2}-(1-w_{m})x_{\psi}$
$\displaystyle+\,2\sqrt{\frac{D-2}{D-1}}[\lambda_{V}x_{V}^{2}+\lambda_{W}f_{1}(u)x_{W}^{2}],$
$\displaystyle\frac{2}{D-1}x_{V}^{\prime}$ $\displaystyle=$
$\displaystyle(1+w_{m})x_{V}+(1-w_{m})x_{\psi}^{2}x_{V}$ (58)
$\displaystyle-\,(1+w_{m})x_{V}^{3}-[f_{2}(u)+w_{m}f_{1}(u)]x_{V}x_{W}^{2}$
$\displaystyle-\,\sqrt{\frac{D-2}{D-1}}\lambda_{V}x_{\psi}x_{V},$
$\displaystyle\frac{2}{D-1}x_{W}^{\prime}$ $\displaystyle=$
$\displaystyle(1+w_{m})x_{W}+(1-w_{m})x_{\psi}^{2}x_{W}$ (59)
$\displaystyle-\,(1+w_{m})x_{V}^{2}x_{W}-[f_{2}(u)+w_{m}f_{1}(u)]x_{W}^{3}$
$\displaystyle-\,\sqrt{\frac{D-2}{D-1}}\lambda_{W}x_{\psi}x_{W},$

where we have used the constraint coming from Eq. (43),

$x_{\rho}^{2}+x_{\psi}^{2}+x_{V}^{2}+f_{1}(u)x_{W}^{2}=1.$ (60)

Note that here the prime symbol denotes the derivative with respect to
$\ln{(a)}$. We also have defined

$\displaystyle f_{1}(u)$ $\displaystyle\equiv$
$\displaystyle\sum_{n=0}^{D-1}c_{n+1}A_{n}u^{D-n-1},$ (61) $\displaystyle
f_{2}(u)$ $\displaystyle\equiv$
$\displaystyle\sum_{n=1}^{D-1}\left(c_{n}+c_{n+1}\right)nA_{n}u^{D-n-1},$ (62)

which are assumed to be bounded functions. As we have seen above, the scale
factor $a(t)$ can be thought of as a parameter in this picture and the
fiducial scale $\tilde{a}(t)$ is only a background in this setup. Therefore,
we conclude that the quantity $u$ in Eq. (40) is not a dynamical variable; it
might be either a function of time, $u(t)$, or a constant. The first case is
called the normal branch, while the latter is called the self-accelerating
branch. Both have been appeared in the context of four-dimensional MTMG [30].
Similar situation also occurs in the cosmological model of dRGT massive
gravity (see, for example, Ref. [47]).

We could also introduce some higher dimensional quantities which are analogous
to the four-dimensional cases. First, the state equation parameter and the
density parameter in this theory are given by

$\displaystyle w_{\text{MG}}$ $\displaystyle\equiv$
$\displaystyle\frac{P_{\text{MG}}}{\rho_{\text{MG}}}=\frac{x_{\psi}^{2}-x_{V}^{2}-f_{2}(u)x_{W}^{2}}{x_{\psi}^{2}+x_{V}^{2}+f_{1}(u)x_{W}^{2}},$
(63) $\displaystyle\Omega_{\text{MG}}$ $\displaystyle\equiv$
$\displaystyle\frac{2\rho_{\text{MG}}}{M_{\text{Pl}}^{D-2}\lambda_{D}^{2}H^{2}}=x_{\psi}^{2}+x_{V}^{2}+f_{1}(u)x_{W}^{2},$
(64)

respectively. Then, the decelerated parameter has the form

$\displaystyle q$ $\displaystyle=$ $\displaystyle-1-\frac{\dot{H}}{H^{2}}$
(65) $\displaystyle=$
$\displaystyle-1+\frac{(D-1)}{2}\Big{[}1+w_{m}+(1-w_{m})x_{\psi}^{2}$
$\displaystyle-\,(1+w_{m})x_{V}^{2}-(f_{2}(u)+w_{m}f_{1}(u))x_{W}^{2}\Big{]},$

where for $q<0$ we have an accelerated universe model. There are five critical
points in the massless sector and three critical points in the massive sector,
which are listed in the Table 1 and 2, respectively, including their
properties. Note that the quantities $\mathcal{A}_{\pm}(D,w_{m},\lambda_{V})$
and $\mathcal{B}_{\pm}(D,w_{m},\lambda_{W})$ mentioned in these tables are
defined as

$\displaystyle\mathcal{A}_{\pm}(D,w_{m},\lambda_{V})\equiv\frac{1+w_{m}}{3-w_{m}}+\frac{(D-2)\lambda_{V}^{2}}{(D-2)(3-w_{m})}$
$\displaystyle\qquad\qquad\times\,\Bigg{\\{}1\pm\Bigg{[}1-\frac{2(D-1)(2-w_{m})(1+w_{m})}{(D-2)\lambda_{V}^{2}}$
$\displaystyle\qquad\qquad+\,\left(\frac{(D-1)(1+w_{m})}{(D-2)\lambda_{V}^{2}}\right)^{2}\Bigg{]}^{1/2}\Bigg{\\}}$
(66)
$\displaystyle\mathcal{B}_{\pm}(D,w_{m},\lambda_{W})\equiv\frac{1+w_{m}}{3-w_{m}}+\frac{(D-2)\lambda_{W}^{2}}{(D-1)(3-w_{m})}$
$\displaystyle\qquad\qquad\times\,\Bigg{\\{}1\pm\Bigg{[}1-\frac{2(D-1)(2-w_{m})(1+w_{m})}{(D-2)\lambda_{W}^{2}}$
$\displaystyle\qquad\qquad+\,\left(\frac{(D-1)(1+w_{m})}{(D-2)\lambda_{W}^{2}}\right)^{2}\Bigg{]}^{1/2}\Bigg{\\}}.$
(67)

Let us first discuss the massless sector in which there are five critical
points, namely, CP1, CP2, CP3, CP4 and CP5 with trivial coupling $W_{0}=0$.
The point CP1 describes a matter-dominated era in which the stability behavior
depends on the state parameter $w_{m}$. We find the late-time attractor
identified by a stable node and the past-time attractor identified by an
unstable node. There are three possible stabilities, namely stable node,
unstable node, and saddle point. The universe is either acceleratedly expanded
for $w_{m}<-{(D-3)}/{(D-1)}$ or unacceleratedly expanded for
$w_{m}\geq-{(D-3)}/{(D-1)}$. The first case is however unphysical since it
contains vacuum with $w_{m}=-1$, while in the latter case the matter of the
universe could be dominated by dust ($w_{m}=0$) or radiation
($w_{m}=\frac{1}{D-1}$).

Critical Points |  $x_{\psi,c}$ |  $x_{V,c}$ |  $x_{W,c}$ |  Existence |  $w_{\text{MG}}$ |  $\Omega_{\text{MG}}$ |  $q$ |  Stability
---|---|---|---|---|---|---|---|---
CP1 |  $0$ |  $0$ |  $0$ |  always |  undefined |  $0$ |  $\frac{D-3+w_{m}(D-1)}{2}$ |  stable node for $w_{m}<-1$, unstable node for $w_{m}>1$, saddle point otherwise
CP2 |  $1$ |  $0$ |  $0$ |  always |  $1$ |  $1$ |  $D-2$ |  stable node for $w_{m}>1$, $\lambda_{V}>2\sqrt{\frac{D-1}{D-2}}$, and $\lambda_{W}>2\sqrt{\frac{D-1}{D-2}}$, unstable node for $w_{m}<1$, $\lambda_{V}<2\sqrt{\frac{D-1}{D-2}}$, and $\lambda_{W}<2\sqrt{\frac{D-1}{D-2}}$, saddle point otherwise
CP3 |  $-1$ |  $0$ |  $0$ |  always |  $1$ |  $1$ |  $D-2$ |  unstable node for $w_{m}<1$, saddle point otherwise
CP4 |  $x_{\psi,c}$ |  $0$ |  $0$ |  $w_{m}=1$ and $0<|x_{\psi,c}|<1$ |  $1$ |  $x_{\psi,c}^{2}$ |  $D-2$ |  unstable for $-1<x_{\psi,c}<0$, or $0<x_{\psi,c}<1$ with at least either $\lambda_{V}<\frac{2}{x_{\psi,c}}\sqrt{\frac{D-1}{D-2}}$ or $\lambda_{W}<\frac{2}{x_{\psi,c}}\sqrt{\frac{D-1}{D-2}}$, nonhyperbolic for $0<x_{\psi,c}<1$
CP5 |  $\frac{\mathcal{A}_{\pm}}{\lambda_{V}}\sqrt{\frac{D-1}{D-2}}$ |  $\frac{1}{\lambda_{V}}\sqrt{\frac{(D-1)(2-\mathcal{A}_{\pm})\mathcal{A}_{\pm}}{2(D-2)}}$ |  $0$ |  $0<\mathcal{A}_{\pm}<2$ |  $\frac{3\mathcal{A}_{\pm}-2}{\mathcal{A}_{\pm}+2}$ |  $\frac{(D-1)(2+\mathcal{A}_{\pm})\mathcal{A}_{\pm}}{2(D-2)\lambda_{V}^{2}}$ |  $\frac{(D-1)\mathcal{A}_{\pm}}{2}-1$ |  see Fig. 1
Table 1:  The properties, the existence, the equation-of-state parameter $w_{\text{MG}}$, the density parameter $\Omega_{\text{MG}}$, the deceleration parameter $q$, and the stability conditions of the critical points of the autonomous system in the massless sector, $W(\psi)=0$. Note that we have introduced the notation $\mathcal{A}_{\pm}$ in Eq. (66). Critical Points |  $x_{\psi,c}$ |  $x_{V,c}$ |  $x_{W,c}$ |  Existence |  $w_{\text{MG}}$ |  $\Omega_{\text{MG}}$ |  $q$ |  Stability
---|---|---|---|---|---|---|---|---
CP6 |  $0$ |  $\sqrt{\frac{\lambda_{W}}{\lambda_{W}-\lambda_{V}}}$ |  $\sqrt{\frac{\lambda_{V}}{|f_{1}(u)|(\lambda_{W}-\lambda_{V})}}$ |  $\lambda_{W}>\lambda_{V}$ and $f_{1}(u)=f_{2}(u)<0$ |  $-1$ |  $1$ |  $-1$ |  stable node for $w_{m}>-1$ and $\lambda_{V}\lambda_{W}<\frac{1}{4}\frac{D-2}{D-1}$, stable spiral for $w_{m}>-1$ and $\lambda_{V}\lambda_{W}>\frac{1}{4}\frac{D-2}{D-1}$, saddle point otherwise
CP7 |  $0$ |  $\sqrt{1-f_{1}(u)x_{W,c}^{2}}$ |  $x_{W,c}$ |  $0\leq x_{W,c}\leq 1$, $\lambda_{V}=\lambda_{W}=0$, and $f_{1}(u)=f_{2}(u)>0$ |  $-1$ |  $1$ |  $-1$ |  unstable for $w_{m}<-1$, nonhyperbolic for $w_{m}>-1$
CP8 |  $\frac{\mathcal{B}_{\pm}}{\lambda_{W}}\sqrt{\frac{D-1}{D-2}}$ |  $0$ |  $\frac{1}{\lambda_{W}}\sqrt{\frac{(D-1)(2-\mathcal{B}_{\pm})\mathcal{B}_{\pm}}{2(D-2)f_{1}(u)}}$ |  $0<\mathcal{B}_{\pm}<2$ |  $\frac{3\mathcal{B}_{\pm}-2}{\mathcal{B}_{\pm}+2}$ |  $\frac{(D-1)(2+\mathcal{B}_{\pm})\mathcal{B}_{\pm}}{2(D-2)\lambda_{W}^{2}}$ |  $\frac{(D-1)\mathcal{B}_{\pm}}{2}-1$ |  see Fig. 2
Table 2:  The properties, the existence, the equation-of-state parameter
$w_{\text{MG}}$, the density parameter $\Omega_{\text{MG}}$, the deceleration
parameter $q$, and the stability conditions of the critical points of the
autonomous system in massive sector, $W(\psi)\neq 0$. Note that we have
introduced the notation $\mathcal{B}_{\pm}$ in Eq. (67).

At CP2 we also find a late-time and past-time attractors according to the
parameter values with also three possible stabilities, namely stable node,
unstable node, and saddle point. At this point, gravitons can be thought of as
nonphantom energy dominated by nonaccelerating expansion process of the
universe ($\Omega_{\text{MG}}=1$). The point CP3 has similar properties as CP2
in terms of its existence and the values of $w_{m}$, $\Omega_{\text{MG}}$, and
$q$. At this point we do not have a late-time attractor (a stable node).

The point CP4 exists on the interval $0<|x_{\psi,c}|<1$ which are unstable on
$-1<x_{\psi,c}<0$, or on $0<x_{\psi,c}<1$ with at least either
$\lambda_{W}<\frac{2}{x_{\psi,c}}\sqrt{\frac{D-1}{D-2}}$ or
$\lambda_{V}<\frac{2}{x_{\psi,c}}\sqrt{\frac{D-1}{D-2}}$. It becomes
nonhyperbolic on $0<x_{\psi,c}<1$ with
$\lambda_{V}>\frac{2}{x_{\psi,c}}\sqrt{\frac{D-1}{D-2}}$ and
$\lambda_{W}>\frac{2}{x_{\psi,c}}\sqrt{\frac{D-1}{D-2}}$. At this point, we
have a class of universe with nonaccelerating expansion since $D\geq 4$. These
universe are filled by the matter and the scalar field with constraint
$x_{\rho,c}^{2}+x_{\psi,c}^{2}=1$ and the equation of state parameters
$w_{m}=w_{\text{MG}}=1$.

Figure 1: The stability conditions of the critical point CP5 obtained by
plotting $\mathcal{A}_{\pm}$ defined in Eq. (66) as a function of
$\lambda_{V}$. The forbidden zone is the area with unphysical properties
$\Omega_{\text{MG}}>1$. This figure consist of two possibilities of the
parameter values $\lambda_{V}$ and $\lambda_{W}$, where (a) for
$\lambda_{V}>\lambda_{W}$, and (b) for $\lambda_{V}<\lambda_{W}$.

The point CP5 exists on the interval $0<\mathcal{A}_{\pm}<2$ which could be
either the late-time attractor, the past-time attractor, or saddle points,
depending on the parameter $\lambda_{V}$ dan $\lambda_{W}$. For
$\mathcal{A}_{\pm}(D,\lambda_{V},\lambda_{W})<\frac{2}{D-1}$, we have a class
of universe which acceleratedly expands and the scalar field plays a role as
quintessencelike. The density parameter are on the interval
$0\leq\Omega_{\text{MG}}\leq 1$ such that
$\lambda_{V}<\sqrt{\frac{(D-1)(2+\mathcal{A}_{\pm})\mathcal{A}_{\pm}}{2(D-2)}}$
is not allowed since it will produce $\Omega_{\text{MG}}>1$. Around
$\lambda_{V}\approx\sqrt{\frac{(D-1)(2+\mathcal{A}_{\pm})\mathcal{A}_{\pm}}{2(D-2)}}$,
the expansion of the universe is dominated by the scalar field
($\Omega_{\text{MG}}\approx 1$) which will be a good candidate for saddle
power-law inflation model. This will be discussed in detail in Sec. VII.

Next, we consider the massive sector which consists of the points CP6, CP7,
and CP8. The points CP6 and CP7 have vanishing kinetic part of the scalar
field, which implies that the scalar becomes very massive. On the other hand,
the kinetic part of the scalar is nonzero at CP8.

The point CP6 could be either the late-time attractor for $w_{m}>-1$ or saddle
point for $w_{m}<-1$ where the graviton mass plays a role as the cosmological
constant which dominates the accelerating expansion of universe. This
cosmological constant has to be positive since the parameter
$\lambda_{W}>\lambda_{V}$. These features give us a good candidate for the
compatible description of the well-known observation result [48].

Figure 2: The stability conditions of the critical point CP8 obtained by
plotting $\mathcal{B}_{\pm}$ defined in Eq. (67) as a function of
$\lambda_{W}$. The forbidden zone is the area with unphysical properties
$\Omega_{\text{MG}}>1$. This figure consists of two possibilities of the
parameter values $\lambda_{V}$ and $\lambda_{W}$, where (a) for
$\lambda_{V}<\lambda_{W}$, and (b) for $\lambda_{V}>\lambda_{W}$.

The point CP7 is defined on the circle $x_{V,c}^{2}+f_{1}(u)x_{W,c}^{2}=1$
which may be either unstable for $w_{m}<-1$ or nonhyperbolic for $w_{m}>-1$.
At this point, the parameter $\lambda_{V}=\lambda_{W}=0$ showing that both the
scalar potential and the graviton mass are constant. They behave like the
cosmological constant which dominates the accelerating expansion of universe.

Finally, the point CP8 which could be either the late-time attractor or the
saddle point depending on the values of $\lambda_{V}$ and $\lambda_{W}$. For
$\mathcal{B}_{\pm}(D,\lambda_{V},\lambda_{W})<\frac{2}{D-1}$, the universe
acceleratedly expands and the mass-varying massive graviton plays a role as
quentessence dark energy. There exists a forbidden zone in CP8; it has
$\Omega_{\text{MG}}>1$ for
$\lambda_{W}<\sqrt{\frac{(D-1)(2+\mathcal{B}_{\pm})\mathcal{B}_{\pm}}{2(D-2)}}$.

##  VI Local-Global Existence of Solutions

In this section we will prove the local-global existence and the uniqueness of
the evolution equations (57), (58), and (59) with constraint (60) using
Picard’s iteration and the contraction mapping properties. We will first
discuss the $f_{1}(u)>0$ case, then continue with the $f_{1}(u)<0$ case.

First of all, we introduce the dynamical variables

$\bm{u}=\begin{pmatrix}x_{\psi}\\\ x_{V}\\\ x_{W}\end{pmatrix},$ (68)

defined on an interval $I\equiv[s,s+\epsilon]$ where
$s\equiv\ln{a}\in\mathbb{R}$ and $\epsilon$ is a small positive constant. The
functions $f_{1}(u)$ and $f_{2}(u)$ in Eqs. (61) and (62)are bounded. In the
first part of this section, we consider the case of $f_{1}(u)>0$ such that
from the constraint (60) we could have

$\left.\begin{array}[]{cc}0\leq|x_{\psi}|\leq 1&\\\ 0\leq|x_{V}|\leq 1&\\\
0\leq|f_{1}(u)|^{1/2}|x_{W}|\leq 1&\end{array}\right\\}$ (69)

In other words, all of the quantities $(x_{\psi},x_{V},f_{1}(u)x_{W})$ are
defined on an open set $U\subset S^{3}$ where $S^{3}$ is the 3-sphere. It is
important to notice that the critical point CP6 is excluded in this setup.

All of the evolution equations (57), (58), and (59) can be simply rewritten
into

$\frac{d\bm{u}}{ds}=\mathcal{J}(\bm{u}),$ (70)

with

$\displaystyle\mathcal{J}(\bm{u})\equiv\frac{1}{2}(D-2)\left(\begin{array}[]{c}(1-w_{m})x_{\psi}^{3}-(1+w_{m})x_{\psi}x_{V}^{2}-(f_{2}(u)+w_{m}f_{1}(u))x_{\psi}x_{W}^{2}\\\
-(1-w_{m})x_{\psi}+2\sqrt{\frac{D-2}{D-1}}\left(\lambda_{V}x_{V}^{2}+\lambda_{W}f_{1}(u)x_{W}^{2}\right)\\\
\\\ (1+w_{m})x_{V}+(1-w_{m})x_{\psi}^{2}x_{V}-(1+w_{m})x_{V}^{3}\\\
-(f_{2}(u)+w_{m}f_{1}(u))x_{V}x_{W}^{2}-\sqrt{\frac{D-2}{D-1}}\lambda_{V}x_{\psi}x_{V}\\\
\\\ (1+w_{m})x_{W}+(1-w_{m})x_{\psi}^{2}x_{W}-(1+w_{m})x_{V}^{2}x_{W}\\\
-(f_{2}(u)+w_{m}f_{1}(u))x_{W}^{3}-\sqrt{\frac{D-2}{D-1}}\lambda_{W}x_{\psi}x_{W}\end{array}\right).$
(79)

######  Lemma VI.1.

The operator $\mathcal{J}(\bm{u})$ in Eq. (70) is locally Lipschitz with
respect to $\bm{u}$.

###### Proof.

We have the following estimate

$\displaystyle|\mathcal{J}|_{U}\leq\frac{1}{2}(D-2)\Bigg{[}|1-w_{m}||x_{\psi}|^{3}+|1+w_{m}||x_{\psi}||x_{V}|^{2}$
$\displaystyle\qquad+\,|f_{2}(u)+w_{m}f_{1}(u)||x_{\psi}||x_{W}|^{2}+|1-w_{m}||x_{\psi}|$
$\displaystyle\qquad+\,2\sqrt{\frac{D-2}{D-1}}\left(\lambda_{V}|x_{V}|^{2}+\lambda_{W}|f_{1}(u)||x_{W}|^{2}\right)$
$\displaystyle\qquad+\,|1+w_{m}||x_{V}|+|1-w_{m}||x_{\psi}|^{2}|x_{V}|+|1+w_{m}||x_{V}|^{3}$
$\displaystyle\qquad+\,|f_{2}(u)+w_{m}f_{1}(u)||x_{V}||x_{W}|^{2}+\sqrt{\frac{D-2}{D-1}}\lambda_{V}|x_{\psi}||x_{V}|$
$\displaystyle\qquad+\,|1+w_{m}||x_{W}|+|1-w_{m}||x_{\psi}|^{2}|x_{W}|$
$\displaystyle\qquad+\,|1+w_{m}||x_{V}|^{2}|x_{W}|+|f_{2}(u)+w_{m}f_{1}(u)||x_{W}|^{3}$
$\displaystyle\qquad+\,\sqrt{\frac{D-2}{D-1}}\lambda_{W}|x_{\psi}||x_{W}|\Bigg{]}.$
(80)

Then, using Eq. (69), we can show that $|\mathcal{J}(\bm{u},x)|_{U}$ is indeed
bounded on $U$.

Moreover, for $\bm{u},\hat{\bm{u}}\in U$ we have

$\displaystyle|\mathcal{J}(\bm{u})-\mathcal{J}(\hat{\bm{u}})|_{U}\leq\frac{1}{2}(D-2)\Bigg{[}|1-w_{m}||x_{\psi}^{3}-\hat{x}_{\psi}^{3}|$
$\displaystyle\qquad+\,|1+w_{m}||x_{\psi}x^{2}_{V}-\hat{x}_{\psi}\hat{x}^{2}_{V}|+|1-w_{m}||x_{\psi}-\hat{x}_{\psi}|$
$\displaystyle\qquad+\,|f_{2}(u)+w_{m}f_{1}(u)||x_{\psi}x^{2}_{W}-\hat{x}_{\psi}\hat{x}^{2}_{W}|$
$\displaystyle\qquad+\,2\sqrt{\frac{D-2}{D-1}}\left(\lambda_{V}|x^{2}_{V}-\hat{x}^{2}_{V}|+\lambda_{W}|f_{1}(u)||x^{2}_{W}-\hat{x}^{2}_{W}|\right)$
$\displaystyle\qquad+\,|1+w_{m}||x_{V}-\hat{x}_{V}|+|1-w_{m}||x^{2}_{\psi}x_{V}-\hat{x}^{2}_{\psi}\hat{x}_{V}|$
$\displaystyle\qquad+\,|1+w_{m}||x^{3}_{V}-\hat{x}^{3}_{V}|$
$\displaystyle\qquad+\,|f_{2}(u)+w_{m}f_{1}(u)||x_{V}x^{2}_{W}-\hat{x}_{V}\hat{x}^{2}_{W}|$
$\displaystyle\qquad+\,\sqrt{\frac{D-2}{D-1}}\lambda_{V}|x_{\psi}x_{V}-\hat{x}_{\psi}\hat{x}_{V}|$
$\displaystyle\qquad+\,|1+w_{m}||x_{W}-\hat{x}_{W}|+|1-w_{m}||x_{\psi}^{2}x_{W}-\hat{x}_{\psi}^{2}\hat{x}_{W}|$
$\displaystyle\qquad+\,|1+w_{m}||x_{V}^{2}x_{W}-\hat{x}_{V}^{2}\hat{x}_{W}|$
$\displaystyle\qquad+\,|f_{2}(u)+w_{m}f_{1}(u)||x^{3}_{W}-\hat{x}^{3}_{W}|$
$\displaystyle\qquad+\,\sqrt{\frac{D-2}{D-1}}\lambda_{W}|x_{\psi}x_{W}-\hat{x}_{\psi}\hat{x}_{W}|\Bigg{]}.$
(81)

After some computations, we obtain

$\left|\mathcal{J}(\bm{u})-\mathcal{J}(\hat{\bm{u}})\right|_{U}\leq
C_{\mathcal{J}}(|\bm{u}|,|\hat{\bm{u}}|)|\bm{u}-\hat{\bm{u}}|,$ (82)

showing that $\mathcal{J}$ is locally Lipshitz with respect to $\bm{u}$. ∎

Next, we rewrite Eq. (70) into the integral form

$\bm{u}(s)=\bm{u}(s_{0})+\int_{s_{0}}^{s}\,\mathcal{J}\left(\bm{u}(\hat{s})\right)d\hat{s}.$
(83)

We define a Banach space

$X\equiv\\{\bm{u}\in
C(I,\mathbb{R}^{2}):\,\bm{u}(x_{0})=\bm{u}_{0},\,\sup_{x\in
I}{|\bm{u}(x)|}\leq L_{0}\\},$ (84)

endowed with the norm

$|\bm{u}|_{X}=\sup_{x\in I}{|\bm{u}(x)|},$ (85)

where $L_{0}>0$. Introducing an operator $\mathcal{K}$

$\mathcal{K}(\bm{u}(x))=\bm{u}_{0}+\int_{x_{0}}^{x}\mathcal{J}\left(\bm{u}(s),s\right)ds,$
(86)

and using Lemma VI.1, we have the following result [49]:

######  Lemma VI.2.

Let $\mathcal{K}$ be an operator defined in Eq. (86). Suppose there exists a
constant $\varepsilon>0$ such that $\mathcal{K}$ is a mapping from $X$ to
itself and $\mathcal{K}$ is a contraction mapping on $I=[x,x+\varepsilon]$
with

$\varepsilon\leq\min\left(\frac{1}{C_{L_{0}}},\frac{1}{C_{L_{0}}L_{0}+\|\mathcal{J}(x)\|}\right).$
(87)

Then, the operator $\mathcal{K}$ is a contraction mapping on $X$.

The above lemma shows that there exists a unique fixed point of Eq. (86)
ensuring a unique local solution of the differential equation (70). We can
further construct a maximal solution by repeating the above arguments of the
local existence with the initial condition $\bm{u}(x-x_{n})$ for some
$x_{0}<x_{n}<x$ and using the uniqueness condition to glue the solutions.

We can now show the existence of global solutions of Eq. (70). Let us consider
the integral form (83) such that

$|\bm{u}(s)|\leq|\bm{u}(s_{0})|+\int_{s_{0}}^{s}|\mathcal{J}\left(\bm{u}(\hat{s})\right)|d\hat{s}.$
(88)

We first consider the self-accelerating branch where the parameter $u$ in Eq.
(40) is constant. Using Eqs. (69) and (80), we get

$\displaystyle|\bm{u}(t)|\leq|\bm{u}(t_{0})|+\frac{1}{2}(D-2)\Bigg{[}3|1-w_{m}|+3|1+w_{m}|$
$\displaystyle\qquad+\,|f_{2}(u)+w_{m}|+\sqrt{\frac{D-2}{D-1}}\left(3\lambda_{V}+2\lambda_{W}\right)$
$\displaystyle\qquad+\,2\left|\frac{f_{2}(u)}{f_{1}(u)}+w_{m}\right|+2\frac{|1+w_{m}|}{|f_{1}(u)|^{1/2}}+\frac{|1-w_{m}|}{|f_{1}(u)|^{1/2}}$
$\displaystyle\qquad+\,\frac{|f_{2}(u)+w_{m}f_{1}(u)|}{|f_{1}(u)|^{3/2}}$
$\displaystyle\qquad+\,\sqrt{\frac{D-2}{D-1}}\frac{\lambda_{W}}{|f_{1}(u)|^{1/2}}\Bigg{]}\ln{\left(\frac{a(t)}{a(t_{0})}\right)}.$
(89)

The second part is to consider the $f_{1}(u)<0$ case. From the constraint (60)
we could have

$\left.\begin{array}[]{rcl}x_{\psi}&=&\cos{\alpha}\\\
x_{V}&=&\sin{\alpha}\cosh{\beta}\\\
|f_{1}(u)|^{1/2}x_{W}&=&\sin{\alpha}\sinh{\beta}\end{array}\right\\}$ (90)

where $\alpha\equiv\alpha(s)$ and $\beta\equiv\beta(s)$. In this case, Lemma
VI.1 and Lemma VI.2 still hold, but we need to modify the estimate to show the
global existence. In the case at hand, using Eq. (90), the estimate (88)
becomes

$\displaystyle|\bm{u}(t)|\leq|\bm{u}(t_{0})|+\frac{1}{2}(D-2)\Bigg{\\{}2|1-w_{m}|\ln{\left(\frac{a(t)}{a(t_{0})}\right)}$
$\displaystyle\qquad+\,\Bigg{[}|1+w_{m}|+|1-w_{m}|+\frac{|1+w_{m}|}{|f_{1}(u)|^{1/2}}+\frac{|1-w_{m}|}{|f_{1}(u)|^{1/2}}$
$\displaystyle\qquad+\,\sqrt{\frac{D-2}{D-1}}\left(\lambda_{V}+\frac{\lambda_{W}}{|f_{1}(u)|^{1/2}}\right)\Bigg{]}\int_{s_{0}}^{s}\cosh{\beta}\,d\hat{s}$
$\displaystyle\qquad+\,\Bigg{[}|1+w_{m}|+\left|\frac{f_{2}(u)}{f_{1}(u)}+w_{m}\right|$
$\displaystyle\qquad+\,2\sqrt{\frac{D-2}{D-1}}\left(\lambda_{V}+\lambda_{W}\right)\Bigg{]}\int_{s_{0}}^{s}\cosh^{2}{\beta}\,d\hat{s}$
$\displaystyle\qquad+\,\Bigg{[}2|1+w_{m}|+\left|\frac{f_{2}(u)}{f_{1}(u)}+w_{m}\right|$
$\displaystyle\qquad+\,\frac{|f_{2}(u)+w_{m}f_{1}(u)|}{|f_{1}(u)|^{3/2}}\Bigg{]}\int_{s_{0}}^{s}\cosh^{3}{\beta}\,d\hat{s}\Bigg{\\}}.$
(91)

For the normal branch with $u(t)$, we employ similar procedure as above and
use the assumptions that $f_{1}(u)$ and $f_{2}(u)$ are bounded. Then, we
obtain the slightly modified forms of Eqs. (89) and (91).

Thus, we have proven

######  Theorem VI.1.

There exists a global solution of the evolution equations (57), (58), and (59)
with constraint (60).

##  VII Cosmological Models

In this section we will discuss some possible cosmological models of the
theory. To simplify the computation, we particularly choose the self-
accelerating branch where the parameter $u$ in Eq. (40) is constant. In the
case of the exponential form of the potentials (51) and (52), we may have an
inflation era in which it can be described by the well-known power law
inflation [50] where $a(t)\varpropto t^{1/\epsilon}$, with the slow-roll
parameter $\epsilon={|\dot{H}|}/{H^{2}}<1$.

If the scalar $\psi$ plays a role as the inflaton field in the early epoch,
then the critical points CP5 and CP8 are the good candidates to describe that
era, with the slow-roll parameter given by

$\displaystyle\epsilon$ $\displaystyle=$
$\displaystyle\frac{(D-2)\lambda_{V}^{2}}{4}<1,$ (92) $\displaystyle\epsilon$
$\displaystyle=$
$\displaystyle\frac{(D-2)\lambda_{W}^{2}}{4}\left(1+\frac{\sum_{n=0}^{D-1}c_{n}A_{n}u^{D-n-1}}{\sum_{n=0}^{D-1}c_{n+1}A_{n}u^{D-n-1}}\right)<1,$

for CP5 and CP8, respectively. In the CP5 case, which is in the massless
sector, we can use Eqs. (64) and (92) to get

$w_{\text{MG}}=\frac{D-2}{D-1}\frac{\lambda_{V}^{2}}{2}-1,$ (94)

which is negative for $D>3$. Hence, from Table 1, we have

$\mathcal{A}_{\pm}=\frac{2(D-2)\lambda_{V}^{2}}{8(D-1)-(D-2)\lambda_{V}^{2}}<\frac{1}{D-1}.$
(95)

Similarly, in the CP8 case, which is in the massive sector, assuming that $u$
is constant, we can use Eqs. (64) and (LABEL:eq:epsiloncp8) to get

$w_{\text{MG}}=\frac{D-2}{D-1}\frac{\lambda_{W}^{2}}{2}\left(1+\frac{\sum_{n=0}^{D-1}c_{n}A_{n}u^{D-n-1}}{\sum_{n=0}^{D-1}c_{n+1}A_{n}u^{D-n-1}}\right)-1,$
(96)

which is negative for $D>3$. Hence, from Table 2, we obtain

$\mathcal{B}_{\pm}=\frac{2(D-2)\lambda_{W}^{2}\sum_{n=0}^{D-1}(c_{n}+c_{n+1})A_{n}u^{D-n-1}}{8(D-1)\sum_{n=0}^{D-1}c_{n+1}A_{n}u^{D-n-1}-(D-2)\lambda_{W}^{2}\sum_{n=0}^{D-1}(c_{n}+c_{n+1})A_{n}u^{D-n-1}}<\frac{1}{D-1}.$
(97)

The scalar $\psi$ in the inflation era has the form

$\psi(t)\propto\frac{\sqrt{M_{\text{Pl}}^{D-2}}}{\lambda_{\alpha}}\ln{\left[\frac{\lambda_{\alpha}^{2}V_{0}\epsilon
t^{2}}{2M_{\text{Pl}}^{D-2}(D-1-\epsilon)}\right]},$ (98)

where $\lambda_{\alpha}=\lambda_{V}$ ($\lambda_{W}$) for the case of CP5
(CP8), respectively.

In order to have a nucleosynthesis era, we have to assume $\epsilon\ll 1$ such
that the reheating phase begins at $t_{i}$ and ends at $t_{f}$. Therefore, we
have the scale factor

$\frac{a(t_{f})}{a(t_{i})}=\left(\frac{t_{f}}{t_{i}}\right)^{1/\epsilon}=\exp{\left[\frac{\lambda_{\alpha}(\psi_{f}-\psi_{i})}{2\epsilon\sqrt{M_{\text{Pl}}^{D-2}}}\right]}.$
(99)

In addition, the scalar potential in $\psi_{i}=\psi(t_{i})$ drops to the
scalar potential in $\psi_{f}=\psi(t_{f})$ implying that the argument of the
exponential is positive and the scale factor increases during this period. At
the end of the inflation (or in the beginning of the reheating at $t=t_{i}$),
the scalar field $\psi$ still dominates the universe with initial density
parameter $\Omega_{\text{MG, i}}\simeq 1$. Then, the phase transition of
$\psi$ occurs such that $\psi$ decays to matter fields in the time interval
$t_{i}<t<t_{f}$. The decay process of the scalar $\psi$ becoming a matter
field occurs relativistically with $P_{m^{\ast}}={\rho_{m^{\ast}}}/{(D-1)}$.
From Eq. (99) and the conservation of energy of relativistic matter
(radiation), we obtain the relation
$\rho_{\text{MG},f}/{\rho_{\text{MG},i}}=({a_{i}}/{a_{f}})^{2\epsilon}$ and
$\rho_{M,f}/\rho_{M,i}=({a_{i}}/{a_{f}})^{D}$, respectively. Since the
inflaton density parameter
$\Omega_{\text{MG}}=2\rho_{\text{MG}}/{(M_{\text{Pl}}^{D-2}\lambda_{D}^{2}H^{2})}$
and the matter field density, we find that $\Omega_{\text{MG}}$ decreases
according to

$\Omega_{\text{MG}}\simeq\left(\frac{1\
\text{MeV}^{4}}{\rho_{M}}\right)^{2\epsilon/D},\qquad(t_{i}\leq t\leq t_{f}),$
(100)

as $\rho_{M}$ increases. Note that here we have assumed that
$\epsilon={|\dot{H}|}/{H^{2}}\ll 1$ in the interval. At the end, we have
$\rho_{M}\sim 10^{4}-10^{8}\text{ MeV}^{4}$ at $t=t_{f}$.

On the other hand, we may also apply the theory to the case of late-time era.
Here, we have at least two interesting possible scenarios. The first is that
the scalar field $\psi$ becomes frozen, after the reheating era, at fixed
value $\psi_{\infty}$ with $w_{\text{MG}}=-1$ at CP6. The value of
$W(\psi_{\infty})$ then becomes the graviton mass in the present time, which
can be determined using observational constraints. Hence, the accelerating
expansion of the universe is due to the constant mass term, as in the dRGT
theory. The second possible scenario is that the scalar field $\psi$ could
play a role as dynamical quintessensial dark energy, either from CP5 in the
massless sector or from CP8 in the massive sector. In the massless sector, the
accelerating expansion of the universe is due solely to the standard
quintessence paradigm. However, in the massive sector, it is due to the
nontrivial interplay between quintessence and massive gravity in the massive
sector.

##  VIII Conclusions

We have constructed higher-dimensional MTMVMG with nonzero scalar potential,
where the graviton mass is varied with respect to the real scalar field
$\psi$. Our construction can be summarized as follows. First, we write the
MVMG action for higher dimensions using the vielbein formulation in the ADM
formalism. We also adopt the vielbein potential in Ref. [35] and couple it to
the masslike scalar potential $W(\psi)$. Inserting the ansatz metric (3) and
(4) into the MVMG action, we then derive the precursor action (9). By imposing
the set of $D$-constraints (19) and (20) to the precursor action (9), we
obtain a theory in which the graviton has $D(D-3)/2$ degrees of freedom for
general $D$ spacetime dimensions, without scalar and vector modes of the
Stückelberg field. The resulting theory admits the Lorentz violation since we
have used the ADM vielbein (3) and (4). Still, the $O(D-1)$ symmetry and the
spatial diffeomorphism are preserved in the MTMVMG action (36). This theory
generalizes the four-dimensional MTMG in Refs. [29, 30] and the four-
dimensional MVMG in Ref. [32].

Then, we derive the Friedmann-Lemaître equations (43) and (46) for the case of
spatially flat spacetimes. To proceed, inspired by Refs. [46, 44, 45], we take
both the scalar potential and the graviton mass couplings to have exponential
forms (51) and (52), such that Eqs. (43) and (46) can be written into a set of
autonomous equations (59) with constraint (60). By performing dynamical system
analysis, we find that in this theory there exists five critical points in the
massless sector, namely CP1-5, whereas in the massive sector we have three
critical points, namely CP6-8. We also discuss their stability, their
existence, and their cosmological aspects related to the state equation
parameter $w_{\text{MG}}$, the density parameter $\Omega_{\text{MG}}$, and the
decelerated parameter $q$. Among them, we may have some critical points that
are suitable to explain either inflation phenomenon or the accelerated
universe in the late-time era.

We also have established the local-global existence and the uniqueness of the
evolution equations (57), (58), and (59) with constraint (60) using Picard’s
iteration and the contraction mapping properties, assuming that the functions
$f_{1}(u)$ and $f_{2}(u)$ are bounded. The discussion is then divided into two
parts: the $f_{1}(u)>0$ case and the $f_{1}(u)<0$ case. Note that our results
apply to all branches, namely the self-accelerating branch and the normal
branch.

Finally, we have particularly discussed some possible cosmological models of
the MTMVMG in the self-accelerating branch. Since both the scalar potential
and the graviton mass couplings have exponential forms (51) and (52), the
theory has a good description of the inflation era in the early universe using
the power-law inflation [50] in which the scale factor $a(t)\varpropto
t^{1/\epsilon}$ with the slow-roll parameter $\epsilon={|\dot{H}|}/{H^{2}}<1$.
This era can be described either by the critical point CP5 or CP8. In other
words, our theory can describe the inflationary era using both the massless
and the massive sectors. Also, we have shown that the MTMVMG could accommodate
the reheating mechanism in this framework, again for both massless and massive
sectors. Perturbative approach needs to be applied, for example, to study the
behavior of primordial gravitational waves based on MTMVMG, as in the case of
four-dimensional MTMG [51, 52]. The detailed construction and the
phenomenological predictions are left for subsequent works. On the other hand,
we have at least two interesting possible scenarios for the late times. The
first scenario is that the dark energy in the present time is due to the
graviton mass which depends on the scalar field $\psi_{\infty}$ that becomes
frozen after the reheating era. The second scenario is that the scalar field
$\psi$ plays role as dynamical quintessential dark energy. Therefore, contrary
to the massless sector where the accelerating expansion is due to the standard
quintessence paradigm, in the massive sector it is due to the nontrivial
interplay between quintessence and massive gravity.

## Acknowledgments

The work in this paper is supported by P2MI FMIPA ITB 2021 and Riset ITB 2021.
The work of H. A. is partially funded by the World Class Research (WCR) grant
from RistekBRIN-IPB 2021.

##  Appendix A Mass Term

In this Appendix we are going to evaluate potential element in Eq. (7) by
taking the ADM vielbein in Eqs. (3) and (4). For $n=0$ we obtain

$\displaystyle\frac{1}{D!}\hat{\epsilon}_{A_{1}A_{2}\cdots
A_{D}}\tensor{\bm{\tilde{E}}}{{}^{A_{1}}}\wedge\tensor{\bm{\tilde{E}}}{{}^{A_{2}}}\wedge\tensor{\bm{\tilde{E}}}{{}^{A_{3}}}\wedge\dots\wedge\tensor{\bm{\tilde{E}}}{{}^{A_{D}}}$
$\displaystyle\quad=\frac{d^{D}x}{D!}\hat{\epsilon}_{A_{1}A_{2}\cdots
A_{D}}\hat{\epsilon}^{\mu_{1}\mu_{2}\cdots\mu_{D}}\tensor{\tilde{E}}{{}^{A_{1}}_{\mu_{1}}}\tensor{\tilde{E}}{{}^{A_{2}}_{\mu_{2}}}\tensor{\tilde{E}}{{}^{A_{3}}_{\mu_{3}}}\cdots\tensor{\tilde{E}}{{}^{A_{D}}_{\mu_{D}}}$
$\displaystyle\quad=\frac{d^{D}x}{D!}\hat{\epsilon}_{A_{1}A_{2}\cdots
A_{D}}\hat{\epsilon}^{B_{1}B_{2}\cdots
B_{D}}\tensor{\delta}{{}^{A_{1}}_{B_{1}}}\tensor{\delta}{{}^{A_{2}}_{B_{2}}}\cdots\tensor{\delta}{{}^{A_{D}}_{B_{D}}}\det{(\tilde{E})}$
$\displaystyle\quad=d^{D}x\,M\det{(\tilde{e})}.$ (101)

Note that we have used the relations $dx^{\mu_{1}}\wedge
dx^{\mu_{2}}\wedge\cdots\wedge
dx^{\mu_{D}}=\hat{\epsilon}^{\mu_{1}\mu_{2}\cdots\mu_{D}}d^{D}x$ in first step
and
$\hat{\epsilon}^{\mu_{1}\mu_{2}\cdots\mu_{D}}=\hat{\epsilon}^{B_{1}B_{2}\cdots
B_{D}}\tensor{\tilde{E}}{{}^{\mu_{1}}_{B_{1}}}\tensor{\tilde{E}}{{}^{\mu_{2}}_{B_{2}}}\cdots\tensor{\tilde{E}}{{}^{\mu_{D}}_{B_{D}}}$
in the second step. Following the same way as above for other $n$, we have

$\displaystyle\hat{\epsilon}_{A_{1}A_{2}\cdots
A_{D}}\tensor{\bm{E}}{{}^{A_{1}}}\wedge\tensor{\bm{\tilde{E}}}{{}^{A_{2}}}\wedge\tensor{\bm{\tilde{E}}}{{}^{A_{3}}}\wedge\dots\wedge\tensor{\bm{\tilde{E}}}{{}^{A_{D}}}$
$\displaystyle=$
$\displaystyle(D-1)!\Big{(}M\det{(\tilde{e})}[\tilde{e}^{-1}e]+N\det{(\tilde{e})}\Big{)},$
$\displaystyle\hat{\epsilon}_{A_{1}A_{2}\cdots
A_{D}}\tensor{\bm{E}}{{}^{A_{1}}}\wedge\tensor{\bm{E}}{{}^{A_{2}}}\wedge\tensor{\bm{\tilde{E}}}{{}^{A_{3}}}\wedge\dots\wedge\tensor{\bm{\tilde{E}}}{{}^{A_{D}}}$
$\displaystyle=$
$\displaystyle(D-2)!\Big{(}M\det{(\tilde{e})}([\tilde{e}^{-1}e]^{2}-[(\tilde{e}^{-1}e)^{2}])+2N\det{(\tilde{e})}[\tilde{e}^{-1}e]\Big{)},$
$\displaystyle\hat{\epsilon}_{A_{1}A_{2}\cdots
A_{D}}\tensor{\bm{E}}{{}^{A_{1}}}\wedge\tensor{\bm{E}}{{}^{A_{2}}}\wedge\tensor{\bm{E}}{{}^{A_{3}}}\wedge\dots\wedge\tensor{\bm{\tilde{E}}}{{}^{A_{D}}}$
$\displaystyle=$
$\displaystyle(D-3)!\Big{(}M\det{(\tilde{e})}\big{(}[\tilde{e}^{-1}e]^{3}-3[\tilde{e}^{-1}e][(\tilde{e}^{-1}e)^{2}]$
$\displaystyle+\,2[(\tilde{e}^{-1}e)^{3}]\big{)}+3N\det{(\tilde{e})}([\tilde{e}^{-1}e]^{2}-[(\tilde{e}^{-1}e)^{2}])\Big{)},$
$\displaystyle\vdots$ $\displaystyle\hat{\epsilon}_{A_{1}A_{2}\cdots
A_{D}}\tensor{\bm{\tilde{E}}}{{}^{A_{1}}}\wedge\tensor{\bm{\tilde{E}}}{{}^{A_{2}}}\wedge\tensor{\bm{E}}{{}^{A_{3}}}\wedge\dots\wedge\tensor{\bm{E}}{{}^{A_{D}}}$
$\displaystyle=$
$\displaystyle(D-2)!\Big{(}N\det{(e)}\left([e^{-1}\tilde{e}]^{2}-[(e^{-1}\tilde{e})^{2}]\right)+2M\det{(e)}[e^{-1}\tilde{e}]\Big{)},$
$\displaystyle\hat{\epsilon}_{A_{1}A_{2}\cdots
A_{D}}\tensor{\bm{\tilde{E}}}{{}^{A_{1}}}\wedge\tensor{\bm{E}}{{}^{A_{2}}}\wedge\tensor{\bm{E}}{{}^{A_{3}}}\wedge\dots\wedge\tensor{\bm{E}}{{}^{A_{D}}}$
$\displaystyle=$
$\displaystyle(D-1)!\left(N\det{(e)}[e^{-1}\tilde{e}]+M\det{(e)}\right),$
$\displaystyle\hat{\epsilon}_{A_{1}A_{2}\cdots
A_{D}}\tensor{\bm{E}}{{}^{A_{1}}}\wedge\tensor{\bm{E}}{{}^{A_{2}}}\wedge\tensor{\bm{E}}{{}^{A_{3}}}\wedge\dots\wedge\tensor{\bm{E}}{{}^{A_{D}}}$
$\displaystyle=$ $\displaystyle D!\,N\det{(e)},$

where $[\cdots]$ denotes the trace. Now we define
$\tensor{X}{{}^{I}_{J}}\equiv\tensor{\tilde{e}}{{}^{I}_{k}}\tensor{e}{{}^{k}_{J}}$
and
$\tensor{Y}{{}^{I}_{J}}\equiv\tensor{e}{{}^{I}_{k}}\tensor{\tilde{e}}{{}^{k}_{J}}$
such that

$\displaystyle\sum_{n=0}^{D}\frac{c_{n}}{n!(D-n)!}\hat{\epsilon}_{A_{1}A_{2}\cdots
A_{D}}\tensor{\bm{E}}{{}^{A_{1}}}\wedge\dots\wedge\tensor{\bm{E}}{{}^{A_{n}}}\wedge\tensor{\bm{\tilde{E}}}{{}^{A_{n+1}}}\wedge\dots\wedge\tensor{\bm{\tilde{E}}}{{}^{A_{D}}}$
$\displaystyle\qquad\qquad=d^{D}x\,N\det{(e)}\bigg{(}|\det{(X)}|\frac{M}{N}\sum_{n=0}^{D-1}c_{n}\mathcal{S}_{n}(Y)+\sum_{n=0}^{D-1}c_{D-n}\mathcal{S}_{n}(X)\bigg{)}.$
(102)

with $\mathcal{S}_{n}$ symmetric polynomial. For example given matrix
$\mathrm{A}$ size $D\times D$,

$\displaystyle\mathcal{S}_{0}(\mathrm{A})$ $\displaystyle=$ $\displaystyle 1,$
$\displaystyle\mathcal{S}_{1}(\mathrm{A})$ $\displaystyle=$
$\displaystyle[\mathrm{A}],$ $\displaystyle\mathcal{S}_{2}(\mathrm{A})$
$\displaystyle=$
$\displaystyle\frac{1}{2!}\Big{(}[\mathrm{A}]^{2}-[\mathrm{A}^{2}]\Big{)},$
$\displaystyle\mathcal{S}_{3}(\mathrm{A})$ $\displaystyle=$
$\displaystyle\frac{1}{3!}\Big{(}[\mathrm{A}]^{3}-3[\mathrm{A}][\mathrm{A}^{2}]+2[\mathrm{A}^{3}]\Big{)},$
$\displaystyle\mathcal{S}_{4}(\mathrm{A})$ $\displaystyle=$
$\displaystyle\frac{1}{4!}\Big{(}[\mathrm{A}]^{4}-6[\mathrm{A}]^{2}[\mathrm{A}^{2}]+8[\mathrm{A}][\mathrm{A}^{3}]$
$\displaystyle+\,3[\mathrm{A}^{2}]^{2}-6[\mathrm{A}^{4}]\Big{)}$
$\displaystyle\vdots$ $\displaystyle\mathcal{S}_{D}(\mathrm{A})$
$\displaystyle=$ $\displaystyle\det{(\mathrm{A})},$
$\displaystyle\mathcal{S}_{n>D}(\mathrm{A})$ $\displaystyle=$ $\displaystyle
0.$

##  Appendix B MTMVMG Action

In this Appendix we fill some gaps in the derivation of the MTMVMG action in
Sec. III. Let us begin with the MTMVMG Hamiltonian (22) which can be expressed
in terms of the time derivative of spatial vielbein and the scalar field,

$\displaystyle\tensor{\dot{e}}{{}^{I}_{i}}$ $\displaystyle\approx$
$\displaystyle\frac{\delta\mathcal{H}_{\text{MTMVMG}}}{\delta\tensor{\pi}{{}_{I}^{i}}}$
(103) $\displaystyle=$
$\displaystyle\frac{1}{M_{\text{Pl}}^{D-2}}\Bigg{[}\frac{N}{\det{(e)}}\left(\tensor{\pi}{{}^{I}_{i}}-\frac{1}{D-2}\tensor{\pi}{{}^{k}_{K}}\tensor{e}{{}^{K}_{k}}\tensor{e}{{}^{I}_{i}}\right)$
$\displaystyle+\,(\nabla_{i}N_{j})\delta^{IJ}\tensor{e}{{}^{j}_{J}}$
$\displaystyle+\,|\det{(X)}|\lambda
MW(\psi)\left(\gamma_{ik}\tensor{\tilde{e}}{{}^{k}_{L}}\delta^{IK}-\frac{1}{D-2}\tensor{Y}{{}^{K}_{L}}\tensor{e}{{}^{I}_{i}}\right)$
$\displaystyle\times\,\sum_{n=1}^{D-1}\sum_{m=1}^{n}(-1)^{m}c_{n}\tensor{\left(Y^{m-1}\right)}{{}^{L}_{K}}\mathcal{S}_{n-m}(Y)\Bigg{]},$
$\displaystyle\dot{\psi}$ $\displaystyle\approx$
$\displaystyle\frac{\delta\mathcal{H}_{\text{MTMVMG}}}{\delta\pi}$ (104)
$\displaystyle=$ $\displaystyle\frac{N\pi}{\det{(e)}}+N^{i}\partial_{i}\psi$
$\displaystyle+\,|\det{(X)}|\lambda
M\frac{dW}{d\psi}\sum_{n=1}^{D-1}c_{n}\mathcal{S}_{n}(Y),$

respectively, such that we have the new conjugate momenta

$\displaystyle\frac{\tensor{\pi}{{}^{i}_{I}}}{\det{(e)}}$
$\displaystyle\equiv$ $\displaystyle
M_{\text{Pl}}^{D-2}(K^{ij}\delta_{IJ}\tensor{e}{{}^{J}_{j}}-K\tensor{e}{{}^{i}_{I}})$
(105)
$\displaystyle-\,\lambda\frac{W(\psi)}{2}\frac{M}{N}\Theta^{ij}\delta_{IJ}\tensor{e}{{}^{J}_{j}},$
$\displaystyle\frac{\pi}{\det{(e)}}$ $\displaystyle\equiv$
$\displaystyle\frac{\dot{\psi}}{N}-\frac{N^{i}}{N}\partial_{i}\psi-\lambda\frac{dW}{d\psi}\frac{M}{N}\Phi.$
(106)

The last terms in Eqs. (105) and (106) arise since the MTMVMG constraints
depend on the conjugate momenta $\tensor{\pi}{{}_{I}^{i}}$ and $\pi$. The
symmetrical property of $\Theta^{ij}$ cancels the antisymmetric terms, namely
$\alpha_{MN}\mathcal{P}^{[MN]}$ and $\beta_{MN}Y^{[MN]}$ out of the MTMVMG
Hamiltonian. Employing the Legendre transformation, we obtain the MTMVMG
action,

$\displaystyle S_{\text{MTMVMG}}$ $\displaystyle\equiv$
$\displaystyle\int_{\mathcal{M}}d^{D}x\,\Big{(}\tensor{\pi}{{}^{i}_{I}}\tensor{\dot{e}}{{}^{I}_{i}}+\pi\dot{\psi}-\mathcal{H}_{\text{MTMVMG}}\big{|}_{\alpha_{MN}=\beta_{MN}=0}\Big{)}$
(109) $\displaystyle=$ $\displaystyle
S_{\text{pre}}-\int_{\mathcal{M}}d^{D}x(\lambda\mathcal{C}_{0}+\lambda^{i}\mathcal{C}_{i})-\frac{2}{M_{\text{Pl}}^{D-2}}\int_{\mathcal{M}}d^{D}x\,N\det{(e)}\left(\lambda\frac{W(\psi)}{4}\frac{M}{N}\right)^{2}\left(\Theta^{ij}\Theta_{ij}-\frac{1}{D-2}\Theta^{2}\right)$
$\displaystyle-\,\frac{1}{2}\int_{\mathcal{M}}d^{D}x\,N\det{(e)}\left(\lambda\frac{dW}{d\psi}\frac{M}{N}\right)^{2}\Phi^{2}+S_{\text{matter}}$
$\displaystyle=$ $\displaystyle
S_{\text{pre}}-\int_{\mathcal{M}}d^{D}x\,(\lambda\mathcal{\bar{C}}_{0}+\lambda^{i}\mathcal{C}_{i})+\frac{2}{M_{\text{Pl}}^{D-2}}\int_{\mathcal{M}}d^{D}x\,N\det{(e)}\left(\lambda\frac{W(\psi)}{4}\frac{M}{N}\right)^{2}\left(\Theta^{ij}\Theta_{ij}-\frac{1}{D-2}\Theta^{2}\right)$
$\displaystyle+\,\frac{1}{2}\int_{\mathcal{M}}d^{D}x\,N\det{(e)}\left(\lambda\frac{dW}{d\psi}\frac{M}{N}\right)^{2}\Phi^{2}+S_{\text{matter}},$

where $S_{\text{pre}}$ and $S_{\text{matter}}$ are the precursor action and
the matter action, respectively, whereas

$\displaystyle\mathcal{C}_{0}$ $\displaystyle=$ $\displaystyle
W(\psi)M\det{(e)}|\det{(X)}|\sum_{n=1}^{D-1}\sum_{m=1}^{n}(-1)^{m}c_{n}\tensor{\left(Y^{m-1}\right)}{{}^{J}_{I}}\mathcal{S}_{n-m}(Y)$
(110)
$\displaystyle\times\left[\left(\gamma_{ik}\tensor{\tilde{e}}{{}^{k}_{J}}\tensor{e}{{}^{I}_{j}}-\frac{1}{D-2}\gamma_{ij}\tensor{Y}{{}^{I}_{J}}\right)\left(K^{ij}-K\gamma^{ij}-\frac{\lambda}{M_{\text{Pl}}^{D-2}}\frac{W(\psi)}{2}\frac{M}{N}\Theta^{ij}\right)-\frac{1}{M}\tensor{\tilde{e}}{{}^{k}_{J}}\frac{\partial}{\partial
t}\tensor{\tilde{e}}{{}^{I}_{k}}\right]$
$\displaystyle+\,M\det{(e)}\frac{dW}{d\psi}\left(\frac{\dot{\psi}}{N}-\frac{N^{i}}{N}\partial_{i}\psi\right)\Phi-\lambda
N\det{(e)}\left(\frac{dW}{d\psi}\frac{M}{N}\right)^{2}\Phi^{2},$
$\displaystyle\mathcal{C}_{i}$ $\displaystyle=$ $\displaystyle
M\det{(e)}\bigg{[}W(\psi)\nabla^{j}\bigg{(}|\det{(X)}|\sum_{n=1}^{D-1}\sum_{m=1}^{n}(-1)^{m}c_{n}\tensor{\left(Y^{m-1}\right)}{{}^{J}_{I}}\mathcal{S}_{n-m}(Y)M\tensor{Y}{{}^{K}_{J}}\delta_{KL}\tensor{e}{{}^{I}_{i}}\tensor{e}{{}^{L}_{j}}\bigg{)}$
(111)
$\displaystyle+\,|\det{(X)}|\partial_{i}\psi\frac{dW}{d\psi}\sum_{n=1}^{D-1}c_{n}\mathcal{S}_{n}(Y)\bigg{]},$

with

$\displaystyle\mathcal{\bar{C}}_{0}\equiv\mathcal{C}_{0}\big{|}_{\lambda=0}$
$\displaystyle=$ $\displaystyle
M|\det{(X)}|\Bigg{\\{}W(\psi)\left[\left(\gamma_{ik}\tensor{\tilde{e}}{{}^{k}_{J}}\tensor{e}{{}^{I}_{j}}-\frac{1}{D-2}\gamma_{ij}\tensor{Y}{{}^{I}_{J}}\right)(K^{ij}-K\gamma^{ij})-\frac{1}{M}\tensor{\tilde{e}}{{}^{k}_{J}}\frac{\partial}{\partial
t}\tensor{\tilde{e}}{{}^{I}_{k}}\right]$ (112)
$\displaystyle\times\sum_{n=1}^{D-1}\sum_{m=1}^{n}(-1)^{m}c_{n}\tensor{\left(Y^{m-1}\right)}{{}^{J}_{I}}\mathcal{S}_{n-m}(Y)$
$\displaystyle-\,\frac{dW}{d\psi}\left(\frac{\dot{\psi}}{N}-\frac{N^{i}}{N}\partial_{i}\psi\right)\sum_{n=1}^{D-1}c_{n}\mathcal{S}_{n}(Y)\Bigg{\\}}.$

##  Appendix C Linear Stability

Eigenvalues |  Stability
---|---
Real eigenvalues
$\mu_{n}<0$, for $n=1,2,3$ |  stable node
$\mu_{n}>0$, for $n=1,2,3$ |  unstable node
$\mu_{n}>0$ and $\mu_{m}\leq 0$, for $n\cup m=1,2,3$ |  unstable
$\mu_{n}>0$ and $\mu_{m}<0$, for $n\cup m=1,2,3$ |  saddle point
$\mu_{n}=0$ and $\mu_{m}<0$, for $n\cup m=1,2,3$ |  nonhyperbolic, linear stability fails to determine and other methods are needed
Complex eigenvalues
$\text{Re}(\mu_{n})<0$, for $n=1,2,3$ |  strongly stable spiral
$\text{Re}(\mu_{n})>0$, for $n=1,2,3$ |  strongly unstable spiral
$\text{Re}(\mu_{n})<0$ and $\text{Re}(\mu_{m})>0$, for $n\cup m=1,2,3$ |  saddle focus
$\text{Re}(\mu_{n})=0$ and $\text{Re}(\mu_{m})\neq 0$, for $n\cup m=1,2,3$ |  weakly center point
$\text{Re}(\mu_{n})=0$, for $n=1,2,3$ |  strongly center point
Real eigenvalues $\mu_{n}$ and complex eigenvalues $\mu_{m}$ for $n\cup
m=1,2,3$
$\mu_{n}<0$ and $\text{Re}(\mu_{m})<0$ |  weakly stable spiral
$\mu_{n}>0$ and $\text{Re}(\mu_{m})>0$ |  weakly unstable spiral
$\mu_{n}<0$ and $\text{Re}(\mu_{m})>0$ |  saddle focus
$\mu_{n}>0$ and $\text{Re}(\mu_{m})<0$ |  saddle focus
Table 3:  Stability properties of the critical point
$(x_{\psi,c},x_{V,c},x_{W,c})$ based on the three eigenvalues
$\mu_{1},\mu_{2},$ and $\mu_{3}$.

In this Appendix we consider the linear perturbation of dynamical equations
(57)-(59) in the critical points $(x_{\psi,c},x_{V,c},x_{W,c})$. First, we
expand autonomous variables around these points

$\displaystyle x_{\psi}$ $\displaystyle=$ $\displaystyle x_{\psi,c}+u_{\psi},$
(113) $\displaystyle x_{V}$ $\displaystyle=$ $\displaystyle x_{V,c}+u_{V},$
(114) $\displaystyle x_{W}$ $\displaystyle=$ $\displaystyle x_{W,c}+u_{W}.$
(115)

The first order equation of motions has the form

$\displaystyle\frac{2u_{\psi}^{\prime}}{D-1}$ $\displaystyle=$
$\displaystyle\Big{[}3(1-w_{m})x_{\psi,c}^{2}-(1+w_{m})x_{V,c}^{2}-(f_{2}(u)+w_{m}f_{1}(u))x_{W,c}^{2}-(1-w_{m})\Big{]}u_{\psi}$
$\displaystyle-\,\Bigg{[}2(1+w_{m})x_{\psi,c}-4\sqrt{\frac{D-2}{D-1}}\lambda_{V}\Bigg{]}x_{V,c}u_{V}-\Bigg{[}2(f_{2}(u)+w_{m}f_{1}(u))x_{\psi,c}-4\sqrt{\frac{D-2}{D-1}}\lambda_{W}f_{1}(u)\Bigg{]}x_{W,c}u_{W},$
$\displaystyle\frac{2u_{V}^{\prime}}{D-1}$ $\displaystyle=$
$\displaystyle\Bigg{[}2(1-w_{m})x_{\psi,c}-\sqrt{\frac{D-2}{D-1}}\lambda_{V}\Bigg{]}x_{V,c}u_{\psi}+\Bigg{[}(1+w_{m})+(1-w_{m})x_{\psi,c}^{2}$
$\displaystyle-\,3(1+w_{m})x_{V,c}^{2}-(f_{2}(u)+w_{m}f_{1}(u))x_{W,c}^{2}-\sqrt{\frac{D-2}{D-1}}\lambda_{V}x_{\psi,c}\Bigg{]}u_{V}-2(f_{2}(u)+w_{m}f_{1}(u))x_{V,c}x_{W,c}u_{W},$
$\displaystyle\frac{2u_{W}^{\prime}}{D-1}$ $\displaystyle=$
$\displaystyle\Bigg{[}2(1-w_{m})x_{\psi,c}-\sqrt{\frac{D-2}{D-1}}\lambda_{W}\Bigg{]}x_{W,c}u_{\psi}-2(1+w_{m})x_{V,c}x_{W,c}u_{V}$
(118)
$\displaystyle+\,\Bigg{[}(1+w_{m})+(1-w_{m})x_{\psi,c}^{2}-(1+w_{m})x_{V,c}^{2}-3(f_{2}(u)+w_{m}f_{1}(u))x_{W,c}^{2}-\sqrt{\frac{D-2}{D-1}}\lambda_{W}x_{\psi,c}\Bigg{]}u_{W},$

such that we can cast the above equations into the matrix form

$\begin{pmatrix}u_{\psi}^{\prime}\\\ u_{V}^{\prime}\\\
u_{W}^{\prime}\end{pmatrix}=\bm{J}\begin{pmatrix}u_{\psi}\\\ u_{V}\\\
u_{W}\end{pmatrix},$ (119)

where $\bm{J}$ is Jacobian matrix. Let the eigenvalues of $\bm{J}$ are
$\mu_{1}$, $\mu_{2}$, and $\mu_{3}$. We write the Jacobian matrix $\bm{J}$ and
its eigenvalues for each of the critical points below. These then can be used
to analyze their stability properties, which we summarize in Table 3.

  * •

The Jacobian matrix for CP1: $(0,0,0)$ is

$\bm{J}=\begin{pmatrix}-(1-w_{m})&0&0\\\ 0&1+w_{m}&0\\\
0&0&1+w_{m}\end{pmatrix},$ (120)

with eigenvalues

$\\{-(1-w_{m}),1+w_{m}\\}.$ (121)

  * •

The Jacobian matrix for CP2: $(1,0,0)$ is

$\bm{J}=\begin{pmatrix}2(1-w_{m})&0&0\\\
0&2-\lambda_{V}\sqrt{\frac{D-2}{D-1}}&0\\\
0&0&2-\lambda_{W}\sqrt{\frac{D-2}{D-1}}\end{pmatrix},$ (122)

with eigenvalues

$\Bigg{\\{}2(1-w_{m}),2-\lambda_{V}\sqrt{\frac{D-2}{D-1}},2-\lambda_{W}\sqrt{\frac{D-2}{D-1}}\Bigg{\\}}.$
(123)

  * •

The Jacobian matrix for CP3: $(-1,0,0)$ is

$\bm{J}=\begin{pmatrix}2(1-w_{m})&0&0\\\
0&2+\lambda_{V}\sqrt{\frac{D-2}{D-1}}&0\\\
0&0&2+\lambda_{W}\sqrt{\frac{D-2}{D-1}}\end{pmatrix},$ (124)

with eigenvalues

$\Bigg{\\{}2(1-w_{m}),2+\sqrt{\frac{D-2}{D-1}},2+\lambda_{W}\sqrt{\frac{D-2}{D-1}}\Bigg{\\}}.$
(125)

  * •

The Jacobian matrix for CP4: $(x_{\psi,c},0,0)$ is

$\bm{J}=\begin{pmatrix}0&0&0\\\
0&2-\lambda_{V}x_{\psi,c}\sqrt{\frac{D-2}{D-1}}&0\\\
0&0&2-\lambda_{W}x_{\psi,c}\sqrt{\frac{D-2}{D-1}}\end{pmatrix},$ (126)

with eigenvalues

$\Bigg{\\{}0,2-\lambda_{V}x_{\psi,c}\sqrt{\frac{D-2}{D-1}},2-\lambda_{W}x_{\psi,c}\sqrt{\frac{D-2}{D-1}}\Bigg{\\}}.$
(127)

  * •

One of the eigenvalues of the Jacobian matrix for CP5:
$\bigg{(}\frac{\mathcal{A}_{\pm}}{\lambda_{V}}\sqrt{\frac{D-1}{D-2}},\frac{1}{\lambda_{V}}\sqrt{\frac{(D-1)(2-\mathcal{A}_{\pm})\mathcal{A}_{\pm}}{2(D-2)}},0\bigg{)}$
is

$\mu_{1}=\left(1-\frac{\lambda_{W}}{\lambda_{V}}\right)\mathcal{A}_{\pm}.$
(128)

  * •

The Jacobian matrix for CP6:
$\bigg{(}0,\sqrt{\frac{\lambda_{W}}{\lambda_{W}-\lambda_{V}}},\sqrt{\frac{\lambda_{V}}{|f_{1}(u)|(\lambda_{W}-\lambda_{V})}}\bigg{)}$
is

$\bm{J}=\begin{pmatrix}-2&-4\lambda_{V}\sqrt{\frac{D-2}{D-1}\frac{\lambda_{W}}{\lambda_{W}-\lambda_{V}}}&4\lambda_{W}\sqrt{\frac{D-2}{D-1}\frac{|f_{1}(u)|\lambda_{V}}{\lambda_{W}-\lambda_{V}}}\\\
-\lambda_{V}\sqrt{\frac{D-2}{D-1}\frac{\lambda_{W}}{\lambda_{W}-\lambda_{V}}}&-\frac{2(1+w_{m})\lambda_{W}}{\lambda_{W}-\lambda_{V}}&-\frac{2(1+w_{m})\sqrt{|f_{1}(u)|\lambda_{V}\lambda_{W}}}{\lambda_{W}-\lambda_{V}}\\\
-\lambda_{W}\sqrt{\frac{D-2}{D-1}\frac{\lambda_{V}}{|f_{1}(u)|(\lambda_{W}-\lambda_{V})}}&-\frac{2(1+w_{m})}{\lambda_{W}-\lambda_{V}}\sqrt{\frac{\lambda_{V}\lambda_{W}}{|f_{1}(u)|}}&\frac{2(1+w_{m})\lambda_{V}}{\lambda_{W}-\lambda_{V}}\end{pmatrix},$
(129)

with eigenvalues

$\Bigg{\\{}-2(1+w_{m}),-1\pm\sqrt{1-\frac{4(D-1)\lambda_{V}\lambda_{W}}{D-2}}\Bigg{\\}}.$
(130)

  * •

The Jacobian matrix for CP7:
$\bigg{(}0,\sqrt{1-f_{1}(u)x_{W,c}^{2}},x_{W,c}\bigg{)}$ is

$\bm{J}=\begin{pmatrix}-2&0&0\\\
0&-2(1+w_{m})(1-f_{1}(u)x_{W,c}^{2})&-2(1+w_{m})f_{1}(u)x_{W,c}\sqrt{1-f_{1}(u)x_{W,c}^{2}}\\\
0&-2(1+w_{m})x_{W,c}\sqrt{1-f_{1}(u)x_{W,c}^{2}}&-2(1+w_{m})f_{1}(u)x_{W,x}^{2}\end{pmatrix},$
(131)

with eigenvalues

$\\{-2,-2(1+w_{m}),0\\}.$ (132)

  * •

One of the eigenvalues of the Jacobian matrix for CP8:
$\bigg{(}\frac{\mathcal{B}_{\pm}}{\lambda_{W}}\sqrt{\frac{D-1}{D-2}},0,\frac{1}{\lambda_{W}}\sqrt{\frac{(D-1)(2-\mathcal{B}_{\pm})\mathcal{B}_{\pm}}{2(D-2)f_{1}(u)}}\bigg{)}$
is

$\mu_{1}=\left(1-\frac{\lambda_{V}}{\lambda_{W}}\right)\mathcal{B}_{\pm}.$
(133)

## References

  * Fierz and Pauli [1939] M. Fierz and W. Pauli, Proc. Roy. Soc. Lond. Sect. A 173, 211 (1939).
  * van Dam and Veltman [1970] H. van Dam and M. J. G. Veltman, Nucl. Phys. B22, 397 (1970).
  * Zakharov [1970] V. I. Zakharov, J. Exp. Theor. Phys. Lett. 12, 312 (1970).
  * Vainshtein [1972] A. I. Vainshtein, Phys. Lett. B 39, 393 (1972).
  * Boulware and Deser [1972] D. G. Boulware and S. Deser, Phys. Rev. D 6, 3368 (1972).
  * Arkani-Hamed  _et al._ [2003] N. Arkani-Hamed, H. Georgi, and M. D. Schwartz, Annals Phys. (Amsterdam) 305, 96 (2003).
  * Creminelli  _et al._ [2005] P. Creminelli, A. Nicolis, M. Papucci, and E. Trincherini, J. High Energy Phys. 09, 003 (2005).
  * de Rham and Gabadadze [2010] C. de Rham and G. Gabadadze, Phys. Rev. D 82, 044020 (2010).
  * de Rham  _et al._ [2011a] C. de Rham, G. Gabadadze, and A. J. Tolley, Phys. Rev. Lett. 106, 231101 (2011a).
  * Hassan and Rosen [2012a] S. F. Hassan and R. A. Rosen, Phys. Rev. Lett. 108, 041101 (2012a).
  * Hassan  _et al._ [2012a] S. F. Hassan, R. A. Rosen, and A. Schmidt-May, J. High Energy Phys. 02, 026 (2012a).
  * de Rham  _et al._ [2012] C. de Rham, G. Gabadadze, and A. J. Tolley, Phys. Lett. B 711, 190 (2012).
  * de Rham  _et al._ [2011b] C. de Rham, G. Gabadadze, and A. J. Tolley, J. High Energy Phys. 11, 093 (2011b).
  * Hassan and Rosen [2012b] S. F. Hassan and R. A. Rosen, J. High Energy Phys. 04, 123 (2012b).
  * Mirbabayi [2012] M. Mirbabayi, Phys. Rev. D 86, 084006 (2012).
  * Golovnev [2012] A. Golovnev, Phys. Lett. B 707, 404 (2012).
  * Hassan  _et al._ [2012b] S. F. Hassan, A. Schmidt-May, and M. von Strauss, Phys. Lett. B 715, 335 (2012b).
  * Kluson [2012] J. Kluson, Phys. Rev. D 86, 044024 (2012).
  * Hinterbichler [2012] K. Hinterbichler, Rev. Mod. Phys. 84, 671 (2012).
  * de Rham [2014] C. de Rham, Living Rev. Relativity 17, 7 (2014).
  * D’Amico  _et al._ [2011] G. D’Amico, C. de Rham, S. Dubovsky, G. Gabadadze, D. Pirtskhalava, and A. J. Tolley, Phys. Rev. D 84, 124046 (2011).
  * Gümrükçüoğlu  _et al._ [2011] A. E. Gümrükçüoğlu, C. Lin, and S. Mukohyama, J. Cosmol. Astropart. Phys. 11, 030 (2011).
  * Gümrükçüoğlu  _et al._ [2012] A. E. Gümrükçüoğlu, C. Lin, and S. Mukohyama, J. Cosmol. Astropart. Phys. 03, 006 (2012).
  * De Felice  _et al._ [2012] A. De Felice, A. E. Gümrükçüoğlu, and S. Mukohyama, Phys. Rev. Lett. 109, 171101 (2012).
  * Higuchi [1987] A. Higuchi, Nucl. Phys. B282, 397 (1987).
  * Fasiello and Tolley [2012] M. Fasiello and A. J. Tolley, J. Cosmol. Astropart. Phys. 11, 035 (2012).
  * Cheung and Remmen [2016] C. Cheung and G. N. Remmen, J. High Energy Phys. 04, 002 (2016).
  * Bellazzini  _et al._ [2018] B. Bellazzini, F. Riva, J. Serra, and F. Sgarlata, Phys. Rev. Lett. 120, 161101 (2018).
  * De Felice and Mukohyama [2016a] A. De Felice and S. Mukohyama, Phys. Lett. B 752, 302 (2016a).
  * De Felice and Mukohyama [2016b] A. De Felice and S. Mukohyama, J. Cosmol. Astropart. Phys. 04, 028 (2016b).
  * De Felice and Mukohyama [2017] A. De Felice and S. Mukohyama, Phys. Rev. Lett. 118, 091104 (2017).
  * Huang  _et al._ [2012] Q.-G. Huang, Y.-S. Piao, and S.-Y. Zhou, Phys. Rev. D 86, 124014 (2012).
  * Saridakis [2013] E. N. Saridakis, Classical Quantum Gravity 30, 075003 (2013).
  * Tannukij and Wongjun [2016] L. Tannukij and P. Wongjun, Eur. Phys. J. C 76, 17 (2016).
  * Hinterbichler and Rosen [2012] K. Hinterbichler and R. A. Rosen, J. High Energy Phys. 07, 047 (2012).
  * Lin and Mukohyama [2017] C. Lin and S. Mukohyama, J. Cosmol. Astropart. Phys. 10, 033 (2017).
  * Mukohyama and Noui [2019] S. Mukohyama and K. Noui, J. Cosmol. Astropart. Phys. 07, 049 (2019).
  * Aoki  _et al._ [2019] K. Aoki, A. De Felice, C. Lin, S. Mukohyama, and M. Oliosi, J. Cosmol. Astropart. Phys. 01, 017 (2019).
  * De Felice  _et al._ [2020] A. De Felice, A. Doll, and S. Mukohyama, J. Cosmol. Astropart. Phys. 09, 034 (2020).
  * Henneaux and Teitelboim [1994] M. Henneaux and C. Teitelboim,  _Quantization of gauge systems_, 2nd ed. (Princeton University Press, Princeton, NJ, 1994).
  * De Felice  _et al._ [2017] A. De Felice, S. Mukohyama, and M. Oliosi, Phys. Rev. D 96, 024032 (2017).
  * Akbar  _et al._ [2019] F. T. Akbar, B. E. Gunara, M. Iqbal, and H. Susanto, Adv. Theor. Math. Phys. 23, 1 (2019).
  * Chatterjee and Bhui [1990] S. Chatterjee and B. Bhui, Mon. Not. R. Astron. Soc 247, 57 (1990).
  * Leon  _et al._ [2013] G. Leon, J. Saavedra, and E. N. Saridakis, Classical Quantum Gravity 30, 135001 (2013).
  * Wu  _et al._ [2013] D.-J. Wu, Y.-S. Piao, and Y.-F. Cai, Phys. Lett. B 721, 7 (2013).
  * Copeland  _et al._ [1998] E. J. Copeland, A. R. Liddle, and D. Wands, Phys. Rev. D 57, 4686 (1998).
  * Alatas  _et al._ [2017] H. Alatas, A. K. Falah, R. Alditia, and B. E. Gunara, Annals Phys. (Amsterdam) 387, 299 (2017).
  * Peebles and Ratra [2003] P. J. E. Peebles and B. Ratra, Rev. Mod. Phys. 75, 559 (2003).
  * Akbar and Gunara [2015] F. T. Akbar and B. E. Gunara, J. Math. Phys. 56, 042502 (2015).
  * Lucchin and Matarrese [1985] F. Lucchin and S. Matarrese, Phys. Rev. D 32, 1316 (1985).
  * Fujita  _et al._ [2019] T. Fujita, S. Kuroyanagi, S. Mizuno, and S. Mukohyama, Phys. Lett. B 789, 215 (2019).
  * Fujita  _et al._ [2020] T. Fujita, S. Mizuno, and S. Mukohyama, J. Cosmol. Astropart. Phys. 01, 023 (2020).