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 ( ) 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 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 10 independent transport coefficients [5,6]. 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 [7] 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 [8,9]. 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 JHEP03(2019)146 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 is an asymptotic expansion with c n+1 c n → −n 1/2 πT (1 − σ 2 ) 8Y α 2 1/2 , as n → ∞ (1.7) 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 non-perturbative effects (the nonhydrodynamic modes in plasma), in the nonlinear elastic theory, the brittle nonperturbative 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 [10][11][12][13][14] 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 [15,16] we study the asymptotic properties of viscoelastic hydrodynamics for the holographic models inspired by [13] and the earlier work [17]. 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 [17]. 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.

JHEP03(2019)146
Our effective action is defined by 2 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 appendix 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' 3 ψ 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. As a definition of the model, we assume the periodicity of ψ I in the field space as in eq. (A. 2) 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 We do not have a string theory embedding of the model. As a result, the boundary interpretation of a higher-derivative bulk coupling λ2 is unclear. 3 The fields ψ I are genuine axions only when λ1 = 1 and λ2 = 0. Thinking about these fields as axions allows for an intuitive understanding why the holographic media shares viscoelastic properties. As we demonstrate explicitly later, the media enjoys nonzero shear elastic modulus for all {λ1, λ2}.

JHEP03(2019)146
see [23]. Following [24,25] we associate the non-equilibrium entropy density S with the Bekenstein-Hawking entropy of the apparent horizon (AH) in the geometry, 4 Taking the derivative of the entropy density and using the holographic equations of motion we find (2.8) 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 (2.10) The corresponding (dynamical) thermal characteristics of the state are 3 (2.11) 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 [26]. To avoid unnecessary cluttering of the formulas we set 8πG N = 1 henceforth.

JHEP03(2019)146
• 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, 5 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 JHEP03(2019)146 invariance breaking (2.2). Following [27], 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), as well as to the bulk couplings {γ, λ 1 , λ 2 }. Indeed, the linearized φ-fluctuations about the RN background satisfy 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 L 2 2 = L 2 12 = 1 12 (3.4) and the effective IR mass of φ: 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 figure 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: JHEP03(2019)146  where we exploited the linearity of (3.2) to fix the normalizable coefficient at the boundary to 1.
Numerical results for the onset of the instability are presented in figure 2 (in both microcanonical and grand canonical ensembles). The onset is pushed to the extremality once s reaches the boundary of the instability region, see (3.6): (3.8)

Spontaneous symmetry breaking and phase diagram
In the previous section we discussed the critical phenomena in the  (a) dominates the symmetric phase both in microcanonical and grand canonical ensembles (as in [27]) (b) does not dominate the symmetric phase in either ensembles (as in some models in [28,29]) (c) dominates the symmetric phase in microcanonical ensemble, but not in a canonical ensemble (as in some models in [30,31]) In our model the case (a) is realized. In figure 3 we present the phase diagram of the model (2.1) with {γ, λ 1 , λ 2 } = {1, 0, 1}, ∆ = 2 and δ 2 = 0 in microcanonical ensemble. 6 The left panel represent the dimensionless ratio E 3 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:

JHEP03(2019)146
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 [32] 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 7 η as 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 [20,33,34]. 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 [35]. For completeness let us just summarize the main results of [35]. Given an external strain source δgxy = t∆, where ∆ is a constant, the corresponding entropy production induced by the work of such external source is obtained as: which very suggestively can be rewritten as: Whenever temperature is the only scale in the system it is natural to expect ∆/T being order unit. Most importantly the previous equation can be used to give a more universal definition of the η/s ratio which does not rely on hydrodynamics. Moreover we can rewrite the KSS bound as a bound on the entropy production of the system: where t pl is the Planckian time t pl ≡ /kBT .

JHEP03(2019)146
4.1 Viscoelastic properties to leading order in Z Z 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.9) one has to solve (4.3) numerically. 8 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 [36]. 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 [37].
We conclude this section with a sample of results obtained using 9 (4.10). 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.10) was obtained to order O(δ 2 ∆ ), large/divergent results are outside the approximation, and the full nonlinear in 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 figures 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 figure 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 figure 2. Note that κ > 0 implying the violation of the KSS bound. 10 8 We will present an example of this computation in section 4.2. 9 Other cases can be analyzed in a similar fashion. 10 The violation of the KSS bound in anisotropic setting has been observed earlier [38].  (4.11) for the parametrization). The color coding corresponds to different values of the ratio of temperature T to the chemical potential µ: T µ = 1 12 (red), T µ = 1 6 (green), T µ = 1 3 (magenta) and T µ = 2 3 (blue). The former two values are below the critical values of T /µ| crit at s = 0, and the latter two are above. Black curves correspond to neutral viscoelastic media, i.e., µ = 0.

Viscoelastic transport of the spontaneously broken Z Z 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 obtained JHEP03(2019)146  and they correspond to a fluid state which saturates the KSS bound and has no shear elastic properties. In figure 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 figure 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 figure 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 JHEP03(2019)146 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 [39,40]: 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 stress-energy 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.14) we immediately obtain E = 3 P + 9 K (4. 16) which in absence of elastic properties is the usual conformal relation E = (d − 1)P.
Using the Smarr relation E + P = ST + µQ and combining it with (4.16) we obtain the final formula for the bulk modulus: At this point is important to notice that the bulk modulus defined in (4.17) is just the "solid" contribution, which indeed vanishes in the pure RN "fluid" solution. On the contrary, the bulk modulus (intended as the total one) is nonzero even in the fluid phase and it gives rise to the common finite speed of longitudinal sound, i.e. c 2 L = K/( + p), which has already been analyzed in the literature. In other words, we can think of K in (4.17) as the additional contribution to the bulk modulus, and therefore the longitudinal sound speed, due to the solid nature of the media. 11 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 figure 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 Z 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 (C.1)-(C.5) as in section 4.2, albeit with the vanishing bulk gauge potential a(z) ≡ 0 . (4.20) The details of the setup and the numerical computations are presented in appendix C. The AdS boundary condition for the scalar which reflects a δ 2 = 0 source term of the explicit Z 2 symmetry breaking are: where we introduced a dimensionless parametrization 12 of δ 2 > 0 as Once again, the background has to be determined numerically. In figure 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 JHEP03(2019)146  result for the shear viscosity [36]. On the contrary, finiteness of p 2 as k T → ∞ suggests that the shear elastic modulus remains finite. These expectations are indeed supported by direct computations. Figure 10 collects the numerical results for 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 gradient expansion) should not feel the spatial lattice of our viscoelastic media. The non-hydrodynamic plasma excitations should know about the translational symmetry breaking. Following [42] we study the quasinormal modes in the shear/scalar channels of the holographic dual. 13 Because of the spatial lattice, the momentum q of the modes must be quantized. All the modes except JHEP03(2019)146 Figure 11. Lowest non-hydrodynamic mode at zero spatial momentum in the shear/scalar channel at fixed δ 2 /T 2 . We defineω = ω/(2πT ) the reduced frequency. those with q = 0 become infinitely heavy in the limit of vanishing lattice spacing, and decouple. Thus, we focus on the q = 0 sector. In this sector there is no distinction between the scalar channel and the shear channel QNMs. Numerical results for the lowest nonhydrodynamic QNM are collected in figure 11. As expected, here (4.24)

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 [15,16]. 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 production (2.8) to order O(δ 2 ∆ ) it is sufficient to study the linearized, order O(δ ∆ ), scalar field dynamics in (2.9):

JHEP03(2019)146
leading to the comoving entropy production rate T (t) is a local temperature, red-shifting as T (t) = T 0 a(t) , T 0 ∼ 1 x h . 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 The recursive relation for T n,∆ can be solved analytically for simple scale factors -e.g., choosing de Sitter expansion, a(t) = e Ht = exp(Hx h τ ), we find T n,∆ = Γ(n + 4 − ∆)(Hx h ) n a n Γ(4 − ∆) .

(5.7)
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 14

JHEP03(2019)146
identifies the bulk viscosity of the viscoelastic medial to order O(δ 2 ∆ ) as The reduced bulk viscosityζ ∆ =ζ ∆ (γ, λ 1 , λ 2 , s, q) is presented in figure 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 figure 9. Much like the shear elastic modulusĜ and the correction to the shear viscosity κ (see figure 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.

Neutral viscoelastic media
Setting λ 1 = q = 0 and γ = λ 2 = 1, we can solve equation (5.5) for F 0,2 and F 1,2 analytically: 15 We solved numerically (5.5) for n = 2 · · · 300, Borel transformed Ω ∆=2 → Ω (B) ∆=2 , Pade approximated the result, extracted the poles and compared them with the appropriate QNMs [15,16]. Figure 13 (left panel) reproduces the result reported in [15] for the agreement between the QNMsω      of-poles are in excellent agreement with the corresponding QNM frequencies. Note that as s increases, the Borel plane singularities to the right of the wall move closer to the imaginary axes (the correspondingω QNM have vanishingly small imaginary part). This is emphasized in the right panel, where the ratio Re[ξ 0 ]/ Im[ξ 0 ] varies from ∼ 10 −19 to ∼ 10 −6 for the poles closest to the origin to the one furthest away. In the limit s → ∞ we expect the 'solid' properties of the viscoelastic media to be enhanced: the QNMs on the branch corresponding to the gravitational bulk scalar fluctuations approach 'normal modes' -the nearly non-dissipative fluctuations characteristic of those of the perfect lattice.

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 figure 2) Since we are below the critical temperature for the spontaneous Z 2 symmetry breaking, the plasma is unstable: there is a branch of the QNMs with Re[ŵ QNM ] = 0. Some of the modes on this branch have Im[ŵ QNM ] > 0, signaling the perturbative instability. Of course, there is also the branch of the QNMs with complex frequencies -this is µ = 0 deformation of the neutral fluid/plasma QNM frequencies computed in [43]. The QNM frequencies at q = 2.24 are represented by the red crosses, while those at q = 0 (µ = 0) are represented by the green crosses in figure 17. The orange lines represent the 'flow' as q ∈ [0, 2.24] of the complex QNM frequencies. As before, solid blue circles represent the Borel plane singularities. There is a good agreement only with the first two QNMs on the Re[ŵ QNM ] = 0 branch. One of these two modes is unstable -we see that the Borel resummation of the hydrodynamic derivative expansion captures the plasma instabilities. . Green crosses represent QNM frequencies at µ = 0, red crosses represent QNM frequencies at corresponding T µ . Orange lines represent the 'flow' of the QNM frequencies as µ changes from zero. between the QNM frequencies and the Borel plane singularities close-to or to the left of the wall. As we turn on s = 0 (at fixed q = 0), viscoelastic media instabilities are suppressed, provided s is sufficiently large (see figure 2) -we observed that in this case all the QNMs on Re[ŵ QNM ] = 0 branch are stable. While we did not study this in details, we believe that the wall-of-poles will move to Re[ξ 0 ] → −∞ as q → 0, revealing better agreement with the Borel plane singularities (for small enough q the QNM branch with Re[ŵ QNM ] = 0 will disappear). It would be interesting to better understand why there is worse agreement between the Borel plane singularities and the QNM frequencies for the charged viscoelastic media. This might be related to the present of additional singularities 16 in (5.5) at q = 0 (also for λ 1 = 0), which renders Pade approximation not reliable for the high-order poles of the Borel transform of Ω ∆ -i.e., there might be additional singularities (besides the poles) of the Borel transform.

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 [7] and the all-order theory of elasticity [8] were argued to be asymptotic expansions (in velocity gradients and strain gradients correspondingly). The nonperturbative effects responsible for the zero radius JHEP03(2019)146 of convergence of the perturbative expansions were argued to be distinct in both cases: these are the QNMs of the fluid and cracks in brittle solids. To analyze the problem in a controlled setting, we embedded it in the framework of holographic gauge theory/gravity correspondence. Specifically, we focused on a class of viscoelastic holographic models introduced in [17]: the condensate of the holographic superconductor is coupled to the axion sector with a flavor symmetry; turning on a source term for the axions necessitates the introduction of the spatial lattice, explicitly breaking the translational invariance; once the axion flavor symmetry is spontaneously (or explicitly) broken, the homogeneous and isotropic states of the viscoelastic media 'feel' the breaking of the translational invariance, leading to the a nonzero shear elastic modulus -an important differentiating feature between fluids and solids; controlling the scale of the source term for the axions (inversely proportional to the spatial lattice spacing ∆x) we can interpolate between the fluid and the solid regimes of the viscoelastic media.
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, 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 .

JHEP03(2019)146
• 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 stressenergy 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 [42]) experience a small but nonzero off-set, see (4.24). The nonhydrodynamic modes of the condensate are profoundly affected: the complex frequencies on this branch have exponentially suppressed imaginary part compare to the real part, see the right panel of figure 16. They represent the nearly nondissipative fluctuations characteristic of those of the perfect lattice.
• We studied the homogeneous and the isotropic expansion of the viscoelastic media to all orders in the velocity gradients, but to leading order in the explicit flavor symmetry breaking. We focused on a particular observable: the non-equilibrium entropy production rate. For neutral viscoelastic media we found a good agreement between the poles of the Pade approximant of the Borel transform of the observable, and the corresponding QNM frequencies. The Borel plane singularities closest to the origin follow the flow of the QNM with ∆x, correctly reproducing the solid limit with nearly nondissipative nonhydrodynamic excitations of the flavor symmetry breaking condensate, see figure 16. We observed a W -shaped wall-of-poles, which presence hampers the agreement between the Borel plane singularities and the QNM frequencies. We observed the lingering effects of the wall in the fluid limit, i.e., ∆x → ∞, see figure 14. Unfortunately we could not establish whether the presence of the wall indicates a new nonperturbative effect contributing to the asymptotic character of the viscoelastic gradient expansion, or is simply an artifact of the Pade approximation. We did establish that the presence of the wall is robust with respect

JHEP03(2019)146
to the truncation order in series expansion (5.6). For charged viscoelastic media there is also a wall-of-poles, albeit a different one -in this case the agreement between the Borel plane singularities and the QNM frequencies is limited to 1-2 modes, closest to the origin. We established that the all-order derivative expansion for the charged fluid correctly identifies nonhydrodynamic instabilities associated with the spontaneous flavor symmetry breaking, see figure 17.
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 figures 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. 17 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 [44] in the presence of the elastic features.
Quenching the holographic media and study of its thermalization, as done for example in [45][46][47], could provide further insights about its viscoelastic features and in particular the possible presence of slow relaxation and aging typical of glassy systems (see [48,49] for some early attempts).
Lastly, it is important to explore different holographic viscoelastic models [14,19] 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' [20,33,35].

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 ∆.

JHEP03(2019)146
The effective action we use throughout the paper is a simple generalization of (A.3). As in [27], 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, Additionally, there are 'constraints': Once the equation (A.9) is satisfied at some time t = 0, by virtue of (A.5)-(A.8) it is satisfied at any time t > 0. Similarly, once equations (A.10)-(A.11) are satisfied at any radial location, e.g., the AdS 5 boundary at x = 0, for all times, they are satisfied at any other radial location. These equations represent the conservation of the energy density (eq. (A.10)) and the charge density (eq. (A.11)).
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, (A.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 18 The term O(ω 2 ) is important to properly set up the incoming wave boundary conditions, but can be neglected otherwise.

C The spontaneously/explicitly broken phases
The spontaneously broken phase. We collect here the details about the results presented in section 4.2 dealing with the spontaneously broken phase and its transport properties.