Nonequilibrium statistical mechanics of crystals
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b Nonequilibrium statistical mechanics of crystals
Joël Mabillard ∗ and Pierre Gaspard † Center for Nonlinear Phenomena and Complex Systems, Université Libre de Bruxelles (U.L.B.),Code Postal 231, Campus Plaine, B-1050 Brussels, Belgium
The local equilibrium approach previously developed by the Authors [J. Mabillard and P. Gas-pard,
J. Stat. Mech. (2020) 103203] for matter with broken symmetries is applied to crystallinesolids. The macroscopic hydrodynamics of crystals and their local thermodynamic and transportproperties are deduced from the microscopic Hamiltonian dynamics. In particular, the Green-Kuboformulas are obtained for all the transport coefficients. The eight hydrodynamic modes and theirdispersion relation are studied for general and cubic crystals. In the same twenty crystallographicclasses as those compatible with piezoelectricity, cross effects coupling transport between linear mo-mentum and heat or crystalline order are shown to split the degeneracy of damping rates for modespropagating in opposite generic directions.
I. INTRODUCTION
Crystals manifest long-range order by the spatial periodicity of their atomic structure, which can beclassified into 14 Bravais lattices, 32 crystallographic point groups, and 230 space groups [1–3]. In contrast,fluids have uniform and isotropic properties and are thus symmetric under continuous spatial translationsand rotations. According to Goldstone’s theorem, the breaking of the three-dimensional continuous groupof 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, andlinear momentum [4–7]. In bulk matter, each mode is characterized by a dispersion relation between itsfrequency (or rate) and its wave number (or wave length). The slow modes, also called hydrodynamicmodes, are characterized by the vanishing of their dispersion relation with the wave number. The mode ispropagative (respectively, diffusive) if the dispersion relation vanishes linearly (respectively, quadratically)with the wave number. As a consequence of continuous symmetry breaking, there exist eight hydrodynamicmodes in crystals: two longitudinal sound modes, four transverse sound modes, one heat mode, and anadditional 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 ofhydrodynamics in solids had identified seven modes only, because the atoms were supposed to move withoutleaving their lattice cell. However, thermal fluctuations may induce the motion of atoms out of their latticecell, yielding point defects called vacancies and interstitials in the crystal [1–3], which is the mechanism ofgenerating the vacancy diffusion mode.If the diffusive parts of the dispersion relations are neglected, the modes are either propagative or static, sothat there is no energy dissipation and entropy is conserved. Instead, the diffusivity of the modes is associatedwith the transport coefficients, which are responsible for damping and irreversible entropy production. Anissue of paramount importance is to deduce these macroscopic properties from the underlying microscopicmotion of atoms, which is the programme of nonequilibrium statistical mechanics.Since the fifties, this programme is carried out for normal fluids, starting from the microscopic expres-sions for the densities of the five fundamentally conserved quantities and their time evolution according toHamiltonian microdynamics [10, 11] and using the local equilibrium approach [12–27] often combined withthe projection-operator method [28, 29]. In this way, Green-Kubo formulas can be deduced for the trans-port coefficients in normal fluids [30–33]. In order to extend this programme to matter phases with brokencontinuous symmetries, the associated order parameter must be identified in the microscopic description ofthese phases [34–36].Since the nineties, the microscopic expression is known for the crystalline order parameter, which isthe displacement vector field, providing the statistical-mechanical formulation of hydrodynamics in crystals ∗ Electronic address: [email protected]; ORCID: 0000-0001-6810-3709. † Electronic address: [email protected]; ORCID: 0000-0003-3804-2110. [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 crystalscan also be investigated. Using the local equilibrium approach, the Authors have recently developed aunified 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 certainnon-centrosymmetric anisotropic phases, this systematic approach combined with Curie’s principle [42] andOnsager-Casimir reciprocal relations [43–45] predicts the existence of cross effects coupling the transport oflinear 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 inorder to obtain their thermodynamic and transport properties. The consequences of the crystallographicsymmetries on the transport properties are investigated. Transport coefficients are identified forming tensorsof rank two, three, and four. Furthermore, the dispersion relations of the eight hydrodynamic modes areobtained 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 Hamilto-nian dynamics and the construction of the crystalline order parameter. Section III presents the main resultsof reference [41] about the local equilibrium approach for the nonequilibrium statistical mechanics of crystalsand the deduction of their hydrodynamic equations. The eight hydrodynamic modes and their dispersionrelation are obtained in section IV for general crystals and in section V for cubic crystals. Section VI isconcluding the paper.
II. MICROSCOPIC DESCRIPTION
In this section, the Hamiltonian microdynamics is formulated for crystals, the densities obeying localconservation laws are introduced at the microscopic level of description, and the crystalline order parameteris 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 positionsand momenta Γ = ( r i , p i ) Ni =1 ∈ R N . They are mutually interacting by the energy potential V ( r ij ) with r ij = k r i − r j k , so that the total energy of the crystal is given by the following Hamiltonian function, H = X i p i m + 12 X i = j V ( r ij ) . (II.1)Consequently, the motion of the atoms is governed by Hamilton’s equations d r ai d t = ∂H∂p ai = p ai m , d p ai d t = − ∂H∂r ai = X j ( = i ) F aij , (II.2)where F aij ≡ − ∂V ( r ij ) ∂r ai (II.3)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 , , . . . , N ). Since the Hamiltonian function is invariant under spatiotemporal translations androtations, this dynamics is known to conserve the total energy E = H , the total momentum P = P i p i , andthe total angular momentum L = P i r i × p i , in addition to the total mass M = mN . We note that, forcomputational simulation with molecular dynamics, the N -body system can be considered on a torus withperiodic boundary conditions [46, 47]. B. Locally conserved quantities
In order to describe the local properties in the crystal, the mass density is defined as ˆ ρ = m ˆ n in terms ofthe particle density ˆ n ( r ; Γ) ≡ X i δ ( r − r i ) , (II.4)the energy density as ˆ e ( r ; Γ) ≡ X i E i δ ( r − r i ) with E i ≡ p i m + 12 X j ( = i ) V ( r ij ) , (II.5)and the momentum density as ˆ g a ( r ; Γ) ≡ X i p ai δ ( r − r i ) . (II.6)Integrating these densities over the volume of the system gives the total mass, energy, and momentum, whichare conserved by the microdynamics.Using Hamilton’s equations (II.2), the aforedefined densities can be shown to obey the following localconservation equations, ∂ t ˆ ρ + ∇ a ˆ J aρ = 0 , (II.7) ∂ t ˆ e + ∇ a ˆ J ae = 0 , (II.8) ∂ t ˆ g b + ∇ a ˆ J ag b = 0 , (II.9)expressed with the corresponding current densities or fluxes, ˆ J aρ ( r ; Γ) ≡ ˆ g a ( r ; Γ) , (II.10) ˆ J ae ( r ; Γ) ≡ X i p ai m E i δ ( r − r i ) + 12 X i In crystals, the continuous symmetry of the microdynamics under spatial translations, T R A [( r i , p i ) Ni =1 ] = A [( r i − R , p i ) Ni =1 ] with R ∈ R (II.14)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 densitybecomes spatially periodic and anisotropic in the equilibrium crystalline phase, although the mass density isuniform and isotropic in the fluid phase. If the translation vector R is infinitesimal, the transformation (II.14)can be expressed as T R A = A + R · { P , A } + O ( R ) (II.15)in terms of the Poisson bracket {· , ·} with the total momentum P = P i p i . A priori , the equilibrium properties could be described by the Gibbs grand canonical distribution p eq (Γ) = 1Ξ ∆Γ e − β ( H − µM ) , (II.16)where H is the Hamiltonian function (II.1), M = mN , β = ( k B T ) − the inverse temperature, µ the chemicalpotential, ∆Γ = h N ( h being Planck’s constant), and Ξ the partition function such that the normalizationcondition R p eq (Γ) dΓ = 1 is satisfied with Z ( · ) dΓ = ∞ X N =0 N ! Z R N ( · ) N Y i =1 d r i d p i . (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 crystallinephase.In order to induce symmetry breaking, an external energy potential V (ext) ( r ) may be introduced, so thatthe Hamiltonian function is modified into H ǫ = H + ǫ H (ext) with H (ext) = Z V (ext) ( r ) ˆ n ( r ; Γ) d r . (II.18)Unless the external potential is uniform, the associated grand canonical probability distribution p eq ,ǫ (Γ , N ) = 1Ξ ǫ ∆Γ e − β ( H ǫ − µM ) (II.19)is no longer invariant under the spatial translations (II.15). Therefore, taking the mean value of equa-tion (II.15) for the particle density A = ˆ n ( r ; Γ) over the distribution (II.19) gives h{ P , ˆ n ( r ; Γ) }i eq ,ǫ = ∂ r h ˆ n ( r ; Γ) i eq ,ǫ = 0 if ǫ = 0 , (II.20)because of the explicit symmetry breaking by the external potential. Under pressure and temperatureconditions where the crystalline phase is thermodynamically stable, the limit ǫ → could be taken in orderto remove the external potential and the periodic equilibrium density of the crystal would be obtained as n eq ( r ) ≡ lim ǫ → h ˆ n ( r ; Γ) i eq ,ǫ . (II.21)In this case, the crystalline phase emerges by spontaneous symmetry breaking and its long-range order ischaracterized by the periodic equilibrium density n eq ( r ) . Accordingly, the spatial structure formed by theatoms becomes periodic in space. Nevertheless, the center of mass R cm ≡ (1 /N ) P Ni =1 r i of the wholecrystal can undergo spatial translations. In this regard, we note that the Hamiltonian function (II.1) canbe expressed as H = P / (2 M ) + H rel in terms of the total kinetic energy of the crystal center of massand the Hamiltonian function H rel ruling the motion of its atoms relative to the center of mass. In thesymmetric grand canonical ensemble (II.16), the center of mass is uniformly distributed in space with aMaxwellian distribution of its velocity, so that its trajectories are free flights, R cm ( t ) = R cm (0) + t P /M ,and the symmetry can only be broken in the frame moving with the center of mass. In the presence ofsuch spontaneous breaking of continuous symmetry, there should exist some external force fields exertingarbitrarily small changes in the total energy of the crystal. In general, an external force field can be describedby the external energy potential V (ext) ( r ) , so that the total energy is changed by the mean external energy h H (ext) i eq , = Z V (ext) ( r ) n eq ( r ) d r (II.22)in the limit ǫ → . Remarkably, this mean external energy is strictly equal to zero if the external potentialis taken as V (ext) ( r ) = C · ∇∇∇ n eq ( r ) (II.23)for some constant vector C . Indeed, for this energy potential, an integration by parts leads to h H (ext) i eq , = Z C · ∇∇∇ n eq ( r ) n eq ( r ) d r = − Z n eq ( r ) C · ∇∇∇ n eq ( r ) d r = −h H (ext) i eq , = 0 , (II.24)hence the result. For such external potentials, the perturbation may induce crystal formation, but does notchange the total energy with respect to the value of the unperturbed Hamiltonian. Since its presence costsno energy on average, the external potential (II.23) is leading to the construction of the crystalline orderparameter, 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 spaceon scales larger than the periodic crystalline structure. This feature can be expressed by requiring thatthe macrofields have no Fourier mode outside the first Brillouin zone B . In general, any function can bedecomposed into Fourier modes according to f ( r ) = Z R d q (2 π ) ˜ f ( q ) e ı q · r , (II.25) ˜ f ( q ) = Z R d r f ( r ) e − ı q · r . (II.26)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 intothe functional space of macrofields. In the position space, this transformation is expressed as f ( r ) → Z R ∆( r − r ′ ) f ( r ′ ) d r ′ (II.28)in terms of the function ∆( r − r ′ ) ≡ Z B d k (2 π ) e ı k · ( r − r ′ ) . (II.29)This function satisfies the property that R R ∆( r − r ′ ) d r ′ = 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 Bravaislattice always includes the inversion. Therefore, we have that ∆( r − r ′ ) ∗ = ∆( r ′ − r ) = ∆( r − r ′ ) . (II.30)As expected, the transformation (II.28) is also a projection because ∆( r − r ′ ) = Z R ∆( r − r ′′ ) ∆( r ′′ − r ′ ) d r ′′ , (II.31)as can be shown using the definition (II.29). According to these considerations, the property (II.27) for themacrofield x becomes x ( r ) = Z R ∆( r − r ′ ) x ( r ′ ) d r ′ . (II.32)In order to identify the local order parameter associated with the continuous symmetry breaking in crystals,we consider the following external perturbation H (ext) = Z V (ext) ( r ) [ˆ n ( r ; Γ) − n eq ( r )] d r = Z C ( r ) · ∇∇∇ n eq ( r ) [ˆ n ( r ; Γ) − n eq ( r )] d r , (II.33)obtained by extending the constant vector C in equation (II.23) into a macrofield C ( r ) and by substractingthe equilibrium value for ǫ = 0 , so that h H (ext) i eq , = 0 . Accordingly, this external perturbation costsarbitrarily 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 bewritten in the following form, H (ext) = Z Z ∆( r − r ′ ) C ( r ′ ) · ∇∇∇ n eq ( r ) [ˆ n ( r ; Γ) − n eq ( r )] d r d r ′ = Z d r C ( r ) · Z d r ′ ∆( r − r ′ ) ∇∇∇ ′ n eq ( r ′ ) [ˆ n ( r ′ ; Γ) − n eq ( r ′ )] , (II.34)using the property (II.30).At the macroscale, such an external perturbation of the crystal would be described as H (ext) = − Z d r φφφ ( r ) : ˆ u ( r ; Γ) (II.35)in terms of the symmetric tensor φφφ = φφφ T and the strain tensor ˆ u ≡ (cid:0) ∇∇∇ ˆ u + ∇∇∇ ˆ u T (cid:1) , (II.36)where ˆ u ( r ; Γ) is the displacement vector field, here supposed to be defined at the microscopic level ofdescription. The tensor φφφ = ( φ ab ) describes the stress applied to the crystal by the external perturbationbecause an integration by parts of equation (II.35) with the definition (II.36) gives H (ext) = Z d r ( ∇∇∇ · φφφ ) · ˆ u = − Z d r f · ˆ u , (II.37)after introducing the force density field f ≡ −∇∇∇ · φφφ . (II.38)Comparing equations (II.34) and (II.37), we see that the vector field C ( r ) should correspond to the forcedensity field f ( r ) and the rest to the microscopic displacement vector field ˆ u ( r ; Γ) . This latter should alsosatisfy the property that, under a homogeneous dilatation r → λ r of the lattice by a factor λ , the meanvalue of the displacement vector field should be given by h ˆ u ( r ; Γ) i = ( λ − r . Accordingly, the microscopicexpression for the displacement vector field should be taken as ˆ u ( r ; Γ) = − NNN − · Z d r ′ ∆( r − r ′ ) ∇∇∇ ′ n eq ( r ′ ) [ˆ n ( r ′ ; Γ) − n eq ( r ′ )] (II.39)with the tensor NNN = ( N ab ) given by N ab ≡ v Z v ∇ a n eq ∇ b n eq d r , (II.40)where v is the volume of a primitive unit cell of the lattice. Since this tensor is symmetric NNN = NNN T , theforce density field can be identified with f ( r ) = − NNN · C ( r ) . (II.41)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: N ab = N δ ab with N ≡ v Z v ( ∇∇∇ n eq ) d r . (II.42)Moreover, the equilibrium density can be expanded into lattice Fourier modes as n eq ( r ) = X G n eq , G e ı G · r (II.43)in terms of the vectors G of the reciprocal lattice. Using the inverse Fourier transform (II.25) and thedefinition (II.29), we find that ˆ u ( r ; Γ) = − ı N X G G n eq , G Z B d k (2 π ) e ı k · r h ˆ˜ n ( k − G ; Γ) − ˜ n eq ( k − G ) i , (II.44)which is precisely the expression given in references [37, 38] for cubic lattices. We note that the displacementvector field is vanishing in fluid phases where the equilibrium density is uniform and ∇∇∇ n eq = 0 , as expectedif the continuous symmetry is not broken.In Cartesian components, the local order parameter (II.39) is given by ˆ u a ( r ; Γ) = − ( NNN − ) ab Z d r ′ ∆( r − r ′ ) ∇ ′ b n eq ( r ′ ) [ˆ n ( r ′ ; Γ) − n eq ( r ′ )] . (II.45)As a consequence of equation (II.7), its time evolution can be expressed as ∂ t ˆ u a + ˆ J u a = 0 (II.46)in terms of the decay rate ˆ J u a ( r ; Γ) = − m ( NNN − ) ab Z d r ′ ∆( r − r ′ ) ∇ ′ b n eq ( r ′ ) ∇ ′ c ˆ g c ( r ′ ; Γ) . (II.47)Now, the time evolution equation of the strain tensor ˆ u ab = ( ∇ a ˆ u b + ∇ b ˆ u a ) / introduced in equation (II.36)takes the form of the local conservation equation ∂ t ˆ u ab + ∇ c ˆ J cu ab = 0 (II.48)with the associated current density ˆ J cu ab ( r ; Γ) = 12 (cid:0) δ ac δ bd + δ ad δ cb (cid:1) ˆ J u d ( r ; Γ) . (II.49)Equation (II.48) along with equations (II.7), (II.8), and (II.9) are ruling the microscopic hydrodynamics ofcrystals and they can be written in general form as ∂ t ˆ c α ( r , t ) + ∇ a ˆ J ac α ( r , t ) = 0 (II.50)with the following densities and current densities, (ˆ c α ) = (ˆ e, ˆ ρ, ˆ g b , ˆ u bc ) and ( ˆ J ac α ) = ( ˆ J ae , ˆ J aρ , ˆ J ag b , ˆ J au bc ) . (II.51)Since equation (II.48) is the direct consequence of equation (II.46) for the displacement vector field com-bined with the definition (II.36) of the strain tensor, equations (II.7), (II.8), (II.9), and (II.46) form theminimal set of eight equations ruling the eight hydrodynamic modes, including the five modes resulting fromthe five fundamentally conserved quantities (i.e., mass, energy, and momentum) and the three additionalNambu-Goldstone modes generated by the spontaneous symmetry breaking of three-dimensional continuousspatial translations in crystals.Accordingly, the methods developed in reference [41] for general continuous symmetry breaking in mattercan here be applied to crystals, where three continuous symmetries are broken. The local order parametersdenoted ˆ x α in reference [41] correspond to the three components of the microscopic displacement vector ˆ u a defined by equation (II.45). The gradients ˆ u aα = ∇ a ˆ x α of the order parameters correspond to the symmetricstrain tensor ˆ u ab = ˆ u 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 Greekindices α, β, . . . should thus be replaced by Latin indices a, b, . . . = 1 , , x, y, z , and the tensors ˆ u aα and φ aα should moreover be symmetrized. III. NONEQUILIBRIUM STATISTICAL MECHANICS In this section, the local equilibrium approach known for normal fluids [12–27] is extended to crystalsusing the results of reference [41]. The Green-Kubo formulas are obtained for all the crystalline transportcoefficients and the vacancy concentration is introduced. A. Local equilibrium distribution In order to deduce the macroscopic equations from the microscopic dynamics, we consider the local equi-librium distribution p leq (Γ; λ ) = 1∆Γ exp [ − λ α ∗ ˆ c α (Γ) − Ω( λ )] , (III.1)expressed in terms of the densities ˆ c = (ˆ c α ) defined in equation (II.51), the conjugated fields λ = ( λ α ) , theintegral f ∗ g ≡ R d r f ( r ) g ( r ) over the volume of the system, and the functional Ω( λ ) = ln Z dΓ∆Γ exp [ − λ α ∗ ˆ c α (Γ)] , (III.2)such that the local equilibrium distribution is normalized to the unit value using equation (II.17). The meanvalues of the densities are thus given by the following functional derivatives with respect to the conjugatedfields, c α ( r ) ≡ h ˆ c α ( r ; Γ) i leq , λ = − δ Ω( λ ) δλ α ( r ) . (III.3)The Legendre transform of the functional (III.2) gives the entropy functional S ( c ) = inf λ [ λ α ∗ c α + Ω( λ )] (III.4)in units where Boltzmann’s constant is equal to k B = 1 [14, 25, 27]. The comparison with the Eulerthermodynamic relation known in crystals [8, 9] allows us to identify (up to possible corrections going as thegradient square) the conjugated fields as λ e = β + O ( ∇ ) , λ ρ = − β µ + O ( ∇ ) , λ g a = − β v a + O ( ∇ ) , λ u ab = − β φ ab + O ( ∇ ) , (III.5)where β = β ( r , t ) is the local inverse temperature, µ ( r , t ) the local chemical potential, v a ( r , t ) the velocityfield, and φ ab ( r , t ) the field introduced in equation (II.35). Accordingly, the entropy (III.4) can be writtenas S ( c ) = R s ( c ) d r + O ( ∇ ) in terms of the entropy density given by Euler’s relation, s = β ( e + p ) − β µ ρ − β v a g a − β φ ab u ab , (III.6)and the functional (III.2) as Ω = R β p d r + O ( ∇ ) with the hydrostatic pressure p ( r , t ) . The formalism isthus 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 = −L p t ,where L ( · ) ≡ {· , H } is the Liouvillian operator defined as the Poisson bracket with the Hamiltonian functionof the microscopic dynamics. Since the Hamiltonian function (II.1) is time independent, the probabilitydensity at time t is given by p t = exp( −L t ) p in terms of the initial density p . This latter is taken as alocal equilibrium distribution (III.1) with some initial conjugated fields λλλ . However, the probability densitydoes 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 −L t p leq (Γ; λ ) = p leq (Γ; λ t ) e Σ t (Γ) , (III.7)by multiplying the local equilibrium distribution corresponding to time-evolved conjugated fields λλλ t with theexponential of the following quantity, Σ t (Γ) ≡ Z t d τ ∂ τ [ λ ατ ∗ ˆ c α (Γ τ − t ) + Ω( λ τ )] . (III.8)The mean value of this quantity represents the entropy (III.4) that is produced during the time interval t h Σ t (Γ) i t = S ( c t ) − S ( c ) ≥ , (III.9)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 d S/ d t =d i S/ d t ≥ , which can thus be calculated using equation (III.8).Now, the macroscopic local conservation equations can be obtained by averaging their microscopic analo-gies (II.50) over the time-evolved probability distribution (III.7), leading to ∂ t c α + ∇ a (cid:0) ¯ J ac α + J ac α (cid:1) = 0 , (III.10)where c α are the mean densities (III.3) at time t and ¯ J ac α ( r , t ) ≡ h ˆ J ac α ( r ; Γ) i leq , λ t , (III.11) J ac α ( r , t ) ≡ h ˆ J ac α ( r ; Γ) (cid:2) e Σ t (Γ) − (cid:3) i leq , λ t (III.12)are respectively the dissipativeless and dissipative current densities [14, 27]. This identification is justifiedbecause they satisfy the following relations, ∇ a λ αt ∗ ¯ J ac α ( t ) = 0 , (III.13) d S d t = ∇ a λ αt ∗ J ac α ( t ) ≥ . (III.14)Indeed, equation (III.13) shows that the mean values of the microscopic current densities over the localequilibrium 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 themean 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 dissipativecurrent densities. These latter are thus expected to provide the transport coefficients as Green-Kubo formulasafter 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 thedensities and current densities.In particular, the local equilibrium mean values of the densities for momentum and energy are respectivelygiven by h ˆ g a i leq = ρ v a , (III.15) h ˆ e i leq = e = e + ρ v / , (III.16)in terms of the mean mass density ρ = h ˆ ρ i leq = m h ˆ n i leq and the velocity field v = ( v a ) , e denoting theinternal 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 displace-ment vector (II.45), the local equilibrium mean value is given by ¯ J u a ( r , t ) ≡ h ˆ J u a ( r ; Γ) i 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 wassupposed to have the following general form, ¯ J u c = − A ac v a − B abc ∇ a v b + O ( ∇ ) (III.18)0with some coefficients A ab and B abc , such that ∇ a A bc = 0 . The comparison with the microscopic re-sult (III.17) shows that, for crystals, these coefficients are equal to A ab = δ ab and B abc = 0 . (III.19)As shown in reference [41], the local equilibrium mean values of the current densities for momentum andenergy are thus given by ¯ J ag b = h ˆ J ag b i leq = ρ v b v a − σ ab + O ( ∇ ) , (III.20) ¯ J ae = h ˆ J ae i leq = e v a − σ ab v b + O ( ∇ ) , (III.21)in terms of the reversible stress tensor σ ab ≡ − p δ ab + φ ab , (III.22)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 ac α = ( J aq , J ag b , J u a ) with equa-tion (III.12) and the heat current density J aq ≡ J ae − v b J ag b − φ ab J u b . (III.23)They can be calculated with the methods of reference [41].At first order in the affinities or thermodynamic forces A ac α = ( A aq , A ag b , A u a ) defined by the followinggradients [48–52], A aq ≡ − T ∇ a T , A ag b ≡ − T ∇ a v b , A u a ≡ − T ∇ b φ ba , (III.24)the dissipative current densities and rates can be expressed as J ac α = X b,β L (cid:0) ac α (cid:12)(cid:12) bc β (cid:1) A bc β , (III.25)in terms of the linear response coefficients given by the following Green-Kubo formulas, L (cid:0) ac α (cid:12)(cid:12) bc β (cid:1) ≡ lim V →∞ k B V Z ∞ d t h δ ˆ J ′ ac α ( t ) δ ˆ J ′ bc β (0) i eq , (III.26)where δ ˆ J ′ ac α ( t ) ≡ Z V δ ˆ J ′ ac α ( r , t ) d r (III.27)are the microscopic global currents defined with δ ˆ J ′ ae ≡ δ ˆ J ae − ρ − ( e + p ) δ ˆ g a , (III.28) δ ˆ J ′ ag b ≡ δ ˆ J ag b + (cid:18) ∂σ ab ∂e (cid:19) ρ, u δ ˆ e + (cid:18) ∂σ ab ∂ρ (cid:19) e , u δ ˆ ρ , (III.29) δ ˆ J ′ u a ≡ δ ˆ J u a + ρ − δ ˆ g a , (III.30)and δ ˆ X ≡ ˆ X − h ˆ X i leq , λ t . We note that the microscopic global energy current (III.28) determines the heatcurrent 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.1 E. Implications of time-reversal symmetry According to the symmetry H (ΘΓ) = H (Γ) of the Hamiltonian function under the time-reversal trans-formation Θ( r i , p i ) = ( r i , − p i ) , the linear response coefficients (III.26) should obey the Onsager-Casimirreciprocal relations [43–45] L (cid:0) ac α (cid:12)(cid:12) bc β (cid:1) = ǫ ac α ǫ bc β L (cid:0) bc β (cid:12)(cid:12) ac α (cid:1) , (III.31)where ǫ ac α = ± if δ ˆ J ac α is even or odd under time reversal (and there is no Einstein summation here). Thequantities ˆ J ae , ˆ J u a , and ˆ g a are odd, while ˆ J ag b , ˆ e , and ˆ ρ are even. As a consequence, the Onsager-Casimirreciprocal relations with ǫ ac α ǫ bc β = +1 are giving L (cid:0) aq (cid:12)(cid:12) bq (cid:1) = L (cid:0) bq (cid:12)(cid:12) aq (cid:1) , L (cid:16) ag b (cid:12)(cid:12)(cid:12) cg d (cid:17) = L (cid:16) cg d (cid:12)(cid:12)(cid:12) ag b (cid:17) , L (cid:0) u a (cid:12)(cid:12) u b (cid:1) = L (cid:0) u b (cid:12)(cid:12) u a (cid:1) , and L (cid:0) aq (cid:12)(cid:12) u b (cid:1) = L (cid:0) u b (cid:12)(cid:12) aq (cid:1) , (III.32)and those with ǫ ac α ǫ bc β = − lead to L (cid:16) ag b (cid:12)(cid:12)(cid:12) cq (cid:17) = −L (cid:16) cq (cid:12)(cid:12)(cid:12) ag b (cid:17) and L (cid:0) u a (cid:12)(cid:12) bg c (cid:1) = −L (cid:0) bg c (cid:12)(cid:12) u a (cid:1) . (III.33)Therefore, time-reversal symmetry is significantly reducing the number of independent linear response coef-ficients. Since the coefficients (III.33) form an antisymmetric linear response matrix, they do not contributeto entropy production [41]. F. Transport properties The independent transport coefficients are thus obtained from the microdynamics with the followingGreen-Kubo formulas, κ ab ≡ T L (cid:0) aq (cid:12)(cid:12) bq (cid:1) = lim V →∞ k B T V Z ∞ d t h δ ˆ J ′ ae ( t ) δ ˆ J ′ be (0) i eq , (III.34) ξ ab ≡ T L (cid:0) aq (cid:12)(cid:12) u b (cid:1) = lim V →∞ k B T V Z ∞ d t h δ ˆ J ′ ae ( t ) δ ˆ J ′ u b (0) i eq , (III.35) ζ ab ≡ T L (cid:0) u a (cid:12)(cid:12) u b (cid:1) = lim V →∞ k B T V Z ∞ d t h δ ˆ J ′ u a ( t ) δ ˆ J ′ u b (0) i eq , (III.36) χ abc ≡ T L (cid:16) aq (cid:12)(cid:12)(cid:12) bg c (cid:17) = lim V →∞ k B T V Z ∞ d t h δ ˆ J ′ ae ( t ) δ ˆ J ′ bg c (0) i eq , (III.37) θ abc ≡ T L (cid:16) ag b (cid:12)(cid:12) u c (cid:17) = lim V →∞ k B T V Z ∞ d t h δ ˆ J ′ ag b ( t ) δ ˆ J ′ u c (0) i eq , (III.38) η abcd ≡ T L (cid:16) ag b (cid:12)(cid:12)(cid:12) cg d (cid:17) = lim V →∞ k B T V Z ∞ d t h δ ˆ J ′ ag b ( t ) δ ˆ J ′ cg d (0) i eq , (III.39)expressed in terms of the microscopic global currents defined by equations (III.27)-(III.30).Because of equations (III.32) and (III.33), the other linear coefficients are given by L (cid:0) u a (cid:12)(cid:12) bq (cid:1) = T ξ ba , L (cid:16) ag b (cid:12)(cid:12)(cid:12) cq (cid:17) = − T χ cab , L (cid:0) u a (cid:12)(cid:12) bg c (cid:1) = − T θ bca , (III.40)and we have the following symmetries κ ab = κ ba , ζ ab = ζ ba , and η abcd = η cdab . The coefficients κ ab are theheat conductivities, η abcd the viscosities, and ζ ab are strain friction coefficients.Consequently, the dissipative current densities and rates take the following forms [60], J aq = − κ ab ∇ b T − χ abc ∇ b v c − ξ ac ∇ b φ bc , (III.41) J ag b = χ cab T ∇ c T − η abcd ∇ c v d − θ abd ∇ c φ cd , (III.42) J u a = − ξ ba T ∇ b T + θ bca ∇ b v c − ζ ac ∇ b φ bc . (III.43)2 G. Implications of crystallographic symmetries According to Curie’s principle [42], the tensorial properties should be symmetric under the transforma-tions 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 crystallineorder. In isotropic phases, rank-three tensors are always vanishing, because such phases have symmetrycenters. However, this is no longer the case in anisotropic phases, as illustrated with the phenomenon ofpiezoelectricity, which is also described by a rank-three tensor [53]. Among the 32 possible crystallographicpoint groups (also called classes), 20 of them are known to allow for non-vanishing rank-three tensors: C , C s , C , C , C , C , C , C , C , C , (III.44) D , D , C , D , D , D , S , D , T , T d . (III.45)The 10 classes (III.44) are compatible with pyroelectricity and the 20 classes (III.44) and (III.45) withpiezoelectrictity [53]. Rank-three tensors are vanishing in the 12 other classes because these point groupsare either centrosymmetric (i.e., they contain the inversion r → − r with respect to a symmetry center),or they contain ◦ 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 inthe 20 crystallographic classes (III.44) and (III.45). 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 tothe periodic equilibrium density n eq ( r ) : n eq , = 1 v Z v n eq ( r ) d r = Z ∆( r − r ′ ) n eq ( r ′ ) d r ′ . (III.46)Under nonequilibrium conditions, the macrofield of particle density may deviate with respect to its equilib-rium 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 byparticles [8, 9, 37, 38]. To describe these latter, a macrofield giving the density of vacancies, also calledvacancy concentration, can be defined as ˆ c ( r ; Γ) ≡ − Z ∆( r − r ′ ) [ˆ n ( r ′ ; Γ) − n eq , + n eq , ∇ a ˆ u a ( r ′ ; Γ)] d r ′ . (III.47)The time evolution of this macrofield is driven by the local conservation equation (II.7) for the mass density ˆ ρ = m ˆ n 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 ) ≡ h ˆ c ( r ; Γ) i leq ,λλλ t = − δn ( r , t ) − n eq , ∇ 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 itsequilibrium 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 ˆ u a . IV. MACROSCOPIC DESCRIPTION Here, the eight hydrodynamic modes of the crystal are deduced from the linearized macroscopic equationsand their dispersion relations are obtained by using expansions in powers of the gradients with respect tothe dissipativeless elastic dynamics.3 A. Macroscopic equations As shown in the previous section III, the macroscopic equations ruling the hydrodynamics of crystals canbe obtained from the microscopic local conservation equations (II.7), (II.8), and (II.9) for mass, energy, andmomentum, 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. Thetransport coefficients are thus given by the Green-Kubo formulas (III.34)-(III.39), entering the expressionsof the dissipative current densities (III.41)-(III.43). Gathering them with the dissipativeless current densi-ties (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, ∂ t ρ + ∇ a ( ρv a ) = 0 , (IV.1) ∂ t e + ∇ a ( e v a − σ ab v b + v b J ag b + φ ab J u b + J aq ) = 0 , (IV.2) ∂ t ( ρv b ) + ∇ a ( ρv a v b − σ ab + J ag b ) = 0 , (IV.3) ∂ t u a − v a + J u a = 0 , (IV.4)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, themacroscopic equations are linearized, leading to the following forms, ∂ t ρ = − ρ ∇ a v a , (IV.5) ∂ t e = − ( e + p ) ∇ a v a + κ ab ∇ a ∇ b T + χ abc ∇ a ∇ b v c + ξ ac ∇ a ∇ b φ bc , (IV.6) ρ ∂ t v b = −∇ b p + ∇ a φ ab + θ abd ∇ a ∇ c φ cd − χ cab T ∇ a ∇ c T + η abcd ∇ a ∇ c v d , (IV.7) ∂ t u a = v a − θ bca ∇ b v c + ξ ba T ∇ b T + ζ ac ∇ b φ bc , (IV.8)respectively, for the mean mass density ρ , the mean internal energy e , the velocity field v b , and the meandisplacement vector u a .The mass density macrofield can be decomposed as ρ = m ( n eq , − n eq , ∇ a u a − c ) (IV.9)in terms of the mean equilibrium density n eq , of the atoms of mass m , the contribution from the trace ofthe strain tensor u aa = ∇ a u a , and the vacancy concentration c . Introducing the fraction of vacancies as y ≡ cn eq , , (IV.10)the mass density ρ is thus determined by the trace of the strain tensor u aa = ∇ a u a and the vacancyfraction y . In particular, the deviations of the mass density with respect to the mean equilibrium massdensity ρ ≃ ρ eq , = mn eq , 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 equationfor the vacancy fraction by using equation (IV.8) for the displacement vector u a : ∂ t y = ∇ a v a − ∇ a ∂ t u a = θ bca ∇ a ∇ b v c − ξ ba T ∇ a ∇ b T − ζ ac ∇ a ∇ b φ bc . (IV.12)Furthermore, we may introduce the entropy per unit mass s ≡ S/M in every element of the crystal, whichis linked to the energy density e , the mass density ρ , and the pressure p by the following Gibbs relation T ρ δ s = δe − e + pρ δρ (IV.13)4for small deviations with respect to equilibrium. Using equations (IV.5) and (IV.6), the evolution equationfor the entropy per unit mass is thus given by T ρ ∂ t s = κ ab ∇ a ∇ b T + χ abc ∇ a ∇ b v c + ξ ac ∇ a ∇ b φ bc . (IV.14)Besides, the velocity field obeys equation (IV.7) and the time evolution of the strain tensor u ab = ( ∇ a u b + ∇ b u a ) / can be deduced from equation (IV.8).In order to close the set of equations (IV.7), (IV.8), (IV.12), and (IV.14), we use δT = (cid:18) ∂T∂ y (cid:19) s , u δ y + (cid:18) ∂T∂ s (cid:19) y , u δ s + (cid:18) ∂T∂u ab (cid:19) s , y δu ab , (IV.15) δp = (cid:18) ∂p∂ y (cid:19) s , u δ y + (cid:18) ∂p∂ s (cid:19) y , u δ s + (cid:18) ∂p∂u ab (cid:19) s , y δu ab , (IV.16) δφ ab = (cid:18) ∂φ ab ∂ y (cid:19) s , u δ y + (cid:18) ∂φ ab ∂ s (cid:19) y , u δ s + (cid:18) ∂φ ab ∂u cd (cid:19) s , y δu cd , (IV.17)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 ∂ t y = D ab y ∇ a ∇ b y + D ab ys ∇ a ∇ b s + F abc y v ∇ a ∇ b v c + D abcd y u ∇ a ∇ b ∇ c u d (IV.18)with the following coefficients, D ab y ≡ − ξ ba T (cid:18) ∂T∂ y (cid:19) s , u − ζ ac (cid:18) ∂φ bc ∂ y (cid:19) s , u , (IV.19) D ab ys ≡ − ξ ba T (cid:18) ∂T∂ s (cid:19) y , u − ζ ac (cid:18) ∂φ bc ∂ s (cid:19) y , u , (IV.20) D abcd y u ≡ − ξ ba T (cid:18) ∂T∂u cd (cid:19) s , y − ζ ae (cid:18) ∂φ be ∂u cd (cid:19) s , y , (IV.21) F abc y v ≡ θ bca . (IV.22)Equation (IV.14) for the evolution of the entropy per unit mass takes the form ∂ t s = D ab sy ∇ a ∇ b y + D ab s ∇ a ∇ b s + F abc s v ∇ a ∇ b v c + D abcd s u ∇ a ∇ b ∇ c u d (IV.23)with D ab sy ≡ κ ab T ρ (cid:18) ∂T∂ y (cid:19) s , u + ξ ac T ρ (cid:18) ∂φ bc ∂ y (cid:19) s , u , (IV.24) D ab s ≡ κ ab T ρ (cid:18) ∂T∂ s (cid:19) y , u + ξ ac T ρ (cid:18) ∂φ bc ∂ s (cid:19) y , u , (IV.25) D abcd s u ≡ κ ab T ρ (cid:18) ∂T∂u cd (cid:19) s , y + ξ ae T ρ (cid:18) ∂φ be ∂u cd (cid:19) s , y , (IV.26) F abc s v ≡ χ abc T ρ . (IV.27)Equation (IV.8) for the displacement vector becomes ∂ t u b = v b + F abc uv ∇ a v c + D ab u y ∇ a y + D ab u s ∇ a s + D abcd u ∇ a ∇ c u d (IV.28)with D ab u y ≡ ξ ba T (cid:18) ∂T∂ y (cid:19) s , u + ζ bc (cid:18) ∂φ ac ∂ y (cid:19) s , u = − D ba y , (IV.29)5 D ab u s ≡ ξ ba T (cid:18) ∂T∂ s (cid:19) y , u + ζ bc (cid:18) ∂φ ac ∂ s (cid:19) y , u = − D ba ys , (IV.30) D abcd u ≡ ξ ba T (cid:18) ∂T∂u cd (cid:19) s , y + ζ be (cid:18) ∂φ ae ∂u cd (cid:19) s , y = − D bacd y u , (IV.31) F abc uv ≡ − θ acb = − F bac y v . (IV.32)Finally, equation (IV.7) for the velocity field gives ∂ t v b = C abcd vu ∇ a ∇ c u d + C ab v y ∇ a y + C ab v s ∇ a s + D abcd v ∇ a ∇ c v d + F abcde vu ∇ a ∇ c ∇ d u e + F abc v y ∇ a ∇ c y + F abc v s ∇ a ∇ c s (IV.33)with C abcd vu ≡ ρ (cid:18) ∂σ ab ∂u cd (cid:19) s , y = − ρ (cid:18) ∂p∂u cd (cid:19) s , y δ ab + 1 ρ (cid:18) ∂φ ab ∂u cd (cid:19) s , y , (IV.34) C ab v y ≡ ρ (cid:18) ∂σ ab ∂ y (cid:19) s , u = − ρ (cid:18) ∂p∂ y (cid:19) s , u δ ab + 1 ρ (cid:18) ∂φ ab ∂ y (cid:19) s , u , (IV.35) C ab v s ≡ ρ (cid:18) ∂σ ab ∂ s (cid:19) y , u = − ρ (cid:18) ∂p∂ s (cid:19) y , u δ ab + 1 ρ (cid:18) ∂φ ab ∂ s (cid:19) y , u , (IV.36) D abcd v ≡ η abcd ρ , (IV.37) F abcde vu ≡ − χ cab T ρ (cid:18) ∂T∂u de (cid:19) s , y + θ abf ρ (cid:18) ∂φ cf ∂u de (cid:19) s , y , (IV.38) F abc v y ≡ − χ cab T ρ (cid:18) ∂T∂ y (cid:19) s , u + θ abd ρ (cid:18) ∂φ cd ∂ y (cid:19) s , u , (IV.39) F abc v s ≡ − χ cab T ρ (cid:18) ∂T∂ s (cid:19) y , u + θ abd ρ (cid:18) ∂φ cd ∂ s (cid:19) y , u , (IV.40)where (IV.34) is the rank-four elasticity tensor divided by the mean mass density ρ = mn eq , .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 couplingproperties.We note that the evolution equations (IV.18), (IV.23), (IV.28), and (IV.33) form a closed set of linearpartial differential equations, which can thus be written in matrix form as ∂ t ψψψ = ˆ L · ψψψ with ψψψ = ( δ y , δ s , δv b , δu b ) T , (IV.41)where ˆ L is a × 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 dispersionrelations z ( q ) of the eight hydrodynamic modes are provided by solving the eigenvalue problem: L · ψψψ = z ψψψ (IV.42)with the × matrix L = − D cd y q c q d − D cd ys q c q d − F cdb y v q c q d − ıD cdeb y u q c q d q e − D cd sy q c q d − D cd s q c q d − F cdb s v q c q d − ıD cdeb s u q c q d q e ı C ca v y q c − F cad v y q c q d ı C ca v s q c − F cad v s 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 ca u y q c ıD ca u s q c δ ab + ıF cab uv q c − D cadb u q c q d . (IV.43)6The eigenvalue problem can be solved perturbatively starting from the dispersion relations of the dissi-pativeless crystal with vanishing coefficients D ’s and F ’s. Such a method corresponds to expanding thedispersion relations in powers of the wave number q . For the dissipativeless crystal, the dispersion relationsare linear in the wave number q = k q k . The next order of the perturbation calculation gives the terms of O ( q ) , providing the damping rates of the modes. Since the dissipativeless and dissipative current densitiesare known at leading order in the gradient expansion of the fields, the corrections of O ( q ) are not relevantto the calculation.Accordingly, we consider the following expansions of the matrix (IV.43), its eigenvalues, and its eigenvec-tors: L = L (0) + L (1) , ψψψ = ψψψ (0) + ψψψ (1) + · · · , z = z (0) + z (1) + · · · (IV.44)with L (0) = ı C ca v y q c ı C ca v s q c − C cadb vu q c q d δ ab (IV.45)and L (1) = − D cd y q c q d − D cd ys q c q d − F cdb y v q c q d − ıD cdeb y u q c q d q e − D cd sy q c q d − D cd s q c q d − F cdb s v q c q d − ıD cdeb s u q c q d q e − F cad v y q c q d − F cad v s q c q d − D cadb v q c q d − ıF cadeb vu q c q d q e ıD ca u y q c ıD ca u s q c ıF cab uv q c − D cadb u q c q d . (IV.46)In order to solve the eigenvalue problem at zeroth order, we introduce the × real symmetric matrix M defined with the elements M ab ≡ C cadb vu q c q d (IV.47)and the three-dimensional vectors N y and N s with the components N a y ≡ C ca v y q c and N a s ≡ C ca v s q c . (IV.48)Since the matrix M is real symmetric, it can be diagonalized by an orthogonal transformation giving the threeeigenvalues λ σ with σ = 1 , , , and the three eigenvectors ǫǫǫ σ , which are supposed to form an orthonormalbasis, ǫǫǫ σ · ǫǫǫ σ ′ = δ σσ ′ [1]: M · ǫǫǫ σ = λ σ ǫǫǫ σ for σ = 1 , , . (IV.49)Since M = O ( q ) , its eigenvalues are also going as λ σ = O ( q ) . The eigenvectors ǫǫǫ σ also depend on the wavevector q , but since they are normalized to the unit value, they only depend on the direction of the wavevector, q /q .At zeroth order, the right eigenvectors ψψψ (0) α and the left eigenvectors ˜ ψ ˜ ψ ˜ ψ (0) α associated with the eigenvalues z (0) α are obtained by solving the following eigenvalue problem, L (0) · ψψψ (0) α = z (0) α ψψψ (0) α , L (0) † · ˜ ψ ˜ ψ ˜ ψ (0) α = z (0) ∗ α ˜ ψ ˜ ψ ˜ ψ (0) α , (IV.50)and they are taken to satisfy the biorthonormality conditions, ˜ ψ ˜ ψ ˜ ψ (0) † α · ψψψ (0) β = δ αβ for α, β = 1 , , . . . , . Theyare given by z (0) α = − ıλ / α , ψψψ (0) α = − ıλ / α ǫǫǫ α ǫǫǫ α , ˜ ψ ˜ ψ ˜ ψ (0) α = 12 ıλ − α N y · ǫǫǫ α ıλ − α N s · ǫǫǫ α − ıλ − / α ǫǫǫ α ǫǫǫ α ( α = 1 , , 3) ; (IV.51) z (0) α +3 = + ıλ / α , ψψψ (0) α +3 = ıλ / α ǫǫǫ α ǫǫǫ α , ˜ ψ ˜ ψ ˜ ψ (0) α +3 = 12 ıλ − α N y · ǫǫǫ α ıλ − α N s · ǫǫǫ α + ıλ − / α ǫǫǫ α ǫǫǫ α ( α = 1 , , 3) ; (IV.52)7 z (0)7 = 0 , ψψψ (0)7 = ı M − · N s , ˜ ψ ˜ ψ ˜ ψ (0)7 = ; (IV.53) z (0)8 = 0 , ψψψ (0)8 = ı M − · N y , ˜ ψ ˜ ψ ˜ ψ (0)8 = . (IV.54)Accordingly, the modes α = 1 - 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 dampingof the modes at this approximation.By standard perturbation theory, the first-order correction to the eigenvalues can be obtained with z (1) α = ˜ ψ ˜ ψ ˜ ψ (0) † α · L (1) · ψψψ (0) α ˜ ψ ˜ ψ ˜ ψ (0) † α · ψψψ (0) α , (IV.55)giving z (1) α = − h − λ − / α (cid:0) C ca v y q c ǫ aα (cid:1) (cid:0) F cdb y v q c q d ǫ bα (cid:1) + λ − α (cid:0) C ca v y q c ǫ aα (cid:1) (cid:0) D cdeb y u q c q d q e ǫ bα (cid:1) − λ − / α ( C ca v s q c ǫ aα ) (cid:0) F cdb s v q c q d ǫ bα (cid:1) + λ − α ( C ca v s q c ǫ aα ) (cid:0) D cdeb s u q c q d q e ǫ bα (cid:1) + D cadb v q c q d ǫ aα ǫ bα − λ − / α F cadeb vu q c q d q e ǫ aα ǫ bα − λ / α F cab uv q c ǫ aα ǫ bα + D cadb u q c q d ǫ aα ǫ bα i ( α = 1 , , 3) ; (IV.56) z (1) α +3 = − h λ − / α (cid:0) C ca v y q c ǫ aα (cid:1) (cid:0) F cdb y v q c q d ǫ bα (cid:1) + λ − α (cid:0) C ca v y q c ǫ aα (cid:1) (cid:0) D cdeb y u q c q d q e ǫ bα (cid:1) + λ − / α ( C ca v s q c ǫ aα ) (cid:0) F cdb s v q c q d ǫ bα (cid:1) + λ − α ( C ca v s q c ǫ aα ) (cid:0) D cdeb s u q c q d q e ǫ bα (cid:1) + D cadb v q c q d ǫ aα ǫ bα + λ − / α F cadeb vu q c q d q e ǫ aα ǫ bα + λ / α F cab uv q c ǫ aα ǫ bα + D cadb u q c q d ǫ aα ǫ bα i ( α = 1 , , 3) ; (IV.57) z (1)7 = − D cd s q c q d + ( C ca v s q c ) (cid:0) M − (cid:1) ab (cid:0) D cdeb s u q c q d q e (cid:1) ; (IV.58) z (1)8 = − D cd y q c q d + (cid:0) C ca v y q c (cid:1) (cid:0) M − (cid:1) ab (cid:0) D cdeb y u q c q d q e (cid:1) . (IV.59)We note that, since λ α = O ( q ) , ǫǫǫ α = O ( q ) , and M = O ( q ) , all the terms of these corrections are going as z (1) α = O ( q ) for α = 1 - . Consequently, the corrections z (1) α and z (1) α +3 with α = 1 , , are giving the dampingrates of the propagative sound modes, while the heat mode and the vacancy diffusion modes are diffusivesince z = z (1)7 = O ( q ) and z = z (1)8 = O ( q ) . Although the perturbative calculation can be continued tofurther corrections of O ( q ) , they are not relevant since the statistical-mechanical theory is limited to theleading terms in the gradient expansion, thus, neglecting the corrections of O ( ∇ ) = O ( q ) .In the case where the coefficients F ’s are vanishing, we note that z (1) α = z (1) α +3 , so that the dampingrate is the same for sound modes propagating in opposite directions. However, in crystals of the samecrystallographic classes as piezoelectric crystals where the coefficients F ’s may be non-vanishing [53], theseresults show that the degeneracy between these damping rates may be split for sound modes propagating inopposite directions, because the terms with the coefficients F ’s have opposite signs in z (1) α and z (1) α +3 . V. CUBIC CRYSTALS Now, the dispersion relations of the eight hydrodynamic modes are analyzed in detail in the case of cubiccrystals. The consequences of the cross effects coupling the transport of momentum to those of heat andcrystalline order are investigated.8 A. Tensors in cubic crystals In the case of cubic crystals, symmetric rank-two tensors are proportional to the unit tensor and symmetricrank-four tensors such as the elasticity and viscosity tensors can be expanded as (cid:18) ∂σ ab ∂u cd (cid:19) s , y = ( c + 2 c )Υ ab Υ cd + ( c − c ) (cid:0) Υ ab Υ cd + Υ ab Υ cd (cid:1) + 2 c (cid:0) Υ ab Υ cd + Υ ab Υ cd + Υ ab Υ cd (cid:1) , (V.1) η abcd = ( η + 2 η )Υ ab Υ cd + ( η − η ) (cid:0) Υ ab Υ cd + Υ ab Υ cd (cid:1) + 2 η (cid:0) Υ ab Υ cd + Υ ab Υ cd + Υ ab Υ cd (cid:1) , (V.2)with [9] Υ ab = √ δ ab , Υ ab = √ (cid:18) δ a δ b − δ ab (cid:19) , Υ ab = √ (cid:0) δ a δ b − δ a δ b (cid:1) , (V.3) Υ ab = √ (cid:0) δ a δ b + δ a δ b (cid:1) , Υ ab = √ (cid:0) δ a δ b + δ a δ b (cid:1) , Υ ab = √ (cid:0) δ a δ b + δ a δ b (cid:1) , (V.4)the three isoentropic elastic constants c , c , and c , and the three viscosity coefficients η , η , and η ,expressed in Voigt’s notations. As a consequence, we have in particular that M = ( C cadb vu q c q d ) = C q x + C ( q y + q z ) ( C + C ) q x q y ( C + C ) q x q z ( C + C ) q x q y C q y + C ( q x + q z ) ( C + C ) q y q z ( C + C ) q x q z ( C + C ) q y q z C q z + C ( q x + q y ) (V.5)with C ≡ c ρ , C ≡ c ρ , C ≡ c ρ , (V.6)and q x ≡ q , q y ≡ q , and q z ≡ q . Similar decompositions hold for other symmetric rank-four tensors. Inparticular, the tensor (IV.37) can be expressed in Voigt’s notations with the three diffusivities D v ≡ η ρ , D v ≡ η ρ , D v ≡ η ρ , (V.7)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 thepoint groups T d and T can accommodate piezoelectricity [53]. The key is that piezoelectricity is describedby a rank-three tensor, which is only non-vanishing for the classes T d and T among cubic crystals. In suchcases, the tensors χ abc and θ abc are specified by a single non-vanishing coefficient because χ xyz = χ yzx = χ zxy = χ xzy = χ zyx = χ yxz ≡ χ , (V.8) θ xyz = θ yzx = θ zxy = θ xzy = θ zyx = θ yxz ≡ θ , (V.9)while the other coefficients are equal to zero [53]. These rank-three tensors can thus be expressed as χ abc = χ Ξ abc and θ abc = θ Ξ abc (V.10)with Ξ abc ≡ δ a δ b δ c + δ a δ b δ c + δ a δ b δ c + δ a δ b δ c + δ a δ b δ c + δ a δ b δ c . (V.11)For cubic crystals in the non-piezoelectric classes O h , O, and T h , these tensors are vanishing because χ = 0 and θ = 0 .9 B. Hydrodynamic modes in cubic crystals For cubic crystals, the matrix (IV.43) thus becomes L = − D y q − D ys q − F cdb y v q c q d − ıD cdeb y u q c q d q e − D sy q − D s q − F cdb s v q c q d − ıD cdeb s u q c q d q e ı C v y q a − F cad v y q c q d ı C v s q a − F cad v s 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 u y q a ıD u s q a δ ab + ıF cab uv q c − D cadb u q c q d (V.12)with the following coefficients, C v y ≡ − ρ − p y + G v y , (V.13) C v s ≡ − ρ − p s + G v s , (V.14) D y ≡ − ˜ ξ T y − ˜ ζ G v y = − D u y , (V.15) D ys ≡ − ˜ ξ T s − ˜ ζ G v s = − D u s , (V.16) D sy ≡ ˜ κ T y + ˜ ξ G v y , (V.17) D s ≡ ˜ κ T s + ˜ ξ G v s , (V.18)rank-four tensors, C abcd vu ≡ − ρ − p u δ ab δ cd + G abcd vu , (V.19) D abcd v ≡ ρ − η abcd , (V.20) D abcd u ≡ ˜ ξ T u δ ab δ cd + ˜ ζ G abcd vu = − D abcd y u , (V.21) D abcd s u ≡ ˜ κ T u δ ab δ cd + ˜ ξ G abcd vu , (V.22)and odd-order tensors, F abc y v ≡ θ Ξ abc = − F bac uv , (V.23) F abc s v ≡ ˜ χ Ξ abc , (V.24) F abc v y ≡ ( − ˜ χ T y + θ G v y ) Ξ abc , (V.25) F abc v s ≡ ( − ˜ χ T s + θ G v s ) Ξ abc , (V.26) F abcef vu ≡ − ˜ χ T u Ξ abc δ ef + θ Ξ abd G cdef vu , (V.27)where p y ≡ (cid:18) ∂p∂ y (cid:19) s , u , p s ≡ (cid:18) ∂p∂ s (cid:19) y , u , p u ≡ 13 tr (cid:18) ∂p∂ u (cid:19) s , y , (V.28) T y ≡ (cid:18) ∂T∂ y (cid:19) s , u , T s ≡ (cid:18) ∂T∂ s (cid:19) y , u , T u ≡ 13 tr (cid:18) ∂T∂ u (cid:19) s , y , (V.29) G v y ≡ ρ tr (cid:18) ∂φφφ∂ y (cid:19) s , u , G v s ≡ ρ tr (cid:18) ∂φφφ∂ s (cid:19) y , u , G abcd vu ≡ ρ (cid:18) ∂φ ab ∂u cd (cid:19) s , y , (V.30) ˜ ξ ≡ ξT , ˜ ζ ≡ ρ ζ , ˜ κ ≡ κT ρ , ˜ χ ≡ χT ρ . (V.31)We note that p u in equation (V.28) can be expressed in terms of the isoentropic compressibility κ s by p u = − /κ s . Since the compressibility is related to the elastic constants c and c according to κ − s =( c + 2 c ) / [54], we have that − ρ − p u = 1 ρκ s = 13 ( C + 2 C ) , (V.32)because of equation (V.6). In this regard, the rank-four tensor G abcd v , u is defined with the two constants G = 23 ( C − C ) = − G and G = C , (V.33)0so that we have C = 1 ρκ s + G , (V.34) C = 1 ρκ s − G , (V.35) C = G , (V.36)as established in the thermodynamics of crystals [54]. The equilibrium thermodynamic stability of thecrystals requires the positivity of the compressibility and the specific heat: κ s = − p u > , c y , u ≡ T (cid:18) ∂ s ∂T (cid:19) y , u = TT s > , (V.37)whereupon p u < and T s > . Furthermore, the Grüneisen parameter γ ≡ ρT (cid:18) ∂T∂ρ (cid:19) s , y , (V.38)which is proportional to the thermal expansivity at constant pressure, is observed to be positive [54]. Underthe same circumstances, the quantity T u defined in equation (V.29) is thus negative, T u = − ρ (cid:18) ∂T∂ρ (cid:19) s , y = − T γ < . (V.39)Vacancies are also known as Schottky defects, in which atoms move to the crystal surface [1–3]. Thenumber of vacancies can be evaluated as N v ≃ N exp[ − β ( ε + pv )] , where N is the number of ions in thecrystals, ε ≃ eV the binding energy of the atoms to their lattice site, v the volume per ion in the lattice,and p the pressure. In such a case, N v ≃ − N at T = 1000 K. Since ε ≫ pv , we have that T y = (cid:18) ∂ y ∂T (cid:19) s , u ≃ εk B T y > and p y = (cid:18) ∂ y ∂p (cid:19) s , u ≃ − v k B T y < , (V.40)so that T y > and p y < . Moreover, the diffusion coefficient of vacancies is known to behave as D y ≃ D y exp( − βE v ) with D y ≃ . cm /s and E v ≃ . eV for the diffusion of Cu in Cu [2]. Therefore, vacancydiffusion coefficients are of the order of magnitude D y ∼ − m /s at T = 1000 K. The vacancy diffusioncoefficient is always positive. Since it is related by D y = − ˜ ξ T y − ˜ ζ G v y to T y > and since the first term isexpected to be larger than the second term, we should have that ˜ ξ < . The heat conductivity is also alwayspositive κ > .On the one hand, experimental data are available about the elastic constants c , c , and c [1–3]. Thesound velocities are typically of the order of a few thousand meters per second in crystals, in agreement withthe values of the corresponding elastic constants and mass density. On the other hand, the experimentalmeasurements of sound attenuation in solids are also available and they show that the diffusivities associatedwith the sound damping rates and due in particular to the viscosities are of the same order of magnitude asin liquids, i.e., D v ∼ − m /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 modesin cubic crystals. The × 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 resultsreported here above. The cases with vanishing and non-vanishing coefficients F ’s are compared. A first setof parameter values is chosen to plot the dispersion relations of the eight modes for vanishing coefficients F ’s: C = 30 , C = 20 , C = 15 ,G v y = 0 . , G v s = − ,p y = − . , p s = 0 . ,T u = − , T y = 1 . , T s = 1 ,D v = 0 . , D v = 0 . , D v = 0 . , ˜ ξ = − . , ˜ ζ = 0 . , ˜ κ = 0 . , (V.41)with ˜ χ = 0 , θ = 0 . (V.42)The units are nm, ps, and K. A second set of parameters is chosen for non-vanishing coefficients F ’s,using (V.41) with ˜ χ = 0 . , θ = 0 . . (V.43) T ( – ) q | I m z | R e z LT , T ( – )( – ) H, V ( – ) LVH ( – ) LT , T ( – )( – ) | I m z | R e z q T ( – ) H, V ( – ) LVH ( – ) L T ( – ) T ( – ) (a) (b) [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, , ; (b) in the direction ( q x , q y , q z ) = ( q, q, q ) / √ .The eight dispersion relations for the parameter values (V.41) and (V.42) are plotted in figure V.1. Themodes are the two longitudinal sound modes L ( ± ) with the largest speeds, the four transverse sound modesT ( ± )1 and T ( ± )2 , the heat mode H, and the vacancy diffusion mode V. The heat mode and the vacancydiffusion modes are always diffusive. The mode with Im z α < are labelled with (+) and those with Im z α > with ( − ) in reference to the direction of propagation with respect to that of the wave vector q . Thedispersion 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, , , but they are distinct in the direction [123] of wave vector ( q x , q y , q z ) = ( q, q, q ) / √ ,which is well known [1–3].Now, for a cubic crystal in the crystallographic classes of piezoelectricity, the eight dispersion relations areplotted in figure V.2 for the parameter values (V.41) and (V.43). The dispersion relations are the same asin figure V.1(a) in the direction [100] of wave vector ( q x , q y , q z ) = ( q, , . However, they differ from thoseof figure V.1(b) in the direction [123] of wave vector ( q x , q y , q z ) = ( q, q, q ) / √ . As seen in figure V.2(b),the damping rates of the propagative sound modes are different for opposite propagation directions. Theresults confirm the expectation from the perturbation calculation of the dispersion relations, showing that thenon-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 thecrystallographic classes (III.44)-(III.45).2 q | I m z | R e z LT , T ( – )( – ) H, V ( – ) LVH ( – ) LT , T ( – )( – ) | I m z | R e z q T ( – ) H, V ( – ) LVH ( - ) L T ( - ) T ( – ) T ( + ) (a) (b) [100] [123] ( + ) L T ( + ) T ( - ) Figure V.2: Dispersion relations of the eight hydrodynamic modes in a crystal with the parameters (V.41)and (V.43): (a) in the direction ( q x , q y , q z ) = ( q, , ; (b) in the direction ( q x , q y , q z ) = ( q, q, q ) / √ . VI. CONCLUSION This paper applies the local equilibrium approach of reference [41] to crystals. Since the three-dimensionalgroup of space translations is broken into one of the 230 crystallographic space groups in these phases, threeNambu-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 con-structed on the basis of the mechanism of continuous symmetry breaking and the result previously obtainedin references [37, 38] for cubic crystals is recovered. The local equilibrium approach to the nonequilibriumstatistical mechanics of crystals provides the microscopic expressions for the thermodynamic and transportproperties 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. Theheat conductivities form a rank-two tensor, as well as the friction coefficients of crystalline order and thecoefficients of coupling between heat transport and crystalline order. The viscosity coefficients form a rank-four tensor. Moreover, in the same 20 crystallographic classes as those compatible with piezoelectricity, theremay exist coefficients coupling the transport of momentum to those of either heat or crystalline order. Thesecross effects are described by rank-three tensors. In the 12 other crystallographic phases, these rank-threetensors are vanishing.The hydrodynamics of crystals around equilibrium is investigated by linearizing the eight macroscopicequations ruling the local conservation laws of mass, energy, and momentum, together with the evolutionequations for the three components of the displacement vector. The eight hydrodynamic modes and theirdispersion relation are obtained by solving the corresponding eigenvalue problem and by using an expansionin powers of the wave number starting from the dissipativeless crystal as a reference. The dispersion relationsare analyzed in detail for cubic crystals.In this way, the damping rates of the hydrodynamic modes are calculated at second order in the wavenumber. The cross effects existing in the 20 crystallographic classes of piezoelectricity and coupling mo-mentum to heat or crystalline order are shown to split the degeneracy of the damping rates for the soundmodes propagating in opposite generic directions. In these crystallographic classes, the damping of soundwaves may thus be stronger in some directions than in the opposite one, as a consequence of this degeneracysplitting. 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 Gaussian3white noise fields δ J ac α ( r , t ) to the dissipative current densities J ac α ( r , t ) . These Gaussian white noise fieldsshould have amplitudes given by the linear response coefficients in terms of the transport coefficients accord-ing to the fluctuation-dissipation theorem [58, 59]. In this way, the effects of fluctuations on the transportproperties 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 microscopicmotion of atoms, but they can be generalized to quantum mechanics. In particular, the classical Green-Kuboformulas can be modified into quantum-mechanical formulas by introducing the Hermitian operators for thedensities and current densities [10, 16–18, 25] and by taking into account their non-commutativity usingimaginary-time integrals up to the inverse temperature [32]. With such considerations, the results can beextended to the dynamics of crystals in quantum regimes. Acknowledgements Financial support from the Université Libre de Bruxelles (ULB) and the Fonds de la Recherche Scientifique- FNRS under the Grant PDR T.0094.16 for the project "SYMSTATPHYS" is acknowledged. Appendix A: Local equilibrium mean decay rate of the order parameter If the crystal is close to equilibrium, the mean particle density n eq ( r ) has the lattice periodicity, so thatit can be expanded in lattice Fourier modes as in equation (II.43). As a consequence, the tensor (II.40) canbe expressed as N ab = X G G a G b | n eq , G | . (A.1)Moreover, we have that X G G a | n eq , G | = 0 , (A.2)because the following integral over a unit cell v of the lattice is vanishing by the periodicity of the density Z v n eq ∇ a n eq d r = 0 . (A.3)Now, since the local equilibrium mean values of the momentum density is given by h ˆ g c i leq = mn eq v c , thedecay rate (II.47) of the crystalline order parameter has the following local equilibrium mean value ¯ J u a ( r ) = h ˆ J u a ( r ; Γ) i leq ,λλλ = − ( NNN − ) ab Z d r ′ ∆( r − r ′ ) ∇ ′ b n eq ( r ′ ) ∇ ′ c [ n eq ( r ′ ) v c ( r ′ )]= − ( NNN − ) ab Z d r ′ ∆( r − r ′ ) ∇ ′ b n eq ( r ′ ) (cid:2) ∇ ′ c n eq ( r ′ ) v c ( r ′ ) + n eq ( r ′ ) ∇ ′ c v c ( r ′ ) (cid:3) ≡ h ˆ J u a i Aleq ,λλλ + h ˆ J u a i Bleq ,λλλ . (A.4)These two terms are calculated separately using the definition (II.29) of the function ∆( r − r ′ ) , the expan-sion (A.1) of the periodic crystal density, and the assumption that the velocity field is a macrofield such thatits Fourier modes belong to the first Brillouin zone B .On the one hand, the first term in (A.4) is obtained as follows, h ˆ J u a i Aleq ,λλλ = − ( NNN − ) ab Z d r ′ ∆( r − r ′ ) ∇ ′ b n eq ( r ′ ) ∇ ′ c n eq ( r ′ ) v c ( r ′ ) (A.5) = ( NNN − ) ab Z d r ′ Z B d k (2 π ) e ı k · ( r − r ′ ) X G , G ′ G b n eq , G e ı G · r ′ G ′ c n eq , G ′ e ı G ′ · r ′ Z B d k ′ (2 π ) e ı k ′ · r ′ ˜ v c ( k ′ )= ( NNN − ) ab X G , G ′ G b G ′ c n eq , G n eq , G ′ Z B d k (2 π ) Z B d k ′ (2 π ) e ı k · r ˜ v c ( k ′ ) (2 π ) δ ( G + G ′ + k ′ − k ) . k and k ′ are both restricted to the first Brillouin zone B , the Diracdelta distribution factorizes as δ ( G + G ′ + k ′ − k ) = δ ( k − k ′ ) δ G , − G ′ . Moreover, because of n eq , − G = n ∗ eq , G and (A.1), we get h ˆ J u a i Aleq ,λλλ = − ( NNN − ) ab X G G b G c | n eq , G | Z B d k (2 π ) e ı k · r ˜ v c ( k ) = − ( NNN − ) ab N bc v c ( r ) = − v a ( r ) , (A.6)since ( NNN − ) ab N bc = δ ac .On the other hand, the second term in (A.4) is vanishing, since h ˆ J u a i Bleq ,λλλ = − ( NNN − ) ab Z d r ′ ∆( r − r ′ ) ∇ ′ b n eq ( r ′ ) n eq ( r ′ ) ∇ ′ c v c ( r ′ )= ( NNN − ) ab Z d r ′ Z B d k (2 π ) e ı k · ( r − r ′ ) X G , G ′ G b n eq , G e ı G · r ′ n eq , G ′ e ı G ′ · r ′ Z B d k ′ (2 π ) e ı k ′ · r ′ k ′ c ˜ v c ( k ′ )= ( NNN − ) ab X G , G ′ G b n eq , G n eq , G ′ Z B d k (2 π ) Z B d k ′ (2 π ) e ı k · r k ′ c ˜ v c ( k ′ ) (2 π ) δ ( G + G ′ + k ′ − k )= ( NNN − ) ab X G G b | n eq , G | Z B d k (2 π ) e ı k · r k c ˜ v c ( k )= − ı ( NNN − ) ab X G G b | n eq , G | ∇ c v c ( r ) = 0 , (A.7)because of (A.2). Replacing the results (A.6) and (A.7) into (A.4), equation (III.17) for the local equilibriummean value of the decay rate is thus proved. [1] N. W. Ashcroft and N. D. Mermin. Solid State Physics . HRW International Editions, Philadelphia, 1976.[2] C. Kittel. Introduction to Solid State Physics . Wiley, New York, 1976.[3] H. Ibach and H. Lüth. Solid-State Physics , 4th edition. Springer, Berlin, 2009.[4] Y. Nambu. Quasiparticles and gauge invariance in the theory of superconductivity. Phys. Rev. , 117:648, 1960.[5] J. Goldstone. Field theories with superconductor solutions. Il Nuovo Cimento (1955-1965) , 19:154, 1961.[6] H. Stern. Broken symmetry, sum rules, and collective modes in many-body systems. Phys. Rev. , 147:94, 1966.[7] P. W. Anderson. Basic Notions of Condensed Matter Physics . W. A. Benjamin, Advanced Book Program, MenloPark CA, 1984.[8] P. C. Martin, O. Parodi, and P. S. Pershan. Unified hydrodynamic theory for crystals, liquid crystals, and normalfluids. Phys. Rev. A , 6:2401, 1972.[9] P. D. Fleming and C. Cohen. Hydrodynamics of solids. Phys. Rev. B , 13:500, 1976.[10] H. Mori. Statistical-mechanical theory of transport in fluids. Phys. Rev. , 112:1829, 1958.[11] L. P. Kadanoff and P. C. Martin. Hydrodynamic equations and correlation functions. Ann. Phys. , 24:419, 1963.[12] J. A. McLennan Jr. Statistical mechanics of transport in fluids. Phys. Fluids , 3:493, 1960.[13] J. A. McLennan Jr. Nonlinear effects in transport theory. Phys. Fluids , 4:1319, 1961.[14] J. A. McLennan Jr. The formal statistical theory of transport processes. Adv. Chem. Phys. , 5:261, 1963.[15] G. P. DeVault and J. A. McLennan Jr. Statistical mechanics of viscoelasticity. Phys. Rev. , 137:724, 1965.[16] D. N. Zubarev. A statistical operator for non stationary processes. Sov. Phys. Doklady , 10:850, 1966.[17] B. Robertson. Equations of motion in nonequilibrium statistical mechanics. Phys. Rev. , 144:151, 1966.[18] B. Robertson. Equations of motion in nonequilibrium statistical mechanics. II. Energy transport. Phys. Rev. ,160:175, 1967.[19] R. A. Piccirelli. Theory of the dynamics of simple fluids for large spatial gradients and long memory. Phys. Rev. ,175:77, 1968.[20] R. C. Desai and R. Kapral. Translational hydrodynamics and light scattering from molecular fluids. Phys. Rev.A , 6:2377, 1972.[21] I. Oppenheim and R. D. Levine. Nonlinear transport processes: Hydrodynamics. Physica A , 99:383, 1979.[22] J. J. Brey, R. Zwanzig, and J. R. Dorfman. Nonlinear transport equations in statistical mechanics. Physica A ,109:425, 1981.[23] T. A. Kavassalis and I. Oppenheim. Derivation of the nonlinear hydrodynamic equations using multi-modetechniques. Physica A , 148:521, 1988. [24] J. P. Boon and S. Yip. Molecular Hydrodynamics . McGraw-Hill, New York, 1980.[25] A. I. Akhiezer and S. V. Peletminskii. Methods of statistical physics . Pergamon Press, Oxford, 1st edition, 1981.Translated by M. Schukin.[26] H. Spohn. Large Scale Dynamics of Interacting Particles . Springer, Berlin, 1991.[27] S.-i. Sasa. Derivation of hydrodynamics from the Hamiltonian description of particle systems. Phys. Rev. Lett. ,112:100602, 2014.[28] R. Zwanzig. Memory effects in irreversible thermodynamics. Phys. Rev. , 124:983, 1961.[29] H. Mori. Transport, collective motion, and Brownian motion. Prog. Theor. Phys. , 33:423, 1965.[30] M. S. Green. Markoff random processes and the statistical mechanics of time-dependent phenomena. J. Chem.Phys. , 20:1281, 1952.[31] M. S. Green. Markoff random processes and the statistical mechanics of time-dependent phenomena. II. Irre-versible processes in fluids. J. Chem. Phys. , 22:398, 1954.[32] R. Kubo. Statistical mechanical theory of irreversible processes. I. General theory and simple applications inmagnetic and conduction problems. J. Phys. Soc. Jpn. , 12:570, 1957.[33] E. Helfand. Transport coefficients from dissipation in a canonical ensemble. Phys. Rev. , 119:1, 1960.[34] D. Forster. Hydrodynamics and correlation functions in ordered systems: Nematic liquid crystals. Ann. Phys. ,84:505, 1974.[35] D. Forster. Hydrodynamic fluctuations, broken symmetry, and correlation functions . W. A. Benjamin, AdvancedBook Program, Reading MA, 1975.[36] P. M. Chaikin and T. C. Lubensky. Principles of Condensed Matter Physics . Cambridge University Press,Cambridge UK, 1995.[37] G. Szamel and M. H. Ernst. Slow modes in crystals: A method to study elastic constants. Phys. Rev. B , 48:112,1993.[38] G. Szamel. Statistical mechanics of dissipative transport in crystals. J. Stat. Phys. , 87:1067, 1997.[39] C. Walz and M. Fuchs. Displacement field and elastic constants in nonideal crystals. Phys. Rev. B , 81:134110,2010.[40] J. M. Häring, C. Walz, G. Szamel, and M. Fuchs. Coarse-grained density and compressibility of nonideal crystals:General theory and an application to cluster crystals. Phys. Rev. B , 92:184103, 2015.[41] J. Mabillard and P. Gaspard. Microscopic approach to the macrodynamics of matter with broken symmetries. J. Stat. Mech. , 2020:103203, 2020.[42] P. Curie. On symmetry in physical phenomena, symmetry of an electric field and of a magnetic field. J. Phys.Théor. Appl. Phys. Rev. , 37:405, 1931.[44] L. Onsager. Reciprocal relations in irreversible processes II. Phys. Rev. , 38:2265, 1931.[45] H. B. G. Casimir. On Onsager’s principle of microscopic reversibility. Rev. Mod. Phys. , 17:343, 1945.[46] D. Frenkel and B. Smit. Understanding Molecular Simulation . Academic Press, San Diego, 2nd edition, 2002.[47] M. P. Allen and D. J. Tildesley. Computer Simulation of Liquids . Oxford University Press, Oxford, 2nd edition,2017.[48] I. Prigogine. Introduction to Thermodynamics of Irreversible Processes . Wiley, New York, 1967.[49] S. R. de Groot and P. Mazur. Nonequilibrium Thermodynamics . Dover, New York, 1984.[50] R. Haase. Thermodynamics of Irreversible Processes . Dover, New York, 1969.[51] G. Nicolis. Irreversible thermodynamics. Rep. Prog. Phys. , 42:225, 1979.[52] H. B. Callen. Thermodynamics and an Introduction to Thermostatistics . Wiley, New York, 2nd edition, 1985.[53] L. D. Landau and E. M. Lifshitz. Electrodynamics of Continuous Media . Pergamon Press, Oxford, 2nd edition,1984.[54] D. C. Wallace. Thermodynamics of Crystals . Dover, Mineola NY, 1998.[55] W. P. Mason. Acoustic Properties of Solids , section 3, pages 98-117. in: D. E. Gray, Editor, AIP Handbook , 3rdedition. McGraw-Hill, New York, 1972.[56] D. R. Lide, Editor. CRC Handbook of Chemistry and Physics , 81st edition. CRC Press, Boca Raton FL, 2000.[57] S. Wolfram. Mathematica . Addison-Wesley Publishing Company, Redwood City CA, 1988.[58] L. D. Landau and E. M. Lifshitz. Statistical Physics, Part 2 . Pergamon Press, Oxford, 1980.[59] J. M. Ortiz de Zárate and J. V. Sengers. Hydrodynamic Fluctuations in Fluids and Fluid Mixtures . Elsevier,Amsterdam, 2006.[60] We note that the coefficient χ abc in equation (6.19) of reference [41] should correctly read χ cabcab