# Black Holes in Multi-Metric Gravity

Kieran Wood<EMAIL_ADDRESS>School of Physics and Astronomy,
University Of Nottingham, Nottingham NG7 2RD, UK  Nottingham Centre Of
Gravity, University Of Nottingham, Nottingham NG7 2RD, UK  Paul M. Saffin
<EMAIL_ADDRESS>School of Physics and Astronomy, University Of
Nottingham, Nottingham NG7 2RD, UK  Nottingham Centre Of Gravity, University
Of Nottingham, Nottingham NG7 2RD, UK  Anastasios Avgoustidis
<EMAIL_ADDRESS>School of Physics and Astronomy,
University Of Nottingham, Nottingham NG7 2RD, UK  Nottingham Centre Of
Gravity, University Of Nottingham, Nottingham NG7 2RD, UK

###### Abstract

We construct a wide class of black hole solutions to the general theory of
ghost free multi-metric gravity in arbitrary spacetime dimension, extending
and generalising the known results in 4-dimensional dRGT massive gravity and
bigravity. The solutions are split into three generic classes based on whether
the metrics can be simultaneously diagonalised – one of which does not exist
in dRGT massive gravity nor bigravity, and is only possible when one has more
than two interacting metric fields. We also linearise the general multi-metric
theory to determine the dynamics of the massive spin-2 modes, including
examples where this can be done analytically, and use the linear theory to
discuss the stability of the 4-dimensional multi-Schwarzchild and multi-Kerr
solutions. We explain how the instabilities that plague these solutions in
dRGT massive gravity and bigravity carry across to the general multi-metric
theory, touching upon ideas of dimensional deconstruction to make sense of the
results.

##  I Introduction

Over the past decade or so, our understanding of the physics of interacting
spin-2 fields has been revolutionised. It has long been known that general
relativity (GR) is the unique local, two-derivative, nonlinear theory that
describes a single, self-interacting, massless spin-2 field – the graviton –
in four spacetime dimensions [1, 2, 3, 4, 5]. For over a century, GR has
remained our leading descriptor of the gravitational interaction, passing most
observational tests to within a remarkable degree of precision. Despite its
numerous successes, viable alternatives to GR have long been sought, as it is
firstly a crucial task to develop theories of gravity for us to test GR
against, and as there secondly remain multiple outstanding problems at the
interface between gravity and particle physics (e.g. the nature of dark matter
and dark energy, the cosmological constant problem, renormalisability etc.).
One such alternative proposal, motivated originally by trying to explain the
origin of the observed late-time accelerated expansion of the universe [6, 7],
but later also by dark matter [8, 9, 10, 11] and the hierarchy problem [12,
13, 14], involves modifying gravity by introducing additional interacting
spin-2 fields over and above the single massless graviton of GR. These
theories then go by the helpful name of ‘multi-metric gravity’, or simply,
‘multi-gravity’. There is, however, a no-go theorem, stating that theories
that include multiple _massless_ interacting spin-2 fields are inconsistent
[15]; in these multi-metric theories, all but one of the spin-2 fields must be
_massive_. Therefore, any consistent multi-metric theory must be built off the
back of a similarly consistent theory of massive gravity, which itself has a
colourful history.

The story of massive gravity dates back to 1939, when Fierz and Pauli first
wrote down the linearised theory describing a massive, self-interacting spin-2
field [16]. For a long while, it was the pervasive view that any nonlinear
completion of the linear Fierz-Pauli theory would be pathological, owing to
the emergence of the so-called Boulware-Deser (BD) ghost – a problematic
scalar mode equipped with wrong-sign kinetic term, signalling an instability
of the vacuum – as soon as nonlinear interactions were taken into account [17,
18]. However, the original BD analysis did not take into account all possible
interaction terms, and a breakthrough finally came much later in 2010 when a
satsifactory nonlinear theory of massive gravity was indeed constructed [19,
20, 21] and subsequently proved to be free of the BD ghost [22, 23, 24, 25,
26, 27, 28]. The theory, built upon groundwork laid earlier in [29, 30], is
now known as dRGT massive gravity, after its progenitors: de Rham, Gabadadze
and Tolley (there were important contributions also by Hassan and Rosen [21,
22, 23, 24]). In four spacetime dimensions, it describes the five degrees of
freedom of a single propagating massive graviton via a framework involving
interactions between the physical spacetime metric and an auxiliary reference
metric, which one inserts by hand (typically taken to be Minkowski, though
this need not be the case). By providing a kinetic term for the reference
metric, thereby promoting it to a second dynamical field, one obtains the
theory of bigravity [31], which, due to the special structure of the dRGT
interactions, is also ghost free [32]. The generalisation to multiple
interacting metric fields followed soon after in [33], although the general
multi-metric theory is only devoid of the BD ghost up to certain conditions,
upon which we shall elaborate in section II. For further details regarding the
development and phenomenology of these theories, we refer the reader to the
excellent and comprehensive reviews [34, 35] on massive gravity, as well as
[36] on bigravity.

Armed with a consistent ghost free framework for multi-metric gravity, the
natural next step is to begin to search for physical solutions that might
describe our world. Much work has focussed on such solutions in the realm of
cosmology, with a lot of initial excitement surrounding the potential of the
theories to address the dark energy problem (see any of the reviews cited
above and references therein for details). However, work has also steadily
been ongoing to construct black hole solutions in these theories and
understand their properties; such an endeavour is of course crucial to test
the theory, since these objects do exist in our universe and provide a natural
arena to look for deviations from GR [37].

In the multi-metric theory, the simplest background solutions are those in
which each metric is proportional, via a constant conformal factor, to some
common GR background [38, 39]. Clearly, these proportional solutions then
contain multi-metric analogues of all of the known GR black hole solutions. In
4-dimensional dRGT massive gravity and bigravity, however, additional non-
proportional solutions are known to exist in which each metric individually is
patterned as a GR solution, but the two metrics are not simultaneously
diagonal [40, 41, 42, 43, 44, 45, 46, 47]. Furthermore, owing to the absence
of the GR no-hair theorems [48, 49, 50] in massive gravity, these theories
contain solutions that are known only numerically and are completely foreign
to GR, in which the black holes are endowed with a cloud of massive graviton
hair [51, 52].

The linear stability of a number of these dRGT and bigravity black holes has
also been studied, with some intriguing results. In [53, 54], it was shown
that the solution where both metrics are simultaneously proportional to the
($D=4$) Schwarzchild(-dS) metric suffers from a radial instability that is
present for certain values of the graviton mass. The instability takes the
same form as the Gregory-Laflamme instability afflicting 5-dimensional black
strings [55, 56, 57]; we shall see later on why this result should not be
surprising. The proportional Kerr solution also suffers from this radial
instability [54], as well as a superradiant one that is again dependent on the
size of the graviton mass [54, 58, 59, 60]. On the other hand, the non-
proportional Schwarzchild solution appears to be linearly stable [61, 62],
although whether or not it excites a (non-BD) ghost is unknown. Likewise, the
final state of any of the aforementioned instabilities of the proportional
solutions remains to be determined.

The above results – which are summarised nicely in the review [63] – as we
stated, are known only for dRGT massive gravity and bigravity, which contain
respectively one and two dynamical metric fields. It remains an open question
as to how the relevant black hole physics carries over to the general multi-
metric theory. In this work, we close this gap somewhat, by explicitly
constructing a wide variety of black hole solutions to the general multi-
metric theory (with both classes of allowed interaction structures – see
section II) in an arbitrary number of spacetime dimensions, generalising and
extending the known dRGT/bigravity results in $D=4$. We also discuss the
stability of the proportional solutions in the $D=4$ multi-metric theory,
determining how the Schwarzchild and Kerr instabilities manifest when there
are multiple massive spin-2 fields, rather than just one. The extensions to
the multi-metric theory are natural, in the sense that many of the bigravity
results carry over in the manner one might naively expect, although there are
some subtleties regarding the non-proportional solutions that do arise in the
multi-metric theory that are not present in dRGT/bigravity. Throughout, we
construct a number of explicit example black hole spacetimes to illustrate
these points, as well as tie in to ideas of dimensional deconstruction [64,
13] to make a clear the link between the instabilities of multi-metric black
holes and higher dimensional black strings.

The structure of the paper, then, is as follows: in section II we review the
general multi-metric theory in arbitrary spacetime dimension, and introduce
its metric and vielbein formulations that will both aid calculations later; in
section III we construct a wide array of arbitrary dimension background black
hole solutions in each of the proportional, non-proportional and (new for
multi-gravity) partially proportional branches; in section IV we linearise the
general multi-metric theory to express explicitly the dynamics of the spin-2
mass modes, providing examples where one may do this analytically; in section
V we use this linear theory to extend the dRGT/bigravity results regarding the
instability of the proportional Schwarzchild and Kerr solutions to the general
multi-metric theory; finally we conclude in section VI.

We work with natural units $c=\hbar=G=1$ throughout, and always use a mostly-
plus metric signature.

##  II Review of multi-metric gravity

The action for multi-metric gravity, living on some $D$-dimensional spacetime
manifold $\mathcal{M}_{D}$, can be written conveniently in the vielbein
formalism as a sum of $N$ Einstein-Hilbert kinetic terms [65, 66, 67] together
with an interaction potential of degree $D$ coupling the various basis 1-forms
(see e.g. [33, 12, 14, 68]):

$\displaystyle S$ $\displaystyle=S_{K}+S_{V}+S_{M}$ (1) $\displaystyle S_{K}$
$\displaystyle=\sum_{i=0}^{N-1}\frac{M_{i}^{D-2}}{2}\int_{\mathcal{M}_{D}}R^{(i)}_{ab}\wedge{\star}^{(i)}e^{(i)ab}$
(2) $\displaystyle S_{V}$ $\displaystyle=-\sum_{i_{1}\ldots
i_{D}=0}^{N-1}\int_{\mathcal{M}_{D}}\varepsilon_{a_{1}\ldots
a_{D}}T_{i_{1}\ldots i_{D}}e^{(i_{1})a_{1}}\wedge\ldots\wedge
e^{(i_{D})a_{D}}\;.$ (3)

The tetrad basis 1-forms are $e^{(i)a}=e_{\mu}^{({i}){a}}\text{d}x^{\mu}$,
where the indices run from 0 to $D-1$, with the vielbeins defined through
$g^{(i)}_{\mu\nu}=e_{\mu}^{({i}){a}}e_{\nu}^{({i}){b}}\eta_{ab}$, and the
shorthand $e^{(i)ab\ldots}$ in the kinetic term means $e^{(i)a}\wedge
e^{(i)b}\wedge\ldots$. We say that the $(i)$ labels refer to a particular
‘site’ within the interaction structure; indices are then raised/lowered site-
wise: Latin indices with $\eta_{ab}$ and Greek indices with
$g^{(i)}_{\mu\nu}$, while we can swap between Latin and Greek indices using
the vielbeins, via changes of basis. $R^{(i)}_{ab}$ is the curvature 2-form
associated with the $i$-th (Levi-Civita) connection, with one index lowered by
$\eta_{ab}$, so that the kinetic term is just $N$ copies of the usual
Einstein-Hilbert action, written in the convenient language of differential
forms. The $T_{i_{1}\ldots i_{D}}=T_{(i_{1}\ldots i_{D})}$ are symmetric
coefficients that characterise the interactions between the tetrads. Finally,
$S_{M}$ is the action for the collective matter fields coupled to the theory.

We note that since we are working in an arbitrary number of dimensions, in
principle the higher dimensional Lovelock invariants (the Euler-Poincaré
forms, in differential form language) may also be included in the kinetic term
$S_{K}$ [69, 70], although we restrict ourselves to only include the Einstein-
Hilbert term for simplicity.

One can also express the theory described by Eq. (1) in the more commonly used
metric formalism [21, 71, 72], where the interaction term is instead written
in terms of the characteristic building block matrices

$S_{i\rightarrow j}=\sqrt{g_{(i)}^{-1}g_{(j)}}\;,$ (4)

with the matrix square root defined in the sense that $(S^{2}_{i\rightarrow
j})^{\mu}_{\;\nu}=g^{(i)\mu\lambda}g^{(j)}_{\lambda\nu}$. Explicitly, the
metric version of the potential term is:

$S_{V}=-\sum_{i,j}\int\differential^{D}x\sqrt{-\det
g_{(i)}}\sum_{m=0}^{D}\beta_{m}^{(i,j)}e_{m}(S_{i\rightarrow j})\;,$ (5)

where the $\beta_{m}^{(i,j)}=\beta_{m}^{(j,i)}$ are arbitrary constants
related to the $T_{i_{1}\ldots i_{D}}$ of the vielbein formalism in a manner
we shall soon specify, and the $e_{m}(S)$ are elementary symmetric polynomials
of the eigenvalues of $S$, given by:

$\displaystyle e_{0}(\lambda_{1},\lambda_{2},...,\lambda_{D})$
$\displaystyle=1$ (6) $\displaystyle
e_{1}(\lambda_{1},\lambda_{2},...,\lambda_{D})$ $\displaystyle=\sum_{1\leq
i\leq D}\lambda_{i}$ (7) $\displaystyle
e_{2}(\lambda_{1},\lambda_{2},...,\lambda_{D})$ $\displaystyle=\sum_{1\leq
i<j\leq D}\lambda_{i}\lambda_{j}$ (8) $\displaystyle\vdots$ $\displaystyle
e_{k}(\lambda_{1},\lambda_{2},...,\lambda_{D})$ $\displaystyle=\sum_{1\leq
j_{1}<j_{2}<...<j_{k}\leq D}\lambda_{j_{1}}...\lambda_{j_{k}}$ (9)
$\displaystyle\vdots$ $\displaystyle
e_{D}(\lambda_{1},\lambda_{2},...,\lambda_{D})$
$\displaystyle=\lambda_{1}\lambda_{2}...\lambda_{D}\,.$ (10)

They can also be constructed iteratively in terms of the trace of $S$, as:

$e_{m}(S)=-\frac{1}{m}\sum_{n=1}^{m}(-1)^{n}\Tr(S^{n})e_{m-n}(S)\;.$ (11)

The structure of the building block matrices means that $S_{i\rightarrow
j}=S_{j\rightarrow i}^{-1}$, so there is a sense in which the interactions are
_oriented_ [39, 73]: we say that a term in the potential, Eq. (5), which
contains $S_{i\rightarrow j}$ (_not_ $S_{j\rightarrow i}$) is positively
oriented with respect to the $i$-th metric and negatively oriented with
respect to the $j$-th metric. The orientation will affect the equations of
motion for the $i$-th and $j$-th metrics differently, as we will soon see, and
it is accounted for in the vielbein formalism within the structure of the
$T_{i_{1}\ldots i_{D}}$.

These metric interactions, as they must, take precisely the special dRGT form
that is required to remove the Boulware-Deser ghost [19, 20, 21]. It was
argued in [74] that the vielbein theory described by the action (1) is
therefore ghost free only if it has an equivalent description in metric form,
which happens whenever the so-called Deser-van Nieuwenhuizen symmetric
vielbein condition,

$e_{\;\;\;\;\;a}^{(i)\mu}e_{\mu}^{(j)b}=e^{(i)\mu b}e^{(j)}_{\mu a}\;,$ (12)

is satisfied. The known vielbein models that satisfy this condition are those
involving only pairwise interactions with no cycles (a cycle is e.g.
$1\rightarrow 2\rightarrow 3\rightarrow 1$), so such models are _manifestly_
ghost free [75, 76, 39], though more recently constructions that evade the
arguments of [74] yet nevertheless remain ghost free were given in [77, 78].

Both the metric and vielbein formalisms of multi-metric gravity prove useful
calculational tools in different situations: indeed, we shall employ both at
different points in this work, depending on which is more appropriate. When
working with the vielbein formalism, we shall choose, however, to restrict to
those ghost free vielbein theories that do permit a metric description, for
simplicity. In particular, we shall take ‘chain’ type interactions [39], where
the $i$-th metric interacts only with its nearest neighbours, and the
interactions are always positively oriented from $i$ to $(i+1)$111In
principle, ‘star’ type interactions: where many metrics all couple to some
common central metric but not to each other, are also allowed. The most
general ghost free interaction one could choose in multi-gravity (at least,
from those vielbein theories permitting a metric description) then consists of
arbitrary combinations of ‘star’ and ‘chain’ type interactions (see [39] for
more details). We will look at how our calculations in the later sections
change for a star type interaction in appendix C.. Such a choice is natural if
one thinks of the theory as arising from some sort of dimensional
deconstruction [64, 13, 14]. In terms of the $T_{i_{1}\ldots i_{D}}$, this
means that one can only permit terms of the form $T_{iiii\ldots}$,
$T_{i+1,iii\ldots}$, $T_{i-1,iii\ldots}$, $T_{i+1,i+1,ii\ldots}$ and so on.

Amongst the general class of theories described by Eq. (1), a particular model
is specified entirely by a choice for both the number of metrics $N$ and the
$T_{i_{1}\ldots i_{D}}$ coefficients. For chain type interactions, one can
always neatly parametrise the $T_{i_{1}\ldots i_{D}}$ in terms of some new
constants $\beta_{m}^{(i,i+1)}$, which characterise the interactions between
the $i$-th and $(i+1)$-th metrics, as [14]

$\displaystyle D!\,T_{iiii\ldots i}$
$\displaystyle=\beta_{0}^{(i,i+1)}+\beta_{D}^{(i-1,i)}$ (13) $\displaystyle
D!\,T_{\\{i+1\\}^{m}\\{i\\}^{D-m}}$ $\displaystyle=\beta_{m}^{(i,i+1)}\;,$
(14)

where $i=0,\ldots,N-1$ and $m=0,\ldots,D$, with
$\beta_{m}^{(-1,0)}=\beta_{m}^{(N-1,N)}=0$ since on each of the two boundaries
of the interaction chain the corresponding metrics have only one nearest
neighbour. The factor of $D!$ is included so that these $\beta_{m}^{(i,i+1)}$
are then precisely the same $\beta_{m}^{(i,j)}$ as in the metric formalism Eq.
(5). The sense of interaction orientation is encoded in $T_{iiii\ldots i}$,
where positively oriented interactions contribute a $\beta_{0}$ term while
negatively oriented interactions contribute a $\beta_{D}$ term. In practice,
we will often choose to write $\beta_{m}^{(i,i+1)}=\alpha_{i}\beta_{m}$, where
$\alpha_{-1}=\alpha_{N-1}=0$ but the rest of the $\alpha_{i}$ are equal,
restricting to the case where the interactions between all of the metrics in
the chain are characterised by the same set of parameters $\beta_{m}$. This is
both to avoid having an abundance of free parameters in the theory and
because, again, such a choice is natural from the deconstruction perspective.
However, we shall keep the generic set of $\beta_{m}^{(i,i+1)}$ in for most of
our expressions in order to be as general as possible.

The equations of motion for the general $D$-dimensional theory are:

$M^{D-2}_{i}G^{(i)\mu}_{\;\;\;\;\;\;\nu}+W^{(i)\mu}_{\;\;\;\;\;\;\nu}=T^{(i)\mu}_{\;\;\;\;\;\;\nu}\;,$
(15)

where $T^{(i)\mu}_{\;\;\;\;\;\;\nu}$ are the energy-momentum tensors for the
various sites, and the new term $W$ characterises the effect of the
interactions over and above the standard GR interactions. Explicitly, in
vielbein form, $W^{(i)\mu}_{\;\;\;\;\;\;\nu}$ reads [14] (see appendix A for
the derivation):

$\begin{split}W^{(i)\mu}_{\;\;\;\;\;\;\nu}&=D!e_{\nu}^{({i}){a}}e_{\;\;\;\;\;{[a}}^{({i}){\mu}}e_{\;\;\;\;\;{b_{1}}}^{({i}){\lambda_{1}}}\ldots
e_{\;\;\;\;\;{b_{D-1}]}}^{({i}){\lambda_{D-1}}}\\\ &\times\sum_{j_{1}\ldots
j_{D-1}}\mathcal{P}(i)T_{ij_{1}\ldots
j_{D-1}}e_{\lambda_{1}}^{({j_{1}}){b_{1}}}\ldots
e_{\lambda_{D-1}}^{({j_{D-1}}){b_{D-1}}}\;,\end{split}$ (16)

with $\mathcal{P}(i)$ counting the number of times the index $(i)$ appears in
the interaction coefficients i.e. a term with $T_{ij_{1}\ldots j_{D-1}}$ has
$\mathcal{P}(i)=1$, a term with $T_{iij_{2}\ldots j_{D-1}}$ has
$\mathcal{P}(i)=2$, and so on. The equivalent metric formalism expression is:

$\begin{split}W^{(i)\mu}_{\;\;\;\;\;\;\nu}&=\sum_{j}\sum_{m=0}^{D}(-1)^{m}\beta_{m}^{(i,j)}Y_{(m)\nu}^{\mu}(S_{i\rightarrow
j})\\\
&+\sum_{k}\sum_{m=0}^{D}(-1)^{m}\beta_{D-m}^{(k,i)}Y_{(m)\nu}^{\mu}(S_{k\rightarrow
i}^{-1})\;,\end{split}$ (17)

where (with respect to the $i$-th metric) $j$ denote positively oriented
interactions, $k$ denote negatively oriented interactions, and we define

$Y_{(m)}(S)=\sum_{n=0}^{m}(-1)^{n}S^{m-n}e_{n}(S)\;.$ (18)

The $T^{(i)\mu}_{\;\;\;\;\;\;\nu}$ are not completely arbitrary: in order to
remain ghost free, the forms of matter coupling that one can permit are
severely restricted. In general, one must couple entirely _separate_ matter
sectors to separate vielbeins, otherwise the BD ghost is resurrected [79, 80,
81]. There is the notable exception, however, where a _single_ matter source
can be coupled to multiple vielbeins in a ghost free manner through the
special ‘effective’ vielbein considered in [80, 82, 83, 84, 77].

The non-interacting theory possesses $N$ copies of diffeomorphism invariance.
Turning on the interactions, these diffeomorphisms are broken to a single
surviving diagonal subgroup, so the theory propagates a single massless
graviton (invariant under transformations of this subgroup) and $N-1$ massive
gravitons, which are linear combinations of the original metric perturbations.

As a result of the Bianchi identities for each Einstein tensor, as well as the
surviving diagonal diffeomorphism invariance, the $W$-tensor is subject to the
_Bianchi constraint_ [14]:

$\sum_{i=0}^{N-1}\absolutevalue{e^{(i)}}\nabla^{(i)\mu}W^{(i)}_{\mu\nu}=0\;.$
(19)

Whenever matter couples to one site only, or when there is no matter coupling
at all, the divergences of each $W$-tensor (i.e. each term in the above sum)
must instead vanish _individually_ , telling us that there can be no flow of
energy-momentum across the chain of interactions.

##  III Background Solutions and Black Holes

General solutions of the multi-metric theory can deviate markedly from
solutions of GR, although it has long been known that the simplest multi-
metric solutions are those where all metrics are proportional to some common
GR background [38, 39]. Such solutions are useful, because they tell us about
those multi-metric solutions that are close to what we already know from GR.
They also admit a sensible perturbative description that allows us to analyse
the mass spectrum (indeed, we shall see this in section IV). However, these
simple solutions are also restrictive: viable cosmological solutions, for
example, necessarily do not lie in this class [85, 14] and, as we shall
demonstrate shortly, there exist numerous black hole solutions where the
metrics are not proportional (a fact that has been known in bigravity for some
time [45, 46, 63]). Nevertheless, both types of solution are important, so we
shall consider each in turn, utilising both metric and vielbein formalisms
where appropriate to aid in their construction.

###  III.1 Proportional Solutions

We look for solutions to the multi-metric equations (15) where the metrics are
conformally related to one another,

$g^{(i)}_{\mu\nu}=a_{i}^{2}\bar{g}_{\mu\nu}\;,$ (20)

where at this stage $\bar{g}_{\mu\nu}$ is some arbitrary fixed metric common
to all sites. Since all of the metrics live on the same manifold
$\mathcal{M}_{D}$, one only has the freedom to rescale the coordinates to fix
_one_ of the $a_{i}$, so their values are physical, up to an overall
normalisation.

With the ansatz (20), the $i$-th vielbein and its inverse are given by:

$\displaystyle e_{\mu}^{({i}){a}}$
$\displaystyle=a_{i}\bar{e}_{\mu}^{\;\;a}\;,$ (21) $\displaystyle
e_{\;\;\;\;\;{a}}^{({i}){\mu}}$
$\displaystyle=\frac{1}{a_{i}}\bar{e}^{\mu}_{\;\;a}\;.$ (22)

Therefore, one may express the $W$-tensor in terms of $\bar{e}_{\mu}^{\;a}$
only, and so remove the vielbeins from the sum over the $j$’s entirely:

$\begin{split}W^{(i)\mu}_{\;\;\;\;\;\;\nu}&=D!\bar{e}_{\nu}^{\;\;{a}}\bar{e}_{\;\;{[a}}^{{\mu}}\bar{e}_{\;\;{b_{1}}}^{{\lambda_{1}}}\ldots\bar{e}_{\;\;{b_{D-1}]}}^{{\lambda_{D-1}}}\bar{e}_{\lambda_{1}}^{\;\;{b_{1}}}\ldots\bar{e}_{\lambda_{D-1}}^{\;\;{b_{D-1}}}\\\
&\times\sum_{j_{1}\ldots j_{D-1}}\mathcal{P}(i)T_{ij_{1}\ldots
j_{D-1}}\frac{a_{j_{1}}\ldots a_{j_{D-1}}}{a_{i}^{D-1}}\;.\end{split}$ (23)

Since all the vielbeins now belong to the same background geometry, they can
be contracted. It is simple enough to show that the vielbeins contract to:

$D!\bar{e}_{\nu}^{\;\;{a}}\bar{e}_{\;\;{[a}}^{{\mu}}\bar{e}_{\;\;{b_{1}}}^{{\lambda_{1}}}\ldots\bar{e}_{\;\;{b_{D-1}]}}^{{\lambda_{D-1}}}\bar{e}_{\lambda_{1}}^{\;\;{b_{1}}}\ldots\bar{e}_{\lambda_{D-1}}^{\;\;{b_{D-1}}}=(D-1)!\delta^{\mu}_{\nu}\;,$
(24)

and so the $W$-tensor takes the following simple, diagonal form, _irrespective
of_ the precise identity of $\bar{g}_{\mu\nu}$:

$W^{(i)\mu}_{\;\;\;\;\;\;\nu}=(D-1)!\delta^{\mu}_{\nu}\sum_{j_{1}\ldots
j_{D-1}}\mathcal{P}(i)T_{ij_{1}\ldots j_{D-1}}\frac{a_{j_{1}}\ldots
a_{j_{D-1}}}{a_{i}^{D-1}}\;.$ (25)

Taking chain type interactions, using the symmetry of the $T_{i_{1}\ldots
i_{D}}$ coefficients, and parametrising them as in Eqs. (13) and (14), one
finds that the components of the $W$-tensor are explicitly given by:

$\begin{split}W^{(i)\mu}_{\;\;\;\;\;\;\nu}=\delta^{\mu}_{\nu}\Bigg{[}&\sum_{m=0}^{D}\beta_{m}^{(i,i+1)}\binom{D-1}{m}a_{i+1}^{m}a_{i}^{-m}\\\
&+\sum_{m=0}^{D}\beta_{D-m}^{(i-1,i)}\binom{D-1}{m}a_{i-1}^{m}a_{i}^{-m}\Bigg{]}\;.\end{split}$
(26)

The Bianchi constraint Eq. (19) forces all the $a_{i}$ to be constant [38].
Lowering an index in Eq. (26) tells us that $W^{(i)}_{\mu\nu}\propto
g^{(i)}_{\mu\nu}$; therefore, we should interpret our $W$-tensors as effective
cosmological constants on each site $(i)$, which arise due to the
interactions. The downstairs index Einstein tensor is scale invariant, that
is, $G_{\mu\nu}(g)=G_{\mu\nu}(ag)$, so we have that
$G^{(i)\mu}_{\;\;\;\;\;\;\nu}=\bar{G}^{\mu}_{\;\nu}/a_{i}^{2}$. We can absorb
these factors of $a_{i}$ into the definitions of the effective cosmological
constants and energy-momentum tensors so that the multi-metric equations
reduce simply to $N$ copies of the standard Einstein equations for
$\bar{g}_{\mu\nu}$:

$\bar{G}^{\mu}_{\;\nu}+\Lambda_{i}\delta^{\mu}_{\nu}=M_{i}^{-(D-2)}\bar{T}^{(i)\mu}_{\;\;\;\;\;\;\nu}\;,$
(27)

where the effective cosmological constants $\Lambda_{i}$ are precisely the
contributions from $W^{(i)\mu}_{\;\;\;\;\;\;\nu}$. Explicitly, the appropriate
rescaled quantities are

$\displaystyle\Lambda_{i}\delta^{\mu}_{\nu}$
$\displaystyle=\frac{a_{i}^{2}}{M_{i}^{D-2}}W^{(i)\mu}_{\;\;\;\;\;\;\nu}\;,$
(28) $\displaystyle\bar{T}^{(i)\mu}_{\;\;\;\;\;\;\nu}$
$\displaystyle=a_{i}^{2}T^{(i)\mu}_{\;\;\;\;\;\;\nu}\;,$ (29)

so that all terms in Eqs. (27) now behave as if they lived on the _common_
background $\bar{g}_{\mu\nu}$ i.e. have their indices manipulated with
$\bar{g}_{\mu\nu}$.

We can take differences of Eqs. (27) to show that the effective cosmological
constants must satisfy:

$\Lambda_{i}=\bar{\Lambda}\;\;\;\forall\;i\;,$ (30)

and similarly the energy-momentum tensors must all be proportional [38],

$\bar{T}^{(i+1)\mu}_{\;\;\;\;\;\;\;\;\;\;\nu}=\left(\frac{M_{i+1}}{M_{i}}\right)^{D-2}\bar{T}^{(i)\mu}_{\;\;\;\;\;\;\nu}\;.$
(31)

This restriction on the matter sources may not be necessarily realistic in
general, but it is certainly true in vacuum, where all the energy-momentum
tensors vanish. Therefore, these proportional solutions describe perfectly
acceptable vacua at the background level in multi-metric gravity. We shall
construct some examples shortly, though before we get there, we note that the
condition on the $\Lambda$’s is not as simple as it may first seem, since on
the ‘boundary’ sites of the interaction chain (i.e. $i=0$ and $i=N-1$) one set
of $\beta$’s vanishes, so only one sum is present in the corresponding
$W$-tensor. If we denote $\Lambda^{(+)}_{i}$ as $a_{i}^{2}/M_{i}^{D-2}$ times
the sum involving $a_{i+1}$ and $\Lambda^{(-)}_{i}$ as $a_{i}^{2}/M_{i}^{D-2}$
times the sum involving $a_{i-1}$, then what we actually have is the
following:

$\displaystyle\Lambda^{(+)}_{i}+\Lambda^{(-)}_{i}$
$\displaystyle=\bar{\Lambda}\;\;\;\;\;\text{all $i$ in bulk}$ (32)
$\displaystyle\Lambda^{(+)}_{0}$ $\displaystyle=\bar{\Lambda}\;\;\;\;\;i=0$
(33) $\displaystyle\Lambda^{(-)}_{N-1}$
$\displaystyle=\bar{\Lambda}\;\;\;\;\;i=N-1\;.$ (34)

As mentioned earlier, one may rescale the coordinates to fix the value of
exactly one of the conformal factors. In vacuum, then, with all
$\bar{T}^{(i)\mu}_{\;\;\;\;\;\;\nu}=0$, the equations of motion (given a
choice for the interaction coefficients and provided that
$\bar{G}^{\mu}_{\;\nu}=-\bar{\Lambda}\delta^{\mu}_{\nu}$) specify a system of
$N$ algebraic equations for $N$ variables: namely, $\bar{\Lambda}$, as well as
the $N-1$ remaining unfixed conformal factors (see Eq. (26)). One may in
principle solve these equations to obtain the vacuum structure of the
corresponding multi-metric theory (coordinate rescaling can always be used to
fix the overall normalisation). In general, there may be multiple solutions,
the physical ones being those where the conformal factors are real; the number
of such physical solutions is in general dependent on both $N$ and
$T_{i_{1}\ldots i_{D}}$.

We further note that the splitting of the $\Lambda$’s given in Eqs. (32)-(34)
means that, unless $\bar{\Lambda}=0$, or unless one includes a bare
cosmological constant _not_ arising from the $W$-tensor on the boundary sites
to account for the missing terms, it is impossible to find solutions where
there is a constant ratio between the conformal factors throughout the chain
of interacting metrics i.e. where $a_{i+1}/a_{i}=C\;\forall\,i$.

These proportional solutions contain multi-metric analogues of all of the
known GR black hole solutions. In particular, in $D$-dimensions, if one takes
any of the Myers-Perry metrics [86] as their $\bar{g}_{\mu\nu}$, the
proportional ansatz Eq. (20) will be a solution of the multi-metric theory,
provided that there exists a physical solution to the equations for the
conformal factors. Within this general class of solutions, there are contained
many special cases that are interesting to consider on their own. We shall
look at a couple of concrete examples in $D=4$ for demonstrative purposes.

####  III.1.1 Deconstructed black string in 4d

As a straightforward example, but one with an interesting physical
interpration, one may take $\bar{g}_{\mu\nu}$ to be the Schwarzchild metric.
One then has the situation where the multi-gravity metrics are:

$\text{d}s^{2}_{(i)}=a_{i}^{2}\left[-\left(1-\frac{r_{s}}{r}\right)\text{d}\eta^{2}+\frac{\text{d}r^{2}}{\left(1-\frac{r_{s}}{r}\right)}+r^{2}\text{d}\Omega_{(2)}^{2}\right]\;,$
(35)

with $r_{s}$ the Schwarzchild radius and $\text{d}\Omega_{(2)}^{2}$ the line
element for the $2$-sphere. This ansatz solves the vacuum equations, provided
that the conformal factors are such that $\bar{\Lambda}=0$.

The physical meaning of the various conformal factors in this scenario is
clear: suppose we make the coordinate changes $r\rightarrow\tilde{r}$ and
$\eta\rightarrow t$ defined by $r=\tilde{r}/a_{i}$ and
$\text{d}\eta=\text{d}t/a_{i}$, then the $i$-th metric becomes:

$\text{d}s^{2}_{(i)}=-\left(1-\frac{r_{s}a_{i}}{\tilde{r}}\right)\text{d}t^{2}+\frac{\text{d}\tilde{r}^{2}}{\left(1-\frac{r_{s}a_{i}}{\tilde{r}}\right)}+\tilde{r}^{2}\text{d}\Omega_{(2)}^{2}\;,$
(36)

we see that an observer minimally coupled to the $i$-th metric would see a
black hole whose Schwarzchild radius is scaled by $a_{i}$. In principle, we
can imagine having separate observers minimally coupled to each metric, who
would each report seeing a black hole with a different sized horizon according
to the vacuum structure dictated by $W^{(i)\mu}_{\;\;\;\;\;\;\nu}=0$.

This situation is reminiscent of the well-known black string solutions in
higher dimensional gravity, where Schwarzchild hypersurfaces are glued
together to form an extra dimension (see [87] for the general $p$-brane
solution). The multi-metric solution where the metrics are all conformally
Schwarzchild is precisely the dimensional deconstruction [64, 13] of these
black strings, where each metric is to be thought of as corresponding to a
discrete location in the extra compact dimension (which must be an _interval_
, rather than an orbifold [14]), and the information regarding the geometry of
the extra dimension is encoded in the structure of the conformal factors.

For example, one may choose to work with a _clockwork_ theory [12, 13, 14],
which is special amongst the general multi-metric constructions in that it
further imposes that the vacuum structure is such that the conformal factors
possess a hierarchy, leading to one end of the chain of metrics being
exponentially suppressed compared to the other i.e. something like:

$a_{i}=\frac{a_{0}}{q^{i}}\,,$ (37)

with $q\gtrsim 1$. The idea is that by coupling matter to the suppressed end
of the chain, we engineer a suppressed coupling to the surviving massless
graviton, since one can show that the structure of the zero-mode is directly
proportional to the vacuum structure [12, 14] – we will see this explicitly in
Section IV . This way, one can imagine a situation whereby the fundamental
scale of the theory is small, but matter interactions with the massless
graviton are still at the Planck scale, with a view to solving the hierarchy
problem (see [88, 89, 90, 91, 92] and references therein for an overview of
this idea in non-gravitational contexts).

One must choose their $\beta$’s in the manner prescribed by Eq. (27) to ensure
that such a vacuum solution exists; two example models were constructed in
[14], one being essentially a deconstruction of the 5D Randall-Sundrum (RS1)
braneworld model [93, 94, 95, 96].

With this vacuum structure, we see that the system is solved by a series of 4D
Schwarzchild metrics whose horizon sizes decrease by a factor $q$ as one moves
along the chain of interacting metrics. This looks a lot like the AdS black
string in 5D [97], where the system is solved by a metric with Schwarzchild
hypersurfaces multiplied by an exponential factor that decays smoothly as one
moves along the extra dimension. Indeed, for the parameter choices
corresponding to the deconstructed RS1 model considered in [14], the solution
(35) is precisely the deconstruction of the RS black string. The general 5D
continuum limit of 4D multi-gravity [14, 13] is a more complicated scalar-
tensor braneworld, but the idea remains the same, with the continuum theory
admitting Ricci-flat hypersurfaces.

####  III.1.2 Kerr-Newman-(Anti-)de Sitter black holes in 4d

The most general black hole solution one is able to write down in 4D GR, as a
consequence of the various no-hair theorems (see e.g. [48, 49, 50]), is the
Kerr-Newman-(A)dS metric, which describes a rotating, charged black hole
living in a universe with non-zero cosmological constant. It is a solution to
the Einstein-Maxwell equations, so for our multi-metric scenario we include
the following matter action:

$S_{M}=\sum_{i=0}^{N-1}\int_{\mathcal{M}_{D}}F^{(i)}\wedge{\star}^{(i)}F^{(i)}\;,$
(38)

where the $F^{(i)}=\text{d}A^{(i)}$ are separate electromagnetic field
strengths on each site, given as exterior derivatives of the corresponding
$U(1)$ connections $A^{(i)}$. In components, the field strengths are
$F^{(i)}_{\mu\nu}=\partial_{\mu}A^{(i)}_{\nu}-\partial_{\nu}A^{(i)}_{\mu}$ and
the corresponding energy-momentum tensors are:

$T^{(i)\mu}_{\;\;\;\;\;\;\nu}=F^{(i)\mu}_{\;\;\;\;\;\alpha}F_{\nu}^{(i)\alpha}-\frac{1}{4}\delta^{\mu}_{\nu}F^{(i)}_{\alpha\beta}F^{(i)\alpha\beta}\;.$
(39)

The common metric $\bar{g}_{\mu\nu}$ is typically written in Boyer-Lindquist
coordinates as [98, 99, 100]:

$\begin{split}\text{d}{\bar{s}}^{2}&=-\frac{\Delta_{r}}{\rho^{2}\Xi^{2}}\left(\text{d}{t}-j\sin^{2}{\theta}\text{d}{\phi}\right)^{2}+\frac{\rho^{2}}{\Delta_{r}}\text{d}{r}^{2}\\\
&+\frac{\Delta_{\theta}}{\rho^{2}}\text{d}{\theta}^{2}+\frac{\Delta_{\theta}\sin^{2}{\theta}}{\rho^{2}\Xi^{2}}\left(j\text{d}{t}-(r^{2}+j^{2})\text{d}{\varphi}\right)^{2}\;,\end{split}$
(40)

where $j$ is the rotation parameter222The rotation parameter is usually
written $a$, but we choose $j$ to avoid confusion with the conformal factors.
and the various functions are defined as

$\displaystyle\Delta_{r}$
$\displaystyle=\left(r^{2}+j^{2}\right)\left(1-\frac{\bar{\Lambda}}{3}r^{2}\right)-r_{s}r+r_{Q}^{2}\;,$
(41) $\displaystyle\Delta_{\theta}$
$\displaystyle=1+\frac{\bar{\Lambda}}{3}j^{2}\cos^{2}{\theta}\;,$ (42)
$\displaystyle\Xi$ $\displaystyle=1+\frac{\bar{\Lambda}}{3}j^{2}\;,$ (43)
$\displaystyle\rho^{2}$ $\displaystyle=r^{2}+j^{2}\cos^{2}{\theta}\;,$ (44)

with $r_{Q}$ a scale related to the electric charge in a manner that will be
fixed below. If $\bar{\Lambda}>0$, the metric is Kerr-Newman-dS, while if
$\bar{\Lambda}<0$ the metric is Kerr-Newman-AdS.

With this choice for $\bar{g}_{\mu\nu}$, the Einstein tensor acquires
contributions from both the cosmological constant and the charge. The former
is accounted for by the $W$-tensor in exactly the manner previously described,
while the latter is supplied by the non-trivial electromagnetic fields:

$A^{(i)}=\frac{Q_{i}r}{\rho^{2}\Xi}\left(\text{d}t+j\sin^{2}{\theta}\text{d}\varphi\right)\;,$
(45)

provided that the charges are expressed in terms of $r_{Q}$ as

$Q_{i}=\sqrt{2}M_{i}a_{i}r_{Q}\;.$ (46)

One can check that the field equations for the electromagnetic fields,
$\nabla^{(i)\mu}F_{\mu\nu}^{(i)}=0$, are also satisfied for this choice of
$A_{\mu}^{(i)}$.

This solution contains all the other proportional black hole solutions of the
$D=4$ multi-metric theory, in the various limits where one takes combinations
of the characteristic parameters ($r_{s},r_{Q},\bar{\Lambda},j$) to 0. For
example, taking $r_{Q}\rightarrow 0$ gives us the Kerr-(A)dS solution, if we
also take $j\rightarrow 0$ we get Schwarzchild-(A)dS, then taking
$\bar{\Lambda}\rightarrow 0$ we recover the Schwarzchild (black string)
solution from earlier. If instead we keep $r_{Q}$ but send $j$ and
$\bar{\Lambda}$ to 0, we get Reissner-Nordstrom. Finally, keeping
$\bar{\Lambda}$ but sending all relevant black hole parameters to 0 we recover
the static coordinate form of the dS vacua that were previously constructed in
FLRW coordinates in [14].

###  III.2 Non-Proportional Solutions

The existence of the charged black hole solutions described above is reliant
on one having $N$ non-interacting copies of the Maxwell action for entirely
separate electromagnetic fields, each minimally coupled to its own separate
metric, yet each taking the same special field configuration that cancels the
Einstein tensor contribution given in Eq. (45). As alluded to earlier, such a
situation is not necessarily realistic, so we are motivated to try and find
solutions where there is, say, a single matter sector coupled to only one
distinguished metric (such a situation, again, is natural from the
deconstruction perspective, where it is analogous to placing matter on a brane
at some distinguished location in the extra dimension).

In this scenario, the proportional ansatz breaks, so the solutions that we
find describe situations where the various metrics cannot be simultaneously
diagonalised. For bigravity ($N=2$), the corresponding _non-bidiagonal_
solutions were constructed in a series of works that established the exact set
of metrics describing static, charged, rotating and asymptotically (A)dS black
holes of this type in $D=4$ dimensions [44, 45, 46, 47]. Here, we generalise
the bigravity results to the multi-metric theory, establishing a wide class of
_non-multidiagonal_ black hole solutions in arbitrary dimension. Unlike the
proportional solutions from before, the field equations do not reduce simply
to equivalent copies of the standard GR field equations; nevertheless, the
metrics are still patterned sitewise as GR solutions, so we refer to these
solutions as ‘GR-adjacent’.

####  III.2.1 Rotating (A)dS black holes in arbitrary dimension

We first look for solutions describing rotating but non-charged
$D$-dimensional black holes that are asymptotically (A)dS. We do not include
charge for the arbitrary $D$ case, as generally this changes the form of the
metric non-trivially whenever there is also rotation, although we shall soon
provide the $D=4$ solution where it is simple enough to include both charge
and rotation, as well as the charged but non-rotating solution in arbitrary
$D$. It proves useful to work in the metric formalism here, and express the
metrics in Kerr-Schild coordinates, where the line elements read [101, 47]:

$g^{(i)}_{\mu\nu}=a_{i}^{2}\left(\bar{g}_{\mu\nu}+2\phi_{i}l_{\mu}l_{\nu}\right)\;.$
(47)

Here, $\bar{g}_{\mu\nu}$ is taken to be the metric of $D$-dimensional (A)dS
space, $\phi_{i}$ are scalar functions whose form will be given shortly, and
$l$ is a null vector that is tangent to a null-geodesic congruence on (A)dS.
Following the conventions of [101], the (A)dS metric is best expressed in
ellipsoidal coordinates, where the line element takes the following form:

$\begin{split}\differential\bar{s}^{2}=&-W\left(1-\lambda
r^{2}\right)\differential t^{2}+F\differential r^{2}\\\
&+\sum_{k=1}^{n}\frac{r^{2}+j_{k}^{2}}{\Xi_{k}}\left(\differential\mu_{k}^{2}+\mu_{k}^{2}\differential\varphi_{k}^{2}\right)\\\
&+\frac{\lambda}{W\left(1-\lambda
r^{2}\right)}\left[\sum_{k=1}^{n}\frac{(r^{2}+j_{k}^{2})\mu_{k}\differential\mu_{k}}{\Xi_{k}}\right]^{2}\;.\end{split}$
(48)

Some elaboration is required here: the coordinate system comprises a time
coordinate $t$, a radial coordinate $r$, $\lfloor(D-1)/2\rfloor$ azimuthal
coordinates $\varphi_{k}$, and $\lfloor D/2\rfloor$ coordinates $\mu_{k}$ that
satisfy $\sum_{k=1}^{\lfloor D/2\rfloor}\mu_{k}^{2}=1$. The sums then run to
$n=\lfloor D/2\rfloor$, although if $D$ is even then there is one fewer
azimuthal coordinate relative to when $D$ is odd, so one should in the even
dimensional case set $\varphi_{n}=\differential\varphi_{n}=0$. The $j_{k}$ are
at this stage simply $\lfloor(D-1)/2\rfloor$ parameters that describe the
ellipticity of the spacetime foliation ($j_{k}=0$ then just gives the (A)dS
metric in spherical coordinates), although they will become genuine rotation
parameters after extending the metrics by the null vector $l$. In the even
dimensional case, one should also therefore set $j_{n}=0$. The various
functions appearing in Eq. (48) are:

$\displaystyle\lambda$ $\displaystyle=\frac{2\Lambda}{(D-2)(D-1)}\;,$ (49)
$\displaystyle\Xi_{k}$ $\displaystyle=1+\lambda j_{k}^{2}\;,$ (50)
$\displaystyle W$ $\displaystyle=\sum_{k=1}^{n}\frac{\mu_{k}^{2}}{\Xi_{k}}\;,$
(51) $\displaystyle F$ $\displaystyle=\frac{r^{2}}{1-\lambda
r^{2}}\sum_{k=1}^{n}\frac{\mu_{k}^{2}}{r^{2}+j_{k}^{2}}\;,$ (52)

and the required null vector (and its dual 1-form) that is tangent to a null-
geodesic congruence in this spacetime is:

$\displaystyle l^{\mu}\partial_{\mu}$ $\displaystyle=-\frac{1}{1-\lambda
r^{2}}\partial_{t}+\partial_{r}-\sum_{k=1}^{n}\frac{j_{k}}{r^{2}+j_{k}^{2}}\partial_{\varphi_{k}}\;,$
(53) $\displaystyle l_{\mu}\differential x^{\mu}$
$\displaystyle=W\differential t+F\differential
r-\sum_{k=1}^{n}\frac{j_{k}\mu_{k}^{2}}{\Xi_{k}}\differential\varphi_{k}\;.$
(54)

Finally, the scalar functions $\phi_{i}$ are given by:

$\phi_{i}=\frac{r_{s,i}}{2U}\;,$ (55)

with each metric now allowed its own independent Schwarzchild radius
$r_{s,i}$, and where the function $U$ differs between the even and odd
dimensional cases [47, 101]:

$\displaystyle U$
$\displaystyle=\sum_{k=1}^{n}\frac{\mu_{k}^{2}}{r^{2}+j_{k}^{2}}\prod_{s=1}^{n}\left(r^{2}+j_{s}^{2}\right)\;\;\;\;\;\;\;\text{(D
is odd)}$ (56) $\displaystyle U$
$\displaystyle=r\sum_{k=1}^{n}\frac{\mu_{k}^{2}}{r^{2}+j_{k}^{2}}\prod_{s=1}^{n-1}\left(r^{2}+j_{s}^{2}\right)\;\;\;\;\;\text{(D
is even)}$ (57)

In principle, one could similarly allow for $(i)$ labels on the rotation
parameters and cosmological constants, though the field equations force these
to be equal on all sites [47], else
$M_{i}^{D-2}G^{(i)\mu}_{\;\;\;\;\;\;\nu}+W^{(i)\mu}_{\;\;\;\;\;\;\nu}=0$ leads
to an inconsistency (unless one has all $j_{k}=0$, as we shall see).

With the ansatz Eq. (47) for the metrics, the Einstein tensors are simply:

$G^{(i)\mu}_{\;\;\;\;\;\;\nu}=-\frac{\Lambda}{a_{i}^{2}}\delta^{\mu}_{\nu}\;,$
(58)

as we had for the proportional solutions.

The utility of writing the metrics in Kerr-Schild form lies in how it
simplifies the calculation of the $W$-tensors: the fact that $l$ is null with
respect to both $\bar{g}_{\mu\nu}$ and $g^{(i)}_{\mu\nu}$ (i.e.
$l^{\mu}l_{\mu}=0$) means that its contribution to the interaction building
block matrices $S_{i\rightarrow j}$ is nilpotent, leading to an early
truncation in the expansion of the matrix square root. In particular, one may
show that, with the metrics given by Eq. (47), the $m$-th power of
$S_{i\rightarrow j}$ is simply [47]:

$(S_{i\rightarrow
j}^{m})^{\mu}_{\;\nu}=a_{j}^{m}a_{i}^{-m}\left[\delta^{\mu}_{\nu}-m\left(\phi_{i}-\phi_{j}\right)l^{\mu}l_{\nu}\right]\;.$
(59)

The next step is to substitute this into Eq. (18) to obtain the form of the
$Y_{(m)}$ matrices that enter the $W$-tensors. The null character of $l$ helps
us here also, since it means that the trace of $S_{i\rightarrow j}$ picks up
only the contribution from $\delta^{\mu}_{\nu}$, so the elementary symmetric
polynomials are:

$e_{n}(S_{i\rightarrow j})=a_{j}^{n}a_{i}^{-n}\binom{D}{n}\;.$ (60)

Using this in Eq. (18), along with the binomial coefficient identities

$\displaystyle\sum_{n=0}^{m}(-1)^{n}\binom{D}{n}$
$\displaystyle=(-1)^{m}\binom{D-1}{m}\;,$ (61)
$\displaystyle\sum_{n=0}^{m}n(-1)^{n}\binom{D}{n}$
$\displaystyle=D(-1)^{m}\binom{D-2}{m-1}\;,$ (62)

one finds that:

$\begin{split}Y_{(m)\nu}^{\mu}(S_{i\rightarrow
j})&=(-1)^{m}a_{j}^{m}a_{i}^{-m}\\\
&\times\left[\binom{D-1}{m}\delta^{\mu}_{\nu}+\binom{D-2}{m-1}\left(\phi_{i}-\phi_{j}\right)l^{\mu}l_{\nu}\right]\;.\end{split}$
(63)

One can then substitute into Eq. (17) to determine the components of the
$W$-tensors. Since for chain type interactions the only building block
matrices present in the potential are the nearest neighbour ones – namely,
$S_{i\rightarrow i+1}$ and $S_{i-1\rightarrow i}$, Eq. (17) involves only two
sums:

$\displaystyle W^{(i)\mu}_{\;\;\;\;\;\;\nu}$
$\displaystyle=\sum_{m=0}^{D}(-1)^{m}\beta_{m}^{(i,i+1)}Y_{(m)\nu}^{\mu}(S_{i\rightarrow
i+1})$
$\displaystyle+\sum_{m=0}^{D}(-1)^{m}\beta_{D-m}^{(i-1,i)}Y_{(m)\nu}^{\mu}(S_{i\rightarrow
i-1})\;,$ (64)

recalling from section II that $S_{i-1\rightarrow i}^{-1}=S_{i\rightarrow
i-1}$. With the $Y_{(m)}$ given by Eq. (63), and the $\phi_{i}$ given by Eq.
(55), the explicit expression for the components becomes:

$\begin{split}W^{(i)\mu}_{\;\;\;\;\;\;\nu}&=\delta^{\mu}_{\nu}\left[\sum_{m=0}^{D}\beta_{m}^{(i,i+1)}\binom{D-1}{m}a_{i+1}^{m}a_{i}^{-m}+\sum_{m=0}^{D}\beta_{D-m}^{(i-1,i)}\binom{D-1}{m}a_{i-1}^{m}a_{i}^{-m}\right]\\\
&+\frac{l^{\mu}l_{\nu}}{2U}\left[\frac{a_{i+1}}{a_{i}}\left(r_{s,i}-r_{s,i+1}\right)\Sigma_{i}^{(+)}+\frac{a_{i-1}}{a_{i}}\left(r_{s,i}-r_{s,i-1}\right)\Sigma_{i}^{(-)}\right]\;,\end{split}$
(65)

where we have defined:

$\displaystyle\Sigma^{(+)}_{i}$
$\displaystyle=\sum_{m=0}^{D}\beta_{m}^{(i,i+1)}\binom{D-2}{m-1}a_{i+1}^{m-1}a_{i}^{1-m}\;,$
(66) $\displaystyle\Sigma^{(-)}_{i}$
$\displaystyle=\sum_{m=0}^{D}\beta_{D-m}^{(i-1,i)}\binom{D-2}{m-1}a_{i-1}^{m-1}a_{i}^{1-m}\;,$
(67)

The part of the $W$-tensor proportional to $\delta^{\mu}_{\nu}$ is exactly as
in Eq. (26) for the proportional solutions, but now there are additional off-
diagonal terms proportional to $l^{\mu}l_{\nu}$. Clearly, for the ansatz Eq.
(47) to be a solution to the multi-metric vacuum equations, since the Einstein
tensor is diagonal, we need these off-diagonal $W$-tensor components to
vanish. There are a few ways of achieving this, and we shall go through them
in order of increasing complexity.

The simplest is to take all of the Schwarzchild radii to be the same,

$r_{s,i}=r_{s}\;\;\;\forall\,i\;,$ (68)

which just recovers the proportional solution from before.

The second is to make all of the $\Sigma_{i}$’s vanish:

$\displaystyle\Sigma_{i}^{(+)}$ $\displaystyle=0\;,$ (69)
$\displaystyle\Sigma_{i}^{(-)}$ $\displaystyle=0\;\;\;\forall\,i\;.$ (70)

These conditions are polynomial equations that fix the ratios of neighbouring
conformal factors. Furthermore, due to the recurrence relation for the
binomial coefficients:

$\binom{n}{k}=\binom{n-1}{k}+\binom{n-1}{k-1}\;,$ (71)

they cause part of the sum defining the diagonal components of the $W$-tensors
in Eq. (65) to vanish also. Thus, the conditions (69) and (70) strengthen Eqs.
(28), which give the value of the cosmological constant, so as to include only
the non-vanishing parts of the diagonal $W$-tensor components:

$\begin{split}&\Bigg{[}\sum_{m=0}^{D}\beta_{m}^{(i,i+1)}\binom{D-2}{m}a_{i+1}^{m}a_{i}^{-m}\\\
&+\sum_{m=0}^{D}\beta_{D-m}^{(i-1,i)}\binom{D-2}{m}a_{i-1}^{m}a_{i}^{-m}\Bigg{]}=\frac{\Lambda
M_{i}^{D-2}}{a_{i}^{2}}\;.\end{split}$ (72)

For a solution to exist, the same ratios of conformal factors that solve Eqs.
(69) and (70) must _also_ satisfy Eq. (72), which can only happen when the
$\beta_{m}^{(i,i+1)}$ are finely tuned to allow for this possibility. That is
to say, Eq. (72) acts as a constraint on which multi-metric theories permit
this class of solutions. If the parameters are chosen such that this
constraint is satisfied, then one has succeeded in constructing their non-
multidiagonal black hole. For example, if one wishes to find a solution with
$\Lambda=0$, one only needs to fix $\beta_{0}$ and $\beta_{D}$ in terms of the
other $\beta$’s and conformal factor ratios to satisfy the constraint (72).

In bigravity, these are the only two available options. We see this by
realising that the sum involving $\beta_{D-m}$ in the $(i+1)$-th off-diagonal
$W$-tensor is proportional to the sum involving $\beta_{m}$ in the $i$-th off-
diagonal $W$-tensor. Explicitly, one has that:

$\Sigma^{(-)}_{i+1}=\left(\frac{a_{i}}{a_{i+1}}\right)^{D-2}\Sigma^{(+)}_{i}\;,$
(73)

so when $N=2$ and the only terms present are $\Sigma^{(+)}_{0}$ and
$\Sigma^{(-)}_{1}$ (since $\beta_{m}^{(-1,0)}=\beta_{m}^{(N-1,N)}=0$), one sum
vanishing implies the other does too. Indeed, our equations recover exactly
the bigravity results in e.g. [46, 47] when one takes $N=2,D=4$ and
$\Lambda=0$333In $D=4$, as we will shortly see, the coordinates of Eq. (48)
can be parametrised by $\mu_{1}=\sin\theta$ and $\mu_{2}=\cos\theta$, leading
to a pair of metrics that are precisely those in [47], and are a coordinate
transformation away from those written in Eddington-Finkelstein coordinates in
[46]..

In the multi-metric scenario, however, whenever one has $N>2$, there are more
$\Sigma^{(+)}_{i}$ to play with, and these do not necessarily all have to be 0
_even if_ some of the others are. More precisely, if one has vanishing
$\Sigma^{(+)}_{I-1}$ for some specific $i=I$, this implies only the vanishing
of $\Sigma^{(-)}_{I}$ and nothing else, so the requirement that the off-
diagonal parts of $W^{(I)\mu}_{\;\;\;\;\;\;\;\nu}$ vanish may still be
satisfied by _either_ $r_{s,I}=r_{s,I+1}$ or $\Sigma^{(+)}_{I}=0$ (which is
independent of $\Sigma^{(+)}_{I-1}$). Therefore, in the multi-metric theory
with $N>2$, the most general way to solve the vacuum equations is to allow for
combinations of both $\Sigma^{(+)}_{k}=0$ and $r_{s,k}=r_{s,k+1}$, for
different $k\subset i$. This corresponds to the situation where some – but not
all – of the metrics can be simultaneously diagonalised. Those sites where
$\Sigma^{(+)}_{k}=0$ then each fix $a_{k+1}/a_{k}$ only, and feed in to modify
only the $k$-th and $(k+1)$-th equations for the cosmological constant in the
manner described by Eq. (72), while the rest of those equations are unchanged
i.e. still involve $\binom{D-1}{m}$ instead of $\binom{D-2}{m}$. If there are
$n$ total sites that have $\Sigma_{k}^{(+)}=0$ ($n$ can therefore be at most
$N-1$), $n+1$ of the conformal factors are a priori fixed before checking
whether these cosmological constant equations are satsified (the additional 1
is fixed by rescaling the coordinates). Then, the $N$ diagonal equations split
into $N-n$ algebraic equations for $\Lambda$ and the remaining $N-n-1$ free
conformal factors, as well as $n$ equations that become constraints on the
$\beta_{m}^{(i,i+1)}$ parameters of the theory. Again, this means that only
for finely tuned parameters can these solutions exist. Indeed, the solutions
form a set of measure zero in the $\beta_{m}^{(i,i+1)}$ parameter space.

This third and most general branch of solutions (which we dub the “partially
proportional” branch) is fiddly and awkward to deal with; in practice, we
shall stick to considering the first two branches of solutions i.e. where
either all the Schwarzchild radii are the same (the proportional solutions) or
all the $\Sigma_{i}$’s vanish (the non-proportional solutions). We note,
however, that a subtlety arises on the overlap of these first two branches,
when one has _both_ $r_{s,i}=r_{s,i+1}$ and $\Sigma^{(+)}_{i}=0$ for every
$i$. In this scenario, one cannot accept the solution as valid: as we will
demonstrate in section IV, the masses of all the linearised perturbations
vanish, and as mentioned in section I, theories involving multiple interacting
_massless_ gravitons are known to be pathological [15].

####  III.2.2 4d Kerr-Newman-(Anti-)de Sitter revisited

In $D=4$, as we mentioned, it is not difficult to extend the above analysis to
include charge, providing a non-proportional generalisation of the Kerr-
Newman-(A)dS metric from section III.1.2. In the $D=4$ case of Eq. (48), there
is only one azimuthal coordinate $\varphi_{1}=\varphi$, only one rotation
parameter $j_{1}=j$ (hence one $\Xi_{1}=\Xi$), and only two $\mu_{k}$
coordinates, which we can parametrise without loss of generality as [101]:

$\mu_{1}=\sin\theta\;,\;\;\;\mu_{2}=\cos\theta\;.$ (74)

With these definitions, the (A)dS metric explicitly takes the form:

$\begin{split}\differential\bar{s}^{2}=&-\frac{\left(1-\lambda
r^{2}\right)\Delta_{\theta}\differential
t^{2}}{\Xi}+\frac{\rho^{2}\differential r^{2}}{\left(1-\lambda
r^{2}\right)\left(r^{2}+j^{2}\right)}\\\
&+\frac{\rho^{2}\differential\theta^{2}}{\Delta_{\theta}}+\frac{\left(r^{2}+j^{2}\right)\sin^{2}\theta\differential\varphi^{2}}{\Xi}\;,\end{split}$
(75)

where $\Delta_{\theta}$ and $\rho^{2}$ are as in Eqs. (42) and (44). The null
vector and 1-form are given by:

$\displaystyle l^{\mu}\partial_{\mu}$ $\displaystyle=-\frac{1}{1-\lambda
r^{2}}\partial_{t}+\partial_{r}-\frac{j}{r^{2}+j^{2}}\partial_{\varphi}\;,$
(76) $\displaystyle l_{\mu}\differential x^{\mu}$
$\displaystyle=\frac{\Delta_{\theta}\differential
t}{\Xi}+\frac{\rho^{2}\differential r}{\left(1-\lambda
r^{2}\right)\left(r^{2}+j^{2}\right)}-\frac{j\sin^{2}\theta\differential\varphi}{\Xi}\;,$
(77)

while the function $U$ becomes:

$U=\frac{\rho^{2}}{r}\;.$ (78)

As in section III.1.2, to incorporate charge we must include copies of the
Maxwell action for the matter sector,

$S_{M}=\sum_{k\subset
i}\int_{\mathcal{M}_{D}}F^{(k)}\wedge{\star}^{(k)}F^{(k)}\;,$ (79)

only now we allow the freedom to couple to only some subset of sites $k\subset
i$, rather than to all of them. This of course breaks the proportional ansatz,
but we are looking for the non-/partially proportional solutions anyway.

In $D=4$, the only effect of the charge is to modify the scalar functions
$\phi_{i}$ to:

$\phi_{i}=\frac{1}{2U}\left(r_{s,i}-\frac{r_{Q,i}}{r}\right)\;,$ (80)

where each metric is now allowed an independent $r_{Q,i}$ as well as
$r_{s,i}$. The extension of $\phi_{i}$ to include charge while keeping
rotation is particularly simple in $D=4$; this is assuredly not the case in
higher dimensions, which is why we treat the charged and rotating cases
separately outside of $D=4$. The coordinate transformation that retrieves the
Boyer-Lindquist form of the metric, Eq. (40), from the Kerr-Schild form is
given in [101].

To account for the additional contributions to the Einstein tensors on the
charged sites, the corresponding electromagnetic fields must take the form:

$A^{(k)}=\frac{Q_{k}r}{\rho^{2}}l_{\mu}\differential x^{\mu}\;,$ (81)

where as before

$Q_{k}=\sqrt{2}M_{k}a_{k}r_{Q,k}\;.$ (82)

Of course, only the sites $k\subset i$ that have an electromagnetic field
coupling have non-vanishing $r_{Q,i}$, so the multi-metric system comprises
Kerr-Newman-(A)dS metrics on those sites with charge, and Kerr-(A)dS metrics
on those without.

Since the only alteration wrought by the inclusion of charge in the system is
the new form of $\phi_{i}$ given by Eq. (80), the only changes to the
$W$-tensor components are in the part multiplying $l^{\mu}l_{\nu}$, where
$r_{s,i}$ is shifted to $r_{s,i}-r_{Q,i}/r$. Therefore, in line with the
discussion of the previous subsection, there are again in principle (depending
on the parameters of the theory) three classes of solutions, corresponding to
the three ways in which one can make these off-diagonal components vanish. The
proportional solutions have the same $r_{s}$ _and_ $r_{Q}$ on every site, the
non-proportional solutions have $\Sigma^{(+)}_{i}=0$ on every site, and the
partially proportional solutions have a combination of the two for different
$k\subset i$. These three classes of Kerr-Newman-(A)dS solutions are then the
most general GR-adjacent black holes in the $D=4$ multi-metric theory.

We note that the proportional solutions only exist in this case when one has a
separate matter sector coupled to each metric in the chain, which we already
argued is unrealistic. If one wishes to solve the system with, say, a single
copy of the Maxwell action coupled to one metric only, then the solution must
lie in either the non-proportional or partially proportional branch, as
$r_{Q,i}-r_{Q,i+1}$ cannot vanish when only one $r_{Q}$ is present. Therefore,
charged and rotating black holes of this type _only_ exist in the $D=4$ multi-
metric theory if the interaction coefficients are finely tuned to allow for
this possibility.

####  III.2.3 Charged (A)dS black holes in arbitrary dimension

Finally, we look for the $D$-dimensional solution for an asymptotically (A)dS,
charged, but _non_ -rotating black hole. As alluded to earlier, when there is
no rotation, one may in principle allow for different cosmological constants
on each site, and so replace $\lambda\rightarrow\lambda_{i}$ in Eq. (48).
However, this then spoils the utility of the Kerr-Schild ansatz, since the
null vectors pick up their own $(i)$ indices, which means that the expressions
for the building block matrices $S_{i\rightarrow j}$ no longer simplify in the
manner they did before. To proceed, we must change tack; once again, the
vielbein formalism comes in handy.

Following [45], it proves useful to now express the metrics in (advanced)
Eddington-Finkelstein coordinates, where the line elements read:

$\displaystyle\differential s^{2}_{(i)}=a_{i}^{2}\Bigg{[}$
$\displaystyle-\left(1-\lambda_{i}r^{2}-\frac{r_{s,i}}{r^{D-3}}+\frac{r_{Q,i}^{2}}{r^{2(D-3)}}\right)\differential
v^{2}$ $\displaystyle+2\differential v\differential
r+r^{2}\left(\differential\theta^{2}+\sin^{2}\theta\differential\varphi^{2}\right)$
$\displaystyle+r^{2}\cos^{2}\theta\differential\Omega^{2}_{(D-4)}\Bigg{]}\;,$
(83)

with $\differential\Omega_{(D-4)}^{2}$ the unit round metric on the
$(D-4)$-sphere, given by:

$\differential\Omega_{(D-4)}^{2}=\differential\psi_{1}^{2}+\sum_{k=2}^{D-4}\left(\prod_{m=1}^{k-1}\sin^{2}\psi_{m}\right)\differential\psi_{k}^{2}\;,$
(84)

As before, to account for the extra contribution to the Einstein tensors on
the sites with non-vanishing $r_{Q,i}$, the corresponding electromagnetic
fields must take the profile:

$A^{(k)}=\sqrt{\frac{D-2}{2(D-3)}}\frac{Q_{k}}{r^{D-3}}\differential t\;,$
(85)

again with:

$Q_{k}=\sqrt{2}M_{k}a_{k}r_{Q,k}\;.$ (86)

With the ansatz Eq. (III.2.3) for the metric, the tetrads have the following
form:

$\displaystyle e^{(i)0}$ $\displaystyle=\differential
r-\frac{1}{2}\left(2-\lambda_{i}r^{2}-\frac{r_{s,i}}{r^{D-3}}+\frac{r_{Q,i}^{2}}{r^{2(D-3)}}\right)\differential
v$ (87) $\displaystyle e^{(i)1}$ $\displaystyle=\differential
r+\frac{1}{2}\left(\lambda_{i}r^{2}+\frac{r_{s,i}}{r^{D-3}}-\frac{r_{Q,i}^{2}}{r^{2(D-3)}}\right)\differential
v$ (88) $\displaystyle e^{(i)2}$ $\displaystyle=r\differential\theta,$ (89)
$\displaystyle e^{(i)3}$ $\displaystyle=r\sin\theta\differential\varphi,$ (90)
$\displaystyle e^{(i)4}$ $\displaystyle=r\cos\theta\differential\psi_{1},$
(91) $\displaystyle e^{(i)5}$
$\displaystyle=r\cos\theta\sin\psi_{1}\differential\psi_{2},$ (92)
$\displaystyle\vdots$

which one can use to determine the components of the $W$-tensors by
substituting into Eq. (16).

The diagonal part of the $i$-th $W$-tensor is precisely as in Eq. (26) and
again gives ($M_{i}^{D-2}/a_{i}^{2}$ times) the $i$-th effective cosmological
constant $\Lambda_{i}$, in the usual manner. The only non-vanishing off-
diagonal term is the $rv$ component, which becomes:

$\begin{split}W^{(i)r}_{\;\;\;\;\;\;v}=\frac{1}{2}&\Bigg{\\{}\frac{a_{i+1}}{a_{i}}\left[\frac{r_{s,i}-r_{s,i+1}}{r^{D-3}}+(\lambda_{i}-\lambda_{i+1})r^{2}-\frac{r_{Q,i}-r_{Q,i+1}}{r^{2(D-3)}}\right]\Sigma_{i}^{(+)}\\\
&+\frac{a_{i-1}}{a_{i}}\left[\frac{r_{s,i}-r_{s,i-1}}{r^{D-3}}+(\lambda_{i}-\lambda_{i-1})r^{2}-\frac{r_{Q,i}-r_{Q,i-1}}{r^{2(D-3)}}\right]\Sigma_{i}^{(-)}\Bigg{\\}}\;.\end{split}$
(93)

Remarkably, both the charge and cosmological constant contributions factor
nicely into the collection of terms multiplying $\Sigma_{i}^{(\pm)}$. The
solutions then split into the usual three branches: the proportional solutions
are those where one has $r_{s,i}=r_{s}\,,\;\Lambda_{i}=\Lambda$ and
$r_{Q,i}=r_{Q}\;\forall\,i$, the non-proportional solutions are those with
$\Sigma_{i}^{(+)}=0\;\forall\,i$, and the partially proportional solutions are
those with combinations of either on different sites.

For the non-proportional case, all of the effective cosmological constants are
fixed by the diagonal part of the vacuum equations once $\Sigma_{i}^{(+)}=0$
has fixed the ratios of all the adjacent conformal factors. However, because
the $\Lambda_{i}$ no longer have to be the same, the $\beta_{m}^{(i,i+1)}$ are
now _unconstrained_.

For the partially proportional case, some, but not all, of the adjacent sites
_do_ still need to have the same $\Lambda$, which, as previously, makes things
more fiddly. The situation is as follows: if there are $n$ total sites
$k\subset i$ which have $\Sigma_{k}^{(+)}=0$, then $n+1$ conformal factors are
a priori fixed before considering the diagonal part of the vacuum equations.
Furthermore, there are only $n+1$ independent $\Lambda_{i}$, as the sites that
do not have $\Sigma_{k}^{(+)}=0$ must have $\Lambda_{k}=\Lambda_{k+1}$. Then,
from the $N$ total diagonal equations, $n$ of them fix $n$ of the independent
$\Lambda_{i}$, leaving only one free, so that one is left with $N-n$ algebraic
equations for the remaining $N-n-1$ free conformal factors and remaining one
free $\Lambda$. Again, the $\beta_{m}^{(i,i+1)}$ are unconstrained here, as
the $n$ equations that constrained them in the rotating case now act instead
to fix the initially free $\Lambda_{i}$.

Because of the lack of constraints on which parameters one must choose for
these non-rotating solutions to exist, one can imagine being able to choose
their model in such a way as to engineer essentially whatever effective
cosmological constants one would like. In particular, taking the limit where
$r_{s,i}$ and $r_{Q,i}$ go to 0, one is faced with the tantalising prospect of
finding a multi-metric theory with a non-proportional dS vacuum whose
effective cosmological constant on the ‘physical’ metric i.e. the one matter
couples to, is small but non-vanishing, with a view to addressing the
cosmological constant problem. However, one expects that in order to do this,
the present fine-tuning problem regarding the size of the cosmological
constant would simply be transferred to a fine-tuning problem regarding the
potential coefficients and/or vacuum structure, so nothing is actually
alleviated. Furthermore, as we have seen, one requires the effective
cosmological constants on each metric to be the same anyway, if one wishes the
theory to admit rotating black holes, which likely comprise the majority of
real, physical black holes that exist in nature. The argument is also a purely
classical one; there is nothing to be said about quelling the contribution of
the QFT vacuum energy. Therefore, it is unlikely that the cosmological
constant problem may be addressed in this manner.

###  III.3 Other Solutions and Remarks on Stability

As we mentioned previously, the wide class of analytic black hole solutions
constructed above were all built starting from metrics that are known to GR.
Even the non- and partially proportional solutions that do differ from GR – in
the sense that the field equations are not simply $N$ copies of the standard
GR equations, as is the case for the proportional solutions – are patterned
sitewise as GR solutions. However, in dRGT massive gravity and bigravity there
also exist additional non-GR-adjacent black hole solutions, which one would
expect to carry over to the multi-metric theory as well.

Firstly, however, there is one remaining GR-adjacent solution we missed that
is worth a mention. Although we did not explicitly write down this solution,
the BTZ black hole in $D=3$ [102] does admit multi-metric analogues in each of
the proportional, non-proportional and partially proportional branches. This
is most easily seen by expressing the BTZ metric in Kerr-(A)dS like
coordinates (see e.g. [103]) and following through exactly as before, the only
changes resulting from the fact that there is no longer a $\theta$ coordinate.
If one wishes to include charge in $D=3$, the electromagnetic field must also
take a logarithmic profile, required to supply the Einstein tensor charge
contribution [104].

As for the non-GR-adjacent solutions, owing to the complexity of the field
equations, solutions must generically be found numerically444Actually, there
is a further class of analytic black hole solutions in $D$-dimensional dRGT
massive gravity, for a particular theory whose non-dynamical reference metric
takes a special degenerate form, and where only the first 4 interaction terms
are included [105, 106, 107] (see also [108, 109, 110]). However, when both
metrics (or more) are dynamical, the degeneracy is of course problematic, so
we shall leave these solutions alone (they are also not completely general
i.e. do not include the full set of ghost free interactions when $D>4$).. An
important such class of numerically determined solutions, at least in
bigravity, describing a family of asymptotically AdS black holes endowed with
a cloud of massive graviton hair was found in [51], in a comprehensive
analysis that also studied (for example) the case where the two metrics are
diagonal but not proportional. This result was extended to asymptotically
_flat_ hairy black holes in [52]. The solutions are found by considering a
generic static, spherically symmetric ansatz for the metrics and utilising the
Bianchi constraint, Eq. (19), to reduce the field equations to a set of
coupled first order ODEs that are then numerically integrated to determine the
appropriate metric functions. A detailed derivation of these equations and
their final form is publicly available in a Mathematica notebook online [111].
Again, one expects that these solutions should extend to the multi-metric
theory, just with the precise form of the ODEs altered by the extra
interactions, though we leave this calculation to future work.

The solutions we constructed earlier, together with the proposed hairy
solutions, comprise the full set of currently known black hole solutions of
multi-metric gravity. It is natural to ask of these solutions the question of
their stability: which (if any) of them are stable to perturbations, and if
so, in what regimes of parameter space are they stable?

Again, significant progress in this direction has already been made in both
$D=4$ dRGT massive gravity and bigravity. There, it is known that the
bidiagonal Schwarzchild(-dS) solution i.e. where both metrics are
proportionally Schwarzchild(-dS), suffers from a spherically symmetric radial
instability for certain values of the graviton mass [53, 54, 61]. The
instability takes precisely the same form as the Gregory-Laflamme (GL)
instability that notoriously plagues black string solutions in higher
dimensions [55, 56, 57]. This should be unsurprising – as we saw in section
III.1.1, the proportional Schwarzchild solutions are just dimensionally
deconstructed black strings. The bidiagonal Kerr solution also possesses this
radial instability [54], as well as a superradiant instability for the
azimuthal modes, again dependent on the graviton mass [54, 58, 59, 60].

The non-bidiagonal Schwarzchild solution, on the other hand, seems to be
linearly stable to metric perturbations [61, 62], although the perturbations
take an unconventional, non-Fierz-Pauli form, and as of yet no analysis exists
as to whether these solutions may contain (non-BD) ghosts. It is interesting
to note that the particular combination of the interaction coefficients,
$\Sigma_{0}^{(+)}=0$, that gives rise to the non-bidiagonal black hole
solutions also shows up when one considers the theory’s cosmological
implications: namely, it is one of two possible ways to satisfy the Bianchi
constraint when one tries to find FLRW solutions in dRGT/bigravity (the multi-
metric extension is in [14], which further admits a partially proportional
branch). The non-proportional branch of cosmological solutions (in
dRGT/bigravity) is littered with pathologies e.g. not all massive graviton
degrees of freedom propagate at linear level, ghosts appear at non-linear
level etc. [36, 112, 113, 114, 115, 116]. However, these pathologies in
cosmology appear to be intimately related to the symmetries of the FLRW
background, which is potentially why they do not show up for the non-
proportional black hole solutions. Also, the scale factors in the cosmological
solutions are time dependent, whereas for our black hole solutions they are
constant, so it is unclear whether one may even identify the
$\Sigma_{i}^{(+)}=0$ branches in either scenario anyway.

One would expect that the dRGT/bigravity results for the non-proportional
black hole solutions extend naturally to the full multi-metric theory, with
the perturbations taking the same generic form as in [61, 62]555Precisely, we
expect that (e.g.) Eq. (11) in [61] has instead of the term written
$\mathcal{A}(r_{S}-r_{f})$ a term that involves both $(r_{s,i}-r_{s,i+1})$ and
$(r_{s,i}-r_{s,i-1})$, accounting for the interactions going in both
directions.. Therefore, the non-proportional multi-metric solutions should
still be linearly stable. Determining the stability of the partially
proportional branch may be more complicated, as on some sites the
perturbations will acquire a standard Fierz-Pauli mass term while on others
they will not. Due to the inherent complexity of these calculations, and the
fact that, as we found, the interaction coefficients generally must be finely
tuned to permit these branches of solutions anyway, an explicit determination
of the stability of the non- and partially proportional multi-metric black
hole solutions is left to future work.

This leaves us with the task of extending the results regarding the stability
of the proportional solutions in bigravity to the general multi-metric
scenario, which we will eventually get to in section V. To begin to address
this point, however, we must first determine the dynamics of the graviton mass
modes, which requires that we consider linear perturbations around the
proportional backgrounds.

##  IV Linearised perturbations and mass modes

The notion of mass arises naturally in Minkowski space as a Casimir invariant
of the Poincaré group relating to spacetime translations. In more general
spacetimes with fewer symmetries, it is not always obvious how this notion
should generalise. However, around the proportional solutions only, the spin-2
fluctuations all acquire a standard Fierz-Pauli mass term, and so give rise to
a well-defined mass spectrum.

The full details of the linearisation procedure around a generic proportional
solution are laid out in appendix B, following the formalism developed in [38,
117] for bigravity (it proves simpler to work in the metric formalism here).
The metrics are expanded as:

$\displaystyle g^{(i)}_{\mu\nu}$
$\displaystyle=a_{i}^{2}\bar{g}_{\mu\nu}+a_{i}M_{i}^{-\frac{D-2}{2}}h^{(i)}_{\mu\nu}\;,$
(94)

where the normalisation is to ensure the kinetic terms are canonical. The
resulting equations of motion for the perturbations read:

$\begin{split}\bar{\mathcal{E}}^{\alpha\beta}_{\;\mu\nu}h^{(i)}_{\alpha\beta}+\bar{\Lambda}h_{\mu\nu}^{(i)}+\frac{\mathcal{M}^{2}_{ij}}{2}&\left(h_{\mu\nu}^{(j)}-\bar{g}_{\mu\nu}h^{(j)}\right)\\\
&=a_{i}M^{-\frac{D-2}{2}}\bar{T}^{(i)}_{\mu\nu}\;,\end{split}$ (95)

where the Lichnerowicz operator is given by [53]:

$\begin{split}\bar{\mathcal{E}}^{\alpha\beta}_{\;\mu\nu}h_{\alpha\beta}=\frac{1}{2}\bigl{[}&-\bar{\Box}h_{\mu\nu}+\bar{\nabla}_{\mu}{\bar{\nabla}_{\alpha}{h^{\alpha}_{\;\nu}}}+\bar{\nabla}_{\nu}{\bar{\nabla}_{\alpha}{h^{\alpha}_{\;\mu}}}\\\
&-\bar{\nabla}_{\mu}{\bar{\nabla}_{\nu}{h}}+\bar{g}_{\mu\nu}\bar{\Box}h-\bar{g}_{\mu\nu}\bar{\nabla}_{\alpha}{\bar{\nabla}_{\beta}{h^{\alpha\beta}}}\\\
&-2\bar{R}^{\alpha\;\beta}_{\;\mu\;\nu}h_{\alpha\beta}\bigr{]}\;,\end{split}$
(96)

and all quantities have been rescaled by the appropriate powers of $a_{i}$, as
in section III.1, so that all tensors in Eq. (95) behave as if they live in
the common background of $\bar{g}_{\mu\nu}$. For chain type interactions, the
mass matrix, $\mathcal{M}^{2}$, for the rescaled perturbations in this common
background, is tridiagonal, with the following non-vanishing components:

$\displaystyle\mathcal{M}_{ii}^{2}$
$\displaystyle=\frac{a_{i}^{2}}{M_{i}^{D-2}}\left(\Sigma_{i}^{(+)}\frac{a_{i+1}}{a_{i}}+\Sigma_{i}^{(-)}\frac{a_{i-1}}{a_{i}}\right)\;,$
(97) $\displaystyle\mathcal{M}_{i+1,i}^{2}$
$\displaystyle=\left(\frac{a_{i+1}}{a_{i}}\right)^{4-D}\mathcal{M}_{i,i,i+1}^{2}=-\frac{a_{i}^{2}\Sigma_{i}^{(+)}}{\left(M_{i+1}M_{i}\right)^{\frac{D-2}{2}}}\;.$
(98)

As per the discussion in Section II, the mass eigenstates are linear
combinations of the $h^{(i)}_{\mu\nu}$. To obtain them, we rotate to the field
basis $\mathbf{H}_{\mu\nu}=(H^{(0)}_{\mu\nu},\ldots,H^{(N-1)}_{\mu\nu})$ in
which the mass matrix is diagonalised, related to the original basis
$\mathbf{h}_{\mu\nu}=(h^{(0)}_{\mu\nu},\ldots,h^{(N-1)}_{\mu\nu})$ via some
$N\times N$ orthogonal matrix $O$ whose columns are the mass eigenvectors;
that is,

$\displaystyle\mathbf{h}_{\mu\nu}$ $\displaystyle=O\mathbf{H}_{\mu\nu}\;,$
(99) $\displaystyle O^{T}\mathcal{M}^{2}O$
$\displaystyle=\text{diag}(m^{2}_{0},\ldots,m^{2}_{N-1})\;.$ (100)

The structure of the potential means that the mass matrix will always possess
one zero eigenvalue, whose associated eigenvector is proportional to the
vacuum structure i.e. $O^{i0}\propto a_{i}$ [39, 12, 14]. Generally, it is not
possible (save for an exceptional case that we will consider shortly) to
determine the higher mass eigenvalues/eigenvectors analytically [118, 119],
although the expressions (97) and (98) for the components allow one to readily
determine them numerically on a case-by-case basis.

Regardless, by substituting Eq. (99) into the linearised Einstein equations
(95) and multiplying by $O^{T}$ from the left, one finds the evolution
equations for the mass modes:

$\begin{split}\bar{\mathcal{E}}^{\alpha\beta}_{\;\mu\nu}H^{(i)}_{\alpha\beta}+\bar{\Lambda}H_{\mu\nu}^{(i)}+\frac{m^{2}_{i}}{2}&\left(H^{(i)}_{\mu\nu}-\bar{g}_{\mu\nu}H^{(i)}\right)\\\
&=O^{ji}a_{j}M_{j}^{-\frac{D-2}{2}}\bar{T}^{(j)}_{\mu\nu}\;.\end{split}$ (101)

Whenever the effective cosmological constant is non-vanishing, it is crucial
that the masses satisfy $m^{2}\geq 2\bar{\Lambda}/(D-1)$. This is the well-
known Higuchi bound, below which the helicity-0 graviton modes become ghost-
like [120, 121, 122, 123]. It is a sufficient condition that the lightest
massive mode exceeds this bound, since all the heavier modes then will do too.

At the point where the Higuchi bound is saturated, the helicity-0 component of
the corresponding mass mode becomes pure gauge and so the number of
propagating degrees of freedom is reduced by 1: this is the linear ‘partially
massless’ (PM) theory [124, 125, 126, 127, 128, 129]. The PM theory has been
the subject of much interest in the context of bigravity, since many of the
linear instabilities that exist when all graviton degrees of freedom propagate
disappear if the one massive mode exhibits PM invariance (see e.g. [130]).
However, there is strong evidence that the PM gauge invariance does not
survive to the full non-linear level within the dRGT framework [131, 132, 133,
134], and in the general multi-metric theory even at the linear level only the
lightest mass mode can exhibit PM invariance anyway, so any instabilities
present still exist for the heavier modes (although it should be stated that a
nonlinear theory of interacting PM spin-2 fields does exist within the realm
of conformal gravity [135, 136, 137]).

As a sanity check, one may compare the results of this section against known
results in bigravity, where the mass matrix is simple enough to diagonalise
explicitly. In bigravity, one has $N=2$ metrics, and usually denotes
$a_{0}=1,a_{1}=c$. Only $\Sigma_{0}^{(+)}$ and $\Sigma_{1}^{(-)}$ are present
in the mass matrix, and by Eq. (73) these are related as
$\Sigma_{1}^{(-)}=c^{2-D}\Sigma_{0}^{(+)}$. Therefore, one has explicitly
that:

$\mathcal{M}^{2}=\begin{bmatrix}M_{0}^{-(D-2)}c\Sigma_{0}^{(+)}&-(M_{0}M_{1})^{-\frac{D-2}{2}}\Sigma_{0}^{(+)}\\\
-(M_{0}M_{1})^{-\frac{D-2}{2}}c^{4-D}\Sigma_{0}^{(+)}&M_{1}^{-(D-2)}c^{3-D}\Sigma_{0}^{(+)}\end{bmatrix}$
(102)

which one can check has eigenvalues:

$\displaystyle m_{0}^{2}$ $\displaystyle=0\;,$ (103) $\displaystyle m_{1}^{2}$
$\displaystyle=\frac{c\Sigma_{0}^{(+)}}{M_{0}^{D-2}}\left[\frac{1+(\gamma
c)^{D-2}}{(\gamma c)^{D-2}}\right]\;,$ (104)

with the corresponding mass modes being:

$\displaystyle H^{(0)}_{\mu\nu}$
$\displaystyle=\frac{1}{\sqrt{1+c^{2}\gamma^{D-2}}}\left(h^{(0)}_{\mu\nu}+c\gamma^{\frac{D-2}{2}}h^{(1)}_{\mu\nu}\right)\;,$
(105) $\displaystyle H^{(1)}_{\mu\nu}$
$\displaystyle=\frac{1}{\sqrt{1+c^{2(D-3)}\gamma^{D-2}}}\left(h^{(1)}_{\mu\nu}-c^{D-3}\gamma^{\frac{D-2}{2}}h^{(0)}_{\mu\nu}\right)\;,$
(106)

defining $\gamma=M_{1}/M_{0}$ as the ratio of the grvaitational couplings. The
equations of motion for the mass modes, as well as the corresponding Fierz-
Pauli mass $m_{1}$, agree precisely with the bigravity results [38, 117], as
they should (though in [117] the perturbations were defined differently to in
our Eq. (94), as $a_{i}^{2}\bar{g}_{\mu\nu}+h^{(i)}_{\mu\nu}$, without the
extra factors of $a_{i}$ and $M_{i}$, so the precise structure of their
$H^{(i)}_{\mu\nu}$ in terms of the $h^{(i)}_{\mu\nu}$ is slightly different to
here – see appendix B for more details).

As before, there are many special cases that are interesting to consider in
their own right, within this general framework. A couple of illustrative
examples in $D=4$ are again useful to demonstrate, before we continue to
investigate the stability of the black hole solutions. First, we return to the
clockwork scenario.

####  IV.0.1 Clockwork gravity in 4d

We defined a clockwork theory back in section III.1.1 as a member of a
particular class of multi-metric theories that admit a vacuum structure with
an exponential hierarchy i.e. $a_{i}=a_{0}/q^{i}$ for $q\gtrsim 1$, the idea
being that we would like to generate a suppressed coupling of matter to the
massless mode without introducing such hierarchies in the parameters of the
theory. Here, we shall see how this works explicitly.

The simplest $D=4$ model one may construct that admits the desired vacuum
structure is the single scale model of [14], which takes all $M_{i}=M$ and has
potential coefficients given by $\beta_{m}^{(i,i+1)}=\alpha_{i}\beta_{m}$,
where $\alpha_{-1}=\alpha_{N-1}=0$, and666Actually, there was a sign error in
the mass matrix in our previous work [14], and we have also been more careful
to match the normalisation of the corrected mass matrix, given by Eqs. (97)
and (98), to that of the kinetic term. This means that the single scale model
we refer to now actually had gravitons with negative mass-squared (the
deconstructed RS model also in [14] was fine), although the work there was
only focussed on the _background_ cosmology so this did not come into play in
the corresponding calculations. Regardless, we shall fix the error here and
flip the sign of the potential coefficients so that the masses are positive.

  * •

$\alpha_{i}=M^{4}\;\;\;\forall\,i\neq-1,N-1$

  * •

$\beta_{0}=-6q^{-1}$

  * •

$\beta_{1}=3$

  * •

$\beta_{2}=-q$

  * •

$\beta_{3}=\beta_{4}=0\;.$

With these choices, the mass matrix has components:

$\displaystyle\mathcal{M}_{ii}^{2}$
$\displaystyle=\frac{a_{0}^{2}q^{-2i}}{M^{2}}\left(\alpha_{i}q^{-1}+\alpha_{i-1}q^{3}\right)\;,$
(107) $\displaystyle\mathcal{M}_{i+1,i}^{2}=\mathcal{M}_{i,i+1}^{2}$
$\displaystyle=-\frac{a_{0}^{2}q^{-2i}}{M^{2}}\alpha_{i}\;,$ (108)

in agreement with [14] (up to the minus sign and numerical factor we have
corrected here). Numerically, the lightest mass mode is found to roughly have
mass $m_{1}\sim qM$, with the heavier modes distributed exponentially above
this. The massless eigenvector is given by $O^{j0}=\mathcal{N}/q^{j}$, where
the normalisation is

$\mathcal{N}=\frac{1}{\sqrt{\sum_{i=0}^{N-1}q^{-2i}}}=\sqrt{\frac{1-q^{-2}}{1-q^{-2N}}}\;.$
(109)

If matter couples only to the metric $g^{(N-1)}_{\mu\nu}$ at the end of the
chain of interactions i.e. only $\bar{T}^{(N-1)}_{\mu\nu}$ is non-vanishing
(to engineer the greatest possible suppression of scales), and one chooses to
fix the overall normalisation of the conformal factors such that $a_{N-1}=1$,
then it follows that the massless mode has the following dynamics:

$\bar{\mathcal{E}}^{\alpha\beta}_{\;\mu\nu}H^{(0)}_{\alpha\beta}=\frac{1}{M_{\text{Pl}}}\bar{T}^{(N-1)}_{\mu\nu}\;,$
(110)

where the effective Planck scale is:

$M_{\text{Pl}}=\sqrt{\frac{1-q^{-2N}}{1-q^{-2}}}q^{N-1}M\;.$ (111)

This can be much larger than $M$ if the number of fields in the chain is big
enough – indeed, generating this scale hierarchy is precisely the purpose of
the clockwork mechanism.

####  IV.0.2 Deconstructed flat extra dimension

Another example, which is interesting to consider because it is simple enough
to analyse analytically, is the multi-metric model arising from the
deconstruction of a flat extra dimension. As we saw in section III.1.1, the 4D
multi-metric (proportional) solutions correspond to dimensionally
deconstructed solutions of some 5D gravitational theory (generically a scalar-
tensor braneworld), with the geometry of the extra dimension encoded in the
vacuum structure of the conformal factors. We saw this explicitly for the
black string solution, but the result is true more generally i.e. for any
$\bar{g}_{\mu\nu}$ satisfying Einstein’s equations. For example, the clockwork
vacuum structure Eq. (37) corresponds to an extra dimension that is warped (as
there is exponential damping through the chain of hypersurfaces), though one
could also consider the simpler situation where all $a_{i}=1$, corresponding
to an extra dimension that is flat.

In particular, one may choose their model parameters in such a way that the 4D
theory, in its continuum limit, becomes pure 5D GR without cosmological
constant. Explicitly, the model in question has again equivalent 4D
gravitational couplings $M_{i}=M_{(4)}$, as well as potential coefficients
given again by $\beta_{m}^{(i,i+1)}=\alpha_{i}\beta_{m}$ with
$\alpha_{-1}=\alpha_{N-1}=0$, except this time [13, 14]:

  * •

$\alpha_{i}=M^{3}_{(5)}/\delta y\;\;\;\forall\,i\neq 0,N-1$

  * •

$\beta_{0}=-6$

  * •

$\beta_{1}=3$

  * •

$\beta_{2}=-1$

  * •

$\beta_{3}=\beta_{4}=0\;,$

where $\delta y$ is the spacing of the hypersurfaces upon which the multi-
gravity metrics live in the extra dimension, and $M_{(5)}$ is the 5D
gravitational coupling, related to the 4D coupling as
$M^{2}_{(4)}=M^{3}_{(5)}\delta y$. With these parameters, one recovers 5D GR
from the 4D multi-metric theory upon taking the limit where $\delta
y\rightarrow 0$ and $N\rightarrow\infty$ while keeping the product $N\delta
y=L$ fixed (corresponding to the size of the extra dimension). One may check
that substituting these parameters into the $D=4$ multi-metric vacuum
equations (27) gives a solution with $\bar{\Lambda}=0$ and
$a_{i}=1\;\forall\,i$.

The form of the mass matrix here reflects the simplicity of the vacuum,
reading:

$\mathcal{M}^{2}=\frac{1}{\delta y^{2}}\begin{bmatrix}1&-1&0&0&\dots&0&0\\\
-1&2&-1&0&\dots&0&0\\\ 0&-1&2&-1&\dots&0&0\\\ 0&0&-1&2&\dots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&0&0&0&\dots&2&-1\\\
0&0&0&0&\dots&-1&1\end{bmatrix}_{(N\times N)}\;,$ (112)

which one can diagonalise using the techniques developed in [118]. The result
is that the mass eigenvalues are given by:

$m_{n}^{2}=\frac{4}{\delta y^{2}}\sin^{2}\left(\frac{n\pi}{2N}\right)\;,$
(113)

the massless mode is simply:

$H^{(0)}_{\mu\nu}=\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}h^{(n)}_{\mu\nu}\;,$ (114)

while the massive modes are:

$H^{(m)}_{\mu\nu}=\mathcal{N}_{m}\sum_{n=0}^{N-1}\left[\sin\left(\frac{(n+1)m\pi}{2N}\right)-\sin\left(\frac{nm\pi}{2N}\right)\right]h^{(n)}_{\mu\nu}\;,$
(115)

with normalisations $\mathcal{N}_{m}$ whose explicit expressions are
irrelevant but can easily be determined.

The mass eigenvalues and structure of the massless mode take the expected
functional forms [64, 73], with the masses of the lightest modes scaling as
$m_{n}\sim n/L$ and the heaviest as $m_{N-1}\sim N/L$. Similarly, if one again
couples matter to one of the metrics, then one recovers the expected relations
between the effective Planck scale of the massless mode and the various
gravitational couplings i.e.
$M_{\text{Pl}}^{2}=NM_{(4)}^{2}=NM_{(5)}^{3}\delta y=M^{3}_{(5)}L$.

##  V Linear stability of Proportional Solutions

Now we have all we need to discuss the linear stability of the proportional
black hole solutions. First, we note that the linearised version of the
Bianchi constraint in vacuum implies that
$\bar{\nabla}^{\nu}H^{(i)}_{\mu\nu}=\bar{\nabla}_{\mu}H^{(i)}$, and so
$\bar{\nabla}^{\nu}\bar{\nabla}^{\mu}H^{(i)}_{\mu\nu}=\bar{\Box}H^{(i)}$.
Using the latter in the trace of Eqs. (95), one is forced into de Donder gauge
for _all_ the mass modes simultaneously:

$\bar{\nabla}^{\nu}H^{(i)}_{\mu\nu}=H^{(i)}=0\;\;\;\forall i\;.$ (116)

With this restriction, the linearised equations for the mass modes become:

$\bar{\Box}H^{(i)}_{\mu\nu}+2\bar{R}^{\alpha\;\beta}_{\;\mu\;\nu}H^{(i)}_{\alpha\beta}=m_{i}^{2}H^{(i)}_{\mu\nu}\;,$
(117)

whose behaviour can then be analysed depending on one’s choice for the
background metric $\bar{g}_{\mu\nu}$.

###  V.1 Multi-Schwarzchild

Eqs. (117), with the common background $\bar{g}_{\mu\nu}$ taken as the $D=4$
Schwarzchild metric, are precisely those equations studied in the context of
the GL instability [55, 56, 57, 53, 61, 54]. In the original work, linear
perturbations to the 5D black string were considered, and split into scalar,
vector and tensor contributions à la Kaluza-Klein. Fourier decomposing around
the extra dimension, it was shown that the 4D tensor perturbations (i.e. the
Fourier coefficients of the tensor contribution) satisfy precisely Eqs. (117),
just with the corresponding Fourier momentum in place of the mass. Spherically
symmetric $s$-wave solutions, regular at the future event horizon, were found
of the form:

$h^{(4)}_{\mu\nu}=e^{\Omega t}\begin{bmatrix}h_{tt}(r)&h_{tr}(r)&0&0\\\
h_{rt}(r)&h_{rr}(r)&0&0\\\ 0&0&K(r)&0\\\
0&0&0&K(r)\sin^{2}{\theta}\end{bmatrix}\;,$ (118)

and were shown to be unstable (i.e. have $\Omega>0$) within the range:

$0<m<\mathcal{O}\left(\frac{1}{r_{s}}\right)\;.$ (119)

In [53], it was argued in bigravity that since the one massive mode satisfies
the same equations as the 4D tensor perturbations in the context of the black
string, the bi-Schwarzchild solutions in bigravity too are unstable if the
mass satisfies the inequality (119). This idea was made more concrete in [54],
who studied the dynamics of the unstable mode in detail for bi-Schwarzchild-
dS, and determined that the instability turns on when $mr_{s}\lesssim 0.86$
(which, interestingly, is precisely the threshold at which the hairy solutions
discussed briefly in section III.3 come into existence [52], signalling that
the hairy black holes may actually be the end point of the instability – for a
nice review of this situation see [63]).

Since the mass modes in the multi-metric theory all independently obey Eqs.
(117) for their respective mass eigenvalues, it is clear how the arguments of
[53, 54] generalise naturally: the multi-metric theory possesses a single
massless mode and a tower of massive modes, whose masses are given in terms of
the parameters of the theory through the mass matrix $\mathcal{M}^{2}$; _all_
of the mass modes are subject to the same stability condition, where the
proportional Schwarzchild solution becomes unstable if, for any of the
$m_{i}$, $m_{i}r_{s}\lesssim 0.86$. Consequently, the solution is stable if
and only if the lightest massive mode $H^{(1)}_{\mu\nu}$ sits above this
inequality. This is for the simple reason that if the lightest mass mode
evades the instability then necessarily so too do all of the others.
Conversely, unstable solutions may have multiple unstable mass modes.

In principle, one could construct a multi-metric theory with whatever masses
one wishes, depending on one’s choice for the potential coefficients. In
practice, however, such choices are tightly constrained, particularly by Solar
System tests of gravity, as the additional degrees of freedom associated with
the mass modes can induce marked deviations from GR even in the weak-field
regime [138]. Typically, the masses have been taken of order the Hubble
parameter today ($m\sim H_{0}\sim 10^{-33}\text{eV}$) with the goal of
addressing the question of dark energy; the Vainshtein mechanism then ensures
that Solar System tests are satisfied, as GR is restored due to nonlinear
self-interactions of the helicity-0 graviton modes [139, 140, 141, 142, 143].
However, this choice is not a necessary one stemming from any sort of
theoretical or observational requirement; in particular, a multi-metric theory
whose masses are all very heavy will also restore GR at the linear level,
without the need for any Vainshtein screening, since the Compton wavelengths
$m^{-1}$ of the mass modes are small enough that their effects would be
invisible to current experimental precision. More precisely, modifications to
weak-field GR solutions (at least in bigravity) are suppressed by a factor
$e^{-mr}$ [52], implying that for large $m$ the large scale behaviour in the
Solar System should be indistinguishable from GR777This point, together with
the fact that the heavy spin-2 still gravitates in the same way as ordinary
matter, was used (in bigravity) in [9, 8] to argue that the heavy mode could
also constitute an interesting dark matter candidate.. This is possible due to
the surviving massless mode in the multi-metric theory – in dRGT massive
gravity, where the only propagating graviton is massive, large masses are
ruled out, as predictions always differ from GR due to the vDVZ discontinuity
[144, 145].

The black hole instability affects the two different mass regimes in different
ways. When the graviton masses are ultra-light, by virtue of Eq. (119) all
astrophysically relevant black holes are unstable: if $m\sim H_{0}$ then any
Schwarzchild black hole weighing less than $10^{22}M_{\odot}$ suffers from the
instability [63]. However, for such light graviton masses, $\Omega$ scales
linearly with $m$ [146, 147, 148], so the characteristic timescale of the
instability $\Omega^{-1}$ is of order the Hubble time: while the instability
may always be present, it is not physically relevant over any observable
timescale. On the other hand, if the graviton masses are very heavy, then
while the instability is far more efficient, it affects only the very lightest
black holes. For example, if the lightest graviton lies at the TeV scale
(which is optimal for a dark matter candidate [9, 8] and is also sensible in
clockwork scenarios [12, 14]), only black holes weighing less than roughly
$10^{-22}M_{\odot}$ are unstable. However, the instability timescale is now
much shorter, and potentially _is_ physically relevant – some primordial black
holes, for example, may exist in the unstable mass range, and as initially
stable black holes evaporate by Hawking radiation, they will inevitably cross
over into the unstable regime at some time. It is difficult to say much more
than this without knowing for definite the end state of the instability (which
requires a full nonlinear analysis, though one possibility is the hairy
solution of [52]). However, it may be possible that such effects possess
observational signatures, or even signal at pathologies within the multi-
metric theory itself.

###  V.2 Multi-Kerr

If one now allows the black holes to rotate, and so instead takes the common
background $\bar{g}_{\mu\nu}$ as the Kerr metric, then the GL monopolar
instability is still present [54], but the value of $m_{i}r_{s}$ for which
modes become stable increases relative to the Schwarzchild case
($m_{i}r_{s}\sim 0.86$) with increasing black hole spin [60] (though it still
always remains order 1). There are also additional instabilities associated
with the azimuthal modes that are not present in the Schwarzchild case, as a
consequence of the superradiant instability of rotating black holes against
massive bosonic excitations [149, 150, 151]. This effect occurs when the
frequency of the perturbation satisfies:

$0<\omega<m_{A}\Omega_{\text{BH}}\;,$ (120)

where $m_{A}$ is the mode’s azimuthal quantum number and $\Omega_{\text{BH}}$
is the angular velocity of the black hole horizon, and is characterised by the
bosons forming a condensate around the black hole, which then spins down and
deposits its rotational energy into the condensate until the above bound
saturates at $\omega=m_{A}\Omega_{\text{BH}}$. The condensate then dissipates
via (almost monchromatic) quadrupolar gravitational wave emission [152].

The superradiant instability described above turns out to be most effective
when the Compton wavelength of the perturbation in question is comparable with
the horizon size of the black hole [58, 59]. Therefore, like the GL monopolar
instability ($m_{A}=0$), it is relevant for black holes that have
$m_{i}r_{s}\sim\mathcal{O}(1)$. Unlike the GL instability, however, this is
not a hard bound at which the superradiant instability switches on, rather a
statement on when its rate is fastest – all black holes with rotation
velocities satisfying Eq. (120) suffer, but the instability rate will only be
non-negligible for a certain range of black hole masses (given a value for
$m_{i}$). Consequently, in the multi-metric theory, a wider range of black
hole masses will be affected than would be with a single massive graviton,
with successively heavier mass modes having potentially relevant superradiant
instabilities for successively lighter black holes ($m_{i}r_{s}$ being order 1
corresponds to $m_{i}\sim 10^{-11}(M_{\odot}/M_{\text{BH}})$ eV, in physical
units).

The rates of the superradiant instabilities for massive spin-2 fields have
been studied semi-analytically for small black hole spins in [54],
analytically in the regime where $m_{i}r_{s}\ll 1$ in [58], and fully
numerically for $m_{i}r_{s}$ up to 0.8 and spins up to $j=0.99$ in [59], where
it was found that the dipole ($m_{A}=1$) mode is the fastest growing – for
certain regions of parameter space it is so fast that it can even affect black
hole ringdown. The authors as a result argued that a measurement of non-zero
rotation in supermassive black holes would rule out large swathes of ultra-
light graviton masses. However, it was demonstrated in [60] that in most of
the parameter space (save for the very fastest spins and largest $m_{i}r_{s}$)
the growth of the superradiant instabilities is always subdominant to that of
the GL monopolar mode, so any such constraints must take into account the
backreaction of the GL mode on the solution. Again, to do this, a nonlinear
analysis with knowledge of the end state of this instability for Kerr black
holes will be necessary. This of course requires one to have a well-posed
dynamical formulation of the multi-metric theory that is suitable for such
simulations, upon which the development is ongoing (see [153] for the case of
dRGT massive gravity with flat reference metric). That said, as an initial toy
model for the nonlinear evolution of the instability, the authors of [60]
considered the linearised system as arising from Einstein-Weyl theory [154,
155], which has a well-posed dynamical formulation [156, 157, 158], but also
contains a ghost. They nevertheless found signatures of the linear theory in
their nonlinear analysis, and postulated that this may be a general feature of
all such nonlinear analyses i.e. even in the ghost free nonlinear theories,
including the multi-metric one.

##  VI Conclusion

In this work, we have sought to extend and generalise numerous results
regarding 4-dimensional black holes in the theories of dRGT massive gravity
and bigravity to the general ghost free multi-metric theory in arbitrary
dimension. To that end, we have explicitly constructed various example black
hole spacetimes that solve the multi-metric equations of motion, including
analogues in the proportional branch of all the higher dimensional Myers-Perry
black holes of GR, as well as multiple additional solutions in which not all
metrics are simultaneously diagonalisable (including a class of solutions –
the partially proportional branch – which is not present in dRGT/bigravity).
The additional solutions we constructed describe, respectively: asymptotically
(A)dS rotating black holes in arbitrary dimension, asymptotically (A)dS
charged black holes in arbitrary dimension, as well as asymptotically (A)dS
charged _and_ rotating black holes in $D=4$ (also in $D=3$, although we did
not write it down explicitly). We suspect that the hairy black hole solutions
in bigravity also carry across to the multi-metric theory in a natural way,
despite not performing the explicit calculation here. Furthermore, we related
these multi-metric black hole solutions to the well-known black string
solutions of higher dimensional GR, with the structure of the conformal
factors that defines the multi-metric vacuum encoding information about the
geometry of the extra dimension (we used the example of clockwork gravity to
represent a warped extra dimension in the dimensional deconstruction limit,
for example).

We later studied the linear stability of these multi-metric black holes. After
linearising the general theory to determine the dynamics of the spin-2 mass
modes, we showed that the GL and superradiant instabilities that plague
4-dimensional proportional black hole solutions in dRGT massive gravity and
bigravity carry over naturally to the multi-metric theory (as they should,
given the relation to black strings). More precisely: in dRGT/bigravity, in
terms of the one graviton mass, it is known that a bound exists at which the
GL instability turns on, as well as a regime for which the superradiant
instability is most efficient, both occuring when $mr_{s}\sim\mathcal{O}(1)$
(differing slightly between the Schwarzchild and Kerr cases). In the multi-
metric theory, this relation holds for _every_ graviton in the spectrum, which
translates to a stability bound on only the lightest massive state for the GL
mode, as well as a wider array of black hole masses that are potentially
affected significantly by the superradiant modes, relative to the situation in
dRGT/bigravity.

Since observations favour either very light or very heavy graviton masses (in
order for the multi-metric theory to agree with GR), the consequences of these
instabilities can be vastly different, and can affect vastly varying sizes of
black hole depending on the particular theory one chooses to work with (i.e.
depending on a particular choice of interaction coefficients and number of
metrics). In order to pin these consequences down, more work is required to
understand how the instabilities saturate. This will inevitably rely on us
having well-posed dynamical simulations within the framework of ghost free
multi-gravity, upon which the development is ongoing. Nevertheless, we hope
that as a result of this work we are now a small step closer to understanding
such interesting questions.

## Acknowledgements

K.W. is supported by a UK Science and Technology Facilities Council
studentship. P.M.S. and A.A. are supported by a STFC Consolidated Grant [Grant
No. ST/T000732/1]. Some calculations involving various vielbein contractions,
and checks thereof, were aided by use of the _xAct_ Mathematica package suite
[159] (specifically _xCoba_). For the purpose of open access, the authors have
applied a Creative Commons Attribution (CC BY) licence to any Author Accepted
Manuscript version arising.

## Data Access Statement

No new data were created or analysed in this study.

##  Appendix A Derivation of the field equations

To derive the vielbein form of the field equations as given in section II, we
determine the variation of the action, Eq. (1), using the differential form
framework prescribed in [160]. For an equivalent derivation in the metric
formalism, we refer the reader to e.g. [117].

The variation of the kinetic term with respect to the $i$-th tetrad gives:

$\displaystyle\delta S_{K}$
$\displaystyle=\frac{M_{i}^{D-2}}{2}\int_{\mathcal{M}_{D}}\delta
e^{(i)a}\wedge R^{(i)}_{bc}\wedge{\star}^{(i)}e^{(i)bc}_{\;\;\;\;\;\;\;a}$
(121) $\displaystyle=-M_{i}^{D-2}\int_{\mathcal{M}_{D}}\delta
e^{(i)a}\wedge{\star}^{(i)}G^{(i)}_{a}\;,$ (122)

where ${\star}G^{a}=-\frac{1}{2}R_{bc}\wedge{\star}e^{abc}$ is the Hodge dual
of the Einstein (vector-valued) 1-form whose components are those of the usual
Einstein tensor i.e. $G^{a}=G^{a}_{\;b}e^{b}$. To see that this is indeed just
the usual Einstein tensor, consider the component expression for the dual (a
$(D-1)$-form):

${\star}G^{a}=-\frac{1}{4(D-3)!}R_{bcmn}\varepsilon^{abc}_{\;\;\;\;d\ldots}e^{mnd\ldots}\;,$
(123)

from which we can use the fact that
$e_{k}\wedge{\star}G^{a}=G^{a}_{\;k}{\star}1$ to extract the components.
First, we have:

$\displaystyle e_{k}\wedge{\star}G^{a}$
$\displaystyle=-\frac{1}{4(D-3)!}R_{bcmn}\varepsilon^{abc}_{\;\;\;\;d\ldots}e_{k}^{\;mnd\ldots}$
$\displaystyle=+\frac{1}{4(D-3)!}R^{mn}_{\;\;\;\;bc}\varepsilon^{abcd\ldots}\varepsilon_{kmnd\ldots}{\star}1$
$\displaystyle=-\frac{1}{4}R^{mn}_{\;\;\;\;bc}\delta^{abc}_{kmn}{\star}1\;,$

using the symmetries of the Riemann tensor and the fact that the epsilon
tensors sum to

$\varepsilon_{\mu_{1}\mu_{2}\mu_{3}\ldots}\varepsilon^{\nu_{1}\nu_{2}\nu_{3}\ldots}=(-1)^{s}\delta^{\nu_{1}\nu_{2}\nu_{3}\ldots}_{\mu_{1}\mu_{2}\mu_{3}\ldots}\;,$
(124)

for metric of signature $(s,D-s)$, where
$\delta^{\nu_{1}\nu_{2}\nu_{3}\ldots}_{\mu_{1}\mu_{2}\mu_{3}\ldots}$ is the
generalised Kronecker delta defined by

$\delta^{\nu_{1}\ldots\nu_{p}}_{\mu_{1}\ldots\mu_{p}}=p!\delta^{\nu_{1}}_{[\mu_{1}}\ldots\delta^{\nu_{p}}_{\mu_{p}]}\;.$
(125)

We see that the Einstein tensor has (tetrad basis) components given by:

$G^{a}_{\;k}=-\frac{1}{4}R^{mn}_{\;\;\;\;bc}\delta^{abc}_{kmn}\;,$ (126)

or in more familiar form:

$\begin{split}G^{a}_{\;k}&=-\frac{1}{4}R^{mn}_{\;\;\;\;bc}3!\delta^{a}_{[k}\delta^{b}_{m}\delta^{c}_{n]}\\\
&=-\frac{1}{4}\left(-4R^{a}_{\;k}+2\delta^{a}_{k}R\right)\\\
&=R^{a}_{\;k}-\frac{1}{2}\delta^{a}_{k}R\;.\end{split}$ (127)

Now, for the potential term, the variation is:

$\displaystyle\delta S_{V}$ $\displaystyle=-\int_{\mathcal{M}_{D}}\delta
e^{(i)a}\wedge{\star}^{(i)}W^{(i)}_{a}\;,$ (128)

in which we define the $W$-tensor 1-form in an analogous way to our Einstein
1-forms above. Explicitly,

$\displaystyle{\star}^{(i)}$ $\displaystyle
W^{(i)}_{a}=\varepsilon_{ab_{1}\ldots b_{D-1}}$ (129)
$\displaystyle\times\sum_{j_{1}\ldots j_{D-1}}\mathcal{P}(i)T_{ij_{1}\ldots
j_{D-1}}e^{(j_{1})b_{1}}\wedge\ldots\wedge e^{(j_{D-1})b_{D-1}}\;.$

To extract the components, we use the same trick as before – namely, that
$\differential
x^{\mu}\wedge{\star}^{(i)}W^{(i)}_{a}=W^{(i)\mu}_{\;\;\;\;\;\;a}{\star}^{(i)}1$,
only this time wedging with the coordinate basis 1-form $\differential
x^{\mu}$, since the tetrads in Eq. (129) do not necessarily belong to the same
geometry (in order to make the calculation tractable, we need to use the
vielbeins to express everything in a common basis in which we can then
identify the volume form ${\star}^{(i)}1$; we make the natural choice of the
coordinate basis). The result is that:

$\displaystyle W^{(i)\mu}_{\;\;\;\;\;\;a}{\star}^{(i)}1$
$\displaystyle=\varepsilon_{ab_{1}\ldots b_{D-1}}\sum_{j_{1}\ldots
j_{D-1}}\mathcal{P}(i)T_{ij_{1}\ldots
j_{D-1}}e_{\lambda_{1}}^{({j_{1}}){b_{1}}}\ldots
e_{\lambda_{D-1}}^{({j_{D-1}}){b_{D-1}}}\differential
x^{\mu\lambda_{1}\ldots\lambda_{D-1}}$
$\displaystyle=-\varepsilon^{\mu\lambda_{1}\ldots\lambda_{D-1}}\varepsilon_{ab_{1}\ldots
b_{D-1}}\sum_{j_{1}\ldots j_{D-1}}\mathcal{P}(i)T_{ij_{1}\ldots
j_{D-1}}e_{\lambda_{1}}^{({j_{1}}){b_{1}}}\ldots
e_{\lambda_{D-1}}^{({j_{D-1}}){b_{D-1}}}{\star}^{(i)}1$
$\displaystyle=D!e_{\;\;\;\;\;{[a}}^{({i}){\mu}}e_{\;\;\;\;\;{b_{1}}}^{({i}){\lambda_{1}}}\ldots
e_{\;\;\;\;\;{b_{D-1}]}}^{({i}){\lambda_{D-1}}}\sum_{j_{1}\ldots
j_{D-1}}\mathcal{P}(i)T_{ij_{1}\ldots
j_{D-1}}e_{\lambda_{1}}^{({j_{1}}){b_{1}}}\ldots
e_{\lambda_{D-1}}^{({j_{D-1}}){b_{D-1}}}{\star}^{(i)}1\;.$

Contracting with $e^{(i)a}_{\nu}$ gives us the appropriate coordinate basis
expression, Eq. (16).

Finally, the variation of the matter action defines the energy-momentum 1-form
as:

$\delta S_{M}=\int_{\mathcal{M}_{D}}\delta
e^{(i)a}\wedge{\star}^{(i)}T^{(i)a}\;.$ (130)

Putting the three variations together and taking the functional derivative
with respect to $\delta e^{(i)a}$ gives us the differential form version of
our Einstein equations:

${\star}^{(i)}\left(-M^{D-2}_{i}G^{(i)}_{a}-W_{a}^{(i)}+T_{a}^{(i)}\right)=0\;,$
(131)

which in coordinate basis components becomes precisely the equations of motion
(15).

##  Appendix B Derivation of the linearised field equations

This time, the calculation is simpler in the metric formalism, and follows
closely the bigravity derivation in [117]. First, recall that the metric form
of the $W$-tensor for chain type interactions reads (C.F. Eq. (64)):

$\begin{split}W^{(i)\mu}_{\;\;\;\;\;\;\nu}&=\sum_{m=0}^{D}(-1)^{m}\beta_{m}^{(i,i+1)}Y_{(m)\nu}^{\mu}(S_{i\rightarrow
i+1})\\\
&+\sum_{m=0}^{D}(-1)^{m}\beta_{D-m}^{(i-1,i)}Y_{(m)\nu}^{\mu}(S_{i\rightarrow
i-1})\;.\end{split}$ (132)

To see that this is entirely equivalent to the correpsonding vielbein
formalism result, Eq. (26), notice that the building block matrices with the
proportional ansatz ($g^{(i)}_{\mu\nu}=a_{i}^{2}\bar{g}_{\mu\nu}$) take the
simple form $S_{i\rightarrow i\pm 1}=(a_{i\pm 1}/a_{i})\mathbbm{1}$, which
means that the eigenvalues of $S_{i\rightarrow i\pm 1}$ are simply $D$ copies
of $(a_{i\pm 1}/a_{i})$. Therefore, the elementary symmetric polynomials are:

$e_{k}(S_{i\rightarrow i\pm 1})=\left(\frac{a_{i\pm
1}}{a_{i}}\right)^{k}\binom{D}{k}\;.$ (133)

Substituting into Eq. (18), and using the binomial coefficient identity
$\sum_{k=0}^{m}(-1)^{k}\binom{D}{k}=(-1)^{m}\binom{D-1}{m}$, one finds that:

$Y_{(m)}(S_{i\rightarrow i\pm 1})=(-1)^{m}a_{i\pm
1}^{m}a_{i}^{-m}\binom{D-1}{m}\mathbbm{1}\;,$ (134)

which recovers Eq. (26) upon substitution back into Eq. (132).

To linearise the system, we perturb around our proportional background. The
metric and its inverse are expanded as:

$\displaystyle g^{(i)}_{\mu\nu}$
$\displaystyle=a_{i}^{2}\bar{g}_{\mu\nu}+\delta g^{(i)}_{\mu\nu}$ (135)
$\displaystyle g^{(i)\mu\nu}$
$\displaystyle=a_{i}^{-2}\bar{g}^{\mu\nu}-a_{i}^{-4}\delta g^{(i)\mu\nu}\;,$
(136)

where we have included factors of $a_{i}$ so that $\delta g^{(i)}_{\mu\nu}$
behaves as if it lived in the _common_ background described by metric
$\bar{g}_{\mu\nu}$ (if $\delta g^{(i)}_{\mu\nu}$ had its indices instead
manipulated with $a_{i}^{2}\bar{g}_{\mu\nu}$, the inverse would simply be
$g^{(i)\mu\nu}=a_{i}^{-2}\bar{g}^{\mu\nu}-\delta g^{(i)\mu\nu}$).

It is well known (see e.g. [34]) that the Einstein tensor linearises to the
Lichnerowicz operator acting on the perturbation, that is:

$\delta G^{(i)}_{\mu\nu}=\bar{\mathcal{E}}^{\alpha\beta}_{\;\mu\nu}\delta
g^{(i)}_{\alpha\beta}\;,$ (137)

so we shall skip this part of the derivation here and focus on the potential.
To that end, the first order variation of the $W$-tensor is:

$\begin{split}\delta W^{(i)\mu}_{\;\;\;\;\;\;\nu}=\delta
g^{(i)\mu}_{\;\;\;\;\;\lambda}&\bar{W}^{(i)\lambda}_{\;\;\;\;\;\;\nu}+\sum_{m=0}^{D}(-1)^{m}\beta_{m}^{(i,i+1)}\delta
Y_{(m)\nu}^{\mu}(S_{i\rightarrow i+1})\\\
&+\sum_{m=0}^{D}(-1)^{m}\beta_{D-m}^{(i-1,i)}\delta
Y_{(m)\nu}^{\mu}(S_{i\rightarrow i-1})\;,\end{split}$ (138)

where, around the proportional background, we know that
$\bar{W}^{(i)\lambda}_{\;\;\;\;\;\;\nu}=\delta^{\lambda}_{\nu}M_{i}^{D-2}\bar{\Lambda}/a_{i}^{2}$
(see section III.1). The variation of the $Y$’s is given by:

$\displaystyle\delta Y_{(m)}(S)=\sum_{k=0}^{m}(-1)^{k}$
$\displaystyle\left[(m-k)S^{m-k-1}e_{k}(S)\delta S\right.$ (139)
$\displaystyle\left.+S^{m-k}\delta e_{k}(S)\right]\;,$

where, by virtue of Eq. (11), we have:

$\delta e_{k}(S)=-\sum_{n=1}^{k}(-1)^{n}\Tr(S^{n-1}\delta S)e_{k-n}(S)\;.$
(140)

Since it is built from the metric variations, $\delta S_{\mu\nu}$ also has its
indices manipulated with the common background metric $\bar{g}_{\mu\nu}$.

Eqs. (138)-(140) hold in general, but around the proportional background
things simplify greatly: as we saw, the building block matrices take the
simple form $S_{i\rightarrow i\pm 1}=(a_{i\pm 1}/a_{i})\mathbbm{1}$. After
substituting this in above (and using some properties of the binomial
coefficients) one can show that the $Y$ variations become [117]:

$\displaystyle\delta Y_{(m)}(S_{i\rightarrow i\pm 1})=$
$\displaystyle(-1)^{m}a_{i\pm 1}^{m-1}a_{i}^{1-m}\binom{D-2}{m-1}$ (141)
$\displaystyle\times\left[\Tr(\delta S_{i\rightarrow i\pm
1})\mathbbm{1}-\delta S_{i\rightarrow i\pm 1}\right]\;,$

and so the variation of the $W$-tensor is:

$\begin{split}\delta
W^{(i)\mu}_{\;\;\;\;\;\;\nu}=\frac{M_{i}^{D-2}\bar{\Lambda}}{a_{i}^{2}}&\delta
g^{(i)\mu}_{\;\;\;\;\;\lambda}+\Sigma_{i}^{(+)}\left[\delta S_{i\rightarrow
i+1}\delta^{\mu}_{\nu}-(\delta S_{i\rightarrow i+1})^{\mu}_{\;\nu}\right]\\\
&+\Sigma_{i}^{(-)}\left[\delta S_{i\rightarrow i-1}\delta^{\mu}_{\nu}-(\delta
S_{i\rightarrow i-1})^{\mu}_{\;\nu}\right]\;,\end{split}$ (142)

where $\Sigma_{i}^{(\pm)}$ are precisely as in Eqs. (66) and (67).

All that remains, in order to determine the linearised equations of motion, is
to calculate the precise form of $\delta S_{i\rightarrow i\pm 1}$. By starting
with $(S_{i\rightarrow i\pm 1}^{2})^{\mu}_{\;\nu}=g^{(i)\mu\lambda}g^{(i\pm
1)}_{\lambda\nu}$ and substituting in Eqs. (135) and (136) for the perturbed
metrics, one can show that the desired variation is given by:

$(\delta S_{i\rightarrow i\pm
1})_{\mu\nu}=\frac{1}{2a_{i}^{2}}\frac{a_{i}}{a_{i\pm 1}}\left[\delta g^{(i\pm
1)}_{\mu\nu}-\left(\frac{a_{i\pm 1}}{a_{i}}\right)^{2}\delta
g^{(i)}_{\mu\nu}\right]\;.$ (143)

With this, the linearised (vacuum) equations take the following form:

$\begin{split}\bar{\mathcal{E}}^{\alpha\beta}_{\;\mu\nu}\delta
g^{(i)}_{\alpha\beta}+\bar{\Lambda}\delta
g_{\mu\nu}^{(i)}+\frac{a_{i}^{2}}{M_{i}^{D-2}}&\Bigg{\\{}\Sigma_{i}^{(+)}\left[\frac{a_{i+1}}{a_{i}}\left(\delta
g^{(i)}_{\mu\nu}-\bar{g}_{\mu\nu}\delta
g^{(i)}\right)-\frac{a_{i}}{a_{i+1}}\left(\delta
g^{(i+1)}_{\mu\nu}-\bar{g}_{\mu\nu}\delta g^{(i+1)}\right)\right]\\\
&+\Sigma_{i}^{(-)}\left[\frac{a_{i-1}}{a_{i}}\left(\delta
g^{(i)}_{\mu\nu}-\bar{g}_{\mu\nu}\delta
g^{(i)}\right)-\frac{a_{i}}{a_{i-1}}\left(\delta
g^{(i-1)}_{\mu\nu}-\bar{g}_{\mu\nu}\delta
g^{(i-1)}\right)\right]\Bigg{\\}}=0\;,\end{split}$ (144)

One can check for bigravity, where only the $i=0$ and $i=1$ terms are present,
that these equations reduce to precisely the linearised equations of [117].

If one instead parametrises their perturbations as

$\delta
g^{(i)}_{\mu\nu}=\frac{a_{i}}{M_{i}^{\frac{D-2}{2}}}h^{(i)}_{\mu\nu}\;,$ (145)

as we did in section IV, then one recovers our linearised equations (95), as
well as our expressions for the mass matrix components, Eqs. (97) and (98).
The mass matrix for the $h^{(i)}_{\mu\nu}$ of course has different components
to the mass matrix for $\delta g^{(i)}_{\mu\nu}$, but since their respective
equations are just related by a simple field rescaling, they share the mass
same eigenvalues. Therefore, the equations for the mass modes that one obtains
after diagonalising are equivalent regardless of whether one initially uses
$\delta g^{(i)}_{\mu\nu}$ or $h^{(i)}_{\mu\nu}$ to express their
perturbations.

Finally, we note that if $\bar{g}_{\mu\nu}$ is the $D$-dimensional (A)dS
metric, then one may additionally pull out the cosmological constant from the
$\bar{R}^{\alpha\;\beta}_{\;\mu\;\nu}$ background curvature contribution to
$\bar{\mathcal{E}}^{\alpha\beta}_{\;\mu\nu}$, to express the first two terms
in Eq. (144) as:

$\tilde{\mathcal{E}}^{\alpha\beta}_{\;\mu\nu}\delta
g^{(i)}_{\mu\nu}-\frac{2\bar{\Lambda}}{D-2}\left(\delta
g^{(i)}_{\mu\nu}-\frac{1}{2}\bar{g}_{\mu\nu}\delta g^{(i)}\right)+\ldots\;,$
(146)

where $\tilde{\mathcal{E}}^{\alpha\beta}$ now contains only the covariant
derivative operators, if one wishes to match up exactly with the corresponding
equations as written in [117, 39] for bigravity. For us, it is important to
keep the background curvature in explicitly, since it feeds into our
discussion of black hole stability, when $\bar{g}_{\mu\nu}$ is either the
Schwarzchild or Kerr metric.

##  Appendix C Star type interactions

Throughout this work, we have consistently worked with chain type
interactions, due to their nice interpretation as arising from dimensional
deconstruction. As mentioned briefly in section II, however, the chain type
interaction is not the only one that is devoid of the Boulware-Deser ghost:
‘star’ type interactions, where many metrics couple to one common central
metric but not to each other, are also allowed, as well as arbitrary
combinations of both stars and chains provided that no interaction cycles form
between the metrics.

We would like to see how our calculations and results are altered for star
type interactions. Of course, the only difference between multi-metric gravity
with chain type interactions and multi-metric gravity with star type
interactions lies in the structure of the interaction coefficients
$T_{i_{1}\ldots i_{D}}$.

For the chain, only terms of the form $T_{iiii\ldots}$, $T_{i+1,iii\ldots}$,
$T_{i-1,iii\ldots}$, $T_{i+1,i+1,ii\ldots}$ etc. were allowed, corresponding
to the scenario where each metric interacts only with its nearest neighbours
on either side. This lead us to paramterise the interaction coefficients
nicely in terms of the $\beta_{m}^{(i,i+1)}$ of Eqs. (13) and (14).

For the star, the situation is slightly different. If we let the index $i=0$
label the common central metric and the indices $j=(1,\ldots,N-1)$ label the
outer metrics, which each couple only to $g^{(0)}_{\mu\nu}$, then the only
permitted interaction coefficients are those that involve combinations of $0$
with any one distinguished $j$ e.g. $T_{0000\ldots}$, $T_{0111\ldots}$,
$T_{0022\ldots}$, $T_{3333\ldots}$ and so on888One cannot allow, for example,
$T_{012\ldots}$, as this would imply an interaction cycle between
$g^{(0)}_{\mu\nu}$, $g^{(1)}_{\mu\nu}$ and $g^{(2)}_{\mu\nu}$, which as we
said in section II renders the theory ghostly.. One may parametrise the
allowed coefficients in the following manner:

$\displaystyle D!\,T_{0000\ldots 0}$
$\displaystyle=\sum_{j=1}^{N-1}\beta_{0}^{(0,j)}\;,$ (147) $\displaystyle
D!\,T_{\\{j\\}^{m}\\{0\\}^{D-m}}$ $\displaystyle=\beta_{m}^{(0,j)}\;,$ (148)

where $m=1,\ldots,D$, the form of the coefficients now reflecting the fact
that the central $0$-th metric interacts with all $N-1$ outer metrics while
the outer $j$-th metrics interact only with the central one, and that all
interactions are oriented outward from the central metric. These
$\beta_{m}^{(0,j)}$ again coincide with those of the metric formalism. With
this parametrisation, one can substitute into Eq. (16) to determine the form
of the $W$-tensors for a given set of vielbeins.

For the proportional metric ansatz, the result is that the $0$-th $W$-tensor
has components:

$W^{(0)\mu}_{\;\;\;\;\;\;\nu}=\delta^{\mu}_{\nu}\sum_{j=1}^{N-1}\sum_{m=0}^{D}\beta_{m}^{(0,j)}\binom{D-1}{m}a_{j}^{m}a_{0}^{-m}\;,$
(149)

while the corresponding expression for the $j$-th $W$-tensor is:

$W^{(j)\mu}_{\;\;\;\;\;\;\;\;\nu}=\delta^{\mu}_{\nu}\sum_{m=0}^{D}\beta_{D-m}^{(0,j)}\binom{D-1}{m}a_{0}^{m}a_{j}^{-m}\;.$
(150)

Again, the same arguments regarding the solvability of the equations of motion
for chain type interactions apply here. Namely, one requires that the
cosmological constant contribution must be the same for each $W$-tensor:

$\displaystyle\sum_{j=1}^{N-1}\Lambda_{j}^{(+)}$
$\displaystyle=\bar{\Lambda}\;,$ (151) $\displaystyle\Lambda^{(-)}_{j}$
$\displaystyle=\bar{\Lambda}\;,$ (152)

defining $\Lambda_{j}^{(\pm)}$ as ($a_{i}^{2}/M_{i}^{D-2}$ times) the sums
involving $\beta^{(0,j)}_{m}$ and $\beta_{D-m}^{(0,j)}$, respectively. The
equations then reduce as before to $N$ copies of Einstein’s equations for the
$\bar{g}_{\mu\nu}$ common background, which one may in principle solve for
$\bar{\Lambda}$ and the $N-1$ free conformal factors, after fixing one of them
via coordinate rescaling.

For the non-proportional metric ansatze, the surviving off-diagonal parts of
the $W$-tensors are also patterned as above, reflecting this new structure of
the interactions. Explicitly, if one defines (C.F. Eqs. (66) and (67)):

$\displaystyle\Sigma_{j}^{(+)}$
$\displaystyle=\sum_{m=0}^{D}\beta_{m}^{(0,j)}\binom{D-2}{m-1}a_{j}^{m-1}a_{0}^{1-m}\;,$
(153) $\displaystyle\Sigma_{j}^{(-)}$
$\displaystyle=\sum_{m=0}^{D}\beta_{D-m}^{(0,j)}\binom{D-2}{m-1}a_{0}^{m-1}a_{j}^{1-m}\;,$
(154)

related as:

$\Sigma_{j}^{(-)}=\left(\frac{a_{0}}{a_{j}}\right)^{D-2}\Sigma_{j}^{(+)}\;,$
(155)

then the star type analogue of (say) Eq. (65) for the chain interactions is:

$\displaystyle\begin{split}W^{(0)\mu}_{\;\;\;\;\;\;\nu}&=\frac{M_{0}^{D-2}\delta^{\mu}_{\nu}}{a_{0}^{2}}\sum_{j=1}^{N-1}\Lambda_{j}^{(+)}\\\
&+\frac{l^{\mu}l_{\nu}}{2U}\sum_{j=1}^{N-1}\frac{a_{j}}{a_{0}}\left(r_{s,0}-r_{s,j}\right)\Sigma_{j}^{(+)}\;,\end{split}$
(156) $\displaystyle W^{(j)\mu}_{\;\;\;\;\;\;\nu}$
$\displaystyle=\frac{M_{j}^{D-2}\delta^{\mu}_{\nu}}{a_{j}^{2}}\Lambda_{j}^{(-)}+\frac{l^{\mu}l_{\nu}}{2U}\frac{a_{0}}{a_{j}}\left(r_{s,j}-r_{s,0}\right)\Sigma_{j}^{(-)}\;,$
(157)

and similarly for the charged variants of the above (C.F. section III.2.3). In
principle, solutions then exist in each of the proportional (all
$r_{s,j}=r_{s,0}$), non-proportional (all $\Sigma^{(+)}_{j}=0$) and partially
proportional (combinations of both) branches.

As for the linearised equations, the star type analogue of Eq. (144) for the
metric perturbations is:

$\displaystyle\bar{\mathcal{E}}^{\alpha\beta}_{\;\mu\nu}\delta
g^{(0)}_{\alpha\beta}+\bar{\Lambda}\delta
g_{\mu\nu}^{(0)}+\frac{a_{0}^{2}}{M_{0}^{D-2}}\sum_{j=1}^{N-1}\Bigg{\\{}\Sigma_{j}^{(+)}\left[\frac{a_{j}}{a_{0}}\left(\delta
g^{(0)}_{\mu\nu}-\bar{g}_{\mu\nu}\delta
g^{(0)}\right)-\frac{a_{0}}{a_{j}}\left(\delta
g^{(j)}_{\mu\nu}-\bar{g}_{\mu\nu}\delta g^{(j)}\right)\right]\Bigg{\\}}$
$\displaystyle=0\;,$ (158)
$\displaystyle\bar{\mathcal{E}}^{\alpha\beta}_{\;\mu\nu}\delta
g^{(j)}_{\alpha\beta}+\bar{\Lambda}\delta
g_{\mu\nu}^{(j)}+\frac{a_{j}^{2}}{M_{j}^{D-2}}\Sigma_{j}^{(-)}\left[\frac{a_{0}}{a_{j}}\left(\delta
g^{(j)}_{\mu\nu}-\bar{g}_{\mu\nu}\delta
g^{(j)}\right)-\frac{a_{j}}{a_{0}}\left(\delta
g^{(0)}_{\mu\nu}-\bar{g}_{\mu\nu}\delta g^{(0)}\right)\right]$
$\displaystyle=0\;.$ (159)

One may check, for example, that in $D=4$ spacetime dimensions, and denoting
$a_{0}=1,a_{j}=c_{j},\gamma_{j}=M_{j}/M_{0}$, these equations imply a mass
matrix for the $\delta g^{(i)}_{\mu\nu}$ that takes the form:

$\mathcal{M}^{2}=\frac{1}{M_{0}^{2}}\begin{bmatrix}\sum_{j=1}^{N-1}\Sigma_{j}^{(+)}c_{j}&-\frac{\Sigma_{1}^{(+)}}{c_{1}}&-\frac{\Sigma_{2}^{(+)}}{c_{2}}&\ldots&-\frac{\Sigma_{N-1}^{(+)}}{c_{N-1}}\\\
-\frac{\Sigma_{1}^{(+)}c_{1}}{\gamma_{1}^{2}}&\frac{\Sigma_{1}^{(+)}}{c_{1}\gamma_{1}^{2}}&0&\ldots&0\\\
-\frac{\Sigma_{2}^{(+)}c_{2}}{\gamma_{2}^{2}}&0&\frac{\Sigma_{2}^{(+)}}{c_{2}\gamma_{2}^{2}}&\ldots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
-\frac{\Sigma_{N-1}^{(+)}c_{N-1}}{\gamma_{N-1}^{2}}&0&0&\ldots&\frac{\Sigma_{N-1}^{(+)}}{c_{N-1}\gamma_{N-1}^{2}}\end{bmatrix}$
(160)

This is precisely the star type interaction mass matrix that was derived for
the $D=4$ multi-metric theory in [39].

## References

  * [1] S. N. Gupta, Gravitation and electromagnetism, Physical Review 96 (1954), no. 6 1683.
  * [2] S. Weinberg, Photons and gravitons in perturbation theory: Derivation of maxwell’s and einstein’s equations, Physical Review 138 (1965), no. 4B B988.
  * [3] S. Deser, Self-interaction and gauge invariance, General Relativity and gravitation 1 (1970), no. 1 9–18.
  * [4] R. Feynman, Feynman lectures on gravitation.  CRC Press, 2018.
  * [5] D. G. Boulware and S. Deser, Classical general relativity derived from quantum gravity, Annals of Physics 89 (1975), no. 1 193–240.
  * [6] A. G. Riess, A. V. Filippenko, P. Challis, A. Clocchiatti, A. Diercks, P. M. Garnavich, R. L. Gilliland, C. J. Hogan, S. Jha, R. P. Kirshner, et al., Observational evidence from supernovae for an accelerating universe and a cosmological constant, The astronomical journal 116 (1998), no. 3 1009, [astro-ph/9805201].
  * [7] S. Perlmutter, G. Aldering, G. Goldhaber, R. Knop, P. Nugent, P. G. Castro, S. Deustua, S. Fabbro, A. Goobar, D. E. Groom, et al., Measurements of $\omega$ and $\lambda$ from 42 high-redshift supernovae, The Astrophysical Journal 517 (1999), no. 2 565, [astro-ph/9812133].
  * [8] E. Babichev, L. Marzola, M. Raidal, A. Schmidt-May, F. Urban, H. Veermäe, and M. von Strauss, Heavy spin-2 dark matter, Journal of Cosmology and Astroparticle Physics 2016 (2016), no. 09 016, [arXiv:1607.03497].
  * [9] E. Babichev, L. Marzola, M. Raidal, A. Schmidt-May, F. Urban, H. Veermäe, and M. von Strauss, Bigravitational origin of dark matter, Physical Review D 94 (2016), no. 8 084055, [arXiv:1604.08564].
  * [10] N. G. Albornoz, A. Schmidt-May, and M. von Strauss, Dark matter scenarios with multiple spin-2 fields, Journal of Cosmology and Astroparticle Physics 2018 (2018), no. 01 014, [arXiv:1709.05128].
  * [11] L. Marzola, M. Raidal, and F. R. Urban, Oscillating spin-2 dark matter, Physical Review D 97 (2018), no. 2 024010, [arXiv:1708.04253].
  * [12] F. Niedermann, A. Padilla, and P. M. Saffin, Higher order clockwork gravity, Physical Review D 98 (2018), no. 10 104014, [arXiv:1805.03523].
  * [13] A. Avgoustidis, F. Niedermann, A. Padilla, and P. M. Saffin, Deconstructing higher order clockwork gravity, Physical Review D 103 (2021), no. 12 124007, [arXiv:2010.10970].
  * [14] K. Wood, P. M. Saffin, and A. Avgoustidis, Clockwork cosmology, Journal of Cosmology and Astroparticle Physics 2023 (2023), no. 7 062, [arXiv:2304.09205].
  * [15] N. Boulanger, T. Damour, L. Gualtieri, and M. Henneaux, Inconsistency of interacting, multi-graviton theories, Nuclear Physics B 597 (2001), no. 1-3 127–171, [hep-th/0007220].
  * [16] M. Fierz and W. E. Pauli, On relativistic wave equations for particles of arbitrary spin in an electromagnetic field, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 173 (1939), no. 953 211–232.
  * [17] D. G. Boulware and S. Deser, Can gravitation have a finite range?, Phys. Rev. D 6 (Dec, 1972) 3368–3382.
  * [18] P. Van Nieuwenhuizen, On ghost-free tensor lagrangians and linearized gravitation, Nuclear Physics B 60 (1973) 478–492.
  * [19] C. de Rham and G. Gabadadze, Generalization of the fierz-pauli action, Phys. Rev. D 82 (Aug, 2010) 044020, [arXiv:1007.0443].
  * [20] C. de Rham, G. Gabadadze, and A. J. Tolley, Resummation of massive gravity, Physical Review Letters 106 (jun, 2011) [arXiv:1011.1232].
  * [21] S. F. Hassan and R. A. Rosen, On non-linear actions for massive gravity, Journal of High Energy Physics 2011 (2011), no. 7 1–22, [arXiv:1103.6055].
  * [22] S. F. Hassan and R. A. Rosen, Resolving the ghost problem in nonlinear massive gravity, Physical review letters 108 (2012), no. 4 041101, [arXiv:1106.3344].
  * [23] S. F. Hassan, R. A. Rosen, and A. Schmidt-May, Ghost-free massive gravity with a general reference metric, Journal of High Energy Physics 2012 (2012), no. 2 1–26, [arXiv:1109.3230].
  * [24] S. Hassan, A. Schmidt-May, and M. von Strauss, Proof of consistency of nonlinear massive gravity in the stückelberg formulation, Physics Letters B 715 (2012), no. 4-5 335–339, [arXiv:1203.5283].
  * [25] A. Golovnev, On the hamiltonian analysis of non-linear massive gravity, Physics Letters B 707 (2012), no. 3-4 404–408, [arXiv:1112.2134].
  * [26] J. Klusoň, Nonlinear massive gravity with additional primary constraint and absence of ghosts, Physical Review D 86 (2012), no. 4 044024, [arXiv:1204.2957].
  * [27] J. Klusoň, Note about hamiltonian formalism for general nonlinear massive gravity action in stückelberg formalism, International Journal of Modern Physics A 28 (2013), no. 30 1350160, [arXiv:1209.3612].
  * [28] T. Kugo and N. Ohta, Covariant approach to the no-ghost theorem in massive gravity, Progress of Theoretical and Experimental Physics 2014 (2014), no. 4 043B04, [arXiv:1401.3873].
  * [29] N. Arkani-Hamed, H. Georgi, and M. D. Schwartz, Effective field theory for massive gravitons and gravity in theory space, Annals of Physics 305 (2003), no. 2 96–118, [hep-th/0210184].
  * [30] P. Creminelli, A. Nicolis, M. Papucci, and E. Trincherini, Ghosts in massive gravity, Journal of High Energy Physics 2005 (2005), no. 09 003, [hep-th/0505147].
  * [31] S. F. Hassan and R. A. Rosen, Bimetric gravity from ghost-free massive gravity, Journal of High Energy Physics 2012 (2012), no. 2 1–12, [arXiv:1109.3515].
  * [32] S. F. Hassan and R. A. Rosen, Confirmation of the secondary constraint and absence of ghost in massive gravity and bimetric gravity, Journal of High Energy Physics 2012 (2012), no. 4 1–16, [arXiv:1111.2070].
  * [33] K. Hinterbichler and R. A. Rosen, Interacting spin-2 fields, Journal of High Energy Physics 2012 (2012), no. 7 1–34, [arXiv:1203.5783].
  * [34] C. de Rham, Massive gravity, Living reviews in relativity 17 (2014), no. 1 1–189, [arXiv:1401.4173].