Black Holes in Multi-Metric Gravity

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 nonproportional 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) Schwarzschild(-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 5dimensional 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 Schwarzschild 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 Schwarzschild 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 [13,64] to make 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 Schwarzschild and Kerr solutions to the general multi-metric theory; finally we conclude in section VI.
We work with natural units c = = 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 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.[12,14,33,68]): The tetrad basis 1-forms are e (i)a = e (i)a µ dx µ , where the indices run from 0 to D − 1, with the vielbeins defined through g ν η ab , and the shorthand e (i)ab... in the kinetic term means e (i)a ∧ e (i)b ∧ . ... 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 η ab and Greek indices with g µν , while we can swap between Latin and Greek indices using the vielbeins, via changes of basis.R (i) ab is the curvature 2form associated with the i-th (Levi-Civita) connection, with one index lowered by η 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 i1...iD = T (i1...iD ) 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 with the matrix square root defined in the sense that λν .Explicitly, the metric version of the potential term is: where the β are arbitrary constants related to the T i1...iD 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: . . .
They can also be constructed iteratively in terms of the trace of S, as: The structure of the building block matrices means that S i→j = S −1 j→i , 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→j (not S j→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 i1...iD .
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, 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 → 2 → 3 → 1), so such models are manifestly ghost free [39,75,76], 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 multimetric 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) 1 .Such a choice is natural if one thinks of the theory as arising from some sort of dimensional deconstruction [13,14,64].In terms of the T i1...iD , this means that one can only permit terms of the form T iiii... , T i+1,iii... , T i−1,iii... , T i+1,i+1,ii... 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 i1...iD coefficients.For chain type interactions, one can always neatly parametrise the T i1...iD in terms of some new constants β (i,i+1) m , which characterise the interactions between the i-th and (i + 1)-th metrics, as [14] where i = 0, . . ., N − 1 and m = 0, . . ., D, with β = 0 since on each of the two boundaries of the 1 In 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.
interaction chain the corresponding metrics have only one nearest neighbour.The factor of D! is included so that these β (i,i+1) m are then precisely the same β (i,j) m as in the metric formalism Eq. ( 5).The sense of interaction orientation is encoded in T iiii...i , where positively oriented interactions contribute a β 0 term while negatively oriented interactions contribute a β D term.In practice, we will often choose to write β (i,i+1) m = α i β m , where α −1 = α N −1 = 0 but the rest of the α 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 β 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 β (i,i+1) m 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: where ν 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)µ ν reads [14] (see appendix A for the derivation): with P(i) counting the number of times the index (i) appears in the interaction coefficients i.e. a term with T ij1...jD−1 has P(i) = 1, a term with T iij2...jD−1 has P(i) = 2, and so on.The equivalent metric formalism expression is: where (with respect to the i-th metric) j denote positively oriented interactions, k denote negatively oriented interactions, and we define The T (i)µ ν 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 [77,80,[82][83][84].
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]: Whenever matter couples to one site only, or when there is no matter coupling at all, the divergences of each Wtensor (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 [14,85] 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.

A. Proportional Solutions
We look for solutions to the multi-metric equations (15) where the metrics are conformally related to one another, where at this stage ḡµν is some arbitrary fixed metric common to all sites.Since all of the metrics live on the same manifold 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: Therefore, one may express the W -tensor in terms of ē a µ only, and so remove the vielbeins from the sum over the j's entirely: 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: and so the W -tensor takes the following simple, diagonal form, irrespective of the precise identity of ḡµν : ) Taking chain type interactions, using the symmetry of the T i1...iD coefficients, and parametrising them as in Eqs. ( 13) and (14), one finds that the components of the W -tensor are explicitly given by: The Bianchi constraint Eq. ( 19) forces all the a i to be constant [38].Lowering an index in Eq. ( 26) tells us that µν ; 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 µν (g) = G µν (ag), so we have that 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 ḡµν : where the effective cosmological constants Λ i are precisely the contributions from W (i)µ ν .Explicitly, the appropriate rescaled quantities are so that all terms in Eqs.(27) now behave as if they lived on the common background ḡµν i.e. have their indices manipulated with ḡµν .We can take differences of Eqs.(27) to show that the effective cosmological constants must satisfy: and similarly the energy-momentum tensors must all be proportional [38], 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 Λ'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 β's vanishes, so only one sum is present in the corresponding W -tensor.If we denote Λ As mentioned earlier, one may rescale the coordinates to fix the value of exactly one of the conformal factors.In vacuum, then, with all T (i)µ ν = 0, the equations of motion (given a choice for the interaction coefficients and provided that Ḡµ ν = − Λδ µ ν ) specify a system of N algebraic equations for N variables: namely, Λ, 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 i1...iD .
We further note that the splitting of the Λ's given in Eqs. ( 32)- (34) means that, unless Λ = 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 ∀ 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 ḡµν , 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.

Deconstructed black string in 4d
As a straightforward example, but one with an interesting physical interpration, one may take ḡµν to be the Schwarzschild metric.One then has the situation where the multi-gravity metrics are: with r s the Schwarzschild radius and dΩ 2  (2) the line element for the 2-sphere.This ansatz solves the vacuum equations, provided that the conformal factors are such that Λ = 0.
The physical meaning of the various conformal factors in this scenario is clear: suppose we make the coordinate changes r → r and η → t defined by r = r/a i and dη = dt/a i , then the i-th metric becomes: we see that an observer minimally coupled to the i-th metric would see a black hole whose Schwarzschild 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)µ ν = 0.This situation is reminiscent of the well-known black string solutions in higher dimensional gravity, where Schwarzschild 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 Schwarzschild is precisely the dimensional deconstruction [13,64] 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: with q 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).
With this vacuum structure, we see that the system is solved by a series of 4D Schwarzschild 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 Schwarzschild 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 [13,14] is a more complicated scalar-tensor braneworld, but the idea remains the same, with the continuum theory admitting Ricci-flat hypersurfaces.

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: where the F (i) = dA (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 µ and the corresponding energy-momentum tensors are: The common metric ḡµν is typically written in Boyer-Lindquist coordinates as [98][99][100]: where j is the rotation parameter2 and the various functions are defined as with r Q a scale related to the electric charge in a manner that will be fixed below.If Λ > 0, the metric is Kerr-Newman-dS, while if Λ < 0 the metric is Kerr-Newman-AdS.
With this choice for ḡµν , 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: provided that the charges are expressed in terms of r Q as One can check that the field equations for the electromagnetic fields, 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 , Λ, j) to 0. For example, taking r Q → 0 gives us the Kerr-(A)dS solution, if we also take j → 0 we get Schwarzschild-(A)dS, then taking Λ → 0 we recover the Schwarzschild (black string) solution from earlier.If instead we keep r Q but send j and Λ to 0, we get Reissner-Nordstrom.Finally, keeping Λ 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].

B. 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'.

Rotating (A)dS black holes in arbitrary dimension
We first look for solutions describing rotating but noncharged 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 [47,101]: Here, ḡµν is taken to be the metric of D-dimensional (A)dS space, φ 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: Some elaboration is required here: the coordinate system comprises a time coordinate t, a radial coordinate r, ⌊(D − 1)/2⌋ azimuthal coordinates ϕ k , and ⌊D/2⌋ coordinates µ k that satisfy The sums then run to n = ⌊D/2⌋, 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 ϕ n = dϕ n = 0.The j k are at this stage simply ⌊(D−1)/2⌋ 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: and the required null vector (and its dual 1-form) that is tangent to a null-geodesic congruence in this spacetime is: Finally, the scalar functions φ i are given by: with each metric now allowed its own independent Schwarzschild radius r s,i , and where the function U differs between the even and odd dimensional cases [47,101]: 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 ν = 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: 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 ḡµν and g (i) µν (i.e.l µ l µ = 0) means that its contribution to the interaction building block matrices S i→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→j is simply [47]: 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→j picks up only the contribution from δ µ ν , so the elementary symmetric polynomials are: Using this in Eq. ( 18), along with the binomial coefficient identities one finds that: 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→i+1 and S i−1→i , Eq. ( 17) involves only two sums: recalling from section II that S −1 i−1→i = S i→i−1 .With the Y (m) given by Eq. (63), and the φ i given by Eq. ( 55), the explicit expression for the components becomes: where we have defined: The part of the W -tensor proportional to δ µ ν is exactly as in Eq. ( 26) for the proportional solutions, but now there are additional off-diagonal terms proportional to l µ l ν .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 Schwarzschild radii to be the same, which just recovers the proportional solution from before.
The second is to make all of the Σ i 's vanish: These conditions are polynomial equations that fix the ratios of neighbouring conformal factors.Furthermore, due to the recurrence relation for the binomial coefficients: 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: 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 β 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 Λ = 0, one only needs to fix β 0 and β D in terms of the other β'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 β D−m in the (i + 1)-th off-diagonal W -tensor is proportional to the sum involving β m in the i-th off-diagonal W -tensor. Explicitly, one has that: so when N = 2 and the only terms present are Σ = 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 Λ = 03 .
In the multi-metric scenario, however, whenever one has N > 2, there are more Σ (+) 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 Σ (+) I−1 for some specific i = I, this implies only the vanishing of Σ (−) I and nothing else, so the requirement that the off-diagonal parts of W (I)µ ν vanish may still be satisfied by either r s,I = r s,I+1 or Σ (+) I = 0 (which is independent of Σ (+) 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 Σ (+) k = 0 and r s,k = r s,k+1 , for different k ⊂ i.This corresponds to the situation where some -but not all -of the metrics can be simultaneously diagonalised.Those sites where Σ (+) 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 If there are n total sites that have Σ (+) 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 Λ and the remaining N − n − 1 free conformal factors, as well as n equations that become constraints on the β (i,i+1) m 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 β (i,i+1) m 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 Schwarzschild radii are the same (the proportional solutions) or all the Σ 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 Σ (+) 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].

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 nonproportional generalisation of the Kerr-Newman-(A)dS metric from section III A 2. In the D = 4 case of Eq. (48), there is only one azimuthal coordinate ϕ 1 = ϕ, only one rotation parameter j 1 = j (hence one Ξ 1 = Ξ), and only two µ k coordinates, which we can parametrise without loss of generality as [101]: With these definitions, the (A)dS metric explicitly takes the form: where ∆ θ and ρ 2 are as in Eqs.(42) and (44).The null vector and 1-form are given by: while the function U becomes: As in section III A 2, to incorporate charge we must include copies of the Maxwell action for the matter sector, only now we allow the freedom to couple to only some subset of sites k ⊂ 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 φ i to [47]: where each metric is now allowed an independent r Q,i as well as r s,i .The coordinate transformation that retrieves the Boyer-Lindquist form of the metric, Eq. ( 40), from the Kerr-Schild form is given in [101].We note that the extension of φ i to include charge while keeping rotation is particularly simple in D = 4; this is assuredly not the case in higher dimensions, where even in standard GR the analogue of the Kerr-Newman metric for D > 4 is not known, so there is no safe starting point to begin to look for the corresponding multi-metric solutions.This is why we treat the charged and rotating cases separately outside of D = 4.
To account for the additional contributions to the Einstein tensors on the charged sites, the corresponding electromagnetic fields must take the form: where as before Of course, only the sites k ⊂ i that have an electromagnetic field coupling have non-vanishing r Q,i , so the multimetric 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 φ i given by Eq. ( 80), the only changes to the W -tensor components are in the part multiplying l µ l ν , 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 Σ (+) i = 0 on every site, and the partially proportional solutions have a combination of the two for different k ⊂ 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.

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 λ → λ 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→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: with dΩ 2 (D−4) the unit round metric on the (D − 4)sphere, given by: 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: again with: With the ansatz Eq. ( 83) for the metric, the tetrads have the following form: r 2(D−3) dv ( 87) r 2(D−3) dv ( 88) . . .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 D−2 i /a 2 i times) the i-th effective cosmological constant Λ i , in the usual manner.The only non-vanishing off-diagonal term is the rv component, which becomes: Remarkably, both the charge and cosmological constant contributions factor nicely into the collection of terms multiplying Σ (±) i .The solutions then split into the usual three branches: the proportional solutions are those where one has r s,i = r s , Λ i = Λ and r Q,i = r Q ∀ i, the non-proportional solutions are those with Σ (+) i = 0 ∀ 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 Σ (+) i = 0 has fixed the ratios of all the adjacent conformal factors.However, because the Λ i no longer have to be the same, the β (i,i+1) m are now unconstrained.
For the partially proportional case, some, but not all, of the adjacent sites do still need to have the same Λ, which, as previously, makes things more fiddly.The situation is as follows: if there are n total sites k ⊂ i which have Σ (+) 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 Λ i , as the sites that do not have Σ (+) k = 0 must have Λ k = Λ k+1 .Then, from the N total diagonal equations, n of them fix n of the independent Λ 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 Λ.Again, the β (i,i+1) m are unconstrained here, as the n equations that constrained them in the rotating case now act instead to fix the initially free Λ 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.

C. 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-GRadjacent 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 θ 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 numerically 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 4 Actually, 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).
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 Schwarzschild(-dS) solution i.e.where both metrics are proportionally Schwarzschild(-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 A 1, the proportional Schwarzschild 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 Schwarzschild 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, Σ (+) 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 [112][113][114] (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,[115][116][117][118][119][120].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 Σ (+) 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] 5 .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,121] for bigravity (it proves simpler to work in the metric formalism here).The metrics are expanded as: where the normalisation is to ensure the kinetic terms are canonical.The resulting equations of motion for the perturbations read: where the Lichnerowicz operator is given by [53]: and all quantities have been rescaled by the appropriate powers of a i , as in section III A, so that all tensors in Eq. ( 95) behave as if they live in the common background of ḡµν .For chain type interactions, the mass matrix, M 2 , for the rescaled perturbations in this common background, is tridiagonal, with the following non-vanishing components: As per the discussion in Section II, the mass eigenstates are linear combinations of the h ) via some N × N orthogonal matrix O whose columns are the mass eigenvectors; that is, 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 ∝ a i [12,14,39].Generally, it is not possible (save for an exceptional case that we will consider shortly) to determine the higher mass eigenvalues/eigenvectors analytically [122,123], 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: Whenever the effective cosmological constant is nonvanishing, it is crucial that the masses satisfy m 2 ≥ 2 Λ/(D−1).This is the well-known Higuchi bound, below which the helicity-0 graviton modes become ghost-like [124][125][126][127].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 [128][129][130][131][132][133].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.[134]).However, there is strong evidence that the PM gauge invariance does not survive to the full non-linear level within the dRGT framework [135][136][137][138], 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 [139][140][141]).
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 de- are present in the mass matrix, and by Eq. ( 73) these are related as Σ 0 .Therefore, one has explicitly that: (102) which one can check has eigenvalues: with the corresponding mass modes being: defining γ = M 1 /M 0 as the ratio of the gravitational 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,121], as they should (though in [121] the perturbations were defined differently to in our Eq.( 94), as µν , without the extra factors of a i and M i , so the precise structure of their µν 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.

Clockwork gravity in 4d
We defined a clockwork theory back in section III A 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 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 β With these choices, the mass matrix has components: 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 ∼ qM , with the heavier modes distributed exponentially above this.The massless eigenvector is given by O j0 = N /q j , where the normalisation is 6 Actually, 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.
If matter couples only to the metric g (N −1) µν at the end of the chain of interactions i.e.only T (N −1) µν is nonvanishing (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: where the effective Planck scale is: 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.

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 A 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 ḡµν 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 β [13,14]: where δ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) δy.With these parameters, one recovers 5D GR from the 4D multi-metric theory upon taking the limit where δy → 0 and N → ∞ while keeping the product N δ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 Λ = 0 and a i = 1 ∀ i.
The form of the mass matrix here reflects the simplicity of the vacuum, reading: which one can diagonalise using the techniques developed in [122].The result is that the mass eigenvalues are given by: the massless mode is simply: while the massive modes are: with normalisations 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 ∼ n/L and the heaviest as m N −1 ∼ 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 2 Pl = N M 2 (4) = N M 3 (5) δ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 ∇ν H (i) µν = ∇µ H (i) , and so ∇ν ∇µ H (i) µν = ¯ H (i) .Using the latter in the trace of Eqs.
(95), one is forced into de Donder gauge for all the mass modes simultaneously: ∇ν With this restriction, the linearised equations for the mass modes become: whose behaviour can then be analysed depending on one's choice for the background metric ḡµν .

A. Multi-Schwarzschild
Eqs. (117), with the common background ḡµν taken as the D = 4 Schwarzschild metric, are precisely those equations studied in the context of the GL instability [53][54][55][56][57]61].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: and were shown to be unstable (i.e. have Ω > 0) within the range: 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-Schwarzschild 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-Schwarzschild-dS, and determined that the instability turns on when mr s 0.86 (which, interestingly, is precisely the threshold at which the hairy solutions discussed briefly in section III C 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 M 2 ; all of the mass modes are subject to the same stability condition, where the proportional Schwarzschild solution becomes unstable if, for any of the m i , m i r s 0.86.Consequently, the solution is stable if and only if the lightest massive mode H (1) µν 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 [142].Typically, the masses have been taken of order the Hubble parameter today (m ∼ H 0 ∼ 10 −33 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 [143][144][145][146][147].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 GR 7 .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 [148,149].
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 ∼ H 0 then any Schwarzschild black hole weighing less than 10 22 M ⊙ suffers from the instability [63].However, for such light graviton masses, Ω scales linearly with m [150][151][152], so the characteristic timescale of the instability Ω −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 [8,9] and is also sensible in clockwork scenarios [12,14]), only black holes weighing less than roughly 10 −22 M ⊙ 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.

B. Multi-Kerr
If one now allows the black holes to rotate, and so instead takes the common background ḡµν 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 Schwarzschild case (m i r s ∼ 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 Schwarzschild case, as a consequence of the superradiant instability of rotating black holes against massive bosonic excitations [153][154][155].This effect occurs when the frequency of the perturbation satisfies: where m A is the mode's azimuthal quantum number and Ω 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 ω = m A Ω BH .The condensate dissipates via (almost monchromatic) quadrupolar gravitational wave emission [156].
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 ∼ 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 ∼ 10 −11 (M ⊙ /M 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 ≪ 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 growingfor 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 [157] 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 [158,159], which has a well-posed dynamical formulation [160][161][162], 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 solutionsthe 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 multimetric 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 ∼ O(1) (differing slightly between the Schwarzschild 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.

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 [164].For an equivalent derivation in the metric formalism, we refer the reader to e.g.[121].
The variation of the kinetic term with respect to the i-th tetrad gives: where ⋆G a = − 1 2 R bc ∧ ⋆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): from which we can use the fact that e k ∧ ⋆G a = G a k ⋆1 to extract the components.First, we have: We see that the Einstein tensor has (tetrad basis) components given by: or in more familiar form: Now, for the potential term, the variation is: in which we define the W -tensor 1-form in an analogous way to our Einstein 1-forms above.Explicitly, To extract the components, we use the same trick as before -namely, that dx µ ∧ ⋆ a ⋆ (i) 1, only this time wedging with the coordinate basis 1-form dx µ , since the tetrads in Eq. (A9) 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 ⋆ (i) 1; we make the natural choice of the coordinate basis).The result is that: Contracting with e (i)a ν gives us the appropriate coordi-nate basis expression, Eq. ( 16).
This time, the calculation is simpler in the metric formalism, and follows closely the bigravity derivation in [121].First, recall that the metric form of the W -tensor for chain type interactions reads (C.F.Eq. ( 64)): 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 ḡµν ) take the simple form S i→i±1 = (a i±1 /a i )½, which means that the eigenvalues of S i→i±1 are simply D copies of (a i±1 /a i ).Therefore, the elementary symmetric polynomials are: which recovers Eq. ( 26) upon substitution back into Eq.(B1).To linearise the system, we perturb around our proportional background.The metric and its inverse are expanded as: µν = a 2 i ḡµν + δg (i) µν (B4) where we have included factors of a i so that δg (i) µν behaves as if it lived in the common background described by metric ḡµν (if δg (i) µν had its indices instead manipulated with a 2 i ḡµν , the inverse would simply be g (i)µν = a −2 i ḡµν − δg (i)µν ).It is well known (see e.g.[34]) that the Einstein tensor linearises to the Lichnerowicz operator acting on the perturbation, that is: αβ , (B6) 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: Since it is built from the metric variations, δS µν also has its indices manipulated with the common background metric ḡµν .Eqs. (B7)-(B9) hold in general, but around the proportional background things simplify greatly: as we saw, the building block matrices take the simple form S i→i±1 = (a i±1 /a i )½.After substituting this in above (and using some properties of the binomial coefficients) one can show that the Y variations become [121]: and so the variation of the W -tensor is: (±) i 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 δS i→i±1 .By starting with (S 2 i→i±1 ) µ ν = g (i)µλ g (i±1) λν and substituting in Eqs.(B4) and (B5) for the perturbed metrics, one can show that the desired variation is given by: (δS i→i±1 ) µν = 1 2a 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 [121].
If one instead parametrises their perturbations as 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 µν of course has different components to the mass matrix for δg µν , 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 δg (i) µν or h (i) µν to express their perturbations.Finally, we note that if ḡµν is the D-dimensional (A)dS metric, then one may additionally pull out the cosmological constant from the Rα β µ ν background curvature contribution to Ēαβ µν , to express the first two terms in Eq. (B13) as: where Ẽαβ now contains only the covariant derivative operators, if one wishes to match up exactly with the corresponding equations as written in [39,121] 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 ḡµν is either the Schwarzschild or Kerr metric.
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 β (0,j) m 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: 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: D−m , respectively.The equations then reduce as before to N copies of Einstein's equations for the ḡµν common background, which one may in principle solve for Λ 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 then the star type analogue of (say) Eq. ( 65) for the chain interactions is: and similarly for the charged variants of the above (C.F.section III B 3).In principle, solutions then exist in each of the proportional (all r s,j = r s,0 ), non-proportional (all Σ (+) j = 0) and partially proportional (combinations of both) branches.
As for the linearised equations, the star type analogue of Eq. (B13) for the metric perturbations is: One may check, for example, that in D = 4 spacetime dimensions, and denoting a 0 = 1, a j = c j , γ j = M j /M 0 , these equations imply a mass matrix for the δg This is precisely the star type interaction mass matrix that was derived for the D = 4 multi-metric theory in [39].

2 i
times the sum involving a i−1 , then what we actually have is the following: µν .To obtain them, we rotate to the field basis H µν = (H (0) µν , . . ., H (N −1) µν ) in which the mass matrix is diagonalised, related to the original basis h µν = (h (0) µν , . . ., h (N −1) µν 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 [163] (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.