Hydrodynamics of dipole-conserving fluids

Dipole-conserving fluids serve as examples of kinematically constrained systems that can be understood on the basis of symmetry. They are known to display various exotic features including glassylike dynamics, subdiffusive transport, and immobile excitations dubbed fractons. Unfortunately, such systems have so far escaped a complete macroscopic formulation as viscous fluids. In this work, we construct a consistent hydrodynamic description for fluids invariant under translation, rotation, and dipole shift symmetry. We use symmetry principles to formulate a thermodynamic theory for dipole-conserving systems at equilibrium and apply irreversible thermodynamics in order to elucidate dissipative effects. Remarkably, we find that the inclusion of the energy conservation not only renders the longitudinal modes diffusive rather than subdiffusive but also diffusion is present even at the lowest order in the derivative expansion. This work paves the way towards an effective description of many-body systems with constrained dynamics such as ensembles of topological defects, fracton phases of matter, and certain models of glasses.


I. INTRODUCTION
Over the years, hydrodynamics has evolved into a universal framework describing the long-wavelength dynamics of many-body systems. It provides a systematic scheme for the evolution of conserved charges, leading to an effective description that is for the most part irrespective of the microscopic details. The structure of a hydrodynamic theory is determined by the symmetries of the system at hand, once the appropriate low-energy degrees of freedom have been identified and a finite set of phenomenological transport coefficients introduced. This approach gives access to the low-energy (long-wavelength) physics of complex many-body systems that are often impossible to be understood from first principles. Over time, the hydrodynamic paradigm has proven to be a robust tool in describing both classical and quantum liquids [1][2][3][4]. Recent effort to apply the formalism of hydrodynamics concentrates, among others, on systems with kinematic constraints. A variant of this problem is particularly important in the field of amorphous solids or glasses [5].
Classical glasses refer to any noncrystalline solid that exhibits a glass transition when heated towards the liquid state. Models of classical glasses are largely stochastic lattice models with imposed kinetic constraints on the allowed transitions between different configurations of the system, while preserving equilibrium conditions. The glassy behavior of such systems is visible in the equilibration timescales that grow exponentially with system size. Unfortunately, the constrained models proposed to elucidate classical glasses are only subject to numerical studies and have not been systematically analyzed in terms of the low-energy hydrodynamic theory.
Systems with mobility restrictions inspired a vast amount of work in theoretical physics, e.g. see [34,35] and references therein. Recent developments revealed that theories conserving multipole moments of conserved charges go beyond conventional field theories [36][37][38][39], requiring a different approach for integrating high-energy modes [40,41], and leading to new geometric structures that appear when coupled to background geometries [42][43][44]. Consequently, transport properties of kinematically constrained many-body systems are fundamentally different from theories without mobility restrictions. This is manifested in the thermalization properties of systems with emergent dipole-conservation that exhibit a characteristic subdiffusive behavior [45][46][47][48]. Such systems have been studied both experimentally [31,32] and numerically [45,46]. These developments have lead to theoretical progress aiming at a consistent hydrodynamic theory of fluids with dipole moment conservation, which is sometimes referred to as fracton hydrodynamics.
First, the hydrodynamic theory for simple systems conserving charge and dipole, but not momentum, has been proposed [49]. Subsequently, ideal hydrodynamics for several classes of kinematically constrained fluids with momentum conservation was formulated [50]. Finally, dissipative effects for fracton fluids without energy conservation were studied in Ref. [51] and other works followed thereafter [52][53][54][55].
In addition to filling in some pedagogical gaps, we strive to generalize previous work by providing a systematic treatment of dissipative dipole-conserving fluids with both, momentum and energy conservation, using the method of irreversible thermodynamics. Within our theory, momentum per particle behaves as a Stüeckelberg field, therefore, in analogy with the superfluid, we propose that constant gradients of such a variable are allowed at thermodynamic equilibrium. In fact, the conjugate variable to the gradient of momentum is understood as the flux of dipoles. Our main results are a set of dynamical equations for zeroth (Eqs. (27)) and first order (Eqs. (46)) linear hydrodynamics, respectively. Contrary to the case of ideal ordinary fluids, the zeroth order equations of motion contain a dissipative transport coefficient interpreted as a thermal conductivity. On the other hand, the first order hydrodynamic equations contain 1+12 transport coefficients. Strikingly, the thermal conductivity α modifies the structure of the hydrodynamic modes predicted in [34,51], upgrading the longitudinal modes to an ordinary diffusive form ω || ∼ −iαk 2 , while allowing for subdiffusion only in the shear sector ω shear ∼ −iηk 4 , where η is the analogue to the shear viscosity of the system. In addition, for some regions in the parameter space, the soundlike modes scale as ±k 2 − ik 2 , whereas in the complementary regions they scale as −ik 2 . Therefore, their propagation will be either strongly attenuated or fully diffusive.
The paper is organized as follows. In Sec. II we introduce the symmetries that characterize the kinematically constrained fluids studied in this work. In Sec. III we propose a thermodynamic description compatible with the symmetries. In Sec. IV the hydrodynamic expansion is systematically derived on the basis of the entropy current formalism, the constitutive relations are determined and the hydrodynamic modes are studied. Finally, in Sec. V we conclude with a brief discussion and outlook.

II. SYMMETRIES
Let us start by introducing the set of symmetries underlying our hydrodynamic construction of kinematically constraint fluids.
To this end, we consider a many-body system enjoying a simultaneous conservation of energy H, momentum P i , U(1) charge Q, dipole moment D i and angular momentum J ij . In the long-wavelength regime, the conserved charges can be expressed in terms of the local densities Furthermore, we will impose parity and time-reversal invariance. The long-wavelength dynamics near thermodynamic equilibrium will be governed by the gapless degrees of freedom. Of course, it is only the locally conserved densities that remain relevant, as non-conserved quantities are expected to be fast-relaxing. Therefore, the hydrodynamic equations will be the local conservation laws Note that neither dipole nor angular conservation lead to additional hydrodynamic equations, since their conservation follows from Eqs. (2), provided that the stress tensor T ij is symmetric.
The set of transformations generating the conserved charges form a Lie group with algebra 1 as follows: Notice that the last bracket in Eq. (3) implies, that momentum is not invariant under the transformation generated by D i . In fact, it transforms as Consequently, the momentum density and the stress tensor must transform under dipole shift in the following manner: This transformation property was also obtained in Ref. [44] by placing the system in curved space and coupling it to Aristotelian background sources as well as appropriate gauge fields. As we see in the following sections, these unusual transformation properties lead to a number of exotic features in the hydrodynamics (and thermodynamics) of dipole-conserving fluids.

III. THERMODYNAMICS AND DIPOLE CONSERVATION
Thermodynamics of dipole-conserving systems cannot be captured by the standard textbook treatment. It re-quires a modified approach that systematically incorporates kinematic constraints arising from the dipole conservation. In this section, we construct a consistent thermodynamic theory with dipole conservation built into it.

A. Dipole-invariant equation of state
The internal energy density of a generic system in equilibrium is a function of the entropy and conserved charge densities, for example ≡ (n, s, p i ). However, owing to the noncommutative structure of the algebra (3), dipole-conserving systems are not generic. In fact, the combination p i /n has a shift symmetry under dipole transformations in analogy to a Nambu-Goldstone mode; therefore, it could enter via the invariant combination V ij = ∂ i (n −1 p j ). It is then necessary to introduce a conjugate variable F ij that, as we will see, can be interpreted as a flux of dipoles. Thus, we infer that different (constant) values of V ij label distinct thermodynamic states. For such systems, we postulate that the first law of thermodynamics reads 2 Besides rotational invariance, we have also assumed that microscopically the system preserves parity. After noticing that under parity transformations V [ij] transform as a pseudo-scalar and as a pseudo vector in two and three dimensions respectively, we conclude the energy density must depend on the symmetric part of V ij only. The pressure of the system is then defined via the standard thermodynamic relation From now on we will refer to the symmetrized tensor as V ij . In our construction, we shall take n, , p i as the hydrodynamic variables. Therefore, we find convenient to use entropy density as the thermodynamic potential. Within this picture, where we have defined thermodynamic quantities The relation s ≡ s( , n, V ij ) is then interpreted as the equation of state and the thermodynamic quantities are understood as functions of the dipole-invariant variables ( , n, V ij ).

B. Thermodynamic relations
In the next section, we will study linearized hydrodynamics around the global equilibrium state (n = n 0 , = 0 , V ij = 0). Therefore, it will be useful to introduce a set of thermodynamic identities that will allow us to relate the variations of ( , n, V ij ) with their corresponding conjugate variables.
To do so, we first expand the entropy density function around the equilibrium state up to the second order in fluctuations where s 0 is the entropy evaluated at the equilibrium state, the traceless symmetrization defined as Therefore, the variations of the thermodynamic quantities Eq. (9) can be expressed as Finally, after using Eq. (12) and Eq. (7) we write the variation of the pressure with respect to the thermodynamic variable as where we have defined IV. DIPOLE-CONSERVING HYDRODYNAMICS In this section, we develop the hydrodynamic framework for dipole-conserving fluids by applying the entropy current formalism. We construct constitutive relations following a derivative expansion and impose the second law of thermodynamics. First, we consider the leading order hydrodynamics, and then we study first-order corrections in a linearized regime. After having the most general constitutive relations, the transport coefficients will be constrained by the entropy production condition.

A. Gradient expansion
Following the canonical paradigm of hydrodynamics, we consider the long-wavelength, near-equilibrium dynamics that is governed by the hydrodynamic variables, that is, the densities n, , p i of the conserved charges. Macroscopic currents are then given by local expressions of the conserved densities organized in a systematic derivative expansion. The explicit form of the currents dubbed as constitutive relations is fixed by the symmetries in Eqs. (1). In writing these constitutive relations a set of unknown parameters will emerge, known as transport coefficients, which are then constrained imposing the laws of thermodynamics, and Onsager relations.
Nonetheless, the non-standard structure the dipole symmetry introduces suggest we should consider the momentum of the system p i to be of order Therefore, our derivative expansion is defined in terms of the order at which the equations of motion are truncated, e.g. we will refer to n−th order hydrodynamics if the set of differential equations is truncated as 3

B. Zeroth order hydrodynamics
To start, we derive the leading order constitutive relations for dipole conserving fluids. Our results are a dissipative completion of the theory constructed in Ref. [34], where the zero temperature ideal constitutive relations were found using the Poisson bracket formalism.
The local form of the first law can be derived from Eq. (8) as follows in the last step we have defined the effective chemical potential and velocitỹ Using the equations of motion Eqs. (2) it is then possible to recast Eq. (16) into a familiar looking equation where we have defined a shifted energy current Throughout the rest of this work, we will often find it convenient to work with the shifted energy current. Using Eq. (18) we can express the entropy production constraint ∂ t s + ∂ i S i ≥ 0 as Thus, after combining the thermodynamic relation Eq. (7) with Eq. (20) and a series of tedious algebraic computations that we show in Appendix A 1, it is possible to rewrite the constraint as Therefore, we conclude that the local version of the second law of thermodynamics will be satisfied provided that the first term in Eq. (21) vanishes and the remainder is semi-positive definite for arbitrary field configurations.
This constraint fixes the zeroth order currents to with α a transport coefficient that can be interpreted as the thermal conductivity of the system, satisfying the inequality α ≥ 0 .
Actually notice that the stress tensor has the required transformation property under dipole shift Eq. (5). In addition, it can be shown that the entropy current reduces to the simple form It is important to emphasize that both of the terms in the entropy current enter with a single spatial derivative of a hydrodynamic variable. Thus, in dipole conserving systems the limit of ideal hydrodynamics, as constructed in Ref. [50], can only be reached by fine-tuning α = 0 rather than neglecting higher order derivative corrections. This happens because the lowest order contributions, in our counting scheme, allow for a dissipative transport coefficient. This is at odds with fluids without kinematic constraints.
We now turn our attention to the study of the hydrodynamic modes. To this aim, we consider linearized perturbations around the equilibrium state ( 0 , n 0 , p 0 = 0). Therefore, the fluctuations read n = n 0 + δn, and the corresponding currents (22) take the form In order to solve for the evolution of the hydrodynamic variables, we express all quantities that appear in the above currents in terms of the variations of the conserved densities δn, δ and δp. This is done using the thermodynamic relations Eqs. (12) and (13), as well as the definition of the effective velocity Eq. (17). After doing so, the zeroth order hydrodynamic equations of motion become . In order to solve the system and find the dispersion relation of the modes, we must Fourier transform the equations with frequencies and momenta (ω, k i ). For this set of equations, the transverse sector (k·δp = 0) contains the non-dispersive mode ω shear = 0.
The solutions to Eq. (28) are with Notice that all modes have a k 2 dependence. Actually, it is possible to expand the solutions for small and large value of the thermal conductivity. In the former case, which corresponds to a 2 1 and a 1 /a 2 fixed, the modes read Whereas in the latter (a 2 1 and a 1 /a 2 fixed) the asymptotic behaviors for of the dispersion relations are The full dependence of the modes as a function of the adimensional thermal conductivity a 2 at fixed a 1 /a 2 is shown in Fig. 1. We notice the existence of two qualitatively distinct regimes corresponding to a 1 /a 2 < 1/9 (Regime I) and a 1 /a 2 > 1/9 (Regime II). Actually, at the critical regime (a 1 /a 2 ) c = 1/9 there is a point for which the three modes are equal and purely imaginary. This situation is shown in the middle plot in Fig. 1. The main difference between Regimes I and II is the existence of a window in the parameter space where the three modes are purely imaginary (Regime I), whereas in Regime II there will always be two modes with a non-vanishing real part.
Finally, it is worth to point out that the thermodynamic constraints a 1 /a 2 < 1 and α ≥ 0 are enough to guarantee that the imaginary part of the modes is negative, and no linear instability will occur.

C. First order hydrodynamics
In the previous section, we have shown the existence of one dissipative transport coefficient α that controls how longitudinal fluctuations diffuse in the system. However, the shear mode ω shear = 0 remained insensitive to the thermal conductivity. Therefore, at that level of the derivative expansion, transverse fluctuations will not diffuse. Although we may think this fact is reminiscent of the fractonic nature of the system, in this section we will prove that this is not the case. Actually, the first order transport coefficients will introduce transverse contributions to the next to leading order hydrodynamic equations of motion and predict a subdiffusive shear mode.
To proceed with the analysis, we decompose the currents into the zeroth and first order contributions according to the derivative counting scheme and plug the decomposition Eq. (34) into Eq. (16) and cancel out the lower order terms (as these satisfy the second law provided that α ≥ 0). Then, the second law requires that Note that for fluctuations around Eq. (25) we have that µ = µ and J i = E i (see Eqs. (17) and (19)). Since our goal is to identify the most general constitutive relations consistent with the above inequality, we find it convenient to rewrite Eq. (35) as follows Ignoring non-linearities, it is possible to express the energy current without loss of generality as E i 1 = ∂ j E ji . With that ansatz at hand, the entropy production constraint can be written as (37) Then, we set the first order correction to the entropy current to be such that the first term in Eq. (37) vanishes. Therefore, the problem is reduced to finding the proper constitutive relations such that the leftover is semi-positive definite. In fact, the most general form for the currents that will allow a positive production read Where we have introduced a set of 18 dissipative transport coefficients. In particular, Onsager reciprocity re-duces the number of off-diagonal coefficients if timereversal invariance is imposed Moreover, the entropy production constraint Eq. (37) can be written in a compact matrix form with Since the two contributions are independent, the second law is then imposed by requiring that matrices A || and A ⊥ are both semi-positive definite. This poses constraints on the transport coefficients, which are summarized below.
In total, we have found 12 independent transport coefficients that we classify in two distinct categories. The first category involves the diagonal coefficients (ζ, η, σ ⊥ , σ || , κ ⊥ , κ || ) satisfying a positivity constraint On the other hand, the second category consists of the off-diagonal terms (α ⊥ , α || , β ⊥ , β || , γ ⊥ , γ || ) obeying inequalities with the coefficients of the previous group The distinction has been motivated by the fact that the value of the off-diagonal transport coefficients is bounded from above by the diagonal ones. The last constraint from semi-positivity corresponds with the positive determinant condition Having the corrections to the zeroth order hydrodynamic constitutive currents, we can plug them into the conservation equations to obtain the first order hydrodynamic equations of motion. In particular, they read ∂ t + (αs ee + e e ∇ 2 )∇ 2 δ + (αs ne + e n ∇ 2 )∇ 2 δn For a detailed derivation of the equations and the relation of the parameters shown in Eqs. (46) with the transport coefficients in Eqs. (39), we refer the reader to the Appendix B. The main output of the first order approach is the conversion of the non-dispersive shear mode into a subdiffusive one On the other hand, the longitudinal modes are not strongly affected by the first order corrections, since their contribution enters at higher order in momentum. In fact, the characteristic polynomial in this case takes the same form as in Eq. (28) with b 0 , b 1 and b 2 derived in Appendix B, and shown in Eq. (B10). Actually, the solutions Eqs. (30) still apply, once we substitute a 1 →ā 1 and a 2 →ā 2 in the Eqs. (31).

V. DISCUSSION
We have presented a hydrodynamic theory of isotropic dipole-conserving fluids up to the first order in the derivative expansion. Our construction gives universal lessons about systems with constrained dynamics.
First, we have shown how to consistently implement the kinematic constraints of dipole-conserving fluids at the level of equilibrium thermodynamics. This thermodynamic state is not the same as for conventional fluids, but requires an additional tensor quantity controlling the flux of dipoles. Secondly, we elucidated the derivative expansion around this equilibrium state and constructed hydrodynamic constitutive relations. Peculiarly, we have found that dipole-conserving fluids are dissipative even if the equations of motion are truncated at the lowest non-trivial order in the derivative expansion.
Our theory gives a conceptually crisp finite temperature description of systems, which preserve charge and its dipole moment, energy and momentum. As a result, it completes previous studies that have given partial answers to this problem for systems without momentum conservation [49], without energy conservation [51] or without dissipative effects [50].
Using the linear response theory, we have found and classified new 1 + 12 transport coefficients that can, in principle, be used as experimental signatures of fracton fluids. In particular, we find non-trivial corrections to both, the transverse and longitudinal modes studied in Ref. [50] for ideal fracton fluids. One dissipative coefficient, which we identify with the thermal conductivity, contributes to the diffusive transport of longitudinal perturbations. The shear mode, which is non-dispersive at the lowest order in derivatives, becomes subdiffusive after the inclusion of the first order corrections.
The generalization of our formalism to generic multipole-conserving systems is straightforward. In fact, we should expect the real part of the longitudinal modes and the shear mode to acquire a higher exponent in momentum. However, the imaginary part of the longitudinal modes should still have the −ik 2 scaling due to its link to the thermal conductivity. Although we have proven a consistent thermodynamic description of the system, its origin from a microscopic perspective is obscure since dxV ij does not correspond to any conserved operator. Therefore, a statistical mechanics picture of the thermodynamic theory proposed here is still an open problem.
Finally, let us comment on the non-equilibrium universality classes recently proposed in [49,51]. In particular, the authors of Ref. [51] have shown that dipoleconserving fluids at rest become unstable against thermal fluctuations in less than four spatial dimensions. However, the dipole-conserving models considered there do not assume energy conservation and consequently display subdiffusive scaling of the hydrodynamic modes. Thus, we suspect that including energy conservation changes the universality class and will render the system stable against thermal fluctuations in three spatial dimensions. It would be interesting to investigate this point in greater detail. In this Appendix we show the algebraic manipulations that were used when going from Eq. (20) to Eq. (21) explicitly.
We start by rearranging Eq. (20) as a total derivative minus the extra terms Driven by intuition from the conventional hydrodynamics, we now incorporate pressure into our construction by adding zero 0 = −∂ i (P Vi T ) + P ∂ i Vi T + Vi T ∂ i P such that (A1) becomes Using the definition of pressure Eq. (7) it is possible to evaluate the gradient We substitute µ =μ + n −1 V i p i and rewrite the last term in Eq. (A3) obtaining Let us now express the last term as a total derivative Using the definition of velocity V i = −n −1 ∂ j F ji and relabelling the dummy indices k ↔ j in the second-to-last term we find Thus, the last term in Eq. (A2) can be written as follows Relabelling the dummy indices i ↔ j in the last two terms and using the definition of pressure Eq. (7) we arrive at Plugging (A8) back into (A2) we obtain an expression that closely resembles Eq. (21) It is then straightforward to reach (21) by adding another zero This final step is necessary in order to obtain a stress tensor that is manifestly symmetric under the exchange of indices.