Holographic Viscoelastic Hydrodynamics

Relativistic fluid hydrodynamics, organized as an effective field theory in the velocity gradients, has zero radius of convergence due to the presence of non-hydrodynamic excitations. Likewise, the theory of elasticity of brittle solids, organized as an effective field theory in the strain gradients, has zero radius of convergence due to the process of the thermal nucleation of cracks. Viscoelastic materials share properties of both fluids and solids. We use holographic gauge theory/gravity correspondence to study all order hydrodynamics of relativistic viscoelastic media.


Introduction
The applications of the gauge theory/string theory correspondence [1,2] towards nonequilibrium dynamics of relativistic plasma led to a modern perspective on hydrodynamics [3][4][5] as an effective field theory (EFT) where the role of higher-dimensional operators and their coupling constants is being played by higher-order gradient combinations of the local fluid velocity field u µ , in irreducible representations of the symmetry group of the theory, and the transport coefficients correspondingly. For example, in the case of uncharged conformal fluids in background metric g µν , the local stress-energy tensor T µν takes the form [3] T µν = ǫ u µ u ν + P (ǫ) ∆ µν + ∞ i=1 Π µν (i) , ∆ µν = g µν + u µ u ν , (1.1) where ǫ and P (ǫ) are the local energy density and the pressure, and the Π µν (i) collects all operators involving a total of i derivatives of the local velocity and the metric. Explicitly, • for i = 1: where η is the shear viscosity; • for i = 2: where τ Π is the shear relaxation time, and κ and λ 1···3 are the four additional second-order transport coefficients. See eq. (3.7) in [3] for the definition of O µν 1···5 .
For non-conformal relativistic fluids there is an additional first-order in the gradients operator, O 1,2 , with the coupling constant being the bulk viscosity ζ, and 15 second-order in the gradients operators with 13 independent transport coefficients [5]. A substantial effort of the community 1 was dedicated to the computation of the transport coefficients from the first principles, using the framework of the holographic correspondence.
While it was clear that from the gauge theory/gravity correspondence it is possible in principle to extend the program of [3,5], to arbitrary order in the velocity gradients, it was not until the work of [6] when this was practically implemented for the boost-invariant flow of N = 4 supersymmetric Yang-Mills (SYM) theory. Using the first 240 terms of the gradient expansion, the authors observed a factorial growth of gradient contributions at large orders, which indicated a zero radius of convergence of the hydrodynamic series. Furthermore, they identified the leading singularity in the Borel transform of the hydrodynamic energy density with the lowest nonhydrodynamic excitation corresponding to a 'nonhydrodynamic' quasinormal mode on the gravity side.
A related conclusion was reached in the analysis of the nonlinear elastic theory of brittle materials [7,8]. Consider a 2D brittle material with a Young's modulus Y , a Poisson ratio σ, and a 'brittle crack surface tension' α. At a finite temperature T , due to thermal nucleation of cracks, the free energy F of the material subject to the external uniform tension P (negative compression) develops an essential singularity, where A is the elastic material area and λ is the ultraviolet cutoff in the theory (roughly, the interatomic distance). Assuming that the free energy is an analytic function in the complex P plane except for a branch cut along P ∈ (−∞, 0], a Cauchy representation for (1.5) implies that the power series for the inverse bulk modulus i.e., , the high-order terms c n in the perturbative expansion for the inverse bulk modulus roughly grow as (n/2)!. For brittle 3D elastic materials the scaling is c n+1 /c n ∼ n 1/4 as n → ∞. We should stress that unlike the boost-invariant N = 4 SYM conformal relativistic hydrodynamics, where the higher-order terms were explicitly calculated and the essential singularity in the Borel-resumed expression was identified with the nonperturbative effects (the nonhydrodynamic modes in plasma), in the nonlinear elastic theory, the brittle non-perturbative effects (cracks) were conjectured to imply the zero radius of convergence of the perturbative in strain expansion.
Fluids and solids are different. Most notably, fluids lack the shear elastic modulus, and as a result the transverse sound modes in fluids are non-propagating (purely dissipative). Viscoelastic materials are rather interesting as they share both the properties of fluids and solids: they flow (and thus one can formulate for them the hydrodynamic theory) and they also exhibit elastic characteristics when undergoing deformation. Recent advances in engineering holographic viscoelastic materials [9][10][11][12][13] opened a possibility to explore large-order in the velocity gradients hydrodynamics of materials with a control parameter that smoothly interpolates between the fluids and solids. This is precisely the focus of the paper: following the set-up of [14,15] we study the asymptotic properties of viscoelastic hydrodynamics for the holographic models inspired by [12] and the earlier work [16]. En passant we analyze, within linear response, the viscoelastic properties of the model and in particular the elastic moduli and the viscosities and we make contact with the results from the QNMs spectrum.
In the next section 2 we introduce the holographic model describing the viscoelastic media. Translational invariance is explicitly broken as in [16]. Typical holographic observables, the one-and two-point correlation functions, 'feel' this explicit symmetry breaking through a messenger sector experiencing explicit or spontaneous flavor symmetry breaking. Next, we discuss the thermodynamics of the model in section 3.
We follow in section 4 highlighting the elastic features of our holographic media. In section 5 we discuss all-order hydrodynamics of the model undergoing homogeneous and isotropic expansion. We conclude in section 6. We provide more details about the model and the computations in the appendices A, B, C.
Our effective action is defined by where ∆ is the scaling dimension of the boundary operator dual to the bulk field φ and the bulk coupling constants satisfy λ i > 0, γ > 0. The action (2.1) can be thought as a generalization of an effective Maxwell-Einstein-Hilbert action with three complex scalar fields with SU(3) symmetry (see App.A for more details about it).
In this paper we discuss states of the holographic viscoelastic media (2.1) with explicit breaking of the translational invariance, i.e., we turn on the non-normalizable components of the 'axions' ψ I : where x J=1···3 are the spatial directions of the boundary gauge theory. We retain the SO(3) invariance in the spatial coordinates.
The ansatz (2.2) implies that all the spatial coordinates are periodically identified: thus, we have a cubic spatial lattice. The remaining gravitational fields are taken to depend only on the AdS radial coordinate x and the time t: The AdS boundary is at x → 0 + , and the background metric of the dual gauge theory is taken as The equations of motions obtained from (2.1) are shown in Appendix A.
Generically, the geometry (2.4) will have an apparent horizon x AH , located at see [22]. Following [23,24] we associate the non-equilibrium entropy density S with the Bekenstein-Hawking entropy of the apparent horizon in the geometry, Taking the derivative of the entropy density and using the holographic equations of motion we find We conclude this section with the following observations.
• Our system admits a general solution describing the homogeneous and isotropic expansion of a plasma in the generic background metric (2.5): From (2.6), the dynamical (apparent) horizon x AH is located at The corresponding (dynamical) thermal characteristics of the state are where µ is the chemical potential, Q the charge density, T the temperature, S the entropy density, E the energy density and finally P the pressure.
Notice that there is no comoving entropy production: This, along with the conformal anomaly contributions to the energy density and the pressure in (2.11), expresses the statement that the state (2.9) is a Weyl transform of the thermal equilibrium state with a(t) = 1, see [25]. To avoid unnecessary cluttering of the formulas we set 8πG N = 1 henceforth.
• In the static limit a(t) = 1 we obtain the familiar AdS 5 Reissner-Nordström black brane solution where x h is the position of the horizon and Q the charge density of the dual theory.
• Despite the fact that translational invariance is broken explicitly due to (2.2), the thermal equilibrium states of our viscoelastic media (2.9) do not feel this breaking because the messenger field φ vanishes identically. In this case the discrete Z 2 symmetry φ ↔ −φ is unbroken. Asymptotically near the boundary 2 , we have where a constant δ ∆ is a source term (coupling constant of the dual operator), and f ∆ is its expectation value. When δ ∆ = 0, the parity Z 2 symmetry is explicitly broken; it can also be spontaneously broken when δ ∆ = 0 provided the energy density of the state with Q = 0 is sufficiently small. Whenever this Z 2 symmetry is broken, the one-point correlation function of the stress-energy tensor is sensitive to the explicit translational invariance breaking (2.2).

Phase diagram of the holographic viscoelastic media
While it is straightforward to analyze the phase diagram of the equilibrium thermal states of the viscoelastic media holographically represented by (2.1), the task is daunting given the dimensionality of the parameter space of the model: the bulk couplings {γ, λ 1 , λ 2 }, the lattice spacing ∆x ≡ 2π/(k √ 3), the scaling dimension ∆ of the relevant operator O φ , and the boundary conformal symmetry breaking coupling constant δ ∆ .
These parameters have to be supplemented by the temperature T and the chemical potential µ (for a grand canonical ensemble equilibrium states), or the energy density E and the charge density Q (for a microcanonical ensemble equilibrium states). In what follows, we mostly restrict the bulk parameters of the holographic dual to the choice The lattice spacing ∆x is allowed to vary to interpolate from the more-fluid-like behavior k ≪ T to the more-solid-like behavior k ≫ T of our viscoelastic media.
We begin in section 3.1 discussing the linearized stability of the Z 2 -invariant RN state ((2.9) with a(t) = 1) at extremality (T = 0). We continue with the linearized stability analysis in the microcanonical and grand canonical ensembles and determine the critical energy density and temperature. We discuss the fully non-linear phase diagram of the ∆ = 2 and δ 2 = 0 model in section 3.2.

BFology and critical phenomena
The Z 2 invariant thermal equilibrium states of the viscoelastic media represented by the bulk geometry (2.9) in the static limit a(t) = 1, do not feel the explicit translational invariance breaking (2.2). Following [26], we expect that close to the extremality, these states with Q = 0 will become unstable to spontaneous Z 2 symmetry breaking due to the condensation of φ. Of course, the onset of the instability is sensitive to ∆x (or k), where we introduced In the next section we determine the precise critical energy and temperature for the presence of the normalizable mode of (3.3). For now, we determine the boundaries of the gravitational bulk parameter space region that would guarantee the onset of the instability at the extremality of the Z 2 symmetric solution. This "BFology" analysis is done applying the Breitenlohner-Freedman (BF) bound to φ scalar in the infrared, i.e., z → 1, of the extremal RN solution which can be achieved by setting q → √ 6 − .
In the latter case the near horizon (IR) geometry is that of AdS 2 × R 3 with the radius of the AdS being The BF bound implying the onset of the instability reads Some examples of regions of the gravitational bulk parameter space leading to spontaneous breaking of Z 2 symmetry are shown in fig. 1. Since λ i > 0 increasing the (dimensionless) lattice spacing s makes the instability more difficult to appear. On the contrary, the coupling γ enhances the instability.
We move now to identify critical points for the onset of the Z 2 symmetry breaking instability in microcanonical and grand canonical ensembles for our benchmark model (3.1). We allow the lattice spacing ∆x to vary. To this end, we need to find the normalizable solutions of (3.2) at the threshold of instability: where we exploited the linearity of (3.2) to fix the normalizable coefficient at the boundary to 1.
The right panel enhances the entropy differences between the symmetry broken and the symmetric phases at fixed energy density E and the charge density Q.

Spontaneous symmetry breaking and phase diagram
In the previous section we discussed the critical phenomena in the Q 4 of the energy density E and the charge density Q, normalized to its extremal value E 3 Q 4 | extremal = 81 32 , versus the dimensionless ratio S Q involving the entropy density S. The red curve is the equilibrium Z 2 symmetric phase, and the {orange,green,purple} curves represent Z 2 -broken phases with s/s 2,max = {0, 1 3 , 1 2 }. The right panel show the entropy differences between the symmetry broken and the symmetric phases. The onset of the instabilities (the location of the bifurcation points) agree to better than 1.5×10 −3 % with computations presented in section 3.1. The analysis used to determined the phase diagram are standard and we will not detail them here.

Viscoelastic properties
It is convenient to study elastic properties of the model (2.1) in Fefferman-Graham (FG) coordinates system. Thus, we assume the following gravitational ansatz for the equilibrium states: where {A, B, C, a 0 } and φ are functions of the radial coordinate only. The corresponding EOMs are given in Appendix A.
To determine the retarded correlation function G R x 1 x 2 ,x 1 x 2 (ω) that encodes the shear elastic modulus and the shear viscosity of the media [31] we study the following timedependent metric fluctuation where C is an unperturbed metric warp factor, see (4.1). It is straightforward to verify that fluctuations (4.2) decouple from all the other fluctuations at the linearized level.
Assuming the usual Fourier expansion h = e −iωt H(x) we find We look for a solution of (4.3) normalized to 1 near the AdS 5 boundary, x → 0, and representing an incoming wave at the horizon, x → x h . From the near-boundary asymptotic expansion of the solution, the corresponding retarded correlation function reads and it determines the shear elastic modulus G and the shear viscosity 4 η as 4.1 Viscoelastic properties to leading order in Z 2 symmetry breaking The Z 2 symmetry breaking of the viscoelastic media state (either explicit or spontaneous) is mediated to the gravitational bulk scalar φ. In general, to find the shear elastic modulus G and the shear viscosity η (4.6) one has to solve (4.3) numerically 5 .
Semi-analytic computations are possible to leading order in the symmetry breaking, i.e., to order O(φ 2 ). We outline the details of these computations in Appendix B. In the main text we just collect the final expressions for {G, η} to leading order in the Z 2 symmetry breaking: where we used the fact that to this order it is consistent to replace the metric warp factors with those of the Z 2 -symmetric RN solution (see (B.1) for details). Notice that with an explicit Z 2 symmetry breaking G is divergent when λ 1 = 0 and ∆ ≥ 3; it is divergent as well with λ 1 = 0 and λ 2 = 0 provided ∆ = 4. When k = 0, the shear elastic modulus vanishes and the shear viscosity attains its universal ratio [35]. In other words for k = 0 the system is a perfect strongly coupled fluid with no elastic properties and the viscosity saturating the KSS bound [36].
We conclude this section with a sample of results obtained using 6 (4.7). We focus on the case of an explicit Z 2 symmetry breaking of the model with {g, λ 1 , λ 2 } = {1, 0, 1} and ∆ = 2. It is important to remember that because (4.7) was obtained to order O(δ 2 ∆ ), large/divergent results are outside the approximation, and the full nonlinear in 4 In presence of the spontaneous/explicit breaking of translational invariance the definition of the viscosity, its Kubo formula and its relation with the momentum diffusion constant become subtle [19,32,33]. Nevertheless the usual expression in terms of Kubo formula still has a clear physical interpretation as the rate of entropy production due to a strain [34]. 5 We will present an example of this computation in section 4.2. 6 Other cases can be analyzed in a similar fashion. the bulk scalar φ analysis must be performed. We define for simplicity The results of the computations (whose details are given in Appendix B) are presented in figs. 4-5. We use T µ = { 1 12 , 1 6 } (red,green) which are below the critical temperature for the Z 2 spontaneous symmetry breaking (see fig. 1), and T µ = { 1 3 , 2 3 } (magenta,blue) which are above the critical temperature. Black curves correspond to neutral viscoelastic media, i.e., µ = 0. Both the shear elastic modulus and the correction to the shear viscosity become large (divergent in O(δ 2 2 ) approximation) at the critical point. This is shown in 5 for T µ = { 1 12 , 1 6 } -the zeros inĜ −1 agree to better than 1% with the results reported in fig. 2.

Viscoelastic transport of the spontaneously broken Z 2 phase
We discuss here the viscoelastic transport (elastic shear and bulk moduli) and the shear viscosity of the spontaneously broken Z 2 phase in the model (2.1). We focus on the parameter set {γ, λ 1 , λ 2 } = {1, 0, 1} and ∆ = 2 (δ 2 = 0, i.e. no explicit source, since we consider the spontaneous symmetry breaking). The technical details about the computations are given in Appendix C.
In the Z 2 symmetric phase the shear modulus and the viscosity can be simply and they correspond to a fluid state which saturates the KSS bound and has no shear elastic properties.
In fig. 6 we collect the numerical results forG ≡ 16πG N G/k 4 andη ≡ (4πη/S − 1) for select values of T µ as a function of k T . Notice that both the shear elastic modulus and the deviation from the universal shear viscosity vanish in the limit k T → 0; they also vanish for large enough k T , when the Z 2 symmetry is restored (see fig. 2). Red curves correspond to T µ = 1 12 at the critical point, while the green curves correspond to T µ = 1 6 at the critical point.
In fig. 7 we study behavior ofG andη in the vicinity of the critical point. Both quantities vanish at criticality with the same critical exponent: So far we focused on the shear elastic modulus. There is a simple way to compute the bulk elastic modulus of the media, provided that the conformal symmetry is unbroken -of course, this is the case for the phase of the model with the spontaneous Z 2 symmetry breaking discussed in this section. To proceed with the computations, notice that the equilibrium stress tensor of our media is compatible with the one of an isotropic crystal [37,38]: where P is the equilibrium thermodynamic pressure, Θ J are the phononic degrees of freedom which can be identified with the Stueckelberg fields of the translational symmetry breaking; is the usual strain tensor. Notice that, contrarily to a perfect fluid, T I crystal,I = 3P, where P = −Ω, related to the grand potential density. In a CFT the full stressenergy tensor vanishes T µ µ = 0. Taking into account the definition of the energy density T t t = E and the form of the spatial components of the stress tensor (4.11) we immediately obtain E = 3 P + 9 K the final formula for the bulk modulus: (4.14) All the physical observables in the previous formulas can be obtained explicitly using holographic renormalization (see Appendix C for details). The numerical results for the reduced bulk elastic modulusK = 16πG N K/k 4 are collected in fig. 8. Left panel presents results as a function of k T for select values of T µ = { 1 12 , 1 6 }, {red,green} curves, at the criticality; right panel demonstrates the critical behavior of the modulus at fixed k T = 5. The critical exponent for the bulk elastic modulus α K in our model is

Viscoelastic transport with explicit Z 2 symmetry breaking
Another interesting regime of the model is the explicit breaking of the Z 2 symmetry at fixed value of the symmetry breaking parameter δ ∆ relative to the equilibrium temperature T , but with varying lattice spacing ∆x (or equivalently k T ). To limit the parameter space, we set the chemical potential µ = 0 and further restrict This regime is represented holographically by the same gravitational bulk equations where we introduced a dimensionless parametrization 7 of δ 2 > 0 as Once again, the background has to be determined numerically. In fig. 9 we present the behavior of the scalar field at the horizon p h,0 , and p 2 (related to the expectation value of the dual O 2 operator) at fixed δ 2 T 2 as a function of k T . Notice that in the limit k T → ∞ (vanishing lattice spacing) the scalar is exponentially suppressed at the horizon, while the expectation value of the corresponding dual operator is finite. As a result, in this limit, the horizon of the gravitation background is 'hairless' Schwarzschild, leading to the universal result for the shear viscosity [35]. On the contrary, finiteness of p 2 as 7 Comparing with [39] δ 2 ∝ +m 2 > 0.
Notice that whileη → 0 as k T → ∞,Ĝ remains finite. The fluid-gravity correspondence [4] identifies the effective theory of small fluctuations of the holographic horizon with the hydrodynamics of the boundary plasma.
In the limit k T → ∞ the horizon fluctuations are sensitive to the translational symmetry breaking only through the mediator -the bulk scalar φ, which is vanishingly small. Thus, we expect that hydrodynamics (at least for the first few orders in the

Hydrodynamics of homogeneous and isotropic flows
In this section we study the all-order hydrodynamics of the viscoelastic media undergoing the homogeneous and isotropic expansion following [14,15]. We focus on the explicit Z 2 flavor symmetry breaking to leading nontrivial order in the source, i.e., O(δ 2 ∆ ). We omit technical details and refer the interested reader to the earlier work.
There are many probes of the all-order hydrodynamics: expectation values of the correlation functions, entanglement entropy, etc. We consider the entropy production of the viscoelastic plasma in the background metric (2.5). To compute the entropy 8 Sound channel fluctuations are more involved and will not be discussed here.
production (2.8) to order O(δ 2 ∆ ) it is sufficient to study the linearized, order O(δ ∆ ), scalar field dynamics in (2.9): leading to the comoving entropy production rate A general solution to (5.1) can be written as a series expansion in the derivatives of the scalar factor a(t), where the functional T n,∆ [a] involves n τ -derivatives of the scale factor, and T n,∆ = aṪ n−1,∆ + (4 − ∆)ȧT n−1 , The all-order entropy production expression (5.3) then takes form Notice that to leading order in the derivative expansion, which when compared with the entropy production for the homogeneous and isotropic flow in the hydrodynamic approximation identifies the bulk viscosity of the viscoelastic medial to order O(δ 2 ∆ ) as The reduced bulk viscosityζ ∆ =ζ ∆ (γ, λ 1 , λ 2 , s, q) is presented in fig. 12 for the viscoelastic models with ∆ = {2, 3} and {γ, λ 1 , λ 2 } = {1, 0, 1} for select values of Notice that the bulk viscosity vanishes in the large-s (small lattice spacing) limit. Again this is a reflection of the fact that the hydrodynamic transport is determined by the near horizon regime of the gravitational dual, and the bulk scalar field is exponentially suppressed at the horizon in this limit, see fig. 9. Much like the shear elastic moduluŝ G and the correction to the shear viscosity κ (see fig. 4), the bulk viscosity diverges in the phase with explicit Z 2 symmetry breaking close to criticality.
In the rest of this section we present results for all-order computation of Ω ∆ for the viscoelastic model with ∆ = 2. We consider {γ, λ 1 , λ 2 } = {1, 0, 1} neutral viscoelastic medial in section 5.1, and charged plasma with spontaneous symmetry breaking in section 5.2.

Charge plasma with spontaneous symmetry breaking
We now consider the 'charged fluid/plasma' regime of the viscoelastic media: γ = 1 and s = 0 but with q = 0. Specifically, we choose q = 2.24, corresponding to (see fig. 2)

Conclusions and discussion
The goal of this project was understanding the large-order perturbative expansion of a viscoelastic media with a control parameter smoothly interpolating between the fluid and the solid. Previously, both the all-order fluid hydrodynamics [6] and the all-order theory of elasticity [7] were argued to be asymptotic expansions (in velocity gradients 10 Similar phenomenon was observed in [15]. Physically, the intuitive picture of our holographic viscoelastic media is that of the crystalline lattice of spacing ∆x, immersed in a fluid: as ∆x → ∞ the lattice is removed and one is left with the fluid only. We presented a rather detailed exploration of the media (2.1): • We studied the equilibrium phase diagram of the model with both explicit and spontaneous flavor symmetry breaking. We identified the critical behavior and the instabilities in the microcanonical and grand canonical ensembles. As common in holographic models, the critical behavior is mean-field.
• We showed that either with explicit or spontaneous flavor symmetry breaking in the model there is a nonzero shear elastic modulus G when the lattice spacing ∆x is finite. Furthermore, the shear elastic modulus vanishes in the 'fluid limit': lim ∆x→∞ G = 0 (6.1) • In case of the spontaneous symmetry breaking, we utilized the underlying conformal symmetry of the model to show that there is a nonzero bulk elastic modulus K as well, when ∆x is finite. Both the shear and the bulk elastic moduli vanish at the criticality, G ∼ (T crit − T ) α G and K ∼ (T crit − T ) α K , with the same critical exponent α G = α K = 1. It is clear that these critical exponents are not universal, and instead depend on the coupling of the condensate to the axion sector of the model. Specifically, modifying the coupling of the bulk scalar φ to the axions ψ I from φ 2 to φ 2n is expected to produce critical exponents α G = α K = n. The reason for this is simple: close to critically φ ∼ (T crit − T ) 1/2 , characteristic of the mean-field behavior, and the elastic moduli respond proportional to the coupling strength, i.e., ∼ φ 2n .
• In case of the explicit flavor symmetry breaking we found that the shear elastic modulus G → const as ∆x → 0, indicative of the 'solid limit'.
• We analyzed the hydrodynamic transport of the model: the shear viscosity (extracted via a Kubo formula from the retarded two-point correlation function of the stress-energy tensor) and the bulk viscosity. The shear viscosity violates the bound η s ≥ 1 4π both for the explicit and the spontaneous flavor symmetry breaking at finite ∆x. This violation vanishes in the fluid limit, and at the criticality (in case of the spontaneous symmetry breaking). Somewhat unexpectedly, both the violation of the shear viscosity bound and the bulk viscosity vanish in the solid limit as well, i.e., when ∆x → 0. The physical reason is unclear to us, but at a technical level this is related to the exponential suppression of the bulk scalar field φ in our model at the horizon of the gravitational dual in the limit ∆x → 0.
Since the hydrodynamic transport is fully determined by the small fluctuations of the horizon, in the solid limit, the horizon is scalar-hair free, leading to the universal shear viscosity and the vanishing bulk viscosity. Note that the bulk scalar φ is not suppressed near the AdS boundary in the solid limit, which is likely responsible for the nonvanishing of the elastic shear modulus in this limit.
• Although the hydrodynamic transport (shear and bulk viscosities) are that of the conformal fluid in the 'solid limit' of our viscoelastic media, the nonhydrodynamic excitations feel the lattice. In the limit ∆x → 0 only the modes with vanishing spatial momentum remain in the spectrum. Modes in the scalar/shear channels (according to classification of [40]) experience a small but nonzero off-set, see There are many open questions left for future analysis: First and foremost, we did not identify the cracks present within the theory of elasticity for brittle materials. We don't know whether our holographic viscoelastic media is brittle. The wall-of-poles on the Borel plane hints to additional singularity structure(s), beyond the poles of the Pade approximant.
We focused on one observable, i.e., the nonequilibrium entropy production rate.
What are the large-order hydrodynamic expansions of one-and higher-point correlation functions and of the entanglement entropy?
To leading order in the explicit symmetry breaking, the elastic moduli and the transport coefficients diverge as one approaches the criticality. The observed divergences (see figs. 4-5 for example) are outside the regime of the validity of the probe approximation and the full nonlinear analysis is necessary.
We are not aware of any discussion of the bulk elastic modulus in holographic nonconformal models 11 .
To facilitate the computations, we restricted the all-order hydrodynamic analysis to a particular flow: the homogeneous and isotropic expansion. It would be interesting to analyze hydrodynamic perturbation theory of the boost-invariant expansion.
We studied nonhydrodynamic excitations in our viscoelastic media in the scalar/shear channels, and for the condensate. The analysis of the QNMs in the sound channel is left for the future work.
It would be interesting to analyze the fate of the bulk viscosity bound [42] in the presence of the elastic features.
Quenching the holographic media and study of its thermalization, as done for ex- 11 With the exception of a formula given, but not studied in any direction, in [12].
ample in [43][44][45], could provide further insights about its viscoelastic features and in particular the possible presence of slow relaxation and aging typical of glassy systems (see [46,47] for some early attempts).
Lastly, it is important to explore different holographic viscoelastic models [13,18] within the context of all-order hydrodynamic expansion. A particularly interesting case is represented by those models where the shear viscosity vanishes in the 'solid limit' [19,32,34].

Acknowledgments
We A Details about the model and its EOMs We start with an effective Maxwell-Enstein-Hilbert action in asymptotically AdS 5 geometry with three complex scalar fields Φ I enjoying SU(3) flavor symmetry: where ∆ is a scaling dimension of the dual boundary operator, and the bulk coupling constant obeys γ > 0. Introducing we rewrite (A.1) with manifest U(1) 3 flavor symmetry Notice that a consistent truncation of (A.3) with ψ I ≡ 0 leads to an effective action S sc of a neutral holographic superconductor: where the operator O φ , dual to the bulk mode φ, has a conformal dimension ∆.
The effective action we use throughout the paper is a simple generalization of (A.3).
As in [26], in the absence of the source term for a real bulk scalar φ, spatially homogeneous and isotropic thermal charged states of (A.4) will become unstable with respect to a φ-condensate (below some critical temperature or the energy density); the role of the bulk coupling constant γ > 0 is to stimulate this condensation process, i.e., increase the critical temperature/energy density.
The equations of motion which follows from the action (2.1) presented in the main text are, d + ≡ ∂ t − x 2 A∂ x and ′ ≡ ∂ x , · ≡ ∂ t , Additionally, there are 'constraints':

EOMs in Fefferman-Graham coordinates
The EOMs in the FG coordinates (4.1) are given by (1 + γφ 2 ) (A. 15) and the first-order constraint, related to reparametrization of the radial coordinate, 16) B Shear elastic modulus and shear viscosity at leading order Notice that to leading order in O(φ 2 ), both the metric warp factors {A, B, C} and the bulk gauge potential a 0 are corrected. However, as we show, the viscoelastic transport is completely determined by the corresponding quantities without the bulk scalar backreaction, i.e., the Z 2 symmetric or the Reissner-Nordström black brane solution: where x h is the location of the horizon, related to temperature T , and µ is the chemical potential.
To extract {G, η} we need to solve (4.3) to order O(ω), subject to the boundary condition (4.4) and an incoming wave boundary condition at the horizon. It is convenient to represent where the two functions satisfy, correspondingly 12 where H 0 and H 1 has to be solved to order O(φ 2 ) -the latter guarantees that (B leading to (see (4.4)-(4.6)) i.e., to this order, elastic shear modulus vanishes, and the shear viscosity is universal, as first established at this level of generality in [48]. Having