Massive and modified gravity as self-gravitating media

We study the effective field theory that describes the low-energy physics of self-gravitating media. The field content consists of four derivatively coupled scalar fields that can be identified with the internal comoving coordinates of the medium. Imposing SO(3) internal spatial invariance, the theory describes supersolids. Stronger symmetry requirements lead to superfluids, solids and perfect fluids, at lowest order in derivatives. In the unitary gauge, massive gravity emerges, being thus the result of a continuous medium propagating in spacetime. Our results can be used to explore systematically the effects and signatures of modifying gravity consistently at large distances. The dark sector is then described as a self-gravitating medium with dynamical and thermodynamic properties dictated by internal symmetries. These results indicate that the divide between dark energy and modified gravity, at large distance scales, is simply a gauge choice.


Introduction and outline
The cosmological data show that the Universe is undergoing a phase of accelerated expansion. See [1] for the first conclusive evidence. Although a simple model of the Universe based on a Friedmann-Lemaître-Robertson-Walker (FLRW) metric in General Relativity (GR) supplemented by a cosmological constant fits well all the observations [2,3], the actual nature of the driving force behind the acceleration remains still unclear. Upcoming probes of the dynamics of the universe [4] will be crucial to shed light on this dark energy (DE). A good deal of theoretical effort has been aimed in the last twenty years to build consistent and compelling models of DE; see [5] for reviews. 1 What the vast majority of them have in common is the addition of some new degrees of freedom to the dynamics of the universe. On the one hand, these hypothetical degrees of freedom have often been interpreted as extra components of the universe beyond baryonic and dark matter, photons and neutrinos. On the other hand, they have also been commonly regarded as part of gravity itself, modifying the behaviour of GR at large distances, in a way that is compatible with the current acceleration of the Universe.
In this work we propose a symmetry-driven approach to DE. By working with an action of four derivatively coupled scalar fields, we show that DE can originate from a medium with concrete mechanical and thermodynamic properties. These four scalars are interpreted as comoving coordinates of the medium and, moreover, they are the degrees of freedom required to describe its low energy physics. At the core of this description of DE lie the symmetries of the medium. They rule the mutual interactions of the four scalars and determine whether the medium behaves as a perfect fluid, a superfluid, an elastic solid or a supersolid.
By using this large class of Lagrangians, we will argue that the distinction that is often made between DE (as some kind of matter) and modified gravity (as a true modification of Einsten's gravity) is purely artificial. Indeed both interpretations should be understood as complementary points of view of a unique (hypothetical) phenomenon. This fact becomes transparent by choosing the adequate coordinates for each perspective, as we explain below.
The four scalars can also be viewed as Stückelberg fields that allow to restore broken diffeomorphisms in models of so-called massive gravity [7,8]. Following common usage, we use the term massive gravity (MG) to refer to any model of gravity for which non-derivative metric perturbations appear squared. We recall, however, that many such models do not contain a massive graviton (of spin two). The introduction of these fields is the key for interpreting massive (and general classes of modified) gravity models as a self-gravitating media. Thanks to diffeomorphism invariance, an adequate choice of the space-time coordinates (called unitary gauge) renders trivial the dynamics of the Stückelberg fields, allowing to draw a direct connection between the theory of self-gravitating media and MG. 2 In particular, we describe in this work the matching to the general Hamiltonian construction of MG models that was studied in [11][12][13].
Most of such MG models are Lorentz breaking in the sense that, around Minkowski spacetime, the metric perturbations are just rotationally invariant in three dimensions. From the point of view of the self-gravitating medium interpretation of MG that we explore in this work, Lorentz breaking is simply triggered by the presence of the medium itself, which generically defines at least one preferred four-vector.
By construction, the general description of self-gravitating media in terms of four derivatively coupled scalars is a well-defined effective field theory (EFT). This allows a systematic exploration of DE in terms of symmetries, with a unified and suggestive interpretation of modified and MG models. This is the take-home message of this paper. The use of an EFT framework to describe modified gravity at large distances an the acceleration of the universe, opens the possibility of confronting simultaneously a wide variety of models with cosmological observations in a simple way. 3 The outline of the paper is the following. In Section 2 we start introducing the EFT of selfgravitating media, emphasizing the role of diffeomorphism invariance. In Section 3 we continue the presentation of the EFT by building the intuition that allows to interpret the low-energy degrees of freedom as comoving coordinates. In Section 4 we write the lowest order scalars of the EFT according to the possible symmetries. In Section 5 we provide a discussion of the different kinds of media according to their symmetries and their field content. In Section 6 the unitary gauge is employed to show how the EFT of self-gravitating media can be put in correspondence with a broad class of MG models. In Section 7 we set the basis for a systematic study of the cosmology of these models by determining, for each kind of medium, the number of propagating degrees of freedom in a FLRW background. We present our conclusions in Section 9. Two appendices are also included. In Appendix A we discuss briefly the conserved currents that appear in the theory. Finally, in Appendix B we give an alternative, shortcut derivation of the number of propagating degrees of freedom.
Through the paper we will use the signature convention (−, +, +, +) and units such that = c = 1.

Derivatively coupled scalar fields and diffeomorphism invariance
Let us consider the action for a set of four scalar fields under diffeomorphisms: Φ A (x µ ), A = 0, 1, 2, 3 , whose mutual interactions obey shift symmetries: At lowest order in derivatives, the action for these fields depends on the kinetic blocks

2)
where g µν is the (four-dimensional) spacetime metric. Including the usual Einstein-Hilbert term the general form of the action at lowest order in derivatives is then given by where we define M 2 pl = 1/(16πG), being G Newton's constant; U is a smooth master function of the 4×4 matrix C AB and L m is the part of the Lagrangian describing matter fields (e.g. baryonic matter). The expression (2.3) is the leading contribution to the action of an EFT where shift symmetries enforce the fields Φ A to be derivatively coupled. 4 The scalars C AB are the only independent ones that can be built out of the four Φ A at leading order (LO) in derivatives, i.e. with a single derivative acting on each Φ A . The next-to-leading order in derivatives (NLO) corresponds to operators that enter in the action with exactly two additional derivatives. 5 These operators must then appear suppressed with respect to the LO ones by some energy scale Λ squared. In principle, one would expect Λ ∼ U 1/4 , because this is the only scale of (2.3) other than M pl (and any masses coming from L m ).
At LO, we are also allowed to multiply the curvature scalar R by an arbitrary function of C AB . However, this function can always be eliminated through a conformal transformation of the metric, 6 see [15]. After such a transformation, we obtain an action with the same structure of (2.3) plus some extra terms which can be disregarded at LO because they are of higher order in derivatives. Similarly, the Lagrangian L m can be multiplied by an arbitrary function of C AB . In fact, denoting by ψ m the matter fields contained in L m , according to the power-counting this part of the action should actually be written as f (C AB )L m [h(C AB )g µν , ψ m ], where f and h denote a set of arbitrary functions [15] . In this paper we will not study the effect of these couplings to matter fields, since our emphasis is on the interactions between the fields Φ A and the metric g µν . Therefore, in what follows we set L m to zero, postponing its study for future work.
Under these assumptions, the LO gravitational energy-momentum tensor (EMT) associated to the scalar fields Φ A , defined via the variation √ −g T µν δg µν = −2δ d 4 x √ −g U , is given by where sums over the repeated indices A and B are implicit. Varying the action (2.3), the equations of motion (EOMs) of the scalars Φ A read where ∇ µ is the covariant derivative defined with the Levi-Civita connection of g µν . Taking the covariant divergence of (2.4), and using that the Levi-Civita connection is torsionless, we obtain Therefore, the covariant conservation of the EMT (∇ µ T µν = 0) can be inferred from the scalar's EOMs (2.5). Clearly, this result would also hold if, N , the total number of scalars Φ A were different from 4. Generically, given some model of scalar fields coupled to gravity, the converse of the previous result is not true: the conservation of the EMT does not imply the EOMs of the scalars. However, in our case the two are equivalent under a specific condition. Let us consider (2.3) with L m = 0 and assume (as we will do in most of this paper) that we have four scalars Φ A . If det(∂ µ Φ A ) = 0, it is possible to choose a system of coordinates such that Then, the equation (2.6) becomes ∇ µ T µν = δ A ν δ B µ ∇ µ ∂U /∂C AB , which is nothing but (2.5) and shows that the EOMs of the scalars now follow from the conservation of the EMT. Since ∇ µ T µν = 0 is ultimately due to diffeomorphism invariance, we can say that the equations of motion are a consequence of it, provided that det(∂ µ Φ A ) = 0. 7 A more direct way of proving it uses that (2.6) can be expressed as the product Q · q, where Q is a 4 × 4 matrix of components ∇ ν Φ A and q is a 4 × 1 matrix. Then, if det Q = 0, the only solution of Q · q = 0 is q = 0, which is (2.5).
It is well-known that diffeomorphism invariance can be restored for any action of metric perturbations around some background. In four dimensions this can be done by introducing four scalars that are called Stückelberg fields. 8 By construction, these scalar fields are coupled among themselves only derivatively. From this point of view, the action (2.3) with L m = 0 is a good candidate for an EFT description of massive (and modified) gravity. The choice of N = 4 scalar fields Φ A is then motivated by the number of spacetime dimensions we assume. As the Stückelberg fields acquire a non trivial 7 An analogous argument can be applied if the number of scalar fields, N , is smaller than the number of spacetime dimensions, D. However, if N > D there would not be enough freedom to eliminate all the scalars from the dynamics through (2.7) and so the EOMs cannot be equivalent the conservation of the EMT in that case. 8 The term Stückelberg field is commonly used for a field that allows to make explicit a (spontaneously broken) gauge symmetry. These fields and the associated mechanism (or "trick"), which we will use later on, are named after E.C.G. Stückelberg, who showed how a (complex) scalar field can make an Abelian vector field massive while preserving gauge invariance [17]. See [18] for a review with various applications.
MG models (like the archetypical one of Fierz and Pauli [19]) are often presented in a way in which general covariance is explicitly broken (and in some cases Lorentz invariance too). In this work we exploit the fact that these models can generically be embedded in a covariant (and Lorentz invariant) action by the introduction of four scalar Stückelberg fields. As we have just explained, the apparent breaking of diffeomorphism invariance is just a consequence of choosing a background solution of the equations of motion for these fields. From this perspective, the action (2.3), together with a background satisfying (2.7), is an appropriate tool for the description of MG models, as we will discuss later in more detail.
The structure of (2.3) makes it convenient to assign dimensions of length to the fields Φ A . With this convention, the LO operators are dimensionless and the NLO ones have dimension of energy squared, which makes explicit the fact that they must appear suppressed by Λ 2 . We can always express the Stückelberg fields Φ A introducing new fields π A as follows: (2.8) By construction, the fields π A inherit the power-counting defined for the Stückelberg fields Φ A . Inserting (2.8) in the action (2.3), assuming that we can neglect NLO and higher order operators, and imposing that |∂ µ π A | ≪ 1, we get an action for the π A that we can expand in powers of ∂ µ or, equivalently, in powers of π A . With the condition |∂ µ π A | ≪ 1, the field redefinition (2.8) suggests that the fields π A parametrize the deviation of the fields Φ A from their background location x A . Therefore, after the interpretation of Φ A as coordinates in a medium, the equation (2.8) tells us that π A may be understood as the degrees of freedom responsible for carrying sound waves in the medium itself, i.e. phonons which can be understood as the Goldstone bosons of broken translations, see e.g. [20,21]. We emphasize that these fields need not be small; it is only their derivatives that have to be under control for the validity of the derivative expansion. The fields π A are the fourdimensional generalization of the standard displacement field used in fluid dynamics, elasticity and strain theory. In the next section we start exploiting this connection, which we will phrase directly using the Stückelberg fields, using the language of the EFTs of solids, fluids and superfluids; different aspects of which have been studied in [9,10,14,15,[22][23][24][25][26][27][28][29][30][31][32]. As a matter of fact, the basic idea of using three comoving coordinates to describe fluid dynamics is nothing but the earlier pull-back formalism, which had already been used to describe the dynamics of continuous media, see e.g. [33][34][35][36][37][38].

Self-gravitating media
We split the scalars Φ A into two distinct sets: {Φ 0 } and {Φ 1 , Φ 2 , Φ 3 }, which are characterized as follows. The elements of {Φ 1 , Φ 2 , Φ 3 }, will be referred to as spatial Stückelberg fields, and we will use small latin letters a, b, · · · to index them. We require the 3 × 3 submatrix of (2.2) whose components are to be positive-definite, so that its three eigenvalues are positive. As we will soon see, the reason for this condition is allowing the inverse of B ab to define an Euclidean metric in IR 3 . As anticipated in the previous sections, we can interpret the fields Φ a as the (spatial) comoving coordinates of an internal space F that describes a continuum medium whose topology we assume to be IR 3 . This picture -that has been known for a long time, see e.g. [34]-allows to write an unconstrained variation principle for a fluid. A pedagogical review of this pull-back formalism is provided in [35] and various applications specific to cosmology can be found in [10,15,27,30,[39][40][41][42][43][44]. Here, we will follow [10] to describe the basics of this construction. We exploit the fact that, given a slicing and a threading of the spacetime (Lorentzian) manifold M, we can formally write M = M T × M S , where M T and M S are respectively the sets of time, {x 0 }, and space, {x 1 , x 2 , x 3 }, coordinates on M. By choice, M S is endowed with a Riemannian metric. Provided that at any given time x 0 , the condition det(∂ i Φ a ) = 0 is satisfied, there exists an invertible map of M S into the medium F, which we assume to be smooth: . While the first of these maps identifies the fluid (or in general, medium) element (labelled by Φ a ) that sits at time x 0 at the point x i ; the second map describes the trajectory of the medium element Φ a in spacetime [10]. This is nothing but the dual description of a continuum in terms of Lagrangian (Φ a ) and Eulerian (x i ) coordinates.
The spacetime metric g µν induces a three-dimensional metric on F, which is given by the matrix B ab defined as follows: The induced metric B ab is the inverse matrix of B ab defined in (3.1), i.e. B ab B bc = δ a c . The condition that B ab is positive definite ensures that the induced spatial metric on the medium F is Riemannian.
Since the fields Φ a can be interpreted as comoving coordinates, they must remain unchanged along the fluid flow. This allows to define a (unique) four-velocity u µ , through the conditions whose only solution is and B B B denotes the 3 × 3 matrix whose components are given by (3.1). In what follows, we will always use boldface capital letters for three-dimensional spatial matrices. When pulled-back into the spacetime manifold M, the induced metric B ab becomes a projector H µν on the embedding of F in M: Only if u [α ∇ µ u ν] = 0 (where the brackets indicate full antisymetrization of the three indices), Frobenius theorem guarantees that u µ is hypersurface orthogonal. In that case there is a one-parameter family of hypersurfaces that are orthogonal (everywhere) to the four-velocity u µ , which thus defines a slicing of the spacetime manifold M. So far we have only dealt with the spatial Stückelberg fields, describing the structure that was used, e.g. in [10] to study (a certain class of) perfect fluids. We will now deal as well, with the temporal Stückelberg, Φ 0 , for which we impose the condition This condition allows to define another four-velocity by means of the derivative ∂ µ Φ 0 , which is thus constrained to be time-like: This defines a hypersurface at each point that is orthogonal to V µ and has an induced metric Since V µ is the gradient of a scalar, its covariant vorticity ω µν = p α µ p β ν ∇ α V β vanishes. In Minkowski and FLRW spacetimes, the four-vectors u µ and V µ coincide, up to a sign. However, if these spaces are slightly distorted, the two velocities generically differ from each other. Given the four scalar fields Φ A , the four-vectors V µ and u µ are the only independent four-vectors that are invariant under the group SO(3) s of internal spatial rotations and under the shift symmetries (2.1).
The condition (3.7) also allows an interesting extension of the maps between M S and F that we mentioned earlier. In particular we can extend the medium space F, spanned by the three spatial Stückelberg fields Φ a , to a space F 4 spanned by the full set {Φ 0 , Φ a }. By construction, and in analogy to (3.2), the matrix C AB , defined through is a Lorentzian metric in the extended medium space F 4 ; and is the inverse of the 4 × 4 matrix of kinetic blocks introduced in (2.2). Hence, the map carries the spacetime manifold onto the extended (self-gravitating) medium space, whereas the inverse gives the location of each element of the extended medium (labelled by Φ A ) on spacetime, as described by the coordinates x µ . In other words, the second of these two maps describes the world hyper-volume of the extended medium F 4 as it propagates on the spacetime M. The role of the temporal Stückelberg Φ 0 is therefore two-fold. First, allows to define an independent time-like four-vector via (3.8). As we will later see, this will give us the possibility of defining two-constituent fluids, as required to describe superfluids (provided that the appropriate symmetries are satisfied) [25]. And second, as we have just described, Φ 0 can also be naturally interpreted as a comoving time coordinate in the medium space, which allows to conceive the self-gravitating media described by (2.3) as four-dimensional objects propagating in spacetime, extending the pull-back formalism to include the temporal dimension as well. Both pictures play important roles in the novel interpretation of MG and DE as self-gravitating media that we explore in this work.
In the next section we focus on the operators that can appear in the action (2.3) at LO in derivatives, and in Section 5 we will classify the different media that arise combining them.

Spatial SO(3) invariance and additional symmetries
For simplicity, we will always impose a global internal spatial SO(3) symmetry in the medium space, that we denote by SO(3) s , so that all the operators in the action are invariant under This basic symmetry will also allow us to make direct contact with the general models of massive gravity classified in [13], as we will discuss in Section 6.
To construct a LO complete set of independent SO(3) s -invariant scalars we start with the traces τ n = Tr(B B B n ) , n = 1, 2, 3 . We recall that the SO(3) s -invariant determinant b = √ det B B B, defined already in (3.5), is not independent of these traces because it can be written as b 2 = (τ 3 1 − 3 τ 1 τ 2 + 2 τ 3 )/6. Introducing the 3 × 3 matrix Z Z Z of components we can define the scalars which are also invariant under SO(3) s . The traces of Z Z Z 2 and Z Z Z 3 (and higher powers) do not need to be considered once y 0 is included because they can be written as powers of the latter. Besides, since the four-vectors ∂ µ Φ 0 and u µ are invariant under SO(3) s , there are two other scalars at LO that are invariant under this group. One of them is X = C 00 , which we already introduced in (3.7), and the other is This closes the set of independent LO SO(3) s -invariant operators, which is therefore composed of nine elements and reads: We point out that the determinant of the 4 × 4 matrix of components C AB , defined in (2.2), is also invariant under SO(3) s , but it is not independent of this set because it can be expressed as What is interesting of this operator is that is closely related to the pull-back of the medium volume element in spacetime. Indeed This expression allows to pass from Eulerian to Lagrangian coordinates (and viceversa) in four dimensions, extending the three-dimensional treatment of these coordinates that was discussed in [10].

Additional symmetries at LO
We can further constrain the form of the action (2.3) by requiring additional symmetries, leading to different types of self-gravitating media. Two examples are elastic solids and perfect fluids, which are derived from symmetries imposed on the spatial sector of the Stückelberg fields. We will see in Section 5 that by imposing symmetries that mix the spatial and temporal Stückelberg fields, new kinds of materials arise. For reference, we summarize in Table 1 the notation for the operators that we use in this paper. Let us now make a list of the additional symmetries of interest for us: • Imposing volume-preserving diffeomorphism invariance of the spatial sector: the scalars τ n have to combine in the action into the determinant operator b and, in addition, the scalars y n are automatically forbidden. Therefore b, X and Y are the only operators allowed at LO. This symmetry was identified in [9] in the context of the EFT formulation of perfect fluids containing only b.
We will also consider a series of symmetries that were proposed in [7]: • If we impose that the action has to be invariant under the transformations for arbitrary f a , we find that at LO the action will generically depend only on X and on the matrix of components Then, the invariance under SO(3) s further enforces W ab to appear in the action only through the traces which are non-linear combinations of the scalars X, τ n and y n : Under the transformation (4.10), Therefore, the operator Y measures the failure of Φ a to remain comoving along u µ when such a transformation is applied. Besides, as in the case of V s Diff, at NLO the same operators with derivatives of u µ are also allowed under (4.10).
• If the action is invariant under the only LO operators that are allowed are Y and the traces τ n . Clearly, u µ is invariant under this transformation and therefore all the NLO operators mentioned above constructed from derivatives fo u µ also respect this transformation.
If we combine the symmetries (4.14) and (4.9) the master function U can only depend on Y and b and the resulting medium has the EMT of a perfect fluid, as we will later see.
• Finally, if the action is invariant under 15) it is constrained at LO to be a function of the traces τ n , w n and the scalars The operators O αβn present some peculiarities with respect to the ones that we have encountered so far. From the point of view of the symmetries, α, β ∈ IR, which implies that there must be an infinite (and uncountable) amount of these operators. However, since all y n vanish at zero order for diagonal metric backgrounds, a well defined perturbative expansion seems to require β to be a positive integer. Even then, there is still an infinite number of these operators at LO, which means that the master function U cannot be arbitrary if the action is to be finite. The same difficulty is also manifest on the associated Noether currents, see Appendix A.2.
If in addition to (4.15), we impose also (4.9), the operators (4.16) get restricted to the subset with β = 0. For later convenience we will denote them as follows: (4.17)

Four-dimensional media
Symmetries of the action LO scalar operators Type of medium  If instead we combine (4.15) and (4.10), the operators w n get selected.
In Table 2 we summarize the different LO scalar operator content according to the internal symmetries of the action and the physical interpretation of the resulting system, which we discuss in the next section. In addition, media with reduced internal dimensionality are possible as well; see Table 3. For instance, if the spatial Stückelberg fields Φ a are not present in the LO action, the only operator at that order would be X.

Lagrangians for self-gravitating media at LO
Given the field content: the metric g µν and the scalar fields Φ A , it is possible to study in a systematic way the models resulting from imposing certain symmetries in addition to SO(3) s and internal shifts. Each class of symmetries corresponds to a medium with specific mechanical and thermodynamic properties, which can be obtained from the EMT. For all the media we consider we require spacetime diffeomorphism invariance. Then, the most general action at LO for a medium described by the four Stückelberg fields Φ A can be constructed in terms of the scalar invariants X, Y , τ n and y n : and the corresponding EMT is given by where O k are the nine scalar LO operators appearing in (5.1) and we use the notation Their partial derivatives with respect to the inverse spacetime metric are where the indices µ and ν are symmetrized in the second term of the last expression. We have included in (5.3) the operator b, that can be written as a combination of the three τ n , for later convenience. The dot (·) represents the standard three-dimensional matrix product and we have introduced the notation In the expressions (5.3), C and ∂ µ Φ have to be understood as a 1 × 3 matrices of components C i and ∂ µ Φ i , respectively. It is also useful to note that In what follows we list the media that arise according to the symmetries that we have discussed and, also, according to the assumptions on the low-energy field content. We focus first on the systems that contain only a partial subset of the Stückelberg fields, which we call temporal or spatial media. These are the solids U (τ n ), and two different types of perfect fluids: U (b) and U (X); the first of which is a special case of the solids. Then we will move on to the media that at low-energy require all the Stückelberg fields. As discussed below, the properties of the most generic class of them have been argued to define (non-relativistic) supersolids [45]. A subclass of these (with an enhanced symmetry group: internal three-volume preserving diffeomorphisms V s Diff) would then describe superfluids. Other media, defined by specific additional symmetries, will also be mentioned, see Table 2. The connection with modified and, specially, massive gravity models will be studied in Section 6. Then, in Section 7 we will discuss the number of propagating degrees of freedom in FLRW, which strongly depends on the symmetries beyond (2.1) and (4.1) that are imposed.

Media with reduced internal dimensionality
We start by describing three special cases for which some of the four Stückelberg fields are missing from the LO action, while the shifts symmetries (and eventually while SO(3) s ) are respected. We refer to these media as temporal or spatial depending on whether the missing fields are the spatial or temporal Stückelbergs, respectively. Concretely, we will discuss two kinds of spatial media (solids and perfect fluids) and the only kind of temporal medium (an irrotational perfect fluid) that exists at LO with our symmetry assumptions. A peculiarity of the spatial and temporal media is that there is no symmetry that forbids one subset of Stückelberg fields at LO, but allows them at a higher order in derivatives. This means that the defining restriction on the low-energy field content of these models must be taken as an assumption, which in principle cannot be justified from the point of view of the symmetries.
Let us first consider the possibility that the only operators relevant at LO are the traces τ n . Then, the action contains only the three spatial Stückelberg fields Φ a . This necessarily assumes that no other fields are needed to describe the system at sufficiently low energies. As we just mentioned, there is no symmetry that forbids Φ 0 at LO but allows it at some higher order in derivatives. Notice that by looking at Table 2, it would seem that combining the operators τ n are selected. However, it is clear that this cannot be considered a true symmetry of an action based on U (τ n ) since Φ 0 is not part of τ n .
In addition, given that the general action (5.1) lacks potential terms and it generically mixes all fields in a democratic way, it is difficult to imagine how Φ 0 could be (on average) frozen at some constant value only at LO, if it is not neglected also at higher orders. Therefore, it seems that if Φ 0 does not appear at LO, it should be completely absent from the action and this must be imposed as an assumption on the physical nature of the system.
The solid is thus the most general medium that we can construct at LO with only the three spatial Stückelberg fields, respecting the basic symmetries that we assumed (SO(3) s and shift invariance). The corresponding EMT can be written as Such a medium can be interpreted as the relativistic generalization of an elastic material, in the sense that this kind of media can support compressing pressures; see e.g. [46]. In standard elasticity theory in a flat space spacetime [47], the stress state is described by the spatial internal inverse metric B ab , so that the medium is relaxed (not stressed) if B ab = δ ab . The object B ab is the covariant generalization of the standard three-dimensional left Cauchy-Green deformation tensor; and thus ∂Φ a /∂x i is the deformation gradient tensor, see e.g. [48]. This type of medium, which has been commonly called solid -see e.g. [20]-or elastic solid -to distinguish it from a perfect fluid-has been applied in astrophysics and cosmology for the description of the dynamics of the interior of neutron stars [49], in a proposal for a dark matter model [50] and also for constructing a model of primordial inflation [27,39].
Solids posses two features, due to the independent operators τ n , that make them different from perfect fluids -which we discuss below-and interesting for certain applications in the context of cosmology. First, they exhibits anisotropic stress, due to the pressure perturbations being in general dependent on the spatial direction. And second, the spatial trace operators τ n are the only ones among the nine LO scalars operators (4.6) that generate a mass for the (spin two) graviton, as we will see in Section 6. Clearly, these two properties are not exclusive of the solids with EMT (5.6) -see Table 2-but these are the simplest systems that display them.
Although it is not immediate to identify a four-velocity from the EMT (5.6), the natural choice is u µ , defined in (3.4), since this is the only SO(3) s -invariant four-vector that is available when only the three spatial Stückelberg fields are present. With this choice of frame, the energy density and the pressure of the solid, defined as usual through the projections ρ = T µν u µ u ν and 3p = T µν H µν = T µν (g µν + u µ u ν ), are given by The four-velocity u µ corresponds to the energy frame -sometimes called rest frame-of the solid, which for an arbitrary EMT is defined (if it exists) as the four-vector U µ (r) that solves g µν + U µ (r) U ν (r) U γ (r) T νγ = 0, see e.g. [10]. Then, the EMT (5.6) can be written as follows: where the non-zero anisotropic stress π µν can be obtained, for instance, from the difference of (5.9) and (5.6). The expression (5.9) is what in cosmology is usually referred to as an imperfect fluid.
The symmetric matrix B ab can be decomposed in terms of its eigenvectors ζ a (n) and (real) eigenvalues λ n , n = 1, 2, 3 as B ab = n λ n ζ a (n) ζ b (n) . Then, we can define the projectors P ab satisfying P (n) · P (m) = δ mn P (n) . The eigenvectors are mutually orthogonal with respect to the three-dimensional metric δ ab , i.e. δ mn = ζ a (m) ζ b (n) δ ab . By inspection, U ab is diagonalized by ζ a (m) as follows: Therefore, the EMT (5.6) can be cast in the following form ν } form an orthonormal tetrad. The elastic solid has principal pressures given by The thermodynamics of solids can be studied following analogous lines to those given below and in [51] for the special case of perfect fluids.

U (b): perfect fluid
If the SO(3) s symmetry of the solid is enlarged to the invariance under volume preserving spatial diffeomorphisms V s Diff, see (4.9), the trace operators τ n must combine into the determinant b 2 defined in (3.5), constraining further the properties of the medium. Indeed, we recall that b 2 can be expressed as a function of τ n using b 2 = (τ 3 1 − 3 τ 1 τ 2 + 2 τ 3 )/6. Then, under the assumption that only the spatial Stückelberg fields are present and imposing V s Diff, the master function U that determines the Lagrangian of the system depends only on b at LO. However, if the restriction on the field content is relaxed, allowing also for the presence of the temporal Stückelberg Φ 0 , the operators Y and X should be included as well under at lowest order in derivatives, see Table 2. In this section we will focus on the special case U (b) and we will later move progressively to the general structure of the media that respect V s Diff.
In the case U (b) the EMT reads with u µ being defined in (3.4). In comparison to the general solid (5.9), the symmetry V s Diff removes the anisotropic stress π µν , thus restricting the medium to be a perfect fluid. In addition, it also sets to zero the (spin two) graviton mass. The application of this type of systems for describing cosmological perfect fluids (and specifically their perturbations) was studied in [10]. For instance, if the fluid has a barotropic equation of state of the form p = wρ, with constant w, the master function is of the form (5.14) Hence, the background dynamics of the ΛCDM universe can be modelled with with Ω r ∼ 10 −4 and Ω m ∼ 0.26, to satisfy the current constraints on the radiation and (cold dark) matter densities [3].
In general, the perturbations of a U (b) perfect fluid can be studied decomposing the spatial phonon fields π i , i = 1, 2, 3 -introduced in (2.8) for the general four-dimensional case-into longitudinal and transverse polarizations: where the transverse fields π i T couple to vector metric perturbations and their evolution is dictated by the conservation of vorticity [10], which is a consequence of the symmetry V s Diff. As we discuss in Appendix A, this can be easily obtained applying Noether's theorem, which implies and infinite set of covariantly conserved currents [9,10].
A peculiar property of the perfect fluids of the kind U (b) is that the transverse modes π i T , which are interpreted as vortices, do not propagate in flat space since the symmetry V s Diff forbids spatial derivatives in their action. Their time evolution is simply given by the equationπ i T = 0. This has lead to the suggestion that these fluids may be strongly coupled at all scales [22]. However, it has been argued that the issues associated to the evolution of these modes (e.g. the appearance of divergences in the scattering of phonons in flat spacetime) can be resolved with a careful choice of the appropriate (physical) observables [31].
It can be easily checked, using a field redefinition, that the inclusion of the NLO operators respecting V s Diff and containing only the three spatial Stückelberg fields does not modify the evolution equation of π i T in Minkowski (with respect to the one at LO) [10]. There are five independent operators at NLO: [10]. Once these operators are are included, the EMT of the system is no longer that of a perfect fluid.
The relativistic formulation of self-gravitating media that we are exploring in this paper allows a straightforward connection to the thermodynamic theory of continuous media. This can be done by constructing a dictionary between two EFT pictures: the action (5.1) and the fundamental equations of thermodynamics, see e.g [23]. A more detailed discussion of this dictionary is beyond the scope of this work and and we provide it, for the case of perfect fluids, in [51]. For the purpose of illustration we use here the perfect fluid U (b). Essentially, we need to match the operator b with an intensive thermodynamic variable (or a combination of them). Choosing the particle number density n and the entropy density s as independent variables, we write b = b(n, s). Then, the chemical potential µ and the temperature are obtained from the first principle of thermodynamics as µ = ∂ρ/∂n and T = ∂ρ/∂s. The pressure enters through the Euler relation µ s + n T = ρ + p , (5.17) which must be valid for any U (b). This can be achieved, for instance, with the identification b = n (and then T = 0), which implies µ = −U b , and is in principle suitable for a degenerate system (i.e. in the limit T → 0). A different option is b = s (with µ = 0) and T = −U b , as for a gas of photons. This second choice is the one that was advocated e.g. in [10,22,23].

U (X): perfect fluid
At the opposite side of the spectrum from solids and perfect fluids lie the self-gravitating media that can be described at LO by introducing only the temporal Stückelberg field Φ 0 . We have not found a symmetry that forbids the spatial Stückelberg fields at LO but reintroduces them at NLO or at higher orders in derivatives. This is analogous to what happens for the solids and perfect fluids studied above, but exchanging the roles of the temporal and spatial fields. By looking at Table 2, it seems that V s Diff and Φ a → Φ a + f a (Φ 0 ) select X, but this combination of transformations is not a symmetry of an action based on U (X), simply because Φ a are not in X. Therefore, the absence of these fields must be assumed. Clearly, with such an assumption, the function U (X) is the only shift-symmetric possibility that exists at LO (with a single scalar field). Further restrictions involve necessarily discrete symmetries such as e.g. X → −X.
Then for a generic U (X), The EMT is given by This type of perfect fluid is fundamentally different from the case U (b) because its four-velocity, V µ , is the derivative of a scalar and thus it has zero vorticity. This is due to the fact that for U (X) the transverse phonons π i T are entirely absent by construction. Actually, since both U (X) and U (b) describe perfect fluids, their phonon expansions (2.8) can be matched setting to zero the transverse modes of the latter and identifying through an equation the divergence of π i with the time derivative of π 0 . This matching generically involves as well a linear combination of the scalar metric perturbation. The correspondence can also be easily established at the level of the background. For instance, a constant barotropic equation of state is obtained choosing U ∝ X (1+w)/(2w) , to be compared with (5.14). For an early study of irrotational perfect fluids where longitudinal fluctuations were already interpreted as phonons and a thermodynamic dictionary was provided, see [52].
In order to extend beyond LO the medium that can be constructed using only the temporal Stückelberg, we must linearly add to the action the operators (∇ µ V µ ) 2 and ∇ µ V ν ∇ µ V ν at NLO. Clearly, the operator ∇ 2 Φ 0 ∼ ∇ µ V µ can be omitted assuming appropriate boundary conditions, since it is a total derivative. Besides, the NLO operators ǫ µναβ ∇ α V µ ∇ β V ν and V µ V ν ∇ µ V α ∇ ν V α need not be included thanks to the fact that V µ is hypersurface orthogonal. The first of the two vanishes and the second one can be written as a combination of (∇ µ V µ ) 2 and ∇ µ V ν ∇ µ V ν .
If we assume that the operator X is absent, and we work only at NLO, the resulting model is a special case of the "Einstein-aether" model [53], see e.g. [54].
Notice that the extension to all orders of U (X) is not a generalization of the ghost condensate [55], but instead the most general modified gravity model based on a single scalar derivatively coupled. Actually, including all the orders in the derivative expansion corresponds to a large class of models contained in the EFT of inflation [56] (or DE [57]) for which the (soft) breaking of the temporal shift symmetry Φ 0 → Φ 0 + c 0 (with constant c 0 ) due to a potential term for the inflation (or, generically, the extra scalar mode) is entirely neglected.

Four-dimensional media
We will now move on to describe the media of Table 2, which pertain to a different class than the cases we have studied so far, because all of them contain the four Stückelberg fields Φ A , which motivates the name we give them. As we will now see, the most general of these media (whose action is invariant under shifts and internal spatial rotations) is the covariant generalization of a non-relativistic model used in [45] to introduce supersolids. If we further constrain these systems by requiring invariance under internal spatial diffeomorphisms that preserve the volume, we get superfluids [23]. Symmetries that mix the temporal and spatial Stückelberg fields -see Table (2)-lead to a set of particular subcases of the general supersolid.
Let us start with the most general media that our symmetries allow. The less symmetric selfgravitating media that we can construct with our assumptions respect only the internal shifts (2.1) and the spatial rotations SO(3) s of (4.1). Therefore, their LO master function U contains all the operators of the set (4.6), that we identified in Section 4, and thus their diffeomorphism invariant action is precisely (5.1). A non-relativistic description of systems characterized precisely by these symmetries was given in [45], where they were interpreted as zero-temperature supersolids. According to [45], the three spatial Φ a would correspond to the comoving coordinates of the medium (as we have been interpreting them), whereas Φ 0 would be the (shift-symmetric) phase field related to a U (1) symmetry associated to particle number conservation.
In spite of some claims, see e.g. [58], whether actual supersolids exist in nature is still unclear. For the purpose of this work we will simply use the term to refer to a continuous medium amenable to a coarse grained description with four degrees of freedom and based on the shift and rotational symmetries (2.1) and (4.1). For a review on supersolids we point the reader to [59].
The same set of symmetries, equations (2.1) and (4.1), was later considered in [29], where the CCWZ method [60] was applied to derive (for Minkowski spacetime) a relativistic version of the LO action given in [45]. However, the action written in [29] for these symmetries depends on a generic function of all our kinetic blocks C AB except X, which is an independent invariant operator under internal translations and SO(3) s rotations. Besides, the concrete operators y n , defined in (4.4) were not given in [29].
Phenomenologically, the media of this kind posses many interesting properties. First of all, they share with the solids U (τ n ) the fact that the graviton acquires a non-vanishing mass that is related to the presence of anisotropic stress and to the speed of propagation of transverse modes. We refer to the group velocity dω/dk, which given e.g. a dispersion relation of the form ω 2 = k 2 + m 2 , clearly involves the mass m. In addition, given that Φ 0 and Φ a lead to two different four-vectors in the medium, the four-velocities (3.4) and (3.8), the energy-momentum tensor exhibits a non-zero "heat-flux".
We point out that the object C AB , defined in (2.2), is the four-dimensional generalization of the left Cauchy-Green tensor, see Section 5.1.1. This agrees with the picture in which the general media described by the EFT (5.1) can be interpreted as four-dimensional media (or hyper-volumes) propagating in spacetime. In Section 6 we will make the connection between these media and massive (and modified) gravity, using the unitary gauge as tool to simplify the matching.

U (X, Y, b): superfluid
Imposing that the action (5.1) has to be invariant under V s Diff, the internal volume preserving diffeomorphisms (4.9), the only invariant LO operators that appear in the master function are X, Y and b. The corresponding EMT is given by In general this is not a perfect fluid, due to the four-velocities V µ and u µ not being parallel. The fact that V µ has zero vorticity has motivated interpreting these media as relativistic superfluids [23,25], following earlier ideas for (non-relativistic) supersolids [45] and superfluids [61]. The key idea of this picture is that the EMT (5.19) can be seen as a two-component fluid whose admixture allows to describe the various phases of a superfluid. For instance, it is straightforward to show that the phonon expansion of the action leads to two kinds of longitudinal waves (propagating with different speeds of sound), which suggests interpreting them as the first and second sound in superfluids [25].
It is worth highlighting some special cases of the general superfluid: • If in addition to the symmetry V s Diff we impose also invariance under Φ 0 → Φ 0 + f (Φ 0 ), the master function has to be a function of the operators O α , which we defined in (4.17), and b. The operators O −1 and b, together with the NLO operators that can be constructed from V µ , were used in [40] to model Lorentz breaking in dark matter.
• The master function U (b, Y ) describes a perfect fluid with four-velocity u µ and energy density and pressure This case is interesting because it comes from imposing not only the symmetry V s Diff but also requiring invariance under Φ 0 → Φ 0 + f (Φ a ), see Table 2. This feature makes it qualitatively different from the other perfect fluids that we have encountered: U (X) and U (b), whose LO structure does not derive from symmetries acting on the four Stückelberg fields. One can obtain a constant barotropic equation of state for U (Y, b) choosing where U is an arbitrary function. Another choice of master function which also leads to a constant barotropic equation of state is U ∝ (b 1+w + Y 1+1/w ).
• An equation of state w = −1 can be obtained from (5.21) by choosing which interestingly corresponds to the enhanced symmetry By looking at equation (4.8), it is now clear why b Y is the operator that allows to switch from Eulerian to Lagrangian coordinates.

Special supersolids
These media are the subclasses of the general action (5.1) listed in Table 2.
• For instance, imposing the symmetry Φ 0 → Φ 0 + f (Φ a ), the leading invariant operators are τ n and Y and the EMT is Given what we have seen so far, a natural interpretation of (5.24) seems to be that it describes a solid coupled to a perfect fluid. An analysis of the thermodynamics of perfect fluids [51] suggests that Y can be interpreted as the temperature so that U (τ n , Y ) may be a possible description for a solid at finite temperature.
• The symmetry Φ a → Φ a +f a (Φ 0 ) enforces the master function to be of the form U = U (w n , X). If in addition we impose the symmetry Φ 0 → Φ 0 + f (Φ 0 ), the operators w n are selected. We will comment on these cases in the next section, where we build the connection between massive gravity and self-gravitating media.

SO(3) spatially symmetric massive gravity
The main motivation for massive and modified gravity in the context of cosmology is the possibility that the observed acceleration of the Universe could be due to a modification of GR that weakens gravity at large distances. In this respect, the idea of massive gravity (MG) is appealing, since endowing the graviton with a mass could effectively make gravity a short range interaction (on cosmological scales). The first attempt to provide a mass to the graviton dates back to 1939; when Fierz and Pauli [19] studied the most general Lorentz invariant mass term for metric perturbations around Minkowski space. It was thought for a long time that any extension of the Fierz-Pauli Lagrangian beyond the the quadratic level was pathological. Indeed, Boulware and Deser argued that at the non-linear level a ghost mode would have necessarily be present [63]. The problem was solved in [64] by finding a ghost free extension of the Fierz-Pauli Lagrangian at the non-linear level [64,65]. This model, which exhibits Lorentz symmetry when expanded around Minkowski, solves the ghost issue at the (heavy) price of lacking spatially flat homogenous FLRW solutions [66]. As shown in [67] the Fierz-Pauli model as an effective theory (and also [64]) has a (very low) cut-off Λ 3 = (m 2 M pl ) 1/3 , where m sets the scale of the graviton mass mass scale. A recent proposal to improve this cut-off has been given in [68,69].
Fierz-Pauli MG fails to reproduce the correct light bending around heavy and dense objects, in sharp contrast with standard GR. This also happens in the limit of vanishing graviton mass, an issue that is known as the vDVZ discontinuity [70]. A possible way out was proposed by Vainshtein in [71], by arguing that non-linear effects restore the correct GR behaviour at short distances from a gravitational source. Whereas this mechanism has been studied for various models [72], it is important to stress that it relies on strong non-linearities; even at the macroscopic scales of the Solar System where the value of the gravitational potential is small.  Table 4: Relationship between material Lagrangians and massive gravity models.
A possible way out of this last difficulty is to abandon Lorentz invariance in the gravitational sector 9 in favour of a simpler symmetry group; in practice, by requiring only spatial rotational invariance [7,74]. All Lorentz-breaking MG models with five DoF having (at least) spatial rotational invariance are free of ghost instabilities and can be classified [11,12]. In addition, they have no vDVZ discontinuity and FLRW solutions do exist for these models [75]. Their range of validity is given by an ultraviolet cut-off Λ, which is dictated by the symmetries that define the EFT. Typically, for Lorentz-breaking MG the cut-off is Λ 2 = (m M pl ) 1/2 ≫ Λ 3 [7,8,74,75]. In order to have any impact on the current acceleration of the Universe, the scale m should thus be of the order of today's Hubble scale, giving Λ −1 2 ∼ 0.1mm. The EFT of four derivatively coupled scalar fields with internal spatial SO(3) s invariance allows to re-interpret Lorentz-breaking MG models as self-gravitating media. From the perspective of this EFT, the breaking of Lorentz symmetry is simply a consequence of a background choice. MG is then a class of models within a broad framework. In other words, the EFT of self-gravitating media leads to a systematic construction of models for the acceleration of the Universe, including MG.
The basic idea that helps to build this connection is the Stückelberg mechanism or trick, which we already mentioned in Section 2. The application of the Stückelberg mechanism for Fierz-Pauli MG was already given in [76] and was later generalized to other models of MG in [67] and [7]. If h µν is a metric perturbation around a (reference) background metric,ḡ µν , such that the full metric is g µν =ḡ µν + h µν , a generic (Lorentz invariant) MG potential may be expressed as function of the perturbation h µν with indices raised using g µν , i.e. √ −g V (h µν , g µν ). Given such a potential, added to the Einstein-Hilbert Lagrangian, the Stückelberg trick can be implemented replacing the background metricḡ µν with a tensor field made out of four scalars: , so that the metric perturbation h µν is replaced with g µν − G µν , see e.g. [77]. Using these replacements to express the action in terms of the spacetime metric g µν and the tensor G µν we obtain a covariant embedding of any Lorentz invariant MG model originally defined with a reference metricḡ µν . This procedure can be adapted to deal also with Lorentz breaking models, as described in [7].
In this work, we started our analysis directly with an explicitly diffeomorphism invariant action (2.3) in a "top-down" approach [8]. This automatically implements the Stückelberg mechanism for the actions obtained imposing the background Φ A = δ A µ x µ . Indeed, assuming that the condition det(∂ µ Φ A ) = 0 holds, it is possible to find a local set of coordinates (unitary gauge) such that Notice that, defined in this way, the unitary gauge is actually a collection of gauges, which are equivalent modulo constant shifts, i.e. Φ A = x A + c A with ∂c A /∂x µ = 0, where for simplicity we do not distinguish between spacetime and internal manifold indices. These gauges are all equivalent since the action is, by assumption, shift invariant in the Stückelberg fields.
The unitary gauge reflects the fact that since we have four scalars Φ A and four spacetime coordinates, we can use the diffeomorphism invariance of the action (5.1) to absorb the dynamics of the scalar fields into the degrees of freedom of the metric. By using the ADM splitting of spacetime to describe the metric g µν and its perturbations [78], it is then straightforward to make a connection with MG in the so-called broken phase. Therefore, this phase (of broken diffeomorphisms) can be simply understood as a way of expressing the dynamics of the comoving coordinates of the selfgravitating medium through their effect on the metric perturbations. The correspondence between self-gravitating media and massive and modified gravity models written explicitly in terms of metric perturbations is depicted in Table 4, where the arrows indicate how to move from one picture to the other.
In the ADM formalism the metric can be written in terms of the lapse N , the shifts N i and the spatial metric γ ij , namely where γ im γ mj = δ i j . In the unitary gauge (6.1), the matrix C AB defined in (2.2) coincides with inverse of the spacetime metric: for convenience, we decompose C AB in the unitary gauge as follows: Then, we can write the building blocks of LO self gravitating media -see (4.6) and Table 1-in terms of the ADM variables: where γ is the 3×3 matrix of components γ ij . In these expressions the dot (·) represents the standard matrix product and we use the following notation: ξ 2 ≡ ξ i γ ij ξ j and (ξ ⊗ ξ) ij ≡ ξ i ξ j . Therefore, the master function U of the action (5.1) becomes an algebraic function of the ADM variables in the unitary gauge: Special cases -see Tables 2 and 3-can be obtained from the remaining SO(3) s invariant scalar operators of Table 1, which are combinations of X, Y , τ n and y n : Notice also that the four-velocities V µ and u µ become in the unitary gauge: If the the shifts are zero, e.g. as in Minkowski and FLRW spacetimes, u µ and V µ coincide, up to a sign; see also Section 3. The action (6.7) was the starting point in [11][12][13] where a large class of non-derivative 10 massive and modified gravity models -those with spatial SO(3) invariance-was studied by means of a Hamiltonian analysis. Models with 6, 5, 3 and 2 degrees of freedom (DoF) were thus found. The condition of SO(3) spatial invariance turns to be important to avoid ghost instabilities. If this restriction is relaxed, 6 DoF generically propagate and one of them (a scalar mode with respect to spatial rotations) is a ghost around flat space, which corresponds to the infamous Boulware-Deser ghost [63].
For example, a SO(3) invariant master function U of the form propagates precisely 5 DoF [12] (provided that U is not a constant). As a matter of fact, this not the most general U that propagates 5 DoF, see [12]. We leave open a detailed exploration of the structure of the most general case, which is anyway contained in our formalism, for a possible future investigation. The functions U and E of (6.10) are generic and thus there exist infinitely many U with 5 DoF. Therefore, according to the discussion in the previous section ghost-free massive gravity can be identified with a specific kind of supersolid. Moreover, setting E = 0 there are still 5 DoF and the medium is a solid. Thus, the construction of ghost-free massive gravity with five DoF needs at least the three spatial Stückelberg fields Φ a . Cases where only 3 or 2 DoF occur also exist [13]. In particular, exactly 3 DoF are present if U has the structure U = U (w n , X) . (6.11) We recall that the above form for U is protected by the symmetry: Φ a → Φ a + f a (Φ 0 ), see Table 2.
It should be stressed that although 3 DoF are present at the non-perturbative level for (6.11), only two DoF are found expanding around flat space [13]. A proposal for a UV completion, where the LO order part of the action is built out of the scalar operators w 1 and w 2 , was put forward in [79] adding the symmetry Φ 0 → Φ 0 + f (Φ 0 ), see Table 2. In that model the fields Φ a are further coupled with a triplet of purely spatial vector fields, whereas the dynamics Φ 0 is dictated by NLO operators. A particular case of (6.11) is given by the Lagrangian where only 2 non-perturbative DoF are found. In this expression λ is a (dimensionless) cosmological constant. This last case is rather interesting since it gives an example of a model of gravity with 2 DoF that is different from GR.

Cosmological perturbations in FLRW
In this section we discuss a basic aspect of cosmological perturbation theory for self-gravitating media.
In particular, we are interested in determining the number of degrees of freedom that propagate in a FLRW background for the different kinds of self-gravitating media that we have considered. A more detailed account of the cosmology of these media will be given elsewhere [98]. Here we present the Lagrangians for perturbations that are required for such a study and discuss their main features. Let us consider a spatially flat 11 FLRW metric written in the following way: where we denote x 0 = t. The derivatives with respect to t will be indicated by primes. For instance, by definition H = a ′ /a is the Hubble function in these coordinates. With this choice of metric, the following field configuration: is compatible with the equations of motion for the Stückelberg fields derived from the action (5.1), giving a spatially homogeneous and isotropic EMT. In what follows we will work in the unitary gauge, setting φ(t) = t. It is important to stress that in general it is not possible to choose N (t) = a(t) or N (t) = 1 at the same time that φ(t) = t. Therefore, N (t) has to be determined using the (background) equations of motion.
We will now study perturbations around the metric (7.1). Choosing the unitary gauge, the dynamics of perturbations for the whole system (metric and Stückelberg fields) is encoded in the metric. 12 Denoting the metric perturbations by h µν , the master function U of the action (5.1) can be expanded up to second order in h µν as follows: where t µν and the masses m 2 i depend on U and its derivatives. Although the operator b is redundant once the τ n are included, it is convenient to single out its effect, since it plays a special role for some symmetry choices. Hence, we write the LO action (5.1) as Then the masses of (7.3) are given by: (n + 1) a 3−2 n U τn , (7.7) n a −2 n U bτn − n a 3−2 n U τn + n 2 a 3−4 n U τnτn , (7.8) These masses can be simplified (on a case by case basis) taking into account the background equations of motion. Concretely, the condition t µν = 0 is required for consistency. This gives two equations: These equations are equivalent to the conservation of the medium's EMT on the background or, equivalently, to the equations of motion of the Stückelberg fields (in the background), see Section 2.
Eliminating H ′ from the second of them we obtain the condition where n a 3−2 n U X τn . (7.13) Therefore, we have the following possibilities: • F 1 = 0 and therefore it is possible to solve for N ′ . As a result, N is dynamical.
• Both F 1 and F 2 are identically zero. In this case N is not fixed by the background (so there exists a residual gauge freedom). We can fix N , for instance, to be equal to a.
Let us now consider the scalar, vector and tensor modes of the metric with respect to spatial rotations. A general perturbation of (7.1) can be written as: where In the scalar sector, one can integrate out the fields v and ψ of (7.15), arriving to the following canonical form for the effective quadratic Lagrangian L (s) where the 2 × 2 matrices K, D and A mass satisfy D t = −D, K t = K and A t mass = A mass , respectively. The propagation of modes is related to the determinant of the matrix K: if Det(K) = 0 there can be at most one propagating scalar DoF. We have that Defining the effective masses we see that if M 2 0 = 0 or M 2 1 = 0, there is a single propagating scalar mode; whereas if M 2 1 = M 2 0 = 0 no scalar mode propagates. In a general case, two different scalar modes can propagate. The study of scalar perturbations is important for cosmological signatures that could help to distinguish a modification of gravity such as the ones discussed here from a pure cosmological constant, see [80] for the prospects with the Euclid satellite. First, given a certain background evolution, the growth of dark matter perturbations is generically altered if there are extra scalar modes. This can be detected by measuring precisely the so-called growth index and its scale and redshift dependence, see e.g. [81][82][83][84][85]. Moreover, a detection of the speed of propagation of DE perturbations would also indicate a deviation from the ΛCDM model, see e.g. [86][87][88][89][90][91]. Besides, a late-time measurement of a non-negligible excess in the gravitational slip (i.e. the difference between the two, gauge invariant, Bardeen potentials [94]) would signal the presence of an anisotropic stress component, see e.g. [83,92,93]. In our case, a gravitational slip appears in solids and supersolids, for which (at linear order in perturbations) where π L is defined via π i L = ∂ i π L and π i = π i L + π i T , being ǫ ijk ∂ j π k L = 0, see (2.8), and we define: n 2 a −2(n−1) U τn . (7.21) Notice also that in a general case, the two scalar DoF that appear from the action (7.4) allow for the propagation of an adiabatic mode and a entropy mode, see also [51]. In the tensor sector there are always two propagating DoF, with quadratic Lagrangian A non-zero M 2 2 changes the propagation speed of gravitational waves with respect to that of light. This could induce, for instance, observable effects on the propagation and lensing of CMB B-modes if the continuous medium is relevant at sufficiently early times, see [95,96]. It is remarkable that the  same quantity, M 2 2 , that is responsible for the gravitational slip also determines the propagation of gravitational waves. The relation between these two effects has also been noted for V Diff systems of three scalars at NLO in derivatives [15] and for other examples of modified gravity models [97].
Finally, in the vector sector the quadratic Lagrangian reads The fields u i have a purely algebraic equation of motion and thus can be integrated out, giving the LagrangianL The vectors s i propagate only if M 2 1 = 0, see Table 5. The dispersion relation is not trivial only when M 2 2 = 0. Thus, M 2 2 , besides controlling the dispersion relation of tensors, also determines the dynamics of the vectors. The fact that M 2 2 = 0 in the unitary gauge is equivalent to the presence of an anisotropic stress in the EMT of the Stückelberg fields. As we mentioned above, a complete study of the cosmology of these media will be given in a separate publication [98]. We summarize in Table  5 the DoF counting for various examples of sets of operators that may appear inside the LO master function.

Dark matter and dark energy signatures
To illustrate the novel phenomenology that arises for the description of dark energy and dark matter from the effective theory of self-gravitating media, we will focus now in certain aspects of cosmological perturbations that are specific of our framework. In particular, we will consider scalar perturbations for non-barotropic fluids and, also, the propagation of gravitational waves.

Scalar perturbations in non-barotropic fluids
As we mentioned in Section 5.2.2, if the symmetry Φ 0 → Φ 0 + f (Φ a ) is imposed on the general superfluid master function U (b, X, Y ), the operator X drops and the resulting U (b, Y ) system is a perfect fluid. This type of perfect fluid is generically non-barotropic because the pressure and the energy density cannot be written one as a function of the other, see (5.20). This kind of action can be used to describe a generic perfect fluid, since b can be identified with the entropy density and Y with the chemical potential. 13 Then, the entropy per particle σ = b/U Y is non-linearly conserved in time [51]. The pressure perturbation can be written as where φ = φ(t) describes the background evolution of the field Φ 0 , see (7.2). In this expression and in all that follows, primes denote derivatives with respect conformal time t. The quantities c 2 s and c 2 b are, respectively, the standard adiabatic sound speed (i.e. the variation of the pressure with respect to the density at fixed σ) and the variation of the pressure with respect to the density at constant b. The latter of the two determines the evolution of φ, which obeys (from the conservation of the EMT): Let us now write the scalar part of the perturbed metric in the conformal Newtonian gauge: We will work in Fourier space, so that the perturbations ψ and τ will be functions of conformal time and the modulus, k, of the comoving momentum k i . Being U (b, Y ) a perfect fluid, there is no anisotropic stress and, in consequence, τ = ψ. The evolution of the metric perturbation is then governed by the equation We recall that the quantity Γ ≡ δp − c 2 s δρ (that is sometimes called "intrinsic entropy perturbation") is gauge invariant and ψ corresponds to one of the (gauge invariant) Bardeen potentials [94]. The equation (8.4) is a linear combination of the 0-0 and the i-i Einstein equations. In the case at hand, it contains a generic time-dependent speed of sound c 2 s , plus a time-dependent source term that is due to the non-barotropic character of U (b, Y ).
Indeed, whereas δσ ≡ δσ(k) is constant in time, with its k-dependence fixed by initial conditions, the factor that multiplies it in equation (8.4) has a non-trivial time dependence determined by c 2 s and c 2 b ; which in turn depend on the function U (b, Y ). Dynamical (and thermodynamic) stability of the perturbations requires c 2 s > 0; see [51]. In the special case of non-zero and constant c 2 b , from (8.2) one has that where φ ′ 0 is a constant and we take a = a 0 = 1 today. Then, the solution of (8.4) during a matter domination period, where H ∝ 1/ √ a, is given by whereψ =ψ(k) is the constant solution of the homogeneous equation and Ω 0 m and H 0 = H 0 are, respectively, the relative matter density and the Hubble parameter in cosmic time, both evaluated today. We have neglected the decreasing solution of (8.4) and the effect of c 2 s , which is assumed to be positive but sufficiently small to be negligible for dark matter.
The energy density perturbation can be written directly from the 0-0 Einstein equation: We get from it the time evolution of the density contrast δ = δρ/ρ : is the non-entropic (or adiabatic) part of the fluctuation and we define the constant With respect to the standard cold dark matter (for which δσ = 0), the entropic term modifies the time evolution of the density contrast. It is interesting to note that as soon as c 2 b < 0, there is an entropic mode that grows with respect to the adiabatic one, both for super-and sub-horizon scales. 14 As an explicit example we can take the function U (b, Y ) corresponding to a perfect non relativistic gas [51]: where g s is the number of spin states and m is the atomic mass of the gas. One gets c 2 b = 2/3 so that φ ′ ∝ a −1 , see (8.5), and Therefore, in this example the entropic corrections decrease in time with respect the adiabatic perturbation (since c 2 b > 0). A more dramatic example of the influence of δσ is obtained, for instance, by taking Notice that in a more general case, where c 2 s is not negligible, there is an interplay between c 2 s and c 2 b that makes the growth condition more complicated that just c 2 b < 1.
In this case, c 2 b = −1 and φ ′ ∝ a 4 , see (8.5), so that b Y = φ ′ 0 is constant in time, and (8.14) In the regime in which λ/a 3 dominates the denominator, w → 0 (as cold dark matter) in a first approximation. In the opposite limit, w → −1 and the fluid describes a phase of dark energy domination. During matter domination, when w ∼ 0, the energy density contrast is given by (8.15) In this case, the entropic perturbation grows faster than the adiabatic one. Although in this example, in order to keep under control the validity of linear perturbation theory, very strong constraints have to be imposed on the initial conditions (specifically onσ), it illustrates the relevance of possible interactions between dark matter and dark energy in the context of the effective theory of selfgravitating media. In general, by splitting one can systematically study the coupling of dark matter and dark energy, being the corresponding EMT tensors not separately conserved.
To conclude this subsection it is also worth hinting how the evolution of perturbations is affected for a superfluid and a solid. In the first of these two cases, described by a function U (b, X, Y ), the function δσ is not constant in time and the conservation equation δσ ′ = 0 will be replaced by a dynamical one. On the other hand, if a solid is considered, an anisotropic stress is present and the two Bardeen potentials are not equal, so that where π L is the longitudinal mode of the spatial Stückelberg perturbations: Φ i = x i + ∂ i π L , see also (7.20). We plan to do a complete analysis of both effects in a forthcoming paper.

Gravitational waves
The dynamics of gravitational waves is governed by (7.22). Without loss of generality, we can choose N (t) = a(t), so that in Fourier space we have the action The resulting propagation equation is The difference with respect to the standard expression is the presence of the time-dependent effective mass squared, M 2 2 , which is non-vanishing for solids and supersolids. As we had anticipated, this mass also controls the difference of the Bardeen potentials, see (8.16). We point out that, at leading order in derivatives (LO), M 2 2 induces a modification of the speed of propagation (the dispersion relation) of gravitational waves, which becomes scale-dependent. Interestingly, it was shown in [15] that at next-to-leading order (NLO), the propagation speed changes in a scale-independent (but time-dependent) way in the case of V Diff invariant media made of just the three spatial Stückelberg fields. The same kind of modification is expected to occur at NLO for a generic medium, leading to a generic propagation equation of the form:

Conclusions and outlook
The main result of this work is that a vast class of modified gravity models, including massive gravity (MG), can be interpreted as self-gravitating continuous media. The low-energy dynamics of these systems is described by the EFT of four scalar fields Φ A that respect shift symmetries. Due to these symmetries, the fields are derivatively coupled between them and are minimally coupled to gravity at leading order in derivatives. At this order, the action is a functional of ten independent scalars encoded in the induced four-dimensional (inverse) metric The four scalar fields can be interpreted as the (comoving) coordinates of the medium. They can be conveniently split into three spatial coordinates Φ a , a = 1, 2, 3, and a temporal coordinate Φ 0 . Using diffeomorphism (diff) invariance, the EFT can be examined in the unitary gauge, where the scalar fields are "frozen" to be coincident with the spacetime coordinates. With this choice of spacetime coordinates, the induced metric is C AB = g AB . This sends all the dynamics of the medium into the gravitational metric and allows to make direct contact with the traditional framework of MG. The inverse path (going from MG to continuous media) can be travelled by using the well-known Stückelberg "trick", which allows to write in a diff invariant way any theory in which diffs appear to be broken. Indeed, the four scalar fields Φ A of our EFT can also be seen as the four Stückelberg fields of a diff invariant formulation of MG. We have discussed in Section 6 how to move from one picture to the other.
We have provided a comprehensive and systematic classification of massive gravity and modified gravity models in terms of symmetries and low-energy degrees of freedom, establishing a direct connection to the theory of (self-gravitating) continuous media and the pull-back formalism. In doing so, we have determined also which are the propagating degrees of freedom for each case (as a function of the symmetries) and obtained the Lagrangians that are necessary for the study of (linear) cosmological perturbations at large scales, highlighting which operators affect the different kinds of observables.
The mechanical and thermodynamic properties of these continuous media (or modified/massive gravity) models depend crucially on extra symmetries that can be imposed on the scalar sector and may restrict drastically the number of allowed operators. The minimal assumption we have used in this work is that the action is invariant under (spatial) SO(3) rotations of the fields Φ a , a = 1, 2, 3. This reduces the number of allowed operators from ten to nine at leading order in derivatives; see Table 1. Imposing further assumptions, such as symmetries that relate the spatial and temporal Stückelberg sectors -see Table 2-, the number of independent operators can be decreased even down to a single one. For instance, volume preserving internal spatial diffs selects just three operators and the resulting action describes superfluids.
Another way to simplify the action is to make it depend (by assumption) only on either the spatial or temporal Stückelberg fields; see Table 3. The models that are allowed in those cases can describe self-gravitating solids, perfect fluids and superfluids. In the general SO(3) spatially symmetric case, which contains the four scalars Φ A , the medium shares some of its properties with both superfluids and solids, and therefore is called a supersolid. Therefore, a general MG model is interpreted a supersolid propagating in spacetime.
Working in the unitary gauge and expanding a generic example of our EFT around Minkowski spacetime, it is straightforward to see that the resulting action for the metric fluctuations will not respect Lorentz symmetry. This is simply a natural consequence of the presence of the medium and the fact that we have assumed the operator content to be dictated by a broad symmetry: internal spatial SO(3) rotations. Nevertheless, it is possible to choose a combination of operators such that the resulting action for metric perturbations is Lorentz invariant, if one wishes to do so [13].
The correspondence between self-gravitating media and massive/modified gravity is intriguing and deserves further study. In the context of cosmology, it can be relevant for a deeper understanding of the properties of dark matter and dark energy. Moreover, the systematic EFT framework presented in this work can aid to devise a program for the interpretation of the data coming from future probes. From this perspective, the key idea is that the symmetry properties of a (coarse-grained) cosmological medium could be imprinted in the the large scale structure of the Universe, including the cosmic microwave background.
We plan to analyze the cosmology of these models in a future publication, by studying in further detail the signatures that would allow to distinguish these models from the standard ΛCDM model [98]. In the present work, we have initiated this program (in Sections 7 and 8) by determining the number of propagating degrees of freedom in FLRW for various representative examples of continuous media (or models of massive/modified gravity), as summarized in Table 5, and by studying analytically some representative examples. In particular, in Section 8 we have considered the modification to the growth function of matter in non-barotropic perfect fluids and the modifications to the propagation of gravitational waves. Other avenues, that are also worthwhile exploring, can be opened if different internal symmetries -such as supersymmetry [26] or Galileon symmetry [99]-are assumed instead of spatial SO(3).
in two different types. First, we have those currents which, by virtue of Noether's theorem, are due to the symmetries that define the Lagrangians of the different models. These currents are grouped in sets of infinite dimension because the defining symmetries are (infinite) subgroups of the internal diffeomorphisms Φ A → Ψ A (Φ B ). Then, we also have currents that are conserved irrespectively of the equations of motion. An example of the latter is where u µ was defined in (3.4). This current is covariantly conserved off-shell and exists for all the models we consider (regardless of their symmetries) that involve the spatial Stückelberg fields Φ a . It has often been interpreted as the entropy current for fluid actions of the form S = d 4 x U (b), see e.g. [10,23], although this is not the only possible interpretation, see e.g. [100].
Independently of whether a current is conserved due to a symmetry of the action or for another reason, an associated charge that is conserved in time can be defined. In general, given a current J µ that satisfies the corresponding time-conserved charge is which can be proven using the divergence theorem and assuming that the fields fall of quickly enough at infinity. For instance, in the particular case of the current (A.1), the charge is and, clearly,V 3 = 0 due to the antisymmetric character of ǫ µαβγ . The physical interpretation of this result is that the flux lines of the medium are neither created nor destroyed, which is a topological statment. In other words, the conservation of V 3 means that these continuous media do not have flux sources nor sinks.

A.1 General volume currents
Let us first consider conserved currents that are independent of the internal symmetries. We just encountered one example of them in (A.1). This type of current can be generalized in a simple way with a permutation of the fields Φ A . Concretely, the four currents are covariantly conserved, and J µ 0 is precisely (A.1). In fact, any four-vector of the form where h is an arbitrary function of the fields Φ A , satisfies This allows to construct further conserved currents by simply multiplying any of the currents J µ D with a function h that does not depend on the field Φ D . The relation (A.7) can be easily proven using the identity √ −g ∇ µ U µ = ∂ µ ( √ −g U µ ), valid for any vector U µ , and the antisymmetry of J µ D under the permutation of two fields Φ A .

A.2 Noether currents
Let us now consider the currents coming from the different symmetries that we have considered in the paper. First of all, we have SO(3) s and the internal shifts (2.1). The first one of these two symmetries leads to the conservation of angular momentum in the internal spatial manifold. The conservation of the currents associated to the shift symmetries is nothing but the equations of motion of the Stückelber fields, see (2.5). We will now enumerate the conserved currents due to the additional symmetries of Section 4.1.

A.2.1 Four-dimensional media
• If the action is invariant under V s Diff, the volume-preserving spatial diffeomorphisms in the space of Stückelberg fields, the operators b, Y and X are selected and the resulting LO medium is a superfluid. The currents associated to this symmetry are known to be related to vorticity conservation, see for instance [9,10]. These currents are of the form: where ε a (Φ j ) satisfies ∂ε a /∂Φ a = 0 . A basis of vorticity charges (in which any other vorticity charge can be expressed) is constructed choosing ε a = ǫ abc α b ∂δ 3 (Φ −Φ)/∂Φ c , with constant α b , see e.g. [10].
• The covariantly conserved Noether currents associated to the symmetry Φ a → Φ a + f a (Φ 0 ), which selects the operators X and w n , are: n U wn W n−1 · j µ a , (A.9) where we define Each set of three functions {f 1 , f 2 , f 3 } of the temporal Stückelberg field defines a different current. By using the unitary gauge expressions (6.5) we interpret Φ a → Φ a + f a (Φ 0 ) as a transformation that does not change the slicing of spacetime but modifies the threading with a time-dependent shift: N → N , N i → N i + N −1 ∂f i /∂Φ 0 . Notice the spatial part of the metric, γ ij , is invariant under Φ a → Φ a + f a (Φ 0 ).
• In the case of the symmetry Φ 0 → Φ 0 + f (Φ a ), the symmetric LO scalar operators are Y and τ n , and the covariantly conserved currents are: It seems natural to interpret this set of currents -recall that there is a current for each choice of the function f -as transporting charges of the fluid in the direction of the four-velocity u µ . Notice that these currents are parallel to the current (A.1). Therefore, if (A.1) is interpreted as the entropy current, the currents (A.11) carry charges that flow with the entropy.
• If the symmetry Φ 0 → Φ 0 + f (Φ 0 ) is imposed, the LO action depends on the operators τ n and O αβn . The operators τ n do not contribute to the conserved currents, since they do not contain Φ 0 . The currents in this case are where the sums extend over all the O αβn operators and As we discussed in Section (4.1), the symmetry Φ 0 → Φ 0 + f (Φ 0 ) allows at LO an infinite number of different scalars O αβn . This is because the exponents α and β appearing in the definition of O αβn are arbitrary real numbers (whereas n can take the values 1, 2 or 3). In consequence the summatory α, β, n has to be understood to extend over all the possible values of these parameters, unless some restriction is enforced on them.

A.2.2 Media with reduced internal dimensionality
• We can also consider the case in which the Lagrangian depends only on the temporal Stückelberg field, although, as we explained is Section (5.1.3), we do not know of any symmetry that forbids the spatial Stückelberg fields at LO (and allows them at higher orders). In this case the master function U depends exclusively on X at LO and the relevant symmetry is the invariance under Φ 0 → Φ 0 + c 0 , with ∂ µ c 0 = 0. This simply tells us that the current is covariantly conserved. The conservation of X µ is precisely the equation of motion for Φ 0 , given the master function U (X).
• Analogously, we can consider the case in which only the spatial Stückelberg fields are present. Since we are assuming always an SO(3) s symmetry this corresponds to the standard solid U (τ n ). Notice, again, that we know of no symmetry that prevents the appearance of Φ 0 at LO but reintroduces this field at higher orders. Aside from the conservation of internal angular momentum, due to SO(3) s , there are also the currents associated to the internal translational invariance Φ a → Φ a + c a , with ∂ µ c a = 0. This symmetry generates in this case the set of covariantly conserved currents t µ = c a 3 n=1 n U τn B n−1 ab ∂ µ Φ b , (A. 15) a base of which is obtained by setting, alternatively, c a to 0 or 1 for a = 1, 2 or 3.
The currents associated to SO(3) s invariance have a similar structure. In fact, they can be formally obtained from the previous ones simply replacing c a by r ab Φ b , where r ab belongs to the Lie algebra of SO(3) s , and so is antisymmetric. These are just the generators of the angular momentum in the internal spatial manifold. Therefore, the currents 16) are also conserved in this case. A basis for this currents can be obtained choosing the matrices r ab as the standard L x , L y , L z generators. The conservation of the currents of this basis are the equations of motion for the spatial Stückelberg fields.

B A shortcut to count degrees of freedom
We can determine the number of DoF by studying directly to the equations of motion of the Stückelberg fields (2.5). From the action (2.3) we know that there are at most six DoF, that is: two from the metric plus those coming from the scalar fields (≤ 4). The equation (2.5) can be rewritten in the following form (using the notation U AB = ∂U/∂C AB ): where K µν AB contains only first derivatives of the fields. From this expression we see that the second order time derivatives of Φ A are always proportional to the matrix kernel K 00 AB , thus the number of Stückelberg propagating DoF is ∆ = Rank[K 00 AB ] .

(B.2)
To simplify the analysis it is convenient to use the phonons (2.8) and study the system perturbatively. The coefficient of second order time derivatives of the perturbations are determined byK 00 AB ≡ K 00 AB |Φ A =x A . From (B.1) we have that (K 00 AB + O(∂π))π B = 0 wherē K 00 AB =L AB g 00 + (L AD,CB +L AD,BC ) g 0C g 0D , where the overbars mean that the corresponding quantities are evaluated on Φ A = x A . It is easy to see thatK 00 AB gives the number of DoF also at the non perturbative level. Indeed, by going to the unitary gauge we have that Rank[K 00 AB ] = ∆ and so the counting of DoF is background-independent. For SO(3) s invariant Lagrangians one finds the same results obtained applying the Hamiltonian formalism [11][12][13], that is: − for U = L(g AB ) we get 4 scalar DoF − for U = L(g ab ) + −g 00L ( γ ab ) we get 3 scalar DoF − for L 1 = L(g 00 , γ ab ) we get 1 scalar DoF − for L 0 = −g 00L ( γ ab ) we get 0 scalar DoF The Hamiltonian analysis applied to field theories allows for a non-integer number of DoF (see [11,13]) while the prescription of (B.2) seems to evade such a problems. For example the naive Hamiltonian counting of DoF for a master function of the form U = a C 0a C 0a /C 00 gives 3+1/2 DoF in the Hamiltonian formalism and 3 DoF when computed following (B.2). This mismatch is being studied.