Black Rubber and the Non-linear Elastic Response of Scale Invariant Solids

We discuss the nonlinear elastic response in scale invariant solids. Following previous work, we split the analysis into two basic options: according to whether scale invariance (SI) is a manifest or a spontaneously broken symmetry. In the latter case, one can employ effective field theory methods, whereas in the former we use holographic methods. We focus on a simple class of holographic models that exhibit elastic behaviour, and obtain their nonlinear stress-strain curves as well as an estimate of the elasticity bounds - the maximum possible deformation in the elastic (reversible) regime. The bounds differ substantially in the manifest or spontaneously broken SI cases, even when the same stress-strain curve is assumed in both cases. Additionally, the hyper-elastic subset of models (that allow for large deformations) is found to have stress-strain curves akin to natural rubber. The holographic instances in this category, which we dub black rubber, display richer stress-strain curves - with two different power-law regimes at different magnitudes of the strain.


Introduction
The response of materials under mechanical (or elastic) deformations is a basic aspect of matter, which is important to understand and characterize. This is an old field of study because of the many practical applications and a large part of it is well understood since long ago [1,2].
The response is best understood when restricted to the 'linear' regime (small deformations), but there are many examples of solids that can undertake large deformations [3]. Common examples of these are the rubbers, but more generically they are referred to as hyper-elastic materials. The non-linear response that these materials exhibit is encoded in the stress-strain relations -the amount of constant stress that must be applied in order to deform by a certain amount the material -see Fig. 1 for a prototypical example. These curves can be easily obtained from experiments, but they are usually difficult to compute from the microscopic ingredients (even within the reversible regime, that is neglecting plasticity and dissipative effects), especially so in strongly coupled materials. Moreover, the nonlinear response is characterized by a rather large number of parameters/observables (e.g., all the derivatives of the stress-strain curve at the origin, the maximum strain that the material can undertake, etc), and it might well be that there exist correlations between them. This motivates the study of the nonlinear mechanical response using effective low energy methods, which on their own might capture some of these correlations and even how these parameters depend on external parameters (temperature, etc).
The 'continuum limit' description of mechanical deformations represents one such effective method that is useful for nonlinear response. This approach embodies already a broad literature, see [4] for a review. As in hydrodynamics, the medium is described by a 3-vector π i (t, x j ), the displacement vector of the solid elements. How the material deforms is encoded in its gradient, the so-called strain tensor ij ∼ ∂ (i π j) . The main difference between solids and fluids in this language is that a solid responds to a constant external stress σ ij with constant ij , whereas a viscous fluid responds with a constant strain rate,˙ ij . The punchline, though, is that the same kind of effective description is possible both for fluids and solids at small frequencies (and momenta) [5].
For small applied stresses the response is linear, σ ij ∝ ij , and the proportionality constants are usually called elastic moduli. Nonlinear elasticity concerns the relation between the stress and strain tensors beyond the linear approximation. The mathematical formalism required for this in the continuum limit was developed long ago, see [4] for a comprehensive review. This results especially relevant for hyper-elastic materials or elastomers, which allow for large reversible deformations. For them, the continuum limit description takes the form of a nonlinear theory for the displacement vector field π i that can be specified by an energy function (how the energy density depends on ij ) or a (nonlinear) constitutive relation.
Symmetries allow to translate nonlinear elasticity into quantum field theory language. Given that condensed matter breaks spontaneously spacetime translations and boosts, it is possible to catch the dynamics for the lightest degrees of freedom using the methods of Effective Field Theory (EFT) applied to the spontaneous breaking of these symmetries. In solids, the Goldstone bosons associated to this spontaneous breaking can be identified as the ('acoustic') phonons [7,8]. These phonons are the fluctuations in the displacement vector π i , they are indeed gapless and therefore belong to the lowest lying excitations, which makes the whole EFT construction consistent.
The possible form of the full nonlinear effective Lagrangian for the phonon fields can Figure 1. Cartoon of a typical stress-strain curve of a hyper-elastic material [6]. In light blue shade the linear regime, in which stress ∝ strain. At large strain deformation, the stress-strain relation can display a power law behaviour, σ ∝ ν , with some exponent ν. The light red area illustrates this behaviour with ν > 1. Materials typically break or fracture after some critical deformation, which translates in the stress-train curves terminating at some point. The red star indicates the breaking point.
be obtained using the coset construction applied to spacetime symmetries, see [9,10]. It was recognized in [11] that, for solids, the resulting phonon effective Lagrangian takes precisely the same form as the continuum limit nonlinear elasticity theory for the displacement vector π i (t, x j ) for hyper-elastic materials [4] (at leading order in derivatives). The crucial advantage is that the solid EFTs are full-fledged effective Lagrangians that include dynamical effects and relativistic corrections among other improvements [11]. Let us emphasize that, even if they are not formulated in terms of the microscopic degrees of freedom, the EFT methods have stringent predictive and constraining power. This point was illustrated in [11], where the correlation between various nonlinear observables was made manifest in the form of elasticity bounds -limits on the strain that a material of certain type can possibly withstand depending on other properties of the material. The simplest example arises by considering a class of materials characterized by power-law scaling in the stress-strain curve, schematically, with some arbitrary exponent ν. Some elastomers in nature follow such a power-law at large deformation [4] with a variety of exponents. For a general analysis of the elastic response, one can treat ν simply as an effective parameter to describe the nonlinear response (at least in a class of materials). Interestingly, assuming this power-law response is enough to place a priori an upper limit on the maximum strain deformation, max , that the material can undertake [11]. This maximum deformation max plays the role of (an upper limit to) the mechanical breaking or failure point of materials.
A nontrivial outcome of the EFT methods is that one can establish a relation between these two nonlinear elasticity parameter, ν and max . Let us emphasize two points in order to highlight the potential value of the EFT methodology for nonlinear elasticity. First, we stress that the correlation between ν and max is entirely based on low energy EFT properties. This illustrates that it is possible to understand some of the properties of the nonlinear response just from the low energy theory, that is, independently from the microscopic details. Second, the constraining power of the EFT methods is expected to apply to many other nonlinear elasticity parameter, beyond the one examined in Ref. [11]. This is especially clear taking into account that the main benefit from the EFT methods is that the full nonlinear structure of the theory is fixed by the symmetries.
This encourages us to continue the analysis to other cases, in particular to the more sophisticated materials that exhibit scale invariance (SI). As elaborated in [12], this case deserves a special treatment, because SI can be realized in several ways and this affects the elastic response -even in the linear response regime. 1 The main distinction concerns whether SI is a broken or a manifest symmetry of the low energy dynamics. The latter case implies that the theory that describes the solid excitations must be akin to a conformal field theory (CFT). In this case, bottom-up AdS/CFT methods [15] provide a useful tool to properly model the material. The opposite case -with SI as a spontaneously broken -can be treated instead within the more conventional EFT methods [7][8][9][10].
The existence of these two types of SI solids gives a 'bonus' of motivation to the present study, as it is interesting to compare how much the low-energy constraints in the nonlinear response differ depending on whether SI is manifest or spontaneously broken. As we will see, there is a significant difference in the relation between the nonlinear parameters introduced above (ν and max ) for the two types of SI materials.
A good fraction of this work is devoted to provide the tools to compute the nonlinear mechanical response for the manifest SI case, by exploiting holographic AdS/CFT methods. The main technical development that we present is the construction of a large family of asymptotically AdS black brane solutions that are subject to finite mechanical deformations 2 and we obtain their stress-strain curves. We shall find that a certain class of models allows for black branes that can be deformed elastically by large amounts without breaking. In these cases, their stress-strain diagrams are similar to that of natural rubber (with O(1) values of the exponents ν, see below), so isn't a great stretch to call these solutions black rubber.
These solutions can be found semi-analytically in the simplest models, which include the displacement vectors π i as new explicit degrees of freedom -also called Stückelberg fields in the previous literature. This paper builds on the recent holographic massive gravity models [18][19][20][21][22][23] which realize in a simple way the spontaneous breaking of translational invariance in 'critical' materials (with manifest SI at low energies).
More recently, several works have improved the framework to accommodate for the spontaneous breaking of translations [21] and the study of the linear elastic response [22], the vibrational modes of the systems [24,25], their viscoelastic nature [26][27][28][29] and their hydrodynamic and physical description [23,30].

Nonlinear Elastic response
In this Section, we review the basic formalism to describe the elastic response under finite ("large") deformations or applied stresses. Linear elasticity theory describes how materials deform in presence of a small ("infinitesimal") external deformation. The mechanical deformation for a solid in d + 1 dimensions can be described by a d-dimensional vector field, the displacement vector, that measures the physical distance from equilibrium position at any given point in the solid. Out of the π i , one builds a rank-2 symmetric tensor, called the strain tensor as Volume-preserving deformations satisfy and are called shear strain. Similarly, strains that change volume but not shape satisfy ij ∝ δ ij (2.4) and they are called pure bulk deformations. In most materials, at small enough deformations (small strains), there is a linear relation between the stress needed to apply on the material and the generated strain. Mathematically, this translates into a linear relation between ij and the stress tensor σ ij of the form, The elastic tensor, C ijkl , is well known to reduce to just two parameters for homogeneous and isotropic materials: the shear and bulk moduli, which encode the linear response to pure shear and pure bulk deformations respectively. Non-linear elasticity concerns the relation between stress and strain beyond linear level -conceptually, the full functional form of σ ij = σ ij ( kl ). In order to extend the relation to the non-linear regime, one must pay attention to how the strain deformations are defined nonlinearly.
Reviewing the logic, one realizes that in materials that are homogeneous and isotropic at long wavelengths there symmetry allows to choose what we call solid elements so that their equilibrium positions coincide with the 'cartesian' coordinates. This suggests to introduce another variable to describe the state of deformation, so that equilibrium corresponds to φ I eq = δ I i x i . This variable is also more amenable to treat homogeneity and isotropy as an internal symmetry for the scalar fields Φ I , as done in [7,8], which is why the index label on π I has been capitalized.
A general state of deformation that is constant along the material is then given by with an arbitrary constant matrix O I j , which is a useful way to parameterize the strain tensor for finite deformations. One can easily convince oneself that in the homogeneous and isotropic limit one can restrict O I j to be a symmetric matrix with no loss of generality. Isotropy forbids the presence of any background shear strain.
The advantage of using (2.6) as a variable is now clear: the natural extension that supersedes (2.3) to the nonlinear regime is This condition extends to non-linear level the requirement that the deformation described by the matrix O I j does not change the volume of the system.
For illustration, in 2 + 1 dimensions, O I j is simply a 2 × 2 symmetric matrix, which contains only 3 free parameters. We shall stick to the following parametrization, where ε serves as a nonlinear version 3 of the shear deformation and α for the pure bulk deformations. We are dropping the 'third' parameter,ε, for deformations of the form O I j = diag(eε, e −ε ) because they only differ from ε deformations in that the shear is introduced in a basis rotated by 45 degrees. Since we are assuming homogeneous and isotropic materials, it suffices to consider one of the two shear 'polarizations'.
The stress-strain curves can now be extracted by computing the stresses in the material that are necessary to support a configuration (2.7). Once the low energy theory for the material is specified, this reduces to just looking at the stresses produced by these configurations. Continuing in the 2 + 1 example above, this can be done by computing the stress as a function of the deformation O I j . Materials that can be deformed by large amounts while in a reversible fashion are generically called hyper-elastic. For these, the stress-strain relation can be obtained from a so-called energy function scalar function [6], that characterizes how energetically 'expensive' every deformation is. From now on, we will consider materials with this property and which realize scale invariance (SI), and we will distinguish between two sub-cases depending on how SI is realized: i) critical solids, that is, which realize SI as a manifest symmetry at low energies (Section 3); ii) solids with spontaneously broken SI, with a gapped spectrum and thus have phonons as the lowest-energy excitations (Section 4).
Case ii) can be dealt with using the solid EFTs, so part of the discussion was already presented in [11]. Here we will extend the analysis with the aim at the SI case and its comparison to the manifest SI case.
As a warm-up, let us remind now how the computation of the stress-strain curve proceeds for a general solid EFT (not necessarily assuming broken SI). Restricting to 2+1 dimensions for simplicity, it can be seen that the most general effective Lagrangian at leading order in derivatives can be written as [9,11] with X and Z defined in terms of the scalar fields matrix 4 I IJ = g µν ∂ µ φ I ∂ ν φ J as X = 1 2 tr I IJ , Z = det I IJ . It is immediately clear that once one restricts to the (strained) static and homogeneous configurations given by (2.7) and (2.9), the action (2.12) plays exactly the same role as the 'energy function'. In other words, we can identify We emphasize that the nontrivial content in the Solid EFT construction is that once one knows the 'energy function' then the whole effective Lagrangian is also known, which can then be used to extract more information such as the elasticity bounds [11]. Instead, for the solids with manifest SI of Section 4 the energy function still exists but it does not correspond directly to the effective Lagrangian -in fact in these cases one expects that a local Lagrangian doesn't exist.
The stress required to support the configurations (2.7) can be read off from the stress tensor associated to these configurations (2.7), which can be easily computed in the EFT. For any time independent scalar field configurations, the stress-energy tensor components are [11] T tt ≡ ρ = V , (2.14) where V X ≡ ∂V /∂X, etc. The deformed field configuration (2.7) introduces both shear and bulk deformation. Setting α = 1, it describes a pure shear strain (i.e. volumepreserving) in the (x, y) directions induced by ε = 0. The full nonlinear stress-strain curve is then found be expressing the stress T xy as a function of the strain ε, These results apply to any solid whose low energy dynamics can be treated with EFT methods. This includes the solids with spontaneously broken scale invariance (SI), which we discuss in more detail in Section 4. However, these steps are not justified for solids which exhibit manifest scale invariance [12].
In principle the procedure is identical for solids with manifest SI, we just want to obtain how much stress is required to support a configuration with given strain ε. However, the main obstacle is that, as in CFTs, scale-invariant solids are expected to lack a local Lagrangian description, therefore the steps after (2.12) do not immediately apply (nor the identification of the 'energy function' with an effective Lagrangian). While this may seem unimportant regarding the response to static and homogeneous strain, it is crucial in order to possibly obtain nontrivial constraints in the nonlinear response (such as the correlations among various nonlinear parameters mentioned in the introduction) because this requires a knowledge of the full theory.
In the next Section, we show how to extract stress-strain curves in (holographic models of) solids with manifest SI, we shall work out the equivalent of Eq.(2.17) for them, and find the constraints and relations among different nonlinear elasticity observables.

Solids with manifest scale invariance
As mentioned in the introduction, our main focus are materials in a critical regimewhich exhibit manifest scale invariance at low energies. We shall model them using the standard holographic dictionary. As usual, it simplifies the analysis to model the scale invariant field theory as a deformation of a CFT. In this case, the AdS/CFT dictionary tells us that the material, which we assume is 2+1 dimensional, is going to be described by asymptotically AdS 4 planar black brane solutions. By assumption, the CFT contains operators that can be identified with the displacement vectors. Their dual incarnation in the AdS 4 , are an identical a set of fields, φ I , which propagate into the holographic dimension too. See [22,31] for more details.
The way to extract the stress-strain curves in these models is simply to find the black brane solutions with nontrivial strain tensor 'emanating' from the horizon. The strain tensor, then, can be thought of as an asymptotic charge of these black branes. Keeping track of the stress tensor for each strain tensor 'charge', one can compute the strain-stress curve.
Given that for every (constant) value of the strain tensor there is a static black brane solution, the process of varying the strain (which is implied in the stress-strain curves) can be assumed to be a reversible process, if done slowly enough. For this reason, we will treat the (static) nonlinear elastic response of these black branes as elastic (reversible). This is, of course, until some instability is reached -and this is the basic guiding principle we shall use to establish elasticity bounds.

Nonlinear response for holographic models
We consider the holographic massive gravity models introduced in [21,22] (see also [24][25][26][32][33][34]), which are obtained by introducing displacement fields Φ I in the AdS bulk, with a generic action Tr(X IJ ) and Z ≡ det(X IJ ). For simplicity, we focus on d = 3 but we the construction can be easily extended to higher dimensions.
For specific choices of the potential W (X, Z), the model (3.1) represents the gravity dual of a CFT at finite temperature and where translational invariance is broken spontaneously. Using the standard AdS/CFT dictionary, this defines for us a CFT that will have non-zero elastic moduli and so it can be interpreted as a model for a solid in a quantum critical regime. We remind the reader that throughout all the manuscript we will only consider standard quantization for the scalar bulk fields Φ I . More precisely, under these assumptions, a well-defined elastic response can be defined for potentials which decay at the boundary as W ∼ u 3 or faster [24] 5 . Moreover, for potentials whose fall-off at the boundary is W ∼ u 5 or faster this elastic response is associated to the presence of massless propagating phonons [25]. The most important point for the moment is that the gravity theory also contains a field Φ J , which is directly linked to the material deformation.
It has been shown before [36] that there exist simple homogeneous asymptotically AdS planar black brane solutions with which from the gravitational perspective acts as a "solid hair " or more technically as magnetically charged 0-forms [37]. Their CFT interpretation fits that of a critical 2+1 planar and homogeneous solid material, with broken translations. How to perturb this solution and read-off the (linear) elastic moduli has already been discussed at length previously [22,24,26]. Our next goal is to find the holographic stress tensor carried by strained configurations (strained solids) with O I j given in (2.9) (with finite ε and α). Since the CFT stress tensor is dual to the metric and we are after the full nonlinear response, we must find the exact holographic stress tensor produced by the deformation 'source' O I j . In practice, this implies that we must find (asymptotically AdS black brane) exact solutions to the Einstein + scalars theory with a nonzero tensor mode -the strain tensor. That is, the spatial part of the metric g ij (with i, j running over x, y) must differ from ∝ δ ij so that it contains a shear (and bulk) deformation.
Fortunately, for deformations that that are constant in time and space it is possible to reduce significantly the equations 6 . Indeed, one can see that the full system of nonlinear Einstein equations can be solved in this case going to the following ansatz where γ ij is a d − 1 dimensional spatial metric with unitary determinant. In d = 3, one can parametrize γ ij in terms of the usual + and × polarizations as where h +, × are functions of u only and σ +, × stand for the Pauli matrices that are usually called σ 3, 1 respectively. The two polarizations h +,× couple to each other at nonlinear level. In order to disentangle them, it is convenient to switch variables to where again h and θ are functions of u only. In these variables, h is the magnitude of the spin-2 mode and θ the direction in the space of polarizations. For physical solutions, h(u) must have a vanishing leading mode -and its subleading mode encodes the stress tensor. As we will see shortly, it follows from the equations of motion that in these solutions the function θ(u) must be a constant. In these variables, then, θ will simply encode the polarization direction of the stress tensor and h the magnitude of the response. We can do a similar representation for the strain matrix O I j : The constant α parametrizes the bulk strain deformation whereas the constants Ω and θ 0 encode the strain magnitude and polarization, and they are related to the nonlinear shear strain parameters ε andε introduced in Section 2. For instance, for θ 0 = 0, one has ε = 2 sinh (Ω/2) .
The magnitude of the shear strain encoded in Ω acts as a source term for the metric in the bulk.
Before showing the equations of motion, note that in homogeneous and isotropic material the elastic response is such that the strain and stress tensors are aligned in the same polarization direction. In our notation, this translates to having θ = θ 0 . We shall see shortly that this is indeed the case in for the physical solutions, but for the moment we keep θ(u) generic in order to see how it is determined by the equations of motion.
In d = 3, the independent equations for the background (3.4) are: where the cosmological constant is fixed to Λ = −3 and we indicate with the subscript (h, θ), the derivative with respect to h and θ. The potential W is evaluated on the background values, which gives Within this ansatz the form of χ(u) and f (u) are completely dictated by h(u) and θ(u). We assume the presence of an event horizon at u = u H defined by f (u H ) = 0 and h(u) reaching a constant value at the horizon, h(u H ) = h H . The associated entropy density is s = 2π/u 2 H and the corresponding temperature reads T = − f (u H ) 4 π e −χ(u H )/2 . In the asymptotic UV region we impose f (0) = 1 and χ(0) = 0.
Assuming that the mass term m 2 W h vanishes sufficiently quickly towards the asymptotic UV region, one finds that the two independents modes of the spin-2 metric deformation are where dots represent higher powers of u. As usual, the subleading term C 3 is identified via the AdS/CFT dictionary with the VEV of the stress tensor T xy , and C 0 with an external spacetime metric source for T xy operator (see below). Throughout all the manuscript, we will consider that the external spacetime deformation source is absent, so that C 0 = 0 . (3.14) With this condition, the only possible source for the stress tensor arises from the mechanical strain deformation that is encoded in the scalars Φ I -ultimately in the parameters α, Ω and θ 0 that parametrize the strain deformation.
We can now look at the equation for θ(u), (3.12). Due to the cross-coupling to h in the second term of (3.12), the two modes of θ(u) near the AdS boundary turn out to depend on the boundary condition assumed for h, i.e. on the choice of C 0 . For C 0 = 0, θ would have a constant mode and a u 3 mode. However, for C 0 = 0, the asymptotic θ modes are the constant mode and a u −3 mode. Regularity at the AdS boundary (rather, consistency with AdS asymptotics) then requires to set the u −3 coefficient to vanish. Given that one is limited to only one free parameter, regularity at the horizon then is then expected to select the θ(u) = θ 0 as the only viable solution at least for generic choices of the potential W . Incidentally, this closes the proof that the elastic response is isotropic, since a strain in a given polarization only sources stress tensor in the same polarization (this would not be true if θ(u) had a nontrivial profile). This was completely expected, but it is easy to show expliticly in the variables (3.6).
Therefore, from now on we will set θ = θ 0 = 0 for the rest of the manuscript. This simplifies the set equations of motion substantially, which makes clear that Ω acts as a 'source' term for h in the bulk. These two equations (3.15) determine uniquely the profiles for f and h and then χ(u) is obtained by integrating (3.9), which reduces to 2 χ = u h 2 .

General results
From the point of view of the gravity theory in the bulk, the solutions to (3.15) can be viewed as black branes with a form of hair encoded in the scalar configurations with nontrivial strain (2.9). To fix ideas, we can think that there are 2 parameters (or charges) that label the solutions: the magnitude of the scalar gradient at zero shear strain, α, which we assume is nonzero (and which can be traded by m for monomial potentials); and the magnitude of the shear strain ε. (We ignore now the angle θ 0 since it only sets a direction.) The novelty of the solutions presented here with respect to the ones previously discussed e.g. in [21] is that we will keep track of how the finite shear strain ε = 0 deforms the solutions. Picturing the shear strain as a standard charge that the black branes can be endowed with is also useful to understand their behaviour and properties. It is clear from (3.4) that we are constructing solutions with a non-zero and static tensor mode of the metric, h(u). This is possible for two reasons: first, because the strain tensor encoded in the scalars acts as a source for the tensor mode h(u); second, because the tensor mode h is a massive graviton. Indeed, the presence of a mass term in the equation of motion grants the possibility to have static response to a static homogeneous source.
At the level of understanding the stress-strain curves that will follow, it is clear that increasing the shear strain one must reach extremal (T = 0) solutions. Also, by changing u H together with ε, it is possible to construct one-parameter family of solutions, say, at constant temperature. Labeling these solutions by the amount of strain ε, and computing the shear stress for each solutions then we can obtain the strain-stress curve. This is the strategy that we follow in this work. Let us now summarize two general results that follow from this prescription.
First of all, it is possible to obtain an approximate expression for the stress-strain curve implied by (3.15) for shear deformations (that is with α = 1). The main observation is that the m 2 term in (3.15) acts as a source term for h for Ω = 0 and that at either m = 0 or Ω = 0, h(u) = 0 is a solution -which means that the stress vanishes in this limit. The profile h(u) is guaranteed to be small then for small m (at least for a class of potentials), and this allows for a perturbative scheme even for large deformations, Ω 1 (equivalently, ε 1). Following [26], we can treat m as a small parameter and find the solution order by order in m 2 . At first order in m 2 this gives This formula is of course only applicable if the integral is finite, which translates into some constraints on the functional form of the potential W (X, Y ). We will elaborate more on the constraints on W (X, Y ) in the coming subsections. At this stage, let us make two remarks. First (3.17) is perfectly consistent with the perturbative expression for the shear modulus found previously in [26]. Second, the formula (3.17) resembles structurally the one arising from EFT methods (2.17) but it also differs substantially in various ways: the bulk potential enters inside an integral which hints at a sense of non-locality; additionally, (3.17) encodes a non-trivial temperature effect, via the dependence on the location of the horizon, u H (which was obviously absent in the EFT formalism). Since the expression (3.17) is obtained perturbatively in m 2 , the factors next to m 2 should be not too large in order to be valid and at large enough ε one expects that this approximation should fail. However, as we will see below the (3.17) still gives a good approximation even for moderately large ε.
The form of the stress-strain curve at very large strain can still be obtained from a different consideration. The behaviour at asymptotically large strain can be understood from the structure of equations of motion (3.15). The key point is that large strain implies h(u) 1 and in this very anisotropic regime the equations (3.15) have an attractor solution (towards the UV), different from AdS 4 . In a subset of models (defined by the choice of the potential W (X, Z)), this UV attractor solution is actually a fixed point that realizes scale invariance with anisotropic scaling in the spatial x , y directions as well as in time. The presence of this additional anisotropic Lifshitz UV fixed point translates into the appearance of a second power-law scaling behaviour at asymptotically large strains. For the sake of clarity, we postpone the discussion of this point to Sec. 3.3 in the context of a specific benchmark model.
Finally, the models (3.1) admit extremal solutions of the form (3.4), whose nearhorizon geometry are then AdS 2 × R 2 . This represents yet another additional emergent scale invariance, this one with with Lifshitz dynamical exponent z → ∞, and isotropic character. This scaling is expected to be manifested in the lightest excitations, governed by the near horizon geometry.

A benchmark model
In order to make further progress, the form of the potential W must be specified. In this work we shall not try to find what is the form that matches the mechanical response of some known materials. Rather, we shall take an approach similar to [12], where we concentrate on potentials W that give rise to power-law stress-strain relation σ ∼ ε ν .
For this purpose we consider a benchmark potential of the form: In order to ensure the consistency of this choice (3.18), and of the model (3.1) in general, we ask the following requirements: absence of ghosts, absence of gradient instabilities, and positivity of the linear elastic moduli in the backgrounds with vanishing shear stress, ε = 0. These conditions constrain the range of the parameters a, b in the benchmark model (3.18) as follows [12], Moreover, restricting the discussion to the case of standard quantization, we need to impose also b > 5 2 to ensure the presence of massless phonons. Notice that the constraints considered are mostly bulk requirements and they represent just necessary but not sufficient conditions for the full consistency of our boundary field theory. In order to have a final verdict a detailed QNMs computation would be needed. For simple enough theories (of the form W (X) = X n or W (Z) = Z m ), that has been done in [23,25,38,39].
At low temperatures, we will find two distinct regimes with a scaling relation σ ∼ ν (both for shear and bulk strain deformations), as seen e.g. in Fig. 2. For the sake of clarity, then, we will introduce the following notation for the corresponding exponents: where the 1, 2 subscript denotes the regime encountered at lower and higher deformation regimes respectively; and S , B stand for the pure shear and pure bulk deformation sectors.

Shear deformations
The nonlinear shear response is encoded in the shear stress strain curve σ(ε), where the stress is given by (see Appendix A): 20) and C 3 is the subleading term in the UV expansion (3.13). The strain on the contrary is produced by the difference between the background configuration (3.7) and the equilibrium one Φ J = x J and it reads: In simple words the shear stress-strain relation is derived by the identification: where C 3 is extracted numerically varying Ω. The solutions can be easily obtained by shooting method, integrating (3.15) from the horizon towards the UV boundary. The boundary conditions that ensures regularity at the horizon are 2 u 2 h f − 4 m 2 W h (X,Z) u H = 0 and h(u H ) = h H , a finite constant that is used as the shooting parameter.
The results for the non linear shear response are shown in Fig.2 for an illustrative choice of potentials. Clearly, at small strains ε 1, the response is linear and the slope is given by the shear modulus studied previously [26]. Moving away from the linear approximation, we can notice that the stress-strain curve exhibits two different scaling regimes. At intermediate strain, a power law behaviour σ ∼ ε ν S 1 appears for ε 1 (but not too large), with ν S 1 = 2a . Additionally, for much larger strains ( ε 10 2 − 10 3 in the examples shown in Fig. 2), the curve again displays a scaling σ ∼ ε ν S 2 with a different exponent The manifestation of two scaling regimes actually happens only at high enough temperature; at low temperature the curve interpolates from the linear regime directly to the ν S 2 scaling as shown in Fig.2. As shown later, both scalings can be obtained analytically (see formula (3.17)).
Let us now supplement the numerical results shown in Fig. 2 with some analytical understanding. Performing the integral in (3.17) for our benchmark potential (3.18), we find the approximate expression Convergence of the integral requires b > 3/2, the same condition of the positivity of the linear bulk modulus. The linear limit ε 1 of formula (3.17) agrees perfectly with what was found previously in [26]. This is an approximate expression at leading order in m 2 , therefore it is valid so long as all factors are of not too large, which allows for finite but moderate ε. Still, the large ε limit of (3.25) already catches the first scaling behaviour with exponent ν S 1 given in (3.23). At larger strain, though, Eq. (3.25) is not expected to hold and more effort is needed to understand what should happen.
At asymptotically large strain, ε 1, it is still possible to obtain the non-linear scaling at large strain analytically by analysing the equations of motion (3.15). Large strain implies Ω 1 and so the tensor mode profile h(u) is expected to perform a large excursion or, equivalently, to have sizeable gradient h (u). It is convenient to introducẽ because the potential depends on h(u) only through this combination. At large Ω, near the boundary this quantity is large, and moreover one can expect thath (u) is also large (say, compared to f (u)) somewhere in the bulk.
The largeh 1 approximation allows to substitute cosh(h) and sinh(h) by eh/2. If the potential W (X, Z) is power-law-like as in our benchmark model, then it is easy to see that in this regime that the solution to the equations of motion (3.15) reaches a constant for f (u) f 0 < 1. In fact, our benchmark potential (3.18) admits a simple approximate solution of the form near the UV (u → 0), with m 2 u 2b 0 = 3 b 2 a /(a 2 + b). The logarithmic shape ofh(u) is clearly seen in Fig. 3. Before discussing further the properties of this approximate solution, let us first see how it determines the elastic response at asymptotically large strain. In this limit, the functionh(u) at the UV is large but finite,h(0) = Ω. Therefore, at some scale u * , there must be a transition between the constant and the logarithmic profile (3.27).
Assuming that (3.27) is correct up to the horizon, we have that approximatelyh(u) Ω − 2 b a log u u * from u * to u H . This allows to identify u * as where in the last step we use that Ω h (u H ) which is reasonable since theh variable is massive and 'tries' to reach 0. Thus, in the large strain limit, Ω 1, we have u * → 0. In other words, the intermediate solution (3.27) extends up to very close to the boundary. Since from the UV viewpoint u * represents the crossover scale wherẽ h changes from constant to logarithmic behaviour, just from dimensional analysis one expects that the subleading term in the AdS UV expansion (3.13), C 3 , must scale like The shear strain is defined as ε = 2 sinh( Ω 2 ), thus the shear stress for a large shear strain will scale with ν S 2 = 3 a b as indeed shown in fig.2.
Let us return now to the solution (3.27), which implies that the bulk metric in this limit takes a special form. Going to the coordinatesx,ỹ where the metric is diagonal, the bulk asymptotic geometry is This geometry represents a new UV fixed point. It is attractive towards the UV, but configurations with large enough strain get close to (3.29) for a long range in log u. Interestingly, it exhibits both a Lifshitz dynamical exponent (between space and time directions) as well as a Lifshitz scaling with respect tox andỹ spatial directions -that is an anisotropic scaling. It is quite natural that for large strain the x and y directions become very anisotropic. One might find more surprising that this occurs via a combined Lifshitz scaling both in the space and time directions. One way to understand this is by recalling that the models (2.12) can be thought of as massive gravity theories. Indeed, in the unitary gauge Φ I = δ I j x j one has X = tr g ij and Z = det g −1 , so the scalar kinetic function W (X, Z) becomes a potential for the (spatial part of the) metric. Recall that the simplest way to obtain a Lifshitz geometry is to support it with a massive spin-1 vector field [40][41][42]. It is not so surprising, then that a massive spin-2 field can do a similar job.
Note that the existence of the anisotropic Lifshitz UV fixed point solution (3.29) depends slightly on the choice of the potential W (X, Z). Looking at the first of Eqns. (3.15), one sees that a necessary condition for f (u) =const. to be a solution is that h(u) is logarithmic. In turn, this gives room for W either to approach a constant or to vanishing towards u = 0. In order to be consistent with the second equation, this requires that W h (which at large h translates into XW ,X ) is constant along the solution. Since in these solutions Z = u 4 , the only way that XW ,X can reach constant is that it depends on X, Z via the special combination X p Z q with and p , q constants which is such that every power of Z can be compensated by the X-dependence. This is precisely what happens in our benchmark model 7 (3.18). However it also happens in much more general choices of the form W X, Z = W 0 X p Z q , with W 0 a free functions of one argument. (An additive function of Z that vanishes at small Z would lead to the same behaviour.) It is also clear that in a more generic potential that does not reduce to this form the anisotropic Lifshitz UV fixed point solution is not present. Still, one expects that some anisotropic solution (not completely scale invariant) should exist and dominate the response in the regime of asymptotically large strains.
In any case, the (near-)Lifshitz form of the geometry is expected to impact be seen also in transport properties like the electric conductivity at finite strain (see for example [43]). Therefore, this suggests an avenue to possibly test whether this anisotropic Lifshitz regime occurs in real materials.
Anisotropic models with a uni-directional scalar field φ = αz have been already used extensively in the literature [44][45][46]. These models are expected to share some of the features of our solutions with large strain, as these can be represented by the scalars in the configuration Φ x = α e Ω and Φ y = α e −Ω in the limit Ω → ∞ keeping α e Ω fixed. 7 Notice that this choice enjoys invariance under scale transformations as a bulk quantum field theory, see [12].

Elasticity bounds
Thus far, we have discussed how to extract and understand the stress-strain curve for the black branes with solid scalar hair (3.3). One can go one step further and give an estimate on the where the stress-strain curve should terminate. As explained in [11], this can be done by studying the stability around the strained configuration. Generally speaking, this translates into a maximum strain max , beyond which the solutions have unstable perturbation modes which would render them unphysical.
The computation of the elasticity bound max for the solids which can be described by EFT methods was presented recently in [11]. We are now ready to perform the analogous computation for our holographic solids that, as emphasized above, model the special case when scale invariance is a manifest -the solid is in a nontrivial fixed point.
To this end, one should study the stability of the strained solutions under small perturbations. In the context of black brane solutions, this proceeds by the computation of (the dispersion relation of) the quasi-normal modes (QNMs) for the solutions (3.4). This effort is beyond the scope of the present work, however we can already gain insight analyzing the perturbations in the decoupling limit.
More precisely, we will exploit the two following approximations. First, rather than working out the QNMs, we will look at the local propagation of fields in the bulk. At this level, it is much easier to identify when the propagation speed c 2 (i) of some mode 'i' develops a wrong sign, c 2 (i) < 0. By the local propagation speeds here we mean the coefficients in the x, y-gradient terms in the bulk equations of motion for the perturbations, which are functions of u. Certainly, the AdS boundary conditions are such that the QNMs might turn out to be stable modes even if some speed c 2 (i) < 0 somewhere in the bulk. However, one expects that the threshold where the QNMs show an instability should be near the threshold where the first mode develops c 2 (i) < 0 somewhere in the bulk. In any case, using this bulk criterion is expected to place a conservative upper bound on max .
On the other hand, we are going to work in the decoupling limit where the mixing between scalar perturbations δΦ I with the metric modes is neglected. In this approximation, the equations for δΦ I are completely parallel to those in a 2 + 1 EFT in flat space -the exercise that was already done in [11]. On the configurations with finite strain (3.3) the phonon sound speeds are anisotropic, but one can still diagonalise the modes and thus obtain two characteristic speeds in the (x, y) plane that depend on the strain magnitude ε, and on the angle ϕ of propagation of the sound wave with respect to one of the the principal axes (eigendirections) of the strain tensor. The subscripts ± refer to the smallest and the largest of the two speeds. (In our holographic model, the speed in the holographic direction u is always 1.) The characteristic speeds are local, in the sense that they are defined at every slice u =const, but their structure in terms of the potential W is identical to the one obtained in [11]. Moreover, for monomial potentials like (3.18), it is easy to see that the u−dependence disappears from characteristic speeds. A conservative implementation of a stability criterium against gradient instabilities, then, is that (for all angles ϕ). In fact, for the model (3.18), the local speeds in the (x, y) plane depend on a, b in the very same way as how they depended on A, B in the first benchmark model of [11]. The condition (3.31) reads exactly the same, however the physical role of the parameters is different. This leads to a maximum strain deformation max which reduces to a certain function of a, b for the model (3.18). In the following, we discuss the obtained results, by referring to physical parameters like the scaling exponents ν S 1, 2 (Eqs. 3.23 and 3.24) instead of a, b.
Note that both ν S 1,2 scaling exponents can be bigger and lower than 1, meaning that both softening and stiffening can appear 8 . Nevertheless the second scaling is always smaller than the first ν S 2 < 2a because of consistency requirements (3.19). Aside from the conditions shown in (3.19) that restrict the values of a and b, we also find consistency conditions on the maximum shear strain that can be applied to the system, as mentioned above. From Fig. 4, we can already see that ν S 2 < ν S 1 , as expected. The upper bound to the maximum deformation max comes from two different consistency conditions depending on the value of ν S 2 : for the region ν S 2 > 3 we will first encounter ghosts (a change of sign in the kinetic term) at a finite value of shear strain, which will determine the value of ε max , while in the region ν S 2 < 3 this value will be determined by a gradient instability (a change of sign in the speed c 2 − ). There is a region with very large max in between these two sectors, i.e. around ν S 2 ∼ 3, where it grows asymptotically. We can find the specific relation between max and the shear scalings close to this area Surprisingly, the expression for the maximum deformation is the same either we are in above or below the value ν S 2 = 3, and it is independent of the value of ν S 1 . The subluminality constraint is also included Fig. 4, in the following way. The local speeds in the bulk (3.30) are allowed to be possibly superluminal, as the speeds in (3.30) do not give directly any physical phonon speeds. A proper limit would require the computation of the QNM dispersion relations at finite strain, which is outside the scope of this work. As a first step, though, we include the known results for the QNMs at vanishing strain [38,39], which certainly places a bound and one expects it gives an idea of the kind of bound one should obtain. We thus cut out the white region, which is where the physical longitudinal phonon QNM would be superluminal at ε = 0. By continuity, the true subluminality bound should make ε max vanish smoothly close to the edge of the white areas in Fig. 4. This has the effect to decrease ε max near the white edge and it should also render the actual plot less blue in the region ν B 1 ∼ ν B 1 ∼ 6, for instance. However, given the different slopes of the 'blue ray' and of the edge of the white area, one expects that the high elasticity region (the blue ray) persists mostly.

Bulk deformations
We can as well consider the bulk response beyond the linear approximation. In this section, we will neglect the finite temperature corrections; as a consequence the results presented here are robust at small temperatures T /m 1 but they will probably get finite T corrections elsewhere. In order to compute the non linear bulk response and avoid mixing the different deformations we set the shear strain ε = 0. In that case, the equilibrium configuration (3.7) is set by Φ I = x I ; this means that a bulk deformation κ = ∂ · φ corresponds to κ = 2(α − 1). Now, in analogy with linear response, we can compute the response in the pressure T xx = T tt /2, i.e. the longitudinal stress, with respect to the longitudinal strain in a full non linear form. We define the bulk stress as σ L = T xx (κ) − T eq xx . In more details for our model (3.18) we obtain: , (3.33) where with u h,κ=0 we mean the value of the BH horizon in absence of any bulk strain, i.e. α = 1. The sign of κ ≡ ∂ · π can be positive, i.e. a compression, or negative, i.e. an expansion.
The results for various potentials are shown in fig.5. At large bulk strain κ 1, we find a universal scaling which is a consequence of conformal invariance and it can be generalized to κ D , with D the number of spacetime boundary dimensions. The scaling can be immediately obtained analytically just realizing that in the limit κ 1 the radius of the horizon scales like u H ∝ 1/κ 9 . Furthermore, this result is in perfect agreement with the EFT computations [11]. Similarly to what happens in shear deformation, at large temperatures, there is an in between scaling that goes as α 2b = 1 + κ 2 2b , as can be seen in fig.5. Notice that there are again two types of deformation: isothermal (T = const.) and adiabatic (s = const.). In this last case, the scaling is always σ L ∝ κ 3 , while the intermediate scaling only shows up for isothermal deformations. Therefore, we can define two different bulk scalings, one that appears at intermediate strains and finite temperature and another that appears at very large strain or low temperatures, just like what happens for the shear deformation. These two bulk scalings are where the first scaling is always equal or larger than the second: So far, we have focused on bulk deformations with κ > 0 but we have not discussed negative values of κ, i.e. expansion. In this case, in order to produce a full expansion we need to go to the limit where α → 0 which corresponds to κ = −2. We can see from (3.33) that for low temperatures the stress goes as σ L ∼ (α 3 − 1) = (( κ 2 + 1) 3 − 1) and for temperatures high enough the scaling would be 2b instead of cubic. However, we do not get to see these non-linear scalings due to the short range of values α and κ have during the expansion. A couple of examples at different temperatures are shown in the right panel of Fig. 5, where we just see a linear scaling that gets saturated at some finite value. There, we also find that the amount of stress needed to produce a full compression is finite. We are not aware of any experimental signature of scale invariance in the linear and non-linear elastic response. Our results suggest that scale invariance should play an important role, providing universal scaling which can be possible tested at quantum critical points or within the quantum critical region.

Solids with spontaneously broken scale invariance
In this section we present the analysis of the nonlinear response, for models of solids that realize scale invariance (SI) as a spontaneously broken symmetry. This case can be treated using EFT methods in [11]. Linear elasticity of this systems has been recently discussed in [12]. The goal now is to extend the analysis to the nonlinear regime properties and to compare it with the manifest SI case studied in the previous section with AdS/CFT techniques.

Nonlinear response from EFT methods
In order to study the nonlinear response of a solid with spontaneously broken scale invariance we are going to employ EFT methods. We have already presented an EFT of a solid in (2.12), which has been studied in the past in [8][9][10][11], but now we want to consider only Lagrangians that are invariant under scale transformations. Therefore, we are going to demand that the Lagrangian that we introduced in (2.12) is invariant under with some 'weight' ∆ for the fields φ I that will depend on the potential, as can be seen in [12]. We find that the most general potential that we can use must be of the form of which is also identified as the equation of state parameter, ω = p/ρ. Notice that, although the potentials (4.2) and (3.1) look alike, their relation is not trivial as the framework where they are used is completely different. The comparison between the two should be done through physical properties of the system they describe such as the non-linear scalings. Moreover, this potential is not describing a system with conformal invariance unless we take the particular case where ω = 1/2.

General results
The non-linear bulk scaling is determined by ω and is completely independent of the function f (x). This scaling is given by Moreover, the bulk strain κ is unconstrained, as opposed to the shear strain ε. Looking at Fig. 6, we can see that the case of a monomial potential the scaling is restricted to satisfy 2 < ν B < 4. For a more general f (x), we can find that the limit is still the same, i.e. we cannot have ω > 1 without gradient instabilities or superluminal modes. Thus, neither the non-linear bulk scaling nor the constraints of it depend on the shape of f (x).

Benchmark models
The results related to the shear response are potential dependent, so we will need to specify what type of potential we want to study. In particular, we will take the two cases considered in [11], which are: 1. The first possibility we will study is the case where f (x) is a monomial, i.e.
with ω as a free parameter. This simple potential realizes a power-law scaling in the stress for large deformations. This potentials will display a power-law behaviour in the stress-strain relation for the shear channel, which will be determined by ν S , i.e. we will have that at large deformations σ S ∼ ν S .
2. The second possibility we are going to consider is taking both f (x) = 1 + v 2 x ν S /2 (4.5) and the limit ω → 0. The advantage of this potential is that the speeds of the phonon modes will be realistic (i.e. much smaller than the speed of light) as long as v 2 1, whereas the monomial potential has relativistic modes for generic values of n and ω. The shear response will only come from the x-dependent term, so the power-scaling of the shear response will be determined by ν S .
The discussion about the shear scaling ν S of the spontaneous broken scale invariance (SBSI) case is of course sensitive to the form of the function f (x). The simplest nontrivial example of potential one can think of is a monomial, ie. f (x) = x ν S /2 , which is shown in Fig.6. In this particular scenario both the bulk and shear scalings are constrained, in particular we find that 2 < ν B < 4 and 0 < ν S < 2. As with the solids with manifest scale invariance, here we can find one region where the maximum strain sustained by the material is significantly large. The most elastic region is found close to ν S ∼ 2 and we find that for this scaling ε max reads as The problem with this kind of potentials is that the speeds of the phonons are excessively large. Typically one would expect that the speeds of these modes are not bigger than ∼ 10 −4 times the light speed in 'earthly' materials (needless to say, this concern does not affect relativistic solids such as neutron star interiors). This constrains very much the possible scalings one can realize in a realistic scenario, specifically we would be restricted to ν S and ν B − 2 not bigger than 10 −8 .
If we do not restrict ourselves to the most simple case, we can describe a solid with slower phonon modes. For this, we are going to take ω → 0 and f (x) = 1 + v 2 x ν S /2 with v 2 1 as proposed in [11]. This form of potential ensures that the speed of the phonon modes are small at least for small shear strain. The shear scaling is forced to be positive ν S > 0 and the maximum shear scaling is thus ν max S 1 for v 2 1. The region with max large is found around ν S 2 and also in the region where ν S > 2 and v 2 1, specifically where we have taken the limit ω → 0 and the last max is valid for ν S > 2 not necessarily close to ν S = 2 as long as v 2 1, as can be seen in Fig.6. Figure 7. Summary of the main differences between the nonlinear elasticity bounds found for materials with spontaneously broken (left) or manifest (right) scale invariance. Both materials present power-law stress-strain relations σ ∼ ε ν (at ε 1) with exponents ν B and ν S for pure-bulk and pure-shear deformations. Both the values of the maximum allowed deformation ε max as well as the range of values of the allowed exponents differ substantially.

Comparison
Let us now compare the nonlinear response obtained in the manifest/spontaneously broken SI cases. In order to make the comparison as meaningful as possible, we will fix the physical observable -the stress-strain curve to be of the type σ ∼ ν for both bulk and shear deformations with exponents ν S, B -and compare the elasticity bounds in the manifest or spontaneously broken case. Our holographic examples present two powerlaw regimes, but we will restrict attention to the first ('intermediate') one labeled with exponents ν S, B 1 because the second one is only realized at extremely large deformations (at finite temperature). 10 The results shown in Sections 3 and 4 make manifest that the universal elasticity bounds that can be obtained from low-energy effective methods differ significantly depending on how scale invariance is realized (as a manifest or a spontaneously broken symmetry). The main differences are basically summarized in Fig. 7, where we repeat here the relevant plots in Figs.4 and 6 but in a comparable scale, for the sake of comparison.
It is clear from Fig. 7, the elasticity bounds present substantial differences, both in the range of allowed values in the ν S −ν B plane where finite deformations can be reached (ε max ∼ 1), as well as in the region where the deformations can be large (ε max 1). In fact, if we restrict the CFT analysis to the models that exhibit massless phonons (which requires ν B 1 > 5) then the areas in the ν S − ν B plane are almost disconnected, with the spontaneously broken case covering lower values of ν B -globally in the 2 < ν B < 4 window.
Similarly, the shear exponent is constrained to 0 < ν S < 2 for the SB case while it is basically unbounded in the manifest case. Another major difference is the presence in the manifest case of a 'very elastic band' (the bluish area) in the region near 11 This is the region that allows for black rubber -like holographic duals. It is worth to emphasize two properties about these solutions: i) the band in parameter space start at ν B, S 1 6 which corresponds to a rather stiff behaviour; ii) these values of ν B, S 1 are far from the free scalar limit (a = b = 1, corresponding to ν B, S 1 = 2). Thus, black rubbers require the presence of scalars with non-canonical (non-linear) kinetic terms in the bulk.
For the EFTs of Sec. 4 instead, the hyper-elastic region collapses down from a strip to basically a 'dot'-like area. Let us remark that this last feature happens in the first EFT benchmark (4.4) but not in the second one (4.5) (which displays small sounds speeds and also an elastic band at ν S 1 ∼ ν B 1 ). However, we prefer to keep the comparison at this level because the benchmark (4.5) depends on an extra parameter.
As mentioned above, the reason to introduce the second benchmark model for in Sec. 4 is to be able to have realistic phonon speeds -much smaller than the speed of 10 Another advantage of focusing on the intermediate scaling is that they allow for a continuum range of bulk exponents ν B 1 whereas the asymptotic one is fixed by conformality to ν B 2 = 3. 11 Interestingly, real-world rubbers are well fitted by power-law stress-strain curves with exponents satisfying (5.1) [6,11]. light c. A fair question, then, is whether this is also possible in the AdS/CFT framework or we are forced to have phonon speeds of the order of c. This issue has been already addressed in [12], where it is claimed that a solid with manifest SI and small phonon speeds is achievable. In the holographic set-up the speeds of the phonons are obtained by finding the spectrum of the quasi-normal modes, which have the form for both transverse and longitudinal modes, as checked in [25,38]. The low energy dynamics of these solids with manifest SI are described by a low energy CFT (and thus an IR fixed point). It is conceivable that the Lorentz group that emerges in the IR has a different light-cone speed c e , and thus the space-time metric is ds 2 e = −c 2 e dt 2 +dx i dx i . This has a crucial impact in (5.2) if we compare the results between a Lorentz invariant theory with a light-cone speed c and one with c e , as we will obtain a rescaling in the dispersion relation as ω → (c/c e ) ω. Therefore, for c e c the phonon speeds would get suppressed as If this is the case, we could also have a "realistic" solid with smaller speeds and, in addition, the white areas -limited by superluminal constraints-in Figure 4 would be enlarged.

Conclusions
We have analyzed the nonlinear elastic response in materials with scale invariance from the low-energy perspective, using effective field theory and holographic methods. The advantage in these effective methods is that they are mainly based on how symmetries are realized and therefore they can help to understand the nonlinear behaviour 'universally', that is, independently of the microscopic details of the material. The constraints imposed by how symmetries are realized imply nontrivial relations among different low-energy observables, especially the ones that encode the non-linear structure of the theory -the interactions in the material. In order to illustrate the appearance of these constraints from the low energy theories (and their dependence on how symmetries are realized), we have focused on the example provided by the elasticity bounds: the maximum deformability that a material may withstand in a reversible form. These were discussed in [11] in the case for materials with no manifest scale invariance. In this work, we have obtained these bounds for scale invariant (SI) materials, which can be of two types: with manifest SI or with spontaneously broken SI. The latter case can be described using the EFT methods and so is a particular case of those discussed in [11]. For the manifest SI case, instead, we have used holographic models. In order to include the two key ingredients (elasticity and manifest scale invariance), we have studied the simple holographic models of 'massive gravity' type. These are well-defined effective field theories in (asymptotically) AdS spacetime and they allow for a straightforward interpretation as a scale invariant field theories.
The main result of this work is to show how the elastic response can be extended to the full non-linear regime by obtaining the full stress-strain curves. The procedure is straightforward, and it gives rise to a very rich phenomenology of nonlinear elasticity behaviours which can easily be extended to other holographic models.
To our knowledge, this is the first time that the full non-linear elastic response of AdS black brane geometries is presented by extracting the corresponding stress-strain curve. Previous studies have discussed the elastic response only in the linear approximation. Ref. [29] discussed the visco-elastic oscillatory (that is, non-static) response and [49] an out-of-equilibrium similar setup.
As in [11], in order to make progress, we have assumed models that, by assumption, display power-law stress-strain curves, σ(ε) ∼ ε ν with some constant exponent ν at large strains ε 1. We have constructed the stress-strain curve and also followed how the appearance of pathologies (gradient instability, ghosts, or superluminal propagation) appears as a function of ε. We have then obtained how the maximum strain ε max depends on the exponents ν S and ν B (for shear and bulk transformations respectively).
We highlight three aspects of our obtained results. First, we find that the elasticity bounds for the solids with manifest SI (the holographic models) differ substantially from the ones for the solids with spontaneously broken SI (the EFT models), as shown clearly in Fig. 7. Both the ranges of allowed exponents and the values of the deformability ε max disagree by O(1) factors. Our interpretation of this discrepancy is that it is physical and due to the fact that SI is realized differently in the two cases, implying that the nonlinear constraints in the theories must differ.
The second aspect to illustrate is that our holographic models present the interesting peculiarity that they exhibit actually two different regimes of power-law stress-strain curves, σ(ε) ∼ ε ν , with different exponents at moderately large ε (say, 1 < ε < 10) and asymptotically large ε. This is due to the fact that these models contain a UV anisotropic Lifshitz fixed point. It is unclear whether this feature is only specific to the present model or whether it should be more generic in (very) elastic solids with manifest SI. Neither it is clear how relevant this feature is for realistic materials as this second scaling regime only appears for extremely large deformations, ε 10, before which the solutions already show some instability.
Finally, we stress that the holographic models that exhibit highest elasticity share features surprisingly similar to familiar real-world elastomers. Indeed, the models that allow black brane solutions which are stable under largest deformations (largest ε max ) turn out to have power-law exponents ν in the stress-strain curves σ ∼ ε ν which are: i) of order a few; ii) similar for pure-shear and pure-bulk ν S = ν B (corresponding to the blue stripe in Fig 7.b). Intriguingly, this is what happens for natural rubber and other elastomers [6] -and it motivates us to call these solutions black rubber.
Very importantly we would like to comment on the possible applications of this framework to realistic systems.
(I) The electric transport properties of quantum critical materials have been subject of lot of recent efforts especially in connection to the "anomalous" scalings found in Strange metals. More recently the question whether also phonons and elastic properties can display surprising and interesting features in quantum critical materials has emerged. More specifically there are preliminary indications that phonons in quantum critical systems can exhibit glassy or viscoelastic features [50]. Moreover the role of these viscoelastic properties has been discussed in connection to the possible implications on the onset of (high-Tc) superconductivity [51]. The framework just presented can provide a useful methods and possible observable predictions in this directions. The utility of these holographic models has been already proven in understanding the glassy features of amorphous solids and ordered crystals in [52][53][54].
(II) Let us also mention the growing interest related to the non-linear mechanical characterization of critical materials (e.g. High-T c superconductors) due to their technological applications [55,56]. Given the absence of robust computational methods, one may not rule out that the holographic methods such as presented here may provide useful insights.
(III) Following [29], a full non-linear characterization of the mechanical response of the holographic homogeneous models [21] considered in this work could definitely shed light on their physical nature, which is still open [23,30].
Finally, we remark a perhaps more technical point that we think is interesting from the gravitational point of view. The computation of the mechanical response to non-linear shear deformations translates in the AdS/CFT dictionary into finding exact black hole solutions with finite shear deformation (and shear stress) in the transverse directions. This implies that these solutions have a non-trivial, nonlinear spin-2 mode, and we devoted some effort to derive them in Sec. 3. These solutions are similar to the nonlinear gravitational waves and 'pp-wave' solutions in General Relativity in that they also have a nonlinear spin-2 (transverse-traceless) mode However, in our models the solutions are static, which is a direct manifestation of the massive character of the metric in these models. Moreover, the spin-2 mode 'sticks out' of a black brane horizon. In this sense, then, these solutions can be thought of as branes with spin-2 hair. We are unaware of solutions of this kind in the literature, but we find it remarkable that they are exist, and are tied to the nonlinear elastic response.
where g uu = u 2 f (u). The extrinsic curvature can be defined as usual θ µν = 1 2 (∇ µ n ν + ∇ ν n µ ) , (A. 4) and it reads θ = u f (u) − f (u) (6 + uχ (u)) 2 f (u) . (A.5) Close to the boundary, we can expand the function h(u) as in (3.13) and f (u) = 1 − M (u)u 3 (where M (0) would correspond to the energy density of the system) and find that the off-diagonal component of the stress-energy tensor is This result is a direct manifestation of the presence of a strain deformation in our background and it will encode the corresponding response, i.e. the shear component of the stress. Interestingly one can also notice that Using the standard holographic renormalization techniques [58] we can also identify µν is the sub-leading term of the induced metric expressed in Fefferman-Graham coordinates. As a first step we have to rewrite our ansatz in the FG form using the coordinate transformation dz 2 z 2 = du 2 u 2 f (u) (A.9) where z will now be the holographic (FG) coordinate. Again, using that the asymptotic behaviour of f (u) is 1 − M (u)u 3 we can find that for small z we have u = z − M (z)z 4 6 . Now we can already look at our metric and derive the stress-energy tensor. For instance, for the off-diagonal term, we are interested in, we have g xy (z) = 1 u(z) 2 sinh(h(u(z))) = (A.10) = 1 z 2 cosh(C 0 ) + (3 C 3 cosh (C 0 ) + M (0) sinh(C 0 )) z 3 , (A. 11) where higher orders in z have been suppressed. We can identify T xy easily in this expression and see that the result is the same we found before in A.6. This result give us a robust definition of the non-linear stress in our system which can be indeed identified as: σ = 1 2 (3 C 3 cosh (C 0 ) + M (0) sinh(C 0 )) (A.12) Since we impose the boundary condition C 0 = 0 the stress, we use in all our computations, is simply defined by σ = 3 2 C 3 .

B Three-phonon interaction terms in Solid EFTs
To further illustrate the predictive power of the low energy methods, we show here another nontrivial powerful statement that follows from the EFT construction which, as explained above, applies when translations (and possibly scale invariance) are broken spontaneously.
The nontrivial statement contained in the EFTs (2.12) is simply that once the stress-strain relations are known then the full Lagrangian is fixed. Let us illustrate now how this impacts for instance in the determination of the cubic phonon interactions. Assuming that the stress-strain relations both for bulk and shear deformations are known, then one can reconstruct the full form of the function V (X, Z) [11], up to an irrelevant additive constant.
Then, the cubic phonon interactions (around the homogeneous, isotropic equilibrium configuration φ I = x I ) are obtained by expanding our Lagrangian around the background solution, i.e. φ I = x I + π I . At third order in π I we obtain V (X, Z) ⇒ C 1 (∂ i π i L ) 3 + C 2 (∂ i π i L ) (π j ) 2 + C 3 (∂ i π i L ) (∂ j π k ) 2 + C 4π iπj ∂ i π j + C 5 (∂ i π i L ) (∂ j π k ) (∂ k π j ) , (B.1) where π I = π I L + π I T . The terms we find are In these expressions, the X, Z-derivatives of V are evaluated on the undeformed configuration. Moreover, by obtaining the relation between V (X, Z) and the stressstrain curve (such as e.g. Eq. (2.17)) one can relate all these V (X, Z) derivatives to derivatives of the stress strain curve at the origin, σ (0), σ (0), etc, which are measurable quantities. For instance, the bulk modulus is K = 4V ZZ + 2V Z + 4V XZ + V XX , G = V X and + p = V X + 2V Z .
This illustrates that the realization of symmetries implies nontrivial relations between distinct low energy observables. In this example, the strength of the phonon cubic interactions are determined by the shape of the stress-strain curves.
We can write the 5 C's as a function of these quantities and we would just need two independent new parameters We can compare our results with Ref. [7]. There he concludes that there are 3 independent new parameters, but we think the difference comes from the fact that he is working in 3 space-dimensions instead of 2. Moreover, in Ref. [7] there are 6 independent operators in the cubic expansion. It is trivial to check that in two dimensions the extra operator can be expressed as a function of the others (∂ i π j ) (∂ i π k ) (∂ j π k ) = (∂ i π i ) (∂ j π k ) 2 + (∂ i π i ) 2 (∂ j π k )(∂ k π j ) − (∂ i π i ) 2 . (B.14) In the case of scale invariance these terms simplify considerably. Let us then take V (X, Z) = Z Therefore cubic interactions are all fixed by these background or linear elasticity quantities. This is true both for general SI as well as conformal solid limit (which is just a particular value of w).