{'fraction_non_alphanumeric': 0.032297591397062855, 'fraction_numerical': 0.009363741107189232, 'mean_word_length': 4.682810226564124, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sum rules for the quasi-static and visco-elastic response of disordered solids at zero temperature

We study exact results concerning the non-affine displacement fields observed by Tanguy et al [Europhys. Lett. {\bf 57}, 423 (2002), Phys. Rev. B {\bf 66}, 174205 (2002)] and their contributions to elasticity. A normal mode analysis permits us to estimate the dominant contributions to the non-affine corrections to elasticity, and relate these corrections to the correlator of a fluctuating force field. We extend this analysis to the visco-elastic dynamical response of the system.

simplifying assumptions-either dealing with an infinite system or valid at zero stress. These simplifications make it difficult to ascertain the domain of validity of various formulae or symmetries. In response to these issues, we wish to attack this problem from two opposite angles: provide an even more systematic and general treatment than before, yet relate this formalism to numerical observations similar to Tanguy et al . [10,11] We thus hope that our treatment could serve, at least, as an introduction to the subject.
From the formalism, we wish to extract information about the contribution of various scales to the non-affine corrections to elasticity: the question we have at heart is whether these corrections originate from small or large scales, or involve a broad distribution of scales as suggested by the observation of vortex-like patterns. [9] This question is directly related to the existence of a continuous limit for elastic properties of amorphous structures. [10,11] To address these questions, we perform a normal mode decomposition of the non-affine displacement field: it permits quantifying the contribution of every frequency shell to the non-affine corrections to elastic constants. In the large size limit, these contributions seem to be self-averaging quantities (in the sense that an ensemble average over subsystems will produce results which are equivalent to a single large system). The existence of this selfaveraging property leads to an expression for the elastic constants in the continuum limit which resembles the sum rules of liquid theory.
Finally, our attention was attracted by several related issues in the recent literature.
Studies of sound propagation and attenuation in granular materials [17,18,19,20] or of the visco-elastic response of dense emulsions and foams, [7,8,21,22,23] emphasize the need for a deep theoretical understanding of the visco-elastic response of disordered solids. Related experiments by McKenna and coworkers indicated that features of the visco-elastic response of amorphous materials are related to measurable changes in their elastic constants. [24] We thus complement the study of static response by a study of dynamical response, and derive a formal expression for the visco-elastic moduli. We shall establish that the relaxation spectrum is directly and simply related to a correlator emerging from the normal mode analysis.
The present work is meant to present the general framework of our analysis, which will be the basis of future numerical studies. Although, here, we will use numerical simulations to illustrate analytical developments, the main core of our numerical study will be presented in a dedicated article. [25] 1

. NON-AFFINE DISPLACEMENT FIELDS
In this work we consider the mechanical response of an amorphous solid quenched at zero temperature. Our formalism permits dealing explicitly with finite size systems: it rests on the idea that, during a quench at zero temperature, any finite size system relaxes toward one of many local minima in the potential energy landscape. [26] Being at zero temperature, the system is then prescribed to lie at this minimum at all times. Small external perturbations are then expected to induce continuous changes in the local minimum. Large external perturbations may induce the vanishing of the local minimum occupied by the system: this vanishing occurs when the basin of attraction of this minimum reduces to a single point, that is when the minimum collides with at least one saddle point. [27,28,29,30] The difference between small, continuous changes of the local configuration and catastrophic events is exemplified figure 1 where the response of an amorphous system is shown as a function of shear. The parameter γ measures the total shear deformation from a quenched state. In this picture, continuous segments are associated with shear-induced changes of local minima and discontinuities to their vanishing. After a minimum has disappeared, the system, coupled to a zero temperature thermal bath, relaxes in search of a new minimum in the potential energy landscape. [30,31] The separation of continuous changes of the local minimum and catastrophic events is, of course, not limited to shear deformation. In a more general context, suppose that we denote by γ the amplitude of any external drive applied to the system. Suppose that the system lies at a minimum which, for γ = γ 0 , is far from any catastrophe. Let us denote δγ = γ − γ 0 the amplitude of a perturbation of the external drive around γ 0 . If we were to slowly increase the amplitude of the perturbation from 0, the system would smoothly follow a trajectory r i (γ) in configuration space for all amplitudes δγ < ǫ such that the basin of the local minimum remains non-vanishing. If the Hessian is non-degenerate, this trajectory r i (δγ) is unique. If we now stop the perturbation at some δγ max < ǫ then revert the drive down to δγ = 0, the system would simply follow the trajectory r i (γ) from γ = γ 0 + δγ max to γ 0 , backward. In this sense we will say that the continuous segments are microscopically reversible.
When the system is driven quasi-statically along these continuous segments, no energy is dissipated. The reason is that at each point along these segments, the system is at mechanical equilibrium: the force applied to any particle is exactly zero and zero forces do no work.
The quasi-static response corresponds to the thermodynamically reversible, elastic, part of the mechanical response. We will see, however, that energy is dissipated when the system is driven at finite deformation rates along these continuous segments. This dissipation results from the fact that finite deformation rates induce non-zero forces which dissipate energy via the coupling with the zero temperature thermal bath.

Notation
We shall let underline and double-underlines respectively indicate vectors and tensors referred to a fixed Cartesian system (x 1 , x 2 , x 3 ). We shall also use the convention that Greek indices refer to Cartesian components of vectors or tensors, while Roman indices refer to the particle numbers. Bold type denotes fields which are defined on every particle in the material: Dots and double dots indicate matrix products and summation convention is always applied on repeated (Greek and Latin) indices. By convention, we also write: The superscript T indicates the transpose of a matrix, −1 its inverse, and −T the inverse of its transpose. The identity matrix is denoted 1, and its components δ αβ .

Displacement fields
In this work, utilizing the formalism underlying the Andersen-Parrinello-Rahman theory [32,33,34,35,36,37] we shall focus on the situation where a system of particles is contained in a periodic simulation cell. The formalism we describe can easily be adapted to study the deformation of a material confined between walls: to do so, it is sufficient to embed the whole system-confined particles plus walls-in the cell and mandate that the particles constituting the walls affinely follow its deformation.
The shape of the simulation cell (which is, by construction, a parallelepiped) is represented by the set of d = 2 or 3 Bravais vectors: h = (a, b) or h = (a, b, c); its volume is V = det(h) = |h|. We consider a system of N particles, with positions r = {r i } in real space.
The interaction potential U(r, h) depends on the positions of the particles but also on the shape of the simulation cell which enforces boundary conditions.
"Macroscopic" deformations of the sample are performed by changing the Bravais vectors. Since we are concerned with variations of the local minimum around some reference configuration, we will often use a reference configurationh of the cell and compare it with a current value h. Following Ray and Rahman, [35,36,37] we introduce a transformation of particle coordinates which maps any vector r onto a cubic reference cell:

Equilibrium trajectories
By definition, elastic constants are second order derivatives of the energy with respect to strain. Since strain is characterized by second order tensors, the elastic constants are fourth order tensors. Before presenting the details of the tensorial formalism, however, it seemed pedagogically sounder to us to first consider the simple situation where the shape of the cell can be parameterized by a single degree of freedom.
Suppose then that we prescribe the tensor h(γ) as a function of a scalar parameter γ. For varying γ, so long as the local minimum does not vanish, the system follows a continuous trajectory in configuration space as illustrated on figure 1. Given a reference cellh = h(γ 0 ) at γ 0 , the energy functional can be written either as a function of r and γ: U(r, γ); or as a function ofr and γ:Ů(r, γ) ≡ U(F (γ).r, γ). We introduce the notationŮ to emphasize that-contrarily to U-this function is defined after a choice of reference cell with Bravais matrixh. Changing γ for fixed {r i } corresponds to performing an affine strain of the whole system-the particles and the boundary.
When particles are constrained to follow deformation-induced changes of a local minimum, their real-space positions {r i (γ)} and the corresponding non-affine displacements {r i (γ)} are now one-parameter functions of γ. An equation of motion for the non-affine displacement fields derives from a straightforward application of the implicit function theorem: The trajectory is specified by the condition that the system is always at mechanical equilibrium: Note that the second equality is a property of the point derivatives of any observable of the and is derived in appendix A. Differentiating the condition ∂Ů ∂r i = 0 once with respect to γ leads to: where the symbol D is introduced to indicate derivatives which are taken under the constraints of mechanical equilibrium. This equation is formally valid for all γ = γ 0 , even though it will primarily be used in the limit γ → γ 0 (i.e. h →h).
Using equation (2), we see that in the limit γ → γ 0 , equation (3) involves the Hessian: and the field of virtual forces, In semi-condensed notation, equation (3) reads: In order to solve equation (6) we need to take care to eliminate the zero modes of the Hessian. Generically there are d zero modes for a system in d dimension: those correspond to translation invariance; their existence indicates that a solution to equation (6) can only The zero modes can thus be eliminated in the standard way by subtracting off the projection onto uniform translations from any particular solution. Up to this invariance, the solution of (6) reads: The numerical practice which consist in, first affinely deforming the simulation cell, then letting the system relax toward the displaced minimum, provides a direct physical interpretation for this equation. Ξ is the field of forces which would result from an elementary affine deformation of all the particles. Equation (7) shows that the non-affine displacements are just the linear response of the system to these extra forces.
Examples for the fields Ξ and Dr/Dγ in simple shear and pure compression are shown figure 2 and 3. (The details of the numerical simulation will be given in section 1.5.) We observe that the short-range randomness of the vector Ξ contrasts with the large vortexlike structures displayed by the non-affine "velocity" fields Dr/Dγ. We will later use the apparent disorder of the field Ξ to construct a statistical treatment of the contributions of these non-affine displacement to the elastic constants. For now, let us move on to study how these non-affine fields enter microscopic equations for the elastic constants.
Upon deformation, the first and second derivative of the energy U(γ) =Ů (r(γ), γ) are related to the components of the stress and the elastic constant respectively. We will later enter a full tensorial derivation of the equations to specify this relation. For now, we can simply differentiate U(γ) with respect to γ. At first order, we obtain a projection of the stress: where the second equality holds because of mechanical equilibrium. The partial derivative which appears at the rightmost side of these equalities corresponds to the projection of stress which would be inferred from the assumption that all the particles undergo affine displacements. Since the total derivative of the energy as a function of γ (while enforcing the condition of mechanical equilibrium) is identical to the partial derivative, it means that non-affine displacements do not contribute to the first derivative of energy with respect to γ, i.e. the projected stress.
The situation, however, is different when the second derivative of the energy is considered.
Indeed, we have: The first term in the right hand side of this equation, ∂ 2Ů ∂γ 2 , corresponds to energy changes associated with strictly affine displacements of the particles. It is the Born approximation for the second derivative of the energy with respect to γ. The last term introduces a correction to the Born approximation which results from the existence of non-affine displacements. Since H −1 is positive definite, the corrective term is negative: non-affine displacements reduce the second order derivative of the energy from its Born approximation. This is another statement of the idea that non-affine displacement are associated with further minimization of the energy of the system from a tentative homogeneous deformation. This property holds for any mode of deformation h(γ). It means, for example, that the non-affine corrections to the shear and compression moduli must be negative.

Tensorial forms
In the preceding section, we considered the non-affine displacements in response to some prescribed mode of deformation parameterized by a scalar parameter, γ. We now proceed to study the tensorial form of the above equations: derivatives of the energy functional with respect to the components of the strain tensor will provide analytical expressions for stress and elastic constants.
The calculations can be simplified by noting that the number of independent components of the strain tensor is reduced by symmetry. In particular, a useful symmetry property holds under the general assumption that [1,2,40,41,42] the interaction potential can be written as a function of the full set of distances between pairs of interacting particles: U(r, h) = U {r ij = |r ij |} , where the index ij runs overs all pairs of particles, and for each pair, r ij = r j − r i modulo the periodic boundary conditions. (Note that for a periodic simulation cell, the distance between two points depends on their positions relative to the boundaries of the cell: the way distances are calculated thus depends on the shape of the cell, hence, on h.) We stress that the ability to write the energy functional in this way is not restricted to potentials involving only pair interactions: it also holds for e.g. three body interactions which depends on bond angles, as used for silicon [43], or embedded atom potentials [44,45], as used for metals. This formalism does not apply, however, to situations when rotational degrees of freedom must be taken into account as, for instance, in Cosserat approaches to granular materials [46,47,48,49] or to model anisotropic objects like nematics. [50] We expect a similar formal treatment to be possible in these situations, but it will require using a slightly more general formalism that we do not wish to address here.
Let us consider a reference configurationh and a current configuration h. Suppose that r andr are the difference between the positions of two particles in both these systems of coordinates. They are related by r = F .r, whence: with the Green-Saint Venant strain tensor: [35,36,37] We want to write the energy functional after the change of variable where any config- Under the general assumption that the potential energy can be written as a function of the distances {r ij } only, equation (10) can be used to write: Where it appears that the energy of the system is a function of {r ij } and η only. [1,2,40,41,42] The notationŮ is introduced here to emphasize that this functional of {r ij } and η is defined after a choice of reference cell with Bravais matrixh (whereas the functional U(r, h) does not depend on this choice). A choice ofh being made, however, the energy functional depends on the current cell coordinates only via the symmetric tensor η, but not via the whole tensor h. We recover the separation of coordinates: changing η for fixed {r i } corresponds to performing an affine strain of the whole system-the particles and the boundary; changing {r i } in the reference cell corresponds to performing non-affine displacements of the particles.
We are now in a position to derive equations of motion for non-affine displacement fields in a fully tensorial form. In any configuration, the force acting on particle i is: and mechanical equilibrium reads: The derivatives are taken here for fixed cell coordinates. The equation of motion for {r i }, is now obtained by differentiation of (13) with respect to the components of η: with the usual summation convention of repeated (Greek and Latin) indices. As above, the symbol D denotes partial derivatives with respect to some tensorial components while enforcing mechanical equilibrium (13). In the limit η → 0, equation (14) involves the Hessian: and the field of third order tensors, In semi-condensed notation, equation (14) reads: Provided standard caution in the inversion of the Hessian, we then get, up to translation modes: which is the non-affine displacement field in response to elementary deformation along η κχ .
To get further insight into the correction term in equation (24), and on the physical interpretation of Ξ, let us write: This equality comes after differentiation of equation (12) with respect to η and use of mechanical equilibrium (equation (13)). As in our pedagogical example of one-parameter strain, the field Ξ κχ , can be interpreted as the forces which would result from an elementary affine displacement of all particles in the strain direction η κχ . (More specifically, ∆η κχ Ξ κχ , with no summation on the greek indices, is the force resulting from a small affine transformation by ∆η κχ ; Ξ κχ is the tangent direction to the changes of forces upon affine transformations.) We see that Ξ i,κχ can be seen as the fluctuation of the force f i in response to an elementary strain. This interpretation emphasizes the random character of field Ξ κχ : local forces are, by definition zero at mechanical equilibrium; however, their variation under an elementary strain depends on the configuration of the particles with which particle i interacts; in particular, Ξ i,κχ should be zero if the conformation of particles surrounding i is symmetric-which is the case of a Bravais crystal. Ξ κχ is thus a measure of the local asymmetry of particle configurations.
Equation (18) further states that the non-affine displacement is nothing but the linear response of the particles to these extra forces.
Let us now derive microscopic equations for the stress and elastic constants. Because the total potential energy is a function of η andr only it is sufficient to characterize the elastic response of a material using the derivatives of the potential energy with respect to the components of η. We note, however, that the usual definition of the Cauchy stress and elastic stiffnesses involves derivatives with respect to the components of the deformation gradient tensor F . Derivatives with respect to the components of η provide the so-called thermodynamic tension and elastic constants (or also thermodynamic stiffnesses). The definitions of thermodynamic tension, Cauchy stress, elastic constants, elastic stiffnesses and their relation are given in appendix A.
By definition, the thermodynamic tension is the derivative of the energy functional with respect to the components of η:

DŮ Dη
Since η is a symmetric tensor, we can conclude without further examination that thermodynamic tension is also symmetric. Using mechanical equilibrium, we next have: As in our pedagogical example of one-parameter strain, we see here that the thermodynamic tension is equal to the partial derivative ofŮ (here, with respect to η). This partial derivative corresponds to changes in energy during a strictly affine displacement of the particles: we see that the existence of non-affine displacement fields does not appear in the expression for the thermodynamic tension. Expression (20) is valid for any finite strain: this will allow us to later differentiate this expression a second time and obtain elastic constants. Often, though, we will use it in the limit η → 0.
Note that this expression for the thermodynamic tension provides another interpretation of the field Ξ. We can write: from the definitions of t and Ξ. This expression is to compare with equation (19). We see that Ξ i,κχ can be seen as a fluctuation of stress (or tension) in response to an elementary displacement of particle i. This interpretation was noted in [14].
The elastic constants are second derivatives of the energy functional with respect to the Green-Saint Venant strain tensor: The derivative is here "total" since this is the second derivative of the energy following deformation-induce changes of a minimum. As we can permute the order of derivatives, elastic constants verify C αβκχ = C κχαβ , and since η is symmetric, To obtain an analytical expression for the elastic constants, it suffices to differentiate expression (20) once and take the limit η → 0 afterwards: Using equation (18) and the definition of Ξ αβ , it then comes: We recognize the first term in equation (24) as the Born approximation C Born αβκχ . The contraction of the inverse of the Hessian on components of Ξ αβ provides the correction terms.

Special case: Isotropic material in 2D
Here, we illustrate these ideas with numerical simulations of a two-dimensional bidisperse mixture of particles interacting through a shifted Lennard-Jones potential. [11] Particle sizes r S = r L sin π 10 /sin π 5 and a number ratio N L /N S = 1+ √ 5 4 are used to prevent crystallization.
The system is prepared via an initial quench from an infinite temperature state. Further details regarding the numerical protocols will be found in our forthcoming dedicated numerical study [25].
The interaction potential is pair-wise additive: U({r ij }) = ij V ij (r ij ). In this case, formal expressions for the fields Ξ i,κχ and the Born approximation for the elastic constants can be related directly to the derivatives of V ij . These expressions ( (71) and (72)) are derived in appendix B. We just recall here equation (72): where we have introduced the normalized vector between pairs of particles: n ij = r ij r ij , as well as the bond tensions and stiffnesses: The material is expected to be isotropic so that the elastic constants take only two independent values under permutations of indices: which define the Lamé constants, λ and µ. Using equation (25) we observe that the Born term is identical for the two Lamé constants λ and µ. Introducing the notation n ij = (cos θ ij , sin θ ij ), we find: To calculate the Lamé constants, it suffices to consider two modes of deformation and the associated non-affine fields. We thus consider Ξ s = Ξ xy = Ξ yx and Ξ c = (Ξ xx + Ξ yy )/2 which are associated with pure shear and pure compression respectively. Plots of Ξ s and Ξ c , and their corresponding displacement fields were given on figure 2 and 3 respectively.
We note that the correction term to the compression modulus appears to be small. This can be understood to arise from a simple property. Let us consider the case when the potential is pairwise additive, and is homogeneous in the sense that the force between each pair of particles is a homogeneous function of their distance: For a compressed Lennard-Jones system, if the total energy is dominated by the pairs of particles which are at close distance, the interaction is thus dominated by the repulsive power-law divergence of the potential. Since this power-law leads to forces which are homogeneous functions of the distances, this repulsive part of the Lennard-Jones potential leads to vanishing corrections to the Born approximation. As a result the overall amplitude of the field Ξ c is small, hence the small correction to K. In particular, this is consistent with the data from Tanguy et al [10,11] and explains why the corrections induce important changes in both λ and µ but smaller changes in their sum.

NORMAL MODE DECOMPOSITION
In expression (24), the correction to elastic constants involves a contraction of the Hessian on the fields Ξ αβ . Since the Hessian is positive definite-we are perturbing around a local minimum-it acts as a scalar product in expressions such as: Ξ αβ .H −1 .Ξ κχ . This form of the corrections to elasticity suggests that further insight may be gained from a normal mode decomposition of the fields Ξ κχ .
We denote the eigenvectors of the Hessian by Ψ p and the associated eigenvalues by λ p .
We introduce a particle mass to correctly scale eigenvalues and eigenfrequencies and write λ p = m ω 2 p . (In our numerical simulations, the mass is taken to be unity.) The vector Ξ κχ is written as: Using this decomposition and equation (18), the non-affine displacement fields read: Dr the non-affine contribution to elasticity reads: As we saw from figures 2 and 3, the fields Ξ display a very random character. The At high frequencies the scatter plots thin out because there are fewer and fewer eigenvectors to sample the density of states close to its upper cut-off.
On the basis of this observation, we make the assumption that the random fields, Ξ κχ , and in particular their projections on normal modes, Ξ p,κχ , are self-averaging quantities. We thus introduce the correlators on frequency shells: In this definition, the average is performed for all the projections of Ξ αβ and Ξ κχ on eigenvectors with eigenfrequency ω p ∈ [ω, ω + dω]. The assumption that Ξ p,κχ are self-averaging means that the quantities Γ αβκχ (ω) are expected to converge toward well-defined functions of ω in the thermodynamic limit: either when a large ensemble of systems of a given size is considered, or when one single system of a large size is considered.
From their definition, we observe that functions Γ αβκχ (ω) verify the same symmetries as the elastic constants, namely: In the thermodynamic limit, we can finally rewrite equation (24) as: where ρ is the density of states [70]. This equation is a sum rule: it relates elastic constants to the integral of a correlator between microscopic fields.
Like the elastic constants, Γ αβκχ are tensorial quantities. Therefore, if the material is isotropic-given the symmetry property (29)-these functions can take only two independent values corresponding to the two Lamé constants λ and µ: In this case, the two correlators, Γ λ and Γ µ provide a complete description of the non-affine corrections to the Lamé constants according to equation (30).
To illustrate this discussion, we return to our two-dimensional numerical study of two modes of deformation, simple shear and pure compression. This protocol suffices to measure the two Lamé constants of our system. The vector fields Ξ s and Ξ c grant direct access to the quantities Γ u = Ξ 2 p,u for u = s, c. Elementary algebra, similar to the calculations of section 1.5, shows that these quantities also verify: Γ µ = Γ s and Γ K = Γ λ + Γ µ = Γ c .
We report, in figure 5, the density of states ρ(ω), the correlators Γ s , Γ c , and the quantities C u (ω) = ρ(ω) Γ u (ω)/(mω 2 ) (with u = s, c) which measure the total contribution of the shell [ω, ω + dω] to the non-affine corrections. Three system sizes have been used in this plot: L = 20, 25, and 30, where L is the length of the square simulation cell in dimensionless units.
For each size, an ensemble of 20 system was used to provide statistical accuracy.
The functions Γ s = Γ µ and Γ c = Γ K measure the variance of the distributions of Ξ p,s and Ξ p,c respectively, which appear on the scatter plot of figure 4. We see that that these functions are roughly increasing as was suggested from the observation of the scatter plot.
Important fluctuations of Γ µ and Γ K at high frequencies result from the lack of statistical representation in the region where the density of states vanishes and where the scatter plots thin out. For their main part, however, the curves lie on top of one another for the three system sizes we have used: this is evidence that the continuous limit has already been reached. This observation is consistent with the numerics of Tanguy et al [10,11] since the values of the elastic constants they measured did not vary much for the different sizes they studied, from L ∼ 20 to L = 150. We believe that the good convergence of these measurements results from the decorrelation of the fluctuating forces Ξ i at very short distances: because the correlation length is short, the finite sizes of the systems considered here provide enough sampling of the distribution of Ξ i .
The density of states, and the Γ's exhibit rather complicated functional forms. In contrast, it is striking to us that for simple shear and pure compression C u (ω) seems rather featureless; particularly in comparison with both the density of states and the functions Γ u (ω). We have observed this for different numerical models and different densities. [25] The monotonous decay of C µ (ω) and C K (ω) seems to indicate that a simple physical mechanism-like a transfer of elasticity from low to high frequency phonons-governs the non-affine corrections. This reminds us of an analysis by Radjai and Roux of quasi-static deformation of granular materials [9] in which the authors suggested that the non-affine velocity field could exhibit statistical properties similar to turbulence. This sound like an appealing suggestion, although, at present, we have no evidence to further test the existence of turbulent-like energy transfers from large to small scales.

VISCO-ELASTICITY
In the above formalism, we assumed that deformation was very slowly forced upon a piece of material, so that the relaxation of the system toward an inherent structure was always faster than the changes induced by the external drive. This guaranteed that the system followed an energy minimum in configuration space. Our analysis, however, is essentially based on a second order expansion of the energy functional: it assumes small displacements, but may, in principle, apply to situations when the timescale of the external drive becomes relevant.
We consider here oscillatory perturbations of small amplitude, and assume that the system remains in the linear response regime around some energy minimum, far away from any catastrophic event. We expect this situation to be of relevance to questions such as e.g.
acoustic damping in granular materials [17,18,19,20] or viscoelasticity in dense emulsions and foams [7,8,21,22,23]. Near a minimum, but not exactly at mechanical equilibrium, the particles are subjected to non-zero forces: an oscillatory external forcing at finite rate thus works on the system of particles. Energy is injected into the system and is lost in dissipative mechanisms at the microscopic scale (friction between grains or viscous dissipation between bubbles).
We will follow the usual practice of characterizing energy transfers in the linear response regime using the components of the complex stress response, that is the storage and loss moduli, G ′ and G ′′ . [51] We insist that we focus here on features of G ′ and G ′′ that arise merely from the existence of non-affine displacement fields around a single minimum. In an experiment-think, for example, of a foam-it is unclear whether sufficiently small strain amplitudes can be achieved so as to stay within the basin of attraction of a single equilibrium configuration. Thus measured values of G ′ and G ′′ would also receive contributions from plasticity in addition to the microscopic visco-elasticity studied here: transitions between distinct inherent structures would become a relevant mechanism of energy dissipation [52,53]. It's only after a careful study of the contribution of linear response that we will be in a position to characterize the relative importance of these various dissipative mechanisms.
The present work is a preliminary step in this direction.

Damped equations of motion
In order to work at finite strain rates, we need to use Newton's equations coupled with the deformation of our simulation cell. In order to prevent the system from heating up, a microscopic mechanism must be specified for energy dissipation. Dissipative mechanisms vary from system to system. In granular packings, dissipation involves the deformation of asperities at contact between grains. In foams, dissipation involves fluid flow inside a liquid phase: it is modeled by Durian [54,55] by viscous terms involving velocity differences between neighboring bubbles. In numerical simulations of structural glasses and supercooled liquids, the coupling with an external bath is modeled either by a Nosé-Hoover [56,57,58] thermostat, by use of a Lagrange multiplier [59,60], or by a Berendsen thermostat [61,62].
It thus appears, after examination of these different systems, that energy dissipation often arises as some sort of viscous damping at the microscopic level. We will thus limit our present discussion to the situation where dissipation enters via a viscous drag force applied on individual particles. Below, we give some rule-of-thumb estimates for what the value of this viscosity might be in the case of dense emulsions.
In the SLLOD algorithm, the viscous term is generally taken to be proportional to the peculiar velocity, [60] that is the velocityṙ i of the particle minus its velocity corresponding to the affine flow u(r i ) =Ḟ .r i =Ḟ .F −1 .r i . Similarly, in foams or granular materials, the viscous damping does not depend on the velocity in real space but on the difference between the particles' velocity and their surrounding environment. We thus first write the equations of motion as: In term of {r i }, it follows that: We see here that a term mF .r i arises which depends on the position of the particle in the simulation cell. Another term of similar flavor would be νḞ .r i : it does not appear in equation (32) because it has already been eliminated by the above assumption that the viscous dissipation applies to the peculiar velocities. Such terms introduce a spatially dependent coupling between particular motion and the cell coordinates: they cannot be allowed to enter the equations of motion if we want to enforce translation invariance. This was originally noted by Andersen [32] and later by Ray and Rahman [36]. Within the formalism of Souza and Martins [63], an equation of motion with no spatially dependent terms can be derived systematically from a Lagrangian. It relies on the idea that the distance traveled by a particle should be calculated while removing a component of rigid rotation.
We do not wish to detail this formalism here, but simply recall that adding appropriate terms in the Lagrangian of the system (particles plus simulation cell) permits eliminating space-dependent terms from equation (32). Following these authors, we eliminate the term mF .r i and focus on translation invariant equations of motion for an athermal system in SLLOD dynamics. These equations read: We want to account for the dynamics of the system when it is submitted to small amounts of strain F − 1 ∼ ǫ << 1. To do so, we write a perturbative expansion of equation (33) around a known equilibrium state {r i }. This expansion is written in terms of the displace- which are also of order ǫ. At first order in ǫ, the perturbed equations of motion read: Using the definition of H and Ξ (equations (15) and (16)) it follows that: which is solved by: The perturbed stress is then obtained by writing at first order in ǫ: Using (37), the complex stress response in the frequency domain thus reads: It can thus be recast in the form: with: Taking the thermodynamic limit, we can furthermore introduce the functions Γ αβκχ (ω) defined in equation (28). Equation (38) can then be written as an integral in frequency domain: It is then an easy task to extract the real and imaginary part of the stress response. Observe that we recover the Born approximation and the true elastic constants in the high and low frequency limits respectively.

Overdamped systems and the relaxation spectrum
The response moduli are often studied in the 'overdamped' limit, when the kinetic energy remains small compared to the elastic energy. [51] In this case, the term m ω 2 is negligible before m ω 2 p and equation (39) can be given a simpler form. To do so, we introduce the timescales: and write equation (39) as: The complex modulus G αβκχ is thus naturally written as a sum of Maxwell elements. By definition, [51] the function H αβκχ is the relaxation spectrum of the system. It is the product of the density of states with the correlator Γ αβκχ . The specific forms of this relation associated with a given experimental geometry -pure shear or pure compression-derive trivially from these tensorial expressions. We plot the shear relaxation spectrum for our model system in figure 6, again for each of the three ensembles of systems with lengths 20,25, and 30, demonstrating the system size independence of the results.
Since the true microscopic mechanisms of dissipation in real structural glasses remain somewhat poorly understood, we focus here on the case of a dense emulsion or foam, where the system is overdamped by construction. To make a semi-quantitative comparison with experiments on such systems, we need to give a rough estimate for the values of the coefficient of viscosity, ν, and the overall stiffness scale,λ, in our model. A truly quantitative connection with experiment is somewhat beyond the scope of the present work and will be explored in our dedicated numerical study [25]. We should also remind the reader that the data presented in this work for illustrative purposes is taken from a compressed system with Lennard-Jones interactions... clearly not sufficient for a fully quantitative comparison. However, the hope is that there is a good deal of qualitative similarity in the structure of the Hessian matrix for these various classes of disordered systems, and we proceed to dimensionalizeλ and ν in order to make a semi-quantitative comparison with experiments on packings of emulsified droplets. Following Durian [55], we ought to choose a value for ν such that the drag force, F drag = νδv, between two droplets sliding past each other with relative velocity δv generates a stress in the film which is equal to the strain rate in the film times the effective viscosity, so we must have: ν = F drag /δv = η eff  figure 6, we expect the viscoelastic effects and the crossover to affine deformation to appear in the range of frequencies from around 1Hz up to about 1kHz in these systems.

CONCLUSION
In this article, we have studied molecular displacements associated with deformationinduced, continuous changes of a minimum in the potential energy landscape of an amorphous solid. For small amounts of deformation, the trajectory followed by the particles is smooth, and the tangent displacements-analogous to a velocity, with strain playing the role of an effective time-involve affine and non-affine fields. Both fields enter at first order in the equation of motion of the particles during quasi-static deformation.
The non-affine fields can be calculated by inverting the Hessian on a fluctuating force field Ξ, which is formally a third order tensor. For every component η αβ of the deformation, Ξ αβ is the field of forces resulting on every particle from the affine motion of its neighborhood.
The non-affine displacement fields and the corrections they induce on elastic constants, can be expressed solely in terms of the Hessian and of various tensorial components of the field Ξ αβ . We next observe that the field Ξ αβ is weakly correlated in space, so that it can essentially be considered as a random vector field. This randomness is confirmed by the scattered values for the projections of Ξ αβ on the eigenmodes of the Hessian.
The normal-mode analysis of the field Ξ αβ grants access to the contribution of each frequency shell to non-affine corrections to elasticity. Our numerical calculations indicate that the correlation function Γ αβκχ (ω) = Ξ p,αβ Ξ p,κχ ωp∈[ω,ω+δω] converges toward a smooth function in the thermodynamic limit: this means that the contribution of each frequency shell to non-affine corrections to elasticity is self-averaging. In our view, the existence of a well-defined limit for the elastic properties of amorphous structures, as observed by Tanguy et al , rest on this self-averaging property.
We moreover observed that the contribution of frequency shells seems to be a rather simple function of the frequency, as opposed to both Γ αβκχ (ω) and the density of states.
This simple form of corrections to elasticity suggests that an elementary mechanism of transfer of energy between frequency shells may be at work and determines corrections to elastic constants. This is reminiscent of the observation by Radjai and Roux of turbulent like features of non-affine displacement fields, [9] although we have not determined whether turbulent-like scalings arise in our non-affine fields. A better understanding of the mechanism underlying this transfer through frequency shells should in principle permit construction of approximations for the elastic constants of amorphous materials.
After studying elastic properties of an amorphous solid in response to quasi-static deformation, we considered the case when the deformation rate is finite. This required us to introduce both a molecular mass and a damping term to provide equations of motions which correspond to a molecular system in contact with a bath at low or zero temperature.
We have shown that the visco-elastic response of the solid can be written in the form of a sum of elementary damped oscillators. The frequency spectrum-that is the distribution of timescales of elementary vibration modes-is directly related to the function Γ αβκχ (ω).
It thus appears, that a broad spectrum of timescales arise at linear order solely from the structure of any particular energy minimum in a high dimensional configuration space. The system remains solid since it resides in the neighborhood of a local minimum in configuration space at all times: this is different from the usual idea that visco-elastic response is associated with transitions between minima in the energy landscape [52,53]. In a real system, this type of dissipation may occur simultaneously with other mechanisms: energy transport via phonons, dissipation resulting from anharmonic terms, or transitions between various energy minima (i.e. true plasticity). The estimation of these various contributions and their interplay seem to us an exciting direction for future research. It will require the construction of models, close to the experiments, where various dissipative mechanisms can be estimated quantitatively.
Finally, we stress that the corrections to elastic constants involve 1/ω 2 p factors as in equation (30): the convergence of the integral in equation (30) thus depends on the low frequency behavior of the functions Γ αβκχ and of the density of states. This becomes a problem when the system develops low frequency eigenmodes which have a non-zero scalar product with the field Ξ: a quick inspection of equation (27) indeed shows that such localized low frequency phonons would lead to diverging terms in the non-affine corrections to the Born approximation. We can identify two situations when this occurs: (i) A localized eigenvector with a vanishing frequency is involved whenever the local minimum occupied by the system reaches a catastrophe. [27,29] This property enabled us to obtain universal scalings for the elastic moduli close to a shear induced catastrophe, [30] and observe the divergence of the non-affine corrections at these points. (ii) Eigenmodes seem to drift toward the low frequency part of the spectrum when the system approaches the unjamming point of Liu and Nagel. [64,65,66,67] A divergence of the non-affine corrections to elasticity could thus arise around this point. In this case, we expect the divergence of the non-affine corrections to elasticity to control the unjamming transition around the J point. We saw the important role played by the Green-Saint Venant strain tensor η: it accounts for the mapping of distances after an affine transformation (see equation (10)). This property permits writing the strain-dependence of the energy functional via η only, whenever the energy is a function of the set of distances {r ij } between the particles. [41] By definition, the thermodynamic tension [12,40,68] is conjugate to the Green-Saint Venant strain tensor: WhereV is the volume of the reference cell. The thermodynamic tensiont is defined after a choice of a reference configuration: it is a Lagrangian derivative. It identifies, in the continuum limit, with the second Piola-Kirchhoff stress tensor. [38,39,69] This object is known, in general, not to have a simple mechanical interpretation: it cannot be interpreted as the tensor generating forces on surface elements.
The Cauchy stress corresponds to the usual definition of a stress: it gives surface forces when contracted on a vector normal to a surface which is prescribed in the physical space.
It can be shown that the Cauchy stress can be written as a derivative of the energy with respect to the deformation gradient tensor: Unlike the second Piola-Kirchoff stress (i.e. the thermodynamic tension) the Cauchy stress does not depend on the choice of a reference state. This equality holds in the limit F → 1: in this sense it is an Eulerian derivative of the energy. As we will show below, the second Piola-Kirchoff stress becomes identical to the Cauchy stress in the limit where the current and reference configuration are identified.
In order to relate these two definitions of the stress tensor, we need to be able to switch between ηand F -derivatives. From the definition (11) of the Green-Saint Venant strain tensor, infinitesimal displacements read: Since η is symmetric, in a d-dimensional problem, it has only d(d + 1)/2 independent components, while F has d 2 . The rotational degrees of freedom are taken into account by considering the antisymmetric infinitesimal displacements: [40,41] which is the generator of infinitesimal rotations, and has d(d−1)/2 independent components.
From equations (43) and (44), it comes: Using this relation, the differential form of an arbitrary function of F , A(F ), reads: Using the property that the contraction of a symmetric with an antisymmetric tensor vanishes and the fact that η is symmetric, we see that the derivative of a function A with respect to η is: Incidentally, we also find from the preceding calculation that iff A only depends on η, hence does not vary under infinitesimal rotations, it verifies: With these formulae in hand, let us come back to the definition (42) of the Cauchy stress.
Taking the limit F → 1 (or η → 0), in equation (47) and looking at equation (42) we see that T is symmetric. Using equation (46), it appears that it equals the thermodynamic tension in this limit, yet in this limit only. We emphasize here that the symmetry of the stress tensor results directly from the invariance of the interaction potential under elementary rotations: the energy functional is expected to be a function of {r i } and η only: its derivatives with respect to the components of ω vanish. This is not, however, equivalent to the energy functional being invariant under global rotations (the Bravais cell-being a parallelepiped-is not). We note, however, that in general ∂W/∂F (evaluated off the reference configuration) need not be symmetric: this may explain the observation of non-symmetric stresses in [10,11].
Barron and Klein [41] provide the following relation between the Cauchy stress and the thermodynamic tension: This relation holds for finite deformations. It can be obtained, following these authors, after a Taylor expansion of the energy versus the different strain tensors. We provide here a different derivation which relies on an interesting property concerning the transport of derivatives.
As usual, a reference configurationh is given and the system is strained to a current configuration h. In order to define the Cauchy stress around the configuration h, we will need to consider h as a new reference configurationh ′ = h and also consider a new current We need to define strain tensors between each pair of these configuration: the deformation gradient tensors are denoted These transformations are illustrated schematically in figure 7. Likewise, the Green-Saint Venant tensors are denoted: η ′′ , η ′ and η. We have the property: F ′ = F ′′ .F . For any arbitrary function A(F ), and for fixed F , we can write: Using this and equation (46) it is then an easy check that: This relation simply means that the tensorial derivative with respect to η transforms as a tensor under a change of reference configuration. We just saw in equation (49) that this is not true of the derivative with respect to F .
Going back to the definition of the Cauchy stress (equation (42)) and taking the limit and using equation (50) in the limit h →h ′ (or F →F ), we recover equation (48).
It is useful to "specialize" the previous expressions for the situations when the energy, or any observable can be parameterized as A({r ij = F .r ij }). We first have: The factor 1/2 results from the fact that each pair is counted twice with the convention of implicit summation over repeated indices. Plugging equation (51) (and its transpose) into equation (46) we immediately obtain, We see from equation (50) that this expression consist of a reference-independent formula, which is transported backward in the reference configuration as in equation (50).
Specializing the functional A to A = W , we can now provide an expression for the Cauchy stress tensor, using either equation (48) and (52) or alternatively (42) and (51). It reads: This is a generalization of the Kirkwood formula to the case of an arbitrary n-body interaction potential. Note that specialization to systems with pairwise interactions immediately yields the standard expression.
Specializing the function A to A = r β ij , we can furthermore obtain a useful formula. Indeed, in the limit F → 1, equation (52) reduces to:

Elastic stiffnesses and constants
By definition, elastic constants are second derivatives of the energy functional with respect to the Green-Saint Venant strain tensor: C αβκχ = 1 V ∂ 2 W ∂η αβ ∂η κχ η→0 (55) In this expression, the energy functional W can be any given function of the tensor, η. For example, it can be the energy functional of an atomistic system W (η) =Ů({r ij }, η) for fixed positions of the particles in a reference configuration. Derivatives of such a functional were denoted as partial derivatives ofŮ throughout the text. But W could also be defined as the energy of an atomistic system following deformation induced changes of a given minimum.
In this case the partial derivatives of W would correspond to total derivatives for the energy functionalŮ . The algebra presented here does not depend on these considerations but only on the existence of a function W (η).
Since we can commute the order of the derivatives, the elastic constants verify C αβκχ = C κχαβ , and since η is symmetric, C αβκχ = C βακχ = C αβχκ . The second order expansion of the energy with respect to the components of η thus reads: whereT αβ =T βα = (1/V ) ∂W/∂F αβ | η→0 is the stress in the reference configuration. Using (41) and (56), we find for the thermodynamic tension: =T αβ + C αβκχ u κχ + O(u 2 κχ ) The elastic constant should not be confounded with the quantities which enter the stressstrain relations-hence, the wave equation. These are the elastic stiffnesses, c αβκχ , and appear in the expansion of the energy with respect to tensor u = F − 1: [41] W −W The relation of the elastic stiffnesses to the elastic constants is obtained by replacing the definition (11) of η in equation (56) and expanding in terms of u: [41] c Note that this expression indicates that elastic stiffnesses do not enjoy the full symmetry of the elastic constants. Instead, they verify: [41] c αβκχ − c κχαβ =T κχ δ αβ −T αβ δ κχ .
The Cauchy stress can now be obtained from (57) by use of equations (48) and (60): [41] T αβ =T αβ + c αβκχ +T αχ δ βκ −T αβ δ κχ u κχ + O(η 2 κχ ) Before proceeding to the derivation of an expression for the Born term, let us derive a general formula for second derivatives with respect to the components of η. We first write equation (46) with indices explicitly expressed: The derivative of this expression with respect to η αβ reads: after application of equation (46). In the limit η → 0, we thus have for the second derivatives, after a further application of equation (46): We can now provide microscopic expressions for the Born approximation to elastic constants. We take the energy functional W to be the energyŮ for some fixed positions of the particles in a reference frame. Using equation (52) and (51) This general formula permits relating Ξ i,κχ to elementary contributions ∂ 2 U ∂r ij ∂r kl of pairs to the Hessian.
Inserting this expression in the sum (69), the second term disappears because the system is at mechanical equilibrium. This yields: It is now an easy task to show that the terms in equation (67) involving contributions from the stress tensor disappear, which leads to the expression: V ij (r ij c ij − t ij ) r ij n α ij n β ij n κ ij n χ ij (72) * Current address:Institut Navier -LMSGC, 2 allée Képler, 77420 Champs-sur-Marne