Nonequilibrium statistical mechanics of crystals

The local equilibrium approach previously developed by the authors (J Mabillard and P Gaspard 2020 J. Stat. Mech. 103203) for matter with broken symmetries is applied to crystalline solids. The macroscopic hydrodynamics of crystals and their local thermodynamic and transport properties are deduced from the microscopic Hamiltonian dynamics. In particular, the Green–Kubo formulas are obtained for all the transport coefficients. The eight hydrodynamic modes and their dispersion relation are studied for general and cubic crystals. In the same twenty crystallographic classes as those compatible with piezoelectricity, cross effects coupling transport between linear momentum and heat or crystalline order are shown to split the degeneracy of damping rates for modes propagating in opposite generic directions.


I. INTRODUCTION
Crystals manifest long-range order by the spatial periodicity of their atomic structure, which can be classified into 14 Bravais lattices, 32 crystallographic point groups, and 230 space groups [1][2][3].In contrast, fluids have uniform and isotropic properties and are thus symmetric under continuous spatial translations and rotations.According to Goldstone's theorem, the breaking of the three-dimensional continuous group of spatial translations in crystals implies the existence of three slow modes, called Nambu-Goldstone modes, in addition to the five slow modes arising from the fundamental conservation laws of mass, energy, and linear momentum [4][5][6][7].In bulk matter, each mode is characterized by a dispersion relation between its frequency (or rate) and its wave number (or wave length).The slow modes, also called hydrodynamic modes, are characterized by the vanishing of their dispersion relation with the wave number.The mode is propagative (respectively, diffusive) if the dispersion relation vanishes linearly (respectively, quadratically) with the wave number.As a consequence of continuous symmetry breaking, there exist eight hydrodynamic modes in crystals: two longitudinal sound modes, four transverse sound modes, one heat mode, and an additional mode, which has been identified as the mode of vacancy diffusion [8].With these considerations, the macroscopic hydrodynamics of crystals is well established since the seventies [8,9].Earlier treatments of hydrodynamics in solids had identified seven modes only, because the atoms were supposed to move without leaving their lattice cell.However, thermal fluctuations may induce the motion of atoms out of their lattice cell, yielding point defects called vacancies and interstitials in the crystal [1][2][3], which is the mechanism of generating the vacancy diffusion mode.
If the diffusive parts of the dispersion relations are neglected, the modes are either propagative or static, so that there is no energy dissipation and entropy is conserved.Instead, the diffusivity of the modes is associated with the transport coefficients, which are responsible for damping and irreversible entropy production.An issue of paramount importance is to deduce these macroscopic properties from the underlying microscopic motion of atoms, which is the programme of nonequilibrium statistical mechanics.
Since the nineties, the microscopic expression is known for the crystalline order parameter, which is the displacement vector field, providing the statistical-mechanical formulation of hydrodynamics in crystals [37,38].Recently, this formulation has led to a systematic study of elastic properties in nonideal crystals, i.e., crystals with point defects [39,40].In addition to elasticity, the transport properties of nonideal crystals can also be investigated.Using the local equilibrium approach, the Authors have recently developed a unified statistical-mechanical theory for hydrodynamics in phases with broken symmetries, in particular, deducing all the Green-Kubo formulas for the transport coefficients from the microdynamics [41].In certain non-centrosymmetric anisotropic phases, this systematic approach combined with Curie's principle [42] and Onsager-Casimir reciprocal relations [43][44][45] predicts the existence of cross effects coupling the transport of linear momentum to those of heat and the order emerging from symmetry breaking.
In the present paper, our purpose is to apply the systematic approach of reference [41] to crystals in order to obtain their thermodynamic and transport properties.The consequences of the crystallographic symmetries on the transport properties are investigated.Transport coefficients are identified forming tensors of rank two, three, and four.Furthermore, the dispersion relations of the eight hydrodynamic modes are obtained for crystallographic classes where rank-three tensors are or are not vanishing.
The plan of the paper is the following.Section II is devoted to the formulation of the microscopic Hamiltonian dynamics and the construction of the crystalline order parameter.Section III presents the main results of reference [41] about the local equilibrium approach for the nonequilibrium statistical mechanics of crystals and the deduction of their hydrodynamic equations.The eight hydrodynamic modes and their dispersion relation are obtained in section IV for general crystals and in section V for cubic crystals.Section VI is concluding the paper.

II. MICROSCOPIC DESCRIPTION
In this section, the Hamiltonian microdynamics is formulated for crystals, the densities obeying local conservation laws are introduced at the microscopic level of description, and the crystalline order parameter is constructed on the basis of symmetry breaking from continuous to discrete spatial translations.

A. The Hamilton equations of motion
At the microscopic scale, the crystal is composed of N atoms of mass m.The atoms have the positions and momenta Γ = (r i , p i ) N i=1 ∈ R 6N .They are mutually interacting by the energy potential V (r ij ) with r ij = r i − r j , so that the total energy of the crystal is given by the following Hamiltonian function, Consequently, the motion of the atoms is governed by Hamilton's equations where is the a th component of the interaction force exerted on the i th atom by the j th atom (for a = x, y, z and i, j = 1, 2, . . ., N ).Since the Hamiltonian function is invariant under spatiotemporal translations and rotations, this dynamics is known to conserve the total energy E = H, the total momentum P = i p i , and the total angular momentum L = i r i × p i , in addition to the total mass M = mN .We note that, for computational simulation with molecular dynamics, the N -body system can be considered on a torus with periodic boundary conditions [46,47].

B. Locally conserved quantities
In order to describe the local properties in the crystal, the mass density is defined as ρ = mn in terms of the particle density the energy density as and the momentum density as Integrating these densities over the volume of the system gives the total mass, energy, and momentum, which are conserved by the microdynamics.Using Hamilton's equations (II.2), the aforedefined densities can be shown to obey the following local conservation equations, expressed with the corresponding current densities or fluxes, Ĵa ρ (r; Γ) ≡ ĝa (r; Γ) , (II.10) where is a lineal distribution joining the positions r i and r j [18,19,26,27].

C. Continuous symmetry breaking into a crystallographic group
In crystals, the continuous symmetry of the microdynamics under spatial translations, acting on some observable quantity A, is broken into one of the 230 discrete crystallographic space groups.This continuous symmetry breaking generates long-range order in the crystal, whereupon the mass density becomes spatially periodic and anisotropic in the equilibrium crystalline phase, although the mass density is uniform and isotropic in the fluid phase.If the translation vector R is infinitesimal, the transformation (II.14) can be expressed as in terms of the Poisson bracket {•, •} with the total momentum P = i p i .
A priori, the equilibrium properties could be described by the Gibbs grand canonical distribution where H is the Hamiltonian function (II.1), M = mN , β = (k B T ) −1 the inverse temperature, µ the chemical potential, ∆Γ = h 3N (h being Planck's constant), and Ξ the partition function such that the normalization condition p eq (Γ) dΓ = 1 is satisfied with (II.17 The issue is that the probability distribution (II.16) is invariant under the spatial translations (II.15) because {P, H} = 0 and {P, M } = 0, although the continuous translational symmetries are broken in the crystalline phase.
In order to induce symmetry breaking, an external energy potential V (ext) (r) may be introduced, so that the Hamiltonian function is modified into with Unless the external potential is uniform, the associated grand canonical probability distribution is no longer invariant under the spatial translations (II.15).Therefore, taking the mean value of equation (II.15) for the particle density A = n(r; Γ) over the distribution (II.19) gives because of the explicit symmetry breaking by the external potential.Under pressure and temperature conditions where the crystalline phase is thermodynamically stable, the limit ǫ → 0 could be taken in order to remove the external potential and the periodic equilibrium density of the crystal would be obtained as n eq (r) ≡ lim ǫ→0 n(r; Γ) eq,ǫ . (II.21) In this case, the crystalline phase emerges by spontaneous symmetry breaking and its long-range order is characterized by the periodic equilibrium density n eq (r).Accordingly, the spatial structure formed by the atoms becomes periodic in space.Nevertheless, the center of mass R cm ≡ (1/N ) N i=1 r i of the whole crystal can undergo spatial translations.In this regard, we note that the Hamiltonian function (II.1) can be expressed as H = P 2 /(2M ) + H rel in terms of the total kinetic energy of the crystal center of mass and the Hamiltonian function H rel ruling the motion of its atoms relative to the center of mass.In the symmetric grand canonical ensemble (II.16), the center of mass is uniformly distributed in space with a Maxwellian distribution of its velocity, so that its trajectories are free flights, R cm (t) = R cm (0) + tP/M , and the symmetry can only be broken in the frame moving with the center of mass.In the presence of such spontaneous breaking of continuous symmetry, there should exist some external force fields exerting arbitrarily small changes in the total energy of the crystal.In general, an external force field can be described by the external energy potential V (ext) (r), so that the total energy is changed by the mean external energy H (ext)  eq,0 = V (ext) (r) n eq (r) dr (II.22) in the limit ǫ → 0. Remarkably, this mean external energy is strictly equal to zero if the external potential is taken as for some constant vector C. Indeed, for this energy potential, an integration by parts leads to H (ext) eq,0 = C • ∇ ∇ ∇n eq (r) n eq (r) dr = − n eq (r) C • ∇ ∇ ∇n eq (r) dr = − H (ext) eq,0 = 0 , (II.24) hence the result.For such external potentials, the perturbation may induce crystal formation, but does not change the total energy with respect to the value of the unperturbed Hamiltonian.Since its presence costs no energy on average, the external potential (II.23) is leading to the construction of the crystalline order parameter, as shown below.

D. Local order parameter in crystals
At the macroscale, the crystal should be described in terms of macrofields that slowly vary in space on scales larger than the periodic crystalline structure.This feature can be expressed by requiring that the macrofields have no Fourier mode outside the first Brillouin zone B. In general, any function can be decomposed into Fourier modes according to For a macrofield x, we thus require that x(q) = I B (q) x(q) , (II.27) where I B (q) is the indicator function of the first Brillouin zone of the crystalline lattice.We note that, for any arbitrary function f (q), the transformation f (q) → I B (q) f (q) is a projection into the functional space of macrofields.In the position space, this transformation is expressed as in terms of the function (II. 29) This function satisfies the property that R 3 ∆(r − r ′ ) dr ′ = 1, as for a Dirac delta distribution.Moreover, the function (II.29) is real because the first Brillouin zone is symmetric under the inversion q → −q.Indeed, the first Brillouin zone is the Wigner-Seitz primitive cell of the reciprocal lattice, which is a Bravais lattice.A Wigner-Seitz primitive cell has the full symmetry of its Bravais lattice and the point group of a Bravais lattice always includes the inversion.Therefore, we have that As expected, the transformation (II.28) is also a projection because as can be shown using the definition (II.29).According to these considerations, the property (II.27) for the macrofield x becomes In order to identify the local order parameter associated with the continuous symmetry breaking in crystals, we consider the following external perturbation obtained by extending the constant vector C in equation (II.23) into a macrofield C(r) and by substracting the equilibrium value for ǫ = 0, so that H (ext) eq,0 = 0. Accordingly, this external perturbation costs arbitrarily small energy in the limit where C(r) becomes constant in space.Since C(r) is a macrofield, it satisfies equation (II.32) with the function (II.29), whereupon the external perturbation (II.33) can be written in the following form, (II.34) using the property (II.30).At the macroscale, such an external perturbation of the crystal would be described as in terms of the symmetric tensor φ φ φ = φ φ φ T and the strain tensor where û(r; Γ) is the displacement vector field, here supposed to be defined at the microscopic level of description.The tensor φ φ φ = (φ ab ) describes the stress applied to the crystal by the external perturbation because an integration by parts of equation (II.35) with the definition (II.36) gives after introducing the force density field Comparing equations (II.34) and (II.37), we see that the vector field C(r) should correspond to the force density field f (r) and the rest to the microscopic displacement vector field û(r; Γ).This latter should also satisfy the property that, under a homogeneous dilatation r → λr of the lattice by a factor λ, the mean value of the displacement vector field should be given by û(r; Γ) = (λ − 1)r.Accordingly, the microscopic expression for the displacement vector field should be taken as with the tensor N N N = (N ab ) given by where v is the volume of a primitive unit cell of the lattice.Since this tensor is symmetric N N N = N N N T , the force density field can be identified with In cubic crystals, the microscopic expression of references [37,38] for the displacement vector is recovered, in which case the tensor (II.40) should be proportional to the identity tensor: Moreover, the equilibrium density can be expanded into lattice Fourier modes as in terms of the vectors G of the reciprocal lattice.Using the inverse Fourier transform (II.25) and the definition (II.29), we find that which is precisely the expression given in references [37,38] for cubic lattices.We note that the displacement vector field is vanishing in fluid phases where the equilibrium density is uniform and ∇ ∇ ∇n eq = 0, as expected if the continuous symmetry is not broken.In Cartesian components, the local order parameter (II.39) is given by As a consequence of equation (II.7), its time evolution can be expressed as in terms of the decay rate Now, the time evolution equation of the strain tensor ûab = (∇ a ûb +∇ b ûa )/2 introduced in equation (II.36) takes the form of the local conservation equation with the associated current density Equation (II.48) along with equations (II.7), (II.8), and (II.9) are ruling the microscopic hydrodynamics of crystals and they can be written in general form as with the following densities and current densities, Since equation (II.48) is the direct consequence of equation (II.46) for the displacement vector field combined with the definition (II.36) of the strain tensor, equations (II.7), (II.8), (II.9), and (II.46) form the minimal set of eight equations ruling the eight hydrodynamic modes, including the five modes resulting from the five fundamentally conserved quantities (i.e., mass, energy, and momentum) and the three additional Nambu-Goldstone modes generated by the spontaneous symmetry breaking of three-dimensional continuous spatial translations in crystals.
Accordingly, the methods developed in reference [41] for general continuous symmetry breaking in matter can here be applied to crystals, where three continuous symmetries are broken.The local order parameters denoted xα in reference [41] correspond to the three components of the microscopic displacement vector ûa defined by equation (II.45).The gradients ûaα = ∇ a xα of the order parameters correspond to the symmetric strain tensor ûab = ûba defined by equation (II.36), and the conjugated fields φ aα to the symmetric tensor φ ab giving the force density (II.38).In order to apply the results of reference [41] to crystals, the Greek indices α, β, . . .should thus be replaced by Latin indices a, b, . . .= 1, 2, 3 = x, y, z, and the tensors ûaα and φ aα should moreover be symmetrized.

A. Local equilibrium distribution
In order to deduce the macroscopic equations from the microscopic dynamics, we consider the local equilibrium distribution expressed in terms of the densities ĉ = (ĉ α ) defined in equation (II.51), the conjugated fields λ = (λ α ), the integral f * g ≡ dr f (r) g(r) over the volume of the system, and the functional such that the local equilibrium distribution is normalized to the unit value using equation (II.17).The mean values of the densities are thus given by the following functional derivatives with respect to the conjugated fields, The Legendre transform of the functional (III.2) gives the entropy functional in units where Boltzmann's constant is equal to k B = 1 [14,25,27].The comparison with the Euler thermodynamic relation known in crystals [8,9] allows us to identify (up to possible corrections going as the gradient square) the conjugated fields as where β = β(r, t) is the local inverse temperature, µ(r, t) the local chemical potential, v a (r, t) the velocity field, and φ ab (r, t) the field introduced in equation (II.35).Accordingly, the entropy (III.4) can be written as S(c) = s(c) dr + O(∇ 2 ) in terms of the entropy density given by Euler's relation, and the functional (III.2) as Ω = β p dr + O(∇ 2 ) with the hydrostatic pressure p(r, t).The formalism is thus consistent with known local equilibrium thermodynamics in crystals.

B. Time evolution
The time evolution of any phase-space probability distribution is ruled by Liouville's equation ∂ t p t = −Lp t , where L(•) ≡ {•, H} is the Liouvillian operator defined as the Poisson bracket with the Hamiltonian function of the microscopic dynamics.Since the Hamiltonian function (II.1) is time independent, the probability density at time t is given by p t = exp(−Lt)p 0 in terms of the initial density p 0 .This latter is taken as a local equilibrium distribution (III.1) with some initial conjugated fields λ λ λ 0 .However, the probability density does not keep the form (III.1) of a local equilibrium distribution during the time evolution.Nevertheless, the phase-space dynamics is point-like, so that it is possible to express the probability density at time t as [14,16,27] p t (Γ) = e −Lt p leq (Γ; λ 0 ) = p leq (Γ; λ t ) e Σt(Γ) , (III.7) by multiplying the local equilibrium distribution corresponding to time-evolved conjugated fields λ λ λ t with the exponential of the following quantity, The mean value of this quantity represents the entropy (III.4) that is produced during the time interval t which can be shown to be always non-negative in agreement with the second law of thermodynamics [27].
Since the system is isolated, the time derivative of the entropy is giving the entropy production rate dS/dt = d i S/dt ≥ 0, which can thus be calculated using equation (III.8).Now, the macroscopic local conservation equations can be obtained by averaging their microscopic analogies (II.50) over the time-evolved probability distribution (III.7),leading to where c α are the mean densities (III.3) at time t and are respectively the dissipativeless and dissipative current densities [14,27].This identification is justified because they satisfy the following relations, Indeed, equation (III.13)shows that the mean values of the microscopic current densities over the local equilibrium distribution at time t do not contribute to the entropy production rate.Therefore, dissipation (i.e., entropy production) is generated according to equation (III.14) by the contributions (III.12) to the mean values of the microscopic current densities over the full probability distribution (III.7) at time t.
In this regard, it is justified to identify equations (III.11) and (III.12) as the dissipativeless and dissipative current densities.These latter are thus expected to provide the transport coefficients as Green-Kubo formulas after their expansion to first order in the gradients of the conjugated fields.

C. Dissipativeless current densities
The dissipativeless time evolution is obtained by considering the local equilibrium mean values of the densities and current densities.
In particular, the local equilibrium mean values of the densities for momentum and energy are respectively given by ĝa leq = ρ v a , (III.15) in terms of the mean mass density ρ = ρ leq = m n leq and the velocity field v = (v a ), e 0 denoting the internal energy in the frame moving with the crystal element.As a consequence of (III.15),equation (II.7) becomes the continuity equation ∂ t ρ + ∇ ∇ ∇ • (ρv) = 0, expressing the local conservation of mass.Furthermore, the microscopic current densities are averaged over the local equilibrium distribution (III.1) to obtain the dissipativeless current densities (III.11).For the decay rate (II.47) of the microscopic displacement vector (II.45), the local equilibrium mean value is given by Ju a (r, t) ≡ Ĵu a (r; Γ) leq,λt = −v a (r, t) , (III.17) as shown in Appendix A. In references [8,41], the local equilibrium mean value of the order parameter was supposed to have the following general form, with some coefficients A ab and B abc , such that ∇ a A bc = 0.The comparison with the microscopic result (III.17)shows that, for crystals, these coefficients are equal to As shown in reference [41], the local equilibrium mean values of the current densities for momentum and energy are thus given by ) in terms of the reversible stress tensor where p is the hydrostatic pressure and φ ab the tensorial field conjugated to the strain tensor (II.36).Therefore, the dissipativeless rate and current densities can be expressed by equations (III.17),(III.20), and (III.21).

D. Dissipative current densities
In crystals, the dissipative current densities and rates are defined by J a c α = (J a q , J a g b , J u a ) with equation (III.12) and the heat current density They can be calculated with the methods of reference [41].At first order in the affinities or thermodynamic forces A a c α = (A a q , A a g b , A u a ) defined by the following gradients [48][49][50][51][52], the dissipative current densities and rates can be expressed as in terms of the linear response coefficients given by the following Green-Kubo formulas, and δ X ≡ X − X leq,λt .We note that the microscopic global energy current (III.28)determines the heat current density (III.23), so that we shall use the notation c α = q in the left-hand side of Eq. (III.26) if c α = e in its right-hand side.

E. Implications of time-reversal symmetry
According to the symmetry H(ΘΓ) = H(Γ) of the Hamiltonian function under the time-reversal transformation Θ(r i , p i ) = (r i , −p i ), the linear response coefficients (III.26)should obey the Onsager-Casimir reciprocal relations [43][44][45]] where ǫ a c α = ±1 if δ Ĵa c α is even or odd under time reversal (and there is no Einstein summation here).The quantities Ĵa e , Ĵu a , and ĝa are odd, while Ĵa g b , ê, and ρ are even.As a consequence, the Onsager-Casimir reciprocal relations with ǫ a c α ǫ b c β = +1 are giving and those with Therefore, time-reversal symmetry is significantly reducing the number of independent linear response coefficients.Since the coefficients (III.33)form an antisymmetric linear response matrix, they do not contribute to entropy production [41].
Because of equations (III.32) and (III.33), the other linear coefficients are given by and we have the following symmetries κ ab = κ ba , ζ ab = ζ ba , and η abcd = η cdab .The coefficients κ ab are the heat conductivities, η abcd the viscosities, and ζ ab are strain friction coefficients.Consequently, the dissipative current densities and rates take the following forms [60],

G. Implications of crystallographic symmetries
According to Curie's principle [42], the tensorial properties should be symmetric under the transformations of the crystallographic group.This principle applies, in particular, to the rank-three tensors (III.37) and (III.38)describing cross effects coupling the transport of momentum to those of heat and crystalline order.In isotropic phases, rank-three tensors are always vanishing, because such phases have symmetry centers.However, this is no longer the case in anisotropic phases, as illustrated with the phenomenon of piezoelectricity, which is also described by a rank-three tensor [53].Among the 32 possible crystallographic point groups (also called classes), 20 of them are known to allow for non-vanishing rank-three tensors:

H. Vacancy concentration
In the crystal at equilibrium, the macrofield of particle density is uniform and equal to the component G = 0 of the lattice Fourier expansion (II.43) and, equivalently, to the macrofield (II.32) corresponding to the periodic equilibrium density n eq (r): Under nonequilibrium conditions, the macrofield of particle density may deviate with respect to its equilibrium value by two possible mechanisms: (1) the lattice dilatation or contraction corresponding to the strain ∇ ∇ ∇ • u = ∇ a u a ; (2) vacancies or interstitials, decreasing or increasing the occupancy of the lattice cells by particles [8,9,37,38].To describe these latter, a macrofield giving the density of vacancies, also called vacancy concentration, can be defined as ĉ(r; Γ) ≡ − ∆(r − r ′ ) [n(r ′ ; Γ) − n eq,0 + n eq,0 ∇ a ûa (r ′ ; Γ)] dr ′ . (III.47) The time evolution of this macrofield is driven by the local conservation equation (II.7) for the mass density ρ = mn and by equation (II.46) for the displacement vector (II.45).The local equilibrium mean value of the vacancy concentration can be expressed as c(r, t) ≡ ĉ(r; Γ) leq,λ λ λt = −δn(r, t) − n eq,0 ∇ a u a (r, t) , (III.48) where δn(r, t) denotes the macrofield giving the deviation of the mean particle density with respect to its equilibrium value (III.46).Accordingly, the dynamics of the vacancy concentration (III.47) is driven by the evolution equations (II.7) and (II.46) for the particle density n = ρ/m and the displacement vector ûa .

IV. MACROSCOPIC DESCRIPTION
Here, the eight hydrodynamic modes of the crystal are deduced from the linearized macroscopic equations and their dispersion relations are obtained by using expansions in powers of the gradients with respect to the dissipativeless elastic dynamics.

A. Macroscopic equations
As shown in the previous section III, the macroscopic equations ruling the hydrodynamics of crystals can be obtained from the microscopic local conservation equations (II.7), (II.8), and (II.9) for mass, energy, and momentum, combined with the evolution equation (II.46) for the microscopic displacement vector (II.45), using the local equilibrium distribution (III.1) and systematic expansions in powers of the gradients.The transport coefficients are thus given by the Green-Kubo formulas (III.34)-(III.39),entering the expressions of the dissipative current densities (III.41)-(III.43).Gathering them with the dissipativeless current densities (III.15),(III.17),(III.20), and (III.21)(respectively, for mass, displacement, momentum, and energy) and the heat current density (III.23), the following macroscopic equations are obtained, with the stress tensor (III.22) and the dissipative current densities (III.41)-(III.43).

B. Linearized macroscopic equations
In order to investigate the time evolution of small deviations in the macrofields around equilibrium, the macroscopic equations are linearized, leading to the following forms, respectively, for the mean mass density ρ, the mean internal energy e 0 , the velocity field v b , and the mean displacement vector u a .The mass density macrofield can be decomposed as ρ = m (n eq,0 − n eq,0 ∇ a u a − c) (IV.9) in terms of the mean equilibrium density n eq,0 of the atoms of mass m, the contribution from the trace of the strain tensor u aa = ∇ a u a , and the vacancy concentration c.Introducing the fraction of vacancies as y ≡ c n eq,0 , (IV.10) the mass density ρ is thus determined by the trace of the strain tensor u aa = ∇ a u a and the vacancy fraction y.In particular, the deviations of the mass density with respect to the mean equilibrium mass density ρ ≃ ρ eq,0 = mn eq,0 can be written as δρ = −ρ (∇ a δu a + δy) , (IV.11) where δ stands for either ∂ t or ∇ ∇ ∇.Accordingly, equation (IV.5) can be replaced by the evolution equation for the vacancy fraction by using equation (IV.8) for the displacement vector u a : Furthermore, we may introduce the entropy per unit mass s ≡ S/M in every element of the crystal, which is linked to the energy density e 0 , the mass density ρ, and the pressure p by the following Gibbs relation for small deviations with respect to equilibrium.Using equations (IV.5) and (IV.6), the evolution equation for the entropy per unit mass is thus given by Besides, the velocity field obeys equation (IV.7) and the time evolution of the strain tensor u ab = (∇ a u b + ∇ b u a )/2 can be deduced from equation (IV.8).
In order to close the set of equations (IV.7), (IV.8), (IV.12), and (IV.and, as a consequence, a similar expansion for the stress tensor (III.22).Accordingly, the linearized evolution equation (IV.12) for the vacancy fraction (IV.10) is transformed into with the following coefficients, Finally, equation (IV.7) for the velocity field gives where (IV.34) is the rank-four elasticity tensor divided by the mean mass density ρ = mn eq,0 .
The coefficients denoted with the letter C are conservative (i.e., adiabatic or dissipativeless) properties, those with the letter D are dissipative properties, and those with the letter F are conservative coupling properties.
We note that the evolution equations (IV.18), (IV.23), (IV.28), and (IV.33) form a closed set of linear partial differential equations, which can thus be written in matrix form as where L is a 8 × 8 matrix with elements involving the gradient operator ∇ ∇ ∇ at the first, second, or third power, and acting on the eight-dimensional vector field ψ ψ ψ(r, t), the superscript T denoting the transpose.

C. Hydrodynamic modes
Supposing that the solution of equation (IV.41) has the form ψ ψ ψ(r, t) ∼ exp(ıq • r + zt), the dispersion relations z(q) of the eight hydrodynamic modes are provided by solving the eigenvalue problem: with the 8 × 8 matrix (IV.43) The eigenvalue problem can be solved perturbatively starting from the dispersion relations of the dissipativeless crystal with vanishing coefficients D's and F 's.Such a method corresponds to expanding the dispersion relations in powers of the wave number q.For the dissipativeless crystal, the dispersion relations are linear in the wave number q = q .The next order of the perturbation calculation gives the terms of O(q 2 ), providing the damping rates of the modes.Since the dissipativeless and dissipative current densities are known at leading order in the gradient expansion of the fields, the corrections of O(q 3 ) are not relevant to the calculation.
Accordingly, we consider the following expansions of the matrix (IV.43), its eigenvalues, and its eigenvectors: with and (IV.46) In order to solve the eigenvalue problem at zeroth order, we introduce the 3 × 3 real symmetric matrix M defined with the elements and the three-dimensional vectors N y and N s with the components N a y ≡ C ca vy q c and N a s ≡ C ca vs q c .(IV.48) Since the matrix M is real symmetric, it can be diagonalized by an orthogonal transformation giving the three eigenvalues λ σ with σ = 1, 2, 3, and the three eigenvectors ǫ ǫ ǫ σ , which are supposed to form an orthonormal basis, Since M = O(q 2 ), its eigenvalues are also going as λ σ = O(q 2 ).The eigenvectors ǫ ǫ ǫ σ also depend on the wave vector q, but since they are normalized to the unit value, they only depend on the direction of the wave vector, q/q.At zeroth order, the right eigenvectors ψ ψ ψ α and the left eigenvectors ψ ψ ψ(0) α associated with the eigenvalues z (0) α are obtained by solving the following eigenvalue problem, and they are taken to satisfy the biorthonormality conditions, β = δ αβ for α, β = 1, 2, . . ., 8.They are given by Accordingly, the modes α = 1-6 are the propagative sound modes with z (0) α = O(q), α = 7 is the heat mode, and α = 8 the vacancy diffusion mode.Since the zeroth order is adiabatic (isoentropic), there is no damping of the modes at this approximation.
By standard perturbation theory, the first-order correction to the eigenvalues can be obtained with = −D cd s q c q d + (C ca vs q c ) M −1 ab D cdeb su q c q d q e ; (IV. = −D cd y q c q d + C ca vy q c M −1 ab D cdeb yu q c q d q e .(IV.59) We note that, since λ α = O(q 2 ), ǫ ǫ ǫ α = O(q 0 ), and M = O(q 2 ), all the terms of these corrections are going as z α = O(q 2 ) for α = 1-8.Consequently, the corrections z α and z α+3 with α = 1, 2, 3 are giving the damping rates of the propagative sound modes, while the heat mode and the vacancy diffusion modes are diffusive since ).Although the perturbative calculation can be continued to further corrections of O(q 3 ), they are not relevant since the statistical-mechanical theory is limited to the leading terms in the gradient expansion, thus, neglecting the corrections of O(∇ 3 ) = O(q 3 ).
In the case where the coefficients F 's are vanishing, we note that z α+3 , so that the damping rate is the same for sound modes propagating in opposite directions.However, in crystals of the same crystallographic classes as piezoelectric crystals where the coefficients F 's may be non-vanishing [53], these results show that the degeneracy between these damping rates may be split for sound modes propagating in opposite directions, because the terms with the coefficients F 's have opposite signs in z V. CUBIC CRYSTALS Now, the dispersion relations of the eight hydrodynamic modes are analyzed in detail in the case of cubic crystals.The consequences of the cross effects coupling the transport of momentum to those of heat and crystalline order are investigated.

A. Tensors in cubic crystals
In the case of cubic crystals, symmetric rank-two tensors are proportional to the unit tensor and symmetric rank-four tensors such as the elasticity and viscosity tensors can be expanded as (V.1) the three isoentropic elastic constants c 11 , c 12 , and c 44 , and the three viscosity coefficients η 11 , η 12 , and η 44 , expressed in Voigt's notations.As a consequence, we have in particular that and q x ≡ q 1 , q y ≡ q 2 , and q z ≡ q 3 .Similar decompositions hold for other symmetric rank-four tensors.In particular, the tensor (IV.37) can be expressed in Voigt's notations with the three diffusivities associated with the corresponding viscosity coefficients.The point groups of cubic crystals are O h , O, T h , T d , and T, among which only the cubic crystals with the point groups T d and T can accommodate piezoelectricity [53].The key is that piezoelectricity is described by a rank-three tensor, which is only non-vanishing for the classes T d and T among cubic crystals.In such cases, the tensors χ abc and θ abc are specified by a single non-vanishing coefficient because while the other coefficients are equal to zero [53].These rank-three tensors can thus be expressed as with For cubic crystals in the non-piezoelectric classes O h , O, and T h , these tensors are vanishing because χ = 0 and θ = 0.

B. Hydrodynamic modes in cubic crystals
For cubic crystals, the matrix (IV.43) thus becomes −D y q 2 −D ys q 2 −F cdb yv q c q d −ıD cdeb yu q c q d q e −D sy q 2 −D s q 2 −F cdb sv q c q d −ıD cdeb su q c q d q e ı C vy q a − F cad vy q c q d ı C vs q a − F cad vs q c q d −D cadb v q c q d −C cadb vu q c q d − ıF cadeb vu q c q d q e ıD uy q a ıD us q a δ ab + ıF cab uv q c −D cadb with the following coefficients, D sy ≡ κ T y + ξ G vy , (V.17) rank-four tensors, and odd-order tensors, as established in the thermodynamics of crystals [54].The equilibrium thermodynamic stability of the crystals requires the positivity of the compressibility and the specific heat: whereupon p u < 0 and T s > 0. Furthermore, the Grüneisen parameter which is proportional to the thermal expansivity at constant pressure, is observed to be positive [54].Under the same circumstances, the quantity T u defined in equation (V.29) is thus negative, Vacancies are also known as Schottky defects, in which atoms move to the crystal surface [1][2][3].The number of vacancies can be evaluated as where N is the number of ions in the crystals, ε ≃ 1 eV the binding energy of the atoms to their lattice site, v 0 the volume per ion in the lattice, and p the pressure.In such a case, N v ≃ 10 −5 N at T = 1000 K. Since ε ≫ pv 0 , we have that so that T y > 0 and p y < 0.Moreover, the diffusion coefficient of vacancies is known to behave as D y ≃ D y0 exp(−βE v ) with D y0 ≃ 0.20 cm 2 /s and E v ≃ 2.0 eV for the diffusion of Cu in Cu [2].Therefore, vacancy diffusion coefficients are of the order of magnitude D y ∼ 10 −15 m 2 /s at T = 1000 K.The vacancy diffusion coefficient is always positive.Since it is related by D y = − ξ T y − ζ G vy to T y > 0 and since the first term is expected to be larger than the second term, we should have that ξ < 0. The heat conductivity is also always positive κ > 0.
On the one hand, experimental data are available about the elastic constants c 11 , c 12 , and c 44 [1][2][3].The sound velocities are typically of the order of a few thousand meters per second in crystals, in agreement with the values of the corresponding elastic constants and mass density.On the other hand, the experimental measurements of sound attenuation in solids are also available and they show that the diffusivities associated with the sound damping rates and due in particular to the viscosities are of the same order of magnitude as in liquids, i.e., D v ∼ 10 −6 m 2 /s [55].Heat diffusivity has a similar order of magnitude in solids [56].Now, our aim is to obtain illustrative examples of dispersion relations for the eight hydrodynamic modes in cubic crystals.The 8 × 8 matrix (V.12) is implemented in Mathematica [57].For illustrative purposes, the value of the vacancy diffusion coefficient is taken larger than expected according to the known results reported here above.The cases with vanishing and non-vanishing coefficients F 's are compared.A first set of parameter values is chosen to plot the dispersion relations of the eight modes for vanishing coefficients F 's: The units are nm, ps, and K.A second set of parameters is chosen for non-vanishing coefficients F 's, using (V. [100] [123] Figure V.1: Dispersion relations of the eight hydrodynamic modes in a crystal with the parameters (V.41) and (V.42):(a) in the direction (q x , q y , q z ) = (q, 0, 0); (b) in the direction (q x , q y , q z ) = (q, 2q, 3q)/ √ 14.
The eight dispersion relations for the parameter values (V.41) and (V.42) are plotted in figure V.1.The modes are the two longitudinal sound modes L (±) with the largest speeds, the four transverse sound modes T (±) 1 and T (±) 2 , the heat mode H, and the vacancy diffusion mode V.The heat mode and the vacancy diffusion modes are always diffusive.The mode with Im z α < 0 are labelled with (+) and those with Im z α > 0 with (−) in reference to the direction of propagation with respect to that of the wave vector q.The dispersion relations of the transverse sound modes T (±) 1 and T (±) 2 coincide in the direction [100] of wave vector (q x , q y , q z ) = (q, 0, 0), but they are distinct in the direction [123] of wave vector (q x , q y , q z ) = (q, 2q, 3q)/ √ 14, which is well known [1][2][3].Now, for a cubic crystal in the crystallographic classes of piezoelectricity, the eight dispersion relations are plotted in figure V.2 for the parameter values (V.41) and (V.43).The dispersion relations are the same as in figure V.1(a) in the direction [100] of wave vector (q x , q y , q z ) = (q, 0, 0).However, they differ from those of figure V.1(b) in the direction [123] of wave vector (q x , q y , q z ) = (q, 2q, 3q)/ √ 14.As seen in figure V.2(b), the damping rates of the propagative sound modes are different for opposite propagation directions.The results confirm the expectation from the perturbation calculation of the dispersion relations, showing that the non-vanishing coefficients F 's are splitting the degeneracy between the damping rates (IV.56) and (IV.57).Accordingly, the damping of the sound modes may be different for opposite propagation directions in the crystallographic classes (III.44)-(III.45).and (V.43):(a) in the direction (q x , q y , q z ) = (q, 0, 0); (b) in the direction (q x , q y , q z ) = (q, 2q, 3q)/ √ 14.

VI. CONCLUSION
This paper applies the local equilibrium approach of reference [41] to crystals.Since the three-dimensional group of space translations is broken into one of the 230 crystallographic space groups in these phases, three Nambu-Goldstone modes emerge in addition to the five slow modes associated with the conservation of mass, energy, and momentum.As a consequence, the hydrodynamics of crystals is governed by eight slow modes.
The microscopic form of the displacement vector field, i.e., the crystalline local order parameter, is constructed on the basis of the mechanism of continuous symmetry breaking and the result previously obtained in references [37,38] for cubic crystals is recovered.The local equilibrium approach to the nonequilibrium statistical mechanics of crystals provides the microscopic expressions for the thermodynamic and transport properties of crystals.In particular, Green-Kubo formulas are obtained for all the transport coefficients.
Because of the crystalline anisotropy, these coefficients form tensors of rank two, three, and four.The heat conductivities form a rank-two tensor, as well as the friction coefficients of crystalline order and the coefficients of coupling between heat transport and crystalline order.The viscosity coefficients form a rankfour tensor.Moreover, in the same 20 crystallographic classes as those compatible with piezoelectricity, there may exist coefficients coupling the transport of momentum to those of either heat or crystalline order.These cross effects are described by rank-three tensors.In the 12 other crystallographic phases, these rank-three tensors are vanishing.
The hydrodynamics of crystals around equilibrium is investigated by linearizing the eight macroscopic equations ruling the local conservation laws of mass, energy, and momentum, together with the evolution equations for the three components of the displacement vector.The eight hydrodynamic modes and their dispersion relation are obtained by solving the corresponding eigenvalue problem and by using an expansion in powers of the wave number starting from the dissipativeless crystal as a reference.The dispersion relations are analyzed in detail for cubic crystals.
In this way, the damping rates of the hydrodynamic modes are calculated at second order in the wave number.The cross effects existing in the 20 crystallographic classes of piezoelectricity and coupling momentum to heat or crystalline order are shown to split the degeneracy of the damping rates for the sound modes propagating in opposite generic directions.In these crystallographic classes, the damping of sound waves may thus be stronger in some directions than in the opposite one, as a consequence of this degeneracy splitting.Instead, in the 12 other crystallographic classes, these damping rates remain degenerate.
To conclude, we note that the fluctuating hydrodynamics of crystals can be established by adding Gaussian white noise fields δJ a c α (r, t) to the dissipative current densities J a c α (r, t).These Gaussian white noise fields should have amplitudes given by the linear response coefficients in terms of the transport coefficients according to the fluctuation-dissipation theorem [58,59].In this way, the effects of fluctuations on the transport properties can be studied in crystals, using the results of the present paper.
Furthermore, these results are here deduced in the framework of classical mechanics for the microscopic motion of atoms, but they can be generalized to quantum mechanics.In particular, the classical Green-Kubo formulas can be modified into quantum-mechanical formulas by introducing the Hermitian operators for the densities and current densities [10,[16][17][18]25] and by taking into account their non-commutativity using imaginary-time integrals up to the inverse temperature [32].With such considerations, the results can be extended to the dynamics of crystals in quantum regimes.
S 4 , D 6 , T , T d .(III.45)The 10 classes (III.44) are compatible with pyroelectricity and the 20 classes (III.44) and (III.45) with piezoelectrictity [53].Rank-three tensors are vanishing in the 12 other classes because these point groups are either centrosymmetric (i.e., they contain the inversion r → −r with respect to a symmetry center), or they contain 90 • rotations around axes perpendicular to the faces as for the cubic (or orthohedral) crystallographic group O. Therefore, the transport coefficients (III.37) and (III.38) may be non-vanishing in the 20 crystallographic classes (III.44) and (III.45).
T u Ξ abc δ ef + θ Ξ abd G cdefWe note that p u in equation (V.28) can be expressed in terms of the isoentropic compressibility κ s by p u = −1/κ s .Since the compressibility is related to the elastic constants c 11 and c 12 according to κ −1 because of equation (V.6).In this regard, the rank-four tensor G abcd v,u is defined with the two constants