Viscoelastic Dynamics in Holography

TL;DR

This paper investigates the viscoelastic response in holographic models, revealing that the traditional Maxwell relaxation time often misestimates the actual relaxation time, especially at low temperatures, and exhibits overshoot phenomena.

Contribution

It provides explicit calculations of relaxation times in holographic models, showing deviations from the Maxwell relaxation time and analyzing temperature-dependent behaviors.

Findings

At high temperatures, $ au_M$ slightly overestimates the actual relaxation time.

At low temperatures, $ au_M$ underestimates the relaxation time significantly.

The real-time response can show overshoot phenomena before relaxation.

Abstract

We study the mechanical response under time-dependent sources of a simple class of holographic models that exhibit viscoelastic features. The ratio of viscosity over elastic modulus defines an intrinsic relaxation time scale -- the so-called Maxwell relaxation time $\tau_M$, which has been identified traditionally with the relaxation time scale. We compute explicitly the relaxation time in our examples and that it differs from $\tau_M$. At high temperatures $\tau_M$ over-estimates the actual relaxation time, although not by much and moreover it still captures reasonably well the temperature behaviour. At sufficiently low temperatures the situation is reversed: $\tau_M$ underestimates the actual relaxation time, in some cases quite drastically. Moreover, when $\tau_M$ under-estimates the real-time response exhibits an overshoot phenomenon before relaxation. We comment on the $T = 0$ limit, where the relaxation is power-law because our models exhibit criticality.

Figures

Equations

$$ϵ(t)=\underset{elastic}{\underbrace{{ϵ}_0\cos (\omega t)\cos \delta }}+\underset{viscous}{\underbrace{{ϵ}_0\sin (\omega t)\sin \delta }}$$

$$W=\pi {\sigma }_0{ϵ}_0\sin \delta$$

$$\sigma (t)={ϵ}_0\left[G^{\prime }\cos (\omega t)+G^{\prime \prime }\sin (\omega t)\right]=G^{\prime }ϵ(t)+\frac{G^{\prime \prime }}{\omega }\dot{ϵ}(t)$$

$$G^{\prime }=\frac{{\sigma }_0}{{ϵ}_0}\cos \delta ,G^{\prime \prime }=\frac{{\sigma }_0}{{ϵ}_0}\sin \delta .$$

$$\sigma (\omega )=G(\omega )ϵ(\omega ),\mathrm{with}G=G^{\prime }+i{G}^{\prime \prime }$$

$$G(\omega )\equiv {\mathcal{G}}_{T_{xy}{T}_{xy}}^{\text{(R)}}(\omega )$$

$$\mu \equiv G^{\prime }(0),\eta \equiv -\underset{\omega \to 0}{\lim }\left(\frac{G^{\prime \prime }(\omega )}{\omega }\right)$$

$$G(\omega )=\mu -i\omega \eta +\mathcal{O}\left({\omega }^2\right).$$

$$\sigma (t)+\frac{\eta }{G_{\infty }}\dot{\sigma }(t)=\eta \dot{ϵ}(t)(\mathrm{Maxwell}\mathrm{model})$$

$$G(\omega )=\frac{\eta \omega }{i+\eta /G_{\infty }\omega }$$

$${\omega }_M=-i\frac{G_{\infty }}{\eta }.$$

$${\tau }_M=\frac{\eta }{G_{\infty }}.(\mathrm{Maxwell}\mathrm{relaxation}\mathrm{time})$$

$${\tau }_{rel}=-1/Im\left({\omega }_{QNM}^{\ast }\right).$$

$$S=M_P^2\int d^{4}x\sqrt{-g}\left[\frac{R}{2}+\frac{3}{{\ell }^2}-m^{2}V(X,Z)\right]$$

$$d{s}^2=\frac{1}{u^2}\left[-f(u)d{t}^2-2dtdu+d{x}^2+d{y}^2\right]$$

$$f(u)=u^3{\int }_u^{u_h}dv\left[\frac{3}{v^4}-\frac{m^2}{v^4}V(v^2,v^4)\right],$$

$$T=-\frac{f^{\prime }(u_h)}{4\pi }=\frac{6-2m^{2}V\left(u_h^2,u_h^4\right)}{8\pi u_h}.$$

$$h_{xy}\left(-\frac{2m^2{V}_X}{f}-\frac{2i\omega }{uf}\right)+h_{xy}^{\prime }\left(\frac{f^{\prime }}{f}+\frac{2i\omega }{f}-\frac{2}{u}\right)+h_{xy}^{\prime \prime }=0$$

$$h_{xy}=h_{xy(l)}(\omega )(1+\dots )+h_{xy(s)}(\omega )u^3(1+\dots )$$

$${\mathcal{G}}_{T_{xy}{T}_{xy}}^{\text{(R)}}(\omega )=\frac{2\Delta -d}{2}\frac{h_{xy(s)}(\omega )}{h_{xy(l)}(\omega )}=\frac{3}{2}\frac{h_{xy(s)}(\omega )}{h_{xy(l)}(\omega )}$$

$$ϵ=\frac{1}{u_h^3}+m^2{\int }_0^{u_h}\frac{V(u^2,u^4)}{u^4}du$$

$$\sigma (t)=Re\left[{\int }_{-\infty }^{\infty }{e}^{-i\omega t}G(\omega )\widetilde{ϵ}(\omega )\right]$$

$$\mathfrak{L}(t,\beta )\equiv \frac{1}{1+e^{2\beta t}}$$

$$\tilde{\mathfrak{L}}(\omega ,\beta )=\sqrt{\frac{\pi }{2}}\delta (\omega )+\frac{i\sqrt{\frac{\pi }{2}}\text{csch}\left(\frac{\pi \omega }{2\beta }\right)}{2\beta }$$

$$\mathfrak{R}(t)\equiv \{\begin{array}{l}t\text{for}t>0\\ 0\text{otherwise}\end{array}$$

$$\tilde{\mathfrak{R}}(\omega )=-\frac{1}{\sqrt{2\pi }{\omega }^2}-i\sqrt{\frac{\pi }{2}}{\delta }^{\prime }(\omega )$$

$$V(X,Z)=Z^2$$

$$\mu =G^{\prime }(0)=0,\frac{\eta }{s}=\frac{1}{4\pi }$$

$$G_{\infty }=\frac{3}{8}ϵ$$

$${\tau }_{rel}=-1/Im\left({\omega }_{QNM}^{\ast }\right)$$

Full text

I Introduction

The difference between a rigid solid and a flowing liquid appears to be very neat and clear in our everyday experience. Despite the apparent simplicity of the question, classifying objects into solids or liquids is less trivial than expected and sometimes it can take considerable time Edgeworth et al. (1984). All natural materials exhibit both viscous and elastic characteristics when undergoing deformation. The interplay of elasticity and viscous behaviour takes the name of viscoelasticity and it represents an important topic not only from the theoretical point of view but especially for applications and technological developments Lakes (1998); Christensen (2003); Gutierrez-Lemini (2013). A description of viscoelasticity from fundamental principles is far from being conclusive. Part of the reason for this is that the hydrodynamic framework is not able to capture the elastic response in a controlled way. On the other hand, there are effective field theory methods which allow to describe elasticity, however dissipation and finite temperature corrections can’t be incorporated in a straightforward way neither. As a matter of fact most of the understanding relies on very empirical and phenomenological models such as the Maxwell model, the Kelvin–Voigt model and extensions thereof (see below).

Holography has recently been established as a valuable tool for hydrodynamics, especially for strongly correlated or strongly coupled fluids Policastro et al. (2002); Son and Starinets (2007); Baier et al. (2008); Hartnoll et al. (2016a). More recently, progress has been made to introduce elastic effects due to the spontaneous breaking of translational invariance and coexisting with the viscous properties implied by the presence of a black hole horizon, i.e. finite temperature. In this direction, a promising direction to describe viscoelastic materials within holography is given by ‘massive gravity’ theories Baggioli and Pujolas (2015); Alberte et al. (2016a)111See Amoretti et al. (2017); Baggioli and Buchel (2018); Andrade et al. (2018); Li and Wu (2018); Grozdanov and Poovuttikul (2018) for other possibilities.. These models break translational invariance while retaining homogeneity and thus allow to deepen the analysis. It has already been shown Alberte et al. (2016a, b, 2018a, 2017) that such theories possess both a finite elastic shear modulus $\mu$ aside from a finite shear viscosity coefficient $\eta$. Additionally, these models contain propagating transverse phonons with a speed of sound that complies with standard elasticity theory Alberte et al. (2018a, 2017).

Beyond the definition and the direct measurement of the viscosity and the elastic modulus, viscoelasticity gives very peculiar imprints on the real time (or finite frequency) response of a system under deformation. A typical experimental test is the so-called dynamic mechanical analysis which consists in applying to the sample a small oscillatory stress $\sigma$ and measuring the resulting strain $\epsilon$. More precisely considering a dynamical load of the form $\sigma(t)=\sigma_{0}\cos(\omega t)$ applied to a viscoelastic sample, the response will be a strain of the form $\epsilon(t)=\epsilon_{0}\cos(\omega t-\delta)$ where $\delta$ is called the loss angle of the material. Expanding the above expression we obtain:

[TABLE]

where the first term is completely in phase with the input while the second term is out of phase. The first term is the elastic non dissipative response; setting $\delta=0$ indeed the stress and the strain are completely in phase. On the contrary the second term is the viscous response and setting $\delta=90$ degrees the stress and the strain are completely out of phase. In order to motivate such statement it is easy to calculate the work done in stressing a material (per unit volume) which is given by $W=\int\sigma d\epsilon$. A short computation shows that the energy loss is:

[TABLE]

which indeed reflects the fact that $\delta$ measures a dissipative contribution, i.e. for $\delta=0$ no energy is dissipated and the process is reversible – the response is purely elastic.

More generically we can say that the elastic response is proportional to the strain itself $\epsilon(t)$ while the viscous response is proportional only to the strain rate $\dot{\epsilon}(t)$. It is indeed standard to assume that a liquid does not react to a static strain deformation.

Sometimes it is convenient to regard the strain as the input and the stress as the output222Notice this is also doable in experimental setups despite it is less under control than fixing the strain and reading the induced stress.. Given this assumption, we can re-write the induced stress as:

[TABLE]

where we defined:

[TABLE]

It is then customary to represent the linear relation between the oscillating stress and strain in Fourier space defining the complex dynamic modulus $G$ as

[TABLE]

The real part $G^{\prime}$ takes the name of storage modulus and the imaginary and dissipative one $G^{\prime\prime}$ of loss modulus333Note that we use the standard notation in condensed matter physics, where prime and double prime on $G$ denote its real and imaginary parts, not derivatives with respect to the argument. Let us also mention that often in the literature, the study of how matter ’flows’ or responds to time dependent strains is referred to as rheology.. In purely elastic materials the response is in phase and therefore ’instantaneous’. In purely viscous fluids, on the contrary, the response has a phase shift of $\delta=90$ degrees which corresponds to $G^{\prime}=0$ (since $\tan\delta=G^{\prime}/G^{\prime\prime}$); viscoelastic materials lie in between. At the same time we can classify the response using the dimensionless ratio $G^{\prime}/|G|=\cos\delta$ between the storage modulus and the absolute value of the complex dynamic modulus. A purely viscous material would have such ratio equal to zero, while in a purely elastic one that number will be maximum, i.e. $G^{\prime}/|G|=1$.

Given the expression (5) it appears clear that the complex moduli $G(\omega)$ can be identified as the (retarded) Green function for the shear stress operator $\sigma\equiv T_{xy}$. Therefore the following identification holds

[TABLE]

and it will be very useful for extracting the complex moduli within the holographic setup. Indeed from the operational point of view the quantity (6) is what we will compute. That said, we can easily identify from the complex modulus the (static) shear elastic modulus $\mu$ and the shear viscosity $\eta$ as:

[TABLE]

which together determine the low frequency behaviour of the complex modulus itself:

[TABLE]

It is illustrative to discuss one of the simplest phenomenological models, the so-called Maxwell viscoelastic fluid model. The response under homogeneous but time dependent source is modelled by assuming that the relation between instantaneous strain and stress follow the linear differential equation

[TABLE]

with $\eta$ the viscosity and ${G_{\infty}}$ a parameter that is immediately recognized as the effective elastic shear modulus at large frequencies, sometimes also referred to as the instantaneous elastic shear modulus. The model then interpolates between fluid behaviour (zero response to a static strain) and solid like behaviour at large frequencies. The fact that liquids behave like solids at frequencies higher than the so-called Frenkel frequency $\omega_{F}\gg 1/\tau_{M}=G_{\infty}/\eta$ is known since long time MacDougall (1947) and it can be easily experienced by diving from high heights.

It is worth showing the complex shear modulus in frequency space,

[TABLE]

which clearly exhibits $G_{\infty}=\lim_{\omega\to\infty}G(\omega)$. Given that this is nothing but a Green function, it is clear that this signals the presence of a quasi-resonance – we shall stick to the naming quasi normal model (QNM) that is more common in the AdS/CFT literature. The Maxwell model contains one single QNM, located at

[TABLE]

Hence this model exhibits relaxation to the steady state with typical time dependence $\sim e^{-t/\tau_{M}}$ with the Maxwell relaxation time

[TABLE]

From our perspective, there are two important lessons to extract from this model. First, that materials (fluid or not) already possess a relaxation time scale as given by (12) and which is important for dynamical viscoelastic response experiments. Let us add that for solids, there is a similar quantity that can be obtained, out of the static elastic modulus $\mu$. Second, the QNMs in the $T_{xy}T_{xy}$ correlator (poles in the complex modulus) are also naturally compared to the inverse Maxwell relaxation time $1/\tau_{M}$. In generic materials, one can expect i) that there are several QNMs and ii) that their locations do not exactly coincide with $-i/\tau_{M}$, not even the lightest one (the one with smallest imaginary part). In the Maxwell model (9) the location of the QNM coincides with the inverse Maxwell time, $\tau_{M}^{-1}=G_{\infty}/\eta$, however this is due to the model simplicity. It is easy to break the coincidence by extending the model (as we show in Sec III.2 this can be done while having still a single QNM). In the holographic models below we will also find that the relaxation rate from the lightest QNM differs from $\tau_{M}^{-1}$, and this has an interesting impact in the viscoelastic response.

Once $G(\omega)$ is known we can also compute the real-time response of a system under an external and time dependent source. Viscoelastic materials exhibit several interesting features in the real-time response to a time dependent deformation. Among these, we shall mainly focus on :

• •

Stress Relaxation.

It is the decrease in stress in response to a constant strain generated in the structure, which can be seen as a delay of the material in reaching the steady-state. Its time dependence can typically be power-law or exponential, and in the latter case one can identify a relaxation time $\tau$. In the phenomenological Maxwell model the relaxation time is simply the ratio between the elastic modulus and the viscosity $\tau_{M}=\eta/G_{\infty}$ where the latter is the value of $G^{\prime}(\omega)$ at infinite frequency.

• •

Stress Overshoot.

It refers to the appearance of a transient excess stress before reaching the (non-zero) steady state. It is typically seen in start-up experiments, (e.g. when the source is turned on at a given time) and has been reported experimentally in a range of materials, from amorphous polymers and metallic glasses to foams and gels (see Falk and Langer (2011) for a review). More precisely the stress overshoot is characterized by the following temporal sequence: the stress $\sigma$ first increases linearly with time, until it reaches a maximum value and then decreases towards it steady-state constant value. Despite several explanatory attempts a complete understanding is still missing, see Divoux et al. (2011).

In the remainder of this paper we shall study the viscoelastic properties of the massive gravity models of Baggioli and Pujolas (2015); Alberte et al. (2017) by is computing the complex modulus and the response under time dependent sources. This opens the possibility of studying viscoelastic effects in strongly coupled materials or systems at quantum criticality. In the following we will show that the holographic setup at hand exhibits standard viscoelastic properties, starting by a clear stress relaxation. As expected, the relaxation time is controlled at finite temperature by the least damped QNM, specifically by the inverse of its imaginary part,

[TABLE]

Moreover, in general this relaxation time does not coincide with the Maxwell relaxation time $\tau_{M}=\eta/G_{\infty}$ (even though in some cases they are within the same order of magnitude, and even the temperature dependence is similar). Additionally, we find that the stress overshoot phenomenon appears when the relaxation time is larger than the Maxwell time, $\tau_{rel}>\tau_{M}$. Finally, at $T=0$ the stress relaxation becomes power-law and governed by the formation of a branchcut on the imaginary axis. The power is simply dictated by the conformal dimension of the stress operator on the near-horizon IR zero temperature geometry which is derived in details in Appendix B.

This paper is organized as follows: in section II we will set the stage presenting the holographic models we will consider, section III contains the main results of the manuscript, finally section IV summarize our findings and propose some future directions.

II Holographic setup

We consider generic viscoelastic holographic massive gravity models Baggioli and Pujolas (2015); Alberte et al. (2016a):

[TABLE]

which are inspired from the original model Andrade and Withers (2014). We define $\mathcal{I}^{IJ}=g^{\mu\nu}\,\partial_{\mu}\phi^{I}\partial_{\nu}\phi^{J}$ and $X=Tr[\mathcal{I}],Z=det[\mathcal{I]}$. The dual field theory has been shown to possess viscoelastic features and was studied in detail in Baggioli and Goykhman (2015); Baggioli and Brattan (2017); Baggioli and Goykhman (2016); Alberte et al. (2016b, 2018a, 2017). We consider 4D asymptotically AdS black hole geometries in Eddington-Filkenstein (EF) coordinates:

[TABLE]

where $u\in[0,u_{h}]$ is the radial holographic direction spanning from the boundary to the horizon, defined through $f(u_{h})=0$. The $\phi^{I}$ are the Stückelberg scalars which admit a radially constant profile $\phi^{I}=x^{I}$ with $I=x,y$ and the emblackening factor takes the simple form:

[TABLE]

The corresponding temperature of the dual QFT reads:

[TABLE]

while the entropy density is simply $s=2\pi/u_{h}^{2}$. The heat capacity can be simply obtained as $c_{v}=Tds/dT$ and was studied in Baggioli and Brattan (2017); Baggioli and Trachenko (2018a).

The deformation and flow properties of the model are encoded in the $T_{xy}$ operator which is dual, by the holographic dictionary, to the bulk deformation $\delta g_{xy}\equiv h_{xy}$. At zero momentum $k=0$ the linearized equation for $h_{xy}$ decouples and in EF coordinates reads:

[TABLE]

where $V_{X}\equiv\partial_{X}V(X,Z)$ and primes denote radial derivatives. The UV asymptotic behaviour of the $h_{xy}$ field is:

[TABLE]

The AdS/CFT dictionary allows us to express the Green’s function of the stress tensor as

[TABLE]

In order to obtain the retarded correlator, one needs to impose ingoing boundary conditions at the horizon $u=u_{h}$, which in EF coordinates corresponds to regularity of $h_{xy}$. We also compute the spectrum of quasi-normal modes by solving the boundary value problem resulting from setting setting the source to zero. They correspond to the poles of the Green’s function of the dual theory.

In the following we will consider potentials of the form $V(X,Z)=X^{n}$, $n>5/2$ and $V(X,Z)=Z^{m}$, $m>5/4$ which realize the spontaneous breaking of translational invariance, and exhibit both elastic and viscous response Alberte et al. (2017). For later use, we recall the expressions for the energy density in these models, this being

[TABLE]

We will refer to the previous class of models as Alberte et al. (2016a, b):

• •

Viscoelastic solid model: $V(X,Z)=X^{n}$, $n>5/2$

• •

Viscoelastic fluid model: $V(X,Z)=Z^{m}$, $m>5/4$

We take as an operational distinction between a fluid and a solid the presence of a finite static shear modulus $\mu\neq 0$. In other words, we will call solids all those configurations which posses a finite response to a zero frequency mechanical deformation. The adjective ”viscoelastic” emphasizes that in both models there is a competition between $G^{\prime}(\omega)$ and $G^{\prime\prime}(\omega)$ at finite $\omega$ due to the non-trivial frequency dependence of the response.

III Holographic Rheology

Rheology is the study of the deformation and flow of materials, both solids, liquids and especially viscoelastic ones. It combines the analysis of the viscous and elastic properties of a given material and it is based, at least at linear level, on the relation between the stress $\sigma(\omega)$ in a material and its consequent mechanical deformation, i.e. the strain $\epsilon(\omega)$. In a more formal way it is based on the relation (5) involving the complex dynamic modulus $G(\omega)$. This can be related to the correlator of particular spatial components of the stress tensor as in (6), which facilitates our holographic treatment. In our language the dynamical modulus is simply encoded in the Green function of the shear operator $T_{xy}$ which can be extracted numerically from (64) using the standard holographic techniques.

Once $G(\omega)$ is obtained we can compute the real-time response of a system under an external and time dependent source $\sigma(t)$. To this end, we will Fourier transform back to the time domain and write

[TABLE]

where $\tilde{\epsilon}(\omega)$ is simply the Fourier transform of the external strain source $\epsilon(t)$. For simplicity we will consider the following external sources:

• •

Logistic function.

[TABLE]

It represents a smooth deformation of a step function and its Fourier transform is:

[TABLE]

• •

Ramp function.

[TABLE]

whose Fourier transform is:

[TABLE]

III.1 Fluids

We start by considering the viscoelastic fluid model $V(X,Z)=Z^{n}$. For concreteness we restrict ourselves to the potential

[TABLE]

which corresponds to consider a dual QFT in a fluid phase Alberte et al. (2016a) which exhibits spontaneous breaking of translational invariance. Other potentials of the form $V(X,Z)=\mathcal{V}(Z)=Z^{n}$, with $n>5/4$ give similar results444Nevertheless notice that this setup can provide also a mechanism for the explicit breaking of translations (for example choosing $\mathcal{V}(Z)=Z$) and therefore a finite DC conductivity still having zero elastic modulus $\mu=0$ Alberte et al. (2016a).. In particular, since $V_{X}=0$, it is easy to prove Alberte et al. (2016b) that for this model:

[TABLE]

which is just a consequence of the fact that the shear component of the graviton remains massless. In other words, this model has no static elastic response but purely a viscous one, which justifies the term “fluid”.

We start by discussing the finite frequency response of the complex dynamical modulus $G(\omega)$. The results are shown in fig.2. Note that by time reflection symmetry, $G^{\prime}(-\omega)=G^{\prime}(\omega)$, $G^{\prime\prime}(-\omega)=-G^{\prime\prime}(\omega)$, so we only show our results for $\omega\geq 0$. The same argument applies to the structure of QNMs which appear in pairs symmetric under ${\rm Im}\omega_{QNM}\to-{\rm Im}\omega_{QNM}$. As already described, at low frequencies we find $\mu=G^{\prime}(0)=0$ while the viscosity saturates the KSS bound $\eta/s=1/4\pi$ Kovtun et al. (2005). The imaginary part $G^{\prime\prime}(\omega)$ goes to zero at large frequency, so no dissipation is present in this regime. Physically, this means that if the system is perturbed at a time scale much shorter than its typical relaxation time, it is not able to dissipate. It is interesting to note that at large frequencies, the real part of $G(\omega)$ takes a constant value which we denote $G_{\infty}\equiv G^{\prime}(\omega=\infty)$. As we show in appendix A, this can be computed analytically in terms of the energy density as:

[TABLE]

independent of the specific form of the potential.

In the intermediate frequency range we observe a non-trivial frequency dependence and, in particular, a clear resonance in $G^{\prime}(\omega)$. Not surprisingly such a resonance is also visible in the QNMs spectrum of the system which is shown for different temperatures in fig.3. Note in particular that the position of the isolated complex pole matches quite well the position of the peak in the real part of the complex modulus.

From the positions of the QNMs, and in particular from the imaginary part of the least damped one, we can read off the relaxation time which is shown in the right panel of fig.3. Interestingly, at a certain value $m/T\approx 2.5$ there is a crossover between the two dominating poles. Upon further lowering the temperature, a series of purely imaginary poles accumulate on the imaginary axes and form a branch cut at exactly zero temperature as usual in models with an extremal horizon at $T=0$, see Faulkner et al. (2011).

Next, we consider the real time mechanical response of the system upon applying an external strain in the form of a ramp and step sources. We display our results in fig. 4. The response to the step source is simple and confirms the fluid nature of the material. In fact, we can see from the right panel of fig.5 that the system reacts only to a gradient of the strain $\dot{\epsilon}(t)$ and only develops a peak close to $t=0$. At late times where the strain is constant, no response is present, as expected for a fluid.

The response to the ramp function is more interesting. At small $m/T$, $\sigma(t)$ follows the expected form of early growth and (exponential) decay similar to the one observed in amorphous solids and disordered systems; see Zaccone et al. (2014) for an example. As we increase $m/T$ the response departs from such a form and develops a clear peak before reaching the late time Newtonian plateau. Such a feature is a distinct signature of the presence of a stress overshoot, and is produced by the interplay between the viscous behaviour and some elastic feature present at finite frequency. Indeed, as shown in fig.2, despite that the shear modulus is zero, i.e. $G^{\prime}(0)=0$, the real part of the complex modulus exhibits a strong peak at intermediate frequency which is governed by the first resonance in the QNM spectrum and whose strength grows with the parameter $m/T$. Correspondingly, we see that the overshoot peak grows with $m/T$.

In both the ramp and the step setups the stress relaxes at late time towards a steady-state value (which turns out to be zero for the step source because of the absence of a finite static shear modulus $\mu$). The stress relaxation, which is exponential at finite temperature, is governed by the imaginary part of the least damped QNM. Indeed, we find good agreement between the real time data at late times and the QNMs, see fig.5, even at considerably small $T$, where the lightest pole is on the imaginary axis. In the $T=0$ limit one expects that the tower of poles in the imaginary axis merge into a branch cut reaching down to $\omega=0$, and this turns the relaxation in real time to be power-law rather than exponential. While this is hard to see explicitly in the real time numerics, we find evidence of the formation of a branch cut on the imaginary axis555However, we do see the power-law scaling in our toy model to be discussed in Sec. III.2.. On general grounds, this power law is controlled by the anomalous dimension that the stress tensor operator can get in the IR, which corresponds to the graviton mass in the $AdS_{2}$ near-horizon geometry. We compute this in appendix B.

Finally we compare the relaxation time at finite temperature which is given by the imaginary part of the least damped QNM:

[TABLE]

to the Maxwell relaxation time:

[TABLE]

Let us emphasize that this parameter is certainly present in the system simply because, just like in the Maxwell model, there is a finite $G_{\infty}\equiv\lim_{\omega\to\infty}G(\omega)$. In contrast with the simple Maxwell model, though, this relaxation time might not coincide with any QNM relaxation time. The comparison between these relaxation time scales can be seen in fig. 3. Quite remarkably, $\tau_{M}$ is rather close to the actual relaxation time. At high $T$, the difference is only a factor $\sim 3$ and the temperature dependence is the correct one. This is so until the pure imaginary QNM becomes dominant at $m/T\gtrsim 2.5$. From then on, the temperature dependence of $\tau_{M}$ and of the actual relaxation time $\tau_{QNM}$ differ. At $m/T\simeq 10$, the Maxwell time $\tau_{M}$ turns from over-estimating to under-estimating the actual relaxation time.666The mismatch between the holographic relaxation time and the Maxwell relaxation time have been already observed in a slightly different context in Baggioli and Trachenko (2018a, b). In hindsight, the failure of $\tau_{M}$ is built in in the fact that at $T=0$ our model has an AdS2 critical point, meaning that the imaginary pole must become arbitrarily light or equivalently the relaxation must be power law, not exponential in time.

The discrepancy between $\tau_{M}$ and $\tau_{QNM}=\tau_{rel}$ also seem to have an impact on the overshoot phenomenon seen in fig. 4. As we see from this figure and fig. 3, overshooting is present when the condition

[TABLE]

is met. This condition will be confirmed in the two more cases studied in the next subsections. Heuristically, it seems reasonable that if a material possesses dynamical modes (resonances) that are light (decay slowly) then effects coming from these resonances are visible. The above condition then represents the criterion that decides whether a resonance is light or not, and it is very natural that this involves the viscoelastic moduli $\eta$, $G_{\infty}$.

III.2 A simple toy model

In this section we construct a simple toy model for the Green function of our system. Despite its simplicity, the model is able to reproduce several features present in the more complex holographic setup. In order to mimic the $T=0$ solution we take a Green’s function which contains a branch cut on the imaginary axis, as in the AdS2 dual geometry, and a series of isolated poles. More precisely we consider:

[TABLE]

where we have set the residues of all the poles equal. By choosing $m_{k}$ and $\Gamma_{k}$ appropriately, one can arrange for an array of poles, such as forming an inverse parabola-like shape like in the holographic models. However, it will be sufficient to consider only the least damped mode of the tower so we shall not enter into specifying this further.

We will specialize to the fluid case by choosing the constant $c$ that controls the shear elastic modulus:

[TABLE]

such that $\mu=0$. The power $\omega^{\alpha}$ produces a branch cut on the imaginary axis starting from the origin. An example of the pole structure is shown in fig.6. To illustrate the capabilities of this toy model, we consider the response to a step function at strictly zero temperature. We display our results for a specific choice of parameters in fig.7. We see that our model captures well the essential aspects of the full holographic model shown in fig.4, showing that the response is dominated by the first QNM.

Moreover, the presence of the branch cut $\sim\omega^{\alpha}$ accounts for the expected power law decay at $T=0$. In fact, considering only the branch cut term, which dominates at low frequencies, and the fact that the step function is constant at late times, we can derive the scaling:

[TABLE]

at late times. This is indeed consistent with the numerical results shown in fig.7. Moreover we see that at $T=0$ the late time dynamics is controlled by the branch cut formed on the imaginary axis, and not by the gapped QNM, which would yield exponential decay.

The toy model is also useful to describe the stress overshoot phenomenon present in the start-up experiment. Indeed, truncating to a single resonance (with non-zero real part),

[TABLE]

we can fix the constant $c$ so as to mimic a fluid model where $\mu=G(0)=0$. This fixes $c=b/m^{2}$. Note that the parameter $\Gamma$ in this model is the imaginary part of the QNM location.

In this model it is very easy to obtain the elastic modulus and viscosity,

[TABLE]

meaning that the Maxwell relaxation time in this case is:

[TABLE]

This now can be compared to the relaxation time encoded in the QNM,

[TABLE]

Therefore in this simple model, the ratio $\tau_{M}/\tau_{QNM}=(\Gamma/m)^{2}$, which can be chosen to be large or small.

Now we can simply compute the response of the system to a ramp strain. The results are shown in fig.8. Clearly, the overshoot is present for $\Gamma\ll m$, that is

[TABLE]

This is the same pattern that is seen in the fluid models of Sec. III.1, and as a matter of fact it will also be found in the solid models below (concerning the overshoot phenonenon): materials that develop a light QNM exhibit the overshoot phenomenon. Here by ‘light’ we mean that the damping scale of the QNM is smaller than the relaxation scale encoded in the material moduli (elasticity, shear).

A rather simple picture seems to emerge, then: in presence of such a light QNM, the approach to steady state is delayed and moreover it can go together with oscillatory features in the response expected from the fact that one is exciting a really dynamical mode in the system. This picture seems conceptually quite different from other models or interpretations that ascribe the overshoot phenomenon to a superposition of elastic and viscous behaviours. It is clear in our models that the feature behind overshoot is a light QNM. It is of course interesting to understand under what conditions does a QNM go below the Maxwell relaxation time. In our solutions, this seems granted because by decreasing temperature the poles on the imaginary axis must get eventually arbitrarily light. In turn, this is a consequence of the fact that our holographic solutions always display criticality at $T=0$ (an AdS2 extremal horizon).

III.3 Viscoelastic solids

In this section we consider a different class of model, which we denote as viscoelastic solid because of the presence of a finite static shear modulus $\mu\neq 0$. More specifically, we study potentials of the form $V(X)=X^{n}$, which realize the spontaneous symmetry breaking of translational invariance and exhibits viscoelastic features along with the presence of damped phonons Alberte et al. (2018a, 2017). As previously analyzed in Alberte et al. (2016a, b, 2017), the shear elastic modulus and the shear viscosity coefficient of a class of holographic models of the form (13), exhibit a viscoelastic behaviour and a smooth glassy transition. More precisely, the shear modulus $\mu$ goes smoothly to zero increasing the dimensionless parameter $T/m$ and it reaches its maximum value at zero temperature. The viscosity behaves exactly in the opposite way. For a schematic representation see fig.9. The holographic systems behave as a perfect dissipationless solid at $T=0$, as a strongly coupled viscous fluid at $T\rightarrow\infty$ and as viscoelastic materials in the intermediate range. As a consequence, dialling the dimensionless parameter $m/T\in[0,\infty]$ we can study our background in the viscoelastic regime. For brevity, concentrate on the case $n=3$. We have checked that higher power-law potentials $V(X)=X^{n}$ behave in a qualitatively similar fashion.

We start by studying the dynamical modulus $G$ in the frequency domain. We show our results for the frequency dependent Green’s function in fig.10.

First, note that our numerics reproduce the generic expansion at low frequencies of the form (8) so that one can readily identify the shear viscosity and elasticity moduli. Note that the main difference with respect to the previous fluid case is that now

[TABLE]

signalling that the model now behaves as a solid in the conventional sense (from the static response). It can be appreciated in fig. 10 that this modulus drops quite quickly to zero at high temperatures, corresponding to a melting-like transition (see Alberte et al. (2017, 2018a)). In fact the whole functional dependence of $G^{\prime}$ and $G^{\prime\prime}$ in the high $T$ limit resemble the one in the fluid limit, fig. 2.

Secondly, at large frequencies $\omega\rightarrow\infty$ the imaginary and dissipative part of the complex moduli goes to zero while the real and elastic part asymptotes a constant value, $G_{\infty}$, again given by the energy density as in (28). This immediately means that there are (at least) two notions now of elastic modulus, and therefore of Maxwell relaxation times:

[TABLE]

On general grounds, it is natural to expect that the relevant one at high temperatures is $\tau^{(\infty)}_{M}$, as in Sect. III.1, because in this limit material melts, meaning that $\mu$ must vanish and the fluid behaviour must be recovered. Conversely, at low temperatures the a static elastic modulus $\mu$ is expected to be sizeable. Then, the late time response is encoded directly in the low frequency expansion of the Green’s function (8), suggesting that the relevant relaxation time is should be given by $\tau^{(0)}_{M}$ (at low $T$). In any case, a glance at the left panel of fig. 10 reveals that in our models at low temperatures $\mu$ and $G_{\infty}$ are comparable in magnitude so there isn’t a big difference in the two Maxwell times in this regime.

Next, we observe that the position of the dominant finite frequency peak in $G^{\prime}$ is approximately given by the real part of the leading gapped QNM, see fig.11. Like in the fluid model, there are two ‘light’ QNMs: a pure imaginary one and an off-axis mode, and the at sufficiently low $T$ the imaginary mode is granted to have lowest imaginary part because $T=0$ the accumulation of poles in the imaginary axis signals the formation of a branch cut. In fact, one can estimate the temperature at which the cross-over between the dominance of the gapped QNM and the branch cut by comparing the imaginary parts of the first propagating and damped QNMs, see right panel of fig. 11.

It is relevant to compare the relaxation times derived from the lowest of these QNM modes, $\tau_{QNM}\,=\,-\,1/Im\left(\omega_{QNM}^{*}\right)$, to the above Maxwell relaxation times $\tau_{M}^{(0)}$, $\tau_{M}^{(\infty)}$. As we can see in the right panel of fig. 11, both $\tau_{M}^{(0)}$ and $\tau_{M}^{(\infty)}$ go below the $\tau_{QNM}$, towards low temperature. In fact, they do so at a surprisingly coincident temperature, roughly $m/T\simeq 4$. Since this crossover takes place at smaller $m/T$ as compared to the fluid case of Section III.1, the QNM responsible for this crossover is the complex one in the $V(X)$ model. At higher $m/T$, there is a second crossover where the imaginary QNM becomes longer lived than the complex one. However, the first crossover at $m/T\simeq 4$ is also seen to correlate well with the appearance of overshoot in the step-function external source, as is especially clear in the right panel of fig. 12.

Another significant difference with respect to the fluid case is that the temperature dependence of the relaxation times from the QNMs $\tau_{QNM}$ or from the Maxwell formulas (38) differ substantially: now both $\tau_{M}^{(0,\infty)}$ decrease much faster towards $T=0$ than $\tau_{QNM}$. It is natural to associate this to the fact that the solid models that we are considering now are known to violate Alberte et al. (2016b); Hartnoll et al. (2016b) the KSS bound Kovtun et al. (2005). In these models $\eta/s$ decreases ‘abnormally’ at low $T$ (while $s$ remains nonzero) so perhaps this is the only reason why the Maxwell time does too. More specifically, since viscosity and entropy density have the same units, it is also possible to construct a relaxation time out of the ratio

[TABLE]

with $G$ denoting any of $G_{\infty}$ or $G(0)$ and the $4\pi$ factor is introduced in order to match correctly the fluid limit (high temperature). As shown in the right panel of fig. 11 the temperature behaviour of $\tau_{s}$ tracks quite well the relaxation times from the lowest QNM. We leave a further study of this for later investigation.

We display our results for the real time response in fig. 12. Let us discuss first the ramp source. In this case, at small $m/T$ the system reacts like a viscous fluids where the strain $\sigma(t)\sim\dot{\epsilon}(t)=const.$; nevertheless before reaching the constant value there is a time delay of the system with a characteristic timescale $\tau_{d}$. Once $m/T$ is increased, a finite shear modulus $\mu\neq 0$ appears and therefore the response continues to grow linearly at late times. Notice that the slope of the late time response indeed grows with $m/T$ as the shear modulus does. The dissipative response provides only a small deviation from the linear growth which is visible in a transient regime. In realistic materials once the strain becomes large, i.e. at late time, non-linear effects become important. As a consequence what usually happens is that the rigid bonds in the material become weaker and the stress stops to grow and it reaches the so-called Newtonian plateau which for example we see at $m/T=0$. Clearly these effects cannot be captured by our analysis which is limited to the linear response approximation. In some cases non-linear effects produce overshooting Zaccone et al. (2014). It would be interesting to see if this effect can be captured by a holographic calculation. We hope to come back to this interesting point in the future.

We now turn to the step function. In that case at $m/T=0$ the response is purely viscous: it develops a peak when the gradient of the strain is large (at the step position) and then it dies off as expected since at late time $\dot{\epsilon}(t)=0$. On the contrary, switching on $m/T$ the shear modulus is non-vanishing and so that the response at late time asymptotes a constant proportional to the shear modulus. Before reaching such a plateau a peak appears in the response due to the competition between the viscous and elastic response. The peak is sharper and higher the stronger is the elastic response, which is governed by the resonance in the QNMs spectrum. The residue of such resonance controls the strength of the peak and its imaginary value the exponential decay that follows it towards the late time plateau. In fact, we have checked that the approach to equilibrium is governed by the first QNM, see fig 13. In the range of temperatures we have considered, we obtain a good fit considering the lightest QNM only.

IV Conclusions

We have studied the real time mechanical response of a class of solid and fluid holographic massive gravity models. For simplicity we have restricted ourselves to the shear sector and to consider the linear response approximation valid only for small external sources. Within this framework we have obtained the complex moduli $G(\omega)$ from the Green’s function of the stress tensor operator $T_{xy}$ at finite frequency. Our results support the classification of the models at hand into a solid and a fluid class accordingly to the value of the shear modulus $\mu\equiv G(0)$.

We have computed, within linear response, the real time response of the systems under an applied external time dependent strain. The results show interesting features which are typically observed in viscoelastic materials: stress overshoot and stress relaxation.

Our results have a number of implications:

They confirm that the holographic massive gravity models Baggioli and Pujolas (2015); Alberte et al. (2016a) as dual bulk descriptions of (strongly coupled) viscoelastic materials, either viscoelastic solids or fluids depending on the choice of potential. 2. 2.

Just like in the phenomenological Maxwell, model these models possess a well defined elastic modulus defined as $G_{\infty}=\lim_{\omega\to\infty}G(\omega)$ that is a key player in the viscoelastic phenomena and which is nonzero even fluids (which have $G(0)=0$). Moreover, in our models we find that this is basically given by the energy density, see Eq. (28). 3. 3.

At finite temperature, the relaxation time is parametrically similar to the Maxwell relaxation time $\tau_{M}=\eta/G_{\infty}$, but in close inspection there are important differences between this and the actual relaxation time. At high temperatures $\tau_{M}$ is only a factor $\sim 3$ different (larger) and moreover it captures correctly the temperature dependence. However, at lower temperatures, the $T$ behaviour starts to differ and at small enough $T$, $\tau_{M}$ severely underestimates $\tau_{rel}$. In our models, this must happen because dynamics at $T=0$ is controlled by an infrared critical point, and the relaxation in time must be power-law. It is tempting to speculate whether a similar behaviour holds more generally – that whenever the mechanical deformation of a material is approaching a critical regime, then the Maxwell relaxation time $\tau_{M}$ must necessarily underestimate the actual relaxation time. 4. 4.

We find that the appearance of a stress overshoot peak preceding relaxation correlates well with the underestimation condition $\tau_{M}<\tau_{rel}$. As we have argued, this should be granted to occur in materials that are close to criticality. It would be interesting to understand whether this can occur in non-critical materials, via other mechanisms Divoux et al. (2011); Laurati et al. (2017); Zaccone et al. (2014), or whether they are two sides of the same reality.

One interesting and natural future direction is extending our computations to the non-linear regime. The latter could produce more realistic features at late time and shed more light on the nature of the stress overshoot in strongly coupled viscoelastic materials. One possibility is to follow the methods used for the non-linear eletric conductivity in Withers (2016); Baggioli and Pujolas (2016) or take inspiration from the models for non-linear elasticity developed in Alberte et al. (2018b); Baggioli et al. (soon).

It is interesting to compare our power law relaxation with the existing literature. A viscoelastic model that predicts a power-law scaling via a fractional relaxation process has been introduced in Pritchard and Terentjev (2017) and analyzed further in Milkus and Zaccone (2017). There is some evidence that there are viscoelastic materials, like polymeric media, biological tissues and cells (see Metzler and Nonnenmacher (2003); Shen et al. (2013) for some examples), where the power law scaling is at work. See La Nave and Phillips (2017); Limtragool et al. (2016); Vlahos et al. (2008) for an incomplete list of references that use AdS/CFT techniques to model phenomena where power-law scaling is important.

As a final and promising direction we have to mention the recent interest around the study of phonons and glassy/viscoelastic features in the context of quantum critical scale invariant systems. In particular two important questions arised: what is the behaviour of the phonons at quantum criticality Ishii et al. (2019)? How do viscoelastic properties affect the onset of high-Tc superconductivity Setty (2019)? Preliminary results suggest a strong correlation between phonons at criticality, viscoelasticity and glass-enhanced high-Tc superconductors. The methods presented here seem useful to address some of these questions. We hope to revisit these issues in the future.

Acknowledgements

We thank Kostya Trachenko and especially Alessio Zaccone for useful discussions and comments about this work and the topics considered. We would like to thank Alessio Zaccone for explaining us all the secrets of viscoelastic materials and the related open problems at any time of the day and the night. The work of TA is supported by the ERC Advanced Grant GravBHs-692951. MB acknowledges the support of the Spanish MINECO’s “Centro de Excelencia Severo Ochoa” Programme under grant SEV-2012-0249. OP acknowledges support by the Spanish Ministry MEC under grant FPA2014-55613-P, FPA2017-88915-P the Severo Ochoa excellence program of MINECO (grant SO-2012-0234, SEV-2016-0588), as well as by the Generalitat de Catalunya under grants 2014-SGR-1450, 2017SGR1069. TA gratefully acknowledges the hospitality of the Lorentz Center during the workshop “Bringing Holography to the Lab: Explaining Strange Metals with Virtual Black Holes”, held in Leiden, 2019, where part of this work was carried out. MB would like to thank Marianna Siouti for the unconditional support. MB would like to thank IFAE and UPC Barcelona for the warm hospitality during the completion of this work. We would also like to acknowledge IFAE for providing the comfortable atmosphere where this project started.

Appendix A The instantaneous elastic modulus and the energy density

In this appendix we obtain the large frequency result for the two point function of the stress tensor. Let us start with the generic equation for the shear mode:

[TABLE]

and let us write down the blackening factor in the general form in terms of the energy density $\epsilon$ (setting the radius of the horizon to $u_{h}=1$) :

[TABLE]

which is always the case for the models we consider. The specific value of $n$ is actually not relevant for the following derivation but importantly $n\neq 3$. If we for example choose a potential of the type $V(X,Z)=X^{p_{1}}Z^{p_{2}}$ then we have $n=4p_{2}+2p_{1}$.

Imposing regularity at the horizon, we find that the series solution near $u=1$ has the form

[TABLE]

On the other hand, the general near boundary expansion has the form

[TABLE]

where

[TABLE]

where we have made the $\omega$ dependence explicit, and $c_{0}$, $c_{1}$ some numerical $\omega$-independent constants, and

[TABLE]

From (41), we observe that an appropriate approximate solution in the limit $\omega\to\infty$ has the form

[TABLE]

Here, we have set $h^{H}=1$ which is always valid by linearity of the equation. In order for (45) to hold approximately everywhere, we choose

[TABLE]

Note that there could be subleading contributions in powers of $u$ in (42), however, these would be suppressed since $0<u<1$. We have verified numerically that indeed, for large $\omega$ our solutions have the form (45). We can read off the Green’s function from this approximate solution, which yields,

[TABLE]

as promised.

We suspect this is a universal relation for holographic systems which we plan to investigate in the nearest future.

Appendix B IR conformal dimension of the shear operator

In this section we compute the conformal dimension of the shear operator $T_{xy}$ at the IR fixed point or in other words at the extremal near horizon geometry AdS${}_{2}\times R^{2}$.

As shown in the main text the UV conformal dimension of the stress tensor does not get modified by the graviton mass. The radial dependent graviton mass indeed can be defined as:

[TABLE]

which indeed is always zero at the UV boundary $u=0$. As a consequence the conformal dimension of the stress tensor is the canonical one, in this case \big{[}T_{xy}\big{]}=3.

At the zero temperature IR fixed point the situation is different. Let us consider the extremal near-horizon geometry:

[TABLE]

where $u_{0}$ is the position of the extremal horizon and $L$ the radius of the AdS4 manifold. It is simple to realize that the above geometry is an AdS${}_{2}\times R^{2}$ spacetime where the radius of the AdS2 geometry reads:

[TABLE]

For the models we consider we find:

[TABLE]

Now we can solve the equation of motion close to the horizon obtaining:

[TABLE]

where:

[TABLE]

and:

[TABLE]

Putting all the results together and assuming a potential of the form $V(X,Z)=X^{a}Z^{(b-a)/2}$ we obtain777The consistency of the model (no instabilities, positive elastic moduli, etc) imposes that the parameters satify $a>0$, $b>3/2$ Alberte et al. (2018b); Baggioli et al. (soon) :

[TABLE]

The conformal dimension $\Delta$ of the $T_{xy}$ operator is then simply:

[TABLE]

In the case of the $X$ dependent potential $a=b$ we obtain $\Delta=2$ which is in agreement with the result for the linear potential presented in Hartnoll et al. (2016b).

The just derived conformal dimension is governing two different phenomena:

• •

The late time power-law relaxation of the stress $\sigma(t)$ which appears at $T=0$. Hereby we sketch the simple argument behind it. The Green function for an operator with dimension $\Delta$ does scale like:

[TABLE]

where $D$ is the number of spatial dimensions and $z$ the Lifshitz scaling parameter. In our case $D=2$ and $z$ is formally infinite because of the AdS2 factor in the near-horizon extremal geometry. In few words the Green function of the stress tensor operator, which controls the linear stress-strain response, scales like $\sim\omega^{2\Delta-1}$ with the conformal dimension obtained above in (56). By dimensional analysis we can then simply deduce that at late time the stress will relax following the power law scaling:

[TABLE]

which is always faster than the linear scaling $t^{-1}$ which appears at zero graviton mass $V_{X}=0$ or equivalently $a=0$. Notice that this implies:

[TABLE]

at zero temperature.

• •

The power-law fall-off of the $\eta/s$ ratio at zero temperature:

[TABLE]

Notice that indeed for $V_{X}=0$ ($a=0$) this gives that the $\eta/s$ ratio is constant in temperature as it must be. Notice also that (61) extends the results of Hartnoll et al. (2016b), which found that $\eta/s$ approaches [math] no faster than $T^{2}$. Indeed, a power larger than 2 can be obtained for $a>b$, which is compatible with consistency of the model Alberte et al. (2018b); Baggioli et al. (soon). Hence, the fall-off of the viscosity/entropy ratio can be arbitrarily fast without appealing to hyperscaling violation or Lifshitz scaling as in Ling et al. (2017).

Appendix C Details about the numerical methods

In this appendix we provide more details about the numerical techniques used in the manuscript and the derivation of the main results.

The Green function

The stress tensor correlator $\langle T_{xy}T_{xy}\rangle$ is extracted by considering the bulk equation for the shear perturbation

[TABLE]

written in EF coordinates. The UV asymptotic behaviour of the $h_{xy}$ field close to the boundary $u=0$ reads:

[TABLE]

and finally the Green function can be read using the standard holographic dictionary as:

[TABLE]

In order to derive this function numerically we use a double shooting matching procedure described in the following steps888For more details about this method see the lectures notes in Baggioli (2019).

We fix a specific form of the potential $V(X,Z)$ in eq.(62). 2. 2.

We solve perturbatively equation (62) close to the horizon $u_{h}=1$ using an ansatz:

[TABLE]

where we impose regularity and we fix $N=5$ to have enough accuracy. The free parameters of the IR expansion are $\{h_{0},\omega,m^{2}\}$. 3. 3.

We construct numerically the solution to equation (62) by shooting from the IR with the boundary conditions obtained in the previous point to an intermediate radial position $0<u_{m}<1$ (we concretely choose $u_{m}=1/2$). 4. 4.

We repeat exactly the same procedure from the UV. In this case the free parameters are $\{h_{xy(l)},h_{xy(s)},\omega,m^{2}\}$. 5. 5.

For a generic choice of parameters the two solutions will be not continuous at the intermediate point $u_{m}=1/2$ as shown in fig.14.

The first step to continue is noticing that one between $h_{0}$ and $h_{xy(l)}$ is redundant and it can be set to $1$ by rescaling the bulk field $h_{xy}(u)$. After doing that, we have simply to impose the continuity of the function $h_{xy}(u)$ and its first derivative at the intermediate position $u_{m}=1/2$. This can be easily done with a FindRoot routine in Mathematica. These two requirement will fix for us $h_{xy(l)},h_{xy(s)}$ and therefore our Green’s function in terms of $\omega$ and $m^{2}$.

Quasinormal modes

The QNMs of the system correspond to the poles of the Green’s function. From (64), we see that this corresponds to setting $h_{xy\,(l)}(\omega)=0$. Imposing ingoing boundary condition at the horizon, (62) turns into an eigenvalue problem for the frequency $\omega$ which is satisfied by a discrete number of complex values for a given set of external parameters. We solve this by discretizing (62) on a Chebyshev grid and solving the resulting linear algebra problem numerically. If one wishes, this problem can also be solved by shooting as explained in the previous section.

Late time relaxation and fits

As it is well-known, the lowest QNMs govern the approach to equilibrium of black holes. In the main text, we have checked that this is the case by employing the late time expansion

[TABLE]

where $N_{QNM}$ is the number of lowest QNMs included in the fit, with frequencies $\omega_{i}$, the values of which we have extracted as explained in the previous section. The quantity $\sigma(t=\infty)$ is the constant late time value of $\sigma$. The complex coefficients $a_{i}$ are obtained by computing the best fit to the numerical value of the left hand side of (66). We have checked the robustness of our results against changing the length of the time interval in which the fits are performed. In this simple setting, over-fitting is signalled by large values of some of the coefficients $a_{i}$. We avoid this by keeping a value of $N_{QNM}$ not higher than 3.