Nonlinear Oscillatory Shear Tests in Viscoelastic Holography

We provide the first characterization of the nonlinear and time dependent rheologic response of viscoelastic bottom-up holographic models. More precisely, we perform oscillatory shear tests in holographic massive gravity theories with finite elastic response, focusing on the large amplitude oscillatory shear (LAOS) regime. The characterization of these systems is done using several techniques: (I) the Lissajous figures, (II) the Fourier analysis of the stress signal, (III) the Pipkin diagram and (IV) the dependence of the storage and loss moduli on the amplitude of the applied strain. We find substantial evidence for a strong strain stiffening mechanism, typical of hyper-elastic materials such as rubbers and complex polymers. This indicates that the holographic models considered are not a good description for rigid metals, where strain stiffening is not commonly observed. Additionally, a crossover between a viscoelastic liquid regime at small graviton mass (compared to the temperature scale), and a viscoelastic solid regime at large values is observed. Finally, we discuss the relevance of our results for soft matter and for the understanding of the widely used homogeneous holographic models with broken translations.

We provide the first characterization of the nonlinear and time dependent rheologic response of viscoelastic bottom-up holographic models. More precisely, we perform oscillatory shear tests in holographic massive gravity theories with finite elastic response, focusing on the large amplitude oscillatory shear (LAOS) regime. The characterization of these systems is done using several techniques: (I) the Lissajous figures, (II) the Fourier analysis of the stress signal, (III) the Pipkin diagram and (IV) the dependence of the storage and loss moduli on the amplitude of the applied strain. We find substantial evidence for a strong strain stiffening mechanism, typical of hyper-elastic materials such as rubbers and complex polymers. This indicates that the holographic models considered are not a good description for rigid metals, where strain stiffening is not commonly observed. Additionally, a crossover between a viscoelastic liquid regime at small graviton mass (compared to the temperature scale), and a viscoelastic solid regime at large values is observed. Finally, we discuss the relevance of our results for soft matter and for the understanding of the widely used homogeneous holographic models with broken translations.

I. INTRODUCTION
In elastic solids, the mechanical stress is proportional to the applied external shear strain [1]. However, in hydrodynamic fluids the stress is proportional to the shear rate [2]. Of course, both of these cases are abstract idealizations, valid only under limiting conditions. In general, all the materials are viscoelastic -they present an interplay between elastic effects and dissipative viscous ones [3]; honey is the most common example. The idea that "everything flows if you wait long enough" lies behind the foundation of a new field of research known as rheology [4,5] -the study of deformation and flow of matter. Even though the first model goes back to Maxwell in 1867 [6], a large part of the theoretical description of viscoelastic materials is still based on phenomenological frameworks (Kelvin-Voigt, generalized Maxwell, Burgers) [7]. The fundamental difficulties are twofold: (I) it is conceptually hard to incorporate dissipation into the effective field theory description of solid materials [8,9] because of the unavoidable requirement of unitarity; (II) elasticity theory can be formulated in the language of standard effective field theory following a well defined action principle [10]. Hydrodynamics, on the contrary, is usually described by a set of conservation equations and constitutive relations [11] and it is not suitable for a description in terms of a local and hermitian action [12]. The problem becomes even more acute when the amplitude of the applied external strain is not small and the linear approximation is of no help anymore -the onset of nonlinear viscoelasticity [13][14][15][16]. For simplicity, in this letter, we focus on oscillatory shear tests in which the external shear strain takes a simple sinusoidal form where γ 0 is the strain amplitude and ω is the characteristic frequency. A convenient characterization of these rheology experiments is defined through the Deborah number De = ω λ and the Weissenberg number Wi ≡ λ γ 0 , where λ is the characteristic relaxation time of the material -in simple words, these two numbers determine how fast and how strong we are probing the viscoelastic system. Small Wi and large De corresponds to linear elasticity, in which the stress output is linearly dependent on the external input. Generally speaking, we can draw the so-called Pipkin diagram [16], picturing the phase space of the system in function of the values of these two numbers. In the largest region of this diagram, strain amplitudes are large and frequencies are neither high nor low; experiments probing that region are called LAOS tests [17,18] and are the subject of this letter.
In the LAOS regime, linear viscoelasticity is not applicable anymore; the response is fully nonlinear, the storage and loss moduli become non-trivial functions of the strain amplitude γ 0 . Very little is known in this regime, quoting Pipkin himself: "Here Be Dragons" [16]. From a totally different perspective, in the last ten years, Holography revealed to be a very useful tool for the development and understanding of Hydrodynamics [19][20][21][22][23][24].
The most famous examples are: (I) the formulation of a universal bound on the viscosity-to-entropy ratio [25] which is so far respected by all known fluids [26]; (II) the discovery of new transport coefficients in anomalous hydrodynamics [27,28] experimentally observed in Weyl semimetals [29]. A fundamental breakthrough in this direction is the observation that black holes (BHs) behave as dissipative hydrodynamic systems [30,31]. From this point of view, the application of an external strain source to the hydrodynamic system corresponds to the perturbation of the BH geometry by dynamical gravitational waves [19]. More recently, a series of works [32][33][34][35] explained how to endow black holes with a solid structure, providing them with a finite elastic response. In these new holographic theories, the BH response is no longer purely hydrodynamic but it becomes viscoelastic [36] in all aspects. Since then, a lot of effort has been devoted to the implementation, the classification and the characterization of these setups and similar ones [37][38][39][40][41][42][43][44][45][46][47][48]. In this letter, we provide the first characterization of the nonlinear and time dependent viscoelastic response of these holographic models, with particular emphasis on the LAOS regime. The relevance of our results is diverse and highly interdisciplinary: (I) to shed light on the challenge of LAOS and in particular the physics of complex fluids (yelding, shear thinning, stress overshoot, dynamical instabilities) [18]; (II) to reach a full characterization and understanding of the homogeneous holographic models with broken translations [32,[49][50][51] and their possible connections with glasses, complex fluids and amorphous systems [52]; (III) to study out of equilibrium processes in strongly coupled field theories and the possible universal evolution after dynamical quenches. Similar studies have been performed in [53], for a CFT driven by an oscillating composite scalar operator, and [54] where a gapped holographic system has been perturbed with a homogeneous gravitational periodic driving. Some qualitative features observed in [53] are totally consistent with our findings.

II. THE HOLOGRAPHIC MODEL
We consider a 4-dimensional holographic massive gravity model [32,33] defined by the following action: with X ≡ 1 2 g µν ∂ µ φ I ∂ ν φ I . The Stückelberg fields admit a radially constant profile φ I = x I which breaks the translational invariance of the dual field theory. For the rest of the manuscript, we focus on the specific potential V (X) = X 3 , which realizes the spontaneous symmetry breaking of translations and it gives rise to a finite elastic response in the dual field theory and to the presence of propagating phonon modes -the corresponding Goldstones [35,39,55,56]. 1 For the potential considered in this work, the sound speeds of transverse and longitudinal phonons are subluminal [35,39,57]. Moreover, the elastic response is accompanied by a viscous dissipative contribution [37], which qualifies the model as viscoelastic [36]. 2 One direct consequence of the competition between elasticity and dissipation is the observation of a sound to diffusion crossover in the spectrum of transverse phonons [40], analogous to the Ioffe-Regel crossover in dissipative systems [58]. In the linear regime -valid when the external deformations are small -we can use linear response theory to obtain the shear correlator from the bulk theory using the holographic dictionary. In the limit of zero momentum, the stress tensor correlator reads and it defines for us the storage modulus G (ω) and the loss modulus G (ω), together with the loss angle (phase shift) tan δ(ω) ≡ G (ω) G (ω) . At low frequency we have: where G 0 and η are the static shear modulus and the shear viscosity, respectively. In a perfect elastic solid, we have G = 0 and δ = 0, while in a purely dissipative fluid G = 0 and δ = π/2. All the materials with 0 < δ < π/2 are by definition viscoelastic. The larger the parameter m -the mass of the gravitonthe stronger the elastic component.

III. NONLINEAR RHEOLOGY
Whenever the amplitude of the applied strain is large, nonlinearities set in and the linear viscoelastic approximation fails. From a gravitational point of view, this problem requires a more complicated time-dependent setup which is explained in detail in the SM , following the seminal work of [59]. 3 Within this regime, the produced stress is no longer linearly proportional to the applied strain but it presents a distorted shape which can be understood as a superposition of different Fourier components. More specifically, in the nonlinear regime, the strain γ and the stress σ can be represented as 4 3 Importantly, out-of-equilibrium, temperature is not a well defined concept [60]. For this reason, we will indicate it everywhere with the symbol T in to specify that such parameter is the temperature of the initial (t = 0) state. 4 The reason why only odd powers appear in the expansion is that the stress response is typically taken independent of the shear direction.
where a 11 , b 11 correspond to the complex moduli G (ω), G (ω) in the linear regime, and the first nonlinear corrections entering at order O(γ 3 0 ). In this letter, we will explore different methods to represent and characterize the nonlinear response at large amplitudes: ( fig.1 that by increasing the amplitude of the applied strain the shape of the stress response gets distorted and it deviates from a simple oscillatory function. This behavior is also displayed in the corresponding Lissajous figures which are no longer a simple oval, as expected in the linear regime. We notice that the shape of the curve after each cycle appears to be slightly modified; this phenomenon emphasizes the complexity of our viscoelastic system. In fig.2 we study the Fourier spectrum of the signal. At small amplitude (orange curve), the spectrum is localized on the first and only harmonic, which is fixed by the frequency of the applied strain signal. This means the system is still in the linear response regime, where the stress is linearly proportional to the applied strain. By increasing the amplitude, higher (odd) harmonics appear in the spectrum confirming the functional structure displayed in eq.(6). The normalized power of the higher harmonics, I n/1 is shown in fig.3. We find preliminary evidence for a power law behavior ∼ γ n 0 which was previously suggested by theoretical arguments in [61]. Continuing along the lines of eq.(6), we can rewrite the stress response as: The normalized intensity I n/1 ≡ P (ωn)/P (ω1) of the first three higher harmonics in function of the strain amplitude γ0. Right: The first complex moduli G 1 (ω, γ0), G 1 (ω, γ0) at fixed frequency, in function of the strain amplitude. We normalized them by their linear value G(0) ≡ G(γ0 1). The colors correspond to m/Tin = 0.01, 1.81, 30 (from yellow to purple). The onset of nonlinearity is roughly independent of the value of m/Tin and it appears around γ0 ∼ 0.3.
The complex moduli are rigorously defined only in the linear regime; however, the measurements of G (γ 0 ) and G (γ 0 ) at a fixed frequency can provide meaningful information. The most common option to calculate the moduli from a non-sinusoidal response consists in looking at the quantities G 1 (ω, γ 0 ), G 1 (ω, γ 0 ), defined as the contributions from the first harmonics sin(2πωt), cos(2πωt) to the expansion in eq.(7). Additionally, the values of G 1 , G 1 are exactly what the commercial rheometers provide in the experiments. Using a simple expansion, we obtain the first term in the sum of eq. (7): where we are neglecting the higher harmonics corrections which naturally appear in eq.(7). Given these notations, the values of G 1 , G 1 at zero strain correspond to the linear response limit in eq.(3). We plot the dependence of the first nonlinear complex moduli G 1 (ω, γ 0 ), G 1 (ω, γ 0 ) at fixed frequency in fig.3. We observe that for small amplitudes the moduli are independent of the strain amplitude. This is not true anymore at large amplitudes where nonlinear effects become important. We find that the onset of nonlinearity is roughly independent of the value of m/T in and it depends solely on the amplitude of the applied strain γ 0 . Notice that in the nonlinear regime both moduli grow in a faster-than-linear fashion. From an operational point of view, this defines the presence of strain stiffening. This behavior is typical of hyperelastic materials such as rubber-like systems or complex polymers and it is in contrast to the so-called strain hardening which is on the contrary a common feature of rigid metals. Let us also notice that at low strain, for small values of the mass m, G 1 > G 1 , indicating that our dual field theory is a viscoelastic liquid. This is reversed at large values of m/T in , where the system becomes a viscoelastic solid with G 1 < G 1 [62] (see fig.(A.2) in the SM). This is totally consistent with the fact that the graviton mass To complete our analysis, we construct the Pipkin diagram of our model in fig.4 by plotting the Lissajous figures at various strain frequencies and amplitudes. We observe a neat transition between a linear viscoelastic regime at low amplitude and frequency to a more complicated large regions where the response becomes highly nonlinear (notice the similarities with [53]). This last result confirms that the regime we investigated cannot be described by linear response and it displays all the main physical properties of LAOS systems.

IV. DISCUSSION
In this letter, we characterize the nonlinear and time dependent mechanical response of viscoelastic (strongly coupled) field theories using holographic techniques. We focus our analysis on oscillatory external strains and on the regime of LAOS. We prove that the viscoelastic response is quite similar to that of complex fluids and hyperelastic materials (e.g. rubber) and in particular it exhibits a very neat strain stiffening phenomenon (see [63] for some concrete experimental data), which suggests that the holographic models at hand are not suitable to describe ordered metallic crystals characterized by strain hardening. We also observe a transition between a viscoelastic liquid behavior at small graviton mass, m/T 1, to a viscoelastic solid regime at large values, which confirms the identification of the graviton mass with the rigidity of the dual field theory. This work opens a new path for the study of complex flu-ids and viscoelastic systems using the holographic methods, which so far have been successfully applied only to strongly coupled liquids with no elastic response. There are several direct and interesting directions to pursue. Firstly, it would be desirable to reach a better theoretical understanding of our numerical data by comparing our results to known phenomenological models such as the multi-mode Giesekus model [61,64]. Secondly, a more extensive exploration of the phase diagram and possibly an extension of the study to include also the Chebyshev analysis [65] are certainly needed to draw universal conclusions. On a more phenomenological perspective, one could consider different types of experiments, i.e. different signals for the applied strain such as building up functions, step functions and quenches. This extension would permit the study of extremely interesting phenomena such as nonlinear relaxation, stress overshoot, yielding, which represent still open challenges for rheology and condensed matter in general. In this respect, the existence of possible universal (relaxation) timescales is certainly a fundamental question to answer. Finally, another relevant question which our work poses is the possibility of having holographic models display-ing strain hardening and being therefore more suitable to describe metallic solids. This nicely connects to a fundamental and still open question: which kind of solids do these holographic systems describe? As shown in this letter, the analysis of nonlinear transport properties is certainly a way of resolving this conundrum.

Linear response
We review the linear viscoelastic analysis of the holographic model used in the main text and defined from the bulk action in eq. (2). More details about this linear regime can be found in [36]. Given an applied external oscillatory strain γ(t) = γ 0 sin(ωt), the linear viscoelastic response is encoded in a time dependent shear stress of the form: where G (ω), G (ω) are the storage and loss moduli, determining the elastic in-phase response and the dissipative out-of-phase one.
Using the Kubo formulae formalism, these two moduli can be read from the shear stress tensor two-point functions as follows: The shear correlator can be derived with standard techniques using the holographic dictionary and considering the bulk dynamics of a gravitational wave perturbation -a geometric perturbation δg xy . Concretely, we consider a four dimensional asymptotically AdS black hole geometry in Eddington-Filkenstein (EF) coordinates: where u ∈ [0, u h ] is the radial holographic direction spanning from the boundary u = 0 to the horizon u = u h . The emblackening factor appearing in the BH geometry takes the simple form: and the corresponding temperature of the dual field theory reads: Additionally, the entropy density is given by s = 2π/u 2 h and the heat capacity can be directly obtained as c v = T ds/dT . As already explained, the linear viscoelastic properties of the model are encoded in the dynamics of the T xy operator which is holographically dual to the bulk perturbation δg xy ≡ h xy . At zero momentum k = 0, the linearized equation for the field h xy decouples and in EF coordinates reads: where V X ≡ ∂ X V (X) and primes denote radial derivatives. The UV asymptotic behavior of the h xy field is: Finally, we can read the shear Green function from G (R) which is the fundamental relation to characterize the linear viscoelastic response.
In the low frequency region: where G 0 is the static shear modulus and η the shear viscosity; in the opposite limit of large frequency: where G ∞ is the instantaneous shear modulus.
At small values of m/T , the static shear modulus can be computed analytically [35,37], using formula: At the same time, the instantaneous shear modulus is universally given by [36] : where ε is the energy density of the system. The only dimensionless and relevant physical parameter of our setup is the ratio m/T ; the physics of the system depends solely on its value. At small values of m/T , the viscosity of the system is larger than its elasticity -we are in a viscoelastic liquid regime. At an intermediate value of m/T ∼ O(1), the situation is inverted. The system becomes more and more rigid and it enters the viscoelastic solid regime. Notice that, in this section, since we are only considering the linear response regime, the parameter T is the well defined temperature of the dual field theory.

Details about the nonlinear computation and the holographic renormalization
The simplest time dependent metric in Eddington-Finkelstein parametrization with g xy = 0 and homogeneous symmetry is cosh (H(u, t))(dx 2 + dy 2 ) + 2 sinh (H(u, t))dxdy To implement the characteristic formulation [59], we use hyperbolic functions such that Volume xy = S 2 . This choice leads to the following equations of motion, where d + A :=Ȧ − A u 2 A is the directional derivative, and prime stands for derivative with respect to radial coordinate u and we use dot for derivative with respect to time t.
The metric functions have the following near boundary expansion 5 accompanied with the Ward identity, Since we want to impose the shear strain as the source for the boundary metric h xy = γ(t) and keep the spatial length scale of the boundary theory fixed h xx = h yy = 1, we use (A.21) Note that the strain function has to be less than one, in our units, to have real metric functions.
Using the standard holographic renormalization, one can read the boundary stress tensor from the near bound-ary expansion of the metric components as 6 The shear stress to the strain function h xy (t) = γ(t) is given by off-diagonal component of the energymomentum tensor T xy (t).

Numerical routine
To impose the boundary conditions, it is convenient to use the following redefined tilde functions In our numerical calculation we use the Chebyshev discretization with 50-100 grid points for integration along the radial coordinate. Following the characteristic feature of the bulk equations of motion in EF coordinates [59], here is the numerical routine to solve the equations of motion (A.14)-(A.18): 1. We start with a static black hole solution withH = 0,S = s 1 (t 0 ) as an initial configuration and choose a strain function γ(t).
2. We check the accuracy of numerical calculation, by pluggingH,S in the constraint equation (A.14).
3. We use the definition of apparent horizon, d + S(u = 1) = 0, as a boundary condition to calculate d + S by solving eq. (A.15). 6. By using the definition of operator d + we find theṠ,Ḣ. Then we integrate in time, by employing fourth-order Runge-Kutta method for the first three time steps and then the fourth order Adams-Bashforth method, to computeH(u, t 0 + δt) and S(t 0 + δt) and repeat the same routine from step 2.
To impose the sinusoidal strain, we turn on the amplitude smoothly (in the spirit of [66,67]) as where parameters t c , w c control how accurate the initial configuration satisfies the constraint equation (A.14) and how fast the maximum strain is reached respectively. In fig.(A.3) we show the stress, the energy density and we check the constraint equation (A.14) for a concrete example.

More about the nonlinear response
In this last section, we give some more details about the nonlinear response. First, we show the shape of the 3D phase space {σ, γ,γ}. In the top panel of fig.A.4, we display it for a large value of the strain, within the LAOS regime. The distorted shape of the curves there are indeed the confirmation that nonlinear corrections are important in size. In the same figure (bottom panel), we emphasize the transition into the nonlinear regime by plotting the 3D curve for small, intermediate and large strain. The difference between the linear green curve and the large amplitude black one is evident. The Lissajous figures shown in the main text are very good indicators of the viscoelastic behavior and they contain important information about the dissipative mechanism in the system. More specifically, the area of the figure coincides indeed with the dissipated energy over the cycle: where C indicates that the integral is performed over a specific cycle. We plot the dissipated energy in function of the strain amplitude in fig.A.5. We observe that increasing the amplitude of the applied strain the dissipated energy grows. Moreover, at low strain, a scaling regime appears, where E ∼ γ 2 0 . Finally, we discuss in more detail the nonlinear response in terms of the complex moduli and the entropy production. Given a Lissajous figure, we can define two impor- tant quantities: where G M is the small strain (tangent) modulus and G L the large strain (secant) modulus. A valid indicator of the nonlinear response is the difference between those two values: Whenever S > 0, the nonlinear response exhibits strain stiffening -a larger elastic response at larger applied strains. We generally observe that our system displays this feature independently of the parameters (see top panel of fig.A.6 for a concrete choice). This confirms the analysis made in the main text regarding G 1 and G 1 . The behavior of the n th complex moduli is repeated in the bottom panel of fig.A.6. As already shown in the main text, at low strain amplitudes G 1 , G 1 are roughly independent of γ 0 and they correspond to the linear moduli G , G . At larger strains they start growing signaling the strain stiffening phenomenon. On the contrary, the higher harmonics moduli G n are zero at small strain amplitude -the response is fully linear. They then grow, rendering the stress response highly nonlinear. To conclude, we investigate the behavior of the entropy density during the nonlinear dynamics. The results are shown in fig. A.7. We observe that the entropy density grows rapidly with time, as indicated by the growing of the BH horizon area. Moreover, we conclude that the entropy density grows in function of the applied strain amplitude. At least for intermediate values of the amplitude, the growth is very well approximated by a quadratic function ∼ γ 2 0 . At this stage, it is not clear to us how to relate this behavior with the strain stiffening that we observe in the nonlinear viscoelastic response.