Fractons in effective field theories for spontaneously broken translations

We study the concomitant breaking of spatial translations and dilatations in Ginzburg-Landau-like models, where the dynamics responsible for the symmetry breaking is described by an effective Mexican hat potential for spatial gradients. We show that there are fractonic modes with either subdimensional propagation or no propagation altogether, namely, immobility. Such class of effective field theories encompasses instances of helical superfluids and meta-fluids, where fractons can be connected to an emergent symmetry under higher moment charges, leading in turns to the trivialization of some elastic coefficients. The introduction of a finite charge density alters the mobility properties of fractons and leads to a competition between the chemical potential and the superfluid velocity in determining the gap of the dilaton. The mobility of fractons can also be altered at zero density upon considering additional higher-derivative terms.


Introduction
An interesting aspect of low-energy effective theories is that of emergent symmetries. In the simplest setup of a complex scalar field with a Mexican hat potential, the U (1) symmetry associated to phase rotations of the scalar is spontaneously broken and the lowenergy effective theory is described by a massless Nambu-Goldstone boson. At sufficiently low energies, the effective action of the theory is that of a massless scalar field, which not only enjoys the original U (1) symmetry in the form of a constant shift of the Nambu-Goldstone field, but it is also conformal invariant and has an infinite set of conserved higher-spin currents associated to coordinate-dependent shifts of the Nambu-Goldstone field. Neither the conformal nor the coordinate-dependent shifts are symmetries of the full theory, and they are broken when higher-derivative corrections to the low-energy action are considered. Nevertheless, they can leave an imprint in the properties of the low-energy effective theory.
Similar emergent symmetries at low energies appear in other contexts like low-energy excitations of a Fermi surface, independent spin and spatial rotation symmetries in nonrelativistic theories, etc. Here we want to explore low-energy effective theories with emergent symmetries that lead to (gapless) fractonic modes. Fractons are excitations that are able to move only along a restricted set of spatial directions, or are even completely immobile [1,2]. Gapless fractons appear in a variety of systems such as spin liquids [3][4][5][6][7], dipole-conserving lattice models [8][9][10][11][12] and quantum elasticity [13][14][15][16][17][18][19][20][21][22][23][24]. Hydrodynamics of fractons has been studied in [25][26][27]. Models with spontaneous breaking of symmetries have also been studied [28,29]. At low energies, the models we are going to discuss have similarities to these last, but with the important difference that it is not necessary to impose any exact coordinate-dependent phase rotation or shift symmetry in order to obtain fractonic dispersion relations.
A second aspect that we want to explore is the effect of spontaneous breaking of spacetime symmetries in the counting of Nambu-Goldstone bosons. It is well known that the naïve counting of a number of gapless modes equal to the number of broken generators does not apply in this case [30][31][32][33]. An interesting case is when time translations are broken by a finite chemical potential. Under these conditions some of the Nambu-Goldstone bosons become gapped when the effective unbroken Hamiltonian does not commute with some of the broken generators [34][35][36][37][38][39][40]. In particular, if scale invariance is spontaneously broken together with a global symmetry, the dilaton will get a gap proportional to the chemical potential [40] since the generator of dilatations does not commute with the Hamiltonian. Integrating out the gapped modes and keeping only the gapless modes would be equivalent to applying the inverse Higgs constraints [30,33].
If, instead of time translations, space translations are broken, we expect to find some qualitative similarities. There will be unbroken generators of space translations of the form P i = P i − k ia Q a , where P i are the ordinary generators of space translations and Q a are the generators of spontaneously broken global symmetries. The generator of dilatations D does not commute with the unbroken generators [D, P i ] = iP i , so this might produce a gap for the dilaton dependent on k ia . However, due to the breaking of spatial symmetries, the dispersion relations of the modes can depend in a non-trivial way on the spatial momenta, so the intuition from the chemical potential does not entirely apply to this more complicated situation.
What we will do in this work is to examine these questions using a simple 2 + 1dimensional model which can be viewed as a generalization of the ordinary Mexican hat model for spatial derivative terms of a complex scalar field. Scale invariance is ensured by introducing an additional real scalar acting as a compensating field. It turns out that there is a large space of possible ground states breaking translation invariance, and the effective theory depends crucially on the symmetry realization of the ground state. We restrict to states leading to homogeneous effective theories. We find emergent symmetries leading to fractonic dispersion relations and a strong dependence on spatial momentum that affects both the dispersion relations and the composition of the modes. We also study generalizations to finite chemical potential and to 3 + 1 dimensions for some cases.
The paper is organized as follows. In Section 2 we introduce the model and discuss its ground states and symmetries. In Section 3 we compute the dispersion relations for linearized fluctuations around the ground states and identify the associated Nambu-Goldstone modes. In Section 4 we extend our results to finite density and in Section 5 we introduce a deformation that removes some of the degeneracy of the simpler model and study its effect on the dispersion relations. In Section 6 we try to compare our results with theorems determining the number of gapless Nambu-Goldstone bosons and finally we conclude in Section 7 with a discussion of possible physical systems where similar symmetry realizations and exotic Nambu-Goldstone modes might be found. We have collected several technical results and generalizations to 3 + 1 dimensions in the Appendices.

Translation-breaking Mexican hat model
We consider a (2+1)-dimensional model with two scalar fields, one complex and one real, governed by the following Lagrangian density The "couplings" A, B and H are all real and positive. The real scalar field Ξ presents a standard kinetic term and plays the role of a "compensator field," introduced in order to ensure scale invariance. The detailed scaling dimensions of the couplings and of the fields (and of the expectation values that we will introduce below) are where we considered natural units of energy. The complex scalar field Φ presents instead a non-standard kinetic term. Specifically, given the positivity of A, the quadratic term with spatial gradients has the opposite sign with respect to the standard relativistic action. This is a key ingredient for triggering the breaking of translation symmetry through configurations with non-vanishing gradients. Intuitively, the "wrong" sign in the gradient term for Φ can be thought in analogy to the negative squared mass term of the standard Mexican hat potential. Thus we say that (2.1) features a "gradient Mexican hat" for Φ [41,42].
The equations of motion are given by 3)

Ground states
There is a large class of possible ground states that break spontaneously translation invariance with different patterns, but it is strongly restricted if we demand that the effective action for perturbations around the ground state is homogeneous, leaving just two possible types (see Appendix A). Following the symmetry breaking pattern they exhibit we dub the first "helical superfluid" and the second "meta-fluid". We will discuss both, pointing out the similarities and differences between the two types of ground states.
• Helical superfluid: We consider the following static ansatz for the solutions where the compensator field is spatially constant, while the complex field configuration corresponds to a plane-wave of amplitude ρ and wave-vector k. All the parameters in the ansatz, ρ, k and v, are assumed to be non zero, and without loss of generality also real and positive.
The equations of motion descending from (2.1), when considered upon the ansatz (2.5) and (2.6) reduce to We can rewrite (2.7) and (2.8) as follows: where we have introduced the dimensionless combination which parameterizes the space of non-trivial static solutions. Positivity (and reality) of ξ implies AB > 0 and HB > 0. This is indeed satisfied by our choice of taking A, B and H all positive. Consistency of all the relations in (2.11) requires the following relation on the Lagrangian coefficients necessary to have non-trivial solutions, i.e. v = 0, k = 0 and ρ = 0; notice that this amounts to a fine-tuning. The significance of the fine-tuning becomes apparent when looking at the energy density for a static configuration. For (2.12) it takes the form of a complete square When evaluated on (2.11), the energy density is zero, so these are minimal energy solutions.
It is easy to see that there are two directions of marginal stability; in fact we are fixing only the combination ξ given in (2.11), but the ansatz (2.5) and (2.6) has three independent parameters. In other words, we have a two-dimensional space of ground states for this particular ansatz.
We will expand for small fluctuations around this ground state using the parameterization (2.15) • Meta-fluid: We still consider model (2.1), but with a different background ansatz, namely where b and v are respectively a complex and a real constant. In principle there can be more complicated solutions of this type where one introduces two complex constants b x and b y such that Φ = b x x + b y y. The main difference with the case we study is that (2.16) keeps a combination of spatial and phase rotations of the complex field unbroken, while the more general solution does not. Since we are mainly interested in the breaking of translation symmetry we keep to the isotropic case in order to avoid further complications. where we have used the condition on the coefficients (2.12). This guarantees that the energy density of the configuration vanishes, so these are also minimal energy solutions of the same model. We will perform an expansion of small fluctuations around the background (2.20) The fluctuations u i can be interpreted as displacement fields in a solid, in the spirit of the effective actions proposed in [43][44][45].
Finally, it is worth mentioning that there is not really an unbroken phase, even for A < 0. Indeed the compensator field Ξ appears in the denominator in the interaction term with coefficient B in (2.1), and hence the limit v → 0 is not well-behaved. We henceforth always keep v > 0.

Symmetries and Ward-Takahashi identities
The action defined by the Lagrangian (2.1) presents the following symmetries

22)
• Dilatation symmetry: Note that the U (1) and complex shift symmetries are not independent, we can always use a U (1) transformation to rotate a complex shift into a real one. The set of independent symmetries we discuss will then be dilatations and either the U (1) and real shift or the complex shifts.
In the helical state the U (1) symmetry is broken together with translations along the x direction to a diagonal combination. Real shifts and dilatations are also broken. The symmetry breaking pattern is quite different in the meta-fluid. In this case, it is the complex shift symmetry the one broken with translations, in both x and y directions, to a diagonal combination. A U (1) symmetry that combines the phase change of the complex field and spatial rotations survives, so this phase is rotationally invariant. As in the previous case, dilatation symmetry is also broken.
The naïve counting of Nambu-Goldstone bosons would give us three gapless modes in each case: the Nambu-Goldstone modes associated to U (1), real shift and dilatations in the helical state and the Nambu-Goldstone modes associated to real and imaginary shifts and dilatations in the meta-fluid state. As we will see the naïve counting fails and a mode becomes gapped. We will return to the issue of this counting in Section 6.
In the meta-fluid state the identification of the fluctuations is more or less evident, u i should be associated to spatial translations/complex shifts while τ should correspond to scale transformations. In the helical state χ is clearly related to U (1) rotations/translations in the x direction, but the role of σ and τ is not so obvious. In order to help with the identification of the modes in the following we will consider the Ward-Takahashi identities associated to symmetries. A more detailed derivation of the identities can be found in Appendix B.
The Ward-Takahashi identities at linear order in the fluctuations return different combinations of the linear equations of motion that we will obtain from the Lagrangian in (3.3), (3.4) and (3.5). The extra information we get from the Ward-Takahashi identities is that, when considering the decoupling or high momentum limit (which we will implement by formally taking k → 0, though of course we keep the premise that k = 0 for symmetry breaking to happen), one can establish a connection between the fluctuation fields χ, σ and τ and the U (1), the real shift and the dilatation symmetries. Similarly, for the meta fluid one can identify the dispersion relations that correspond to each mode at high momentum. Note that in the perspective where (2.1) is already an effective theory, the dispersion relations at high momentum would in principle be modified by putative higher derivative terms not included in the Mexican hat model we are studying (scale invariance would be explicitly broken by such corrections). However those would come suppressed by a mass scale that we assume to be much larger than any of the scales in the model so it is still sensible to discuss a high momentum limit.

U (1) symmetry
The U (1) current corresponding to the Lagrangian ∂ µ Φ * ∂ µ Φ has the form Thus, for the model (2.1) we have whose conservation is encoded in the continuity equation 1 Expanding to linear order in the fluctuations of the helical superfluid we have In the k → 0 limit one finds 29) indicating that at large frequency and momentum compared to k, the perturbation χ maps to the Nambu-Goldstone boson of the U (1) symmetry, with a dispersion relation This mode has an unusual dispersion relation, and we will refer to it as a 'lineon' since it moves on a line. We will discuss this in more detail when we introduce the connection to fractons.

Shift symmetry
The (complex) shift current corresponding to the Lagrangian ∂ µ Φ * ∂ µ Φ is given by where the s label stands for "shift". The current is linear in the field because the field variation is a constant. Using (2.31), the current for the model (2.1) can be expressed as follows: The associated continuity equation is which, at linear level in the fluctuations of the helical superfluid, gives two linearly independent equations, (2.28) and In the k → 0 limit, we get Therefore, at large frequencies and momenta compared to k, the perturbation σ can be identified with the Nambu-Goldstone mode of (real) shifts. Again, the unusual dispersion relation ω 2 0 will be discussed later on. For the meta-fluid it is convenient to study only the Ward-Takahashi identity of complex shifts. To linear order in the fluctuations the conservation of the complex shift current produces the equations v ∂ 2 At high momentum, τ is decoupled and the displacements u i combine in two modes with dispersion relations where the trivial mode corresponds to the transverse component ∂ k u k = 0 and the propagating mode to the longitudinal component.

Dilatation symmetry
The Lagrangian (2.1) being scale invariant ensures us that we can improve the energymomentum tensor such that the dilatation conserved current takes the form where V µ is called the virial current. Therefore, the conservation equation is equivalent to say that, on-shell, the trace of the improved energy-momentum tensor is zero up to a total divergence of the virial current The improved energy-momentum tensor contains the following terms By injecting the equations of motion in the trace of (2.42), we have that the virial current is given by We now have an explicit expression for (2.41), which at linear order in the fluctuations of the helical superfluid gives In the k → 0 limit (assuming v, ρ can be kept fixed), one obtains where we have used (2.36). Then, for large values of frequency and momenta, τ can be identified with the Nambu-Goldstone boson for dilatations. In this case the dispersion relation is the usual one for a relativistic massless mode For the meta-fluid, the dilatation Ward-Takahashi identity produces the following equation At high momentum the displacement fields decouple and the dilaton has a relativistic dispersion relation (2.49) as in the helical superfluid.

Connection to fractons
The unusual dispersion relations we have found in (2.30) and (2.36) are not just a peculiarity of the decoupling limit but they are also observed at small frequency and momentum, as we will show in the next sections. A possible way to understand their origin is through emergent symmetries of linearized perturbations around the translation-breaking ground states. These symmetries involve coordinate-dependent shifts of the fields similar to those introduced in fracton models [1,2] and are linked to excitations that are immobile or restricted to move in a subdimensional space. In order to identify the emergent symmetry more easily, we will proceed by studying the quadratic Lagrangian of the perturbations and integrating out the gapped mode. The resulting effective Lagrangian admits a derivative expansion where the symmetry becomes manifest.

Helical superfluid
The action to quadratic order in the fluctuations is In this form, we already observe emergent coordinate-dependent shift symmetries, namely The emergence of these symmetries may explain in part the fractonic behavior observed from the analysis of the Ward-Takahashi identities. The dilaton τ has the symmetry of a massless field, the symmetry under the transformation β can be used to introduce an arbitrary dependence on both x and y in σ while the remaining transformation α allows an arbitrary dependence on y in χ. In this way, the identification of σ as a fracton and χ as a lineon appears naturally. Note that these are not symmetries of the full action, so it is expected that terms of higher order in the fluctuations will not be invariant under them, however this only affects indirectly the dispersion relations by radiative corrections. We can diagonalize the mass terms by performing a rotation of the fields The action becomes The mass of η equals to We can now group terms linear in η inside a squared term (after integrating by parts) and subtract the appropriate η-independent terms (2.57) Next, we integrate out η expanding its solution in derivatives, starting at lowest order with Then, up to the fourth order in derivatives, we get (2.59) Both χ and ϕ are gapless and have constant shift symmetries so there are corresponding conserved charges. Furthermore, up to total derivatives in the Lagrangian (2.59), χ can be shifted by a term depending on the coordinates Symmetry under shifts by linear terms imply that the dipole moment of the charge is conserved, while shifts under quadratic terms imply the conservation of quadrupole and second radial moment. This is characteristic of models of fractons that are immobile. Although higher derivative terms might spoil the shift symmetries, this would only affect the dispersion relations at higher order.
To quadratic order in momentum, the dispersion relations of the gapless fluctuations are

Meta-fluid
To linear order, the spatial derivatives of Φ are Then, expanding the action to quadratic order in the fields, we find We can also write this action in the following way The coefficients C ijkl are the components of the elasticity tensor, that in this case only has a bulk component with K the bulk modulus, which also enters in the coupling between the dilaton and the bulk strain. The mass of the dilaton is (2.66) A large v limit would make the dilaton massless and decoupled from the elastic theory at low energies, this latter remaining otherwise unaffected. Roughly, if there is a big hierarchy between the spontaneous breaking of dilatations and that of translations, one does not expect the low-energy elastic theory to be sensitive to the dilaton physics.
Since the shear modulus vanishes, any deformation with τ = 0, ∂ i u i = 0 has zero energy. Then, the elastic part describes a fluid or a meta-fluid. Note that constant changes in volume can be compensated with a shift of the dilaton, so scale invariance is preserved in this sense. This implies that there is a zero mode associated to the dilatation symmetry and a massive mode which corresponds to the combination squared in (2.63).
We can separate the gapped and gapless modes by doing the shift We can further complete the square Integrating out σ implies solving order by order in derivatives with the leading term To sixth order in derivatives in the action, we are left with (2.71) In this form, we also observe that the shear strain has zero energy and that the action is symmetric under constant changes of the bulk strain. This implies that there are linear and quadratic shift symmetries Then, we have that, not only the charges associated to the constant shifts, but also their dipole and second moments are conserved, this is characteristic of fractonic models. The larger symmetry associated to arbitrary shear and rotational strains corresponds to transverse transformations

Dispersion relations
By a standard pertubation analysis of the model (2.1) around the respective backgrounds of the helical superfluid and the meta-fluid, we will compute the dispersion relations of the fluctuations. This will support and refine some of the results and interpretations which we already derived in the preceding sections.

Helical superfluid
As stated in Subsection 2.1, we perform a fluctuation of the model (2.1) around a planewave background where we consider the parameterization (3. 2) The equations of motion at linear order for the fluctuations are 2 Going to Fourier space, we obtain a homogeneous algebraic system determined by the kinetic matrix: where the first row corresponds to σ, the second one to χ and the third one to τ . In order to have non-trivial solutions for the fluctuations, the determinant for the kinetic matrix should vanish, This leads to a set of conditions for the frequency and momenta that determine the dispersion relations. The fluctuation determinant (3.7) has a ω 2 factor, producing a gapless mode whose dispersion relation is trivial, i.e. identically zero, ω = 0. Apart from such a trivial mode, the spectrum features a gapless and a gapped mode Let us remark that the mass (3.9) agrees with the mass coming from the effective field theory analysis (2.56), m 3 ≡ m η . Proceeding to compute the dispersion relations, we obtain: The expansion at low momenta provides 14) To recover the results obtained from the Ward-Takahashi identities of Section 2.2, and in particular the fractonic behavior, we consider the large-momentum behavior of (3.11).
Of course, the exact trivial mode will remain trivial in any q limit. In order to take the large-momentum limit, we simply take q x , q y m η . Then, we find and Thus, ω 2 and ω 3 swap their roles depending on the sign of (2A−1)q 2 x −q 2 y , which in general depends on the direction in the momentum plane. Note that if A ≤ 1/2 this quantity is always negative, so in that case ω 2 and ω 3 do not change with direction.

Identification of the modes
In order to study the Nambu-Goldstone nature of (3.13), (3.14) and (3.15), we need to determine how they relate to a local spacetime modulation of the various symmetryoriginated zero modes. The study of the Ward-Takahashi identities gave us already a glance into such associations at large momentum. At low momentum, instead, one can get useful information from the effective action, which we have already derived to establish the connection to fractons. Accordingly, we will identify χ with the U (1) Nambu-Goldstone mode and τ and σ with the dilaton and shifton respectively.
iii Gapped mode: ω 3 , mixture of τ and σ orthogonal to the gapless mode. ii Lineon: iii Relativistic mode: These identifications unveil a strong change in the nature of the modes as a function of momentum, this being a reflection of the mixing induced by the breaking of translation symmetry. For A > 1/2, the transmutation does not only occur in the transition from low to large momentum but also depending on the direction in the momentum plane.
We would like now to pause a moment to comment on the relation with the non-scale invariant model of [41], where precisely the helical ground state was considered. The model is basically the same as the present one, where however the fluctuation τ is frozen. The spectrum is easily obtained from the determinant of the upper-left 2-by-2 submatrix of (3.6). It consists of a trivial fractonic mode, and a gapped lineon. Hence we see that the compensator field enforcing scale invariance is a highly non-trivial addition to the model, yielding non-trivial mixing among the modes, and their identification.

Meta-fluid
In order to get the dispersion relations for the meta-fluid, we proceed in a similar fashion as for the helical superfluid. By referring to Section 2.1, we study the fluctuations of the model (2.1) around the background given in (2.16) and (2.17) with the following parametrization: At first order in the fluctuations, the equations of motion are given by The quadratic fluctuation matrix in Fourier space is where the two first lines correspond to u x and u y while the last line is associated to τ . The determinant of this matrix is given by where we have used q 2 = q 2 x + q 2 y . Take notice that the above expression is completely isotropic. The dispersion relations, given by the roots of the determinant, are the following where we have used (2.66). At low momenta, we find the following dispersion relations: We obtain a similar mass spectrum as for the helical superfluid: two massless modes, one being exactly trivial, and a gapped mode. However, we have some qualitative differences in the dispersion relations. Indeed, the non-trivial massless mode has a quadratic dispersion relation while in the helical background it has a linear behavior. This can be traced back to the effective theory for gapless modes (2.71), where there are terms with two time derivatives and four spatial derivatives but there are no terms with just two spatial derivatives. An additional qualitative difference with the helical case is that all the dispersion relations are isotropic in the meta-fluid case. As it was mentioned in Section 2.2, the metafluid background preserves an effective rotation symmetry from the diagonal breaking of U (1) and spatial rotations. This is actually due to the particular ansatz (2.16) for the solution. It is possible to choose a more general solution which will lead to anisotropies in the determinant of the kinetic matrix and hence in the spectrum. The other features of the latter would however be unchanged. Hence, we prefer to deal with the isotropic meta-fluid, for a better clarity of the resulting expressions. On the other hand, in the plane wave background spatial rotations are necessarily broken by the choice of a preferred direction in the solution.
At small momentum, the massive mode (3.29) is associated to the fluctuation τ , as it can be observed in the diagonalization of (3.23) in the q i → 0 limit. We notice that the association to τ matches the effective study of Section 2.3 and that we recover the mass (2.66).
For large momenta, the non-trivial modes have the following dispersion relations so that when A ≤ 1 we have that ω 2 √ A|q| and ω 3 |q| (with |q| = q 2 ), while when A > 1 we have the opposite, ω 2 |q| and ω 3 √ A|q|. By looking at the kinetic matrix in the large q limit, i.e. neglecting all non-leading terms in ω or q, we find that the mode with ω |q| is always aligned with τ , while the mode with ω √ A|q| is aligned with the longitudinal combination of the u i (the other being always the trivial immobile mode). Hence, when A > 1, the modes ω 2 and ω 3 switch nature when going from low to high momenta. On the other hand, note that for the meta-fluid, the trivial mode is always the transverse part of u i , for all momenta.
Again, let us comment briefly on the possibility to have a meta-fluid ground state in the model of [41], where τ is frozen. In this case, we can easily see that there is no scale in the spectrum. Eventually, the spectrum consists of a trivial fractonic mode, and a gapless isotropic mode with linear dispersion relations. Hence we notice that also in this case, the addition of the compensator field changes quite radically the spectrum, due to non-trivial mixing.

Finite density
As it has already been established for internal compact symmetries, working at finite density modifies the spectrum of Nambu-Goldstone modes associated to spontaneous symmetry breaking [36][37][38]40]. We would like to probe how these results extend in our specific spacetime symmetry breaking pattern. To do so, we switch on a chemical potential µ for the U (1) symmetry of the theory (2.1) in the framework of the helical superfluid background. We do not extend the analysis to the meta-fluid because in the presence of a chemical potential the effective action is no longer homogeneous.
The chemical potential introduces a new term in the effective potential ∼ −µ 2 ρ 2 , that makes it unbounded from below and would produce a runaway behavior. Something similar occurs in the model of [40] for the simultaneous breaking of scale invariance and an internal symmetry. In that simpler case the issue was solved by introducing a small deformation of the model that lifts the space of minimal energy states at zero chemical potential and stabilizes it at finite chemical potential. The results at zero density can be recovered by simultaneously sending the chemical potential and the deformation to zero. Following the same reasoning we introduce a new term with coupling λ 2 that preserves the U (1) and dilatation symmetries The additional term breaks explicitly the shift symmetry, and would introduce an explicit dependence on the coordinates in the effective action of the meta-fluid. The equations of motion are given by where the parameters v, ρ, k and µ are all real and non vanishing, and assumed to be positive for simplicity. The equations of motion are We keep the same relation between the coefficients of the action so that the relation (2.11) remains unchanged, but there is an additional condition Therefore ρ is fixed in terms of µ/λ. The zero density limit can be taken keeping ρ fixed if both µ and λ are taken to zero at the same rate. The chemical potential µ is seen as an external parameter that fixes the ensemble. Therefore, ρ, k and v are parameters of the solution that should be solved in terms of A, B and µ. This can alternatively be achieved by minimizing the effective potential: 3 In the present case the ratio v/ρ is fixed by k/µ, more precisely v 6 /ρ 6 ∼ k 2 /µ 2 . Then, if k µ we expect the results to be quite similar to the µ = 0 case with v/ρ 1, in which case the gapless mode would be mostly τ . On the other hand, for µ k they are expected to be closer to the case v/ρ 1, where the mode with a gap proportional to k is mostly τ .

Dispersion relations
We are now ready to perform the fluctuations around our background. The fluctuations are parameterized as follows: (4.12) The linearized equations of motion are 4 Notice that the term 4µ 2 σ in (4.13) and (4.14) spoils the space-modulated shift symmetry we had in the case µ = 0 = λ. We therefore do not expect a trivial mode in the spectrum. In Fourier space, the kinetic matrix for the fluctuations is where, as before, the first line corresponds to σ, the second one to χ and the third one to τ . Its determinant is given by Setting the momenta to zero, one gets whose zeros give the mass spectrum. One thus finds a gapless mode m 2 1 = 0 and two gapped modes, whose gaps are The reduction of the number of massless modes compared to the zero-density case is expected due to the explicit breaking of the shift symmetry by the pair µ and λ. Intuitively, such breaking leads to one less flat direction and hence, to one fewer gapless mode. If we take µ k while keeping ρ and v fixed, one gets (4.20) where we recall that m 2 η as given in (3.9) is of O(k 2 ). Note that since as we already noticed, we have that µ/k ∼ (ρ/v) 3 , the leading term in (4.20) goes to zero really as m 2 2 ∼ µ 8/3 k −2/3 . In any case, the zero density limit returns the spectrum computed in Section 3.1 as expected.
In the opposite limit, k µ, we have Again, note that taking into account the behaviour of ρ/v in the limit, we have that m 2 2 ∼ k 4/3 µ 2/3 , still very much suppressed with respect to m 2 3 ∼ µ 2 . The upshot is thus that in both limits there is a large separation between the larger and smaller gap m 3 m 2 . The columns refer to the case of zero and non-zero chemical potential -to make it more visual, the zero chemical plots are the solid curves while the non-zero chemical ones are dashed. All plots are done with A = 0.125, k = 1.5; the left column is obtained with µ = 0 = λ while the right column is obtained with µ = 1 and λ = 0.5. The VEV value ρ is fixed in the µ = 0 case by the preceding cited parameters but it is not so in the zero chemical potential case. For practicality, we took the same value for ρ in both cases. Since k > µ, it means that v > ρ. Hence, at low momentum, the green curve is mostly shiftonic while the red curve is mostly dilatonic.
It is also possible to compute analytically the dispersion relation for ω 1 , by taking the limit ω, q µ, k in (4.18). Actually, one can solve the resulting equation by further postulating ω q, so that we get: Note that this expression is valid for momenta smaller than the chemical potential. Taking the zero density limit for any fixed momentum one recovers that ω 1 is the trivial mode. At large momentum, in any non-zero q x direction, we can again solve for ω, k, µ q x and get

Identification of the modes
The presence of a chemical potential does not fundamentally alter the equations of motion at very large momentum discussed in Section 3.1. Therefore, we expect that the high momentum identification of the modes remains unchanged at finite density, with the trivial mode ω 1 being mostly σ, the lineon ω 2 or ω 3 being mostly χ (depending on the value of A and the direction in the momentum plane) and the relativistic mode ω 3 or ω 2 being mostly τ . This is confirmed by Figure 1 where we have the same trends at large momentum for µ = 0 and µ = 0, the only difference being the non-zero plateau of the blue curve, which formally is produced by a subleading contribution at high momentum.
Similarly, in the regime k µ we do not expect the identification at low momenta of the modes to be significantly altered. Therefore ω 1 would be mostly χ and, given that v/ρ 1 in this regime, the lower gapped mode ω 2 with gap m 2 ∼ (µ 4 /k) 1/3 would be mostly τ and the higher gapped mode ω 3 with gap m 3 ∼ k would be mostly σ.
Finally, when µ k the most relevant terms producing the mixing of different modes is changed. Taking the matrix (4.16) at zero momentum and µ k leads to In principle both χ and τ become gapless in this limit (they are eigenvectors of M with zero eigenvalue for ω = 0). However, taking into account the O(k 2 ) corrections we see that τ acquires a gap proportional to k while χ remains as the true gapless mode to leading order. Finally, the gapped mode with ω 2 √ 2µ is a linear combination ∼ σ − i √ 2 χ. Summarizing, the identification of the modes is

Removing the degeneracy
The Mexican hat model we have studied in the previous sections has a large emergent symmetry that results in the presence of trivial modes in the spectrum. We can remove partially the emergent symmetry and generate non-trivial dispersion relations for all the modes by introducing additional terms to the action, while at the same time keeping the same symmetry breaking pattern. At fourth order in spatial derivatives and fields, there are two possible extensions 5 However they do not produce qualitatively different results. For simplicity we will set F = 0 in the following. We will then study the extended model where the G-term is the additional part. In general the new term will change the energy of the solutions. If we want to ensure that the helical superfluid background (2.5), (2.6) is a minimal energy solution we have to modify the relation between the coefficients to With this choice the dimensionless combination This breaks partially the symmetry characterized by the transformation β in (2.52). With the new term, β is restricted to be a function at most linear in y, but yet arbitrary in x.
So the emergent symmetry with nonzero G is Following the same reasoning we did previously, we expect χ and σ to be both lineons, with χ moving along the x direction and σ along the y direction.
For the meta-fluid the ansatz (2.19) introduces a term in the action for fluctuations where u ij = ∂ i u j + ∂ j u i − δ ij ∂ k u k is the shear strain, and the second equality is obtained up to total derivatives. Note that, in contrast to the helical superfluid, we do not have to change the relation of H with the other coefficients since this term does not give a contribution to the energy density of the background. Furthermore, the meta-fluid is stable for G < 0, while the helical fluid is stable for G > 0, since this gives the right sign to the kinetic terms in (5.6) and (5.8). So, with the new term, only one of the two states would be realized depending on the values of the coefficients we choose to extend the model. In the meta-fluid the new term introduces a shear modulus G = −4G|b| 2 /v 6 , that removes most of the symmetries in (2.72) and (2.73), leaving just the symmetries for massless fields.
Integrating out the massive dilaton as before will remove the zero-momentum bulk modulus, but the higher derivative terms only affect the dispersion relation at higher order in momentum. Then, the effective low-energy theory is almost the same as ordinary elasticity, the dispersion relation for the fluctuations u i is at lowest order in momentum (5.9)

Ward-Takahashi identities Helical superfluid
For the helical superfluid the U (1) Ward-Takahashi identity (2.28) does not change when G is introduced, so the U (1) Nambu-Goldstone mode has the same dispersion relation at high momentum. The real shift symmetry Ward-Takahashi identity becomes The dilatation Ward-Takahashi identity also acquires a new contribution Since ξ is a fixed quantity, in the high momentum limit k → 0, the dispersion relation of the shifton is modified to On the other hand the dilaton keeps a relativistic dispersion relation in this limit. This confirms our analysis of the emergent symmetries where we predicted that σ would behave as a lineon moving along the y direction.

Meta-fluid
When G is introduced, the dilatation Ward-Takahashi identity for the meta-fluid does not change, but there is a new term in the complex shift Ward-Takahashi identity At high momentum this gives two modes with dispersion relations with the first mode corresponding to the transverse and the second to the longitudinal components of u i .

Dispersion relations for the helical superfluid
We obtain the following equations of motion at linear order for the fluctuations: Sending the parameter G to zero (keeping the parameter k, ρ and v fixed) permits to recover the vacuum as well as the equations of motion of the G = 0 model (2.1). Hence, in this specific limit, we expect to recover smoothly the original spectrum. Going to Fourier space, we obtain a homogeneous algebraic system for the equations of motion driven by the kinetic matrix: The determinant evaluates to where det(M ) G=0 is given by (3.7).
If we specifically look for a trivial root of the determinant, we do not find one: This immediately tells us that there is no longer a trivial mode. This is consistent with the analysis we made based on the emergent shift symmetries. The spectrum features two gapless modes and one gapped mode Notice that the difference with the case G = 0 is hidden in the relation among the parameters, where the correction is given by a factor (B − G) instead of simply B. So sending G to zero smoothly provides the masses of the G = 0 case. Proceeding to compute the dispersion relations at small momenta, we obtain: We recover smoothly the G = 0 case in the limit of zero G. Note that the expression for the mode ω 1 has an unusual non-analytic dependence with momentum but overall it goes like ω 1 ∼ q 2 , while the other gapless mode is linear ω 2 ∼ q. The dispersion relation of ω 1 is confirmed by the numerical study of Figure 3. In fact, the plots display the trivialization of the dispersion relations for q x = 0 and q y = 0, and the non trivial bump in the quadrant in between. Also, the analytic expression predicts the changes in slope we observed in the 3D plot. Indeed, according to (5.24), the starting slope of the dispersion relation at fixed q x > 0 is larger than the starting slope at fixed q y > 0. Hence, in order for the dispersion relations to join continuously, the fixed q x > 0 dispersion relation should bend downwards.
A final comment is that at large momentum we observe that (5.18) diagonalizes. In particular, the modes ω 1 and ω 2 are respectively transverse and longitudinal lineons with the dispersion relations 2A q x when q x , q y k . (5.27) This is in agreement with the analysis of Ward-Takahashi identities we did previously.

Identification of the modes
Since the new term proportional to G adds a contribution to the kinetic matrix proportional to q 2 y , the separation between gapless and gapped modes at low momentum is the same as for G = 0. Indeed, the gap in (3.15)  gapless mode (3.14) are the same with nonzero G (5.24). The nature of the high momentum modes is easily identified with the help of the Ward-Takahashi identities. So the identification of the modes is essentially the same as in section 3.2, except ω 1 has a non-analytic behavior at low momentum and becomes a lineon propagating in the y direction at high momentum.

Dispersion relations for the meta-fluid
Let us consider directly the kinetic matrix where we recall that G = −4G|b| 2 /v 6 > 0. We notice that only the first two diagonal terms are modified compared to the G = 0 case. Therefore, we expect that only two of the three dispersion relations will be more significantly affected by the correction, namely the lightest modes. The determinant of the kinetic matrix reads as follows It rightly reduces to (3.24) when G = 0. From this expression one can immediately see that what was formerly the immobile fracton, acquires isotropic and linear dispersion relations which are valid for any momenta, and are entirely controlled by G. One can further find the exact analytical expression for the other two modes, which will depend non-trivially both on G and m τ . At low-momentum ω, q m τ , one can see that the condition det(M ) = 0 gets an additional factor of (ω 2 − Gq 2 ), giving the two gapless modes expected from the low-energy effective theory.
In more detail, at low momentum we have the expansions We recover the expected results from the effective analysis as well as the idea that two of the three modes are more substantially affected by G. Looking at (5.30), (5.31) and (5.32) we get back the dispersion relation of the original model when we send G to zero. At large momentum we can drop the last term in (5.29), so that the determinant completely factorizes and the modes will behave as When G → 0, we recover the original large momentum behavior. According to our previous analysis, at low momentum we can identify ω 1 with the transverse component of the displacement u i , while ω 2 corresponds to the longitudinal part. The gapped mode is mostly the dilaton. At high momentum the Ward-Takahashi identities keep the same identification for the modes as for low momentum.

Extended model at finite density
We generalize the helical superfluid solutions to finite density, as already done in Section 4 for G = 0. In order to stabilize the ground states, we add a shift symmetry breaking λ-term in the Lagrangian, Setting G to zero, we recover the background solution of the finite density G = 0 model. The linearized equations of motion are The kinetic matrix associated to the equations of motion is (5.39) Since G only contributes by terms proportional to the momentum, there are no qualitative  Figure 1). The array of plots is such that each lines corresponds respectively to the longitudinal direction (q y = 0) and the transverse direction (q x = 0). The columns refer to the case of zero and non-zero chemical potential -to make it more visual, the zero chemical plots are the solid curves while the non-zero chemical ones are dashed. All plots are done with A = 0.125, G = 0.25, k = 1.5; the left column is obtained with µ = 0 = λ while the right column is obtained with µ = 1 and λ = 1. The VEV value ρ is fixed in the µ = 0 case by the preceding cited parameters but it is not so in the zero chemical potential case. We took the same value for ρ in both cases for practical reasons. We have that v > ρ, hence, at low momentum, the green curve is mostly shiftonic while the red curve is mostly dilatonic. Notice that since µ < k, the identification of the modes for the finite density case matches the one with zero chemical potential.
differences in the gaps, it is enough to replace B → B − G in the expressions found in Section 4. The high momentum behavior will once more be the same as for zero density. For low and intermediate momenta, we resort to numerics, our results are plotted in Figure 4 and in Figure 5. Comparing with Figure 3, we observe that the ω 1 mode is lifted at q y = 0 when µ = 0, as also happened at G = 0. On the other hand, comparing Figure  2 and Figure 5, the effect of G is to introduce a change in the slope of the dispersion relation in the q x direction. The identification of the modes will be the same as that made at G = 0 in Section 4.

Counting the Nambu-Goldstone modes
The counting of Nambu-Goldstone modes for internal symmetries has been systematically studied through the years, also in relation to their dispersion characteristics [35,[46][47][48][49][50][51][52][53][54], and -in some cases-encompassing translational symmetries too [55] (see [56] for a review). A comparison of the specific results found above with the general knowledge on Nambu-Goldstone counting is interesting because non-trivial. To this purpose, we recapitulate in Table 1 the Nambu-Goldstone modes found explicitly from the study of the helical fluid and meta-fluid fluctuation Lagrangians, as well as their dispersion and analyticity properties. In particular, we stress that we found in general two gapless and a gapped mode.
For internal symmetries, the number of Nambu-Goldstone modes n NG is generically bounded by the number of spontaneously broken symmetries n BS . If there are no terms with single time derivatives in the effective action, however, all the Nambu-Goldstone modes are of type A in the classification of [50,53] and we have n NG = n BS . A priori, that would be the case for the effective actions we found (we set µ = 0 for the moment). However, we are not dealing with internal symmetries only, hence the notion of n BS has to be qualified, as we will do shortly.
An alternative classification, perhaps more pertinent to our situation, is provided by counting theorems which split the total number of Nambu-Goldstone modes n NG according to specific dispersion properties. Defining as type I/type II the modes with an odd/even dispersion relation, respectively, [46] established that n I + 2n II ≥ n BS . (6.1) Turning to broken spacetime symmetries, there are no general counting rules, yet it is known that the number of independent modes can be reduced [30-33, 37, 57, 58]. In essence, if Q a denotes the generators of broken symmetries, P i denotes the unbroken translations and Φ(x) denotes the expectation value of the order parameter, then the following set of identities allows to reduce the number of independent fields where c iab indicates the relevant structure constant of the symmetry algebra. This introduces a constraint such that would-be Nambu-Goldstone bosons appearing on each side of the identity are not independent in general. Let us discuss how the explicit results found in the previous sections relate to the general counting rules. To the counting purposes, we have to consider the following symmetries: • Spacetime: Translations P 1 , P 2 , rotations R and dilatations D; • Internal: U (1) transformations (phase rotations) Q, complex (i.e. real + imaginary) shifts S R , S I .
Let us report the symmetry content of the two kinds of ground states separately.
Helical superfluid: the unbroken symmetries are two translations which leaves in principle five broken symmetries n BS = 5. However, the commutation relation of unbroken translations with the broken generators result in additional conditions This would imply that rotations, dilatations and broken translations are described by a single mode, and there would be a single mode associated to both real and imaginary shifts (indeed, we had already commented earlier on about this). Effectively we would be left with a number of independent broken symmetries n  Table 2, which requires some discussion. For G = 0 in the helical superfluid, there is a trivial gapless mode and a type I mode.  The trivial mode could be counted as type I or II. Still in the helical case, when G = 0, we have almost the same type of modes except that the trivial mode has a non-analytic dispersion relation, but we will still consider it as type II, since at low momenta ω ∼ O(q 2 ). In so doing, the counting rule (6.1) is satisfied if we use n (in) BS as the number of broken symmetries.
Turning to the meta-fluid case, for G = 0 there is an analytic/type II mode and the trivial mode, which could be counted either as type I or II. On the other hand, if G = 0, the two gapless modes are type I/non-analytic. Also in this case we observe that the counting rule (6.1) is satisfied with n (in) BS as the number of broken symmetries. Finally, switching on a chemical potential in the helical superfluid gaps one of the two massless modes. This could be seen as a consequence of introducing a mixing through a single time-derivative term in the effective Lagrangian, so effectively we would be left with a single type B Nambu-Goldstone mode in the classification of [50,53], which turns out to be also type II if we just count the power of the momentum and ignore the breaking of rotational invariance. Note that in going to finite density we introduced an additional coupling that breaks the shift symmetry, so that in this case the remaining gapped mode is a pseudo-Goldstone mode and not a true gapped Namubu-Goldstone. A final comment is the following. A finite chemical potential can be implemented as a linear time dependence in the phase of the charged field. In so doing, however, time-translations would be broken by our choice of ensemble and not, strictly speaking, by a dynamical feature described at the level of the Lagrangian. In contrast, the space modulations that we consider in the helical case, although being formally similar in some aspects, are determined by the gradient Mexican hat potential. Some further comments on this are given in relation to ghost condensation in the next section.

Discussion
The concomitant breaking of dilatations and spatial translations constitutes the main focus of the present study, which adopts the framework of effective field theory. Together with dilatations and translations, an internal Abelian symmetry is broken, too. This symmetry serves two purposes, one physical and another technical. The former consists in modeling a conserved current, thus providing the possibility of considering finite density circumstances; the latter consists in the fact that the breaking of translations and a U (1) symmetry to their diagonal subgroup allows for homogeneous symmetry breaking, in the spirit of Q-lattice models [59].
We have utilized a simple non-relativistic field theoretical setup which allows one to characterize the low-energy modes and derive their effective description in different regimes. The breaking of translations is dynamically induced by a gradient Mexican hat mechanism, namely the competition among a quadratic gradient term driving towards instability then stabilized by higher terms [41]. The gradient Mexican hat, when discretized, connects to lattice models with frustrated interactions [60].
By means of a neat particular example, we clarified the generic fact that the analysis of the modes revolves about three relevant bases: the basis given by the fluctuations of the fields appearing in the Lagrangian; the basis of the fluctuations which diagonalize the lowest-order dynamics at low energy; and the basis associated to the symmetries of the model (this latter is possibly incomplete). Whenever the connection among such bases is non-trivial, we have mixing phenomena. For instance, the Nambu-Goldstone mode associated to a specific symmetry can results from different combinations of the UV or the IR modes, as a function of momenta.
We showed the presence of two degenerate classes of vacua, one associated to a plane wave configuration and possessing a helical structure (i.e. a global phase rotation can be compensated by a suitable translation along the wave-vector of the plane wave), the other associated to complex field configurations which are linear in the coordinates. The latter class admits a specific subclass of isotropic solutions, where a global phase rotation of the background can be compensated by a suitable spatial rotation. We referred to the latter subclass as meta-fluids, because they show a trivial shear elastic response alongside isotropy.
An important feature of the model studied here is the presence of low-energy modes with reduced propagation properties. This fractonic behavior can be associated to enhanced polynomial shift symmetries of the low-energy, linear effective theory and translates into the trivialization of some elastic coefficients. More specifically, we have encountered both completely immobile modes and subdimensional modes, like lineons propagating only along one spatial direction. The immobile fractons can be thought intuitively as plastic deformations which cost zero energy, corresponding to an enlarged vacuum degeneracy, a property which can be compared to the diverging zero-temperature entropy of some fractonic lattice models in their continuum limit [61].
Despite describing an elastic effective field theory with fractonic excitations, the models studied in the present paper differ from the setup where fracton-elasticity duality has been demonstrated [16,17]. There, the two gapless modes of a symmetric gauge field represent the dual encoding of the transverse and longitudinal phonons. Our models lack a gauge field and encode the phonons as Nambu-Goldstone modes, dynamically generated by the interactions. Furthermore, the fractons described in the present paper do not correspond to defects or non-perturbative configurations. They are the low-energy encoding of a trivial (or partially trivial) elastic response.
Immobile fractons cannot move when in isolation but can move due to interactions [16]. Analogously, we expect that the flat fractonic dispersion relations that we encountered are in general "bent" when considering higher non-linear terms in the effective theory, leading to non-trivial propagation.
To conclude, we briefly comment the relation of the present study with some applications in condensed matter and cosmology.
1. LOFF (or FFLO) superconductors: The Larkin-Ovchinnikov-Fulde-Ferrell (LOFF) state [62,63] is characterized by a spatially modulated order parameter and a gapless phonon associated to the breaking of translation invariance, and it might be realized in some unconventional superconductors [64,65]. A similar color superconducting state can also arise at high density in QCD [66,67]. The Mexican hat model is similar to the Ginzburg-Landau functional used to describe superconducting states in the particle-hole symmetric case [68,69]. It could be interesting to revisit the description of the LOFF state and other translation-breaking superconducting states to explore possible emergent symmetries and the fractonic nature of the associated modes. 6 2. Wigner crystals: In the low electron density regime, Wigner crystals can be treated according to classical elastic theory [73,74]. There, the chemical potential µ can possibly be sufficiently low as to allow for a µ-gapped dilaton to become relevant for the low-energy collective-mode description of the crystal response. Despite the non-vanishing µ, in a clean Wigner crystal, all the phonons (either longitudinal or transverse) are gapless. This matches with what we observed in the meta-fluid model, which however suggests that the standard elastic description for Wigner crystals could lack an extra (gapped) dilatonic degree of freedom.
3. Charge density waves: Optical conductivity measurements show a rich structure of peaks [75,76]. In the presence of disorder, the general pattern is characterized by a low-frequency gapped mode corresponding to a pinned collective sliding mode of the density wave condensate. At the opposite end of the spectrum, there is a highfrequency mode associated to the single excitation through the density wave gap.

4.
Ghost condensates: In the search of possible infrared modifications of General Relativity, a mechanism similar to the gradient Mexican hat has been proposed in the time derivative sector, this is usually referred to as ghost condensation [82]. 8 A relativistic generalization of (2.1) is possible, yet it leads to trivial results. Specifically, one can consider the model For λ = 0 the equations of motion imply Φ = 0. Whereas, for λ = 0, the resulting low-energy effective theory features just a relativistic gapless mode, the other two degrees of freedom in (7.1) being associated to an emergent gauge redundancy at the quadratic level in the fluctuations.
Let us first consider the term with the time derivatives: Homogeneity is achieved requiring |f | 2 to be spacetime independent. Hence, f can be considered to have only a space-dependent phase.
Consider now the expansion of the expression squaring to the 'gradient Mexican hat,' to linear order in the fluctuations: The quadratic Lagrangian involves the square of the above expression, and will be homogeneous if and only if each coefficient of the four terms above is itself space-independent (or zero). Taking into account that they come in complex pairs, we have the two conditions relating f and Φ 0 : where a i and b are generic complex space-independent constants. Let us now implement the fact that f must have all its space-dependence in a real phase: From the first of (A.5) we get The second of (A.5) gives now These last two sets of equations imply that ∂ i θ are both constant (and must be real for consistency).
If at least one of the constants is not zero (i.e. b = 0), then we can write θ = k i x i , and we have f ∝ Φ 0 = ρ e ik i x i , i.e. the helical solution (rotated towards a generic direction).
If on the other hand both constants are zero (i.e. b = 0), then θ is a constant that can be reabsorbed in f 0 , the constant value of f . Then Φ 0 is linear, Φ 0 = b i x i , i.e. we have the meta-fluid solution (generalized to the non-isotropic case).

B Ward-Takahashi identities
The model in 2 + 1 dimensions has a Lagrangian density The Noether energy-momentum tensor is One can check that it is conserved on-shell ∂ µ T µ ν = 0. The spatial components are symmetric so the complete stress tensor is The T 00 component is Then, the trace is T µ µ = T 00 + δ ij T ij . (B.7) The traces are All together (B.9) Using the equation of motion for Ξ, we can write this as (B.10) Using now the equations of motion for Φ, Φ *

B.1 Conserved current
The current is The ordinary current is In this model, the components of the conserved current are And the current conservation equation is

B.2 Shift symmetries
The Lagrangian has additional shift symmetry (we consider here only real shifts, to avoid overcounting) Φ → Φ + α, Φ * → Φ * + α. (B.21) The Noether currents associated to this symmetry is If the action only depends on derivatives of Φ, then the current is conserved, since it is a combination of the equations of motion for Φ, Φ * ∂ µ J s µ = 0. (B.23) Note that adding this equation makes the system of equations from the Ward-Takahashi identities equal to the system of equations from the Lagrangian, we have to solve for all the modes. The components are

C Generalizations to 3+1 dimensions
In this appendix, we briefly outline generalizations to 3 + 1 dimensional systems, to show that the essential features of both the helical superfluid and the meta-fluid are unchanged. The only difference is that we have to use different models to generalize the helical superfluid and the meta-fluid, respectively. We will keep the analysis of both models to a minimum, since it turns out that they are very similar to their 2+1 dimensional cousins.

C.1 3+1 dimensional helical superfluid
In order to generalize the helical superfluid, we keep the field content to be a complex scalar Φ and a real scalar Ξ. Only the scaling dimensions of the scalars changes, and hence the compensating powers of Ξ. We thus start with the following Lagrangian, where we have already implemented a condition like (2.12): (C.1) The equations of motion are solved for The expansion is exactly as before Φ = ρ e ikx (1 + σ + iχ) , so that the effective quadratic Lagrangian for the fluctuations about the helical vacuum is which is exactly similar to (2.51) except for a numerical coefficient. The spectrum will then be exactly the same: there is an immobile fracton, a gapless mode which has linear and isotropic dispersion relations at low momentum, but becomes a lineon at high momentum (propagating along x, now one out of three directions), and a gapped mode which has relativistic dispersion relations at high momentum (this is the spectrum for A ≤ 1/2; if A > 1/2 then as before the lineon and the relativistic mode switch roles according to the direction of propagation). The gap is given by

C.2 3+1 dimensional meta-fluid
In order to generalize the meta-fluid, the model has to contain as many (real) scalar fields as space directions, plus the compensator scalar field. Hence we start with 3 real scalar fields Φ i , to which we add Ξ. The Lagrangian is now (C.6) The solution to the equations of motion is We take the fluctuations to be and the quadratic Lagrangian becomes Again, this is very similar to (2.63), up to some numerical coefficients. However, now it involves 4 modes instead of three. But we can immediately see that the only modes that will have non-trivial dispersion relations are the mixtures of τ and the longitudinal component of u i . Then, both transverse modes of u i will be immobile fractons.
As for the non-trivial modes, one is gapped with gap given by and relativistic dispersion relation at high momentum, while the other is gapless with quadratic dispersion relation at low momentum at high momentum. If A > 3/4, the high momentum behavior is switched between the two modes.
To summarize, we see that the generalization to 3+1 dimensions yields physics very similar to the 2+1 dimensional case that we have analyzed in detail, so that we expect that the latter transposes to 3+1 dimensions straightforwardly.