Solid Holography and Massive Gravity

Momentum dissipation is an important ingredient in condensed matter physics that requires a translation breaking sector. In the bottom-up gauge/gravity duality, this implies that the gravity dual is massive. We start here a systematic analysis of holographic massive gravity (HMG) theories, which admit field theory dual interpretations and which, therefore, might store interesting condensed matter applications. We show that there are many phases of HMG that are fully consistent effective field theories and which have been left overlooked in the literature. The most important distinction between the different HMG phases is that they can be clearly separated into solids and fluids. This can be done both at the level of the unbroken spacetime symmetries as well as concerning the elastic properties of the dual materials. We extract the modulus of rigidity of the solid HMG black brane solutions and show how it relates to the graviton mass term. We also consider the implications of the different HMGs on the electric response. We show that the types of response that can be consistently described within this framework is much wider than what is captured by the narrow class of models mostly considered so far.


Contents
Recently, a somewhat surprising connection between Condensed Matter (CM) and Massive Gravity (MG) has been unveiled: MG theories are especially relevant for CM applications by way of the AdS/CFT correspondence [1][2][3][4]. This is quite appealing because CM contains a large variety of poorly understood systems, whereas both MG and AdS/CFT have enjoyed a vast development in the last decade. The fact that CM and MG must be connected somehow is actually very natural: any distribution of matter gives rise to a plasma mass for the graviton, much in the same way the photon gets a mass in a medium with a finite density of freely moving charges. Since a plasma mass is linked to the presence of a material, it is inevitably non-Lorentz invariant and so it comes as no surprise that it can be connected with Lorentz Violating Massive Gravity (LVMG). It is well known now that there are many types or phases of LVMG that are free from the consistency problems which plague the Lorentz invariant versions of massive gravity [5,6]. Moreover, the LVMGs are believed to be better behaved also in the UV region in comparison to the LI massive gravity theories. 1 Thus, it seems to be possible to match the great diversity of types of condensed matter materials with the great diversity of consistent phases of LVMG.
This connection, surely quite obvious to most MG experts, has always been clear but seemed of little consequence because the graviton plasma mass inside a material with an energy density ρ can be estimated as m 2 g ∼ ρ M −2 Pl , which is much smaller than the inverse size of human-scale materials. Cosmology is the area where one can have, in principle, more significant effects, even if to date it is unclear whether any form of MG may help with the open problems in cosmology (for a recent review see, e.g. [8]).
In AdS/CFT correspondence [9], though, this connection opens an entirely new dimension. MGs in AdS-like spaces are dual to strongly coupled materials that incorporate a crucial aspect of CM: a sector that breaks spatial translation invariance. From the CM perspective, including this sector in the low-energy description is of central importance to capture e.g. the effects from phonons and from disorder (see below).
In our view, however, the way how MG plays this central role has not been clearly stated yet in the AdS/CMT literature. To make contact with the AdS/CMT applications, we shall stick to the common terminology that refers to the MG theories in AdS (or more general spaces that allow a holographic interpretation) as Holographic Massive Gravity (HMG). The main aim of this paper is to clarify the role, consistency, and classification of HMG theories. The key messages will be: 1. HMGs can be broadly and clearly separated into solids and fluids 2 . This follows both from the analysis of the residual diffeomorphisms that are preserved as well as from the presence/absence of an elastic shear response.
2. The large family of HMGs considered below are very healthy: they are fully consistent as Effective Field Theories (EFTs) in that they are free of ghosts and other pathologies at linear and nonlinear level, even if they are not necessarily of dRGT type [10]. In addition, and in contrast to the application to cosmology, in holography the graviton mass does not need to be so much smaller than the Planck mass: it suffices that it is around the AdS curvature scale 1/L wich can be taken to be, say, one to two orders of magnitude below M Pl . In this case, the strong coupling scale (which is usually a geometric mean of the form Λ n ∼ (m n−1 g M Pl ) 1/n for n = 2, 3, 5) associated with the graviton mass sector does not differ significantly from the Planck scale, and in any case is not below 1/L. This is quite important for the holographic application because it allows the gravitational description to continue to be weakly coupled (at least at low frequencies and long wavelengths) even in the presence of a fully nonlinear mass sector. 3. The class of HMGs relevant for CM is much more general than the dRGT models [10] that have been considered mostly so far (see e.g. [1,2,4]), as was already argued in [11]. The family of theories that are free from various pathologies (as discussed e.g. in [11,12]) and have the same symmetries is parameterized by free functions that are related to the physical properties of the materials. This has some practical implications because there are many properties of the dRGT-like massive gravities that are not generic in the AdS/CMT context, and therefore can lead to misguided implications for CM. For instance, the full set of HMG theories that can be applied to modelling the low energy limit of CM systems is not only parameterized by two real numbers.
All in all, this suggests that holographic massive gravities might be very relevant for condensed matter applications.
In this paper we shall use both the presentation of massive gravity in terms of broken diffeomorphisms and in its covariantized form relying on the Stückelberg fields φ A . In the latter language MG can be seen as a theory of general relativity coupled to a number of scalar fields. The application of the scalar fields formalism to holographic massive gravity was initiated in [13,14] and was instrumental in identifying the additional degrees of freedom as phonons [11].
The most meaningful way to characterize the various different phases of massive gravity is phrased in terms of the symmetries that are left unbroken in each phase. The situation is very similar to the organization of the EFTs for various materials according to the spontaneous symmetry breaking pattern of the time and/or space translations, [15][16][17][18]. The classification of the phases of MG proceeds in very similar terms because on physical grounds the LV mass terms are the possible forms of the plasma-mass that is generated in different types of materials. In other words, the EFTs of solids/fluids/etc. are related to the phases of LVMG by gauging the spacetime symmetries, that is, by introducing the coupling to the dynamical metric by the usual covariantization prescription. Thus it is not at all surprising that one can speak of, e.g., solid and fluid phases of MG. For the sake of simplicity, we shall restrict ourselves to these two cases and defer for future work a more thorough analysis of other phases of MG.
As is nowadays relatively well-understood, the EFT for fluids and solids in flat space involves a set of phonon scalars φ I (in 2+1 dimensions, I = 1, 2) that enjoy internal shift and rotation symmetries for homogeneous and isotropic materials [15][16][17][18]. The internal symmetry group for solids is the two-dimensional Euclidean group of translations and rotations. For fluids, the internal group is much bigger and includes also volume preserving diffeomorphisms (VPDiffs). The scalars acquire an expectation value φ I = δ I i x i and break the product of the (space transformations) ⊗ (internal transformations) to the diagonal subgroup. For fluids, the preserved symmetry includes a volume preserving diagonal subgroup.
The effective Lagrangian at the lowest order in derivatives in the two cases can be written as [18][19][20][21] L (solids) = V s (X, Z) and where X = tr I IJ and Z = det I IJ with I IJ ≡ η µν ∂ µ φ I ∂ ν φ J . 3 The functions V s,f encode the linear and nonlinear properties of the solid and fluid, and they are free functions subject to mild consistency constraints. It is easy to realize that gauging these theories leads to graviton mass terms around the solution with φ I = δ I i x i . The simplest way to see this is to replace η µν with g µν and to go to the unitary gauge where the scalar fields are fixed to be equal to their background configuration. The above solid/fluid Lagrangians then become nonlinear potential terms for the metric where g ij denotes the spatial part of the inverse metric. We note that exactly the same procedure was followed in the so-called 'solid inflation', ref. [21]. At this point, it is quite clear that there must exist a fully equivalent nonlinear formulation of solid/fluid MGs that is phrased entirely in terms of a unitary-gauge metric variable. The form of this action is dictated by requiring it to be invariant under certain subset of the diffeomorphisms, that do not include the spatial diffeomorphisms x i →x i (t, x j ). The preserved diffeomorphisms are the ones enjoyed by the potential terms in (1.2). Both for solid and fluid MGs, these include the time-reparametrizations t → f (t) plus global translations and rotations that force the potential not to depend explicitly on x i and to contract the spatial indices with Kronecker delta δ ij . 4 For fluid MG the potential is also invariant under the spatial VPDiffs, that forces it to be a function of det g ij only. Importantly, as we shall see below, the VPDiff symmetry protects the vanishing of the physical mass parameter of the metric tensor modes.
In order to link the defining symmetries of solid/fluid MGs to the structure of mass parameters, it is important to identify the relevant notion of the mass terms, which is not entirely obvious since we need to work in curved backgrounds. The most practical definition is to expand the action around the solutions of interest and look at the nonderivative quadratic terms once the kinetic terms have a canonical form. Limiting the discussion here to homogeneous and isotropic backgrounds, one can follow [5,6] and parameterize the possible mass terms by five constant mass parameters, 5 For fluids/solids, the preserved symmetries at the level of the Lagrangian in the unitary gauge formulation, and the allowed mass terms are as follows: Solids Preserved symmetries: time reparametrizations and the diagonal subgroup of space translations and rotations. Expressed as infinitesimal diffeomorphisms ( The ⊕ label denotes that it is a combination of a spatial translation/rotation and a corresponding internal transformation that leaves the background configuration invariant.The allowed mass terms are m 1,2,3 = 0 Fluids Preserved symmetries: time reparametrizations and the diagonal subgroup of translations, rotations, and volume-preserving diffeomorphisms. In infinitesimal form, these look like ξ t (t), ξ i V P ⊕ (x j ) (VP = volume preserving, ∂ i ξ i = 0). The allowed mass terms are m 1,3 = 0. Importantly, for fluids m 2 = 0. 3 In higher dimensions there are more invariants. For instance, in 3+1 dimensions tr I I J I J K gives an independent invariant. 4 For simplicity, we shall assume that he kinetic part of the action is given by the standard Einstein-Hilbert term, which is not the most general one compatible with these symmetries. 5 The precise link between these mass terms and the nonlinear Lagrangian (1.2) will be done in full detail in sections 2 and 5 in a class of models that covers fluids and solids.
Let us insist that in both cases the spatial translations are broken (or non-linearly realized), which we emphasize with the ⊕ label that should remind that this is not a standard translation. In the realization of these phases as HMGs presented in the main body of the paper, both cases lead to a finite DC conductivity in the electric response. Also, both the fluid/solid EFTs and the fluid/solid MGs are very healthy in the sense that the Lagrangian can be always chosen in order to avoid ghosts or other pathologies at linear and nonlinear levels. We show that explicitly in Section 5.
The main differences between the two types of theories are i) that the solid phases exhibit propagating transverse phonons -the Goldstone modes of the broken space translations, inhomogenenous in spatial coordinates -whereas the fluids do not; and ii) that the tensor modes are massive/massless for solid/fluid phases respectively. In the HMG constructions below, both statements will hold at the level of local propagating fields in AdS. However, since local degrees of freedom in the dual field theory (such as the transverse phonons) can be over-damped in black brane backgrounds, the difference in this respect between the solid/fluid HMGs becomes slightly blurred. Instead, a clear distinction survives at the level of the tensor modes being massive/massless. In turn, as we shall show below, this allows to characterize the solid types of HMGs as the ones with rigidity, a nonzero static elastic shear response, encoded in the non-vanishing m 2 .
After this quite generic review of the physical distinctions between MG phases, let us introduce a bit more the specifics of the holographic application (for some recent developments concerning the importance of the translation breaking sectors see [1][2][3][4][22][23][24][25][26][27] and references therein). As already mentioned, we will focus on asymptotically AdS 4 backgrounds, since they admit simple and well-understood 2+1 CFT interpretations. We shall denote by t and x i the coordinates along the boundary directions and r the holographic coordinate dual to the renormalization group scale. Since the holographic map already gives r a special role, we consider gravity theories with an anisotropic mass term. The crucial ingredient to make contact with realistic CM is to break translations in the x i directions. Thus, in the non-unitary gauge description we are going to need at least two scalars φ I for I = 1, 2 with vacuum expectation values φ I = δ I i x i . 6 One can also consider additional radial and temporal fields φ r and φ t . The role of a φ t ∝ t scalar is to introduce super-fluid/solid versions of the dual material, and the role of φ r is basically the same as of the commonly-used dilaton field. This identification especially holds when φ r does not enjoy an internal translational symmetry, and thus provides a rich variety of radial dependence in black brane as well as in vacuum solutions (zero temperature and density). For the sake of simplicity, we shall restrict here mostly to the simplest case that provides momentum dissipation in the x i directions -a two-field model with the φ I scalars only.
These two-field HMG models are already quite rich (as they include solid and fluid behaviours), simple and interesting. Because of their dual field theory interpretation, it seems appropriate to call them solid/fluid CFTs. Indeed, the asymptotically AdS black brane solutions in the solid/fluid MG admit an interpretation as CFTs with solid/fluid behaviour that are able to dissipate momentum, and which are deformed only by a finite temperature and finite density. The fact that momentum is not conserved can be seen as a coupling of a standard momentum-conserving CFT to a translation-breaking sector. The latter is realized in the black brane solutions via the φ I -hair. The magnitude of the translation-breaking is proportional to the local energy density in the φ I sector and is position-(i.e., scale-) dependent. Importantly, depending on the form of the potential V (X, Z) this energy density can be located closer or farther from the black brane horizon. In the CFT, this characterizes the typical energy scale of the translation breaking. At least for the solids it seems suggestive and consistent to identify the explicit breaking (that is localized near the AdS boundary) as some type of lattice disorder in the sense that it mimics a random distribution of lattice defects in an otherwise perfect crystal. In any case, the presence of defects seems unrelated to the fact that momentum dissipates, as is clear a posteriori since the DC conductivities involve horizon data only. We shall not deepen much more into the nature of this disorder in this work (for the role that holographic version of disorder can play in various contexts, see [28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44] and references therein). Still, having identified a clear distinction between solid and fluid materials might help identifying the nature of disorder, simply because the notion of disorder itself seems much better defined for solids than fluids in the first place.
The focus of our work is to analyze the full set of mass terms that are available in HMG and to clarify how the various mass terms control different transport properties of the dual materials. Obviously, there are many more mass terms if one restricts to homogeneity and isotropy of the CFT in the spatial directions, x i , but not in the holographic direction. Hence, the full classification of the phases of HMG is rather complicated. We shall therefore initiate the analysis of the effects from the full range of possible mass terms, but will also restrict to the cases that are simpler to interpret. Along the way, we will perform various checks such as that the DC conductivity is indeed finite in both solids and fluids, and that the solids, in addition, enjoy a nontrivial elastic response.
The rest of this paper is organized as follows. In section 2 we set up our notations and introduce a classification of the massive gravity theories according to the preserved subgroup of the diffeomorphism symmetry. We describe two particular cases, on which we focus in the present paper: the general theory with negligible background stress energy tensor and the theory written in terms of two Stückelberg fields -the simplest theory describing the solid phase of HMG. Section 3 is devoted to the electric response in the dual field theory. We derive the expression for the DC conductivity in terms of the mass parameters and discuss the phenomenology of the representative models. We proceed in section 4 with discussing the elastic properties of the duals to the solid and fluid phases of HMG. Section 5 contains the analysis of the massive gravity with two scalar fields. Section 6 is devoted to the conclusions and outlook. A detailed analysis of the HMGs with four scalar fields and zero stress-energy tensor is presented in appendices A and B.

Phases of Holographic Massive Gravity
For holographic applications in condensed matter theory we are interested in massive gravity theories that allow for asymptotically AdS charged black brane solutions. The action that will be considered in this paper is the Einstein-Maxwell action with a negative cosmological constant and a graviton mass term: 7 The graviton mass term L φ can be written as a Lagrangian for the Stückelberg scalar fields φ A and will be specified in the following sections. We shall concentrate on the probe limit and neglect the backreaction from the mass term on the background metric everywhere apart from section 5. In this limit the action admits an asymptotically AdS Reissner-Nordström black brane solution: with the emblackening factor f (r) given by Here the mass M is a constant of integration that can be determined by demanding that f (r) vanishes on the horizon r = r h to be The solution for the Stückelberg scalars and the Maxwell field takes the form In the following sections we shall discuss the different phenomenological consequences in the dual theory induced by the presence of the graviton mass in the bulk and investigate the bulk stability of the perturbations around the black brane background. To do so we consider the perturbations of the metric and Maxwell fields defined as 6) and the perturbations of the Stückelberg scalar fields π A defined as Since we are interested in homogenous and isotropic condensed matter systems, we shall limit ourselves to the mass terms that preserve the constant shifts and rotations of the transverse coordinates x i = {x, y}. It allows us to classify the perturbations according to the scalar, vector, and tensor representations of the transverse O(2) rotation group. The split depends on whether we consider homogeneous modes, i.e. independent on the transverse coordinates x i , or inhomogeneous modes. For simplicity we shall consider only the homogeneous case when the perturbations can be classified as (2.10)

General mass terms
A straightforward way of analyzing the stability and phenomenology of a general theory of massive gravity is to start with the non-covariant form of the massive gravity action with the mass term written it terms of the metric perturbation h µν and then restore the diffeomorphism invariance by the Stückelberg trick. A similar analysis of the stability of the Lorentz violating massive gravity around Minkowski background has been performed previously in refs. [5,6].
We define the inverse metric perturbations 8 as and write the most general quadratic mass term that preserves the rotations of the transverse coordinates in the following form: where all the masses m 2 i are functions of the radial coordinate. The repeated transverse coordinate indices are contracted with δ ij . The numbering of the mass parameters is chosen to be consistent with the notations of refs. [5,6]. The novelty introduced by the planar AdS black brane background is the special role of the holographic coordinate r. For future reference we note that the mass parameters defined in (2.12) can be classified with respect to the perturbations that they affect as: We note that this classification is valid only for the homogeneous modes.

Symmetries
The mass term (2.12) explicitly breaks the spacetime diffeomorphisms since the metric perturbation h µν is not invariant under the coordinate transformations x µ →x µ (x ν ). In fact, as we have discussed in the introduction, this can be considered as the defining property of massive gravity -a generic massive gravity theory is a theory that breaks some subset of the diffeomorphism invariance of the Einstein-Hilbert gravity. In order to reveal the symmetry breaking pattern of a given massive gravity theory it is useful to restore the diffeomorphism invariance by introducing four Stückelberg scalar fields φ A , which represent the physical coordinates. 9 One can then write the mass term using two ingredients: a gauge invariant version of the space-time metric and a reference metric in the configuration space of the scalar fields which can be a function of the Stückelberg fields themselves. The metric f AB is used to raise and lower the scalar fields space indices. 10 Any mass term written in terms of I AB and f AB will be manifestly invariant under the general coordinate transformations. In addition, depending on the exact form of the scalar fields Lagrangian, it can be invariant under certain internal symmetries of the scalar fields φ A → ψ A (φ B ). The field configuration spontaneously breaks the internal symmetries and the spacetime diffeomorphisms to a diagonal subgroup -it is left invariant only by a combination of simultaneous internal field redefinitions and spacetime diffeomorphisms. For example, under the transformations This is the diagonal subgroup of the internal and spacetime symmetries.
In general, the scalar fields Lagrangian does not admit the full reparameterisation invariance of the scalar fields due to the fact that the reference metric f AB is fixed and, hence, non-dynamic. Thus, the scalar fields background spontaneously breaks the spacetime diffeomorphisms to the product of the residual internal symmetries and their spacetime counterparts. For example, if the reference metric does not depend on the scalar fields, the Lagrangian is invariant under constant shifts φ A → φ A + c A and the diagonal subgroup contains the spacetime translations. An important case is when all the internal scalar field indices are contracted in the Lagrangian. Such theories are invariant under the reparameterisations of the scalar fields that are isometries of the background metric f AB . This is the case in the Lorentz invariant massive gravity theories where the Minkowski reference metric leads to the invariance under internal Poincaré transformations. In the context of effective field theories of solids and fluids (see [17,18] and references therein), the Lagrangian is invariant under rigid rotations and volume preserving diffeomorphisms of the scalar fields, respectively [19]. The background configuration (2.18) breaks these symmetries down to the diagonal subgroup of internal and spacetime translations and rotations, or volume preserving diffeomorphisms. In the present work we are making the first steps towards using the effective theories of solids and fluids in the holographic context and applying these theories to the condensed matter systems.
Leaving some of the scalar field indices not contracted (e.g. when the Lagrangian may contain a sole I AB as opposed to I AB f AB ) breaks the internal scalar field metric isometries to some residual subgroup. In general, this leads to Lorentz violating massive gravity theories which are of the main interest of this paper. In particular, we choose to preserve the internal O(2) rotations of the transverse scalar fields φ I . This allows us to use I AB as the building block of our effective Lagrangian with the only condition that the transverse indices have to be contracted with a reference metric proportional to δ IJ . As described above, the full spacetime diffeomorphisms are preserved by the action and are broken spontaneously on the background configuration (2.18).

Probe limit Lagrangian with T µν = 0
A straightforward but not unique way of constructing a gauge invariant form of the most generic quadratic mass term (2.12) is to set the scalar fields reference metric to be equal to the background black brane metricĝ µν : We can then introduce diffeomorphism invariant metric perturbations as In the unitary gauge, when the scalar fields are set to be equal to their background solution (2.18), the field H AB equals to the inverse metric perturbations: The covariant form of the action (2.12) is then obtained by replacing the different components of the metric perturbations h µν with the gauge invariant fields H AB and the explicit r-dependence of the masses with an explicit dependence on the radial Stückelberg field φ r . The resulting covariant mass term L φ (H AB , φ r ) gives a stress-energy tensor that vanishes on the background and, thus, does not contribute to the background equations for the metric. Although the exact vanishing of the backreaction on the metric is a very particular case it can be considered as a probe limit approximation for a generic massive gravity theory. In particular, the number of propagating degrees of freedom as well as the stability requirements of the theory should not be altered in this limit. The exact form of the covariant mass term is given in eqn. (A.1) in appendix A, where also the linear stability analysis of the theory is performed. The phenomenological consequences of the mass term (2.12) relevant for the electric response of the dual field theory will be discussed in section 3. We note here the obvious fact that the number of scalar fields used in our Lagrangian is a matter of choice. However, it limits the available components of H AB and, consequently, the subset of mass terms out of the ones presented in (2.12) that arise if not all four scalar fields φ A with A = t, r, x, y are used. Of particular interest for us is the case of only two scalar fields, i.e. φ I with I = x, y, that is commonly used in the context of holographic massive gravity and will be discussed in great detail in section 5. In this case, only the H IJ components can be constructed which implies that all the quadratic mass terms in Moreover, since the black brane reference metric f AB (φ r ) explicitly depends on the absent field φ r then even H IJ can only be present in a particular combination that includes its traceless part only, i.e. (H IJ ) 2 − 1 2 H 2 . This leads to a quadratic mass term of type (2.12) with The same conclusion about the available mass terms could have been reached by noting that the absence of fields φ t , φ r in the mass term building blocks I AB and f AB leads to a residual symmetry that is left unbroken by the scalar fields background configuration φ I = x i . These are the spacetime diffeomorphisms x a →x a (x µ ) with a = t, r. Hence, all mass terms in (2.12) that transform non-trivially under these diffeomorphisms are forbidden.
As a result, the above conditions prohibit the appearance of any components of the metric perturbations that are scalars and vectors under the transverse rotations in quadratic action. From the phenomenological point of view it is a case of no particular interest. We note, however, that such strict constraints only arise if we consider generally covariant mass terms with zero background stress-energy tensor, i.e. with zero backreaction on the metric. In practice, the mass Lagrangian L φ (H AB , φ r ) is not the only diffeomorphism invariant mass term that reduces to the quadratic action (2.12) once expanded up to the second order in perturbation theory. The most general covariant mass terms can be obtained by writing the action in terms of the matrix I AB and arbitrary functions of the scalar field φ r . The quadratic perturbative mass term obtained by these other actions will also be of the form (2.12). However, a general mass Lagrangian of this sort will give rise to a non-zero stress-energy tensor that will backreact on the spacetime metric. As a result, some of the quadratic non-derivative terms arising in the perturbative expansion of L φ will vanish on the background equation of motion once combined with the Einstein-Maxwell part in the full action (2.1). Hence, there is no one-to-one correspondence of the quadratic non-derivative term of metric perturbations in L φ and the actual mass terms of the metric perturbations in the case of generally covariant mass potentials with non-zero backreaction.

Two fields case: solids and fluids
A particular example of a massive gravity action with non-zero stress-energy tensor we would like to consider is a generic mass term that depends on only the two transverse Stückelberg fields φ I . Such model is a direct generalisation of the two fields dRGT theory, and provides the simplest way to model momentum dissipation in the dual theory. Since the action does not depend on the φ a scalar fields, the background solution φ I = x i δ I i leaves the {t, r} spacetime diffeomorphisms x a →x a (x µ ) unbroken in this case. As before, the most general graviton mass Lagrangian can be written in terms of the two building blocks I IJ ≡ g µν ∂ µ φ I ∂ ν φ J and some reference metric f IJ . Since we wish to preserve the internal rotations in the (I, J) plane we set the reference metric to be the identity matrix f IJ = δ IJ . Hence, the mass Lagrangian is an arbitrary function of the powers of I with all indices contracted with the identity matrix δ IJ . For a 2 × 2 matrix there are only two algebraically independent contractions. We choose them to be Since any other contraction of I is a function of X and Z, the most general two fields mass term takes the form: The dRGT theory with two fields considered in ref. [1] is a particular case of the two fields massive gravity with the Lagrangian given by The Lagrangian (2.27) of the scalar fields is similar to the effective field theories used to describe perfect fluids [19] and solids [20,21]. In this context the scalar fields φ I denote the comoving Lagrangian coordinates of the material. In order to describe perfect fluids the scalar fields action has to be invariant under field space volume preserving diffeomorphisms For solids, the action is only invariant under constant shifts and rotations We see that, in general, the scalar fields Lagrangian (2.27) describes solids. In the special case when it depends only on the function Z, its symmetry is enhanced, and it describes perfect fluids. It is the presence of the field X that "turns" the material into a solid. In this sense, the two fields Lorentz violating massive gravity is an effective theory that describes gravity in the presence of a solid or a fluid. We shall present a detailed analysis of the two fields massive gravity in section 5.

Electric response
In this section we investigate the phenomenological consequences of allowing for the broad class of the graviton mass terms as in (2.12). We shall concentrate on the dynamics of the vector modes in particular and shall discuss the implications for momentum relaxation and electric response in the dual CFT. Following the discussion of section 2.3 we shall consider the diffeomorphism invariant mass term relevant for vector perturbations For simplicity we assume that the φ r dependence of the masses is given by is the emblackening factor written in a gauge invariant form, and M 2 (φ r ) is a universal mass function of dimension mass squared that is regular and non-vanishing at the horizon when φ r = r h . The masses m 2 i and the powers α i are dimensionless constants. A detailed stability analysis of the quadratic Lagrangian in appendices A and B leads to conditions α 8 = 0 and α 1 = −α 7 = 1. We also find in appendix B that for arbitrary mass parameters m i the vector modes of the metric propagate on an effective acoustic background metric that is different from the Reissner-Nordström background of the Maxwell field. The conditions on the mass parameters for the two effective light cones to coincide read A priori these conditions do not need to be imposed as long as the acoustic metric describes a causally stable spacetime [46]. In appendix B we find no indications that this would not be the case here. Importantly, the mass condition (3.3) is automatically satisfied in the two fields massive gravity described in 2.4. This means that in these theories both dynamical vector modes propagate on the same effective metric. Since the two fields dRGT theory is a subclass of these theories same conclusion applies. Nevertheless, in what follows we shall leave the mass parameters m 2 i unconstrained unless otherwise specified. We also note that the mass term (3.1) gives rise to a vanishing stress-energy tensor and as such can be considered as particular. As was already discussed in 2.3 both the stability and phenomenology of mass terms with non-zero backreaction coincide with the results presented below in the probe limit of vanishing graviton mass.
In appendix A we find that the dynamics of the vector sector of our model can be described in terms of two gauge invariant vector fields a i and λ i . The first one is the perturbations of the Maxwell field while the second field, λ i , was introduced as a Lagrange multiplier and can be expressed in terms of the original fields as Henceforth we shall drop the index i. The resulting equations of motion for the Fourier modes then read

DC conductivity
In this subsection we make use of the analytic approach of ref. [2] for calculating the DC conductivity (for an alternative approach, see [29]). The method relies on the existence of a 'massless mode' in the bulk that is a linear combination of the Maxwell field a and the graviton, which in our calculations is represented by the field λ. In zero frequency limit, the existence of this massless mode implies the conservation of a certain quantity Π in the radial direction. In [2] it was shown that the quantity Π encodes the universal behaviour of the DC conductivity and ensures that it remains constant as one moves from the horizon to the boundary.
Here we shall repeat the analysis for the equations of motion (3.6), (3.7) with the mass parameters set to m 2 8 = 0, m 2 7 = −2m 2 1 so that to satisfy the condition (3.3). In this case the equations of motion (3.6) and (3.7) can be written in the form (3.9) The equations (3.8) can be compared to the corresponding equations in [2] for the dRGT massive gravity. We find that the dRGT mass can be expressed in terms of our parameters m 2 1 and M 2 (r) as: In particular, the β 1 term corresponds to the mass function with ν = −3 and the β 2 term is equivalent to ν = −2.
It is straightforward to check that the mass matrix M has a vanishing determinant and, thus, one zero eigenvalue. The eigenvectors α 1 , α 2 corresponding to the zero and non-zero eigenvalues respectively can be taken to be (3.12) The matrix M can then be diagonalized as M = U DU −1 where D = diag(µ 1 = 0, µ 2 ) and U = ( α 1 α 2 ) is a matrix with its column vectors given by the eigenvectors α i . The massless and massive modes, δλ 1 and δλ 2 , are then: The equations of motion (3.8) then take the form (3.14) We emphasize that unless the mass function M 2 (r) = const the matrix U (r) is a function of r. In the zero frequency limit the massless equation can be written in the form of a radial conservation law A few remarks are in order. First, we observe that the conserved quantity Π can be rewritten in terms of the massless and massive modes as This form is very similar to the conserved Π found in [2] for the dRGT massive gravity. Second, we notice that in the unitary gauge the field λ is related to the metric components as (see equations (A.12) in appendix A) relates to Hence, the conserved quantity Π, as defined in (3.15), reduces to This coincides exactly with the radial momentum conjugated to the Maxwell field that can be obtained by varying the on-shell boundary action with respect to the boundary value of the Maxwell field A i(r) [22,24]: The electric conductivity is then defined as the functional derivative The fact that in the zero frequency limit, the canonical momentum with respect to the r-foliation of the Maxwell field, Π, is conserved in the radial direction, was used by Iqbal and Liu in [24] as a motivation to introduce a fictitious membrane DC conductivity for each constant r slice as a response to the massless modē .
We note that this definition of conductivity is equivalent to computing the linear response of the boundary theory to the massless mode δλ 1 . In practice, we are interested in the linear response to the Maxwell field perturbations δA. However, the membrane DC conductivity has nice properties which we would like to exploit in order to extract information about the actual DC conductivity defined as (3.22) In the limit r → 0 the two conductivities are related to each other as where we have used (3.13) to express the fields a, λ in terms of δλ 1 , δλ 2 . For the two DC conductivities to coincide we have to demand that the second term in the above expression vanishes near the boundary, i.e. that In the case of M 2 (r) = const this means imposing the condition λ (0) = 0 on the boundary. For other choices of mass function M 2 (r) one has to similarly find the appropriate condition on the field λ at r = 0. We shall discuss this in more detail in the next section. The claim of [24] is that the membrane DC conductivity does not evolve in the radial direction and can be evaluated at arbitrary r. In particular, it can be evaluated at the horizon, thus, emphasizing the fact that the DC conductivity of the boundary theory can be entirely determined in terms of the horizon quantities. These ideas were used in the context of massive holography by Blake and Tong in [2]. They have shown that also in the case of massive gravity the membrane DC conductivity (3.21) is conserved. Thus, the electric DC conductivity on the boundary can be evaluated as the membrane DC conductivity at the horizon. Near the horizon the fields a and λ are proportional to f (r) −iω/(4πT ) , and the mass function M 2 (r h ) is regular. Hence f δλ 1 = iωδλ 1 while f δλ 2 vanishes. Using the equation (3.16) we find the DC conductivity: With the identification (3.11) this result coincides with the expression for the DC conductivity in the dRGT massive gravity derived in [2].

Phenomenology
A phenomenologically interesting question to investigate is the different types of materials that can be described within the framework of holographic massive gravity and, in particular, their ability to conduct an electric current. The distinction between different classes of materials is well captured by the temperature dependence of their DC electric conductivity. In the models of massive holography proposed in this paper, it is controlled by the radial dependence of the mass function M 2 (r). This determines the dependence of the DC conductivity on the horizon temperature through T = |f (r h )| /(4π). For a mass function of the form the DC conductivity becomes . (3.27) This defines the distiction between a metallic behaviour (dσ/dT < 0) and an insulating behaviour (dσ/dT > 0). The parameter ν determines the nature of the dual CFT as shown in fig. 1. For ν < −2 the behaviour is metallic while for the opposite case, ν > −2, we are in the presence of an insulator. 11 Since the mass function (3.10) in the dRGT theory corresponds to the cases ν = {−2, −3}, the dual materials exhibit a metallic behaviour there. The possibility to mimic metallic and insulating behaviour (and a transition between them) in the context of holographic massive gravity has been already pointed out in [11]. At last, we remark that the special temperature at which all the models with different values of ν have the same electric conductivity is given by (3.28) and corresponds to the temperature of a black brane with the horizon radius r h = L. We further exploit the freedom offered by the generic power ν to investigate the possibility of having a linear T resistivity ρ = 1/σ DC ∝ T −1 . This is a special feature of strange metals (see [25] and references therein) which evade the usual scaling predicted by the Fermi liquid theory (σ DC ∝ T −2 ). Within our class of models we observe a linear scaling in the resistivity at low temperature only for the power ν = −3 as shown in fig. 2. As mentioned above, the ν = −3 mass function coincides with the β 1 term of dRGT massive gravity. The linear scaling of the DC resistivity for this particular case in the low temperature regime 12 has been already observed earlier in [14]. In [26], the same effect has been seen for a more complicated dilatonic model.

AC conductivity
In order to find the electric conductivity in the dual field theory we need to numerically solve the equations (3.6) and (3.7) with appropriate boundary conditions on the black brane horizon r = r h and on the AdS boundary r = 0.
Boundary conditions on r = 0 In order to see what are the appropriate boundary conditions for the field λ near the AdS boundary, we expand the equations (3.6) and (3.7) in the limit r → 0. We use the following ansatz for the radial dependence of the fields: a = r γ , λ = r β , M 2 (r) = r ν (3.29) and find that the leading order expansion of the fields a and λ near the boundary is a = a 0 + r L a 1 + . . . ti (see e.g. [22]). For M 2 (r) = r ν with ν + 1 < 0 the second term in (3.31) diverges and and the λ 0 term becomes subleading with respect to the λ 1 term. The correct near boundary expansion in this case would be For the regularity of the field λ near the AdS boundary we shall thus demand that In section 3.1 it was shown that the DC conductivity that is finite and conserved in the radial direction is due to the response of the theory to the mode that is massless in the high frequency limit and is given by (3.21). In distinction from the Einstein-Maxwell theory the massless mode is a linear combination of both fields, a and λ. This means that if we wish to consider only the response to the Maxwell field then the additional condition (3.24) needs to be imposed. For mass function M 2 (r) = r ν this reads: (3.34) This ensures that the contribution of the field λ to the membrane DC conductivity vanishes. Without this condition, the usual DC conductivity as defined in (3.22) is infinite. In the case when ν + 1 < 0 the expression above automatically vanishes at r = 0. We notice that the two conditions (3.33) and (3.34) complement each other and impose a general condition that λ (0) = 0 for all values of ν.

Boundary conditions on r = r h
In order to determine the asymptotic behaviour of the two fields a and λ in the vicinity of the horizon we use the ansatz where the functions A(r) and Λ(r) can be expanded as A(r) = a 0 + a 1 (r − r h ) + . . . and Λ(r) = l 0 + l 1 (r − r h ) + . . . . Expanding the equations of motion (3.6) and (3.7) leads to There are three cases when these equations can be satisfied. Case 1: The first case occurs when γ = β so that the right hand sides of the above equations become subleading and the equations reduce to The first equation gives γ = ±iω/|f (r h )| and we impose the usual ingoing boundary condition for the Maxwell field a: where T = |f (r h )|/(4π) is the black brane temperature. Given the initial assumption γ = β, the equation (3.39) can only be satisfied when the graviton masses obey the relation This is exactly the condition (B.14), discussed in appendices A.1 and B, ensuring that one of the sides of the light cone of the effective acoustic metric g ab coincides with one side of the gravitational metricĝ ab . For this particular choice of the masses we can thus impose the ingoing boundary condition also for the field λ: Case 2: Another possible choice of the coefficients β and γ is γ − β + 1 = 0. In this case, the power γ is again determined by the equation (3.38) and, hence, we can use the following boundary conditions for the fields a and λ: The additional dependence on f (r) in the above ansatz for λ(r) arises due to the relation β = γ + 1. Effectively, such an ansatz results in demanding that the field λ(r) vanishes on the horizon. We therefore expect that the effect of the additional degree of freedom introduced by massive gravity is negligible in this case and, thus, cannot help in rendering the DC conductivity finite. Nevertheless, it is interesting to perform a more detailed study of this case due to the fact that, unlike in Case 1, there are no constraints on the mass parameters m 2 i . This allows us to probe the parameter space of massive gravity that admits superluminal propagation velocities.
Case 3: The last possibility to satisfy the near horizon equations for integer β and γ is to set β − γ + 1 = 0. Then the boundary conditions for the two fields read However, this amounts to demanding that the Maxwell field vanishes on the horizon which makes this an irrelevant case for the discussion of the electric conductivity and we shall not consider it in the remaining of the paper. Finally, we note that due to the fact that the mass function M 2 (r) is chosen to be regular on the horizon, the particular form of it has no impact on the near horizon boundary conditions.

Results and phenomenology
The AC conductivity can be found by numerically solving equations (3.6) and (3.7) and using the relation (3.20): The numerical results depend on the mass parameters m 2 i and on the functional dependence of the universal mass function M 2 (r). Here we consider two different cases depending on whether the relation (3.41) is satisfied or not. These correspond to the Case 1 and Case 2 in our above discussion on the near-horizon boundary conditions. When the mass relation (3.41) is satisfied, we use the ingoing boundary conditions (3.40) and (3.42) of Case 1 together with the additional condition (3.33)/(3.34) for the derivative of λ on the AdS boundary r = 0. The latter is needed for the finiteness of the DC conductivity and regularity of the solution, as explained above. We also normalize the Maxwell field on the horizon which completes the set of four boundary conditions needed to numerically solve a system of two second order equations for two fields.
If the mass relation is not satisfied, we use the ingoing boundary conditions (3.43) of Case 2. In this case the near horizon expansion imposes an additional relation between the horizon values of the fields A(r h ) = a 0 and Λ(r h ) = l 0 , given in (3.37). It is also necessary to impose the boundary condition (3.33)/(3.34) on the field λ at r = 0. This comprises the set of four boundary conditions needed for solving the system of equations. In distinction from previous case, we cannot anymore normalize the Maxwell field on the horizon since this would overdetermine the system. Nevertheless, the numerical solution for both the Maxwell field and the field λ is everywhere regular. However, the DC conductivity has a divergent imaginary part indicating that the real part is infinite. Hence, the values of mass parameters away from the condition (3.41) do not provide a well behaved dual description to a field theory with momentum dissipation and we shall not consider this case in what follows. One can solve numerically equations (3.6), (3.7) for a generic function M 2 (r) = r ν and extract the behaviour of the optical conductivity of the system. As already shown in fig. 1 the power ν determines whether the model is dual to a metallic or to an insulating state. The choice of M 2 (r) influences also the AC response of the dual CFT as shown in fig. 3. The value of the real part of the conductivity at zero frequency coincides with the analytic value for the DC conductivity given in (3.25). For ν 0 we observe a formation of a peak, which becomes sharply localized for higher values of ν and moves towards zero frequency. This suggests the presence of a localized excitation whose width decreases with ν very fast. This resonance shows up only in the insulating state and it is reminescent of what has been found in [11]. It would also be interesting to proceed with an analysis of the quasinormal modes spectrum as in [32] to better understand the nature of this pinned response.
We associate the appearance of the peak at these values of ν with the fact that the mass functions that are more localized near the AdS boundary give rise to an explicit source of disorder that seems equivalent to the lattice disorder discussed in introduction. We find this a reasonable interpretation because such a disorder amounts to an additional source of an explicit breaking of translations that makes the otherwise exactly massless phonons become pseudo-Goldstone bosons. This is indeed what can be seen in the solid HMGs: one can identify almost massless transverse phonon poles in cases where the breaking is localized near the horizon, while as the profile of the mass function is moved towards the AdS boundary the phonons become massive [11]. 13 Hence, the physical phenomenon related to the appearance of peaks in the electric conductivity for ν 0 is expected to be the small collective field excitations -the phonons -due to the spontaneous breaking of the translational symmetry. In the presence of charge density this can be interpreted as a polaron formation as first suggested and observed in [11] (for more on polarons, see [49] and references therein).

Elastic response
In addition to the classification of the materials into solids and fluids according to the residual symmetries that are preserved in the low energy EFT, there is another perhaps more intuitive way to distinguish them, namely, according to the type of response they exhibit under the shape deformation. In standard mechanical response theory [50,51], the magnitude that encodes the material deformation is a displacement vector u i , and the direct measure of the deformation is the linearized strain tensor, A clear distinction between fluids and solids is that they respond very differently to a constant applied shear stress (given by a traceless stress tensor T ij leaving the volume of the material unchanged). For fluids, this leads to a constant shear velocity gradient (traceless part ofu ij ) and the corresponding response parameter is the shear viscosity. For solids, instead, a small constant applied shear stress leads to a constant shear strain (traceless part of u ij ), and the response parameter is the elasticity. Thus, a practical distinction between solids and fluids is that the static shear elasticity (or rigidity) modulus is nonzero for solids and vanishes for fluids -fluids do not offer resistance to a constant shear deformation. 14 This distinction between the solid and fluid phases is exactly reproduced in HMG: it is encoded exclusively in the m 2 (r) mass parameter, which vanishes for the fluid HMGs but not for the solid HMGs. Of course, the response of different materials (and also black branes) is more complex than the simplified picture above. Materials can, for instance, exhibit both elastic and viscous (i.e., viscoelastic) response. Given that the elastic properties of HMG black branes have not yet been unveiled, we shall restrict here to the elastic response by considering only static applied stress and static deformation or strain. The full viscoelastic response can be studied by considering time dependent applied stresses, but we defer it to a separate publication [52].
In homogeneous and isotropic materials, the (static) elastic shear modulus or modulus of rigidity can be defined as the stress/strain ratio where the superscript T stands for the traceless part. Equivalently, one can extract the modulus of rigidity, G, from a Kubo-like formula that relates it to the Fourier transform at zero frequency and wavenumber of the retarded correlator as: We also note that it differs from the Kubo formula for the shear viscosity given by η = − lim ω,k→0 Im T (T ) xy T (T ) xy R /ω, cf. [53].
Once we have G expressed in terms of the stress tensor two-point functions, it is easy to apply the holographic prescription to extract it. The bulk field dual to T (T ) ij is the 14 Let us note also that fluids do have a non-zero elastic compression response. Therefore this response does not allow to distinguish between solids an fluids, which is why we shall not consider it here. Let us mention though that the compression response is encoded in HMG in the m 3 (r) mass parameter, which is indeed present both for fluids and solids (see section 1). Again, then HMG correctly reproduces that both fluids and solids have an elastic compression modulus. traceless tensor mode of the metric perturbation h T ij (r, t, x k ). Its equation of motion reads (it can be derived from (5.14) in the next section) For the constant and homogeneous (ω = k = 0) response, this equation depends exclusively on the m 2 (r) mass parameter. The retarded Greens function is extracted as usual by solving (4.3) with ingoing boundary conditions and taking the ratio of the subleading to the leading mode. The resulting response vanishes only for m 2 = 0, i.e. for the fluid HMGs. In fig. 4 we show the dependence on temperature of the shear modulus for some representative model choice which we take to be r 2 m 2 2 (r) = 1. This shows that there is a well-defined sense in which the HMG black branes enjoy a nontrivial elastic shear response and thus behave as solids. We observe that the rigidity modulus G increases with temperature, whereas most ordinary solids display the opposite dependence and G decreases with increasing T . However, in ordinary solids this behaviour occurs at roughly constant energy density while in the middle panel of fig. 4 the energy density is strongly increasing with T (which relates to the fact that the CFT degrees of freedom are inevitably in a relativistic regime). The ratio of G to the energy density , instead, does display the more standard decreasing in T form, so in this sense the result seems to be consistent with expectations from ordinary solids. We defer a more complete study of how G depends on other parameters for future work [52].

Solid and fluid CFTs
In this section we study in detail the model with two scalar fields, which provides the simplest effective description of holographic solids and perfect fluids. The two fields mass term (2.27) is the most general diffeomorphism invariant Lagrangian that one can write using two scalar fields only. All the stability requirements and holographic predictions can, therefore, be translated directly to the form of the function V . This is a great advantage of the two fields massive gravity. However, we shall see below that the phenomenology of (2.27) is less rich than that arising from the generic mass terms (A.1) that we consider in appendix A.

Stress-energy tensor and the background solution
Let us start by asking when such a theory admits an asymptotically AdS black brane solution. In distinction from the mass terms (A.1), the scalar fields with action (2.27) provide a non-vanishing contribution to the background stress-energy tensor even on the solution (2.5). The scalar fields stress-energy tensor is given by Here and in what follows we drop the distinction between the indices I and i, use V with subscripts to denote partial derivatives, e.g., V X ≡ ∂V ∂X , and set L ≡ 1. On the solution (2.2), (2.5) the quantities X and Z take valuesX = r 2 andẐ = r 4 , and the stress-energy tensor takes the form This structure exactly matches the structure of the Einstein tensor for the black brane metric, i.e. G i i = G a a −r 2 d dr 2 G a a . It thus follows that the two fields massive gravity with the action (2.27) admits an AdS black brane solution with the metric (2.2) for an arbitrary choice of the function V (X, Z). The background solution of the total action (2.1) is then completely determined by the r-dependence of the background valueV (r), with the emblackening factor given by The integration constant parametrises the position of the black brane horizon as well as its mass density and temperature. The fact that the resulting function f (r) is influenced by the graviton mass term only through the background valueV (r) means that there are infinitely many different actions that lead to the same background metric. From equation (5.5) it is easy to see that, in principle, it is possible to design a two fields Lagrangian that does not affect the background metric. However, as we show below such theories do not contain any propagating vector modes and cannot reproduce the generic mass term of the form (2.12). This is different from the four fields case with zero backreaction studied in appendix A. In order to finish the discussion of the background solution we also note that the configuration (2.5) is a solution to the equation of motion of the scalar fields for any choice of the function V (X, Z). Indeed, the equation of motion has the form: On the configuration φ i = x i , g µν =ĝ µν the equation reads and is satisfied since neither the background metric nor the functionV (r) depend on x i .

Quadratic mass terms
In order to connect to the generic graviton mass term in (2.12) let us find the quadratic mass terms for the metric perturbations around the black brane. The expansion of the scalar fields Lagrangian up to the second order in inverse metric perturbations in the unitary gauge is given by where we have used the relations It is straightforward to check that the above Lagrangian transforms as a scalar under the {t, r} diffeomorphisms which is expected from the fact that the φ a fields are absent. However, in situation when the scalar fields provide a non-vanishing background stress-energy tensor, the quadratic expansion of the scalar fields Lagrangian L φ alone is insufficient to describe the graviton mass term. Instead we shall consider the quadratic expansion of the full action (2.1) including the Einstein-Maxwell part. As expected, the vector, tensor, and scalar modes come in the quadratic expansion separately.
In a generic two fields massive gravity the r-dependence of this mass parameter can be more general than in the dRGT case (3.10) and is constrained only by the stability requirements. In particular, from the quadratic action of the vector sector (A.18) one can 15 According to these definitions, the dimensionless mass parameters m 2 i take the values 2m 2 1 = −m 2 7 = 1 see that for the mode λ to have a healthy kinetic term the vector modes mass parameter has to be positive, dV (r) dr 2 > 0 . (5.13) We also note that the functions V (X, Z) with constant background valuesV (r) are special. In this case the scalar field background contribution to the stress-energy tensor takes the form of a cosmological constant and can be reabsorbed in the definition of 3/L 2 . Therefore, the scalar fields have vanishing stress-energy tensor and do not affect the background metric. However, it also leads to the vanishing of the vector mass parameter (5.11). In the vector mode analysis of section A.1 it corresponds to the case when all the mass parameters relevant for the vector sector are vanishing, m 2 1 (r) = m 2 7 (r) = m 2 8 (r) = 0. The vector sector of such theory at quadratic level is identical to the one of pure Einstein-Maxwell theory and does not contain any propagating gravitational degrees of freedom. On the CFT side it implies the absence of momentum dissipation and infinite value of electric DC conductivity. We conclude that the two fields action cannot provide a healthy graviton vector mode without contributing non-trivially to the background stress-energy tensor. It also implies, that in order to have a holographic massive gravity with propagating vectors and finite DC conductivity, which in the meantime does not contribute to the background stress-energy tensor, one has to add more than two scalar fields.
Let us now turn to the helicity-two tensor mode of the metric. In 3+1 dimensional bulk this mode exists only for perturbations homogenous in the transverse directions and encodes the viscoelastic response of the boundary theory in the holographic description. The quadratic action for the traceless tensor mode takes the form where h T stands for any of the two helicity-two components, which we parameterise as h + ≡ 1 2 L 2 r 2 h xx − h yy and h × ≡ L 2 r 2 h xy . This action is equivalent to the transverse traceless sector of the generic quadratic mass term (2.12) with the mass parameter m 2 2 (r) given by The mass of the tensor mode (5.15) is not completely determined by the background behaviour ofV (r), and is thus independent of the mass of the vector mode and of the value of the DC conductivity. We would like to note, that the tensor mass vanishes for the two field actions that only depend on Z, i.e., for the theories describing perfect fluids. This is to be expected, since such theories possess an additional symmetry of transverse space diffeomorphisms that forbids the appearance of the mass term for the helicity two mode. This observation also agrees with an earlier finding in the two fields dRGT massive gravity in [12] that the β 2 term in the mass term (2.28) gives zero contribution to the mass of the transverse graviton. Hence, in dRGT theory the β 1 term describes solids, while β 2 term corresponds to a fluid. In accordance to the symmetry considerations of the section 2.4, the scalar sector of the theory does not contain any dynamical degrees of freedom. In particular, the mass parameters m 10 and m 12 of the quadratic action vanish, and as shown in appendix A.2 all the scalar degrees of freedom can be eliminated by using constraints and are thus non-dynamical.

Transverse phonons
In order to analyse further the stability of the two fields theory we proceed with the study of the quadratic action for the inhomogeneous vector modes in the decoupling limit and derive their propagation speeds. For this we consider the quadratic action for the scalar field perturbations π i = φ i − x i while keeping the metric to be fixed, i.e. g µν =ĝ µν . In the presence of x i dependence, the two component vector π i can be split in the longitudinal and transverse parts in the following way The quadratic Lagrangian for these components then reads where the repeated a indices are contracted with the metricĝ ab r −2 . In accordance with the full analysis of the homogeneous vector modes, the absence of ghost requires that dV dr 2 is positive. We can also read off the propagation speeds of both vector modes: The absence of gradient instabilities puts additional constraints on the function V (X, Z) One notices that the speed of propagation of the transverse phonons c T is proportional to the m 2 mass parameter, which is the one that is forced to vanish by the residual VPDiffs in the fluid case. Hence, exactly as it happens for fluids in flat space, the VPDiffs are responsible for decoupling the transverse phonons, π T i , by forcing them to be nonpropagating (the dynamics of these degrees of freedom is related to vortices in a quite interesting way, see [54]). That the vanishing of c T is equivalent to the vanishing of the tensor mode mass m 2 can be easily understood in terms of the Stückelberg trick. The reason is that, in the space part of the metric, h ij is accompanied by ∂ (i π T j) to form a gauge invariant quantity, and so the mass term for the transverse traceless mode and the spatial gradient term for the transverse phonons π T i are in fact the same thing. 16 Notice, however, that the fact that the tensor mode mass must vanish for fluids is not completely obvious from the nonlinear form of the metric potential (1.2), which is a function of det g ij . Naively, it seems that it can still give rise to a mass term for tensor modes. However, this cannot happen and the simplest physical reason is that the transverse phonons in a fluid must have vanishing spatial gradient term. In turn, this further implies that black brane solutions acquire rigid properties in solid HMGs but not in fluid HMGs.

Linear and nonlinear consistency
Let us stress that the two fields massive gravity is free from the Boulware-Deser ghost [55], i.e. that at full nonlinear level the number of unconstrained degrees of freedom is the same as at the linear level. Showing this is almost trivial if one does not go to the unitary gauge, where φ i are frozen to their vev φ i = x i , but considers it as a theory of two scalar fields minimally coupled to gravity, cf. [56]. For the most general Lagrangian that depends on the first derivatives of the scalar fields, V (X, Z), the equation of motion of the scalar fields is given in (5.6) and schematically it looks like Since this equation is of the second order in time and space derivatives even at full nonlinear level, the full dynamics of the model is completely determined by specifying the same initial data as for the linearized problem. This means that there are no additional degrees of freedom showing up at nonlinear level. The other dynamical equations of the problem, the Einstein equations, do not alter this conclusion, because the stress tensor from φ i is of first order in derivatives. It is crucial that the scalar fields φ i are minimally coupled to gravity, so that no additional degrees of freedom could arise in the gravity sector.
Note that the same conclusion is not so easy to reach if one starts directly in the unitary gauge. In that gauge, the two fields theory looks like adding a potential term for the metric V (tr(g ij ) , det(g ij )) and the only dynamical equations are the Einstein equations. Nevertheless, because of the above argument, this theory is completely free from the BD problem for arbitrary choice of V (X, Z). The theory propagates four degrees of freedom at non linear level in total. There are separate conditions to be imposed in order to insure that the two additional degrees of freedom are healthy. But this is different from the BD ghost.

Conclusions
In this paper we have analyzed theoretical consistency and phenomenological implications of a wide class of holographic massive gravity theories. These encompass all possible Lorentz violating graviton mass terms that are compatible with symmetries of the AdS black brane background solution and preserve the homogeneity and isotropy in the two spatial directions of the dual field theory.
Such theories are known to provide a holographic framework for momentum relaxation and are used as a phenomenologically viable description of transport properties in various condensed matter systems. However, only very particular massive gravity model, i.e. the dRGT theory, has been mostly exploited for this purpose so far. Such restriction is justified in the Lorentz invariant massive gravity where it is well-known that a generic massive gravity is plagued by a ghost instability. In the Lorentz violating case, however, there are many more graviton mass terms that are allowed by theoretical consistency. For the first time, in this work we have spelled out all the mass terms that are permitted by the symmetry of the black brane metric in the spirit of a similar analysis that has been carried out earlier for the Minkowski background in [5,6]. We have shown how generally covariant form of these masses can be reconstructed by the usual Stückelberg trick by using a set of four scalar fields.
A large part of the paper is devoted to the discussion of massive gravity theories where only two transverse Stückelberg fields are employed. We argue that such theories are the covariant version of the effective field theories that describe solids and fluids. As such, also the two fields holographic massive gravity can be broadly split into theories describing solids and fluids. Due to the enhanced symmetry in fluids, the set of allowed mass terms is more constrained in this case. Of particular importance for this work is the observation that the m 2 mass term is absent for fluids. We further show that this mass term implies important phenomenological consequences for the response in the dual field theory. We find that theories with m 2 mass term enjoy a non-zero response to static shear deformations thus signifying a sharp distinction between holographic solids and fluids.
We have also investigated the phenomenological consequences of the new classes of the holographic massive gravities on the electric response in the dual theory. We find that it can be parametrised by a universal mass function M 2 (r) that encodes the electric properties of the dual system. We show that this function is proportional to the m 2 (r) mass parameter in the dRGT theory. However, in our case the form of this function is free. We generalize the universal expression for the DC conductivity found for the dRGT massive gravity in [2] for this case. By assuming the particular dependence M 2 (r) ∝ r ν we find that the parameter ν determines whether the dual material has the properties of an insulator or a metal. It also regulates the appearance of the phonons in the dual theory. Moreover, we also show that for ν = −3 the DC resistivity at low temperatures scales linear with temperature. The same phenomenon has been previously observed in the dRGT theory in [14]. Let us add here that the metric potential V (X, Z) that gives rise to this behaviour scales like V ∼ √ X (or equivalently V ∼ Z 1/4 ) which is similar to the square root structure in the scalar DBI Lagrangian. Since the shape of the latter is protected by a specific reparametrization-like symmetry it might be that also the universal linear in T resistivity behaviour follows from a symmetry principle.
The case of holographic massive gravity with all four Stückelberg fields employed has never been considered elsewhere. Here we perform a thorough analysis of the homogeneous scalar, vector and tensor perturbations of these theories. We find that the problematic scalar sector can be healthy and propagate up to one helicity-0 degree of freedom if the mass parameters, as defined in (2.12), satisfy The vector sector is in general stable and contains two dynamical vector fields with two degrees of freedom each. However, we find that the effective light cones on which these two fields propagate are, in general different from each other. While the Maxwell field perturbations propagate on the background Reissner-Nordström geometry, the second vector field, λ propagates on an effective geometry determined by the graviton mass parameters. The two light cones partially coincide only in the case when the mass parameters satisfy The light cone of the gravitational vector mode, although apparently causally stable, can be wider than that of the Maxwell field and thus allows superluminal propagation velocities. Although bearing no apparent pathologies we find that the superluminal sector of the holographic massive gravities is not suitable for describing momentum relaxation in holographic framework. The reason for this is that the ingoing boundary conditions for the field λ cannot be imposed in the usual manner in this case. In the theories, where the light cones of both vector modes coincide, the electric response of the dual field theory can also be parametrized by a single mass function M 2 (r) and is identical to the two fields case.

Acknowledgments
We

A Second order action -Linear consistency
In this section we shall study the linear stability of the generally covariant graviton mass term with four Stückelberg fields described in section 2.3 and explicitly given by: For the analysis of quadratic Lagrangian it is enough to expand the gauge invariant fields H AB defined in (2.21) up to first order in perturbations: where the last term arises from the perturbation of the reference metric f AB . We shall classify the perturbations according to the scalar, vector, and tensor representations of the transverse O(2) rotation group as in (2.8)-(2.10). At quadratic level the Lagrangians of the three sectors decouple from each other and will therefore be studied separately.

A.1 Vector modes
We consider the mass term relevant for vector perturbations and assume for simplicity that the φ r dependence of the masses is given by Here i = 1, 7, 8, f (φ r ) is the emblackening factor written in a gauge invariant form, and M 2 (φ r ) is a universal mass function that is regular and non-vanishing at the horizon φ r = r h . For concreteness, we also assume that the mass function M 2 (φ r ) ≥ 0 everywhere and has the dimension of mass squared. The masses m 2 i on the right hand side of (A.4) and the powers α i are dimensionles constants.
As was clarified in section 2, vector perturbations are described by the components Since the mass term above is diffeomorphism invariant by construction it is possible to rewrite it in terms of gauge invariant fields. Indeed, we find that up to the first order in perturbations the fields H ir and H it can be written as whereπ i r andπ i t are the gauge invariant fields introduced in [12] and are defined as The Maxwell perturbation a i is diffeomorphism invariant by itself, and there is one more useful gauge invariant quantitȳ that is, however, related to the other variables. With these definitions at hand, it is straightforward to rewrite the total action (2.1) up to the second order in perturbations in a diffeomorphism invariant form (A.10) The above action depends on the Maxwell field a i and the three fieldsh i ,π i t ,π i r that are related to each other via the constraint (A.9). The vector indices of the fields are always contracted trivially, and in what follows we shall suppress them in order to simplify the appearance of the expressions.
Since one expects at most two dynamical vector fields in the theory -the Maxwell field and the helicity-one part of the massive graviton -it should be possible to integrate out two of the four vector fields. In order to consistently impose the constraint (A.9) we add it to the Lagrangian with a Lagrange multiplier λ: The action (A.10) together with the above constraint contains five vector fieldsπ t ,π r ,h, a, and λ, all of which should be treated independently. The corresponding equations of motion read The equations (A.12) and (A.13) can be used to eliminate the fieldsπ t ,π r andh from the action (A.10) and rewrite it in terms of a and λ. We notice that the system of equations (A.12) has a non-degenerate solution only in the case when M 2 (r) = 0 and 16) In order to satisfy the above condition at the horizon for a generic choice of the masses m 2 i one has to impose α 8 = 0 , Under these assumptions the det P is constant and non-vanishing everywhere. The case when det P = 0 marks one particular phase of the Lorentz violating massive gravity and has to be considered separately. When det P = 0 we can solve for the fieldsπ t ,π r andh algebraically. After substituting the solutions back in the Lagrangian (A.10) with the constraint term added we obtain: When det P = 0 we find that the Maxwell field a and the field λ decouple from each other. Moreover, λ = 0 by its equation of motion and the Lagrangian for a reduces to the case of pure Einstein-Maxwell theory. We shall not considered this case anymore in what follows.
The Lagrangian (A.18) can be used to study the stability of the theory. The first term in the above expression coincides with the corresponding term found in the holographic dRGT massive gravity in [12]. The equations of motion found by varying this action with respect to the fields a and λ can be easily recast in terms of the metric perturbations by expressing λ and its derivatives from the equations (A.12) and (A.13). The equations of motion that follow from the action (A.18) are: (A.20) In appendix B we analyse the stability of the vector modes in the high-frequency approximation. We find that the local causal structure for the field λ is determined by an effective metric g ab given in (B.2). In general, this emergent light cone is different from the light cone determined by the background metricĝ ab . In particular, if the effective light cone is wider than that of the Reissner-Nordström background, the perturbations λ propagate superluminally with respect toĝ ab . A similar behaviour can be seen when discussing the perturbations of the k-essence [46,57] and Galileons [58]. For the emergent geometry to be causally stable it has to satisfy several consistency requirements investigated in detail in appendix B.

A.2 Scalar modes
We further consider the stability of the homogeneous scalar sector modes (2.8). The relevant mass terms include where we allow for an arbitrary φ r dependence of the masses.
are two relations between the seven gauge invariant fieldsh,h t ,h r ,ā,π t ,π r ,π tr . We find these to be: These expressions allow one to eliminate the fieldsh r andh t from the action (A.27). Moreover, from (A.34) we see directly how the presence of the Stückelberg fields π t , π r introduces dynamics in the scalar sector of the theory. The previously non-derivative termsh t andh r turn into terms with second order time derivatives in the action.
We thus obtain the final second order action for perturbations in two steps. We first replace the fieldsh t andh r in the Einstein-Maxwell action (A.27) by the corresponding equations (A.34). We then add the mass term (A.21) where we replace the different components of the matrix H AB by the corresponding expressions (A.32) and (A.33). The resulting action depends on the five gauge invariant fieldsā,h,π r ,π t ,π tr . Three of them are non-dynamical and can be integrated out or act as Lagrange multipliers and impose constraints among the remaining fields. The non-dynamical fields are thē a,π t ,π tr while the two dynamical degrees of freedom are carried byh andπ r . Hence, there are at most two propagating degrees of freedom in the scalar sector of the Einstein-Maxwell massive gravity. It is a standard result in general theories of massive gravity. Around Minkowski background it is known that one of the two scalar degrees of freedom is unhealthy and propagates a ghost. This happens in both Lorentz invariant [55,59] and Lorentz violating [5,6] massive gravity theories. In Lorentz violating massive gravity theories the number of degrees of freedom in the scalar sector depends on the choice of the graviton masses: by setting some of the m i 's to zero, one can propagate either two, one or none of the helicity-0 fields. Below we identify the mass terms that are responsible for the degree of freedom count and present the main choices that lead to distinct number of propagating fields.
We shall not present the full Lagrangian here. The expression is quite long and becomes even longer after we start integrating out the non-dynamical fields. We list the different steps we take in manipulating the action and discuss the different results depending on the mass parameters m i . The actual calculations are straightforward. The first step is to integrate out the Maxwell fieldā. It enters the action linearly and quadratically without derivatives, which allows use to use its equation of motion to eliminate it from the action. The next step is to integrate out the fieldπ tr . It enters the action as: Lπ tr = 1 2 m 2 10 (r)π 2 tr −π tr F 1 ḣ ,ḣ,π r ,h,π t ,π r .
The result depends crucially on whether the mass parameter m 10 in front of the term quadratic inπ tr is vanishing or not.
In the case of a non-vanishing m 10 the equation of motion ofπ tr allows to express it in terms of the other remaining fieldsh,π t ,π r and eliminates it from the action. We then proceed by integrating out the field π t which appears in the action as Lπ t = 1 8m 2 10 (r)f 4 4m 2 0 (r)m 2 10 (r) − m 4 12 (r) π 2 t −π t F 2 ḣ ,ḣ,π r ,h ,h ,π r ,h,π r .
Further results depend on whether the combination 4m 2 0 m 2 10 − m 4 12 is vanishing or not.
• 4m 2 0 m 2 10 − m 4 12 = 0 In this case we use the equation of motion for the fieldπ t to eliminate it from the action. The remaining action contains the two fieldsh andπ r . They both appear with quadratic time derivatives and are kinetically mixed. We thus see two propagating degrees of freedom. In analogy with the usual behaviour of the helicity-0 fields in massive gravity we expect that one of them is a ghost. We therefore conclude, that the massive gravity theory (2.1) with the mass term (A.21) is unstable in this case.
• 4m 2 0 m 2 10 − m 4 12 = 0 In this case the quadratic termπ 2 t vanishes from the action and the fieldπ t becomes a Lagrange multiplier imposing the constraint F 2 = 0 on the fieldsh andπ r . This particular constraint allows to eliminate all the spatial derivatives ofπ r from the action. It means that in a high frequency limit one of the degrees of freedom has the dispersion relation ω 2 = 0 and does not propagate. The other degree of freedom is a linear combination of the fieldsh ,h,π r and is dynamical. Hence, we conclude that there is one propagating helicity-0 degree of freedom in the high energy limit and this phase can be stable.
When the mass parameter m 10 vanishes, the fieldπ tr turns into a Lagrange multiplier and its equation of motion imposes a constraint F 1 on the three fieldsh,π t ,π r . The constraint has the form Further results depend on the choice of m 12 .
• m 12 = 0 For non-zero mass parameter m 12 we can use the constraint F 1 = 0 to eliminateπ t from the action. The resulting action contains second order time derivatives of the both remaining fieldsh andπ r . This shows that there are two propagating degrees of freedom. As before, we conclude that this case is unstable since one of the degrees of freedom is a ghost.
• m 12 = 0 and m 0 = 0 In this case we cannot use the constraint F 1 to eliminate the fieldπ t . Instead the constraint provides a relation between the fieldsh andπ r . We use it in the action in order to remove the termπ tr F 1 . The fieldπ t enters the remaining action only linearly and quadratically without derivatives and can thus be integrated out by its equation of motion. In the final action there are no second order time derivatives. We conclude that there are no propagating helicity-0 degrees of freedom in this case.
• m 12 = 0 and m 0 = 0 As before we use the constraint F 1 to remove the termπ tr from the action. We further use it to express the spatial derivativeπ r and eliminate it from the action. As a result the only remaining terms that involve the fieldπ r are:π rḣ ,π 2 r ,π rh . The highest derivative term of theh field that arises in the quadratic action is (ḣ) 2 , (h ) 2 . In the high frequency limit this leads to the fact that both fields have a dispersion relation ω 2 = 0. We thus conclude that there are no propagating degrees of freedom in this case.

Conclusion about the stability
To summarise we conclude that the theory propagates two degrees of freedom and is thus necessarily unstable if 4m 2 0 m 2 10 − m 4 12 = 0 . (A.36) The high energy limit allows us to further conclude that there is at most one dynamical scalar degree of freedom in the parameter region: We note that from the above cases, (A.37) and (A.38) were obtained in the high frequency limit.
Our results presented here coincide with the earlier results for the Lorentz violating massive gravity around the Minkowski background in [5,6]. In particular, in the high energy limit it was found that there is one propagating scalar degree of freedom when m 2 0 = 0, m 2 10 = 0 and no propagating scalars when m 2 10 = 0, wherease the case with m 2 0 = 0 and m 2 10 = 0 is unstable [6] 17 . These findings coincide with our conclusions above.

B Stability analysis in eikonal approximation
Here we investigate the short wavelength dynamics of the helicity-1 fields a and λ that describe the dynamics of the vector perturbations of the Maxwell field and the metric and were introduced in appendix A.1. We analyze the action (A.18) and the corresponding equations of motion (A. 19) and (A.20) by using the eikonal approximation [60]. For this we expand the fields as a(t, r) ∝ exp(iωŜ(t, r)) and λ(t, r) ∝ exp(iωS(t, r)) and take the high frequency limit ω → ∞. In this limit only the terms quadratic in ω survive and the equations readĝ ab ∂ aŜ ∂ bŜ = 0 , and g ab ∂ a S ∂ b S = 0 , (B.1) where the indices a, b = {t, r} andĝ ab is the background metric (2.2) defining the metric on which the perturbations of the Maxwell field propagate. In turn, the metric g ab is a new effective metric defining the geometry on which the perturbations λ propagate: The physical meaning of the equations (B.1) is that at the leading order in the eikonal limit the fields a and λ can be approximated as waves with a slowly changing phaseŜ and S respectively. The surface of constant phase is the corresponding wavefront and its normal is given by the gradient ∂ a S. The effective acoustic cone of the wavefront propagation is determined by the vectors tangential to the wavefront n a ≡ g ab ∂ b S. That the vector n a is indeed tangential to the characteristic surface follows from equation (B.1). By further rewriting of the equation (B.1) we find another equation g ab n a n b = 0 , where g ab is the inverse of g ab . Hence, the wavefronts of the field λ moves as if it were light in the spacetime g ab . The effective metric g ab arises due to the scalar field condensate that fills the Reissner-Nordström spacetime and defines the local causal structure for the field λ. In general, this emergent light cone is different from the light cone determined by the background metriĉ g ab . In particular, if the effective light cone is wider than that of the Reissner-Nordström background, the perturbations λ propagate superluminally with respect toĝ ab . In such situation, there are several conditions the emergent geometry has to satisfy in order for it to be causally stable.

Hyperbolicity
Obviously, the equations of motion (B.1) have to by hyperbollic, i.e. the metricg ab has to have a Lorentzian signature:

Closed time-like curves
Due to the fact that, in general, superluminal propagation might be possible in the metric g ab one has to investigate whether closed time-like curves can possibly form. To prove the causal stability of the emergent metric g ab it is sufficient to find a differentiable function τ so that ∂ µ τ is a future directed time-like vector field (see [46,61] and references therein). The function τ then serves the role of a global time function of the spacetime under consideration. We choose the function τ to be the time t of the original Reissner-Nordström spacetime. Then the following condition has to be satisfied g ab ∂ a t ∂ b t = g tt = −2m 2 1 f α−2 < 0 .

(B.5)
Depending on the choice of α this expression can become infinite on the horizon when f (r h ) = 0 stressing the fact that the coordinate t is not a good choice for the time direction near the horizon. This is the case also for the Reissner-Nordström metric and thus we will disregard this issue here. The above condition together with (B.4) implies the following constraints on the masses: In the special case when m 2 8 = 0 this leads to m 2 7 < 0.

Light cones
In order, to compare the light cone defined by g ab with the one defined byĝ ab we consider a vector v a that is null with respect to the metric g, i.e. g ab v a v b = 0. We then investigate whether this vector is space-like or time-like with respect to the background metricĝ ab . By knowing that the vector v a is a null vector with respect to g we find v r ± = 1 2 where the plus and minus indicate the two boundaries of the light cone. It is interesting to notice, that only in the case m 2 8 = 0 the light cone of g is symmetric with respect to the time axis. One can further check that The two light cones coincide only if bothĝ ab v a determines whether the light cone of g is wider or narrower than the light cone ofĝ. We note that the mutual alignment of the two light cones for some fixed values of the masses m 2 i can flip while the radial coordinate changes from 0 to r h due to the factor of f 2(1−α) inside the square brackets. In order for this not to happen we demand that A general condition for the two metricsĝ ab and g ab to have the same causal structure is that they are related to each other by a conformal transformation Ω(r) so that g ab = Ωĝ ab . Since the light cone ofĝ ab is symmetric with respect to the time axis, i.e.ĝ ab is a diagonal metric, then this is possible only when m 2 8 = 0. The conditionĝ ab v a ± v b ± = 0 then simplifies to m 2 7 + 2m 2 1 = 0 , when m 2 8 = 0 . To analyse the mutual alignment of the two light cones for a particular choice of the masses m 2 i we have to check the sign of the norms g ab v a ± v b ± . If both g ab v a ± v b ± > 0, then the null vectors v a ± with respect to the metric g ab are spacelike with respect to the Reissner-Nordström metricĝ ab meaning that the light cone determined by the new effective metricg is wider than the light cone ofĝ. Hence, the perturbations λ can propagate superluminally.

Summary
From the comparison of the two light cones defined by the background metric and the effective metric, we first conclude that for the mutual alignment of the two cones not to flip while the radial coordinate changes from 0 to r h we need to demand that α = 1. We further find that a necessary condition for the two light cones to coincide is given in (B.11) whereas the requirements for only one of the sides of the light cones to coincide read: m 2 8 = −m 2 7 − 2m 2 1 , or m 2 8 = m 2 7 + 2m 2 1 . (B.14) In general, whenever m 2 8 = 0 superluminal propagation of the perturbations λ becomes possible. This result has been obtained in the eikonal approxmation that corresponds to taking the high frequency limit. In this limit, the phase velocity determined by the effective metric g approaches the front velocity of the propagation of λ and describes the speed of the signal propagation. Such a situation for perturbations propagating on a Minkowski background would be pathological and should be avoided. However, as has been argued previously in [46], having a superluminal propagation of perturbations on an emergent effective geometry due to a condensate of scalar fields does not lead to any causal paradoxes as long as certain conditions are met. These include the requirement of the causal stability of the emergent spacetime and of a well-posed Cauchy problem on an initial hypersurface that is spacelike with respect to both the effective metric g and the gravitational metricĝ [46]. From our analysis above we find that, the spacetime defined by g ab is causally stable if m 2 1 > 0 whereas the hyperbolicity condition is satisfied when m 2 7 < m 4 8 /(8m 2 1 ). The issue of correctly posed Cauchy problem is more involved and has not be discussed here.
We also note that, strictly speaking, the eikonal approximation assumes that the underlying theory is UV complete. The theory at hand is, however, an effective theory and is known to have a very low UV cutoff scale [62]. It might therefore be that the superluminal propagation found in the eikonal limit is an artifact of the high frequency limit and would not be present in the hypothethic UV complete massive gravity theory (for a recent proposal of such a theory, see [7]). Whether or not such a completion is possible in the presence of superluminalities in the low energy effective field theory is another subtle question that has been studied in a number of earlier works [63,64] and is beyond the scope of this work.
We therefore conclude that restricting the allowed range of the graviton masses m 2 i basing on the fact that the theory permits superluminal propagation of field perturbations would be premature. Instead we shall accept the point of view that a superluminality in an effective field theory is allowed and investigate whether this causes any inconsistencies in the dual field theory.