Dipole superfluid hydrodynamics

We construct a theory of hydrodynamic transport for systems with conserved dipole moment, U(1) charge, energy, and momentum. These models have been considered in the context of fractons, since their elementary and isolated charges are immobile by symmetry, and have two known translation-invariant gapless phases: a “p-wave dipole superfluid” phase where the dipole symmetry is spontaneously broken and a “s-wave dipole superfluid” phase where both the U(1) and dipole symmetries are spontaneously broken. We argue on grounds of symmetry and thermodynamics that there is no transitionally-invariant gapless fluid with unbroken dipole symmetry. In this work, we primarily focus on the hydrodynamic description of p-wave dipole superfluids, including leading dissipative corrections. That theory has, in a sense, a dynamical scaling exponent z = 2, and its spectrum of fluctuations includes novel subdiffusive modes ω ∼ −ik4 in the shear sector and magnon-like sound mode ω ∼ ±k2 − ik2. By coupling the fluid to background fields, we find response functions of the various symmetry currents. We also present a preliminary generalization of our work to s-wave dipole superfluids, which resemble z = 1 fluids and feature sound waves and diffusive shear modes, as in an ordinary fluid. However, the spectrum also contains a magnon-like second-sound mode ω ∼ ±k2 ± k4 − ik4 with subdiffusive attenuation.


Introduction
This work weaves together two threads of recent interest, namely fractons and hydrodynamics. Fractons arise in exotic lattice models of interest in both condensed matter and high energy physics [1][2][3][4][5][6]. Generally speaking, these models do not possess a simple continuum effective field theory description at long distances and low energies. For example, they have excitations of restricted mobility, propagating only along sublattices of the microscopic lattice; sometimes have a non-topological, UV-sensitive ground state degeneracy, which, when present, is an example of UV/IR mixing; can have continuum limits that do not commute with the thermodynamic limit; and have exotic symmetries that are sensitive to the details of the underlying lattice.
This last fact usefully organizes many of the unconventional features of these systems [7][8][9]. An instructive prototype is provided by the X-Cube model [7], a soluble theory of quantum mechanical spins on a three-dimensional planar lattice. This model has a subsystem symmetry that flips all the spins along any given plane of the lattice, characterized by a large symmetry group (Z 2 ) Lx+Ly+Lz with L i being the number of lattice sites in the i th direction. The ground states of this model are not invariant under this lattice-sensitive symmetry, leading to an exponentially large ground state degeneracy. 1 Despite how far these systems seem from an effective field theory (EFT) framework, it is sometimes possible to find a near-EFT description; see e.g. [6,10,11]. Focusing on the X-Cube model, one can construct a field theory with exotic charge-2 bosonic matter featuring a U(1) subsystem symmetry and tune the potential so that the preferred ground state spontaneously breaks the U(1) symmetry to a Z 2 subsystem symmetry. Upon coupling this gapless phase to a tensor gauge theory, thereby Higgsing the subsystem symmetry, one finds a massive phase with a description resembling that of a topological quantum field theory (TQFT). This TQFT-like description matches several features of the X-Cube model, including its global symmetry structure and ground state degeneracy.
Another, somewhat simpler, exotic spacetime symmetry considered in the context of fractons is dipole invariance, which, roughly speaking, concerns models where a U(1) charge and the associated dipole moment are conserved; see e.g. [12]. We will give a precise definition in the next section. These models have previously been considered in the literature as a simpler alternative to those with subsystem symmetry, as their elementary and isolated charges are also immobile by symmetry. There are also good reasons to expect that systems with an approximate dipole symmetry are realized in nature. Elastic media in two spatial dimensions [13], vortices in superfluid helium [14], and tilted optical lattices [15] are all expected to have dipole-invariant low-energy descriptions. We will focus on such systems in this work.
Like the subsystem symmetry, dipole symmetry also depends on the details of the underlying lattice. A complex scalar field Ψ(t, x) with U(1) charge q transforms under U(1) dipole transformations along a phase Λ and vector ψ i as Ψ(t, x) → exp(iqΛ − iqψ i x i )Ψ(t, x) [12]. In particular, dipole transformations are valued in the space of single-particle momenta. When the microscopic theory is defined on a lattice in finite volume, the dipole symmetry group is discrete and isomorphic to the dual lattice [16]. In the limit that the lattice spacing vanishes and the volume is infinite, the dipole symmetry becomes R d with d being the number of spatial dimensions.
Pretko has shown how to write down field theories of bosonic matter invariant under this exotic dipole symmetry [12,17]; for similar models on a lattice see [16]. Unlike the case of subsystem symmetries, these models can be made rotationally invariant. In both cases, these bosonic theories are generically strongly correlated as terms with spatial derivatives in the effective action have four or more powers of the fundamental fields. This makes it a challenging task to map out their phase diagram. By tuning the tree-level potential to have a very deep well away from the origin, one might hope to spontaneously break the U(1) symmetry, in which case there is a candidate weakly coupled sigma model description of the phase. However, it is not a priori clear that interactions do not strongly renormalize the potential and restore the symmetry. More generally, one might wonder what other phases exist and with what symmetry breaking patterns.
The best answer to this question to date comes from an analysis of large N versions of models with conserved dipole moment, where one can access the physics at finite interaction strength by studying models with large number of degrees of freedom [16]. Similar techniques may also be applied to study the large N versions of models invariant under subsystem symmetries [18]. One considers models with an N -component charged scalar field Ψ = (Ψ 1 , . . . , Ψ N ), for a large number of fields N ≫ 1, with U(N )-invariant interactions. The U(N ) transformations act on Ψ in the usual way, Ψ → M Ψ, where M ∈ U(N ) is a constant matrix. The diagonal subgroup U(1) ⊂ U(N ), associated with M = exp(iqΛ) for some constant parameter Λ and charge q, is the large N version of the U(1) symmetry. For models with dipole symmetry, we further require invariance under M = exp(−iqψ i x i ). These models reveal three translation-invariant phases of matter with dipole and U(1) symmetries: 1. A gapped dipole-symmetric phase in which both the dipole and U(1) symmetries are preserved.
In the large N models of [16], this phase is found at lattice-scale temperatures or when the U(1) density vanishes.
2. A gapless p-wave dipole superfluid phase in which the dipole symmetry is spontaneously broken but the U(1) symmetry is preserved. This phase may be characterized by a neutral vector condensate iΨ † ↔ ∂ i Ψ . 3. A gapless s-wave dipole superfluid phase in which the U(1) symmetry is spontaneously broken, which in turn leads to the spontaneous breaking of the dipole symmetry as well. 2 This phase is characterized by a charged scalar condensate Ψ . 2 We will have more to say about this shortly, but, for now, note that the discrete dipole symmetry is continuous in infinite volume, leading to a conserved Noether charge, i.e. the dipole moment. That dipole moment D i does not commute with the momenta P i with a commutator [P i , D j ] = iQδ ij , where Q is the U(1) charge [19]. This commutator simply expresses the fact that translating a charge q by a displacement ℓ changes the dipole moment by q ℓ. This commutator implies that if the U(1) symmetry is spontaneously broken but translations are unbroken, then dipole is spontaneously broken as well.
In infinite volume, long-wavelength fluctuations may restore the dipole and U(1) symmetries, leading to a quasi-ordered rather than an ordered phase as in the two-dimensional XY model. At finite temperature (using arguments like that in [20]), this is expected to happen for the p-wave dipole superfluid phase in d = 1, 2 spatial dimensions, and in the s-wave dipole superfluid phase for d ≤ 4. There may also be other phases in which the translation symmetry is spontaneously broken.
The gapped dipole-symmetric phase has immobile charged excitations resembling fractons, albeit with a trivial ground state degeneracy. The n-point functions of charged operators act as an order parameter for spontaneous breaking of dipole symmetry. Consider the two-point function Ψ † (t, x)Ψ(0, 0) , which is invariant under dipole transformations if and only if it vanishes when x i = 0. More generally, an n-point function with charges q m inserted at positions x i m is dipoleinvariant only when it vanishes when n m=1 q m x i m = 0. 3 In this sense, the elementary and isolated charges cannot be transported in a dipole-symmetric phase. By contrast, dipole superfluid phases admit a "dipole condensate," allowing the charges to move in pairs. One can also obtain a gapped phase with immobile excitations by coupling the N = 1 version of the s-wave dipole superfluid phase to a tensor gauge theory [21,22]. Without coupling to a tensor gauge theory, the dipole superfluid phases have a non-trivial ground state degeneracy equal to the number of lattice sites mandated by the spontaneous breaking of the dipole symmetry [16]. In the thermodynamic limit in infinite volume, where the dipole symmetry is continuous, there is a candidate weakly coupled sigma model description of the p-wave dipole superfluid phase in terms of a vector Goldstone φ i [16], shifting under the dipole transformations as φ i → φ i − ψ i . There is also a candidate sigma model description of the s-wave dipole superfluid phase in terms of a scalar Goldstone φ, shifting under the U(1) transformations as φ → φ − Λ. The low-energy behaviours of these two superfluid phases are quite distinct from each other as we shall explore in this work.
In infinite volume and at low energies, many-body systems with conserved dipole moment can be universally described using the framework hydrodynamics, agnostic of the microscopic realizations of the model under consideration. Hydrodynamics is a relatively ancient subject, which has undergone a modern revolution thanks to insights from the AdS/CFT correspondence [23,24] and effective field theory [25][26][27][28][29]. This modern perspective provides us with a rigorous prescription for building a dissipative hydrodynamic theory of fluids with spontaneously broken dipole symmetry. In this work, we will primarily focus on p-wave dipole superfluid hydrodynamics, where the U(1) symmetry is spontaneously unbroken. We will briefly discuss the generalization of our framework to s-wave dipole superfluid hydrodynamics, where both the U(1) and dipole symmetries are spontaneously broken, and will return with full details in an imminent publication [30]. Our main motivation in constructing these hydrodynamic models is to make predictions for future experiments. While there is good reason to expect that models with dipole symmetry can be realized in the lab, a conclusive discovery of these models in nature is yet lacking. Our theory of transport leads to model-independent predictions for the spectrum of fluctuations and the functional form of the response functions of conserved charges. We hope that these results will aid in the future discovery of fractonic matter.
Before we delve into the details of our formalism and results, we note that our work is not the first to tackle the hydrodynamic description for dipole-invariant systems; see [31][32][33][34][35][36]. Nonetheless, there are several distinctive features of our work. We introduce a consistent derivative counting appropriate for dipole superfluids that had not been identified before, allowing us to conclusively enumerate all transport coefficients at leading order in derivatives, as well as the first subleading derivative corrections. We discuss in detail how to couple these hydrodynamic theories to curved spacetime background, allowing us to give a complete account of hydrostatic equilibrium and so the thermodynamics of these phases. In particular, this allows us to transparently infer that for there to be non-trivial charge transport in a gapless transitionally-invariant phase of a system with dipole symmetry, the said dipole symmetry must be spontaneously broken. Relatedly, we can solve the hydrodynamic equations in the presence of background fields and thereby obtain the response functions of various conserved operators. Since we are in a superfluid phase, the system also admits equilibrium states with dissipationless dipole "superflow", characterized by a profile for the dipole Goldstones φ i = x j ξ ji for some constant matrix ξ ij . While we do not focus on these in this work, the framework we setup can be also used to describe transport in such more general states. Finally, while most of the previous work on the subject was focused on p-wave dipole superfluids, our work is the first to identify the qualitative features of s-wave dipole superfluids with conserved momentum. We provide a more detailed comparison to previous work in Appendix D.
The remainder of this manuscript is organized as follows. In Section 2, we give an overview of the rest of the paper, orienting the reader to our main results and approach. We set the stage for our analysis in Section 3, reviewing models with dipole symmetry, coupling to curved spacetime background, spontaneous breaking of dipole symmetry, and subtleties regarding the derivative counting scheme. We go on to construct the hydrostatic effective actions for p-wave dipole superfluids in Section 4 and the full theory of p-wave superfluid hydrodynamics in Section 5. We use this theory to compute the spectrum of linearized fluctuations and response functions in Section 6. We sketch the construction of s-wave superfluid hydrodynamics in Section 7, and wrap up with a Discussion in Section 8. Various technical results and a comparison with previous work are relegated to the Appendix.
Note: While this work was near completion, we were made aware of the forthcoming work of [37] that has significant overlap with our own and is set to appear on the same day.

Overview of hydrodynamics with dipole symmetry
In this Section we provide a summary of the rest of the manuscript. We intend for it to serve as a guide for the technical results that follow.
The basic algorithm towards a hydrodynamic theory we follow in this work is determined almost entirely by symmetries, the symmetry-breaking pattern, and, crucially, on the coupling of conserved currents to slowly varying background fields. This algorithm is quite useful when it comes to constructing the hydrodynamic constitutive relations. It has the added advantage that we can thereby vary the constitutive relations with respect to the background fields and obtain the linear response functions of conserved currents. To avoid getting lost in technical details, here we summarize the overarching formalism and the main results of our construction. We shall focus on p-wave dipole superfluids for concreteness and return to s-wave case towards the end of our summary.
Hydrodynamic fields and conservation laws.-Our starting point is the Ward identities associated with the symmetries under consideration, together with the symmetry-breaking pattern, which in the present context refers to the spontaneous breaking of dipole and/or U(1) symmetry. We then postulate a hydrodynamic description in terms of a set of slowly varying hydrodynamic fields, one for each conserved current, and the Goldstone fields mandated by the symmetry-breaking pattern. Focusing on p-wave dipole superfluids for concreteness, the hydrodynamic fields are a local temperature T , chemical potential µ, velocity u i , and a dipole Goldstone φ i . The relevant symmetry currents are the energy density and flux ǫ t , ǫ i , momentum density π i , symmetric spatial stress tensor τ ij , and U(1) density and flux J t , J i . The dipole symmetry implies a Ward identity, whereby the U(1) flux J i is expressible as the divergence of a symmetric dipole flux J ij through Since the dipole symmetry is discrete in finite volume, one might expect there to be no associated Ward identity. However, as we discuss in Section 3.2, global dipole symmetry can be promoted to a continuous spurionic symmetry, so that dipole transformations imply the Ward identity (2.1) in both finite and infinite volume and, suitably generalized, in curved spacetime as well. The remaining Ward identities are the standard conservation equations for a non-relativistic system with translational and rotational symmetry, i.e.
The curved space versions of these equations are slightly more complicated and will be outlined in Section 3. We should emphasize that systems with dipole-invariance do not feature any boost symmetry, therefore the momentum density π i is not fixed in terms of the U(1) flux J i or the energy flux ǫ i , as it would be in Galilean and Lorentzian theories respectively. In any case, we impose these Ward identities as equations of motion in the hydrodynamic description. Note that there are as many equations as hydrodynamic variables, so that we have a well-defined differential system given that we specify the constitutive relations, i.e. how various conserved currents depend on the hydrodynamic fields. In particular, the dipole Ward identity can be understood as the Josephson condition for the dipole Goldstone.
invariant under global U(1) transformations. Global dipole transformations, however, do act nontrivially on π i and τ ij as One can check that this is a symmetry of the conservation equations (2.2). As we shall explain in our subsequent discussion, ǫ t , ǫ i also transform under dipole transformations on a generic curved spacetime background, while J t , J i , J ij are left invariant; see [38] for more discussion.
Note that due to these transformation properties, it is not possible to construct hydrodynamic constitutive relations for π i and τ ij in the dipole-symmetric phase where the dipole symmetry is not spontaneously broken, provided that J t , J i , and J ij are non-vanishing. This is because in the absence of the dipole Goldstone φ i , all the hydrodynamic variables are dipole-invariant and so there are no low-energy fields available to reproduce the requisite transformation laws of π i and τ ij . This implies that J t , J i , and J ij must vanish in the dipole-symmetric phase, in which case π i and τ ij are dipole-invariant themselves, and hence the low-energy description is effectively that of ordinary neutral hydrodynamics.
When the dipole symmetry is spontaneously broken, we can use the dipole Goldstone φ i to construct the constitutive relations with correct transformation properties, e.g. π i ∼ J t φ i + . . ., where ellipsis characterize further dipole-invariant terms in the constitutive relations. Using suitable redefinitions involving the Goldstone field φ i , we can construct the dipole-invariant versions of various conserved currents and the respective Ward identities, denoted by "tilde" throughout this work. This allows us to formulate a manifestly dipole-invariant version of dipole superfluid dynamics, which can always be undone by performing a local dipole transformation with the parameter ψ i = φ i . For example,π i = π i − J t φ i , etc. This is analogous to constructing manifestly boost-invariant hydrodynamics in the Galilean setting [39].
Since the dipole symmetry is lattice sensitive, our description exhibits UV/IR mixing characteristic of fracton models. However, the UV-sensitivity is relatively modest for dipole symmetries compared to subsystem symmetries. Despite being a low-energy field, the dipole Goldstone field φ i is UV-sensitive, in the sense that it is compact with periodicity ∼ 2π/a, with a being the lattice spacing [16]. This means that n-point functions of φ i are UV-sensitive, but derivatives of φ i are not. Thus, π i and τ ij are UV-sensitive, even in the course-grained hydrodynamic limit, while the other conserved currents are not.
Derivative counting.-An important ingredient in the hydrodynamic description is the derivative counting scheme. Since hydrodynamics describes small fluctuations away from thermal equilibrium, we can organize the respective constitutive relations as a perturbative expansion in derivatives. However, for a consistent truncation of this expansion to any given order in derivatives, we must agree on the relative derivative ordering of various constituent fields. In textbook hydrodynamic treatments [40,41], the local temperature T , chemical potential µ, and velocity u i are treated as O(∂ 0 ) quantities. One then classifies terms in the constitutive relations according to how many space and time derivatives they contain, treating both space derivatives ∂ i and time derivatives ∂ t as O(∂ 1 ). The terms in the constitutive relations with no derivatives are part of "ideal hydrodynamics," while those with a single derivative are the first-order corrections, and so on.
The derivative counting scheme for dipole superfluids is considerably more intricate and is one of the central results of our work. Firstly, the dipole Goldstone field φ i is counted as O(∂ −1 ), as is conventional for a Goldstone field in superfluid hydrodynamics. The temperature T and chemical potential µ are still counted as O(∂ 0 ) quantities. However, the spatial velocity u i is counted as O(∂ 1 ). The underlying physical reason is that the dipole symmetry forbids homogeneous equilibrium states with nonzero charge and nonzero velocity. Therefore, velocity is always a fluctuation in a state of nonzero density. A consistent counting scheme then requires that time derivatives ∂ t scale as O(∂ 2 ), while the spatial derivatives ∂ i are still O(∂ 1 ). In other words, these systems are characterized by a dynamical scaling exponent z = 2.
Due to the disparate counting of time and space derivatives, the density and flux components of conserved currents also have different power countings. The energy density ǫ t , charge density J t , spatial stress τ ij , and dipole flux J ij are O(∂ 0 ); the energy flux ǫ i and charge flux J i are O(∂ 1 ); and the momentum density π i is O(∂ −1 ). The unconventional scaling of π i is due to the fact that it transforms almost like a dipole Goldstone as in (2.3). With these scalings, the hydrodynamic equations are homogeneous in derivatives, and one can classify terms in the constitutive relations by how many derivatives they have relative to the leading scaling.
Note that not all allowed terms in the constitutive relations that can be written down at a particular derivative-order are physical. While the hydrodynamic variables have well-defined physical meaning in thermodynamic equilibrium, they do not have fixed definitions out-of-equilibrium in the hydrodynamic regime. One may use this freedom together with the leading order hydrodynamic equations to eliminate various unphysical terms from the constitutive relations. These steps are analogous to using field redefinitions and equations of motion to remove unphysical coupling constants in a field theory Lagrangian.
Hydrostatic partition function.-By coupling our hydrodynamic theory to a time-independent but slowly spatially-varying background spacetime, the fluid can be placed in hydrostatic equilibrium. These hydrostatic states can be simply described using Euclidean thermal field theory. When the state has finite correlation length, one can dimensionally reduce the theory on the thermal circle and obtain a gapped long wavelength description [42,43]. These states are then described by a local partition function expressed in terms of the time-independent background fields and their derivatives. Varying the hydrostatic partition function with respect to background fields must reproduce the hydrostatic part of the conserved currents under consideration and thus imposes strong consistency constraints on the hydrodynamic constitutive relations. For systems with infinite correlation length, as in a superfluid, the hydrostatic partition function can be constructed using a Euclidean effective action for the Goldstone field coupled to time-independent background fields, which again yields strong hydrostatic constraints on the constitutive relations [44].
Hydrostatic constraints provide us another argument for why transitionally-invariant systems with conserved dipole moment must have their dipole symmetry spontaneously broken for there to be non-trivial charge transport at low-energies. In the dipole-symmetric phase, the correlation length is finite and the hydrostatic equilibrium states ought to be characterized by a local partition function constructed out of time-independent background fields. Let us consider the background sources A t , A i , and a ij coupled to J t , J i , and J ij respectively. If the dipole symmetry was absent, one could have a U(1)-invariant term in the partition function like 1 2 χ(µ 0 + A t ) 2 , with µ 0 the equilibrium chemical potential and χ the charge susceptibility, leading to nonzero equilibrium charge density J t ∼ χ(mu 0 + A t ). However, in the presence of a background frame velocity v i coupled to the momentum density π i , the dipole transformations act on A t as A t → A t − ψ i v i . Therefore, such a term is not compatible with dipole symmetry in the presence of momentum sources and forces the charge density to vanish. More generally, it is an easy exercise to show that there are no terms compatible with dipole symmetry that depend on A t , A i , a ij and transport is effectively neutral. On the other hand, when the dipole symmetry is spontaneously broken, one can write down a dipole-invariant term in the Goldstone effective action as 1 2 , resulting in non-trivial charge density in equilibrium.
Second Law of thermodynamics.-Out-of-equilibrium, the hydrodynamic constitutive relations must satisfy a local version of the Second Law of thermodynamics. That is, we require the existence of an entropy density and flux s t , s i , whose divergence is positive semi-definite everywhere when evaluated on any solution of the hydrodynamic equations. In a homogeneous equilibrium state, these take the form s t = s, s i = 0, with s being the thermodynamic entropy density. Historically, the local Second Law was imposed on hydrodynamic descriptions as a physically motivated constraint, but is now understood to follow from more fundamental considerations in Schwinger-Keldysh effective field theory [45][46][47].
As it turns out, the Second Law requirement exactly reproduces all the "equality-type" constraints arising from hydrostatic consistency mentioned above. In addition, it also results in a set of "inequality-type" constraints, mandating certain dissipative transport coefficients at low derivative orders to be non-negative. In textbook hydrodynamics, these lead to the non-negativity of viscosity and conductivity. Fundamentally, these inequality constraints can be understood as arising from unitarity in a microscopic description viz. the non-negativity of spectral functions. The higher derivative corrections are, roughly speaking, always consistent with the positivity of entropy production provided that the leading corrections are, thanks to an argument by S. Bhattacharyya [48].
Constitutive relations.-Let us briefly look at the constitutive relations allowed by the Second Law of thermodynamics in the absence of various background fields. Since the dipole Goldstone is counted as O(∂ −1 ), there is a spatial tensor ∂ i φ j at zeroth order in derivatives. This tensor is constant in a superflow state with φ i = x j ξ ji , and has a number of independent components that grows rapidly with the dimensionality of space. This leads to a dimension-dependent number of zeroth order scalars that can appear in the thermodynamic equation of state of a dipole superfluid. The complete theory of nonlinear hydrodynamics in general dimension is then quite complicated, analogous in some ways to nonlinear magnetohydrodynamics [49][50][51]. Motivated by the desire to predict future observations, we content ourselves to work at the level of linearized constitutive relations around a zero superflow state where the analysis is more tractable. However, we emphasize that our methods allow for a fully nonlinear hydrodynamic description of dipole superfluids.
Turning off the background fields and working to leading order in derivatives and up to linearized order in fluctuations around a zero superflow state, we find the constitutive relations for the conserved densities where ǫ and q are the thermodynamic energy and charge densities, respectively. The conventional kinetic term in the momentum density, π i ∼ ρ u i where ρ is the kinetic mass density, appears as a leading derivative correction in a p-wave dipole superfluid. The associated fluxes are given as where we have identified the thermodynamic pressure p, dipole pressure p d , and susceptibilities χ m , χ mn , and G d . Furthermore, κ is the dissipative thermal conductivity transport coefficient, which, crucially, affects the constitutive relations at the same order as thermodynamics. This is a qualitatively novel feature of the dipole-invariant hydrodynamic description.
We have also enumerated the most general leading-derivative corrections to the constitutive relations. Accounting for the dependence on background fields, there are 28 hydrostatic coefficients including the kinetic mass density ρ at this order, 12 dissipative coefficients including the traditional shear viscosity η and bulk viscosity ζ, and 7 new non-dissipative non-equilibrium coefficients. More details can be found in Sections 5.2 and 5.3.
Modes and correlators.-We can then plug the constitutive relations above into the hydrodynamic equations and solve for the spectrum of linearized hydrodynamic fluctuations. In the longitudinal sector, the longitudinal velocity u is effectively gapped due to the dipole Ward identity, leaving the longitudinal component of the dipole Goldstone φ , temperature T , and chemical potential µ as the effective low-energy degrees of freedom. These mix to produce a pair of magnonlike modes with dispersion relations ω ∼ ±k 2 − ik 2 + O(k 4 ) and a diffusive mode with dispersion relation ω ∼ −ik 2 + O(k 4 ). The attenuation of the magnon-like modes and the diffusion constant are both controlled only by the thermal conductivity κ, meaning that if the system under consideration is a poor thermal conductor, these modes will become subdiffusive. This feature does not exist in ordinary fluids where other dissipative coefficients, e.g. the bulk viscosity ζ, appear at the same order as κ. Similarly in the transverse sector, the transverse velocity u ⊥ is effectively gapped, while the transverse fluctuations of the dipole Goldstone φ ⊥ give rise to a subdiffusive mode with dispersion relation ω ∼ −ik 4 . The subdiffusion coefficient is fixed in terms of the shear viscosity η of the dipole superfluid.
We also compute the response functions of conserved densities and fluxes using our hydrodynamic theory in Section 6.3. The full suite of hydrodynamic correlation functions is quite complicated, therefore we focus on response functions in the transverse sector for illustrative purposes. We have included a supplementary Mathematica notebook with our arXiv submission that enables users to compute all longitudinal and transverse response functions depending on the particular application in mind. We also present the frequency-dependent optical responses at zero wavevector in Section 6.3. These allow us to find simple Kubo formulae for κ and 8 out of 19 out-of-equilibrium transport coefficients at the order of leading-derivative corrections. The Kubo formulae for the remaining transport coefficients require response functions at nonzero wavevector and frequency and can be obtained using our supplementary Mathematica notebook.
p-wave vs s-wave dipole superfluids.-For concreteness and clarity of presentation, we primarily focus on p-wave dipole superfluids in this work, where the dipole symmetry is spontaneously broken but the U(1) symmetry is preserved. In Section 7, we briefly discuss the s-wave dipole superfluids, where the U(1) symmetry is also spontaneously broken, highlighting some important distinctions with the p-wave case. This phase is characterized by a scalar U(1) Goldstone φ and the vector Goldstone is no longer independent, φ i = −∂ i φ − A i . Following a preliminary analysis, we find evidence that the mode spectrum of this phase has ordinary sound waves ω ∼ ±k − ik 2 + O(k 3 ) in the longitudinal sector and ordinary shear diffusion ω ∼ −ik 2 + O(k 4 ) in the transverse sector. However, as a consequence of dipole symmetry, we do find additional magnon-like second sound waves ω ∼ ±k 2 ± k 4 − ik 4 + O(k 6 ) in the longitudinal sector with subdiffusive attenuation. We leave a rigorous analysis of s-wave dipole superfluids to an imminent publication [30].

Preliminaries
In this Section we review some background we require in order to construct a hydrodynamic theory for systems with conserved dipole moment. We start with a quick introduction to field theories with a dipole symmetry, followed by a review of how to couple such theories to a fixed curved spacetime background. We then discuss the spontaneous breaking of U(1) and dipole symmetry in these models. There we argue that the U(1) breaking leads to dipole breaking as well, and show how to use a dipole Goldstone to reformulate the currents and Ward identities in a dipole-invariant way. We end with an argument for the unconventional derivative counting in the hydrodynamic description of these models.

Dipole symmetry
We start optimistically with the goal of writing down a continuum quantum field theory of a unit-charge scalar field Ψ(t, x) which, in infinite volume, is invariant under spacetime translations, rotations, U(1) phase rotation, and dipole transformations. Under a U(1) phase rotation Λ and dipole transformation ψ i , the scalar field transforms as Pretko has shown how to write down theories of this sort [12]. A representative Lagrangian is The bi-local differential operator D ij involves two derivatives and acts on two fundamental fields of charges q 1 , q 2 as By construction, this operator transforms covariantly under the aforementioned U(1) and dipole ; see [38]. In particular, the usual spatial kinetic term ∂ i Ψ † ∂ i Ψ is forbidden by dipole symmetries, and the simplest terms with spatial derivatives involve at least four powers of the charged fields written above.
The dipole symmetry is simpler to understand in momentum space. Dipole transformations act on a momentum space field Ψ(t, k) of charge q as i.e. as translations in momentum space. Dipole-invariant interactions are those that are, in a sense, invariant under translations in momentum space. For example the interaction D ij (Ψ † , Ψ † )D ij (Ψ, Ψ) can be expressed in momentum space as where the vertex is defined as together with k mn = k m − k n . Dipole invariance here is made manifest as the interaction only depends on the invariant differences of momenta. By going to momentum space, it is almost immediate to write down dipole-invariant theories on a hypercubic spatial lattice with lattice spacing a. The simplest lattice implementation of this interaction takes the same form as above, provided that we now interpret the integral over momenta as a sum over lattice momenta and the interaction is replaced with with the delta function understood as the appropriate one for lattice momenta.
What is the dipole symmetry group? As mentioned in the introduction, dipole transformations are valued in the space of single-particle momenta. This space depends on whether we are in finite or infinite volume, in the continuum or on a lattice. If by finite volume we mean a torus, then the four options are [16,21,22] 1. Infinite volume, continuum: Continuous non-compact dipole symmetry R d .
3. Infinite volume, lattice: Continuous compact dipole symmetry U(1) d . 4. Finite volume, lattice: Discrete compact dipole symmetry Z L 1 × . . . × Z L d , where the lattice has L i sites in the i th direction.
In the lattice cases, the dipole transformations are valued in the Brillouin zone. For a lattice Γ, the dipole group is isomorphic to the dual lattice Γ. Observe that lattice translations are valued in Γ, so that lattice translations and dipole transformations are dual to each other. Equivalently, as stated above, dipole transformations act as translations in momentum space.
In the continuum, there are conservation laws associated with spacetime translations and U(1) transformations given in (2.2), which become Ward identities in the quantum theory. Furthermore, from the action (3.2), one can compute the U(1) and dipole fluxes and check that they satisfy the relation given in (2.1). In the next subsection, we will see that one can define the current J ij by varying the action with respect to a suitable external field, and that this equation may also be understood as a Ward identity associated with a spurionic continuous dipole symmetry.
The Hamiltonian H, total momenta P i , and U(1) charge Q can be constructed from the conserved densities appearing in (2.2) as and, in infinite volume, so that we have rotational symmetry, we have angular momenta M ij = −M ji given by M ij = d d x (π j x i − π i x j ). Importantly, we do not require invariance under any kind of boost symmetry, Lorentzian or Galilean. Thus, we do not have a conserved boost generator. In infinite volume, so that the dipole symmetry is continuous, we also have the dipole moment which is conserved thanks to (2.1) and (2.2), provided that J t falls off suitably fast at infinity. Note that this operator is only well-defined in infinite volume, 4 and that the flux associated with D i is almost but not quite J ij . These charges generate an algebra with nonzero commutators [19] [M ij , (3.10) As a consequence of the symmetry algebra, the momentum operator shifts under dipole transfor- The conserved currents transform as usual under spacetime symmetries and are invariant under U(1) transformations. Interestingly, there is a local version of the transformation law for momentum above, whereby the momentum density π i and spatial stress τ ij are not invariant under dipole symmetry transformations [38]. The transformation laws for a constant dipole transformation in flat space are given in (2.3). As we will see in the next subsection, on a generic curved spacetime background, ǫ t , ǫ i also transform under dipole transformations, while J t , J i , J ij are invariant.

Coupling to curved spacetime background
In this subsection, we discuss how to couple field theories with dipole symmetries to curved spacetime background. This is a necessary ingredient to compute the response functions of conserved currents in the theory. It also makes manifest the transformation properties of these currents under global dipole transformations and will be indispensable for the construction of hydrostatic partition functions for dipole-invariant fluids. The discussion here is a digest of [38]; see also [52]. Some of these matters were considered in [53].

Aristotelian backgrounds
We would like to couple dipole-symmetric field theories to a fixed curved spacetime. Systems of this sort do not respect any boost symmetry and generically couple to an "Aristotelian" spacetime [29,54,55], characterized by a nowhere-vanishing clock-form n µ and a degenerate spatial metric tensor h µν . Together, these background fields couple to the energy current, momentum density, and spatial stress tensor; we will return to the precise coupling structure in the following. In flat space, these sources take the trivial form n µ = δ t µ and h µν = δ i µ δ j ν δ ij . The zero eigenvector of h µν is called the frame velocity v µ and is normalized such that n µ v µ = 1. We can also define an inverse spatial metric h µν such that h µν n µ = 0 and h µρ h ρν = δ µ ν − v µ n ν ≡ h µ ν . Generally, we use h µν to raise indices. We can define a connection on this geometry where the covariant derivative ∇ µ acts on tensors as The covariant derivatives of the spacetime background reads Here £ X denotes the Lie derivative along a vector field X µ . Note that this connection is torsional It is also helpful to define a volume measure γ = det(n µ n ν + h µν ). We can check that (3.14) The curvature tensor associated with this connection is In addition to the spacetime background, we introduce a background U(1) gauge field A µ coupled to the conserved U(1) current. Finally, we need to introduce a symmetric tensor gauge field a µν , satisfying a µν v ν = 0, to couple to the symmetric dipole flux J µν , satisfying J µν n ν = 0 [38,52].
Due to the normalization conditions, the only independent components of the sources are n t , n i , v i , g ij , A t , A i , and a ij . We can express them as (3.16) Here non-covariant indices are lowered and raised with g ij and its inverse g ij respectively.
To define the dipole symmetry covariantly we promote it to a spurionic symmetry that can be imposed locally, in finite or infinite volume, and even in curved space. In flat space it acts by with charged fields Ψ inert. This spurionic symmetry implies the dipole Ward identity J i = ∂ j J ij which must be accounted for in the low-energy hydrodynamic description we obtain later.
With the background sources in place, the curved space version of the scalar theory in (3.2) reads with the derivative operators defined as The definition of D µ contains a Aristotelian covariant derivative when acting on tensor fields. This action, so constructed, is invariant under finite diffeomorphisms, local U(1) transformations, and local dipole transformations, provided that we take the charged scalar to transform as It is also helpful to define an effective dipole gauge field as This field contains information about both the U(1) gauge field A µ and the symmetric tensor gauge field a µν . It is a useful quantity because in the absence of Aristotelian sources, it transforms as a gauge field under dipole transformations, i.e. A λ µ → A λ µ + ∂ µ ψ λ . More generally, its transformation is given as We can also construct the associated field strength This is invariant under dipole transformations in flat space, but transforms non-trivially when coupled to curved Aristotelian spacetime [38].

Conservation equations in curved space
Having identified the appropriate background sources, we can define the conjugate symmetry currents. Given a field theory described by a generating functional W = −i ln Z, we define expectation values of the currents by Here ǫ µ is the energy current, π µ (normalized as π µ n µ = 0) is the momentum density, τ µν (normalized as τ µν n ν = 0 and τ µν = τ νµ ) is the stress tensor, J µ is the charge current, and J µν (normalized as J µν n ν = 0 and J µν = J νµ ) is the dipole flux. In terms of the non-covariant currents appearing in (2.2), these are defined as Note that we have defined J µν using the coupling to the effective dipole gauge field A λ µ . If we had instead coupled to the symmetric tensor gauge field δa µν , various currents would need to be slightly adjusted on account of The benefit of using A λ µ is that the transformation properties of various currents under the dipole shift symmetry is "nicer" [38], i.e. (3.28) We have introduced the notation where we have defined the effective Lorentz force f µ and dipole stress tensor τ µν d as These boil down to the conservation equations in (2.2) in the absence of background sources. After some work, one can verify that the conservation equations are invariant under dipole shifts of the conserved currents given in (3.28).
Were we to interpret A λ µ as an independent field transforming like (3.23), we could introduce an independent operator J [µν] = 0 which would amount to adding intrinsic dipoles. However, in this work, we will always take symmetric dipole gauge field a µν to be fundamental, with A λ µ defined via (3.22). Therefore, we don't have access to the operator J [µν] = 0 in our description.

Spontaneous breaking of dipole symmetry
We would like to obtain a thermodynamic and ultimately hydrodynamic description of many-body systems with dipole symmetry. A consequence of the dipole shift transformations (2.3) is that the homogeneous equilibrium states of such a system must either have no charge or no momentum. A uniform slab of moving charge does not conserve dipole moment and hence is disallowed by dipole symmetry. A more formal argument using the representation theory of the dipole and translation symmetries leads to the same result [16]. Furthermore, as we explained in the Introduction, there is no transport of charge in a dipole-symmetric phase and the hydrodynamic description is effectively that of a neutral fluid.
Therefore, dipole-invariant systems with charge transport must have their dipole symmetry spontaneously broken. 5 Spontaneous breaking of the dipole symmetry gives rise to a UV-sensitive degeneracy of ground states in dipole-invariant systems, as one finds in gapped phases of fracton order [16]. The spontaneous breaking of dipole symmetry is characterized by a vector Goldstone φ i in the low-energy effective theory that transforms nonlinearly under dipole transformations ψ i as a shift, i.e.
More generally, we can also consider systems where the U(1) symmetry is spontaneously broken, giving rise to a scalar Goldstone φ that transforms as φ → φ − Λ under local U(1) transformations. Spontaneously breaking the U(1) symmetry automatically spontaneously breaks the associated dipole symmetry, but the would be dipole Goldstone φ i is generally massive. 6 We will revisit this possibility in Section 7.
The identification of a Goldstone boson implementing the dipole shift symmetry in the hydrodynamics of dipole-invariant systems has also been described in [35,36], where such a Goldstone boson is implicit in that the momentum density π i ∝ J t φ i . This is a necessary consequence of the symmetry algebra (3.10). However, this is only the case when ∂ i φ j = 0, in the absence of external sources, and to leading order in gradients. A non-vanishing expectation value for this operator is analogous to a background superflow in a U(1) superfluid [40,56]. In a forthcoming work [57], we investigate the consequences of such a superflow. Nevertheless, it is worthwhile to consider the simpler and perhaps more physical case of ∂ i φ j = 0 first, which will be the concern of this manuscript. We emphasize that our construction differs somewhat from [35,36] in that we are able to explicitly identify the relationship between φ i and fluctuating background sources and in particular are able to write down all response functions in this limit. Furthermore, we are able to include non-linear interactions among components of the dipole Goldstone. Our formalism also naturally extends to scenarios with nonzero ∂ i φ j .
We can introduce a curved space version of the dipole Goldstone φ µ , normalized such that The generating functional W is given by a path integral over possible configurations of φ µ , i.e.
). The off-shell variations of the effective action S with respect to the background and dynamical fields can be parametrized similar to (3.25) as where X µ = 0, normalized as X µ n µ = 0, is the equation of motion for the dipole Goldstone. We note that the conservation equations for off-shell configurations of φ µ modify (3.29) to Note that the U(1) conservation equation remains the same; the energy and momentum conservation equations also remain the same in form, but the Lorentz force is now modified with The dipole invariance of the action implies so the equations of motion for the dipole Goldstone implement dipole conservation on-shell.
We are welcome to work with the Goldstone as an independent field in the low-energy description, but we find that coupling to background sources is greatly simplified if we instead work with dipoleinvariant combinations of the fields and currents. Using the transformation properties of the dipole Goldstone φ µ , we can define dipole-invariant versions of the U(1) gauge field, symmetric tensor gauge field, and auxiliary dipole gauge field Similarly, we can define the field strength associated with the effective dipole gauge field We can also define dipole-invariant versions of the symmetry currents which are conjugate to variations of the dipole-invariant sources, (3.40) All dependence on φ µ in this expression appears implicitly via the dipole-invariant combinations of various objects. The dipole-invariant currents satisfy an analogous set of conservation equations where the dipole-invariant Lorentz forcef µ and dipole stress tensorτ µν d are now given as We should note that the dipole conservation equation in this interpretation is a consequence of the dipole Goldstone equation of motion, since all fields in (3.40) are manifestly dipole-invariant.
In (3.39), we see that when π µ = 0, the charge density n ν J ν = q is nonzero, and the clockform torsion F n µν vanishes, then the momentum density is proportional to the Goldstone boson, It turns out that this identity holds in a homogeneous equilibrium state, but only under the assumptions mentioned above. A novel feature of our construction is that it allows us to discuss the momentum current and dipole Goldstone separately. This is also crucial when there is an equilibrium superflow, ∂ i φ j = 0. For our construction of superfluid hydrodynamics to be consistent, the system must be linearly stable over some range of finite superflow, even in the absence of an underlying boost symmetry (an argument for standard U(1) superfluids is in [56]). We will explore this feature in an upcoming work [57]. Now that we have discussed coupling the dipole Goldstone to a background spacetime, we can move ahead with our hydrodynamic construction.

Derivative counting
A distinctive feature of hydrodynamics for dipole-symmetric systems is the derivative counting of the hydrodynamic variables, symmetry currents, and background fields. For example, p-wave dipole superfluids have a dynamical scaling exponent z = 2, with time derivatives counting as second order in spatial derivatives . This derivative counting will be crucial in the remainder of this manuscript, therefore we spend some time to arrive at it systematically.
To arrive a consistent derivative counting, we need to postulate the spectrum of hydrodynamic variables for these fluids. As we mentioned in the Introduction, we postulate a variable for each conserved quantity, a local temperature T , chemical potential µ, and velocity u i , in addition to a dipole Goldstone φ i . We are interested in describing thermodynamic systems with nonzero charge density, energy density, and pressure, so we enforce the leading order terms in J t , ǫ t , and τ ij to be O(∂ 0 ). We take the local temperature T and chemical potential µ to similarly be O(∂ 0 ). In dipole-symmetric systems there are no homogeneous equilibrium states with nonzero charge density and finite velocity [16]. Therefore, velocity is always a fluctuation around an equilibrium state, and for this reason we take the spatial velocity u i to be O(∂ 1 ). As with the Goldstone mode of an ordinary superfluid, we take the dipole Goldstone to be of O(∂ −1 ). 7 This counting for the dipole Goldstone field is the same as in recent works [35,36].
The form of the dipole-invariant momentum densityπ µ in (3.39) implies that π µ contains a term proportional to φ µ , meaning that π µ scales as O(∂ −1 ). The source-free momentum conservation equation ∂ t π i +∂ j τ ij = 0, together with the fact that τ ij ∼ O(∂ 0 ), then implies that time derivatives count twice as much as spatial derivatives. The energy and U(1) conservation equations further imply that the fluxes ǫ i , J i are both of O(∂ 1 ). To summarize, we arrive at a counting for the leading order symmetry currents, Given that the generating functional W in (3.25) should be counted as O(∂ 0 ), this in turn induces a counting prescription for the background fields conjugate to the currents, This appropriate derivative counting of background spacetime fields is novel to our work.
With this in place, one can check that the dipole-invariant quantitiesǫ t ,ǫ i ,π i ,τ ij ,Ã t ,Ã i , and a ij have the same derivative counting as their physical counterparts. It is also useful to note the derivative orders of different components of the connection Different derivative counting between space and time components of various objects makes this derivative counting scheme a bit unnatural to implement covariantly. Therefore, it is convenient to In the end, this covariant scheme gives rise to the same physics as the non-covariant scheme outlined above. For the conserved currents, this prescribes Note that the covariant charge and energy density n µ J µ and n µ ǫ µ are still counted as O(∂ 0 ) in this covariant scheme. For the Goldstone field and thermodynamic parameters, we have where u µ is the covariant version of fluid velocity defined such that u µ n µ = 1. Note that under the covariant scheme, the connection has a uniform derivative order The different countings of different time components can then be understood as arising from different contractions with n µ and v µ , respectively.

Hot dipoles in thermal equilibrium
In this section, we discuss fluids with dipole symmetry in hydrostatic equilibrium. We start with a formal discussion of hydrostatic partition functions [42][43][44]58] (see [28,29,59] for discussion specialised to non-relativistic fluids), and again demonstrate the necessity of spontaneously breaking the dipole symmetry. We discuss ideal fluids with dipole symmetry in hydrostatic equilibrium and outline the roadmap to including derivative corrections.

Hydrostatic partition function
Let us warm up with a formal discussion of systems with conserved dipole moment at finite temperature. The primary object of interest in thermal equilibrium is the partition function(al) that can be used to obtain the thermodynamic charge densities and fluxes in equilibrium. In the grand canonical ensemble, in flat space, the equilibrium thermal state of a system can be described by the partition function where β 0 = 1/T 0 is the inverse temperature of the thermal state and H K is the Hamiltonian operator defined as Here K = (K µ , Λ K ) collectively labels the equilibrium thermodynamic temperature T 0 , velocity u i 0 , and chemical potential µ 0 , via K µ = β 0 (1, u i 0 ) and Λ K = β 0 µ 0 . In the context of field theory, the trace operation above can be understood as a path integral over all configurations of the dynamical fields on a Euclidean space where Euclidean time has periodicity β 0 . Note that, while the definition of the thermal partition function makes explicit reference to a preferred time-coordinate, both H K and W eqb are actually time-independent due to the conservation equations.
As we discussed in the previous subsection, for a system with conserved dipole moment, the momentum density π i shifts under dipole transformations as π i → π i − J t ψ i . Consequently, the chemical potential µ 0 must also shift as µ 0 → µ 0 +u i 0 ψ i while keeping u i 0 unchanged, so as to keep the thermal state invariant. This non-invariance of the chemical potential has far reaching consequences for the thermodynamics of the system. In particular, in the "ordinary" phase where none of the symmetries are spontaneously broken, the grand canonical equation of state p 0 = p(T 0 , µ 0 , u 2 0 ) is no longer allowed to have any dependence on µ 0 as long as u i 0 = 0. Consequently, the associated number density n 0 = ∂p 0 /∂µ 0 is identically zero in equilibrium. This is to be expected of fractonic systems with conserved dipole moment, because a state with nonzero density cannot have nonzero momentum and vice-versa [16]. Therefore, to have non-trivial thermodynamic phases in a system with conserved dipole moment and momentum, the dipole symmetry must be spontaneously broken.
The definition of the thermal state can, in principle, be generalized on any time-independent background. The rationale is that, given enough time, any physical system coupled to time-independent sources will tend to relax back to its thermal equilibrium. To this end, consider coupling the physical system of interest to a set of background sources n µ , h µν , A µ , and in the present case A λ µ , that admit a thermal isometry K = (K µ , Λ K ), i.e.
The isometry requirements also hold for v µ and h µν . The thermal isometry K characterizes the observer with respect to which the system is in thermodynamic equilibrium, and the requirements above are the statement that background sources must be static in the reference frame of the equilibrium observer. Without loss of generality, we can take K µ = δ µ t /T 0 and Λ K = µ 0 /T 0 , however it is convenient to stick to the covariant notation. On a time-independent background, we can define a thermal state via the partition function (4.1) as before, but with the definition of H K generalized to include background sources as Here dΣ µ is the spatial volume element associated with a Cauchy slice transverse to the timelike vector K µ . In the second line, we have written the same expression in non-covariant notation for the benefit of the reader. It can be explicitly checked that H K is invariant under "time-independent" symmetry transformations χ µ , Λ, ψ µ , satisfying Similarly, H K is invariant under the choice of Cauchy slice used to define it, i.e. is time-independent, provided that the Ward identities are satisfied. This final property follows from the fact that the integrand in H K above is conserved, i.e.
The partition function W eqb is now a functional of the background fields and can be used to compute the equilibrium expectation values and correlators of various currents.
When the dipole symmetry is spontaneously broken, the hydrostatic generating functional [42,43] can be expressed as a Euclidean path integral over the hydrostatic configurations of the Goldstone φ µ (as for an ordinary superfluid [44] Here, F is the free energy density, which is a scalar under diffeomorphisms and is required to be invariant under U(1) transformations. We have expressed the free energy in terms of dipoleinvariant fields, so it is manifestly invariant under dipole shift symmetry. For later use, let us identify the local gauge-and dipole-invariant versions of thermodynamic parameters: local temperature T , velocity field u µ , and U(1) chemical potential µ as which are all invariant under time-independent symmetry transformations. With these definitions, the velocity is normalized as n µ u µ = 1. Formally, we take K µ ∼ O(∂ 1 ), so that the thermodynamic variables have the right derivative counting as suggested in section 3.4. Note that, in non-covariant notation, this means K t ∼ O(∂ 0 ) and K i ∼ O(∂ 1 ). Note that the original thermodynamic definition of µ in (4.2) is not dipole-invariant. We have compensated for this here by defining µ in terms of the dipole-invariant gauge fieldÃ µ , and hence has dependence on the dipole Goldstone φ µ .

Ideal p-wave dipole superfluids in equilibrium
As an example of the construction above, let us look at ideal p-wave dipole superfluids with a spontaneously broken dipole symmetry. These fluids are qualitatively distinct from ordinary ideal fluids due to the disparate derivative counting. In particular, the kinetic term proportional to u 2 = h µν u µ u ν in the free energy gets suppressed to second order in derivatives, which is a nod to the fact that the kinetic energy of a fluid with conserved dipole moment is parametrically small. On the other hand, there are also terms like F 2 n = F µν n F n µν in the free energy that would ordinarily be counted as second order in derivatives, but get promoted to zeroth order in the presence of dipole symmetry.

Thermodynamics and constitutive relations
Ideal fluids are characterized by a free energy density F expressed in terms of zeroth order scalars made out of n µ , h µν ,Ã µ , andÃ λ µ . Because of the existence of zeroth order in derivative tensor fields a µν , F µν , and F n µν , there are in general many scalars one can write down by mutually contracting chains of these tensors. For simplicity, let us focus on terms in the free energy that are at most quadratic in fields. The ensuing constitutive currents will be at least linear in fields. These are the fields that will be relevant for linear response. The free energy density of an ideal fluid is characterized by the thermodynamic pressure where we have identified the zeroth order scalars Note that there is no kinetic energy term proportional to u 2 = h µν u µ u ν in the free energy at ideal order, as there would be in a non-boost-invariant fluid without dipole symmetry [29,54,55]. 8 This is because unlike ordinary fluids, the scalar u 2 is counted as O(∂ 2 ) according to our derivative counting scheme outlined in subsection 3.4.
Using the thermodynamic derivatives of the pressure p with respect to its arguments, we can define the thermodynamic relations dp = s dT + q dµ + 1 2 (4.11) We can identify the usual thermodynamic observables: the entropy density s, charge density q, and the energy density ǫ. The coefficient χ m can be understood as an inverse magnetic permeability, and we have similar coefficients χ n and χ mn related to the clock torsion F n µν . These coefficients normally show up at second derivative order in an ordinary fluid, but have been promoted to ideal order due to the derivative counting scheme of dipole fluids. In addition, we also have new thermodynamic coefficients specific to systems with conserved dipole moment: the dipole pressure p d and dipole "shear modulus" G d . 9 Eventually, we will be interested in linearized hydrodynamics which will be sensitive to at most quadratic fluctuations in the free energy. To this end, we can also define the thermodynamic derivatives of s, q, and p d leading to the susceptibilities (4.12) Varying the partition function with respect to the dipole-invariant background sources, we can read out the equilibrium constitutive relations for the dipole-invariant quantities (4.13) 9 An analogy can be drawn here with the theory of elasticity which also features a vector degree of freedom similar to δφi, i.e. the displacement field. The analogue of the quantityãij = aij + 2∂ (i φ j) (in the absence of the spacetime sources) is the strain tensor. The coefficients p d , B d , and G d map to lattice pressure, bulk modulus, and shear modulus respectively; see e.g. [60,61]. Similarly, χ sd /B d and χ qd /B d are identified as thermal and charge expansion coefficients. However, unlike the theory of classical elasticity, we also have dependence on the analogue of antisymmetric strain tensorFij = Fij + 2∂ [i φ j] . It will be interesting to further explore this analogy in the context of fracton-elasticity duality [13,62,63], which we leave for future work. Table 1: Transformation properties of various quantities under charge conjugation C, parity P, and time-reversal T discrete transformations. We have only listed the independent components here; the transformation properties of remaining components can be obtained from here using normalization conditions.
Here, the indices onã µν andF µν are raised using h µν , and angular brackets denote a symmetrictraceless combinationã µν =ã µν − 1 d h µν trã. While deriving the above constitutive relations, we have used the fact that δφ µ = h ν µ δφ ν + n µ φ ν v ρ δh νρ due to the normalization condition. The grayed out terms above are non-linear in fields. Since we decided to ignore cubic and higher non-linear terms in the free energy density (4.9), the non-linear contributions to the constitutive relations are not complete.
These constitutive relations generically depend on the equilibrium configurations of φ µ implicitly throughã µν andF µν . The respective classical configuration equations can be obtained by extremizing the free energy, leading to (4.14) These relations validate our derivative counting scheme where it is observed that u µ is suppressed to O(∂).
It is interesting to note the transformation properties of the constitutive relations under discrete symmetries that the system might possess; see table 1. In particular, for systems respecting T symmetry, the coefficients p d , χ sd , and χ qd must vanish in equilibrium, while B d = −2∂p d /∂ trã can be non-vanishing. Another interesting case is CT symmetry. Since the equilibrium chemical potential µ 0 flips sign under CT symmetry, this requires that q, χ sq , χ qd , and χ mn must be odd functions of µ, while all the remaining objects must be even functions of µ. Note that unlike T symmetry, CT symmetry does not require any thermodynamic parameters in the theory to vanish. These discrete transformation properties will have important physical consequences in our subsequent hydrodynamic discussion.

Equilibrium configurations
To gain some intuition into the constitutive relations outlined above, let us turn off the background fields and look at the configuration equations for φ i . From (4.14) it follows that The scalar coefficients in this equation are further functions of We are interested in equilibrium states with vanishing equilibrium velocity u i 0 = 0 and non-vanishing charge density q = q 0 . In this case, the generic ground state solution for φ i takes a linear form for arbitrary constant ξ ij and π i 0 . The respective constitutive relations are where ξ S ij = 2ξ (ij) and ξ A ij = 2ξ [ij] . The coefficients appearing here are understood to be evaluated on the solution for φ i . Note that even though the dipole-invariant momentum density is zero, the physical momentum density is nonzero due to nonzero expectation value of φ i . We find the physical momentum density and stress tensor The energy density/flux remain the same as their invariant counterparts in flat space. It is also possible to find more general solutions, but we will not explore these in this work.
Looking at the explicit equilibrium solutions allows us a further opportunity to verify our derivative counting scheme. Since the generic solution for φ i in (4.17) involves a linear term in the spatial coordinates x i , it follows that φ i ∼ O(∂ −1 ).
Lastly, we note that not all equilibrium configurations in (4.17) are thermodynamically equivalent. The thermodynamically favored solution is one that minimizes the free energy. Throughout this work, we will mostly be interested in the zero superflow state φ i = π 0 i /q, ∂ i φ j = 0. This state is thermodynamically stable if when evaluated on the solution φ i = π 0 i /q. When these stability conditions are not satisfied, the system will prefer to transition to a nonzero superflow state ∂ i φ j = 0.

Derivative counting and higher-derivative corrections
The ideal fluid free energy can generically admit higher-derivative corrections. The most physically relevant among these is the kinetic term that appears at second order in derivatives, i.e.
where ρ is the kinetic mass density. In a Galilean fluid, ρ is proportional to the charge density q, while in a relativistic fluid, it is proportional to ǫ + p. Since a fluid with conserved dipole moment does not have any boost symmetries, ρ is generically an independent thermodynamic coefficient.
Since we are only interested in the free energy up to quadratic order in fluctuations, we can take ρ to be a constant. Including the kinetic energy term modifies the momentum density to take more conventional formπ The contributions to all other currents is at least quadratic in fluctuations. Interestingly, the kinetic energy term is not the only hydrostatic scalar we can introduce in the free energy density at second order in derivatives. There are many more hydrostatic coefficients at this order which are as important as ρ that must be taken into account for consistent truncation of the derivative expansion. These contributions are represented by ellipses in the expression above.
More generally, we can use the derivative counting scheme outlined in section 3.4 to arrange the free energy density F in a perturbative series. Schematically, we have where S (n) I denotes the collection of all independent hydrostatic scalars at O(∂ n ), while the coefficients p and λ (n) I are arbitrary functions of O(∂ 0 ) scalars given in (4.10). As it turns out, there are no hydrostatic scalars at O(∂ 1 ) for a parity-preserving fluid. There are a total of 27 independent hydrostatic scalars at O(∂ 2 ) in addition to u 2 that can affect the linearized constitutive relations; we have worked these out in Appendix A. These lead to corrections to the ideal fluid constitutive relations, which enter at the same order as those coming from the kinetic term 1 2 ρ u 2 .
As an illustration, let us look at all the hydrostatic scalars that affect the momentum density. There are seven such scalars in addition to u 2 , leading to (4.24) Here, R = R λ µλν h µν is the Ricci scalar and we have dropped the superscript "(2)" from the transport coefficients for ease of reading. This results in the full hydrostatic momentum density up to second order in derivatives and linear order in fluctuations . These and the 21 other hydrostatic coefficients also contribute to the remaining conserved currents. However, we will not report these expressions here. Instead, in the next section we will depart from hydrostatic equilibrium and construct a hydrodynamic theory with conserved dipole moment.

Hydrodynamics of p-wave dipole superfluids
Using the ingredients outlined in sections 3 and 4, we will now develop a theory of hydrodynamics for p-wave dipole superfluids using the techniques of [25,27,29,42,43,54,55,[64][65][66]. We will explicitly write down the constitutive relations for our hydrodynamic theory up to subleading orders in derivatives and up to linearized order in fields. We will use this theory to obtain the dispersion relations and response in section 6.

The setup and the adiabaticity equation
Thermal fluctuations of a many-body system around thermal equilibrium are described by the framework of hydrodynamics. The starting point of hydrodynamics are the conservation equations associated with the global symmetries of the system under consideration, which can be used to obtain the dynamical evolution of various thermodynamic observables. For our case of interest, the conservation equations are given by (3.29) and the relevant hydrodynamic fields are given by a set of symmetry parameters:B = (β µ , Λ β , ψ β µ ). In thermal equilibrium, the hydrodynamic fields are aligned along the thermal isometry asK = (K µ , Λ K , 0), meaning that the fluid flow is stationary with respect to the thermal observer. However, out-of-equilibrium, the hydrodynamic fields take a meaning of their own. Note that we have introduced one hydrodynamic field per conservation equation, so the system of equations is closed provided that we specify the hydrodynamic constitutive relations, i.e how the conserved currents of the theory are fixed in terms of the dynamical and background field content. The hydrodynamic fields are valued in the Lie algebra of the symmetry group and transform under an infinitesimal symmetry transformationX = (χ µ , Λ, ψ) as More details regarding these symmetry transformations can be found in [38]. Since the dipole symmetry is spontaneously broken in a p-wave dipole superfluid, we also have the associated Goldstone φ µ in the hydrodynamic theory. This will need to be supplied with its own equations of motion, known as the Josephson equation. For now, we take this to have a schematic form where X µ is some vector field satisfying X µ n µ = 0. In hydrostatic equilibrium, X µ could be obtained by varying the free energy with respect to φ µ . We do not have access to a free energy description out-of-equilibrium, so we will need to implement a different way to find X µ . We will get to it by the end of this subsection; let us keep X µ as an arbitrary operator for now. Together, the dynamical fields can be used to construct dipole and gauge-invariant hydrodynamic fields identified as fluid temperature T , velocity u µ , U(1) chemical potential µ, and dipole chemical potential ̟ µ respectively. We will shortly see that ̟ µ can be integrated out of the hydrodynamic theory using the Josephson equation for φ µ . The situation here is similar to ordinary superfluids where the U(1) chemical potential µ can be integrated out using the Josephson equation for the U(1) Goldstone φ.
To complete the hydrodynamic set of equations, we need to provide a set of constitutive relations for how ǫ µ , π µ , τ µν , J µ , J µν , and X µ are fixed in terms of the dynamical fields β µ , Λ β , ψ β µ , φ µ and the background fields n µ , h µν , A µ , A λ µ . The constitutive relations are required to satisfy the local version of the Second Law of thermodynamics, i.e. there must exist an entropy current s µ whose divergence is positive semi-definite, i.e.
for all solutions of the conservation equations (3.34), when the dipole Goldstone field φ µ is kept offshell. Note that the Second Law is just an inequality statement, however the fact that it is imposed for every possible hydrodynamic configuration results in strong constraints on the functional form of the constitutive relations. We emphasize that the Second Law is imposed for off-shell configurations of the Goldstone φ µ , which appears to be a stronger requirement than the conventional statement of the Second Law of thermodynamics. Ultimately, this generality serves to fix the Josephson equation (5.2). See [67] for more discussion in this regard.
In conventional hydrodynamics, one writes down the most generic expressions for the constitutive relations and entropy current allowed by symmetries, truncated to some desired order in derivatives, and investigates the constraints resulting from imposing (5.4); see e.g. [40]. However, an alternative and simpler route to allowed constitutive relations involves the construction of an off-shell formulation of the Second Law. The details of this procedure for relativistic and Galilean hydrodynamics can be found in [39,[68][69][70][71], and for the boost-agnostic hydrodynamics relevant to this work in [29]. Here we find an analogue of this off-shell constraint for systems with dipole symmetry, whose final form is given in (5.14). To derive this, consider an off-shell generalization of the Second Law statement by adding combinations of conservation equations A priori, the coefficients multiplying the conservation equations in this expression can be arbitrary. Earlier, we introduced β µ , Λ β , and ψ β µ as a means to solve the conservation equations out-ofequilibrium, though we emphasize that they have no intrinsic physical meaning. However, the off-shell Second Law can now be seen as a particular definition for the hydrodynamic fields, which turns out to be equivalent to the requirement that the hydrodynamic fields are aligned with the thermal isometry K in hydrostatic equilibrium, i.e. β µ = K µ and Λ β = Λ K , and ψ β µ = 0. We will return to the redefinition freedom of hydrodynamic fields in section 5.3. Having fixed these definitions, the off-shell statement (5.5) is entirely equivalent to the on-shell statement of the Second Law of thermodynamics in (5.4). We can also massage the Second Law into a more useful form by defining the free energy current which converts the Second Law statement into the so-called adiabaticity equation Here δB denotes a diffeomorphism, gauge transformation, and dipole transformation along the hydrodynamic fields B = (β µ , Λ β , ψ β µ ). Note that there is no entropy production ∆ in hydrostatic equilibrium and δB = δ K acting on background fields vanishes. Consequently, the right-hand-side of (5.7) is trivially zero in equilibrium and the free energy current N µ is identically conserved. In fact, the equilibrium value of N µ is precisely −FK µ , where F is the hydrostatic free energy density introduced in section 4.
For simplicity and convenience, it is useful to recast above equations in terms of dipole-invariant objects. First, the off-shell Second Law of thermodynamics in (5.5) can be re-expressed as where the dipole-invariant entropy current is defined as The transformation vanishes on-shell, so the two entropy currents are physically indistinguishable. Similarly, we define the dipole-invariant free energy current Using this, the adiabaticity equation becomes where δ B denotes a Lie derivative along β µ together with a gauge transformation along Λ β . This equation encodes the same information as the original adiabaticity equation in (5.7), only expressed in terms of dipole-invariant observables.
Let us return to the Josephson equation for the dipole Goldstone field φ µ . There is a trivial solution of the adiabaticity equation (5.11) given by for some positive coefficient σ φ . Upon imposing X µ = 0, this implies that

(5.13)
In principle, this equation can get further derivative corrections compatible with the adiabaticity equation, denoted by ellipses. This equation is supposed to determine the dynamics of φ µ . In practice, however, we can view this as defining ̟ µ or ψ β µ in terms of φ µ and other hydrodynamic fields. In turn, the dipole Ward identity can be used to obtain the dynamics of φ µ . Upon taking this view, the final term in the adiabaticity equation (5.11) drops out, and we are left with a simpler equation Thus, to obtain the hydrodynamic constitutive relations allowed by the Second Law of thermodynamics, we need to find the most general expressions for the dipole-invariant currentsǫ µ ,π µ ,τ µν , J µ , and J µν in terms of the hydrodynamic fields β µ , Λ β , φ µ and the background fields n µ , h µν , A µ , A λ µ , for some free energy currentÑ µ that satisfy the (5.14) for some positive semi-definite quadratic form ∆ ≥ 0.

Constitutive relations: hydrostatic sector
The constitutive relations allowed by the adiabaticity equation (5.14) can be broadly classified into hydrostatic and non-hydrostatic sectors. The hydrostatic sector concerns the part of the constitutive relations that remains non-trivial in equilibrium after the hydrostatic constraints have been implemented. The hydrostatic constitutive relations can be characterized by the free energy density F introduced in Section 4 and does not cause entropy production. On the other hand, the non-hydrostatic sector concerns part of the constitutive relations that vanish in equilibrium. They are further classified into the dissipative and non-hydrostatic non-dissipative sectors, depending on whether they contribute to entropy production or not.
Ideal fluids are part of the hydrostatic sector. It can be checked that the ideal fluid constitutive relations presented in (4.13) satisfy the adiabaticity equation with the free energy current and ∆ ideal = 0. Note that the additional terms in the first line vanish in equilibrium. The terms in the second line have been added for convenience; they are a total derivative and trivially drop out of the adiabaticity equation. The corresponding entropy current is given bỹ Therefore entropy flow in an ideal fluid is purely in the direction of the fluid velocity.
To derive the hydrostatic constitutive relations associated with more general corrections to the hydrostatic free energy density F, we note the identity and ∆ hs = 0. Therefore, transport in the hydrostatic sector is non-dissipative.

Constitutive relations: non-hydrostatic sector
As discussed earlier, equilibrium hydrodynamics is subject to the hydrostatic constraint δ K (. . .) = 0. Out-of-equilibrium, these constraints do not apply and the constitutive relations pick up new contributions, called non-hydrostatic constitutive relations. Due to the construction of the adiabaticity equation, it is natural to write the corrections in terms of δ B acting on background sources, since they identically drop out in hydrostatic equilibrium when δ B = δ K is promoted to an isometry.
Before we get to the non-hydrostatic constitutive relations, we need to address the issue of hydrodynamic frames. As we have stated, β µ , Λ β are arbitrary fields that we chose to describe the hydrodynamic system and we can always redefine them to suit our needs. These induce redefinitions of temperature, chemical potential, and fluid velocity where n µ δu µ = 0. Incidentally, in writing the off-shell Second Law of thermodynamics (5.5), or equivalently the adiabaticity equation (5.14), and choosing the hydrodynamic fields to be the Lagrange multipliers for equations of motion, we further fixed part of the redefinition freedom available to us. Such partial fixing is known as the "thermodynamic frame" of hydrodynamics. However, this still leaves some redefinition freedom where we only shift the hydrodynamic fields with non-hydrostatic variables; see [29]. Heuristically, since the energy, momentum, and charge conservation equations become trivial in equilibrium, they can be seen as relating certain linear combinations of non-hydrostatic data out-of-equilibrium. Therefore, we can eliminate certain nonhydrostatic data from the constitutive relations by suitably redefining the hydrodynamic variables. A suitable choice involves disallowing any non-hydrostatic corrections to the energy, momentum, and charge densities, i.e. choosing n µǫ µ nhs =π µ nhs = n µ J µ nhs = 0. (5.21) This is the so-called "thermodynamic density frame". Conveniently, this frame is known not to exhibit the unphysical instabilities in the linearized spectrum that can appear in conventional choices like the Landau and Eckart frames [29,[72][73][74]. In addition, we have a similar redefinition freedom in the dipole Goldstone field where v µ δφ µ = 0. This will allows us to additionally fix Having fixed the hydrodynamic frame, non-hydrostatic corrections can only appear in the energy flux h µ νǫ ν , stress tensor τ µν , and the dipole flux J µν . The non-hydrostatic constitutive relations satisfy a simpler adiabaticity equation There are two kinds of solutions to this equation. First, we have non-hydrostatic non-dissipative constitutive relations, where contributions to the adiabaticity equation cancel among themselves and lead to no entropy production. They are naturally left unconstrained by the Second Law requirement. Second, we have dissipative constitutive relations, where terms arrange themselves into a quadratic form that is required to be sign constrained. Note that the Second Law does not disallow any non-hydrostatic transport coefficients; at most it imposes inequality constraints on them.
Let us note the derivative order of the non-hydrostatic data appearing in (5.24). We find . (5.25) This disparate derivative counting between vector and other sectors has interesting consequences for non-hydrostatic transport as we discuss below.

Zero-derivative order
An interesting consequence of the derivative ordering (5.25) is that we obtain a dissipative transport coefficient at the same order as ideal fluids discussed in section 4.2. This is given bỹ where κ is thermal conductivity. The entropy production quadratic form and the associated free energy current are given as The positivity of entropy production requires that the thermal conductivity is non-negative All other non-hydrostatic transport is at least second order in derivatives

Two-derivative order
The non-hydrostatic sector is more interesting at subleading orders in derivative expansion. For systems that respect parity invariance, there are no possible contributions at one derivative order. However, at two derivative order, keeping only terms that contribute to linearized hydrodynamics around a state with u i = 0 and φ i = π i 0 /q, the full set of non-hydrostatic constitutive relations are given by  where the objects appearing above are differential operators There are 19 new coefficients at second order. The coefficients ζ and η are bulk and shear viscosities respectively, which have now been suppressed to second order due to the derivative suppression of fluid velocity. Five coefficients Ω ǫ , σ ǫτ ,σ ǫτ , γ ǫτ , andγ ǫτ appear at second order in ordinary hydrodynamics and still show up at second order in the presence of dipole conservation. 10 On the other hand, there are two coefficients σ ǫ , γ ǫ that would ordinarily show up at third order, but which are now promoted to second order in a dipole-conserving fluid. The remaining ten coefficients σ ǫd ,σ ǫd , γ ǫd ,γ ǫd , ζ τ d ,ζ τ d η τ d ,η τ d , ζ d , and η d are specific to the conserved dipole moment. We emphasize that when expanding around an equilibrium state with u i = 0 or ∂ i φ j = 0, there would be many more transport coefficients.
After some computation, we can check that these constitutive relations satisfy the adibaticity equation with the entropy production quadratic form where we have defined together with the free energy current .

(5.33)
We have ignored the terms that are cubic or higher order in fields while deriving these expressions. Furthermore, to manifest the quadratic form structure of the entropy production ∆, we have included some O(∂ 6 ) terms in ∆, arising when O(∂ 3 ) terms in V µ multiply each other. Formally, we can think of these as some some contributions to ∆ coming from four-derivative order constitutive relations that are not relevant at our derivative order of interest. See [48,75] for more details on this procedure. The first line of (5.31) refers to dissipation in the vector sector of the constitutive relations, the second to dissipation in the scalar sector, and the third to dissipation in the tensor sector.
We see that all 7 barred coefficients drop out from entropy production and can be identified as non-hydrostatic non-dissipative. The remaining 12 second order coefficients, together with the zeroth order thermal conductivity κ, are dissipative and lead to entropy production. The positivity of ∆ requires that 11 The non-hydrostatic constitutive relations are subject to additional constraints in cases where the microscopic system obeys microscopic time-reversal invariance. In the presence of dissipation, the macroscopic constitutive relations are typically not expected to be time-reversal invariant. Nonetheless, provided the microscopic theory under consideration is time-reversal invariant, the two-point response functions obtained from the dissipative hydrodynamic construction must obey certain discrete symmetry properties, known as the Onsager relations. We will discuss these in detail in the next section. For now, let us record the implications of the Onsager's relations. Under T invariance, the following non-hydrostatic transport coefficients must vanish 12 On the other hand, under CT symmetry the following non-hydrostatic coefficients must be odd functions of the chemical potential µ, i.e.
This concludes our discussion of hydrodynamics with spontaneously broken dipole symmetry up to second order in derivatives and linear order in fluctuations. In the next section, we will use this theory to derive the hydrodynamic predictions for mode spectrum and linearized response functions.

Dispersion relations and linear response functions
In this section, we solve the linearized hydrodynamic equations and use these to obtain the mode spectrum and linear response functions of various hydrodynamic operators. In the longitudinal sector, we find a thermal diffusion mode ω = −iD k 2 and pair of "magnon-like" propagating sound modes ω ∼ ±v s k 2 − iΓ s k 2 . On the other hand in the transverse sector, we find a "subdiffusive" shear mode ω ∼ −iD ⊥ k 4 . Subdiffusion is a common feature of kinematically constrained systems as has been emphasized in [32,35,36,76].

Response functions from hydrodynamics
Let us start with a quick overview of response functions in hydrodynamics; for more details see [41]. Within linear response, turning on a small external source δs(t ′ ) for an operator O(t ′ ) at some time t ′ , leads to a small deformation at later times of all operators O(t) which do not commute with O(t ′ ), i.e. [O(t), O(t ′ )] = 0. The deformation can be expressed as where . . . s denotes expectation values in the presence of background sources and the object Ordinarily, imposing T symmetry switches off all non-hydrostatic non-dissipative transport coefficients in a parity-preserving hydrodynamic theory. This is because, ordinarily, all odd/even rank tensors (in spatial indices) in a parity-preserving theory are odd/even under T. In the presence of dipole symmetry, however, we have even rank objects J ij , aij that are odd under T. Consequently, Onsager's relation end up switching off certain dissipative coefficients appearing in J ij . is known as the response function or a retarded correlation function. It is often convenient to talk about the Fourier transform of the response function instead The generalization to multiple operators is straightforward. The poles of the Fourier transformed response function can be used to read off the dispersion relations of the hydrodynamic system. The fact that the response function is only nonzero for t > t ′ , means that poles always lie in the lowerhalf complex ω plane. Given a set of hydrodynamic equations coupled to background sources, we can use these to solve for hydrodynamic observables in terms of the respective background sources and use these to explicitly compute the response functions [41]. This method goes beyond traditional Kadanoff-Martin approach [77], as it allows us to compute response functions of not just conserved densities but also their respective fluxes, while keeping track of the necessary contact terms required to satisfy the Ward identities.
We are interested in response functions for a dipole-invariant system in flat spacetime in the absence of external sources. To this end, we introduce plane wave fluctuations of the hydrodynamic fields Furthermore, we are interested in states with zero superflow, so we introduce fluctuations of the Goldstone field as 13 We also introduce fluctuations of the background sources on top of the flat background, i.e.
The set of sources and operators for our hydrodynamic theory can be read off using (3.25) together 13 The states with nonzero superflow are spatially inhomogeneous, so a simple linearized analysis in terms of Fourier modes is not accessible. Note that the case for spontaneously broken dipole symmetry is qualitatively distinct from the states with nonzero superflow in usual superfluids with spontaneously broken U(1) symmetry, because fluctuations in the U(1)-invariant operators are still homogeneous. In the present case, however, due to covariant derivatives in the definition, the fluctuations of U(1) and dipole-invariant operatorãij is non-homogeneous.
with (3.16) and (3.26) as where the operators are evaluated on the solutions of the equations of motion. The operators presented here are only correct up to linear order in background fluctuations; these are sufficient to compute two-point response functions. Our goal is to solve the hydrodynamic equations to determine δT , δu i , δµ, and δφ i in terms of the background field variations δn t , δn i , δg ij , δv i , δA t , δA i , and δA ij . Having done that, we can substitute these to find operator expectation values O s in terms of background fields and use these to obtain the response functions and dispersion relations of our hydrodynamic model.
If we require the underlying microscopic description of the theory to feature certain discrete time-reversal symmetry, e.g. T or CT, the response functions satisfy the Onsager's reciprocity relations where Θ is a diagonal matrix denoting the eigenvalues of operators under the said discrete transformation given in table 1; see [41] for more details. Note that the overall sign of π i 0 and µ 0 must be flipped under T and CT respectively.

Dispersion relations
Let us first obtain the dispersion relations for our hydrodynamic theory. To get some physical intuition of the subsequent results, it is useful to first ignore the dipole Ward identity and leave the dipole Goldstone φ µ off-shell. The ensuing mode spectrum will be that of a boost-agnostic fluid without dipole symmetry; see [29]. Under this assumption, we can solve the energy-momentum and charge conservation equations for the hydrodynamic fields β µ and Λ β in terms of the dipoleinvariant background sources δs. The equations can be written schematically as Here we have denoted δβ =k i δβ i and δβ i ⊥ = (δ i j −k ik j )δβ j to be the longitudinal and transverse projections of β i in the Fourier basis, wherek i = k i /|k|. Furthermore, δs denotes the dipoleinvariant version of background fields δs, combined with the dipole Goldstone φ µ . HereS andM are matrices of suitable dimensionality. Setting δs = 0, the φ µ -off-shell modes are found by solving detM = detM ⊥ = 0. We find a sound and thermal diffusion mode in the longitudinal sector where Ξ = χ ss χ qq − χ 2 sq . The notationω above is to remind us that these are not the physical modes of our theory. Similarly, we find a shear mode in the transverse sector This mode spectrum should be qualitatively familiar from ordinary hydrodynamics.
Physical dispersion relations in our model are obtained when we put the Goldstone φ µ on-shell by solving the dipole Ward identity. In other words, we solve for δβ µ , δΛ β , δφ µ in terms of the physical background sources δs. The equations of motion can be schematically written as where we have introduced the decomposition for φ i similar to the one introduced for β i above. Setting δs = 0, the physical mode spectrum can be obtained by setting det M = det M ⊥ = 0. In the longitudinal sector, we find the linearly propagating sound modes are suppressed to magnon-like modes and the thermal diffusion mode stays intact, i.e.
In the transverse sector, we find that the shear diffusion mode is suppressed to become sub-diffusive, i.e. ω = −iD ⊥ k 4 . (6.14) Here the magnon velocity and attenuation, and the diffusion constants are given by Note that we do not find any real k 4 contribution to the subdiffusive shear mode in the transverse sector. This fact can be understood as a generic consequence of rotational invariance and reality of the hydrodynamic equations of motion. Hydrodynamic response functions can schematically be expressed as G R (ω, k) ∼ P (ω, k)/Q(ω, k), where P and Q are polynomials in ω and k. We can always rescale these polynomials in a way that the reality and rotation-invariance of hydrodynamic equations implies P * (ω, k) = P (−ω, −k) = P (−ω, k) and Q * (ω, k) = Q(−ω, −k) = Q(−ω, k). In particular, this implies that both P and Q are real polynomials of iω and k 2 , and thus their roots for iω are either real or arise in conjugate complex pairs. In particular, for the roots of Q of the form ω = ω n (k) representing the mode spectrum of the theory, this means that either Re ω n (k) = 0 or that there exists another mode ω = ω m (k) such that Re ω n (k) = − Re ω m (k). The transverse sector of our theory has precisely one mode ω = ω ⊥ (k) and, hence, it follows that Re ω ⊥ (k) = 0, explaining the absence of the real k 4 contribution. It also explains why the pair of magnon-like sound modes in the longitudinal sector have oppositely signed real parts and the thermal diffusion mode is purely imaginary.
The fact that dissipation in the longitudinal sector appears at O(k 2 ), the same order as the thermodynamic contributions, emphasizes the importance of thermal fluctuations in a dipole superfluid. As we have seen, dissipation at this order is controlled by a single transport coefficient κ from the vector sector. Other dissipative transport coefficients from the vector sector, whose signs are not fixed by the positivity of entropy production when κ = 0, will appear at O(k 4 ) in the dispersion relations and contribute to fracton-like subdiffusive dissipation. Setting κ = 0 for illustrative purposes, we indeed find that thermal diffusion mode becomes subdiffusive When κ = 0, the coefficients γ ǫ , σ ǫ are constrained by entropy production to be positive.

Response functions
In this section, we look at the hydrodynamic response functions in the transverse sector, involving the operators π ⊥ , J ⊥ , and ǫ ⊥ . The response functions for other operators in the transverse sector, namely τ ⊥ and J ⊥ , can be obtained from these using Ward identities. The response functions in the longitudinal sector can also be easily obtained using our formalism, but the respective expressions are too involved to report here. We have provided a supplementary Mathematica notebook with this submission for the readers who might want to give it a go. We will also present the results for optical response functions obtained in the limit k → 0, in which limit the theory becomes isotropic and the transverse and longitudinal sectors coincide.
Our formalism can predict accurate results for response functions up to first non-trivial correction in the gradient expansion, which are k 2 corrections on top the leading order result. Further taking into account the isotropy of the hydrodynamic model, this means that we can infer results for G R π ⊥ π ⊥ that are accurate up to O(k 2 ) corrections, those for G R J ⊥ π ⊥ , G R ǫ ⊥ π ⊥ up to O(k 4 ) corrections, corrections. However, due to the sheer number of transport coefficients, the general results are quite complicated to present here even in the transverse sector. These have been included in the supplementary Mathematica notebook.
Let us first look at the response function of dipole-invariant momentaπ ⊥ , obtained by keeping the dipole Goldstone φ µ off-shell by ignoring the dipole Ward identity. We find This resembles the respective result from ordinary hydrodynamics without dipole symmetries. We notice that the pole of this correlator is given by the transverse φ µ -off-shell shear mode in (6.11). Instead, if we take φ µ on-shell by solving the dipole Ward identity, we find the response function of physical momenta π ⊥ , leading to where λ ⊥ = λ 5 + λ 6 and we have defined the denominator There are two poles of this correlator: one given by the subdiffusive hydrodynamic shear mode in (6.14) and another non-hydrodynamic mode at ω ∼ −i/k 2 + O(k 0 ) that is not relevant for the lowenergy spectrum. The dictionary between the dipole-invariant and physical response functions is explained in appendix B. Notice that, in the limit q → 0, the physical response function reduces to its dipole-invariant value, meaning that dipole symmetry has no consequence for response functions in chargeless states. It is also interesting to note that as ω → 0, we recover the hydrostatic response function (6.20) justifying the role of π i as a dipole Goldstone.
For the remaining response functions, we only look at the respective physical versions where the Goldstone φ µ has been taken on-shell. We find the diagonal response functions Our formalism can also predict the O(k 4 ) corrections to these correlators, but these expressions are quite involved and not illuminating. We can also consider the off-diagonal response functions.
For the off-diagonal response involving the momentum density, we find We have used the notation G R+ AB = G R AB and G R− AB = G R BA . Here we see that the off-diagonal terms are symmetric up to three signs of p d and certain dissipative coefficients. Imposing Onsager reciprocity relations as in (5.35) sets these coefficients to zero, leading to symmetric off-diagonal response. Similarly, for the energy-charge off-diagonal correlator we find Again, we have dropped the O(k 4 ) corrections to this correlators for clarity.
Now, let us look at the frequency-dependent correlators at zero wavevector. In this limit, all correlators involving the conserved densities ǫ t , J t , and π i and the U(1) flux J i are zero due to the Ward identities (up to contact terms). The optical correlators then read We can use these to find the Kubo formulae The Kubo formulae for other transport coefficients require us to probe the response functions at nonzero wavevector. We shall not report these in the manuscript.

Towards s-wave dipole superfluids
In our discussion thus far, we have presented a detailed analysis of p-wave dipole superfluids where only the dipole symmetry is spontaneously broken, but the U(1) symmetry is left intact. In this section, we briefly outline how our results can be generalized to s-wave dipole superfluids where the U(1) symmetry is spontaneously broken as well. At nonzero temperature these phases are only expected to exist in spatial dimension d > 4, below which infrared fluctuations of the Goldstone restore the U(1) symmetry at long distance.
We leave a comprehensive analysis to a forthcoming manuscript [30], but record the most salient features here. In particular, we find that s-wave dipole superfluids admit ordinary sound modes and so resemble a model with a dynamical scaling exponent z = 1 instead of z = 2 as for p-wave dipole superfluids. This raises the prospect that the natural derivative counting scheme for this phase is different than in the p-wave phase, which we comment on below and will address in [30].
Hydrodynamic equations.-Instead of an independent vector Goldstone φ µ , s-wave dipole superfluids are expressed in terms of a scalar U(1) Goldstone φ transforming according to The dipole symmetry is automatically spontaneously broken as a result of the spontaneously broken U(1) symmetry, and the associated vector Goldstone is given by Introducing independent Goldstones for U(1) and dipole symmetries is redundant from the perspective of low-energy effective field theory; see footnotes 2 and 6.
Most of the discussion of p-wave dipole superfluids can be extended to s-wave by simply replacing φ µ with its definition above in terms of φ. In particular, after substituting for φ µ , the U(1)-and dipole-invariant versions of the background gauge fields becomẽ where we have defined Φ = v µ (A µ + ∂ µ φ). We note that, unlike the p-wave case,ã µν here mixes both a µν and A µ . Furthermore,Ã µ is only sensitive to the temporal part of the U(1) gauge field v µ A µ . The resulting chemical potential is given as meaning that, in hydrostatic equilibrium when δ B φ = 0, the chemical potential is just Φ.
Taking these observations into account, one can construct a theory of hydrodynamics for s-wave dipole superfluids similar to our discussion in Section 5.1. The hydrodynamical variables can be taken to be a local temperature T , chemical potential µ, velocity u µ (satisfying u µ n µ = 1), and scalar Goldstone φ. The hydrodynamical equations are the Ward identities for energy, momentum, and U(1) current conservation, together with a Josephson condition. At the level of constitutive relations the vector current h mu ν J ν is exactly fixed to be the gradient of the dipole current, h µ ν J ν = ∇ ′ ν J µν , and we supply constitutive relations for the energy current, momentum current, spatial stress density, charge density n µ J µ , and dipole current J µν .
The derivation of the adiabaticity equation also follows exactly the same as before in Section 5.1, culminating into its final form in terms of dipole-invariant objects Recall from above that the spatial part ofÃ µ is not an independent field in the s-wave phase and thus has been substituted in terms of φ in this expression.
An argument similar to (5.13) gives a trivial solution of the adiabaticity equation for some positive coefficient σ φ . Upon imposing the Josephson equation X = 0, we find where ellipses represent further corrections that can be added consistent with the adiabaticity equation. Noting that Φ = v µ ∂ µ φ in the absence of background fields, this equation is reminiscent of the Josephson equation from ordinary U(1) superfluids [78]. We could proceed as we did following (5.13), by treating the Josephson equation as the defining equation for µ and drop the δ B φ in adiabaticity equation (7.5). However, it is customary to fix the definition of µ by forbidding the U(1) density J µ n µ to obtain any non-hydrostatic corrections. In this interpretation, (7.7) should be seen as determining the time-derivatives of the scalar Goldstone v µ ∂ µ φ in terms of µ − v µ A µ and further derivative corrections.
Derivative counting and constitutive relations.-With the hydrodynamic variables and equations out of the way, we require a derivative counting scheme to construct constitutive relations. Given that our derivative counting for p-wave dipole superfluids required that φ µ ∼ O(∂ −1 ), we anticipate φ to be O(∂ −2 ) along with the p-wave scalings for the fields and currents. However, as we shall see below, the linearized mode spectrum of the theory suggests another potential derivative counting scheme where we take φ ∼ O(∂ −1 ) and φ µ ∼ O(∂ 0 ), together with dynamical scaling exponent z = 1, meaning that v µ ∼ O(∂ 0 ) and v µ ∂ µ ∼ O(∂ 1 ). In this counting, all the hydrodynamic fields u µ , T , µ are O(∂ 0 ), similar to ordinary hydrodynamics. Following through the consistency arguments similar to those in Section 3.4, we find that, under this alternative counting, the conserved currents ǫ µ , π µ , τ µν , J µ are all O(∂ 0 ), except J µν that is formally O(∂ −1 ). Correspondingly, the background fields . The same countings would apply for the dipole-invariant versions of these quantities as well.
At this time it is unclear which of these counting schemes is the one that organizes the s-wave hydrodynamics. We will address this question in [30]; for now we proceed agnostically in order to draw out some of the basic physical features of these fluids.
To illustrate the qualitative features of s-wave dipole superfluids, let us write down a representative set of constitutive relations consistent with the adiabaticity equation (7.5). We leave a more thorough analysis for future work [30]. The hydrostatic sector of the constitutive relations is characterized by the free energy density F. Let us include a pressure term, p 0 (T, µ), a "kinetic energy" term ∝ u 2 , a term ∝ trã that contributes to the dipole current one-point function, and terms ∼ã 2 that, in hydrostatics, serve as kinetic terms for the Goldstone φ: We emphasize that this does not characterize the most general free energy up to second-order in derivatives. Following the procedure outlined in Section 5.2, we are led the linearized hydrostatic constitutive relations ǫ µ hs = ǫ u µ + p u µ + p d v ρF ρ µ + . . . , π µ hs = ρ u µ + . . . , τ µν hs = p h µν + . . . , n µ J µ hs = q + . . . , where p = p 0 + 1 2 ρ u 2 and dp = sdT + qdµ Note that we have only written down constitutive relations for the density part of J µ , since its spatial part is identically fixed by the dipole Ward identity. We have also written down "constitutive relations" for the Goldstone equation of motion X in agreement with the adiabaticity equation, as we decided not to impose it identically as a definition of µ.
In flat spacetime and in the absence of background sources, the non-linear configuration equation for the scalar Goldstone is given by This is solved by the generic hydrostatic configuration of the scalar Goldstone which give rise to the same configurations for φ i as (4.17). Note that, in contrast to p-waves, there is no equation in hydrostatic equilibrium enforcing a derivative suppression of the fluid velocity u i .
Similarly, we can also write down non-hydrostatic corrections to the constitutive relations following our discussion in Section 5.3. A representative but non-generic subset of these arẽ

(7.13)
Mode spectrum and order-mixing.-We can use the hydrodynamic equations outlined above in the absence of background sources to obtain the linearized mode spectrum of s-wave dipole superfluids. In contrast to the p-wave phase, we find a normal sound mode in the longitudinal sector with linear dispersion relations ω = ± sχ qq (ǫ + p) ΞρT k + O(k 2 ). (7.14) In the transverse sector, we similarly find the normal shear diffusion mode In particular, we see that the shear mode is not subdiffusive for s-wave dipole superfluids. The above modes suggest the alternative derivative counting scheme we mentioned earlier with z = 1. Note that both the normal modes above crucially depend on the kinetic mass density ρ.
In addition to the modes above, we also find a pair of magnon-like second-sound modes in the longitudinal sector These are reminiscent of the second sound mode in ordinary U(1) superfluids, but has now been suppressed in k due to the presence of dipole symmetry. We have gone beyond what is presented here and included many more derivative corrections to the constitutive relations. While we do not have a conclusive proof at this time, we comment that so far we have only found further corrections to this mode at O(k 4 ). In particular, in contrast to the magnon-like mode observed in the p-wave phase, the attenuation correction to this mode goes as ω ∼ −ik 4 and is purely subdiffusive. In this sense, this mode is "purely fractonic" in origin, and has no qualitative contamination from the non-fractonic hydrodynamic degrees of freedom.
The mode spectrum presented here is consistent with two derivative counting schemes. One is the z = 1 scheme we suggested earlier in this Section, where ω ∼ k, so that the "kinetic energy density" ρ is O(∂ 0 ), the viscosity η is O(∂ 1 ), but the coefficients B d and G d are O(∂ 2 ). However, this scheme is in some tension with a more general hydrostatic analysis, which allows for equilibrium states (7.12) with superflow with φ ∼ x 2 .
The other scheme is that for the p-wave phase together with φ ∼ O(∂ −2 ). In this counting ω ∼ k 2 , B d and G d are O(∂ 0 ), but η and ρ are O(∂ 2 ). While this counting is seemingly at odds with the existence of a sound mode, we note that (1) it is consistent with the existence of equilibrium states with superflow and (2) further is consistent with the linearized spectrum. However, with this counting the kinetic energy density ρ is analogous to a dangerously irrelevant operator in quantum field theory, appearing in the denominator of the dispersion ω = ω(k) and thereby changing the scalings of infrared dispersion relations.
Both of these countings are consistent with the linearized mode spectrum in the absence of superflow. We expect that only one is consistent in the presence of interactions, or equivalently in a superflow, and will report on the matter soon [30].

Discussion
In this work, we have systematically constructed the hydrodynamic description of translationallyand rotationally-invariant fluids with a spontaneously broken dipole symmetry. We have explained our algorithm in great detail already in the Introduction and Section 2. In this Discussion, we wrap up by drawing attention to a few interesting aspects of our analysis and point out a few natural questions and future directions suggested by our work.
First, thanks to the unconventional derivative counting appropriate to fluids with dipole symmetry, dipole superfluids have the unique feature that there is dissipative transport already at the leading order in derivatives. The thermal conductivity κ is formally as important in the hydrodynamic equations as the "perfect fluid" part of the constitutive relations. This leading order dissipative transport is in the vector sector. There is dissipative transport in the scalar and tensor sectors as well, although these appear exclusively through derivative corrections to the leading order constitutive relations. In our analysis, we have classified all 47 transport coefficients consistent with the Second Law of thermodynamics up to first subleading derivative order. In doing so, we are able to deduce the full suite of inequality-type constraints on the transport coefficients enforced by the positivity of entropy production. As in relativistic hydrodynamics, we expect that out-ofequilibrium transport with even more derivatives is automatically consistent with the Second Law once the constraints presented here are satisfied.
In order to implement the Second Law constraints up to subleading order in derivatives, we built upon some existing techniques in the literature. We appropriately generalized the "adiabaticity equation" and the off-shell formulation of the Second Law of thermodynamics from [69] to fluids with dipole symmetry. We illustrated that the constitutive relations can be consistently classified into a hydrostatic part, obtained by variations of a hydrostatic effective action, and a non-hydrostatic part. In particular, using redefinitions of the hydrodynamic variables, we argued that it is possible to choose the non-hydrostatic corrections to only appear in ǫ i , τ ij , and J ij . In doing so, we landed on a relatively simple parametrization of out-of-equilibrium transport given in (5.29) and showed that in this form it is straightforward to impose the positivity of entropy production.
Given our results, there are some natural directions to follow upon. A particularly interesting point to note is that hydrostatic equilibrium allows for the dipole flux J ij to have a nontrivial one-point function J ij = p d δ ij in the absence of dipole superflow. Indeed, the large N models of [16] feature such a thermal one-point function. However, we expect that this is an indication that the thermodynamic state of the fluid is unstable and the system would prefer to transition to a state with nonzero superflow. Furthermore, we have constructed a hydrodynamic description of p-wave dipole superfluids at the level of linearized constitutive relations. Our methods can be generalized to produce the full theory of nonlinear dipole superfluid dynamics, at the expense of keeping track of a dimensiondependent number of zeroth order scalars constructed from the zeroth order tensorsã µν ,F µν , and F n µν . Constructing such nonlinear hydrodynamics has the useful corollary of allowing us to describe transport in a state with superflow characterized by ∂ i φ j = 0. These questions will be explored in more detail in [57].
As another immediate follow-up to this work, we intend to construct a hydrodynamic description of dissipative s-wave dipole superfluids in an imminent publication [30]. The initial steps towards this construction were outlined in Section 7. In particular, we found that the linearized spectrum and natural derivative-counting in this phase is quite different from its p-wave counterpart. It will be valuable to systematically revisit the constitutive relations of s-wave dipole superfluids and inspect the constraints imposed by the second law of thermodynamics.
We noted in our introduction that the dipole Goldstone field φ i is UV-sensitive, in the sense that it is compact with a periodicity ∼ 2π/a, with a being the lattice spacing [43]. The U(1) density and flux J t , J i , the dipole flux J ij , and the energy density and flux ǫ t , ǫ i (in flat space) only depend on the derivatives of φ i , and are therefore UV-insensitive. On the other hand, the momentum density π i and the spatial stress tensor τ ij depend directly on φ i , making these UV-sensitive. It would be interesting to explore what consequences does this sensitivity imply for the two-point functions of π i and τ ij , accounting for one-loop fluctuations of the dipole Goldstone.
In two spatial dimensions, the low-energy limit of elasticity theory is expected to be described by a dipole-symmetric phase, with disclination defects in the crystal playing the role of elementary charges [13]; see also [62]. As such, we also expect our hydrodynamic description to describe a phase of liquid crystals with dislocations, seen as condensed disclination-antidisclination pairs. It would be interesting to see how this works out in detail and whether the low-energy description of these physical systems can indeed be reinterpreted as dipole-invariant field theories. More generally, it would be interesting to explore the relations between our results and transport in elastic media.
In nature, we expect any dipole symmetry to be emergent at low-energies, broken by irrelevant operators in the IR description. In a regime where the dipole symmetry is approximately preserved, such systems should still possess a hydrodynamic description at small temperatures, but with the explicit breaking of the dipole symmetry generating a small mass for the dipole Goldstone. It would interesting to investigate whether a hydrodynamic theory of a fluid with a light dipole Goldstone along the lines of [79,80] can be turned into a diagnostic of these phases.
Lastly, in this work we have only considered the hydrodynamic description of systems with dipole symmetry. In the context of fractons, dipole symmetry is only a toy model for subsystem symmetries realized in lattice models with fractonic excitations. It would be interesting to explore whether subsystem-invariant theories also have phases which admit a hydrodynamic description like those considered in this work, and if so, how do these qualitatively differ from the ones with dipole symmetry.

A Hydrostatic scalars
In this appendix, we will enlist all hydrostatic scalars up to second order in derivatives that can appear in the free energy density F and characterize the hydrostatic constitutive relations. In particular, we will only enlist scalars that are at most quadratic in fluctuations. With these considerations, we find a total of seven hydrostatic scalars at zero-derivative order All transport coefficients in the theory can be seen as arbitrary functions of these zero-derivative scalars. The first three of these are linear in fluctuations, while the remaining four are quadratic. For a fully non-linear theory, there are other zero-derivative hydrostatic scalars one can construct combining chains of zero-derivative tensors, e.g.F µνF νρF ρσF σµ . Furthermore, for a parity-violating theory in even spatial dimensions, we can additionally construct zero derivative scalars as exterior products ofF µν and F n µν , such as where ε µ 1 µ 2 µ 2 ...µ d is the spatial Levi-Civita symbol satisfying n µ 1 ε µ 1 µ 2 µ 2 ...µ d = 0. This term is quadratic in fluctuations in d = 4 and linear in d = 2. In odd spatial dimensions, we instead have one-derivative terms such as These are at least quadratic in fluctuations in d = 3 and linear in d = 1. We are focusing on parity-preserving fluids in this work, so we have not considered these parity-violating terms in the main text. Furthermore, due to our derivative counting scheme, there are no hydrostatic scalars at one-derivative order for a parity-preserving theory. So, the first corrections to the ideal order free energy appears at second order in derivatives.
We will now perform the counting of parity-preserving hydrostatic scalars at second order in derivatives and at most quadratic in fluctuations. First, we have one hydrostatic scalar that does not involve explicit derivatives Next, let us look at the hydrostatic scalars involving one explicit derivative. Hydrostatic identities δ B (. . .) = 0 require that u µ ∂ µ S (0) = 0. In addition to these, we have one more scalar hydrostatic identity relevant at second order in derivatives, ∇ ′ µ β µ = 0. We have 3 hydrostatic scalars of the kind u µ ∂ µ S (0) , i.e.
We have not counted scalars likeū µ ∂ µã 2 etc., because these are cubic in fluctuations. In principle, we also have a scalar like ∇ ′ µ u µ . However, note that within an integral, for some transport coefficient λ, we find using integration by parts that Therefore, it is not an independent hydrostatic scalar for the purposes of the hydrostatic partition function. We will use similar arguments to kill other "total-derivative" scalars in our hydrostatic counting. Similarly, we have scalars like u µ ∇ ′ ν S µν , where S µν is a placeholder for zero-derivative tensorsã µν ,F µν , and F µν n . Using the same procedure as above, these can be replaced with S µν ∇ µ u ν , resulting in 3 scalars Next, let us consider hydrostatic scalars with two explicit derivatives. There are scalars of the kind ∇ µ S (0) ∇ µ S (0) , where we have defined ∇ µ = h µν ∇ ν . Up to quadratic order in fluctuations, this yields 6 possibilities Using the integration by parts argument outlined above, we can skip the total-derivative hydrostatic scalars of the kind h µν ∇ µ ∇ ν S (0) . Likewise, we can skip scalars of the kind ∇ µ ∇ ν S µν in favor of ∇ µ S µν ∇ ν S (0) , which in turn can be skipped in favor of S µν ∇ µ ∇ ν S (0) . This results in three more hydrostatic scalars µν F n µν u λ ∂ λ S (0) are quartic in fluctuations and hence are not important for our purposes. In addition to these, there are 4 hydrostatic scalars of the kind ∇ µ S νρ ∇ µ S νρ , i.e.
24 ≡ ∇ µãνρ ∇ νF µρ , S 26 ≡ ∇ µ F n νρ ∇ νF µρ . (A.11) We have skipped scalars of the kind S µν ∇ ρ ∇ ρ S µν , S µν ∇ ρ ∇ µ S ρν , and ∇ ρ S ρν ∇ µ S µν using integration by parts. Finally, there are two scalars constructed out of the Riemann tensor S This exhausts the list of potential hydrostatic scalars. In summary, at second order in derivatives, there are 28 hydrostatic scalars which must be added to the ideal order pressure to capture linearized hydrodynamics. Interestingly, with the exception of the Ricci scalar S (2) 27 , all other two-derivative hydrostatic scalars in our counting are quadratic in fluctuations. This means that we can take the respective transport coefficients to be constants for our purposes.

B Computation of response functions
In this appendix, we will outline an analytic path to computing hydrodynamic response functions of our model. We will break this problem into two steps. We recall that the conservation equations of the system can be written in terms of dipole-invariant objects in (3.41). Therefore, we will first use the dipole-invariant energy, momentum, and charge conservation equations in (3.41) to determine δT , δu i , and δµ in terms of the spacetime background fields δn t , δn i , δg ij , and δv i , and the dipole-invariant versions of the gauge fields δÃ t , δÃ i , and δã ij . Pretending for the moment that δÃ t , δÃ i , and δã ij are background sources, we can use these solutions to compute the response functionsG R OÕ between dipole-invariant objects through (B.1) The answers thus obtained will, of course, not be the physical response functions because we still need to solve for δφ i . Nonetheless, the dipole-invariant response functions are more analytically tractable and can be used as building blocks to construct the full physical response functions. The second part of the problem will be to solve the dipole conservation equation in (3.41) for δφ i in terms of the physical background fields. These can be used to convert the dipole-invariant response functions into their physical counterparts and also obtain the hydrodynamic mode spectrum.
Let us imagine that we have obtained the dipole-invariant response functionsG R OÕ . We will like to relate this to the physical response functions that will actually be observable in an experiment. Firstly, let us note the transformation between dipole-invariant and physical background fields Note that the equilibrium value ofÃ i is π 0 i /q, while that ofÃ t andã ij is zero. Using these together with (6.7) and (B.1), we can obtain the dictionary between physical and invariant operators up to linear order in fluctuations, i.e.
These transformation rules can be summarized as Here " †" denotes the transpose conjugate operation and we have identified the transformation coefficients as The free row indices are denoted by i, j while the free column indices by k, l. Note that C † = C. Note also that for an equilibrium state with zero momentum density, i.e. π i 0 = 0, these simplify to R = 1 and C = 0, while X i and Y i are independent of π i 0 . Using (B.3) and noting that Õ s only depends on δs, we can write down the relation between physical and dipole-invariant response functions All we need now to complete the mapping are the solutions for φ i using the dipole conservation Ward identity. At linearized order in momentum space, it is given as Since both J ij and J i are dipole-invariant operators, we can express this in terms of dipole-invariant response functions as where the coefficient matrix M ij is defined as the inverse of (M −1 ) ij = X †i ·G R OÕ · X j . (B.8) Plugging this back into (B.5), we find Note that the terms inside the square brackets are independent of π i 0 , which are nothing but the physical correlation functions at zero momentum density. The physical correlation functions at nonzero momentum density can be obtained by performing a global dipole transformation and eventually result in an operator rotation R and a contact term C.

Onsager relations
When the underlying microscopic description has some discrete time-reversal symmetry, the dipoleinvariant response functions are also required to satisfy the same Onsager relations as the physical ones in (6.8), i.e. G R OÕ (ω, k) T = Θ · G R OÕ (ω, − k) · Θ, (B.10) together with the Onsager constraints for the coefficient matrices (B.11) required for consistency of the transformation rules (B.9).

C Transformation properties of conservation equations
In this appendix, we give some details regarding the equations of motion. Let us express the conservation equations for off-shell configurations of the dipole Goldstone field φ µ as where f µ and τ µν d are defined in (3.35). Using (3.28), we note that f µ and τ µν d transform under dipole transformations as The conservation equations E q and E µ d are manifestly invariant under dipole transformations. For the other two, after some non-trivial algebra, we can show that they transform as In deriving these, we have used the curvature identities (C.4) Therefore the full system of conservation equations is invariant under dipole transformations.
As we motivated in the main text, the dipole-invariant combinations of energy current, momentum density, and stress tensor themselves satisfy a set of conservation equations wheref µ andτ µν d are defined in (3.42). The dipole-invariant quantities are obtained by performing a dipole transformation along φ µ . Therefore, using (C.3) we can read out the relations (C.6) We can express the off-shell Second Law of thermodynamics in (5.5) in terms of the quantities outlined above, i.e.
It can be checked that the Second Law statement is invariant under dipole transformations, provided that we also shift the entropy current as Note that this transformation vanishes on-shell, so the physical entropy is still invariant. We can also note the transformation properties of the free energy current defined in (5.6), i.e.
We have used the fact that 1/ √ γ ∂ ν ( √ γ X µν ) = ∇ ′ ν X µν + v µ 1 2 F n ρσ X ρσ . The second term in the transformation is irrelevant because its gradient vanishes and does not contribute to the adiabaticity equation. However, the free energy current still transforms non-trivially under dipole transformations even for on-shell configurations, and the transformation law is similar to the conserved energy current ǫ µ as given in (3.28).

D Comparison with previous works
Our results for p-wave superfluids have some overlap with previous works of [35,36], and in this Appendix we would like to emphasize the commonalities and distinctions between our works.
We have argued at length in this work that the global dipole symmetry of a physical system must be spontaneously broken in the low-energy description for there to exist non-trivial charge transport, leading to a vector Goldstone φ i . In the absence of dipole superfluid, i.e. when ∂ i φ j = 0, we can identify the one-point function of momentum density with that of the dipole Goldstone, i.e. π i = φ i /q 0 . In [35], this identification is posited even out-of-equilibrium, though, as we show, this is not true on curved backgrounds, or whenπ i = 0, for instance in the presence of a kinetic mass density ρ. By utilizing the dipole-invariant construction of hydrodynamics, we are able to relax this assumption and explicitly identify the role of the Goldstone field in the momentum response beyond equilibrium thermodynamics. A tangible feature of this is the appearance of the kinetic mass density ρ in momentum response function G R π ⊥ π ⊥ . As we emphasize in the main text, a consistent formulation of p-wave dipole superfluid hydrodynamics involves a novel derivative counting scheme in the presence of background sources. Fundamental to this counting is the requirement that the hydrodynamic fields µ and T scale as O(∂ 0 ). The existence of superflow indicates the Goldstone φ i scales as O(∂ −1 ), from which it follows via the homogeneity of hydrodynamic equations that A i ∼ O(∂ −1 ) and u i ∼ O(∂). Note that while u i is a hydrodynamic field, the dipole Ward identity requires that it vanishes in equilibrium, as we demonstrate in (4.15). In fact, this counting scheme is implicitly present in the work of [35] as well, as illustrated via the inclusion of thermal conduction at the same derivative order as the ideal order hydrostatic terms in the dipole-invariant energy fluxǫ i .
Next, the authors of [35] noticed that the use of a shifted energy current and effective chemical potential led to a more familiar expression for the first law and allowed for more straightforward Second Law constraints. In fact, these are none other than the dipole-invariant versions of the energy fluxǫ i and our dipole-invariant definition of the chemical potential. That such a rewriting is useful then becomes transparent as the hydrodynamic theory written in terms of dipole-invariant fields is effectively that of a charged boost-agnostic fluid without dipole symmetry.
Conveniently, we find that with a certain consistent truncation, our results match with those of [35]. For the sake of matching, we will work in flat spacetime background. To compare, consider the ideal constitutive relations for the U(1) and dipole currents from (2.5). The linearized Goldstone equations of motion then imply eq. (4.15). So, denoting the notation of [35] with a hat, we identifŷ and becauseV [ij] = 0 in [35], to match, we then must set χ m = 0. This is the equivalent of settinĝ F ij = −J ij , an equality obtained from the ideal Second Law in their work. Continuing, we see that dp = sdT +ndµ − p d δ ij − G d a ij dV ij , (D.2) so that δF ij = δ ij B d δV k k + χ sd δT + χ qd δµ + 2G d δV ij . (D. 3) which tells us to identifyf = B d andf ⊥ = 2G d , as well as set χ sd = χ qd = 0. These are reasonable if in our model, we choose p d = − 1 2 B d trã which then vanishes in equilibrium; see (4.12). We can also see thatf =n −1 0 (f + d−1 df ⊥ ) = B ′ d /q 0 in our notation. Having made these identifications, we can consider constitutive relations. In the absence of background sources our hydrostatic dipole-invariant energy current is agreeing with [35]. Importantly, to match the two expressions, we must include certain non-linear terms inǫ i . Including non-hydrostatic corrections, in flat spacetime, we find This allows us to identifyα = T 2 κ.
Turning to higher order non-hydrostatic terms, whereas we predict 7 non-trivial non-dissipative terms and [35] predicts 6, both our works predict 12 non-trivial transport coefficients which contribute to entropy production. In fact, once we have integrated by parts, the expressions for entropy production are quite similar, being expressible in terms of ∂ i ∂ j T , ∂ i ∂ j µ, and ∂ iVj , though the relation between transport coefficients is somewhat obscured. The extra coefficient we find isΩ ǫ which explicitly contains a time derivative. It seems reasonable that using equations of motion, this can be replaced with spatial derivatives, but whether or not this term can be eliminated altogether at second order was not evident to us.
After imposing time reversal invariance, we are left with only 6 dissipative coefficients, as opposed to 12, and we retain 6 non-dissipative coefficients at second order (now,Ω ǫ = 0). Hence, our total counts after Onsager's relations are imposed are the same, but our accounting of entropy production differs slightly. This is partially related to our choice of hydrodynamic frame, where we impose h µν J µ nhs = 0. Nevertheless, this distinction does not affect the leading behavior of the hydrodynamic modes in k. Given the matching discussed above, where we turn off many coefficients, we find the same expressions for the longitudinal modes to O(k 2 ) and the transverse mode to O(k 4 ).
In summary, the work of [35], can be seen as a consistent truncation of the work presented here, subject to the caveat that the relationship between the momentum π i and the Goldstone φ i is required to hold out-of-equilibrium. In this respect, our method can be used to find response functions for the model presented in that work, and can be used as a diagnostic to identify relevant contributions to the hydrodynamics. For instance, when ρ, p d = 0, then lim ω→0 −1 ω 2 Re G R ǫ i ǫ j (ω) = ρp 2 d q 2 , (D. 6) whereas for [35] would find a vanishing result. Likewise, in the limit that q → 0, we would recover the response of a neutral fluid G R π ⊥ π ⊥ = −ρ + ωρ 2 ρω + iηk 2 . (D.7) We would also like to mention the work of [36] which found similar results to our work and to the results of [35] in the context of p-wave dipole superfluids. The work considers dissipation up to O(∂ 2 ) in the limit that energy is not conserved, as well as hydrodynamics with energy conservation up to O(∂) via inclusion of the thermal conductivity term. Overall, the schematic structure of the hydrodynamic modes is similar to what we find, but because of the non-trivial mixing between energy, charge, and the dipole current, we found that it was necessary to consider all currents simultaneously within a consistent derivative scheme. We nevertheless note a curious observation that the subdiffusive constant D ⊥ is oblivous to explicit derivative corrections appearing in the non-hydrostatic sector in the sense that the only relevant dissipative parameter is η, but not terms like η d which appear at the same order in the tensorial sector. In fact, the subdiffusion constant matches the result of [36] given in the absence of energy conservation, suggesting that the effects of energy conservation are subleading to the coupling of stress to dipole fluctuations in the transverse sector, leading to a universal propagation behavior. This is connected to the fact that the dissipative vector terms appear O(∂ −2 ) enhanced relative to the scalar and tensor terms and the energy flux only sees the vector corrections.