Viscoelastic hydrodynamics and holography

We formulate the theory of nonlinear viscoelastic hydrodynamics of anisotropic crystals in terms of dynamical Goldstone scalars of spontaneously broken translational symmetries, under the assumption of homogeneous lattices and absence of plastic deformations. We reformulate classical elasticity effective field theory using surface calculus in which the Goldstone scalars naturally define the position of higher-dimensional crystal cores, covering both elastic and smectic crystal phases. We systematically incorporate all dissipative effects in viscoelastic hydrodynamics at first order in a long-wavelength expansion and study the resulting rheology equations. In the process, we find the necessary conditions for equilibrium states of viscoelastic materials. In the linear regime and for isotropic crystals, the theory includes the description of Kelvin-Voigt materials. Furthermore, we provide an entirely equivalent description of viscoelastic hydrodynamics as a novel theory of higher-form superfluids in arbitrary dimensions where the Goldstone scalars of partially broken generalised global symmetries play an essential role. An exact map between the two formulations of viscoelastic hydrodynamics is given. Finally, we study holographic models dual to both these formulations and map them one-to-one via a careful analysis of boundary conditions. We propose a new simple holographic model of viscoelastic hydrodynamics by adopting an alternative quantisation for the scalar fields.


| Introduction
When undergoing deformations, most observable materials are known to exhibit both elastic and viscous responses. Due to the coupling between fluid and elastic behaviour, such materials are said to be viscoelastic. Despite being a century old subject and an active research field with multiple technological applications [1], the understanding of viscoelasticity has been mostly based on phenomenological models that assume linear strain responses such as the Kelvin-Voigt model, Maxwell model and the Zener model as well as on some nonlinear generalisations thereof (see e.g. [2]). Several efforts have been made in order to formulate viscoelasticity from general principles. In particular, the work of Eckart [3] was the key towards the geometrisation of strain and the introduction of the notion of a dynamical material reference state. More recently, works based on non-equilibrium thermodynamics have brought some of these aspects to covariant form in the non-relativistic context [4,5] and in the relativistic context [6][7][8] and recovered the stresses and rheology of a few of the above mentioned viscoelastic models. However, while significant, these works have not characterised the full interplay between fluid and elastic behaviour due to defining assumptions. The aim of this paper is to provide such a full characterisation, under relaxable conditions, in the hydrodynamic regime.
According to Maxwell, the defining property of viscoelasticity is the capacity for continuous media to exhibit elasticity at short time scales and fluidity at long time scales compared to the strain relaxation time [9,10]. If the relaxation time is very large, fluid and elastic behaviour can coexist in the hydrodynamic limit. This is the realm of (liquid) crystal theory. We are thus interested in a long-wavelength long-distance effective description of crystals. A crystal is characterised by a regularly ordered lattice of points (atoms or molecules) discretely distributed over space. More generally, the lattice cores that constitute the crystal can be higher dimensional, such as strings and surfaces, where the atoms/molecules have no positional ordering within the cores and can move freely like a "liquid". These crystals are called liquid crystals. Crystals may be present in different phases, such as elastic (solid) phase, smectic or nematic, among others (see e.g. [11,12]). In the non-relativistic context, the hydrodynamics of (liquid) crystals has been considered in several works [11,13,14] but these treatments assume isotropy, no external currents and do not explicitly derive the constitutive relations and stress/strain relations that couple the fluid and elastic degrees of freedom. In this paper we focus on describing the elastic and smectic phases within a modern framework of hydrodynamics, which includes effective field theory [15], offshell adiabatic analysis [16,17], and hydrostatic partition functions [18,19].
to formulate a hydrodynamic description of anisotropic (liquid) crystals with nonlinear strains of arbitrary strength under the assumptions of (i) absence of dislocations and disclinations, (ii) lattice homogeneity (i.e. invariant under φ I → φ I + a I with a I being a constant translation), and (iii) non-dynamical intrinsic crystal metric (i.e. the crystal reference state does not change in time and thus we do not consider plastic deformations). Within this formulation, we describe the structure of first order elastic responses and transport properties of viscoelastic fluids in a long-wavelength expansion. In the elastic phase of isotropic crystals and under the assumption of linear strain responses, we uncover 5 extra transport coefficients that have not been considered in the literature.
Recently, it has been argued that viscoelastic hydrodynamics can be recast as a theory of higherform hydrodynamics making its global symmetries manifest and avoiding the need to introduce microscopic dynamical fields [22]. This follows the recent line of research where hydrodynamic systems with dynamical fields, such as magnetohydrodynamics with dynamical gauge fields, are recast in terms of dual hydrodynamic systems with higher-form symmetries [23][24][25][26][27][28]. Building up on this idea, here we provide a completely equivalent description of viscoelastic hydrodynamics of anisotropic (liquid) crystals in arbitrary spacetime dimensions by identifying the correct degrees of freedom of higher-form hydrodynamics. The resulting theory describes higher-form superfluidity in which the higher-form symmetries are partially broken (as in the context of magnetohydrodynamics [26,27]). The usefulness of this formulation resides in the fact that the hydrodynamic description can be recast entirely in terms of symmetry principles where the trivial conservation of the higher-form current J I = dφ I , with being the Hodge operator in d + 1 dimensions, is a consequence of the absence of defects (e.g. disclinations or dislocations). We provide an exact map between the two formulations.
The lack of control over viscoelastic theories and the absence of complete formulations have prompt a series of works where holographic methods are employed to study putative strongly coupled viscoelastic theories and probe regimes of elasticity and fluidity (e.g. [22,[29][30][31][32][33][34][35][36]). Two types of models have been considered in the literature: gravity coupled to a set of scalar fields Φ I [30,31,33] (and with additional fields [31,32]) and gravity minimally coupled to a set of higher-form gauge fields B I [22]. The former is supposed to describe the dynamics of viscoelastic materials with spontaneously broken translation symmetries while the latter is supposed to describe viscoelastic theories with higher-form currents (see also [37,38]). However, the establishment of a precise map between the two hydrodynamic formulations has prompt us to investigate whether such a map exists at the level of holographic models. Indeed, a careful analysis of boundary conditions has led us to propose a simple model of viscoelasticity consisting of a set of scalars Φ I minimally coupled to gravity but with an alternative quantisation for the scalar fields and a double trace deformation of the boundary theory. The model is thus the linear axion model of [39] (see also e.g. [40][41][42]) used in the context of momentum relaxation [43,44] but which does not treat the fields Φ I as background fields (as in the setting of forced fluid dynamics [45,46]). Instead, the scalar fields are dynamical fields, as in a viscoelastic theory where the Goldstone scalars have inherent dynamics.
Summarising, in this paper we make the following advancements and solve the following problems/issues: • We provide a complete formulation of nonlinear viscoelastic hydrodynamics of anisotropic (liquid) crystals in terms of Goldstone scalars of spontaneously broken symmetries up to first order in a long-wavelength expansion. We derive the Josephson equations for the Goldstone modes, akin to that found in the context of superfluids [47]. We provide a classification of the response and transport of linear isotropic materials and recover Kelvin-Voigt materials as a 2. Broken translations and elasticity | 5 special case.
• We formulate the same theory of nonlinear viscoelastic hydrodynamic in terms of a higher-form superfluid by identifying the correct hydrodynamic degrees of freedom. This formulation, based on generalised global symmetries, provides an organisation principle and a first principle derivation of viscoelastic hydrodynamics that does not involve additional microscopic dynamical fields.
• We provide holographic models for both these formulations in D = 4, 5 bulk spacetime dimensions as well as a map between the two models corresponding to each of the formulations. We identify the boundary conditions and boundary action necessary for obtaining holographic viscoelastic dynamics with the simple model of [39].
This paper is organised as follows. In section 2 we review the classical elasticity effective field theory in terms of the Goldstones of broken translational symmetries. However, we reformulate it in terms of surface calculus, which besides aiding in understanding the appropriate degrees of freedom of viscoelasticity, also leads to a precise covariant geometrization of elastic strain. In section 3 we formulate viscoelastic hydrodynamics in terms of Goldstone scalars up to first order in a derivative expansion. We obtain the Josephson conditions and construct a hydrostatic effective action that characterises the equilibrium viscoelastic states in the theory. We also study the rheology equations and comment on phenomenological viscoelastic models. In section 4 we formulate the same theory as a novel theory of higher-form superfluidity, generalising [25][26][27] to arbitrary d-dimensional higher-form currents. In this section, we also provide a detailed map between the two formulations. Section 5 is devoted to the construction of holographic models of viscoelastic dynamics and provides appropriate holographic renormalisation procedures. It also contains a study of conformal viscoelastic fluids. In sec. 6 we conclude with a summary of the results obtained in this paper together with interesting future research directions. We also provide further details on the geometry of crystals in appendix A, while in appendix B we give the details of hydrostatic constitutive relations. Finally, in app. C we provide precise comparisons between our different formulations and earlier ones in the literature.

| Broken translations and elasticity
Utilising elements from earlier formulations (e.g. [20,21,48]), where Goldstones of spontaneously broken translational symmetries play a key role, we introduce a classical effective field theory for crystals exhibiting solid and smectic phases. As mentioned in the introduction, crystals arrange themselves into a structured lattice of points, strings, or surfaces (generically called lattice cores). In order to deal with this wide range of higher-dimensional objects, we present a new reformulation of classical elasticity effective field theory in terms of surface calculus, which proves to be useful in later sections for tackling the hydrodynamic regime of liquid crystals. In particular, this formulation provides a simple and covariant notion of strain and allows us to cover solid and smectic phases simultaneously. We begin with zero temperature considerations, moving on to finite temperature effects towards the end of this section. We consider (d + 1)-dimensional spacetimes where d is the number of spatial dimensions. In the continuum limit, valid at long distances and low energies, the worldsheets of (d − k)-dimensional crystal cores can be parametrised by a set of k spacetime dependent one-forms e I µ (x) normal to the cores, with I = 1, 2, . . . k ≤ d. Point-like cores correspond to k = d, string-like to k = d − 1, and so on. In general, these normal one-forms have an inherent spacetime dependent GL(k) ambiguity due to arbitrary normalisation: e I µ → M I J e J µ with M I J ∈ GL(k). We keep this redundancy unfixed for now by allowing for a local GL(k) symmetry in the effective theory.
Given that the background spacetime is equipped with a metric g µν (which can be set to η µν = diag(−1, 1, 1, . . .) for crystals in flat space), the physical distance between the cores is determined using the crystal metric The metric h IJ is the transverse metric to the crystal cores, obtained by projecting the spacetime metric along the normal one-forms. The indices I, J, . . . can be raised/lowered using h IJ and h IJ . For later convenience, we also define a pair of spacetime projectors transverse and along the crystal cores by pushing forward the crystal metric whereh µν is the longitudinal projector and h µν the transverse projector. On the other hand, the crystal also carries an intrinsic reference metric that captures the lattice structure of the crystal and determines the "preferred" distance between the cores when no external factors are at play. We define this reference metric as where h IJ is an arbitrary non-singular symmetric matrix. The difference between the two metric tensors on the crystal defines the strain tensor which captures distortions of the crystal away from its reference configuration. Subjecting a crystal to a strain, i.e. distorting the crystal, causes stress depending on the physical and chemical properties of the material that constitute the crystal. Within an effective field theory framework, we will attempt to characterise the most generic such responses, given the symmetries of the crystal.

Crystal fields
It is a known result in differential geometry of surfaces that a generic set of one-forms e I µ does not have to be surface forming, i.e. there might not exist a foliation of crystal core worldsheets normal to all the e I µ . For this to be the case, one needs to invoke the Frobenius theorem, ensuring that there must exist a set of spacetime one-forms a I νJ such that ∂ [µ e I ν] = −a I [µJ e J ν] . As a consequence, the variations of the one-forms e I µ are not independent and we cannot use them or the strain tensor u IJ directly as fundamental degrees of freedom in the effective theory of crystals. To get around this nuisance, we assume that the normal one-forms can locally be spanned by a set of k closed one-forms, is an arbitrary invertible matrix and φ I (x) are possibly multi-valued smooth scalar fields. This choice corresponds to the crystal core worldsheets being level surfaces of the functions φ I (x) and satisfies the Frobenius condition as . The fields φ I (x), which we refer to as crystal fields, describe the position of the crystal structure in the ambient spacetime. These crystal fields can be physically understood as Goldstones of spontaneously broken translations. If the crystal does not have any topological defects such as dislocations or disclinations, the fields φ I (x) can be taken to be single-valued and well-behaved (see e.g. [49]).
Recall that we had an arbitrary GL(k) renormalisation freedom in e I µ , which we can now fix by setting Λ I J (x) to the identity matrix. Consequently, 5) and the physical and reference metrics of the crystal can be expressed as We should note that as such, like in any field theory, there is an arbitrary redefinition freedom in the choice of the fundamental crystal fields φ I . This will be useful later.

Plasticity, homogeneity, and isotropy
Generically, the reference metric h IJ (x) of the crystal is a dynamical field and can evolve independently with time. This physically describes "plastic materials" for which the applied strain can permanently deform the internal structure of the crystal over time. In this work, however, we will focus on "elastic materials", wherein we assume that h IJ (x) = h IJ (φ(x)), i.e. the reference metric is an intrinsic property of the crystalline structure and is not dependent on a particular embedding of the crystal into the spacetime. The functional form of h IJ (φ) is a property of the physical system under observation and needs to be provided as input into the theory.
Furthermore, the crystals we wish to describe using this effective field theory are homogeneous in space at macroscopic scales. Therefore, there exists a choice of crystal fields φ I such that the reference metric h IJ is constant and the theory is invariant under a constant shift φ I → φ I + a I . In fact, for homogeneous crystals, we can utilise the φ I redefinition freedom to set the reference metric to be the Kronecker delta, that is This leaves just a global SO(k) rotation freedom among φ I , i.e. φ I → Ω I J φ J where Ω I J is a constant matrix valued in SO(k). As long as we properly contract the I, J, . . . indices, we do not need to worry about this redundancy while constructing the effective field theory.
Finally, if we wish to describe a crystal that is isotropic at macroscopic scales (possibly due to randomly oriented crystal domains), we can impose the aforementioned global SO(k) freedom of φ I as an invariance of the theory. Along with the constant shift invariance due to homogeneity, this results in a Poincaré invariance on the field space. Practically, it means that besides e I µ and K ext I , the field space indices I, J, . . . in the theory can only enter via the reference metric δ IJ . If, instead, the crystal under consideration has long-range order, the parameters of the effective theory can be arbitrary rank tensors on the field space. 1 We provide further details on the geometry of crystals in appendix A. In the bulk of this work we will assume the crystals being described to be "elastic" and "homogeneous", while no assumption is made regarding isotropy except in some explicit examples.

Elasticity at zero temperature
Previously we have introduced the geometric notions required to describe crystals but we have not yet attributed dynamics to the crystal fields. Here we consider classical elasticity field theory at zero temperature, for which the dynamics follows from an action principle that is written in terms of the appropriate crystal fields (determined earlier to be φ I ).

Effective action
We posit that our theory of interest is described by an effective action with functional form S[φ I ; g µν ], where φ I are the dynamical Goldstones of broken translations and g µν is taken to be a background metric field. We focus on homogeneous crystals, for which the action is invariant under a constant translation of the crystal fields φ I (x) → φ I (x) + a I and all the dependence on φ I appears via its derivatives e I µ = ∂ µ φ I . The action can then be parametrised as We define the crystal momentum currents by varying the action with respect to e I µ = ∂ µ φ I , that is (2.9) Given homogeneity, the equations of motion for φ I simply imply the conservation of crystal momentum currents ∇ µ σ µ I + K ext where K ext I is a background field, which can be understood as an external force sourcing the crystal fields. 2 The conservation eq. (2.10) is not protected by any fundamental symmetry and will in general be violated by thermal corrections as we will see in the next section. If we further assume that all the dependence on e I µ comes via h IJ or equivalently u IJ , the crystal currents can also be obtained by varying the action with respect to the strain tensor 3 (2.11) 1 This structure only pertains to the geometric structure of the crystal itself. In general the atoms/molecules occupying the lattice sites can also carry other preferred vectors like spin or dipole moment which will need to be considered independently. 2  Finally, we can obtain the energy-momentum tensor of the theory by varying the action with respect to the background metric Given that the action as constructed is invariant under background diffeomorphisms, the energymomentum tensor is conserved, modulo background sources 4 (2.13)

Linear isotropic materials at zero temperature
As an illustrative example, we consider classical elasticity theory where all the dependence of e I µ comes via the strain h IJ , and expand the Lagrangian in a small strain expansion. In the case of homogeneous crystals, the strain is given by u IJ = 1 2 (h IJ − δ IJ ) and the most generic such effective action for an isotropic crystal at zero-derivative order and quadratic in strain is given by 5 (2.14) Here C IJKL is the elasticity tensor of the crystal 15) and the coefficients B and G are the bulk modulus and shear modulus of the crystal respectively, whereas P does not have a standard physical interpretation in the literature. In eq. (2.15), we have chosen to express the coefficients using h IJ for convenience, but we could equivalently have used δ IJ . This choice leads to the same physical currents up to linear order in strain. By varying the action with respect to h IJ , we can read out the crystal momenta On the other hand, the energy-momentum tensor is given by The bulk modulus B and the shear modulus G couple to the trace and traceless parts of strain, respectively, in the energy-momentum tensor. The coefficient P, on the other hand, gives a constant pressure contribution along the field directions in the energy-momentum tensor, modelling a repulsion between lattice points. Such crystals cannot be supported without non-trivial boundary conditions on their surface. For most phenomenological applications, the lattice points are effectively neutral and the coefficient P can be dropped. However, as we will see in section 5, this coefficient appears naturally in holographic models of elasticity. 4 At zero derivative order, where all the dependence on the metric in the Lagrangian comes via hIJ as well, the energy momentum tensor reads T µν = L g µν + σIJ e Iµ e Jν . This leads to the well known definition of stress tensor as the conjugate to strain T IJ = T µν e I µ e J ν = L h IJ + ∂L/∂uIJ (see for example section 6.3.3 of [11]). 5 Note that 1 2 log det(hIJ ) = 1 2 log det(δIJ + 2uIJ ) = h IJ uIJ + h I(K h L)J uIJ uKL + O(u 2 ). So, generically, the P term here can also be replaced by the term linear in strain h IJ uIJ by redefining B and G.

Heating up the crystals
So far we have focused on the effective field theory describing crystals at zero temperature. However, for the phenomenological applications that we have in mind, we need to take into account the effects of finite temperature. In this section we discuss crystals in thermodynamic equilibrium using the Matsubara formalism of finite temperature field theory and introduce the equilibrium effective action. Towards the end we motivate the hydrodynamic formulation of crystals seen as small dynamical perturbations around thermodynamic equilibrium, which we later elaborate in section 3.

Equilibrium effective action
The fundamental entity of interest at finite temperature is the thermal partition function written in a given statistical ensemble. However, our understanding of a complete partition function describing arbitrary non-equilibrium thermal processes in a quantum field theory is still very limited. Nevertheless, if we focus on just equilibrium (time-independent states) in the theory, the grand canonical partition function can be computed using the Matsubara imaginary time formalism Here S eqb is the equilibrium effective action of the theory that is, naively, obtained by Wick rotating the Lorentzian action S. It should be noted that defining the partition function above requires us to pick a preferred time coordinate with respect to which the equilibrium is defined and with respect to which the Wick rotation is to be performed. Consequently, in an effective field theory approach, the equilibrium effective action S eqb can contain many new terms dependent on the preferred timelike vector that have no analogue in the original zero temperature effective action S. To make this precise, let us define K µ = δ µ t /T 0 to be the preferred timelike vector, with T 0 being the inverse radius of the Euclidean time circle interpreted as the global temperature of the thermal state under consideration. The requirement of equilibrium implies that the Lie derivative of the constituent fields g µν and φ I along K µ is zero, leading to The equilibrium effective action and the resulting thermal partition function for a crystal can be schematically represented as (2.20) The integral in the first line is performed over a constant-time Cauchy slice Σ with the respective differential volume-element denoted by dσ µ . The free-energy current (N µ ) eqb is conserved rendering the effective action independent of the choice of Cauchy slice. The finite temperature action, for instance, can have dependence on the scalar K 2 = K µ K ν g µν = g tt /T 2 0 , which has no analogue in the zero temperature effective action. This scalar is related to the local observable temperature T eqb in the field theory (as opposed to the global thermodynamic temperature T 0 ) as Once the effective action is at hand, we can work out the finite temperature version of the φ I equations of motion (2.10) with the crystal momentum currents (2.23) We can also read out the energy-momentum tensor of the theory in thermal equilibrium to be which satisfies the conservation equation (2.13) owing to the background diffeomorphism invariance of the equilibrium effective action.

Linear isotropic materials at thermal equilibrium
Focusing on the model of linear elasticity from section 2.2.2, it is possible to heat it up to finite temperature while keeping it in equilibrium. At zero-derivative order, the equilibrium effective action has a form similar to eq. (2.14), except that here we also need to take into account the dependence on T eqb . We find where Here P f is interpreted as the thermodynamic pressure of the crystal, which is purely a finite temperature effect. Note that at finite temperature, the elastic modulii B and G of the crystal as well as the crystal pressure P are functions of the local temperature. Varying the resulting effective action, we can read out which enter (2.10) and (2.13) to give the φ I equations of motion and energy-momentum conservation in thermal equilibrium respectively. These can be directly compared to their zero temperature counterparts in section 2.2.2. The crystal momentum currents remain similar in form except for the temperature dependence of the coefficients, while the energy-momentum tensor has a few novel terms. The first of these terms corresponds to the thermodynamic energy density while the second to the thermodynamic pressure, as promised.

Leaving equilibrium -hydrodynamics
Although generic non-equilibrium processes in a thermal field theory are not accessible with the machinery at hand, we can leave equilibrium perturbatively using the framework of hydrodynamics. The basic premise of hydrodynamics is that we can describe slight departures from thermal equilibrium by replacing the isometry K µ with a slowly varying dynamical field β µ . The time-evolution of these fields is governed by the energy-momentum conservation (2.13), which in the out of equilibrium context is not trivially satisfied as a mathematical identity. It is customary to isolate the normalisation piece and re-express β µ as Here u µ is the fluid velocity and T is the fluid temperature out of equilibrium. Note that in equilibrium, obtained by setting β µ = K µ , the temperature T reverts back to its equilibrium value T eqb in eq. (2.22), while the fluid velocity u µ is just a unit vector along δ µ t describing a fluid at rest. Out of equilibrium, we no longer have the luxury to derive the φ I equation of motion or the conserved energy-momentum tensor using an effective action. Instead we assume the existence of these as the starting point of hydrodynamics. To wit Note that we have replaced ∇ µ σ µ I from eq. (2.10) by an arbitrary operator K I out of equilibrium, making contact with our previous comment that there is no fundamental symmetry at play to enforce the φ I equations of motion to take the form of a conservation law. As such, K I and T µν can be arbitrary operators made out of the constituent fields in the theory. However, the existence of a partition function in thermal equilibrium implies that the theory must admit a free energy current N µ which reduces to the conserved free-energy current (N µ ) eqb in equilibrium (upon setting β µ = K µ ). Performing a time-dependent deformation of the equilibrium effective action in eq. (2.20), it is not hard to convince oneself that where ∆ is at least quadratic in δ B . Here δ B is defined similar to eq. (2.19) and denotes Lie derivatives along β µ . This is commonly referred to as the adiabaticity equation and determines the allowed terms in K I and T µν in agreement with the equilibrium partition function.
Generally, hydrodynamic systems are required to satisfy a slightly stronger constraint: the second law of thermodynamics. It is possible to define an entropy current S µ = N µ − T µν β ν , which upon using eq. (2.30) and eq. (2.29) satisfies (2.31) The second law requires that the divergence of the entropy current should be locally positive semi-definite, forcing ∆ in eq. (2.30) to be a positive semi-definite quadratic form.
Having motivated hydrodynamics of crystals from the viewpoint of thermal field theories, in the next section we will reintroduce hydrodynamics as its own framework based on the second law of thermodynamics. We will revisit the hydrodynamic elements discussed here and work out the equations governing a crystal in the hydrodynamic regime up to first order in derivatives in agreement with the second law.

| Viscoelastic hydrodynamics
In this section we formulate viscoelastic hydrodynamics as a theory of viscous fluids with broken translation invariance, analogous to [13,14]. This is done by introducing the set of crystal (scalar) fields, one for each spatial dimension along the crystal, as in the previous section, which can be seen as Goldstones of broken momenta. Contrary to previously studied cases of forced fluid dynamics [45] and models of momentum relaxation [44] where the scalar fields are background fields, these Goldstone fields are dynamical. Their dynamics is governed by a Josephson-type condition similar to that encountered in the context of superfluids. We formulate viscoelastic fluids with one-derivative corrections in arbitrary dimensions and study carefully the case of isotropic crystals with linear responses in strain. Attention is given to the resulting rheology equations and a linearised fluctuation analysis is carried out, identifying dispersion relations for phonons and sound modes.

The setup
As discussed in detail in the previous section, a crystal can be characterised by a set of normal one-forms e I µ , with I = 1, 2, . . . k. In the hydrodynamic regime, the dynamics of a viscoelastic crystal is governed by the conservation of energy-momentum tensor along with the "no topological defect" constraint that requires the normal one-forms describing the crystal to be closed and the crystal evolution equations The conservation of the energy-momentum tensor T µν is being sourced by the external sources K ext I coupled to e I µ . It governs the time evolution of a set of hydrodynamic fields: fluid velocity u µ (normalised such that u µ u µ = −1) and fluid temperature T . The constraint in eq. (3.1b) can be identically solved by introducing a set of crystal fields φ I such that e I µ = ∂ µ φ I . Physically, the crystal fields can be understood as Goldstones of broken momentum generators. 6 The dynamics of these crystal fields themselves is governed by eq. (3.1c), where K I is an effective macroscopic operator composed of the field content of the theory. A priori, we do not have any knowledge of the form of this operator. However, much like T µν , within the hydrodynamic derivative expansion, constitutive relations for K I can be fixed using the second law of thermodynamics [47]. Eq. (3.1c) is the finite temperature counterpart of eq. (2.10) but since effective actions for viscoelastic fluids describing dissipative dynamics have not yet been constructed, there is no first principle derivation of K I . 7 6 A closely related formulation of viscoelastic fluids is found in [7], which models viscoelasticity as a sigma model given by d + 1 scalar fields, seen as coordinates on an internal worldsheet. We provide a discussion of the similarities and distinctions between the two formulations in appendix C.1. 7 A nice parallel can be made with the theory of magnetohydrodynamics where the crystal fields φ I are replaced by the photon Aµ and the normal one-forms e I µ by the field strength Fµν . The three equations in (3.1) find their respective analogues in energy-momentum conservation ∇µT µν = −F νρ J ext ρ , Bianchi identity ∂ [µ F νρ] = 0, and Maxwell's equations J µ + J µ ext = 0. See [26,27] for more details.
Hydrodynamics is an effective theory where the most generic constitutive relations for T µν and K I are obtained order-by-order in a derivative expansion in terms of the constituent fields u µ , T , φ I and background field g µν . These constitutive relations are required to satisfy the second law of thermodynamics, which states that there must exist an entropy current S µ , whose divergence is positive semi-definite in an arbitrary φ I -offshell configuration. To wit where β µ is an arbitrary multiplier that can be chosen to be u µ /T using the inherent redefinition freedom in the hydrodynamic fields. A more helpful version of the second law is obtained by defining a free energy current which converts eq. (3.2) into the adiabaticity equation where we have denoted the Lie derivatives of g µν and φ I along β µ as To obtain the hydrodynamic constitutive relations allowed by the second law of thermodynamics, we need to find the most generic expressions for T µν and K I , within a derivative expansion, which satisfy eq. (3.4) for some N µ elastic and ∆. It is thus required to establish a derivative counting scheme. Following usual hydrodynamic treatments, we consider u µ , T , and g µν to be O(1) in the derivative expansion. On the other hand, we treat the derivatives of the scalars φ I as O(1), formally pushing the scalars themselves to O(∂ −1 ). We also treat the sources K ext I coupled to the scalars to be O(∂). This counting scheme is reminiscent of the one employed in the context of superfluids, and guarantees that the crystal cores composing the lattice, which are responsible for the elastic behaviour, appear at ideal order in the constitutive relations. Thus, we will be describing viscoelastic fluids with arbitrary strains, avoiding working in the restrictive regime of small strains as in [7].
Similar to the case of magnetohydrodynamics with dynamical gauge fields, not all terms in the adiabaticity equation (3.4) appear at the same derivative order. In particular, . This leads to order mixing in the constitutive relations, that is, the same transport coefficients can appear across derivative orders, forcing the analysis of the constitutive relations to consider multiple derivative orders simultaneously -an expression of one of the fallbacks of hydrodynamic formulations with dynamical fields. In sec. 4, we show that this problem can be avoided by working instead with formulations in terms of higher-form symmetries.

Ideal viscoelastic fluids
Given the establishment of a derivative counting scheme, we can use the adiabaticity equation (3.4) in order to find the constitutive relations of a viscoelastic fluid at ideal order. It is possible to infer that at leading order in derivatives, the adiabaticity equation has the solution The coefficient matrix σ IJ can be arbitrary except that its eigenvalues are constrained to be positive semi-definite. 8 Noting that This is the equivalent of the Josephson equation for superfluids and implies that the crystal fields are stationary at ideal order in derivatives. 9 In practice, this equation algebraically determines the time-derivatives of the crystal fields. It is useful to define the independent spatial derivatives of the crystal fields as where P µν = g µν + u µ u ν is the projector orthogonal to the fluid velocity. The spatial derivatives (3.8) capture all the onshell independent information contained in φ I .
In order to proceed further, we consider eq. (3.4) at one-derivative order, i.e. N µ elastic and T µν appear at ideal order in derivatives order while K I only appears at one-derivative order. The most generic constitutive relations are characterised by a free-energy current of the form where the fluid pressure P (T, h IJ ; h IJ ) is an arbitrary function of all the zero-derivative scalar fields in the theory, namely, the temperature T and the crystal metric h IJ = g µν e I µ e J ν . In particular we have allowed for an independent dependence on each component of h IJ . Additionally, h IJ labels the reference state of the material but has no inherent dynamics so hereafter we omit it for simplicity. Introducing (3.9) in the adiabaticity equation (3.4) and noting that ∇ µ (P β µ ) = δ B P + 1 2 P g µν δ B g µν along with we find the ideal viscoelastic fluid constitutive relations with ∆ remaining the same as eq. (3.6). In writing (3.11), we have defined the thermodynamic relations dP = sdT + 1 2 Thus, we can identify P as the thermodynamic pressure, as the energy density, s as the entropy density, and r IJ as the thermodynamic stress that models elastic responses. The φ I equation of motion now becomes The symbol σ has been used to draw a parallel with the respective term in magnetohydrodynamics, where higher-form fluids find another useful application [26,27]. There, the non-hydrodynamic field is the electromagnetic photon Aµ with the respective equation of motion given schematically as J µ = . . . − T σP µν δ B Aν + . . . = −J µ ext . 9 Note that, unlike superfluids, we do not have a chemical potential whose redefinition freedom could be used to absorb the plausible derivative corrections in eq. (3.7). Technically, the fluid velocity itself serves as a chemical potential along spontaneously broken translations. To see this, one can expand the Goldstones along a reference position as φ I = δ I i x i + δφ I and note that eq. (3.7) becomes u 0 ∂0δφ One can in principle absorb the derivative corrections into the redefinitions of u i , but such redefinitions will be incompatible with the manifest Lorentz covariance of the theory.
The constitutive relations (3.11) are quite general at this point but we will specialise to the case of an isotropic viscoelastic fluid later in section 3.4 leading to more familiar expressions. It is worth noticing that the same transport coefficient P (T, h IJ ) that is introduced at zero-derivative order in N µ elastic , appears at zero-derivative order in T µν but at one-derivative order in K I (via thermodynamic relations). However, as we will see in the next subsection, both T µν and K I get further corrections at one-derivative order. Hence, the constitutive relations for K I mix different derivative orders. This is the manifestation of order-mixing that we alluded to above.
For later use, it is helpful to explicitly write the energy-momentum conservation equations eq. (3.1a), given the constitutive relations (3.11). In particular, we find 14) or equivalently Formally, these equations can be used to eliminate u µ δ B g µν at one-derivative order in favour of P Iµ P Jν δ B g µν and δ B φ I .

One derivative corrections
The philosophy implemented for ideal viscoelastic fluids can also be extended to include one-derivative corrections to the constitutive relations. For simplicity, we focus on the elastic phase of crystals (as opposed to liquid crystals) for which k = d. The derivative corrections can naturally be classified into hydrostatic and non-hydrostatic constitutive relations: those that do not vanish when promoting β µ = u µ /T to an isometry and those that do vanish, respectively (see [47]).
In order to characterise the hydrostatic sector, we need all the one-derivative hydrostatic scalars that will make up the respective hydrostatic free-energy current. For this purpose, we list all the hydrostatic one-derivative structures The presence of the vectors u µ and e Iµ in the theory completely breaks the Poincaré invariance, so we can convert all of these into independent scalars 10 When k = d, this is no longer true and the counting of independent scalars needs to be more carefully implemented. Supplementing with arbitrary transport coefficients f 1 I , f 2 [IJ] , f 3 I(JK) as functions of T and h IJ , we construct the hydrostatic free energy density at first order in derivatives as with N µ elastic,hs = N β µ . Noting that ∇ µ (N β µ ) = δ B N + 1 2 N g µν δ B g µν and using the adiabaticity equation (3.4), we can read off the respective modified hydrostatic constitutive relations (see app. B for details). The free energy density is defined up to total derivative terms. Hence, it is possible to use the total derivative term ∇ µ e Iµ to eliminate the trace part of f 3 I(JK) ∼ f 3T I h JK and take f 3 to be traceless in the JK indices without loss of generality.
In the non-hydrostatic sector, the constitutive relations are the most generic expressions that involve δ B g µν and δ B φ I . At one-derivative order, the contribution to the respective free energy density happens to be zero, while the actual constitutive relations are We have defined T µν nhs = P Iµ P Jν T nhs IJ and have used the first order conservation equations (3.15) to eliminate u µ δ B g µν as well as to set the Landau frame condition T µν nhs u ν = 0. The associated quadratic form is given as The second law (3.4) requires that all the eigenvalues of the coefficient matrix are non-negative.
To summarise, the constitutive relations of a viscoelastic fluid, including the most generic one derivative corrections, are given by where the contributions T µν f i are given in app. (B.3). The φ I equations of motion modify to which is now correct up to two derivative terms. These constitutive relations describe the dynamics of a viscoelastic fluid fully non-linearly in strain. In the next subsection we focus on the linear regime.

Linear isotropic materials
For concreteness, we study the constitutive relations of an isotropic viscoelastic fluid. In this case, all the I, J, . . . indices appear due to the crystal metric h IJ and the reference metric h IJ . For simplicity, we work linearly in strain u IJ = 1 2 (h IJ − h IJ ), though the formalism introduced previously is sufficient to handle any possible non-linearities.

Constitutive relations
Firstly, we note that we cannot construct an odd-rank tensor or an antisymmetric 2-tensor (in field space) using just h IJ and h IJ . Therefore we are forced to set and hence the hydrostatic sector (3.18) at one-derivative order is rendered trivial. The ideal order pressure P can be expanded up to quadratic terms in strain as where This should be contrasted with the zero-temperature Lagrangian density in eq. (2.14). We have expanded the pressure up to quadratic terms because their derivatives can generically contribute to the constitutive relations with terms linear in strain via thermodynamics (see (3.12)). Thus In the non-hydrostatic sector, we can expand the coefficients η IJKL and σ IJ linearly in strain and obtain together with the associated quadratic form The ellipsis denote terms quadratic or higher order in strain. For positive semi-definiteness, the leading order transport coefficients η(T ), ζ(T ), and σ(T ) must be all non-negative, while the remaining ones are unconstrained. It should be noted that the transport coefficientζ u does not cause any dissipation, and is an example of non-dissipative non-hydrostatic transport in hydrodynamics. 11 In the end, the complete set of constitutive relations for an isotropic viscoelastic fluid up to first order in derivatives and linear in strain is given by where we have defined the expansion and shear of the fluid according to and used eq. (2.7). Using eq. (3.22), the φ I equation of motion takes the form The first line in eq. (3.29) contains the usual constitutive relations of an isotropic fluid, with η(T ) being the shear viscosity and ζ(T ) being the bulk viscosity. The terms in the second line correspond to lattice pressure P(T ), shear modulus G(T ), and bulk modulus B(T ), decoupled from the fluid except for the temperature dependence of the coefficients, which are present at zero temperature as well (see section 2.2). When P = 0, then the second line describes the well-known stresses of Hookean materials. The terms in the third and fourth lines denote one-derivative corrections that are linear in strain and correspond to the true coupling between fluid and elastic degrees of freedom. Such terms have not been explicitly considered in traditional treatments [11,13,14] neither in recent ones [4][5][6][7] and represent types of sliding frictional elements in rheology analyses. 12

Rheology and phenomenological models
Rheology is the study of stress/strain relations in flowing viscoelastic matter and is traditionally based on phenomenological models composed of mechanical building blocks designed for the purpose of describing observed properties of matter. The dynamics of viscoelastic materials studied in this paper is governed by energy-momentum conservation (3.1a) and the Goldstone equations (3.1c). In order to recast the equations in a more suitable form for comparison with rheology studies, it is useful to consider the implications of the Josephson condition (3.31), namely where £ β denotes the Lie derivative along β µ . These are the rheology equations. The first equation in (3.32) expresses the relation between the time-evolution of strains and viscous stresses while the second is a consequence of one of the basic assumptions in this work, namely, that the reference crystal metric is non-dynamical (i.e. absence of plastic deformations). This corresponds to the elastic limit in the language of [7] (see also app. C.1). Given the rheology equations (3.32), one can compare the constitutive relations found here with existent viscoelastic models. First of all, it should be noted that the last two lines in (3.29) describe several couplings between fluid and elastic degrees of freedom and a proper account of them in material models has not been considered in generality. Doing so requires introducing many new mechanical building blocks of the sliding frictional type. For simplicity, we consider the case in which P = η u 1 = η u 2 = ζ u 1 = ζ u 2 =ζ u = 0 which leads to the energy-momentum tensor This form of the stress tensor, together with (3.32), is known as the Kelvin-Voigt model and usually represented as in figure 1(a), where we have focused on the bulk stresses sector (i.e. we have depicted the effect of bulk viscosity and bulk elastic modulus) and ignored the ideal fluid part.
Another model that illustrates the use of coupling terms between elastic and fluid degrees of freedom is the Bingham-Kelvin model for which P = η u 1 = η u 2 = ζ u 2 =ζ u = 0 and the energymomentum tensor becomes where we have ignored the shear contribution to the stresses. The last term in (3.34) is the term responsible for the frictional slide element (black box) depicted in figure 1(b).
The entire possibility of linear responses (3.29) allows for more intricate and rich material diagrams. The full nonlinear theory of section 3.3 also allows for nonlinear responses to strain and hence for the description of non-Newtonian fluids. However, it is not capable of describing Maxwell-type models or Zener models as these violate the second condition in (3.32). Such models allow for plastic deformations and require that we consider a dynamical reference crystal metric h IJ as in [7]. We intend to pursue this generalisation in the future.

Linearised fluctuations
In this section we study linearised fluctuations of equilibrium states of isotropic crystals. We consider crystals coupled to a flat background g µν = η µν and vanishing external sources K I ext = 0, with static equilibrium configurations given by Note that in equilibrium we have h IJ = δ IJ , corresponding to a crystal subjected to no strain. In general the system also admits solutions of the type φ I = α x I corresponding to a uniform strain. Such configurations are allowed for a space filling crystal, but will need to be supplied with appropriate boundary conditions if the crystal was finite in extent. Since such configurations can be obtained by a trivial rescaling of φ I 's, we do not consider them here.

Modes
Let us consider small perturbation of the equilibrium state parametrised by Plugging in a plane wave ansatz and solving the equations of motion (3.15) linearly in the perturbations, we can find the solutions We have suppressed the arguments of various transport coefficients but they are understood to be evaluated on the equilibrium configuration. In addition, we have omitted the effect of viscosities but we consider it explicitly below. Primes denote a derivative with respect to temperature. We have also used the isotropy of the system to decompose using the same transport coefficients introduced in sec. 3.4. 13 The symbols A and A I ⊥ (with A I ⊥ k I = 0) denote arbitrary amplitudes corresponding to "longitudinal" and "transverse" modes | 22 respectively. The respective dispersion relations, in small momentum and frequency regime, are given by Solving these equations, we find that in the longitudinal sector we have the usual sound mode along with a new diffusion mode characteristic of a lattice 14 where On the other hand, in the transverse sector we have another sound mode where We see that the transverse sound mode is controlled by the shear modulus G. On the other hand, the new diffusive mode is controlled by the transport coefficient σ and η. In the absence of the lattice pressure P, these expressions simplify which might be more familiar to some readers.
For linear stability of the system, the imaginary part of ω(k) must be non-negative. This leads to the constraints v 2 , v 2 ⊥ , Γ , Γ ⊥ , D > 0. In terms of coefficients, assuming that T s + P > 0 and the second law constraints η, ζ, σ ≥ 0, we land on the parameter space On the other hand, for causality, we require v 2 , v 2 ⊥ < 1. This leads to (s + P ) 2 This gives the allowed range of parameters for a sensible evolution of the dynamical equations.

Linear response functions and Kubo formulas
We can extend the analysis above to read out the linear response functions of the theory by switching on plane wave background fluctuations. Let us start by setting σ → ∞ and η, ζ → 0 turning off the dissipative corrections for simplicity. Let us take a perturbation of the background sources We can read out the solution of the equations of motion (3.15) by a straightforward computation Without loss of generality, we can choose the momentum to be in k I = δ I x and denote the remaining spatial indices by a, b, . . .. Defining we read out the respective two point functions; in the longitudinal sector we have Similarly, in the transverse sector Finally, we have non-zero contributions in the cross sector Note that all the correlators with odd number of transverse indices vanish due to isotropy G tt,ta T T,R = G tx,ta T T,R = G tt,xa T T,R = G tx,xa T T,R = G xx,xa T T,R = G xx,ta T T,R = G ta,bc T T,R = G xa,bc T T,R = 0 , G x,ta φT,R = G x,xa φT,R = G a,tt φT,R = G a,tx φT,R = G a,xx φT,R = G a,bc φT,R = 0 , G x,a φφ,R = 0 . Upon turning on the dissipative transport coefficients, the two-point functions become much more involved. However, we report the respective Kubo formulas These can be used to read out the transport coefficients in terms of linear responses functions. 15 This finishes our quite detailed discussion of viscoelastic fluids. We have written down the most generic constitutive relations determining the dynamics of a viscoelastic fluid up to first order in the derivative expansion. In particular, we specialised to linear isotropic materials and obtained the respective constitutive relations, modes, and linear response functions. In the next section we present an equivalent formulation of viscoelasticity in terms of a fluid with partially broken higher-form symmetries.

| Viscoelastic fluids as higher-form superfluids
In this section we present a formulation of hydrodynamics with partially broken generalised global symmetries and show their relation to the theory of viscoelastic fluids formulated in the previous section. Generalised global symmetries are an extension of ordinary global symmetries with one-form (vector) conserved currents and point-like conserved charges to higher-form conserved currents and higher dimensional conserved charges such as strings and branes [52]. It has been observed that when a fluid with a one-form symmetry 16 has its symmetry partially broken along the direction of the fluid flow, it implements a symmetry-based reformulation of magnetohydrodynamics [24,26,27,53]. In this section, we extend this partial symmetry breaking to hydrodynamics with multiple higher-form symmetries and show that the resultant theory is a dual description of viscoelastic fluids with translation broken symmetries in arbitrary dimensions. In d spatial dimensions, one requires d number of partially broken (d − 1)-form symmetries in order to describe viscoelasticity. The case of d = 2 involving two one-form symmetries was considered in [22], albeit in a very restrictive case and ignoring the issues that require partial symmetry breaking. The understanding of partial symmetry breaking is essential for consistency of higher-form hydrodynamics with thermal equilibrium partition functions, as has been previously observed in [25,26].

A dual formulation
In viscoelastic fluids with translation broken symmetries, all the dependence on the crystal field φ I comes via its derivatives e I µ . As we have argued in section 3, the φ I equations of motion can be used to eliminate u µ e I µ in favour of P Iµ = P µν e I ν and other constituent fields in the theory. Thus, it should be possible to reformulate the physics of viscoelastic fluids purely in terms of P Iµ and the hydrodynamic fields u µ and T , without referring to the microscopic fields φ I .
To make this precise, we formally define a set of d-form currents associated with the viscoelastic fluid by Hodge-dualising the derivatives of φ I as It is understood here that the φ I equations of motion have been taken onshell to eliminate δ B φ I in terms of the remaining fields and background sources. Due to the symmetry of partial derivatives, these currents are conserved by construction A priori, these conservation equations have d(d + 1)/2 independent components for every value of I but only d of these contain a time-derivative and hence govern dynamical evolution, while the remaining d(d − 1)/2 components are constraints on an initial Cauchy slice. Conspicuously, these are the exact number of dynamical equations required to evolve the d physical components in P Iµ for every value of I (note that u µ P Iµ = 0).
The conservation equation (4.2) implies that there is a set of d topological conserved charges Q I of the form where Σ 1 is a given one-dimensional surface and is the Hodge operator in d+1 dimensional spacetime. The charges Q I [Σ 1 ] count the number of lattice hyperplanes that intersect the one-dimensional surface Σ 1 .
The current (4.1) couples to the field strength of a higher-form gauge field. More precisely, we can replace the external currents K ext I by the field strength such that Since H Iµ 1 ...µ d+1 is a full-rank form, locally it can be re-expressed as an exact form Using this definition, the energy-momentum conservation equation (3.1a) takes the form The dynamics of the hydrodynamic fields u µ and T , and P Iµ is governed by energy-momentum conservation (4.7) and d-form conservation equations (4.2). 17

Ordinary higher-form hydrodynamics
Having motivated a dual formulation of viscelastic fluids in terms of higher-form symmetries, we consider higher-form hydrodynamics in its own right, following [25,26]. Consider a fluid living in (d + 1)-dimensions that carries a conserved energy-momentum tensor T µν and a k number of conserved d-form currents J Iµ 1 ...µ d where I = 1, 2, . . . , k. When coupled to a background metric g µν and background d-form gauge fields b Iµ 1 ...µ d , the associated conservation equations are given as 17 Note that, by the definition of the d-form currents, we have the following relation which can be understood as a frame choice from the higher-form hydrodynamic perspective. In general, we can choose a different set of fields in the hydrodynamic description which are aligned with P µ I in this particular frame, but can be arbitrarily redefined otherwise.

| 27
In a generic number of dimensions, the conservation equations lead to (d+1+kd) dynamical equations and kd(d − 1)/2 constraints. From this counting procedure, it can be checked that eqs. (4.9) can provide dynamics for a set of symmetry parameters (4.10) Under an infinitesimal symmetry transformation parametrised by X = (χ µ , Λ χ Iµ 1 ...µ d−1 ), they transform according to where £ χ denotes the Lie derivative with respect to χ µ . Let us repackage these fields into the fluid velocity u µ , temperature T , and (d − 1)-form chemical potentials µ Iµ 1 ...µ d−1 according to such that u µ u µ = −1. Interestingly, while the fields u µ and T are gauge invariant, µ Iµ 1 ...µ d−1 transform akin to (d − 1)-form gauge fields Hence, they have the required kd physical degrees of freedom, which, along with u µ and T , match the number of dynamical components of the conservation equations.
Similar to our discussion in section 3.1, higher-form fluids need to obey a version of the second law of thermodynamics. In the current context, this statement translates into the existence of an entropy current S µ that satisfies Compared to eq. (3.2), here we have also taken into account the higher-form conservation equation and the respective multiplier has been chosen to be µ Iµ 1 ...µ d−1 /T using the inherent redefinition freedom in the higher-form chemical potential. It is straightforward to formulate a higher-form analogue of the adiabaticity equation to ease the derivative of the constitutive relations. Defining the free energy current the adiabaticity equation for higher-form fluids reads We have identified the variations of the various background fields according to To obtain the constitutive relations of a higher-form fluid, it is required to find the most generic expressions for T µν and J Iµ 1 ...µ d in terms of the dynamical fields u µ , T , and µ Iµ 1 ...µ d−1 , as well as background fields g µν and b Iµ 1 ...µ d , arranged in a derivative expansion, that satisfy eq. (4.16) for some N µ and ∆.

Partial symmetry breaking
When a higher-form symmetry is partially broken in its ground state, the hydrodynamic description should include the associated Goldstone modes ϕ Iµ 1 ...µ d−2 with u µ 1 ϕ Iµ 1 ...µ d−2 = 0, that transform according to 18 This allows us to define a gauge-invariant version of the (d − 1)-form chemical potentials The Goldstone fields are reminiscent of the crystal fields φ I from sec. 3.1 and are accompanied by their own equations of motion, which can be schematically represented as (4.20) As in the previous formulation, the operator K Iµ 1 ...µ d−2 is not predetermined without the knowledge of the microscopics but we can fix its form up to certain transport coefficients by imposing the second law of thermodynamics. In this context, this translates into the requirement that the fluid must admit an entropy current S µ whose divergence must be positive-semi-definite in any arbitrary ϕ I -offshell configuration. Fixing the Lagrange multipliers associated with the respective conservation equations to be β µ and ζ Iµ 1 ...µ d−1 , the statement of the second can be written as which is different than the corresponding statement in the symmetry unbroken phase given in (4.14).
While ζ Iµ 1 ...µ d−1 is fundamentally more transparent, for most of the explicit computations, it will be helpful to work with a Hodge-dualised version where we note that u µ ψ µ I = 0. This representation makes it clear that ψ µ I , and hence ζ Iµ 1 ...µ d−1 , have the same degrees of freedom as P Iµ and can be used as a fundamental hydrodynamic field for viscoelastic fluids instead. The relation between the two is typically non-trivial and needs to be obtained order-by-order in the derivative expansion.
Compared to eq. (3.4), the adiabaticity equation (4.16) in the dual formulation does not exhibit order mixing, as both δ B g µν and δ B b Iµ 1 ...µ d are O(∂). This allows for a more transparent analysis of the constitutive relations, as we shall see in the next section. Another benefit of working in the dual formulation is that we directly obtain the constitutive relations for (the Hodge dual of) the physically observable crystal momenta, rather than for the equations of motion of the crystal fields. This can considerably simplify the computation of the respective correlation functions and Kubo formulae, as in the case of magnetohydrodynamics [25].

Revisiting ideal viscoelastic fluids
The constitutive relations of an ideal viscoelastic fluid in higher-form language are characterised by an ideal order free energy current (4.27) Here p(T, γ IJ ) is an arbitrary function of all the available ideal order scalars in the theory, namely the temperature T and matrix Introducing this into eq. (4.16) and noting that 19 we obtain the ideal viscoelastic fluid constitutive relations 20 Here we have defined These are the most generic constitutive relations of an ideal viscoelastic fluid in higher-form language. As promised earlier, there is no order-mixing across derivative orders in this formulation. It is useful to explicitly work out the equations of motion for ideal viscoelastic fluids, which result in In order to find the relation between the two formalisms, we need to perform the identification according to eq. (4.1), which results in the following map between formulations Additionally, introducing this into the energy-momentum tensor, we find while temperature, entropy density and energy density agree in both formulations.

One derivative corrections
Following the same arguments as the previous subsection, we consider the hydrostatic free energy density in the dual picture The hydrostatic constitutive relations can be obtained using the variations given above and refer the reader to appendix B for details. In turn, in the non-hydrostatic sector, we can expand The most general non-hydrostatic constitutive relations are correspondingly The transport coefficient matrices have necessary symmetry properties. The positivity constraint requires that the symmetric part of the 5 × 5 transport coefficient matrix is positive semi-definite.
In order to provide the map of transport coefficients at first order between the two formulations, we first use the identification in eq. (4.1) in order to obtain This provides a definition of P µ I and δ B φ I in the conventional formulation in terms of the dual formulation variables We wish to begin the comparison with the free energy currents in the two formulations. Using eq. (4.15), we know that Using the map (4.39) and the results from appendix B, we infer the map between the hydrostatic transport coefficients To obtain the map in the non-hydrostatic sector, we need to compare the energy-momentum tensors and δ B φ I in the two formulations. In hindsight, we allow for a relative field redefinition of the fluid velocity between the two formulations, i.e u µ elastic = u µ + δu µ with u µ δu µ = 0. The ϕ I -equation of motion (3.22), upon using the said field redefinition, implies that On the other hand, using the results from appendix B, it is straight-forward, albeit cumbersome, to obtain that T µν elastic,hs = T µν hs + 2T u (µ sδu ν) − ψ ν) Given that all the non-hydrostatic corrections are in the Landau frame, we must choose the velocity field-redefinition to be δu µ = ψ µ I δ B φ I /s, mapping the hydrostatic sectors of the two formulations to each other. Comparing δ B φ J from eq. (4.42) to eq. (4.39) and comparing the non-hydrostatic energy-momentum tensors in the two formulations, we obtain the map This completes the formulation of viscoelastic hydrodynamics in terms of generalised global symmetries and shows that it can exactly accommodate viscoelastic hydrodynamics with broken translation invariance. In the next section we look at particular realisations of both these formulations in the context of holography.

| Conformal viscoelastic fluids and holography
In this section we provide, and study the properties of, holographic models in D = 4, 5 bulk dimensions (i.e. d = 2, 3 spatial dimensional fluids). The models we consider in general break conformal symmetry due to double trace deformations but conformal symmetry can be recovered in a specific case. Thus, in the beginning of this section we consider conformal fluids. In connection with viscoelastic holography, we consider two classes of models that have been considered in the literature. The first class of models has translation broken symmetries involving a set of (D − 2) scalar fields Φ I minimally coupled to gravity. The second class is formulated in the context of generalised global symmetries and involve a set of (D − 2) gauge fields B Ia 1 ...a D−2 minimally coupled to gravity [22]. The latter class describes particular equilibrium states of the higher-form hydrodynamics described in sec. 4, which we explicitly show by generalising the work of [22] to D = 5. Given that in this case the dual fluid is governed by conservation equations alone (i.e. no dynamical fields), which has been one of the motivations in the holographic digressions of [22,54,55], we consider it first. Later, generalising aspects of [56], we "dualise" the model with higher-form symmetries and obtain the class of viscoelastic models with translation broken symmetries, which consist of the model of [39] but with an alternative quantisation of the scalar fields and a double trace deformation of the boundary theory. We show that this process results in the dual fluid given in section 3.

Conformal viscoelastic fluids
A viscoelastic fluid is said to be conformal if it is invariant under the conformal rescaling of the background metric g µν → Ω 2 g µν for some arbitrary function Ω(x). In practice, it implies that the energy-momentum tensor of the theory is traceless (modulo conformal anomalies) and the constitutive relations are only constructed out of the conformal covariants.

Constitutive relations
Focusing on the non-anomalous case, setting the trace of the energy-momentum tensor (3.21) to zero, we get certain constraints on the respective transport coefficients, namely Furthermore, requiring T µν and K I to only involve conformal covariants requires The first equation in eq. (5.1) determines the energy density in conformal fluids as expected. At the one-derivative hydrostatic order, we see that we are only left with f 2 [IJ] which is the only conformal covariant term in the free-energy current. In the non-hydrostatic sector, we essentially just eliminate the conformal-non-invariant ∇ µ u µ term from the constitutive relations. Consequently we get and We will now focus on a special case of these constitutive relations.

Linear conformal isotropic materials
For conformal viscoelastic fluids truncated to linear order in strain, we need to additionally impose the constraints (5.1) on top of the constitutive relations eq. (3.29), leading to This gives the following constitutive relations of a linear conformal viscoelastic fluid while the φ I equation of motion is still given by eq. (3.31). In the special case that the internal pressure of the lattice P = 0, we infer that the bulk modulus B = 0 and these constitutive relations, along with the φ I equations of motion (3.31), simplify to The first line of the energy-momentum tensor represents the constitutive relations for an ordinary uncharged conformal fluid. The second line has the expected shear modulus term along with the variants of shear viscosities η u 1 , η u 2 representing the coupling of the conformal fluid to the strain of the crystal.
We would like to note that the constitutive relations (5.5) are, in principle, different from the ones obtained in [8]. The strain u µν defined in eq. (2.4) transforms inhomogeneously under a conformal transformation, i.e. u µν → Ω 2 u µν + 1 2 Ω 2 − 1 h µν , because conformal rescaling only acts on the physical distances and not on the reference distances between the crystal cores. It follows that the Using the renormalisation procedure of (5.18) we evaluate the on-shell Euclidean action in order to find the pressure and extract the temperature and entropy from the black brane geometry whereX given by (5.29) but where the contraction is performed with γ µν . Variation of the total action with respect to Φ I , upon using the bulk equations (5.31) leads to the boundary term for some scalar operator O I . It is clear from here that if the boundary action (5.32) is considered then the set of Φ I are taken as sources in the boundary theory. In the hydrodynamic limit, this is the setting of forced fluid dynamics [45] and of momentum relaxation [44] in which case the scalar fields Φ I are background sources to which the fluid couples to but not dynamical fields as in viscoelasticity. 25 In order to use the holographic model (5.29) for describing viscoelastic materials, another type of boundary conditions is necessary.

Dualising the holographic model
It is well known that in D = 3, the dynamics of a U (1) gauge field A µ is equivalent to the dynamics of a scalar field Φ since dA ∼ dΦ. It is also known that a massless scalar field in D = 3 AdS admits two possible quantisations corresponding to different dimensions of the boundary theory operator (∆ = 1 ± 1) [58]. This was exploited in [56] to show that the correct boundary conditions that describe the dynamics of the gauge field A µ are not those that fix the scalar Φ to be the source but instead those that fix its conjugate momentum. This corresponds to the quantisation with ∆ = 0.
In this section we generalise the analysis of [56] to higher-form fields in order to dualise the model of sec. 5.2. The difference between the cases considered here and that of [56] is that the dualisation takes a theory with mixed boundary conditions (in the sense of [59]) and maps to another theory with mixed boundary conditions. Naturally, in D bulk dimensions, the dynamics of the bulk gauge field B Iµ 1 ...µ D−2 is equivalent to the dynamics of a scalar field Φ I since dB I ∼ dΦ I . This can be seen directly at the level of the path integral. Consider the action for the B I field as in (5.12) but integrate over the field strength H I = dB I instead of over B. In such case, one needs to enforce the Bianchi identity d H I = 0 by introducing a Lagrange multiplier Φ I such that such that the path integral becomes 36) which is that of a massless scalar field in D dimensions. Having established the duality at the level of the bulk path integrals, one may include the boundary action. Focusing just on the boundary Acting with the covariant derivative on the boundary stress tensor, using the Codazzi-Mainardi equation n a R aµ = −∇ µ K + ∇ ν K ν µ (see e.g. [61]) and the bulk equations (5.31), one obtains the Ward identity Comparing this with (3.1a) we identify K ext I = Π I and O I = φ I . Thus fixing the boundary value of the source Π I provides dynamics for the Goldstone scalars φ I and has the interpretation of applying external forces to the crystal lattice. We will now study thermal states within the model (5.29) with boundary action (5.40).

Thermal state in D = 4
The bulk metric in D = 4 was considered in [37] but the thermodynamic properties have not been properly evaluated. The metric takes the form (5.14) but with metric function and scalar fields where r = r 0 is the location of the horizon. We now wish to determine the thermodynamics of this black brane and the holographic stress tensor and scalar currents for later interpretation in terms of a viscoelastic fluid. Noting that the free energy of the viscoelastic fluid is (minus) the pressure P as in (3.9), we obtain the pressure from the onshell action while the entropy and temperature are extracted from the black brane where we have set M p = 1. In order to obtain the stress tensor we use (5.44) and for the scalar operators we use (5.41), finding while the sources in this case vanish, i.e. Π I = 0. 26 We now identify the thermodynamic properties of the viscoelastic fluid by comparison with (3.11). We find These quantities satisfy the thermodynamic relations (3.12). Introducing these quantities in the map (4.3), we obtain exactly the same thermodynamic properties as in section 5.2.3 provided we identify r h = r 0 and m = m. In the case M = 0, these thermodynamic quantities describe a conformal fluid as in sec. 5.1.2.
We would like to note that the strain of the viscoelastic fluid is given by The trace of the energy-momentum tensor is non-vanishing except if M = 0, in which case the theory is conformal.

| 42
Therefore m, in some sense, controls the strength of the elastic strain and hence holography can provide models of viscoelasticity with arbitrary strains. 27 Let us focus on the linear regime to make contact with section 5.1.2. Expanding and noting section 3.4.1, we can read out the lattice pressure and bulk modulus The fluid pressure P f can be read out trivially from P by setting strain to zero. We cannot comment on the shear modulus G because we are working in a state with diagonal strain. We find that the coefficients P and B do not scale homogeneously under conformal transformations. This is in accordance with the comments made on inhomogeneous conformal scaling of strain in section 5.1.2.
Since the thermodynamic properties derived in (5.47) describe an equilibrium state with vanishing elastic shear tensor, it is not possible to extract from it the shear modulus. Thus, we do not have a complete knowledge of the transverse phonon and the longitudinal sound dispersion relations. However, [22] showed that for small m the dispersion relations under the assumption of vanishing shear modulus agree with numerical results. 28 According to the analysis of [22], it is expected that a Gregory-Laflamme like instability [62] is present for specific values of the parameters, including when M = 0. However, a more in-depth analysis is necessary in order to make definite statements.

Thermal state in D = 5
Bulk metrics dual to viscoelastic fluids in D = 5 have not been studied in depth but they straightforwardly generalise their lower dimensional counterpart. The bulk metric function and scalar fields that solve (5.31) are given by The onshell action (i.e. pressure), the temperature and entropy are given by 54) 27 In previous considerations of viscoelastic holography (see e.g. [35]), the bulk field ΦI has been related to the crystal displacement field δφ I = φ I − x I at the boundary and not with φ I itself. Unlike our model, where the strainless limit is given by m = 1/ √ 2, this alternative choice places the strainless limit at m = 0. Realising that the theory in the bulk becomes an ordinary charged black brane at m = 0 that is known to be dual to a pure fluid at the boundary and not an unstrained crystal, we do not make this choice. Furthermore, it is unclear if this choice can be implemented at a non-linear level in strain. A similar choice has been made in the higher-form setup of [55], but the authors there approached it as fluctuations around a state without "dynamical defects" (no crystal cores), distinct from a crystal phase where such defects are obviously present. 28 We believe that the discrepancy between hydrodynamic and numerical results identified in [22] is due to the fact that [22] assumed a vanishing shear modulus in their hydrodynamic calculations.
while the boundary stress tensor takes the form Comparing with the constitutive relations (3.12), we identify Again, these quantities satisfy the thermodynamic relations (3.11) and using the map (4.3) they lead to the constitutive relations and thermodynamic relations of sec. 5.2.4 provided r h = r 0 and m = m.
Similar to D = 4, we can perform a small strain expansion. We find that and corresponds to a conformal fluid when M = 0. Once again, the transport coefficients appearing here are not homogeneous under conformal scalings.

| Outlook
In this paper we introduced two novel formulations of relativistic viscoelastic hydrodynamics capable of dealing with elastic and smectic crystals phases. The first formulation of sec. 3 follows traditional treatments [11,13,14] where the elastic degrees of freedom are described by the dynamics of Goldstones of translational broken symmetries. However, it generalises earlier literature by considering the effect of external currents, anisotropy, and nonlinearities. When applied to the case of linear isotropic materials we uncovered 6 new transport coefficients (5 dissipative and 1 non-dissipative) in section 3.4 that characterise the coupling between elastic and fluid degrees of freedom, constituting sliding frictional forces in viscoelastic material diagrams. The second formulation of section 4 uses the framework of generalised global symmetries in order to recast traditional viscoelastic treatments as higher-form superfluidity. This provides a fully symmetry-based approach to viscoelastic hydrodynamics and we show how the two formulations map one-to-one.
In section 5 we studied holographic models to both these formulations and proposed a new and simple model for viscoelasticity, consisting of the model of [39] with an alternative quantisation for the scalar fields and a double trace deformation. We also classified conformal linear isotropic viscoelastic materials in section 5.1 and shown that they correctly reproduce holographic results when there is no double trace deformation. We also identified new holographic transport coefficients, which have not been considered in earlier holographic works [22,[29][30][31][32][33][34][35][36]. Namely, by expanding the equation of state in a small strain expansion, we notice that there is a linear term in strain in the free energy, which is the lattice pressure and usually ignored in classical elasticity treatments (see sec. 5.3).
The work presented here naturally opens up the possibility for various extensions and generalisations which we now describe.
Non-homogeneity and dynamical reference metric: The hydrodynamic formulations considered here assumed homogeneity of the crystal lattice and a non-dynamical reference metric. As such, it was shown that phenomenological viscoelastic models such as the Kelvin-Voigt and Bingham-Voigt models are special cases of the general constitutive relations obtained here. However, other existent models such as the Maxwell model are not captured within this approach. In order to do so, it is required to consider dynamical reference metrics, allowing for the possibility of plastic deformations. This has been implemented in [7] but only for linear strains. 29 Disclinations and dislocations: We have assumed that the fields φ I are surface forming, that is, the Bianchi-type condition in (3.1a) is satisfied. This implies that no defects (e.g. disclinations or dislocations) are present. It would be interesting to understand viscoelastic hydrodynamics in the presence of defects. In terms of generalised global symmetries this implies that the higher-form currents J Iµ 1 ...µ d−1 are not conserved, which makes it harder to understand from this dual point of view. However, it may be the case that in some cases, the violation of current conservation can be understood as an anomaly in quantum field theories with generalised global symmetries.
Other crystal phases: This paper was mostly concerned with elastic (solid) crystal phases, although some of the results are valid for smectic-A phases. However, liquid crystals can exhibit many other types of phases such as other smectic, nematic, and hexatic phases. The hydrodynamics of these phases have been consider in traditional treatments [11,13,14,67], though without a careful analysis of the constitutive relations. It would be interesting to revisit these works using modern hydrodynamics and to develop equivalent models in terms of higher-form symmetries.
Charged lattices, charge density waves, and holography: Charged Wigner crystals and charge density waves are charged generalisations of elastic and nematic phases of neutral crystals. It would be interesting to consider such extensions as it can aid in the understanding of recent holographic studies [50,68,69]. It is natural to consider the works [29,32,33,37,38] and charged generalisations [50,68,69] within the framework of generalised global symmetries. We would also like to understand whether sec. 3 is sufficient for modelling the viscoelastic fluids encountered in [29][30][31][32][33][34][35][36]).
Finite relaxation time: As mentioned in the introduction to this work, Maxwell's original idea of viscoelasticity consisted of materials that exhibited elasticity at short time scales and fluidity at long time scales. In this paper we focused on situations in which elasticity and fluidity coexist at long time scales and long wavelengths by assuming the strain relaxation times to be very large. It would be interesting to consider the case of arbitrary finite relaxation times as in [7] in such a way that deviations away from the hydrodynamic regime are under control. In these situations, the framework of quasi-hydrodynamics will most likely be useful, as in [70].
Fluid/gravity of viscoelastic fluids: We explored holographic models to viscoelastic hydrodynamics in section 5 but only at ideal order in a long wavelength expansion. It is clear from the analysis of section 5 that the holographic models describe constitutive relations nonlinear in strain. As it is non-trivial to categorise all possible materials nonlinearly in strain, it would be interesting to continue the expansion one order higher and to uncover viscoelastic transport coefficients that are present in gravity, ultimately obtaining a possible phenomenological model of viscoelasticity. In this context, it will be possible to uncover the conformal fluid structure of section 5 for small strains.
Shear elastic modulus and instabilities: In section 3.5, we have performed a thorough analysis of the dispersion relations of viscoelastic fluids, identified a longitudinal sound, transverse sound and a diffusive mode as well as generalised the relation between longitudinal and sound modes in a conformal solid [31] to the case of finite temperature and in the presence of lattice pressure. In order to apply these results to the holographic models of section 5, it is necessary to obtain the shear elastic modulus. This transport coefficient does not follow from the equation of state in equilibrium since the equilibrium state that we have considered has vanishing elastic shear tensor. The Kubo formulae of section 3.5 can be used in the holographic models studied here to obtain the shear modulus. Nevertheless, the results of [22] suggest the existence of an instability for certain values of parameters in the model of [39] with alternative boundary conditions. It would be interesting to obtain the shear modulus precisely and study instabilities in these models more thoroughly.
Finally, the work presented here shows that formulating hydrodynamics in terms of generalised global symmetries can be extremely useful, not only because it allows to rewrite hydrodynamic theories with dynamical fields just based on symmetries (and their spontaneous breaking), but also because it allows for a cleaner understanding of potential holographic models and their boundary conditions. It should be noted that, as in the case of magnetohydrodynamics, viscoelasticity when written in the language of generalised global symmetries has global symmetries partially spontaneously broken (along the fluid flows). In the context of condensed matter systems, broken global symmetries are only natural [71]. At this point, we are not aware of a physical hydrodynamic system with unbroken generalised global symmetries at the boundary as described in [27]. This fact has repercussions to several other constructions studied in [70], where the Goldstone modes of spontaneous broken global symmetries have not been taken into account. We wish to study these constructions more carefully in the near future. connection the acts on the crystal indices as 30 C I µJ = −e λ J ∇ µ e I λ = −e λ J ∂ µ e I λ + e λ J Γ σ µλ e I σ , (A.1) where ∇ µ denotes the spacetime covariant derivative compatible with g µν and associated with the Christoffel connection Γ λ µν . Additionally, we introduce the covariant derivative D µ , associated with Γ λ µν and C I µJ and compatible with both g µν and h IJ . The covariant derivative D µ transforms as a tensor under GL(k) transformations of the normal one-forms. With this definition at hand, it is easy to check that h λν D µ e I ν = h λν D µ e ν I = D µ h IJ = D µ h IJ = 0 . (A. 2) The projection of the structures D µ e I λ , D µ e λ I along the crystal directions vanishes as in (A.2) but in general D µ e I ν =h ν λ ∇ µ e I λ . (A.3) Finally, the curvature associated with the connection C I µJ is not independent and is related to projections of the spacetime Riemann curvature tensor R µν ρ σ , that is which is a generalised Ricci-Voss equation (see e.g. [74]).
The crystal metric h IJ and the reference metric h IJ contain all the information about the internal geometry of the crystal and, in particular, about the deformations of the crystalline structure. The one-forms e I µ , however, contain additional information about the shape and orientation of the crystal as embedded into the spacetime. If the spacetime does not have broken rotational invariance except for the existence of the crystal itself, this extra information can only be accessed via derivatives. 31 Given the structures introduced above, we see that all the information about the derivatives of e I µ is stored in D µ e I ν and C I µJ . In turn, the derivatives of h IJ are all captured by Hence, D µ e I ν and C [IJ] µ can be seen as containing the additional information about the crystal embedding. C [IJ] µ captures the differential of the local SO(k) orientation of the crystal with respect to the reference coordinate system, while D µ e I ν captures the shape of the crystal cores via the tangential extrinsic curvature tensor K I µν =h µ λ D λ e I ν and that of the crystal via the normal extrinsic curvature tensor L IJ ν = −e Jλ ∇ λ e I ν . It is well known that a generic set of one-forms e I µ does not have to be surface forming, i.e. there might not exist a foliation of crystal core worldsheets normal to all the e I µ . For this to be the case, one needs to invoke the Frobenius theorem, ensuring that there must exist a set of spacetime one-forms a I νJ such that ∇ [µ e I ν] = −a I [µJ e J ν] . This is equivalent to the condition that the extrinsic curvature tensor is symmetric K I µν = K I νµ . Introducing the normal one-forms as in sec. C | Comparison with previous works

C.1 Comparison with Fukuma-Sakatani formulation
In this section we compare our work with that of [7]. The work of [7] differs from ours in several ways: it assumes small (linear) strains; isotropy of the material; absence of external forces; introduces an extra scalar field related to breaking of time translations; it considers the possibility of a dynamical reference metric; and introduces a finite strain relaxation time. By contrast, the work presented here works fully non-linear in strain; incorporates external forces; includes anisotropy; assumes a non-dynamical reference metric; does not break time translations; and considers very large relaxation times. Within this regime, it is capable of capturing many viscoelastic effects which [7] cannot.
As a result of not breaking time translations, the strain tensor in the work presented here is purely spatial. The (d + 1)-dimensional spacetime strain tensor u µν is defined via the pullback of the d dimensional crystal strain tensor u µν = e I µ e J ν u IJ , (C.1) which implies that there exists a timelike vector v µ such that v µ u µν = 0. Physically, the vector v µ can be understood as a reference frame in which the stress tensor is purely spatial. In equilibrium, this frame is the same as the fluid frame of reference. We see that the formulation of viscoelastic fluids given in this paper necessitates the strain tensors to admit such a reference frame. 32 In the notation of [7] we identify E FS µν = u µν where the upperscript FS denotes quantities defined in [7]. We can split the strain tensor into components according to ε µν FS = P µρ P νσ E FS ρσ = P Iµ P Jν u IJ , tr ε FS = P µν E FS µν = h IJ u IJ + T 2 u IJ δ B φ I δ B φ J , ε µ FS = −2P µρ u ν E FS ρν = −2T P Iµ u IJ δ B φ J , θ FS = u µ u ν E FS µν = T 2 u IJ δ B φ I δ B φ J . (C. 2) The extrinsic curvature defined in [7] is given as We first examine the φ I equations of motion. In sec. 3 we deduced that and hence, in the formulation presented in this paper, ε µ FS is pushed to O(∂) and θ FS to O(∂ 2 ) and are both algebraically determined. On the other hand, noting that we also have one spatial direction. We can define (d − 1)-and d-dimensional volume forms We can also define a projector In terms of these, the ideal order constitutive relations (4.30) reduce to T µν = ( + p) u µ u ν + pg µν − q 11 γ 11 Π µν + O(∂) , This can be compared directly to section 2 of [25] with Q = q 11 √ γ 11 and µ = √ γ 11 . In principle, we can extend this comparison to include one-derivative corrections as well. However for technical simplicity, in this work we have only focused on one-derivative corrections for elastic fluids (k = d), so such a comparison is beyond the scope of the current analysis.