Microscopic approach to the macrodynamics of matter with broken symmetries

A unified set of hydrodynamic equations describing condensed phases of matter with broken continuous symmetries is derived using a generalization of the statistical-mechanical approach based on the local equilibrium distribution. The dissipativeless and dissipative parts of the current densities and the entropy production are systematically deduced in this approach by expanding in powers of the gradients of the macrofields. Green-Kubo formulas are obtained for all the transport coefficients. The results apply to both crystalline solids and liquid crystals. The consequences of microreversibility and spatial symmetries are investigated, leading to the prediction of cross effects resulting from Onsager-Casimir reciprocal relations.


I. INTRODUCTION
The spontaneous breaking of continuous symmetries is a ubiquitous phenomenon in Nature. It manifests itself in the quantum vacuum of spacetime [1][2][3], in condensed phases of matter at equilibrium [4], or in dissipative structures existing far from equilibrium [5,6]. In nonrelativistic condensed matter at equilibrium, different kinds of continuous symmetries may be broken including internal gauge symmetries in superfluids or superconductors [7][8][9][10], and spatial symmetries of translations or rotations in crystals, liquid crystals, or magnetic materials [11,12]. The breaking of continuous symmetries generates long-range order, rigidity (i.e., the possibility to drag the whole condensed phase from its boundaries), as well as Nambu-Goldstone modes [8,9] beside the hydrodynamic modes resulting from the five fundamental conservation laws for mass, energy, and linear momentum. All these soft modes have frequencies vanishing with their wave number, since they represent perturbations with respect to equilibrium becoming slower and slower as their wave length increases. These modes may be propagative (e.g., the sound modes) or diffusive (e.g., the heat mode). In any case, they are damped because of energy dissipation coming from thermal fluctuations at positive temperature.
All these effects can be described in terms of macroscopic equations ruling the time evolution of these modes. The macroscale formulation of hydrodynamics in matter characterized by broken symmetries has been achieved, in particular, for crystalline solids and liquid crystals [13,14]. A basic issue is to deduce the macroscopic equations from the underlying microscopic dynamics of atoms and molecules composing matter. For this purpose, sta-tistical mechanics is required, not only for matter at equilibrium, but also away from equilibrium to obtain the time-dependent properties including the transport coefficients associated with energy dissipation [15][16][17][18][19][20][21]. In the regime of relaxation towards global equilibrium, linear response theory combined with projection-operator method has been much developed to deduce the transport properties of condensed phases with broken symmetries [11,22]. This approach has also been formulated for crystalline solids in Refs. [23,24]. However, nonlinear effects often arise because of advection induced by the velocity of the system, as it is the case in fluid turbulence. In order to deduce the Eulerian terms in hydrodynamics, a more systematic approach consists of using local equilibrium probability distributions instead of global equilibrium distributions. This method, which has been developed since the sixties [25][26][27][28][29][30][31][32][33][34][35], provides not only the microscopic expressions for the dissipativeless fluxes of Eulerian type, but also local thermodynamics and the statistical-mechanical expression for entropy production in terms of the dissipative fluxes. Furthermore, these latter fluxes can be obtained at leading order in the gradients of the macrofields together with the transport coefficients given by Green-Kubo formulas. In this framework, the consequences of microreversibility and spatial symmetries can also be investigated. However, this method has been defined and used only for normal fluids and its generalization to a system with broken symmetry is still lacking.
Our aim is here to use this systematic approach to derive a unified set of macroscopic equations applicable to both crystalline solids and liquid crystals. The full derivation relies on the identification of the hydrodynamic variables and local order parameters associated with the broken continuous symmetries at the microscopic level of description. These variables obey balance equations in the form of local conservation laws, which can be obtained from the microscopic Hamiltonian dynamics, as presented in Sec. II. On this basis, the local equilibrium probabil-ity distribution is introduced and its time evolution can be investigated using the microscopic Hamiltonian dynamics. In this way, the entropy functional can be defined, as well as the Massieu functional given by its Legendre transform. The microscopic expression for entropy production is deduced, allowing us to identify the dissipativeless and dissipative fluxes (also called current densities) as summarized in Sec. III. The calculations are carried out by expanding the macrofields in powers of their gradients around local equilibrium. The local thermodynamic relations for matter with broken symmetries are obtained at leading order in the gradients in Sec. IV. The dissipativeless current densities are determined at next order in Sec. V, including the extra contributions arising from the broken symmetries. In Sec. VI, the dissipative current densities are derived and the Green-Kubo formulas are found giving all the possible transport coefficients. Even if the microscopic expressions of the quantities of interest are essential for a concrete evaluation of the dissipative coefficients, the full derivation can be performed without their explicit knowledge. The applicability to crystalline solids is discussed in Sec. VII and the case of liquid crystals in Sec. VIII. Finally, our concluding remarks are presented in Sec. IX.

Notations and conventions:
Latin letters a, b, c, ... = x, y, z correspond to spatial coordinates and Greek letters α, β, γ, . . . label the hydrodynamic variables and the order parameters. The indices i, j, k, · · · = 1, 2, . . . , N label the particles, where N is the total number of particles. We take, for simplicity, a single species of particles. The microscopic expression of a field such as a densityĉ α or a current densityĴ a c α is denoted with a hat to make the distinction with respect to its expectation value denoted without the hat. Einstein's convention of summation over repeated indices is adopted.

A. Hamiltonian dynamics
We consider a system of N particles of mass m, all supposed of the same species. The particle i has the position r i , the velocityṙ i = dr i /dt, and the momentum p i = mṙ i . Their microscopic dynamics is ruled by the following Hamiltonian function where V is the interaction potential energy and r ij = ||r i − r j || the distance between the parti-cles i and j (with i, j = 1, 2, . . . , N ). The positions and the momenta are three-dimensional vectors with the components r i = (r a i ) and p i = (p a i ) (with a = x, y, z).
The time evolution of this system is represented in the phase space Γ = (r i , p i ) N i=1 of dimension 6N , where the trajectories Γ t = Γ t (Γ 0 ) are uniquely determined by their initial condition Γ 0 . Any observable function A(Γ) defined in phase space is thus evolving in time according is the Poisson bracket between the phase-space functions A and B. For time-independent observable functions such that ∂ t A = 0, time integration gives A t = A(Γ t ), since dΓ t /dt = {Γ t , H}. This time evolution can be expressed in terms of the Liouvillian operator LA ≡ {A, H} according to dΓ t /dt = LΓ.
The observable quantities may be global, i.e., extensive, such as the total energy given by the Hamiltonian function, or local, i.e., intensive, such as densities. In particular, the densities of energyê, masŝ ρ, and momentumĝ = (ĝ a ) are the variables associated with five fundamental conservation laws in the system. The time variations of these densities can be expressed in terms of the divergence of a current density or flux, leading to slow modes called hydrodynamic modes.
Beyond, there may exist other densities, such as the densities of kinetic energy and potential energy, which are not directly associated with conservation laws. Their time variations are ruled by a rate, instead of a flux divergence, so that they most often generate fast modes, called kinetic modes. Nevertheless, upon continuous symmetry breaking, some of these fast local quantities may turn into slow modes, called Nambu-Goldstone modes [8,9], which thus behave as hydrodynamic modes. According to the Goldstone theorem [9], there are as many such slow modes as continuous symmetries that are broken. Therefore, the total number of hydrodynamic variables is the sum of the five conservation laws and the number of broken continuous symmetries (one, two or three depending on the type of system in consideration, i.e., a liquid crystal or a crystalline solid). For the rest of this section, we introduce the microscopic definitions of the hydrodynamic variables and their current densities.

B. Conserved quantities
The Hamiltonian system ruled by Eq. (II.1) has several conserved quantities of fundamental origin. Since the Hamiltonian function is time independent ∂ t H = 0, this function is itself a conserved quantity because {H, H} = 0, which represents the total energy E = H. Besides, the total mass M = mN is conserved. Moreover, the Hamiltonian function is invariant under spatial translations r i → r i + a with a ∈ R 3 , so that the total momentum P = i p i is conserved. Because of the symmetry of the Hamiltonian function (II.1) under rotations O ∈ SO(3), the angular momentum L = i r i ×p i is also conserved.
In order to obtain the local conservation laws, we introduce the densities associated with the total mass, energy, and linear momentum according to mass density: energy density: wheren denotes the particle density. The extensive quantities are given by integrating these densities over space: M = ρ dr, E = ê dr, and P = ĝ dr. We note that a density of angular momentum could be introduced similarly. However, the density of angular momentum is not strictly local since it is defined with respect to a coordinate origin, so that its status is different from the five aforedefined densities [11]. It is known [19] that the densities (II.2)-(II.4) obey the following local conservation equations, given in terms of the following current densities or fluxes, is the force exerted on the particle i by the particle j, and is a uniform linear density distributed on the straight line joining the positions of the particles i and j [35]. As expected, there are five hydrodynamic variables coming from the conservation laws.

C. Breaking of continuous symmetries
A low enough temperature, phase transitions happen from normal fluids to liquid crystals or crystalline solids. In these new phases, continuous symmetries are broken in the structure of matter at equilibrium. In nematic liquid crystals, continuous rotational symmetry is broken by the emergence of a special orientation of the molecules, while the continuous translational symmetry is broken in crystals, where only discrete translational symmetry remains. In these phases of matter, the continuous symmetries of uniform and isotropic normal fluids are thus broken. Such phenomena are not described by Gibbsian statistical distributions based on the Hamiltonian function (II.1) since this latter is invariant under continuous translations and rotations. Therefore, an external potential energy should be added to the Hamiltonian in order to break explicitly the symmetries, The symmetric Hamiltonian function (II.1) is recovered in the limit ǫ → 0. At the inverse temperature β = (k B T ) −1 and chemical potential µ, the equilibrium phase of matter can be described by the following Gibbsian grand canonical probability distribution p eq (Γ) = 1 Ξ ∆Γ e −β(Hǫ−µM) , (II. 14) expressed in terms of the total mass M = mN and depending on the random variables Γ = (Γ N , N ) since the particle number N is also a random variable. In Eq. (II.14), Ξ is the partition function given by the normalization condition (II.15) and ∆Γ = h 3N is the elementary phase-space volume, where h is Planck's constant. In the following, Boltzmann's constant is set equal to unity, k B = 1, except if it is explicitly written. The statistical average with respect to the probability distribution is denoted A ≡ A(Γ) p(Γ) dΓ.
Because of the external potential V (ext) , continuous symmetries are explicitly broken. For instance in crystals, the particle density can now have an equilibrium mean value n eq (r) = n(r; Γ) eq that is a periodic function in space, which is invariant under the discrete crystalline group of symmetry, but no longer under the continuous group of threedimensional translations and rotations.
We note that the system can undergo a phenomenon of spontaneous symmetry breaking without the help of the external potential, i.e., for the Hamiltonian function (II.13) in the limit ǫ → 0. For ǫ = 0, all the continuous symmetries are restored for the probability distribution (II.14). However, the system manifests long-range order. For instance, in crystals, the particle density is periodic with respect to the center of mass of the crystal. Accordingly, the symmetric Hamiltonian function (II.1) should be split into the part describing the motion of the center of mass and the other part ruling the dynamics in the frame moving with the center of mass. This latter part is no longer symmetric under continuous translations and can itself define a grand canonical probability distribution with the broken symmetry. This reasoning shows that including the center of mass in the Hamiltonian function restores the continuous symmetry in Gibbsian equilibrium probability distributions. Similarly, in nematic liquid crystals, the part of the Hamiltonian ruling the rotation of the system around the orientation selected by spontaneous symmetry breaking should be separated from the rest of the Hamiltonian function in order to define equilibrium probability distributions describing the properties of the phase with broken symmetry.
Phases with broken symmetries can be characterized by local order parameters denotedx α , where the index α runs over the subset of variables originating from symmetry breaking. The decay of these variables is given by a rateĴ x α such that ∂ tx α (r; Γ t ) +Ĵ x α (r; Γ t ) = 0 . (II. 16) We may also introduce the variablesû bα ≡ ∇ bxα that obey the following equations similar to the local conservation equations (II.5)-(II.7), ∂ tû bα (r; Γ t ) + ∇ aĴ a u bα (r; Γ t ) = 0 , (II. 17) where the corresponding current density is defined asĴ The microscopic expression forx α depends on the phase in consideration, as discussed in Secs. VII and VIII for crystalline solids and liquid crystals. However, it is possible to proceed with a general derivation of the macroscopic equations without using any explicit microscopic expression forx α . Because of long-range order, the equilibrium correlation functions of the variablesx α decay in space as (II. 19) or for their Fourier transforms, In this regard, these order parameters have a singular behavior, contrary to regular conserved quantities that have correlations of medium or short range. Nevertheless, the gradientsû aα ≡ ∇ axα of the order parameters have medium-or short-ranged correlations because δû aα (q) δû bβ (−q) eq ∼ q 0 , (II. 22) so that they are regular as for the conserved variables.

D. Nambu-Goldstone modes
As a consequence of the emergence of long-range order, there exist Nambu-Goldstone modes behaving as the conserved modes with vanishing dispersion relations for long enough wave length [6,8,9]. The equations ruling all the conserved and kinetic modes ψ β = ψβ could be written in the following form (II.23) Let us suppose that there exists an equilibrium solution Now, we may consider small perturbations of different kinds with respect to this equilibrium solution. On the one hand, an additive perturbation may be considered giving solutions of the following form, Substituting into Eq. (II.23) and linearizing, we obtain the following evolution equations, Such a mode is decaying exponentially in time with a rate that is not vanishing for long enough wave length, because ∂F β /∂ψ γ eq defines a matrix with non-vanishing elements in general. On the other hand, another perturbation may be considered where the order parameters are locally varying in space and time as For this solution, we have that Now, substituting into Eq. (II.23) and linearizing, we find where the term with the non-vanishing coefficients ∂F β /∂ψ γ eq no longer appears. If we multiply Eq. (II.31) by ∂ψ β eq /∂x γ , sum over β, integrate over the volume V of the system to average out the local variations of the symmetry-breaking equilibrium solution, and relabel the quantities, we get In Eq. (II.32), the dots denote terms with higher spatial derivatives for x β (r, t). Since the solution (II.28) is breaking continuous symmetries, we have that (∂ψ β eq /∂x α ) = 0 and thus the matrix N N N = (N αβ ) is non vanishing and can be inverted to define the quantity V V V = N N N −1 · M M M. If we suppose that the local order parameters behave as x x x(r, t) ∼ exp(ıq·r−ıωt), the frequency ω is related to the wave vector q according to showing that the dispersion relations of these solutions are vanishing with the wave number as lim q→0 ω(q) = 0, as in the case of locally conserved quantities. The eigenvalues of the matrix V V V · q are giving the propagation speed of the modes. The mode is diffusive if its propagation speed is equal to zero. The existence of the Nambu-Goldstone modes is thus a consequence of continuous symmetry breaking. There are as many such modes as components of the vector x x x = (x α ), i.e., as continuous symmetries that are broken, which is the statement of the Goldstone theorem [6,9]. In order to investigate the effects of the Nambu-Goldstone modes, we should thus consider Eq. (II.17) on the same footing as the equations (II.5)-(II.7) for the locally conserved quantities.

A. Time evolution
On the one hand, Eqs. (II.5)-(II.7) for the locally conserved quantities, as well Eq. (II.17) for the gradients of the order parameters can all be written as are respectively the densities and the corresponding current densities or fluxes. At time t, the densities are given in terms of the Liouvillian operator L or the trajectories Γ t bŷ with similar expressions for the current densities.
On the other hand, any phase-space probability density p t (Γ) at time t is given by in terms of the initial probability density p 0 (Γ) and the reversed trajectory Γ −t going from the current phase-space point Γ back to the initial conditions Γ 0 of the trajectory. Consequently, the macroscopic densities can be obtained by taking the mean value of the time-independent densities over the timeevolved probability distribution p t (Γ) or, equivalently, the mean values of the time-dependent densities over the initial probability distribution p 0 (Γ), because of Liouville's theorem dΓ 0 = dΓ t , the expectation value with respect to p t (Γ) being denoted as · t . Similar results hold for the current densities and other fields.

B. Local equilibrium distribution
The key assumption of the formalism is the initial condition being the local equilibrium distribution where λ = (λ α ) = (λ c α ) are inhomogeneous fields conjugated to the density fieldsĉ = (ĉ α ), and the asterisk * corresponds to the integration over space The normalization condition (II.15) for the local equilibrium distribution (III.7) gives the functional The expectation value with respect to the local equilibrium distribution (III.7) is denoted by · leq,λ . In this formalism, the expectation values of the densities can be obtained by taking the functional derivative of the functional (III.9) with respect to the conjugated fields as follows, (III.10) The entropy is defined as leading for the local equilibrium distribution (III.7) to the entropy functional which is the Legendre transform of the previously introduced functional (III.9). The conjugated fields are thus given by the following functional derivatives, this relation being called the second identity in Ref. [35]. Vice versa, the Legendre transform of the entropy functional (III.12) gives back the functional (III.9). We note that the equilibrium grand canonical distribution (II.14) with ǫ = 0 is recovered if the conjugated fields λ are uniform with λ e = β, λ ρ = −βµ, λ g a = 0, λ u aα = 0, and Ω = ln Ξ.

C. Time evolution of the local equilibrium distribution
The basic idea of the formalism is that the timeevolved probability density (III.5) should remain close to the local equilibrium distribution (III.7) with the conjugated fields λ t considered at time t. In this regard, these latter should be determined by the conditions ĉ α (r; Γ) t = ĉ α (r; Γ) leq,λt ≡ c α (r, t) , (III. 14) according to which the expectation values (III.6) of the densities with respect to the probability distribution p t (Γ) are equal to their expectation values with respect to the local equilibrium distribution with the conjugated fields λ t at time t. The conditions (III.14) are thus defining the macroscopic densities c α (r, t), which are given by the functional derivatives (III.10) of the functional (III.9) with respect to the conjugated fields λ = λ t at time t. Now, we consider the time evolution of the probability density starting from the initial condition given by the local equilibrium distribution (III.7) with λ = λ 0 , Since the normalization of this probability distribution should be preserved during the time evolution, we must have that (d/dt) p t (Γ) dΓ = 0. The calculation using Eq. (III.1) leads to the following relation, which should hold for any conjugated field λ 0 that may thus be replaced by λ to get This relation has been obtained notably in Refs. [31,35] and is called the first identity in Ref. [35].
Remarkably, we have that with the quantity as shown in Refs. [25,35]. Therefore, the expectation value of any observable A(Γ) with respect to the time-evolved probability distribution p t (Γ) is thus given in terms of the expectation value with respect to the local equilibrium distribution according to which is called the third identity in Ref. [35]. In particular, the conditions (III.14) are equivalent to the following relations,

D. Entropy production and dissipative current densities
The identity (III. 19) is a universal relation, which is reminiscent of the nonequilibrium work and integral fluctuation theorems [35][36][37][38][39][40][41][42][43]. Choosing A(Γ) = e −Σt(Γ) in this third identity gives [35] e −Σt(Γ) t = 1 , (III. 21) which implies by Jensen's inequality [44] that In open systems, the entropy S changes in time due to the exchanges d e S with the environment and its production d i S inside the system: dS = d e S + d i S.
Since the system is here isolated, there is no exchange with the environment d e S = 0, so that the change in time of the entropy is equal to the entropy production dS = d i S. In this regard, the result (III.22) may be interpreted as the nonnegativity of the entropy production, in agreement with the second law of thermodynamics. Using Eqs. (III.9) and (III.10), we have that According to the definition (III.12) of the entropy functional, the relation (III.14), Eq. (III.1), and integrations by parts, the entropy production is thus given by Now, using the identity (III. 19) with A taken aŝ J a c α , the expectation values of the current densities with respect to the phase-space probability distribution (III.5) can be decomposed as [25] The entropy production is thus given by where the first term vanishes because of the identity (III. 16). This term is thus expressing the conservation of entropy in adiabatic (isoentropic) processes induced by the dissipativeless current densities defined by Eqs. (III.27). The second term in Eq. (III.29) is in general non vanishing and related to the production of entropy, leading to the definition of the dissipative current densities by Eqs. (III.28). Accordingly, the entropy production can be expressed as in terms of the dissipative current densities (III.28) and the gradients of the conjugated fields, which play the role of thermodynamic forces, also called the affinities, which is in accordance with macroscopic nonequilibrium thermodynamics [45][46][47][48][49].
As we will show explicitly below, the three identities (III.16), (III.13), and (III.19) allow us to fully deduce the macroscopic equations and identify the dissipative coefficients. The local conservation equations for the mean values (III.14) can indeed be obtained from the expectation values of Eq. (III.1), giving in terms of the dissipativeless (III.27) and dissipative (III.28) current densities. In summary, the method is carried out as follows using expansions in powers of the gradients: 1. First, the slow modes of the system are identified and the local thermodynamic relations between these variables are established at leading order in the gradients.
2. Once the hydrodynamic variables are identified, the identity (III.13) is used to obtain the conjugate fields λ.
3. The dissipativeless current densities are computed by taking the expectation values of the microscopic current densities over the local equilibrium distribution, according to Eq. (III.27).
4. The dissipative current densities arise from Eq. (III.28) as a direct consequence of the third identity (III.19).
5. The Green-Kubo relations giving the linear response coefficients between the dissipative current densities and the gradients of the conjugated fields can thus be obtained.
We now proceed with the explicit computation for systems with broken symmetries, such as crystalline solids and liquid crystals.

IV. LOCAL THERMODYNAMICS
A. The local Euler, Gibbs, and Gibbs-Duhem relations The identity (III. 16) shows that the dissipativeless current densities (III.27) are leaving the entropy constant in time. All the processes involved by the dissipativeless current densities may thus be considered as reversible (i.e., adiabatic or isoentropic). The approach of nonequilibrium statistical mechanics based on the local equilibrium distribution (III.7) is therefore providing the local thermodynamic relations in every element of matter, as shown here below.
At leading order in the expansion in the gradients, the functionals (III.9) and (III.12) may be supposed of the forms defined by introducing the densities ω(λ λ λ) and s(c), which are respectively functions of the conjugated fields and mean densities. Since both functionals are interrelated by Legendre transforms, we have that giving the local relations s = λ α c α + ω , ds = λ α dc α , dω = −c α dλ α , (IV.3) up to terms of second order in the gradients. In Eq. (IV.3), the first relation can be identified as the local Euler relation, the second as the local Gibbs relation for the entropy density, and the third as the associated Gibbs-Duhem relation.
At this stage, a comparison becomes possible with previous works on the thermodynamics of matter with broken symmetry [13,14], where the relevant thermodynamic relations are given in the laboratory frame by Euler relation: e = T s + µρ + v a g a + φ aα u aα − p , (IV.4) Gibbs relation: de = T ds + µdρ + v a dg a + φ aα du aα , (IV.5) Gibbs-Duhem relation: dp = sdT + ρdµ + g a dv a + u aα dφ aα , (IV. 6) where e ≡ E/V is the mean energy density, T the temperature, s ≡ S/V the entropy density, µ the chemical potential, ρ ≡ M/V the mean mass density, v a the velocity, g a = ρv a the mean momentum density, u aα the gradients of the mean order parameters x α , φ aα the fields thermodynamically conjugated to u aα , and p the hydrostatic pressure.

B. The conjugated fields
The Euler relation (IV.4) can be written to give the entropy density instead of the energy density. After substitution into the entropy functional (III.12) and taking the functional derivatives (III.13), the conjugated fields are obtained as in terms of the local inverse temperature β = (k B T ) −1 . Furthermore, the comparison between the phenomenological Euler relation (IV.4) and the theoretical expression given by the third relation in Eq. (IV.3) allows us to identify the Legendre transform of the entropy density as the local thermodynamic potential ω = βp, proportional to the hydrostatic pressure p and referred to as a Massieu function (here, per unit volume) [49]. Now, the Gibbs relation (IV.5) written for the entropy density gives the relations holding in the laboratory frame where the center of mass of the matter element moves at the velocity v = (v x , v y , v z ). In the frame moving with the element where v = 0, the energy density and the chemical potential are given by (IV. 19) In the frame moving with the element where v = 0, we get where the subscript 0 denote the rest-frame quantities (IV.12).

V. DISSIPATIVELESS CURRENT DENSITIES
The expression (III.30) for the entropy production is showing that entropy is conserved if the dissipative parts of the current densities are vanishing. This is the case for the leading parts of the current densities giving the dissipativeless current densities (III.27). In normal fluids, these parts lead to Euler's equations. Our purpose is now to obtain their expressions in phases with broken symmetries. With this purpose, we first consider Galilean transformations from the laboratory frame to the frame moving at the velocity v of the system. In this way, the expectation values of the current densities with respect to the local equilibrium distribution (III.7) can be calculated according to the definition (III.27). In phases with broken symmetries, extra contributions are expected, which need to be determined, in particular, using the microscopic expressions of the corresponding decay ratesĴ x α introduced in Eq. (II.16). These microscopic expressions will be given in Secs. VII and VIII for crystals and liquid crystals. However, the most general form of their expectation value is known on the basis of time-reversal symmetry [13]. Once this general form is fixed, the contributions of broken symmetry to the dissipativeless current densities of momentum and energy can be obtained using the identity (III.16) expressing the conservation of entropy, which finds its origin in the adiabaticity of the reversible processes ruled by the dissipativeless current densities (III.27).

A. Galilean transformation
The particle momenta p i in the laboratory frame and those p i0 in the frame moving at the local velocity v(r) are related to each other according to p i = p i0 + mv(r i ) by the principle of Galilean relativity. Carrying out this change of variables in the microscopic expressions for the energy densityê, the momentum densityĝ a , and their current densities, we findê where the quantities with the subscript 0 are those with the momenta p i0 replacing p i , and [35] Similar expressions hold for the current densitieŝ J a u bα of the gradients of order parameters.
In Eq. (V.3), the contribution (V.5) would vanish if the velocity field v was uniform. Furthermore, we note that∆ a goes as the square of gradients, in particular, the square of the gradients of the velocity field, so that this term can be dropped if the corrections of O(∇ 2 ) are neglected. Now, the terms with odd powers of the momenta p i0 are vanishing after averaging over local equilibrium in the frame moving with the element of matter. Consequently, we get ĝ a 0 leq = 0, Ĵ a e0 leq = 0, and ρ leq = ρ. Moreover, the internal energy density is given by ê 0 leq = e 0 . The velocity field is thus defined as the ratio of the momentum to the mass densities: v a ≡ ĝ a leq / ρ leq .
As aforementioned, we need to consider the possibility that symmetry breaking may introduce extra contributions in the local equilibrium expectation values of the currents densities of energy, momentum, and the gradients u aα of the order parameter, so that we obtain the following forms for these expectation values, where e = e 0 + ρv 2 /2, the hydrostatic pressure p is separated from the contributionsJ a e BS andJ a g b BS of broken symmetries, and Since the order parametersx α are usually even under time-reversal symmetry, their rateĴ x α should be odd. Therefore, the most general form of the dissipativeless mean rates is given bȳ where the coefficients A aα and B abα are defined with the sign convention of Ref. [13]. This form will be justified on the basis of the microscopic dynamics in Secs. VII and VIII. Here, we note that B abα = 0 in crystals and A aα = 0 in nematic liquid crystals. In Eqs. (V.6)-(V.10), the terms vanishing with the velocity are due to the advection of the corresponding quantity by the motion of the element of matter at the velocity v = (v a ).
The next issue is to determine the contributions J a e BS andJ a g b BS to the dissipativeless current densities of energy and momentum by using the identity (III.16).

B. Dissipativeless current densities of energy and momentum
The dissipativeless currents are satisfying the identity (III. 16). As shown in Eq. (III.29), this identity is equivalent to the requirement that the dissipativeless terms conserve the entropy. This identity can be used to deduce the energy and momentum current densities from the expression (V.12), as follows.
Using the leading-order contributions of the conjugate fields λ α , derived in Eqs. (IV.7)-(IV.10), and writing the dissipativeless currents asJ a c α ≡ Ĵ a c α leq,λ , the identity (III.16) is giving up to higher-order corrections. We note thatJ a ρ = ρv a because of Eq. (V.7). Using the Gibbs-Duhem relation (IV.18) with the differential d replaced by the gradient ∇ a , we get Besides, the scalar product of the term ∇ a (βp) with the velocity v a can be transformed according to  16) because the integral (III.8) of the divergence ∇ a (βp v a ) is vanishing. Replacing the dissipativeless current densities by their expressions (V.8)-(V.10) with g b = ρv b given by Eq. (V.7), we find that the contributions of broken symmetry to the dissipativeless current densities should satisfy the following identity after using Eq. (V.11) with the mean decay rate of the order parameter given by Eq. (V.12). In order to satisfy this identity, we first require that ∇ a v bJ a g b BS + ∇ a φ aαJ x α can be expressed as a divergence, which leads tō Now, the divergence ∇ a (φ aα A bα v b ) resulting from the last two terms in Eq. (V.18) is combined with the factor β in front of them to give an extra contribution to the first term in ∇ a β, so that we finally obtain The dissipativeless parts of the momentum and energy current densities are thus determined as with the reversible stress tensor defined by (V.24) These results are consistent with those obtained in Refs. [13,14].
To summarize, the dissipativeless current densities are obtained from the condition that they should conserve the entropy, which is equivalent to the identity (III.16) of the formalism. Galilean invariance and the behavior of the variables under timereversal symmetry provide the identification of the most general form of the dissipativeless current densities. This approach, based on the conservation of entropy, does not require the knowledge of the explicit microscopic expressions for the variablesx α .

C. Dissipativeless local conservation equations
Now, the expressions obtained for the dissipativeless current densities can be substituted back into the local conservation equations (III.31), here neglecting the dissipative part of the current densities. We introduce the total time derivative of any field f (r, t) along a streamline of matter advected by the velocity v = (v a ) according to Expanding the Eulerian local conservation equations, we obtain their following Lagrangian forms, in terms of the reversible stress tensor (V.24) and up to higher-order corrections. These equations of motion thus describe adiabatic processes leaving constant the entropy.

D. Dissipativeless equations for the conjugated fields
Since the conjugate fields β, βµ 0 , and βφ aα are functions of the variables e 0 , ρ, and u aα , their time evolution can be obtained from the dissipativeless equations (V.26)-(V.29). As shown in detail in App. A using the Maxwell relations (IV.14)-(IV.17), they obey the following equations, up to corrections of order ∇ 2 . Furthermore, we also have which is obtained from Eq. (V.28) with the reversible stress tensor (V.24), using Eq. (IV.20) and the fact that u aα = ∇ a x α = O(∇) by definition.

A. Heat current density
Once the dissipativeless current densities are identified, we turn to the derivation of the dissipative current densities defined by Eq. (III.28) and contributing to the entropy production (III.30). This latter can be expressed in terms of the gradients of the conjugated field given by Eqs. (IV.7)-(IV.10) at leading order in the gradients. Using the fact that the dissipative current density of mass is equal to zero J a ρ = 0, the relation J a u bα = δ ab J x α , and the expansions ∇ a (βf ) = f ∇ a β + β∇ a f for any field f , we obtain in terms of the heat current density defined as This result shows that the dissipative current density of energy is reduced to the heat current density under the conditions where v a = 0 and φ aα = 0, i.e., in the frame moving with matter and in the absence of the contribution from φ aα to the stress.

B. Deduction at leading order
According to Eq. (III.28), the leading-order term in the gradient expansion of the dissipative current densities is given by Therefore, their derivation relies on the evaluation of the quantity (III.18), which can be expressed as The detailed calculations are carried out in App. B. Under the conditions v a = 0 and φ aα = 0 where the dissipative current density of energy coincides with the heat current density, we find at first order in the gradients that with the following definitions, δĴ ′a e ≡ δĴ a e − ρ −1 (e 0 + p) δĝ a , (VI.8) Substituting Eq. (VI.7) back into Eq. (VI.4), the dissipative current densities are given at leading order by Eq. (VI.3), giving up to higher-order corrections. Assuming that the characteristic length and time scales of the conjugated fields λ are much larger than the correlation length and time of the current densities (VI.8)-(VI.10), we can replace ∇ ′c λ γ (r ′ , τ ) by ∇ c λ γ (r, t) in the previous equation to find Since the conjugated fields evolve in time on a longer time scale than the correlation time of the current densities, the local equilibrium distribution at time t may be considered as the equilibrium distribution at the local values of the conjugated fields in the frame where v a = 0 and φ aα = 0, i.e., for given local values of temperature and chemical potential. Since the equilibrium distribution is stationary, we have that δâ(r, 0) δb(r ′ , τ − t) = δâ(r, t − τ ) δb(r ′ , 0) . Replacing t − τ by τ , the integral over time becomes t 0 dτ δâ(r, 0) δb(r ′ , τ − t) = t 0 dτ δâ(r, τ ) δb(r ′ , 0) , where the limit t → ∞ can be taken since the time scale t of the conjugated fields is longer than the correlation time of the current densities. Furthermore, the material properties over spatial scales larger than the microscopic structure of the phase can be defined by averaging over space. The microscopic currents are thus introduced aŝ δĴ ′a e (t) = δĴ a e (t) − ρ −1 (e 0 + p) δP a (t) , (VI.14) where δP a = δĝ a dr, δÊ = δê dr, and δM = δρ dr are the deviations in the total momentum, energy, and mass with respect to their local equilibrium value. Since the total momentum, energy, and mass are conserved, they are constants of motion [i.e., they do not fluctuate in time, but they are still random variables because the initial conditions in phase space may give different values to these constants of motion, for instance, in the grand canonical ensemble of distribution (II.14)]. They may thus be added into the equilibrium time correlation functions because δĴ a c α eq = 0. Consequently, the unprime quantities can be replaced by the prime ones in the time correlation functions. Therefore, Eq. (VI.12) becomes which holds for the dissipative current densities (J a c α ) = (J a e , J a g b , J x β ) under the conditions v a = 0 and φ aα = 0, where the dissipative energy current density is equal to the heat current density J a e = J a q .

C. Green-Kubo formulas for the transport coefficients
As a consequence, Eq. (VI.17) leads to the following expressions for the dissipative current densities of heat, momentum, and order parameters, with the transport coefficients given by the following Green-Kubo formulas, Here, the limit V → ∞ is taken at constant chemical potential, in order to define the bulk transport properties in arbitrarily large systems, removing in this way possible finite-size effects encountered in molecular dynamics simulation [50]. We note that, in isotropic phases of matter, the rank-three tensors χ abc and θ abα are vanishing, so that they are expected to be small in general.
Consequently, we have the Onsager-Casimir reciprocal relations κ ab = κ ba , η abcd = η cdab , ξ aα = ξ αa , and ζ αβ = ζ βα . This latter relation explains how the coefficient in front of ∇ b φ bα in Eq. (VI.18) is the same as the one in front of ∇ a T /T in Eq. (VI.20). Moreover, we have the Onsager-Casimir reciprocal relations χ abc = −χ bca and θ abα = −θ αab , which explains the changes of sign in front of ∇ c T /T in Eq. (VI. 19) and in front of ∇ a v b in Eq. (VI.20) with respect to the corresponding terms with the same coefficients. The coefficients χ abc generate a coupling between momentum transport and temperature gradients in a similar way as at the interface between two phases [54], which is the mechanism inducing the phenomenon of thermophoresis [55].

E. Entropy production
Because of the antisymmetry of the coefficients χ abc and θ abα , the associated terms do not contribute to entropy production and thus dissipation. In order to confirm this result, the entropy production (VI.1) is calculated from Eqs. (VI.18)-(VI.20), giving where, indeed, the terms with the coefficients χ abc and θ abα do not appear. In this regard, these terms may be considered as dissipativeless contributions to the current densities. In particular, the terms with the coefficients θ abα in Eq. (VI.20) are similar to those with the coefficients B abα in Eq. (V.12). The non-negativity of the entropy production results in particular from the conditions η abab ≥ 0, κ ab ≥ 0, ζ αα ≥ 0, and κ aa ζ αα ≥ (ξ aα ) 2 /T [47].

F. Macroscopic equations
Finally, the macroscopic equations read with the dissipativeless and dissipative current densities and rates given by Eqs.
which is exactly Eq. (6.1) of Ref. [13] in the case where χ abc = 0 and θ abα = 0. The coefficients κ ab can be interpreted as the heat conductivities and η abcd as the viscosities. Finally, note that the microscopic expressions for the order parameters and its associated current densities were not used in the derivation. However, these expressions are essential in practice and in order to evaluate the transport coefficients with the Green-Kubo formulas.

VII. CRYSTALLINE SOLIDS
The method applies in particular to crystals where the continuous symmetry under the group of spatial translations is broken into the discrete symmetry of one of the 230 crystallographic space groups. Accordingly, crystals have eight hydrodynamic modes with dispersion relations vanishing with the wave number according to the Goldstone theorem. These modes are the two longitudinal sound modes, the four transverse sound modes, the mode of heat conduction, and the mode of vacancy diffusion. These modes are damped because of energy dissipation and the present results give the Green-Kubo formulas for the coefficients ruling their damping.

A. Order parameters
For crystals, the variablesx α associated with continuous symmetry breaking are the componentsû a of the displacement field. Since the broken symmetry is a spatial symmetry in three dimensions, we can replace the Greek index α by a Latin index a. The microscopic expression for the displacement field is known [23,24,56,57]. In cubic crystals, it readŝ ∂r ′an (r ′ ; Γ) , (VII.1) given in terms of the particle densityn =ρ/m, the equilibrium particle density n eq (r), an integral over the Brillouin zone (BZ), and a normalization factor N . The equilibrium particle density is a periodic function of space, which has the symmetry of the crystallographic group. In a uniform phase where ∂n eq /∂r a = 0, the displacement field would be vanishing, as expected for order parameters. The associated rate defined by Eq. (II.16) is thus given bŷ In crystals, the microscopic strain tensor is defined as the symmetric rank-two tensor Accordingly, the associated current density introduced in Eq. (II.17) is here given bŷ According to Eq. (III.27), the dissipativeless part of the rate can be obtained by computing the expectation value of the microscopic expression (VII.2) over the local equilibrium distribution. For conjugated fields slowly varying in space over length scales much larger than the size of the lattice unit cell, we have that Ĵ u a (r; Γ) leq,λt = −v a (r, t) . (VII.5) Consequently, the comparison with Eq. (V.12) gives A ab = δ ab and B abc = 0 . (VII. 6) We note that the reversible stress tensor (V.24) is thus also symmetric. In a crystal, it is not possible to assign a particle to a lattice site due to the presence of defects, which are vacancies and interstitials. As a consequence, there is an associated phenomenon of diffusion and a corresponding mode [14,23,24,56,57]. To describe this phenomenon, the density of vacancies is defined asĉ where n eq,0 ≡ 1 v v dr n eq (r) , (VII.8) is the mean equilibrium density at equilibrium, v being the volume of the unit cell. In this way, the eight hydrodynamic modes of crystals can be obtained with the methods of Ref. [14].

C. Green-Kubo formulas
As shown in Sec. VI, the dissipative current densities are given by Eqs. (VI.18)-(VI.20) with the coefficients given by the Green-Kubo formulas (VI.21)-(VI.26). The formulas for the heat conductivities κ ab , the viscosities η abcd , the coefficients ζ ab related to vacancy diffusion, and ξ ab to the cross effect of vacancy thermal diffusion [14] are consistent with the results of Ref. [24].
Moreover, there also exist dissipativeless cross effects described by the rank-three tensors χ abc = χ acb and θ abc = θ bac . In isotropic media, such rank-three tensors are known to vanish according to Curie's principle based on space rotational symmetries. Such rank-three tensors may be non-vanishing only for 20 among the 32 crystallographic structures. These 20 crystallographic structures are the same as those selected to allow the possibility of a nonvanishing piezoelectric tensor, which is also of rank three [58]. Nevertheless, the cross effects described by the coefficients χ abc and θ abc often play a negligible role.

VIII. LIQUID CRYSTALS
Liquid crystals are composed of nonspherical molecules, which interact with different types of intermolecular forces. Rotational or translational symmetries may be broken in liquid crystals, because of the emergence of a privileged orientation, e.g., in nematics, or two-dimensional columnar order, e.g., in some phases of discotic liquid crystals [12,59].
For rotational symmetry breaking, we note that the total angular momentum can be decomposed as L = L 0 + k L k in terms of the angular momentum L 0 with respect to the origin of the laboratory frame and the angular momenta L k of the molecules k with respect to their center of mass (or any other property). Such angular momenta L k allow us to carry out local rotations in the system. The corresponding Nambu-Goldstone modes can be defined as soft modes associated with such local rotations, as discussed in Subsec. II D.

A. Order parameters
For apolar nematogens (i.e., nematic molecules), an external electric field E a (r) will explicitly break the rotation symmetry. In this case, the total external potential energy in the Hamiltonian function (II.13) is given by with the local traceless polarizability tensor q ab (r).
In general, this local order parameter can be taken as the quadrupolar contribution to the density of some property associated with the nematogenŝ where q i is the relevant property attached to the atom i ∈ k in the molecule k, and r k = (r a k ) is a position at the center of the molecule k [12].

B. Dissipativeless current densities
If the variables (VIII.2) are taken as the order parametersx α , the corresponding microscopic rates are given by Eq. (II. 16). Carrying out the change of variables p i = p i0 + m i v(r i ) where v(r) is the velocity field, the expectation values of those rates over the local equilibrium distribution give the dissipativeless parts For nematics, we have that A abc = −∇ c q ab leq = 0, since they are uniform at equilibrium. Therefore, Eq. (V.12) can be also justified for such liquid crystals on the basis of the microscopic approach.

C. Green-Kubo formulas
Again, the transport coefficients will be given by the Green-Kubo formulas (VI.21)-(VI.26). The coefficients χ abc and θ abα describing the dissipativeless cross effects may be expected to be equal to zero in most phases of liquid crystals. However, since piezoelectricity also exists in some liquid crystals [60,61], it is possible that such cross effects exist here as well, although being small.

IX. CONCLUSION
In this paper, we have shown that the macroscopic equations ruling the time evolution of matter with broken continuous symmetries can be derived in a unified microscopic approach based on the local equilibrium distribution, extending to crystalline solids and liquid crystals results previously obtained for normal fluids.
In the presence of broken continuous symmetries, the description should be extended to include the microscopic expressions of the order parameters beside the microscopic densities of mass, energy, and momentum, which are locally conserved. The time evolution of these variables is generated at the fundamental level of description by the underlying Hamiltonian microdynamics. The manifestation of spontaneous symmetry breaking is the emergence of as many Nambu-Goldstone modes as there are broken continuous symmetries. Those modes have frequencies ω(q) vanishing with their wave number q, as for the hydrodynamic modes associated with the local conservation laws. All these modes are damped because of their interaction with the other degrees of freedom. This damping is determined by the transport coefficients responsible for energy dissipation and entropy production.
Here, we have deduced these properties using the microscopic approach based on the local equilibrium distribution and its time evolution ruled by the Hamiltonian microdynamics. For every density, we have systematically obtained the dissipativeless and dissipative parts of the corresponding current density using expansions in powers of the gradients of the macrofields. With this approach, the dissipativeless part of each current density can be inferred including their nonlinear dependence on the velocity field, and their dissipative part can be identified in direct relation with entropy production. In this way, Green-Kubo formulas have been derived for all the possible transport coefficients, giving them a microscopic foundation.
The symmetries under time reversal and the point group of unbroken spatial rotations have been used to reduce the number of transport coefficients. Time-reversal symmetry leads to the Onsager-Casimir reciprocal relations between the coefficients coupling different transport processes. These latter may contribute to entropy production if their coupling is symmetric under time reversal, or be dissipativeless if their coupling is antisymmetric. Among the former, there are the heat conductivities, the viscosities, the coefficients associated with the order parameters, and the cross effects between heat transport and the order parameters. Among the latter, cross effects are furthermore predicted between heat and momentum transport, and between momentum transport and the order parameters. In isotropic phases of matter, these latter cross effects are absent according to Curie's principle based on the full continuous group of three-dimensional rotations. However, such cross effects become possible in the presence of anisotropies, such as planar interfaces between two isotropic phases [54]. Here, we have found that, in some classes of crystalline solids and liquid crystals, properties may be coupled together with a rank-three tensor of non-vanishing coefficients because of anisotropy, as it is the case for piezoelectricity [58]. The microscopic expressions of all these transport coefficients are here given by Green-Kubo formulas. For crystalline solids, previously obtained results are recovered [23,24]. Moreover, the unified approach also provides the microscopic expressions of transport properties in liquid crystals, depending on their broken continuous symmetries.
We note that the approach can be extended to quantum systems [17,18,[26][27][28]30]. In future work, we hope to use the methods developed in the present paper, in particular, to investigate the hydrodynamic and Nambu-Goldstone modes and their damping in crystalline solids and liquid crystals.