# Continuum mechanics for the elastic properties of crystals: Microscopic
approach based on projection-operator formalism

Florian Miserez  Saswati Ganguly  Rudolf Haussmann  Matthias Fuchs
Fachbereich Physik, Universität Konstanz, D-78457 Konstanz, Germany

###### Abstract

We present a microscopic derivation of the laws of continuum mechanics of
nonideal ordered solids including dissipation, defect diffusion, and heat
transport. Starting point is the classical many-body Hamiltonian. The approach
relies on the Zwanzig-Mori projection operator formalism to connect
microscopic fluctuations to thermodynamic derivatives and transport
coefficients. Conservation laws and spontaneous symmetry breaking, implemented
via Bogoliubov’s inequality, determine the selection of the slow variables.
Density fluctuations in reciprocal space encode the displacement field and the
defect concentration. Isothermal and adiabatic elastic constants are obtained
from equilibrium correlations, while transport coefficients are given as
Green-Kubo formulae, providing the basis for their measurement in atomistic
simulations or colloidal experiments. The approach and results are compared to
others from the literature.

††preprint: APS/123-QED

##  I Introduction

Continuum mechanics describe the physical properties of condensed matter at
large spatial and temporal scales. For simple liquids, the theory is called
hydrodynamics where there are five relevant variables for the five physical
degrees of freedom : the mass density $\rho_{m}$, three components of the
momentum density $j_{x}$, $j_{y}$, $j_{z}$, and the energy density $e$. In
crystalline solids, spontaneous symmetry breaking takes place and a periodic
structure arises for all physical microscopic variables. Therefore, on
macroscopic scales, three more relevant variables arise and these are the
three components of the displacement field $u_{x}$, $u_{y}$, and $u_{z}$.
These additional physical degrees of freedom means there are a total of eight
relevant variables for a crystalline solid.

However, in conventional elasticity theory there are only seven degrees of
freedom [1]. In case of defect-free ideal crystals with perfect lattice
structures, each lattice site is occupied by exactly one particle. This
constrains the density to be the divergence of the displacement field ${\bf
u}$. In this way the physical degrees of freedom are reduced by one, from
eight to seven.

On the other hand, a one-component crystalline solid at a non-zero temperature
is expected to have a finite concentration of point-defects. These point
defects may be vacant lattice sites or particles at interstitial places. As a
consequence, there is an eighth physical degree of freedom in nonideal
crystals best described by the defect density $c$. This eighth degree of
freedom, describing the diffusion of point-defects, has been introduced by
Martin, Parodi and Pershan [2]. Later, Fleming and Cohen [3] further developed
and elaborated a continuum-mechanics description based on the phenomenology of
this idea.

The first microscopic approach, to understand reversible mechanical response
and dissipative transport in crystals with local-defects, has been theorised
by Szamel and Ernst [4, 5]. They investigate the microscopic density
$n_{\mathbf{g}}(\mathbf{q})$ in Fourier representation where the wave vector
$\mathbf{k}=\mathbf{g}+\mathbf{q}$ is decomposed into the discrete reciprocal
lattice vector $\mathbf{g}$ and the continuous wave vector $\mathbf{q}$
restricted to the first Brillouin zone. Szamel and Ernst [4] suggest a
microscopic formula for the displacement field
$\mathbf{u}=u_{x}\mathbf{e}_{x}+u_{y}\mathbf{e}_{y}+u_{z}\mathbf{e}_{z}$
associated with the linear elastic response due to small deformations of the
crystal. From microscopic principles, they derive the dynamic equations for
the eight degrees of freedom governing the continuum mechanics of the system.
As a result they obtain explicit expressions for the linear elastic constants.
Later, Szamel’s [5] extension of the theory includes dissipative effects. He
applies the concepts of the projection-operator formalism [6] and derives the
Green-Kubo [7, 8, 9] relations for the transport coefficients.

The microscopic approach of Szamel and Ernst [4] has been extended by Walz and
Fuchs [10]. Their representation of the microscopic particle density
$n_{\mathbf{g}}(\mathbf{q})$ in terms of the displacement field $\mathbf{u}$
and the defect density $c$ also identifies previously ignored correlations
between displacements and defect density fields. These general theoretical
frameworks were first implemented to a specific model of cluster crystals by
Häring _et al._ [11]. Cluster crystals are defect-rich crystals, where an
inhomogeneously distributed number of soft particles occupies lattice sites
[12, 13]. An extensive examination of the elastic properties of the same
cluster crystals, taken up by Ganguly _et al._ [14], infers how local disorder
quantitatively impacts mechanical response. The perspective for the elasticity
of hard-sphere crystals, presented by Lin _et al._ [15] highlights the
influence of microscopic interactions and direct correlations on the thermo-
mechanics. Ras _et al._ [16], through their study of disordered binary
crystals, further extends the scope of these theories. An extension to include
nonlinear effects and fluctuations has been provided by Haussmann [17]. This
theory uses projection operators defined [18] for ensembles far from
equilibrium but with an assumption of local entropy maximisation. Finally, an
alternative and equivalent approach has been presented by Mabillard and
Gaspard [19, 20]. Their approach, unlike [4, 5, 10, 11, 14, 15, 16, 17],
avoids the explicit use of projection operators.

In this paper we build on the microscopic theories of Walz-Fuchs [10] and
Häring _et al._ [11]. These previous theories are restricted to the reversible
isothermal case where the temperature is constant and dissipative processes
are neglected. In this paper we consider non-constant temperatures which imply
heat transport by diffusion. Furthermore, we consider dissipative
contributions to calculate the transport coefficients of diffusive processes
and attenuation in propagative modes. We focus on general concepts and the
derivation of explicit formulas from microscopic principles.

The paper is organised as follows: Section II is devoted to developing and
explaining the microscopic basis of our theoretical framework. In it, section
II.1 briefly summarises the conceptual premise of the Mori-Zwanzig projection
operators and introduces the set of microscopic dynamical variables relevant
for this paper. Section II.2, then derives the equations of motion for these
relevant variables after defining the static and dynamic correlations
responsible for the reactive and dissipative couplings in the equations. The
main focus of section III is to obtain the coarse-grained fields of elasticity
from the microscopic fluctuating fields through an appropriate ansatz. Once
the equations of motion are obtained in the reduced space of the coarse-
grained relevant variables, section IV derives their connection to macroscopic
thermodynamic properties of the system. This is achieved by the consideration
of thermodynamic identities and expansions which allow us to interpret our
theoretical perspective in the broader context of material properties in
different thermodynamic ensembles. Finally, in section V we conclude giving
outlines of future directions.

##  II The microscopic theory

###  II.1 Densities of the relevant variables and their projected dynamics

The thermodynamics of macroscopic systems predominantly depend on the dynamics
of a few relevant variables. Specialised projection tools allow the derivation
of the dynamics of the few relevant variables from the microscopic degrees of
freedom governed by the Hamiltonian dynamics. In the crystalline phase, the
reversible parts of the equations of motion of the slowly relaxing relevant
variables govern the macroscopic mechanical response, while the dissipative
parts describe the coefficients associated with heat or momentum transport.
The Mori-Zwanzig projection operator formalism [21, 6] provides a way to
connect these measurable equilibrium and non-equilibrium thermodynamic
properties to the underlying large number of microscopic degrees of freedom.
The success of this formalism leading to an useful representation of a
material phase relies on a good choice of the relevant variables. They will be
called ’hydrodynamic variables’ in the following and will be identified in the
’hydrodynamic limit’ of small frequencies and long wavelengths. Conservation
laws and long-ranged elastic correlations arising from spontaneous symmetry
breaking provide the basis for their selection.

The equations of motion for a chosen set of relevant variables
$\\{A_{i}(t)\\}$ within a linear response framework is the eventual outcome of
the Mori-Zwanzig formalism [21, 6]. Averages
$\langle\delta\hat{A}_{k}(t)\rangle^{\text{lr}}$ of interest for time $t>0$
describe the relaxation of small initial perturbations
$\langle\delta\hat{A}_{k}(t=0)\rangle^{\text{lr}}$ of the relevant variables
generated by external fields in the past ($t<0$). The general form of these
equations and the expressions for the matrices involved are summarised here.

$\displaystyle\partial_{t}\langle\delta\hat{A}_{k}(t)\rangle^{\text{lr}}=$
$\displaystyle
i\sum_{i,j}\chi^{-1}_{ij}\omega_{jk}\langle\delta\hat{A}_{i}(t)\rangle^{\text{lr}}$
$\displaystyle-$
$\displaystyle\sum_{i,j}\int_{0}^{t}d\tau\chi^{-1}_{ij}m_{jk}(t-\tau)\langle\delta\hat{A}_{i}(\tau)\rangle^{\text{lr}}+F_{k}(t)$
(1)

The projection operation splits the time evolution of the relevant variables
into reversible (first term in Eq. (II.1)) and dissipative (second term in Eq.
(II.1)) parts. Besides these two contributions arising from the present and
earlier values of the relevant variables, there exists a random force term
$F_{k}(t)$. As we aim to derive linear elasticity theory along with the
attenuation in the elastic waves due to the dissipative couplings, only the
motion of the averaged fields are of interest. Including the fluctuating
forces would lead to stochastic equations introduced by Langevin. Here that
would correspond to fluctuating elasticity theory which is not our aim as it
would provide a far too detailed description. So we neglect the fluctuating
force $F_{k}(t)$ in our calculations. The angular brackets in Eq. (II.1)
denotes ensemble averages and the relation between the small fluctuation of
the relevant variable and the average correlation functions within the linear
response theory is highlighted by the use of index $``\text{lr}"$. Because of
the closeness to equilibrium, the dynamical parameters entering Eq. (II.1) can
be then evaluated by (grand) canonical averaging. Yet, before defining the
matrices of static susceptibility $\chi_{ij}$, frequency $\omega_{ij}$, and
memory $m_{ij}$, that appear in the dynamical equations of the averaged
fields, one needs to define the Liouville operator $\mathcal{L}$ governing the
dynamics of the microscopic variables. The Liouville operator $\mathcal{L}$,
acting on the dynamical variables of a system with a conserved phase space
volume, is defined as the Poisson bracket of an arbitrary dynamical variable
$\hat{A}$ and the Hamiltonian $\mathcal{H}$ of the system. The time evolution
of a microscopic variable $\hat{A}$ is given by [22]

$\displaystyle\partial_{t}\hat{A}(t)=\\{\hat{A}(t),\mathcal{H}\\}=i\mathcal{L}\hat{A}(t)$
(2a) $\displaystyle\hat{A}(t)=e^{i\mathcal{L}t}\hat{A}(0).$ (2b)

With the Liouville operator, the definitions of the matrices in Eq. (II.1) can
be given:

$\displaystyle\chi_{ij}=\beta\langle\delta\hat{A}^{*}_{i}\delta\hat{A}_{j}\rangle$
(3a)
$\displaystyle\omega_{jk}=\beta\langle\delta\hat{A}^{*}_{j}\mathcal{L}\delta\hat{A}_{k}\rangle$
(3b) $\displaystyle
m_{jk}(t)=\beta\langle\delta\hat{A}^{*}_{j}\mathcal{L}\mathcal{Q}e^{-i\mathcal{Q}\mathcal{L}\mathcal{Q}t}\mathcal{Q}\mathcal{L}\delta\hat{A}_{k}\rangle.$
(3c)

Here, $\beta^{-1}=k_{B}T$ $i.e.$ the temperature $T$ multiplied to the
Boltzmann constant $k_{B}$, represents the scale of energy in the system. As
$k_{B}$ is a constant, its value is set to one without any loss of generality.
The generalised Langevin equation (Eq. (II.1)) is derived [21] by splitting
the Liouville dynamics onto two orthogonal subspaces described by the
projectors

$\displaystyle\mathcal{P}=\sum_{i,j}\;\delta\hat{A}_{i}\langle\delta\hat{A}^{*}_{i}\delta\hat{A}_{j}\rangle^{-1}\langle\delta\hat{A}^{*}_{j}=\beta\sum_{i,j}\;\delta\hat{A}_{i}\rangle\chi^{-1}_{ij}\langle\delta\hat{A}^{*}_{j}$
(4a) $\displaystyle\mathcal{Q}=1-\mathcal{P}$ (4b)

The latter projector $\mathcal{Q}$ enters the matrix of memory functions
$m_{ij}(t)$.

Guided by hydrodynamic description of fluids or magnetic systems [6, 23],
conserved quantities and variables associated with spontaneously broken
continuous symmetries are chosen as the hydrodynamic variables for an ordered
solid. This includes non-ideal crystals that possess finite concentrations of
point-defects. Mass, three components of linear momentum and energy are the
conserved variables in the system. The fluctuations in the densities of the
following set of relevant variables

$\displaystyle\\{\delta(\hat{A}_{i}/V)\\}=\left(\delta\hat{\rho}_{{\bf
g}}({\bf q}),\delta\hat{j}_{\alpha=1,2,3}({\bf q}),\delta\hat{e}({\bf
q})\right)$ (5)

are chosen to describe the mechanical response and transport processes in a
three dimensional crystalline solid. Here, ${\bf q}$ is a wave vector
restricted to lie in the first Brillouin zone of reciprocal space, and the
limit ${\bf q}\to 0$ is of interest for deriving hydrodynamic theories. For a
thermodynamic system of volume $V$, the first set of variables
$\delta\hat{\rho}_{{\bf g}}({\bf q})$ are the density fluctuations with almost
the periodicity of the lattice; see their definition in Eq. (10) below. They
are indexed by the reciprocal lattice vectors ${\bf g}$, whose number we call
$N$. While $N\to\infty$ is the relevant limit, considering $N$ finite helps in
interpreting some algebraic manipulations later on. The finite values ${\bf
g}\neq 0$ reflect the broken translational symmetry of an ordered phase. The
relation of $\delta\hat{\rho}_{{\bf g}}({\bf q})$ to the microscopic particle
density (Eq. (II.1)) justifies its use as the relevant variable associated
with the conservation of mass as well. The second set
$\delta\hat{j}_{\alpha=1,2,3}$ and the third variable $\delta\hat{e}$ are the
fluctuations in the three components of linear-momentum-densities and the
energy density. The respective conservation laws are given in Eq. (13) and Eq.
(14).

All the conservation laws in this section have the general form of the
continuity equation

$\displaystyle\partial_{t}\hat{\rho}^{A}({\bf
q},t)+iq_{\alpha}\hat{j}^{A}_{\alpha}({\bf q},t)=0$ (6)

for the density $\hat{\rho}^{A}$ and current $\hat{j}^{A}$ of a conserved
variable $A$. Here, the equations are given in reciprocal space. In real
space, these equations state that for a conserved quantity, like the total
energy of the system, any small temporal change in its density at some spatial
point ${\bf r}$, will be due to a resultant current of that quantity to or
from a small volume element around that spatial point. We denote microscopic
operators with a hat on the symbol of the variable to distinguish it from the
averaged quantity of the same variable (for example see Eq. (12)). The use of
microscopic operators further implies that these conservation laws are obeyed
locally for any microscopically defined phase space function. Macroscopic
conservation laws relating the average densities to averaged currents, in the
complex dynamical system, can be obtained through ensemble or coarse-grained
averages. This is the conceptual basis for deriving the generalised Langevin
equations (Eq. (II.1)).

With particles of unit mass, the microscopic particle density operator is
$\hat{\rho}({\bf r},t)=\sum_{i=1}^{N_{p}}\delta\left({\bf r}-{\bf
r}_{i}(t)\right)$ where the sum runs over all $N_{p}$ particles, indexed by
$i$, in the system. In the Fourier space, it is

$\displaystyle\hat{\rho}({\bf k},t)$
$\displaystyle=\int\mathrm{d}^{3}re^{-i{\bf k}\cdot{\bf r}}\hat{\rho}({\bf
r},t)$ $\displaystyle=\sum_{i=1}^{N_{p}}e^{-i{\bf k}\cdot{\bf
r}_{i}(t)}=\sum_{i=1}^{N_{p}}e^{-i({\bf g}+{\bf q})\cdot{\bf r}_{i}(t)}.$ (7)

The total wave vector ${\bf k}={\bf g}+{\bf q}$ is written as a sum of the
reciprocal lattice vector ${\bf g}$ and the wave vector ${\bf q}$ of the first
Brillouin zone. This separation is possible because, in equilibrium crystals
and in the thermodynamic limit, the averages of $\hat{\rho}({\bf g}+{\bf
q},t)$ have only contributions at reciprocal lattice vectors and thus the
ensemble averaged Bragg peak amplitudes at ${\bf g}$ are

$n_{{\bf g}}=\dfrac{1}{V}\langle\hat{\rho}({\bf
g},t)\rangle=\dfrac{1}{V}\left\langle\sum_{i=1}^{N_{p}}e^{-i{\bf g}\cdot{\bf
r}_{i}(t)}\right\rangle.$ (8)

The reciprocal lattice vectors ${\bf g}$ form a Bravais lattice [24]. In the
absence of thermal fluctuations for an ideal crystal, with rigidly fixed
particles at each lattice site, $n_{{\bf g}}$ simplifies to the inverse of the
volume of the lattice unit cells. But, in any other scenario, the deviations
in $n_{{\bf g}}$ results from particle motions due to thermal fluctuations or
defect diffusion. In ordered structures like crystalline solids, the
Bogoliubov inequality [6, 25] indicates a long range correlation of the
density fluctuations $\delta\hat{\rho}_{{\bf g}}({\bf q})$, whose correlation
function diverges as $\propto q^{-2}$ for wave vectors close to all non-zero
reciprocal lattice vectors ${\bf g}\neq 0$. This is the argument for the
inclusion of the density fluctuation close to a reciprocal lattice vector,
$\delta\hat{\rho}_{{\bf g}}({\bf q})$, in the set of slow variables (Eq. (5)).
With the continuity equation for the mass density or number density for
particles of unit mass

$\displaystyle\partial_{t}\hat{\rho}({\bf g}+{\bf
q},t)+i(g+q)_{\alpha}\hat{j}_{\alpha}({\bf g}+{\bf q},t)=0,$ (9)

the amplitude at ${\bf g}=0$ still serves as the slow hydrodynamic variable
associated with the conservation of mass.

In Eq. (9) and subsequent equations the Einstein convention of summation over
repeated indices is used. Here, $\hat{j}_{\alpha}$ is a component of the
density of the linear momentum. The Fourier representation of the space and
time dependent number densities, allows us to circumvent the necessity of
defining the fields of elasticity attached to a reference lattice structure.
Instead, the microscopic operator corresponding to the fluctuations in the
density close to Bragg peaks

$\displaystyle\delta\hat{\rho}_{{\bf g}}({\bf q},t)=\hat{\rho}({\bf g}+{\bf
q},t)-n_{{\bf g}}V\delta_{{\bf q}0}$ (10)

are defined as deviation from the equilibrium ensemble averaged Bragg peak
amplitudes $n_{{\bf g}}$ (Eq. (8)) for reciprocal lattice vectors ${\bf g}$.
Now averaging $\delta\hat{\rho}_{{\bf g}}({\bf q},t)$ over the linear response
many-body distribution yields the observable microscopic field for the density
fluctuations measured within the linear hydrodynamic regime[10] when evaluated
in the small ${\bf q}$ limit

$\delta n_{{\bf g}}({\bf q},t)=\langle\delta\hat{\rho}_{{\bf g}}({\bf
q},t)\rangle^{\text{lr}}$ (11)

Here we want to reiterate the distinct notations used for microscopic
operators and averaged variables in Eq. (11). In all our calculations, we use
symbols with hat for microscopic density operators like $\hat{\rho}_{{\bf
g}}({\bf q},t)$ to distinguish them from their averaged counterparts
$\langle\hat{\rho}_{{\bf g}}({\bf q},t)\rangle^{\text{lr}}$. For notational
clarity, we choose to represent all average fluctuations by latin letters and
thus use the symbol $\delta n_{{\bf g}}({\bf q},t)$ instead of
$\langle\hat{\rho}_{{\bf g}}({\bf q},t)\rangle^{\text{lr}}$ (this connects to
the notation in Eq. (8)). These conventions are consistently followed for all
variable notations in this paper.

The next set of slow, relevant fluctuations are given by

$\delta j_{\alpha}({\bf q},t)=\langle\hat{j}_{\alpha}({\bf
q},t)\rangle^{\text{lr}}.$ (12)

justified by the conservation of linear momentum

$\displaystyle\partial_{t}\hat{j}_{\alpha}({\bf
q},t)+iq_{\beta}\hat{\sigma}_{\alpha\beta}({\bf q},t)=0,$ (13)

derived from the operator for the momentum density components
$\hat{j}_{\alpha}({\bf r},t)=\sum_{i=1}^{N_{p}}p_{\alpha}({\bf
r}_{i})\delta\left({\bf r}-{\bf r}_{i}(t)\right)$ and their corresponding
spatial Fourier transforms $\hat{j}_{\alpha}({\bf q},t)=\int d^{3}re^{-i{\bf
q}\cdot{\bf r}}\hat{j}_{\alpha}({\bf r},t)$. The term
$\hat{\sigma}_{\alpha\beta}$ is the stress tensor; see Appendix A.1 for its
definition. The conservation of energy

$\displaystyle\partial_{t}\hat{e}({\bf q},t)+i{\bf
q}\cdot{\bf\hat{j}}^{e}({\bf q},t)=0$ (14)

in an isolated system also provides a relevant variable and its coupling with
number density is related to heat transport and thermal expansion in the
system. The spatial Fourier transform for the microscopic energy density

$\displaystyle\hat{e}(r,t)$ $\displaystyle=\sum_{i=1}^{N_{p}}E({\bf
r}_{i})\delta({\bf r}-{\bf r}_{i}(t))$ (15a)
$\displaystyle=\sum_{i}\dfrac{\hat{j}_{\alpha}({\bf
r}_{i})\hat{j}_{\alpha}({\bf r}_{i})}{2m}\delta({\bf r}-{\bf
r}_{i})+\dfrac{1}{2}\sum_{i\neq j}V(r_{ij})\delta({\bf r}-{\bf r}_{i}),$ (15b)

is given by

$\displaystyle\hat{e}({\bf q},t)=\int d^{3}re^{-i{\bf q}\cdot{\bf
r}}\hat{e}({\bf r},t)$ (16)

similar to the mass and momentum densities. The microscopic definition of the
energy current $\hat{{\bf j}}^{e}$ is given in the appendix A.1. Following
arguments similar to the momentum conservation law, the relevant fluctuations
for the energy is denoted by

$\displaystyle\delta e({\bf q},t)=\langle\delta\hat{e}({\bf q},t)\rangle^{lr}$
(17)

This concludes the introduction of the slow variables enlisted in Eq. (5) and
required in providing a microscopic basis for the hydrodynamic equations of a
crystalline solid including all dissipative processes such as heat and defect
diffusion.

###  II.2 Microscopic basis for the equations of motion in the hydrodynamic
limit

In order to understand the dynamics of the relevant variables, that impact the
macroscopic properties like the elastic constants and the different transport
coefficients, we first need to focus on the microscopic definitions of the
three key quantities $\chi_{ij},\omega_{ij}$ and $m_{ij}$ (Eq. (3)) of the
generalised Langevin equation (Eq. (II.1)). In the following paragraphs
II.2.1,II.2.2,II.2.3, we define and discuss each of these quantities. Finally,
in paragraph II.2.4 we present the microscopic time evolution equations for
the set of $N+4$ relevant variables.

####  II.2.1 The static susceptibility matrix and intensive variables

The static susceptibility matrix contains the direct and the cross
correlations, measured at equilibrium, between the different relevant
variables (see Eq. (3)a). In the linear response framework, the equilibrium
susceptibility quantifies the small change in a system property, for example
density, on being subjected to an external field, for example a changing
chemical potential. Reading this relation in the opposite direction, intensive
variables $\tilde{A_{i}}$ conjugate to the selected slow variables can be
introduced. The set of these conjugate variables $\\{\tilde{A_{i}}\\}$ is
built from the general form

$\langle\delta\tilde{A}_{i}\rangle^{\text{lr}}=\sum_{j}\chi^{-1}_{ij}\langle\delta\hat{A}_{j}\rangle^{\text{lr}}.$
(18)

The generalised Langevin equation in Eq. (II.1) translates to a simpler one
written in terms of these conjugate variables; see Eq. (II.2.1) below.
Equation (18) also helps in interpreting $\chi_{ij}$. The comparison with
phenomenological approaches, e.g. the one in Ref. [3], can take place on the
level of equations combining conjugate variable sets. Moreover, the intensive
variables are used in nonlinear projection operator formalism extending our
linear response study [17]. Therefore, we use conjugate sets of variables in
the following presentation and derive explicit expressions for their equations
of motion. Starting with the relation between the three sets of relevant
variables (see Eq. (5), Eq. (11), Eq. (12), Eq. (17)), their respective
thermodynamic conjugates are defined using the static susceptibility matrix,

$\displaystyle\begin{bmatrix}\delta\mathbf{a}_{{\bf g}}({\bf q},t)\\\ \delta
b({\bf q},t)\\\ \delta\mathbf{v}({\bf q},t)\\\
\end{bmatrix}=V\mathbf{\chi^{-1}}\begin{bmatrix}\delta\mathbf{n_{{\bf
g}}}({\bf q},t)\\\ \delta e({\bf q},t)\\\ \delta\mathbf{j}({\bf
q},t)\;\end{bmatrix}.$ (19)

Here $\mathbf{n_{{\bf g}}}({\bf q},t)$ and its conjugate $\mathbf{a}_{{\bf
g}}$ denote N-dimensional column vectors whose components $n_{{\bf g}}$ (Eq.
(11)), $a_{{\bf g}}$ and their fluctuations are indexed by ${\bf g}$ (the set
of N reciprocal lattice vectors) ordered in some fixed but arbitrary way. The
$3$ dimensional vectors $\mathbf{j}$ (Eq. (12)) and $\bm{v}$ comprise of the
three Cartesian components of the linear momentum density and their respective
conjugates. Finally, $b$ is the thermodynamic conjugate to the internal energy
density $e$ (Eq. (17)) of the system. Consequently, the static susceptibility
matrix is a $(N+4)\times(N+4)$ dimensional matrix denoted by $\bm{\chi}$. The
matrix $\bm{\chi}$ contains blocks representative of self and cross
correlations between the different sets of hydrodynamic variables.

$\displaystyle\bm{\chi}({\bf q})$
$\displaystyle=\begin{bmatrix}\bm{\chi}^{\rho\rho}_{(N\times
N)}&\bm{\chi}^{\rho e}_{(N\times 1)}&\bm{\chi}^{\rho j}_{(N\times 3)}\\\
\bm{\chi}^{e\rho}_{(1\times N)}&\bm{\chi}^{ee}_{(1\times
1)}&\bm{\chi}^{ej}_{(1\times 3)}\\\ \bm{\chi}^{j\rho}_{(3\times
N)}&\bm{\chi}^{je}_{(3\times 1)}&\bm{\chi}^{jj}_{(3\times 3)}\end{bmatrix}$
(20a) $\displaystyle=\begin{bmatrix}\bm{\chi}^{\rho\rho}_{(N\times
N)}&\bm{\chi}^{\rho e}_{(N\times 1)}&0\\\ \bm{\chi}^{e\rho}_{(1\times
N)}&\bm{\chi}^{ee}_{(1\times 1)}&0\\\ 0&0&\bm{\chi}^{jj}_{(3\times
3)}\end{bmatrix}$ (20b)

The general form of $\bm{\chi}$ can be ascertained using arguments pertaining
to the time reversal symmetries of density, energy and momentum [6]. While
density and energy are even, momentum is odd with respect to time reversal.
Quantities with opposite parity under time reversal cannot have non-zero
static correlations i.e. the blocks $\bm{\chi}^{\rho j}$, $\bm{\chi}^{ej}$ and
their complex conjugates are null matrices. Constituent components of the
matrices $\bm{\chi}^{\rho\rho}_{(N\times N)}$, $\bm{\chi}^{\rho e}_{(N\times
1)}$ and $\bm{\chi}^{ee}_{(1\times 1)}$ are $\beta\langle\delta n^{*}_{{\bf
g}}({\bf q})\delta n_{{\bf g}^{\prime}}({\bf q})\rangle$, $\beta\langle\delta
n^{*}_{{\bf g}}({\bf q})\delta e({\bf q})\rangle$ and $\beta\langle\delta
e^{*}({\bf q})\delta e({\bf q})\rangle$ respectively.

Acquiring explicit expressions of the variables $\delta a_{{\bf g}}$, $\delta
b$ and $\delta v_{\alpha}$ in terms of the relevant variables (see Eq. (19)),
requires $\bm{\chi}^{-1}$ and some shorthand definitions for the constituent
matrix blocks of $\bm{\chi}$. The simple block diagonal form of the matrix
(Eq. (20)) allows the independent inversion of the diagonal blocks (see Eq.
(85)a). The diagonal block $\bm{\chi}^{jj}_{(3\times 3)}$ involving the
correlations between the momentum density fluctuations, can be inverted easily
using the classical equipartition theorem $\langle
p^{i}_{\alpha}p^{j}_{\beta}\rangle=mk_{B}T\delta_{ij}\delta_{\alpha\beta}$.
With the mass per particle $m$ set to one, the correlations between the
fluctuations of the densities of different components of linear momentum

$\displaystyle\bm{\chi}^{jj}({\bf
q})=\beta\langle\delta\hat{j}^{*}_{\alpha}({\bf q})\delta\hat{j}_{\beta}({\bf
q})\rangle=n_{0}V\delta_{{\alpha\beta}}$ (21)

leads to the identification of the field conjugate to momentum; it is the
velocity as given in Eq. (22)c (see Appendix B). The block diagonal structure
and the matrix identities in Eq. (85) are used to perform the inversion of the
$(N\times N)$ submatrix in $\bm{\chi}$ and to derive the relation between the
remaining pairs of thermodynamic conjugate variables

$\displaystyle\delta a_{\mathbf{g}}({\bf q},t)=\sum_{{\bf
g}^{\prime}}J^{*}_{{\bf g}{\bf g}^{\prime}}\delta n_{{\bf g}^{\prime}}({\bf
q},t)-U_{\bf g}({\bf q})\delta b({\bf q},t)$ (22a) $\displaystyle\delta b({\bf
q},t)=-L^{-1}({\bf q})\sum_{{\bf g}^{\prime}}U^{*}_{{\bf g}^{\prime}}({\bf
q})\delta n_{{\bf g}^{\prime}}({\bf q},t)+L^{-1}({\bf q})\delta e({\bf q},t)$
(22b) $\displaystyle\delta v_{\alpha}({\bf q},t)=n^{-1}_{0}\delta
j_{\alpha}({\bf q},t).$ (22c)

Each of the terms $J_{{\bf g}{\bf g}^{\prime}}$, $U_{{\bf g}}$ and $L$ needs
further interpretation. Given the general form (Eq. (3)a) for the components
of the matrix $\bm{\chi}$, first we focus on the block
$\bm{\chi}^{\rho\rho}_{(N\times N)}$. Let us define a $(N\times N)$
dimensional matrix ${\bf J}^{\rho\rho}$ such that ${\bf
J}^{\rho\rho}=(\bm{\chi}^{\rho\rho})^{-1}$. This inverse density correlation
matrix ${\bf J}^{\rho\rho}$ has components $J_{{\bf g}{\bf g}^{\prime}}$ and
is Hermitian. The Ornstein-Zernike relation (Eq. (23)a) provides a connection
between the density correlations and the inverse density correlation matrix.
Moreover, the components $J_{{\bf g}{\bf g}^{\prime}}$ can be obtained from
the direct correlation function $c({\bf r}_{1},{\bf r}_{2})$. Previous works
[14, 11] by the authors, explored this connection in great detail and here we
present these relations for the sake of completeness.

$\displaystyle k_{B}TV\delta_{{\bf g}{\bf g}^{\prime\prime}}=$
$\displaystyle\sum_{{\bf g}^{\prime}}\langle\delta\hat{\rho}^{*}_{{\bf
g}}({\bf q},t)\delta\hat{\rho}_{{\bf g}^{\prime}}({\bf q},t)\rangle J_{{\bf
g}^{\prime}{\bf g}^{\prime\prime}({\bf q})}$ (23a) $\displaystyle J_{{\bf
g}{\bf g}^{{}^{\prime}}}({\bf q})=$ $\displaystyle\frac{k_{B}T}{V}\int
d^{3}r_{1}\int d^{3}r_{2}e^{i{\bf g}.{\bf r}_{1}}e^{-i{\bf
g}^{{}^{\prime}}.{\bf r}_{2}}e^{i{\bf q}.({\bf r}_{1}-{\bf r}_{2})}$
$\displaystyle\left[\frac{\delta({\bf r}_{1}-{\bf r}_{2})}{n({\bf
r}_{1})}-c({\bf r}_{1},{\bf r}_{2})\right]$ (23b)

Next, we define $L$ and $U_{{\bf g}}$ through the introduction of some short
hand notations to represent components of specific correlations in the matrix
$\bm{\chi}$ (Eq. (20))

$\displaystyle K_{{\bf g}}({\bf q})=\beta\langle\delta\hat{e}^{*}({\bf
q})\delta\hat{\rho}_{{\bf g}}({\bf q})\rangle$ (24a) $\displaystyle K({\bf
q})=\beta\langle\delta\hat{e}^{*}({\bf q})\delta\hat{e}({\bf q})\rangle$ (24b)
$\displaystyle L({\bf q})=K({\bf q})-\sum_{{\bf g}{\bf g}^{\prime}}K_{{\bf
g}}({\bf q})J_{{\bf g}{\bf g}^{\prime}}({\bf q})K^{*}_{{\bf g}^{\prime}}({\bf
q})$ (24c) $\displaystyle U_{{\bf g}}({\bf q})=\sum_{{\bf
g}^{\prime}}J^{*}_{{\bf g}{\bf g}^{\prime}}({\bf q})K_{{\bf g}^{\prime}}({\bf
q})$ (24d)

The term $K_{{\bf g}}({\bf q})$ is the correlation between the fluctuations in
energy density and Bragg peak amplitude around a reciprocal lattice vector
${\bf g}$. $K({\bf q})$ is the length (${\bf q}$) dependent second moment of
the energy density. The inverse density correlation components $J_{{\bf g}{\bf
g}^{\prime}}({\bf q})$ has been defined already in Eq. (23). Multiplicative
combinations of these three terms show up as $L({\bf q})$ (Eq. (24)c) and
$U_{{\bf g}}({\bf q})$ (Eq. (24)d) in the components of the matrix
$\bm{\chi}^{-1}$. In matrix representation, these constitute the matrix
$\bm{\chi}^{-1}$ (see Appendix B) finally leading to the relations given in
Eq. (22). These microscopically defined quantities $L({\bf q})$ and $U_{{\bf
g}}({\bf q})$ will be revisited (see Eq. (40)), their small wave vector $q$
limits will be examined (see Eq. (41)) and the implications of their symmetry
properties will be discussed in section III.2.

Incorporating the thermodynamic conjugates of $\hat{A}_{i}$ (see Eq. (18) and
Eq.(22)), the generalised Langevin equations Eq. (II.1) can be re-written as

$\displaystyle\langle\partial_{t}\hat{A}_{k}(t)\rangle^{\text{lr}}$
$\displaystyle=i\sum_{j}\omega^{*}_{kj}\langle\delta\tilde{A}_{j}(t)\rangle^{\text{lr}}$
$\displaystyle-\sum_{j}\int_{0}^{t}d\tau
m^{*}_{kj}(\tau-t)\langle\delta\tilde{A}_{j}(\tau)\rangle^{\text{lr}}.$ (25)

Note the use of notations $\hat{A}_{i}$ and $\tilde{A_{i}}$ to denote
microscopic dynamical variables and their respective thermodynamic conjugates.

$\displaystyle\mathcal{P}$ $\displaystyle=\beta V^{-1}\sum_{{\bf
q}}\left(\sum_{{\bf g},{\bf g}^{\prime}}\delta n_{{\bf g}}({\bf q})\rangle
J_{{\bf g}{\bf g}^{\prime}}\langle\delta n^{*}_{{\bf g}^{\prime}}({\bf
q})+\delta e({\bf q})\rangle L^{-1}({\bf q})\langle\delta e^{*}({\bf
q})-\delta e({\bf q})\rangle L^{-1}({\bf q})\sum_{{\bf
g}^{\prime\prime}}U^{*}_{{\bf g}^{\prime\prime}}({\bf q})\langle\delta
n^{*}_{{\bf g}^{\prime\prime}}({\bf q})+n^{-1}_{0}\delta j_{\alpha}({\bf
q})\rangle\langle\delta j^{*}_{\alpha}({\bf q})\right)$ (26a)
$\displaystyle=\beta V^{-1}\sum_{{\bf q}}\left(\sum_{{\bf g}^{\prime}}\delta
a_{{\bf g}^{\prime}}({\bf q})\rangle\langle\delta n^{*}_{{\bf
g}^{\prime}}({\bf q})+\delta b({\bf q})\rangle\langle\delta e^{*}({\bf
q})+\delta v_{\alpha}({\bf q})\rangle\langle\delta j^{*}_{\alpha}({\bf
q})\right).$ (26b)

The projectors in Eq. (4), constructed from the set of relevant variables can
now be given a more specific form for our system of interest. The Eq. (26)
highlights the relation between the projection operator and the pairs of
conjugate variables in Eq. (19).

Before we present Eq. (II.2.1) for the crystalline solid with local defects,
we simplify its terms further, in section II.2.2and section II.2.3 to finally
derive the equations of motion in section II.2.4.

####  II.2.2 The frequency matrix

Once we have defined and simplified the components of the matrix $\bm{\chi}$,
we attempt to first simplify and then evaluate the components of the
$(N+4)\times(N+4)$ dimensional frequency matrix $\bm{\omega}$. This, as
defined in Eq. (3)b, governs the reversible dynamical response of a system,
pushed slightly out of equilibrium by an external perturbing field. Given the
expression Eq. (3)b for the components of the frequency matrix, the symmetry
of the relevant variables with respect to time-reversal renders the components
of $\bm{\omega}^{\rho\rho}_{(N\times N)}$, $\bm{\omega}^{\rho e}_{(N\times
1)}$, $\bm{\omega}^{ee}_{(1\times 1)}$ and $\bm{\omega}^{jj}_{(3\times 3)}$ to
zero [6]. The rest of the components have been evaluated, from their
microscopic definitions, in the Appendix A.3 (Eq. (82) and Eq. (83)).

$\displaystyle\bm{\omega}({\bf q})$
$\displaystyle=\begin{bmatrix}\bm{\omega}^{\rho\rho}_{(N\times
N)}&\bm{\omega}^{\rho e}_{(N\times 1)}&\bm{\omega}^{\rho j}_{(N\times 3)}\\\
\bm{\omega}^{e\rho}_{(1\times N)}&\bm{\omega}^{ee}_{(1\times
1)}&\bm{\omega}^{ej}_{(1\times 3)}\\\ \bm{\omega}^{j\rho}_{(3\times
N)}&\bm{\omega}^{je}_{(3\times 1)}&\bm{\omega}^{jj}_{(3\times
3)}\end{bmatrix}$ (27a) $\displaystyle=\begin{bmatrix}0&0&\bm{\omega}^{\rho
j}_{(N\times 3)}\\\ 0&0&\bm{\omega}^{ej}_{(1\times 3)}\\\
\bm{\omega}^{j\rho}_{(3\times N)}&\bm{\omega}^{je}_{(3\times
1)}&0\end{bmatrix}$ (27b)

Since our calculations, seek to provide connections between microscopically
derived equations of motion to well defined macroscopic thermodynamic
variables, we need to interpret the different correlation functions in the
small ${\bf q}$ limit. To this end, we present here the definitions of the
non-zero components $\bm{\omega}^{\rho j}_{(N\times 3)}$,
$\bm{\omega}^{j\rho}_{(3\times N)}$, $\bm{\omega}^{ej}_{(1\times 3)}$ and
$\bm{\omega}^{je}_{(3\times 1)}$ of the frequency matrix at ${\bf
q}\rightarrow 0$

$\displaystyle\bm{\omega}^{\rho j}_{{\bf g}\alpha}({\bf q})$
$\displaystyle=\beta\langle\delta\hat{\rho}^{*}_{{\bf g}}({\bf
q})\mathcal{L}\delta\hat{j}_{\alpha}({\bf q})\rangle$
$\displaystyle=-V(g+q)_{\alpha}n^{*}_{{\bf g}}+\mathcal{O}(q^{2})$ (28a)
$\displaystyle\bm{\omega}^{j\rho}_{\alpha{\bf g}}({\bf q})$
$\displaystyle=\beta\langle\delta\hat{j}^{*}_{\alpha}({\bf
q})\mathcal{L}\delta\hat{\rho}_{{\bf g}}({\bf q})\rangle$
$\displaystyle=-V(g+q)_{\alpha}n_{{\bf g}}+\mathcal{O}(q^{2})$ (28b)
$\displaystyle\bm{\omega}^{ej}_{\alpha}({\bf q})$
$\displaystyle=\beta\langle\delta\hat{e}^{*}({\bf
q})\mathcal{L}\delta\hat{j}_{\alpha}({\bf q})\rangle$
$\displaystyle=-q_{\alpha}V\left(e_{0}+p_{0}\right)+\mathcal{O}(q^{2})$ (28c)
$\displaystyle\bm{\omega}^{je}_{\alpha}({\bf q})$
$\displaystyle=\beta\langle\delta\hat{j}^{*}_{\alpha}({\bf
q})\mathcal{L}\delta\hat{e}({\bf q})\rangle$
$\displaystyle=-q_{\alpha}V\left(e_{0}+p_{0}\right)+\mathcal{O}(q^{2})$ (28d)

We evaluate them in terms of familiar thermodynamic parameters like energy and
pressure. The explicit derivations of these terms from the microscopic
expressions of the correlations are given in the Appendix A.3 (Eq. (82) and
Eq. (83)). It is important to note that the expressions for $e_{0}$ and
$p_{0}$ are achieved by taking the small wave vector ${\bf q}$ limit of the
microscopically defined, spatially varying energy and pressure, as depicted in
the Appendix A.2. The section III.1 will once again take up the discussion of
this frequency matrix $\bm{\omega}({\bf q})$ after identifying, through
coarse-graining, the connection between these definitions and the hydrodynamic
definition of the analogous quantities in terms of the elastic fields like
macroscopic density and displacements.

####  II.2.3 The memory matrix

Next we focus on the memory matrix $\bm{m}$ introduced in Eq. (II.1) and
defined in Eq. (3)c. Invoking the general continuity equation (Eq. (6)) and
the definition of the Liouville operator (Eq. (2)), it becomes evident that
the memory terms deal with the overlap between currents of the relevant
variables. Therefore, these terms are expected to have finite correlation
times. This again is a $(N+4)\times(N+4)$ dimensional Hermitian matrix
separated into blocks similar to the susceptibility and the frequency matrix.

$\displaystyle\bm{m}({\bf q},t)=\begin{bmatrix}{\bf m}^{\rho\rho}_{(N\times
N)}&{\bf m}^{\rho e}_{(N\times 1)}&{\bf m}^{\rho j}_{(N\times 3)}\\\ {\bf
m}^{*\rho e}_{(1\times N)}&{\bf m}^{ee}_{(1\times 1)}&{\bf m}^{ej}_{(1\times
3)}\\\ {\bf m}^{*\rho j}_{(3\times N)}&{\bf m}^{*ej}_{(3\times 1)}&{\bf
m}^{jj}_{(3\times 3)}\end{bmatrix}$ (29)

The components of the constituent matrix blocks are given as follows

$\displaystyle m^{\rho e}_{{\bf g}}({\bf
q},t)=\beta\langle\delta\hat{\rho}^{*}_{{\bf g}}({\bf
q})\mathcal{L}\mathcal{Q}e^{-i\mathcal{Q}\mathcal{L}\mathcal{Q}t}\mathcal{Q}\mathcal{L}\delta\hat{e}({\bf
q})\rangle.$ (30)

The explicit expression for $m^{\rho e}_{{\bf g}}({\bf q},t)$ is provided here
and the rest of the memory terms $m^{\rho\rho}_{{\bf g}{\bf g}^{\prime}}({\bf
q},t),m^{\rho j}_{{\bf g}\beta}({\bf q},t),m^{ee}({\bf
q},t),m^{ej}_{\beta}({\bf q},t),m^{jj}_{\alpha\beta}({\bf q},t)$ are analogous
in their expressions. Since the memory terms represent the dissipative
dynamics of the relevant variables and are interpreted as the overlaps between
currents of conserved quantities, they are expected to have finite relaxation
time scales. This we present as justification for our Markovian approximation
[21] to simplify our equations.

$\displaystyle\int_{0}^{t}d\tau m_{jk}({\bf
q},t-\tau)\Delta\tilde{A}_{j}(\tau)=\Gamma_{jk}({\bf
q})\Delta\tilde{A}_{j}(t)$ (31)

In the small ${\bf q}$ limit, the reduced dynamics of
$\mathcal{Q}\mathcal{L}\mathcal{Q}$ can be replaced (for proof see references
[6, 26]) with the dynamics of $\mathcal{L}$. The Markovian approximation (Eq.
(31)) in the hydrodynamic limit $({\bf q}\rightarrow 0)$ gives the matrix
comprising of Onsager transport coefficients presented in Eq. (32). In the
time evolution equation of the densities of the relevant variables, the
reactive parts of the currents results from coupling (see Eq. (3)b) of the
time derivative of one variable to another variable with opposite sign under
time reversal. The dissipative parts, on the other hand, result from coupling
of time derivatives of variables with same parity under time reversal. The
Onsager transport coefficients $\Gamma_{jk}({\bf q})$ follow from the
Markovian approximation of the memory terms and are explicitly defined as
follows

$\displaystyle\frac{\Gamma^{\rho e}_{{\bf g}}({\bf q})}{V}$
$\displaystyle=\frac{\beta}{V}\int_{0}^{\infty}\mathrm{d}t\,\langle\delta\hat{\rho}_{{\bf
g}}^{*}({\bf
q})\mathcal{L}\mathcal{Q}\mathrm{e}^{-\mathrm{i}\mathcal{L}t}\mathcal{Q}\mathcal{L}\delta\hat{e}({\bf
q})\rangle$ (32)

Here we only give the example of $\Gamma^{\rho e}_{{\bf g}}({\bf q})$ as the
other elements, $\Gamma^{\rho\rho}_{{\bf g}{\bf g}^{\prime}}({\bf
q}),\Gamma^{jj}_{{\alpha\beta}}({\bf q}),\Gamma^{ee}({\bf q})$ are built
completely in analogy. Physical significance of these transport coefficients
will be examined in section III.1 by seeking their relations to familiar
quantities like viscosity, heat conductivity etc.

####  II.2.4 Equations of motion

The equations of motion for the relevant variables, introduced in Section II.1
Eq. (5), for the crystalline solid with a finite concentration of point-
defects can now be obtained. Here, in addition to the simplifying
approximations in the small wave vector limit, we invoke the thermodynamic
conjugates, derived in Eq. (22), of the relevant variables. The relatively,
simplified equations of motion for the microscopic relevant variables are
given by

$\displaystyle\partial_{t}\delta n_{{\bf g}}({\bf q},t)=$
$\displaystyle\left(\dfrac{\omega^{\rho j}_{{\bf g}\alpha}({\bf
q})}{V}-\frac{\Gamma^{\rho j}_{{\bf g}\alpha}({\bf q})}{V}\right)\delta
v_{\alpha}({\bf q},t)$
$\displaystyle-\sum_{\mathbf{g}^{\prime}}\frac{\Gamma^{\rho\rho}_{{\bf
g}{\mathbf{g}^{\prime}}}({\bf q})}{V}\delta a_{\mathbf{g}^{\prime}}({\bf
q},t)-\frac{\Gamma^{\rho e}_{{\bf g}}({\bf q})}{V}\delta b({\bf q},t)$ (33a)
$\displaystyle\partial_{t}\delta e({\bf q},t)=$
$\displaystyle\left(\dfrac{\omega^{ej}_{\alpha}}{V}-\frac{\Gamma^{ej}_{\alpha}({\bf
q})}{V}\right)\delta v_{\alpha}({\bf q},t)$ $\displaystyle-\sum_{\bf
g}\frac{\Gamma^{e\rho}_{{\bf g}}({\bf q})}{V}\delta a_{\mathbf{g}}({\bf
q},t)-\dfrac{\Gamma^{ee}}{V}\delta b({\bf q},t)$ (33b)
$\displaystyle\partial_{t}\delta j_{\alpha}({\bf q},t)=$
$\displaystyle-\dfrac{\Gamma^{jj}_{\alpha\beta}}{V}\delta v_{\beta}({\bf
q},t)$ $\displaystyle+\left(\sum_{{\bf g}}\dfrac{\omega^{j\rho}_{\alpha{\bf
g}}}{V}-\sum_{{\bf g}}\frac{\Gamma^{j\rho}_{\alpha{\bf g}}({\bf
q})}{V}\right)\delta a_{\mathbf{g}}({\bf q},t)$
$\displaystyle+\left(\dfrac{\omega^{je}_{\alpha}}{V}-\frac{\Gamma^{je}_{\alpha}({\bf
q})}{V}\right)\delta b({\bf q},t).$ (33c)

These equations convey little intuition about the actual system without the
frame of the system’s thermodynamic properties to provide context. This
connection is achieved by coarse-graining (see section III.1) the high
dimensional space of the microscopic dynamics to the space of relevant
thermodynamic variables. Besides substituting the components of the frequency
matrix using Eq. (28), we also show how the Onsager transport coefficients
(Eq. (32)) relate to constants (Eq. (49)) obtained from the dissipative
dynamics in the system.

##  III The coarse-grained fields of elasticity

###  III.1 Coarse graining procedure

In order to derive connections between the microscopic equations of motions
and the fields associated with the elasticity of a crystalline solid, the
following ansatz was introduced in reference [10].

$\displaystyle\delta n_{{\bf g}}({\bf
q},t)=-\mathrm{i}n_{\mathbf{g}}g_{\alpha}\delta u_{\alpha}({\bf
q},t)+\frac{n_{\mathbf{g}}}{n_{0}}\delta n({\bf q},t).$ (34)

In an ideal crystal, the density fluctuation is equal to the divergence of the
displacement field defined with respect to a fixed reference lattice. So the
Bragg peak amplitude fluctuations of an ideal crystal at finite temperature
would be given by the density fluctuations or the divergence of the
displacement fields alone. This description breaks down in presence of mobile
point defects like vacancies and interstitials. The ansatz in Eq. (34) assumes
that the fluctuations in the Bragg peak amplitudes in a crystal with local
defects will have contributions from both divergence of displacements as well
as density fluctuations. This description in the Fourier space also does not
require any reference lattice for the definition of the displacement fields.

Consistent with the ansatz Eq. (34), reference [10] proposed two linear
combinations that lead to the standard thermodynamic variables, i. e. the
number density and the displacement fields from the Bragg peak amplitudes.

$\displaystyle\delta n({\bf q},t)$
$\displaystyle=\frac{n_{0}}{\mathcal{N}_{0}}\sum_{\bf
g}n_{\mathbf{g}}^{*}\delta n_{{\bf g}}({\bf q},t),$ (35a) $\displaystyle\delta
u_{\alpha}({\bf q},t)$
$\displaystyle=\mathrm{i}\mathcal{N}^{-1}_{\alpha\beta}\sum_{\bf
g}n_{\mathbf{g}}^{*}g_{\beta}\delta n_{{\bf g}}({\bf q},t),$ (35b)

The two normalizations in Eq.(35) are $\mathcal{N}_{0}=\sum_{\bf g}\lvert
n_{\mathbf{g}}\rvert^{2}$ and $\mathcal{N}_{\alpha\beta}=\sum_{\bf g}\lvert
n_{\mathbf{g}}\rvert^{2}g_{\alpha}g_{\beta}$. We get accordingly $\langle
n({\bf q},t)\rangle=N\delta_{{\bf q},0}$ and $\langle u_{\alpha}({\bf
q},t)\rangle=0$. Also the symmetry argument, $\sum_{g}\lvert
n_{\mathbf{g}}\rvert^{2}g_{\beta}=0$, essentially ensures orthogonality of
$\delta n({\bf q},t)$ and $\delta u_{\alpha}({\bf q},t)$. Here we want to
point out how our definition of the density field (Eq. (35a)) is different
from alternative approaches addressing the question of mechanical response in
crystals with defects. The coarse-grained density fields in the theories
developed by Haussmann [17] and Mabillard _et al._ [20] arise from density
fluctuations only around the centre of reciprocal space (viz. at ${\bf g}=0$)
while we define the density field(Eq. 35a) as a sum over fluctuations around
all reciprocal lattice vectors in the system. Implications of this difference
will become discernible only through future implementation of these approaches
to study specific crystalline solids.

The large set of variables $\delta n_{{\bf g}}({\bf q})$ (Eq. (11)), having
been reduced to four coarse-grained variables $\delta n({\bf q})$ and $\delta
u_{\alpha=1,2,3}({\bf q})$ through the ansatz Eq. (34) leaves us with eight
coarse-grained hydrodynamic variables. These variables, paired with their
respective thermodynamic conjugates, give us a coarse-grained version of the
Eq. (19)

$\displaystyle\begin{bmatrix}\delta a({\bf q},t)\\\ \delta\mathbf{y}({\bf
q},t)\\\ \delta b({\bf q},t)\\\ \delta\mathbf{v}({\bf q},t)\\\
\end{bmatrix}=V\mathbf{\chi^{-1}}\begin{bmatrix}\delta n({\bf q},t)\\\
\delta\mathbf{u}({\bf q},t)\\\ \delta e({\bf q},t)\\\ \delta\mathbf{j}({\bf
q},t)\;\end{bmatrix}.$ (36)

Each pair of conjugate variables, along with reference to linear response
relations, are derived and discussed in section III.2 as well as appendix D.
Attention is given to the vector $\delta\mathbf{u}$ constituted of $\delta
u_{\alpha=1,2,3}$ in three dimensions since its relation to its conjugate
$\delta\bm{y}$ with constituents $\delta y_{\alpha=1,2,3}$ leads to the
coefficients of stiffness within linear elasticity.

###  III.2 The coarse-grained projectors and coarse-grained intensive
variables

Now that we have established the ground work by deriving the important
connections between the microscopic and the coarse-grained framework, we can
go on to derive the equations of motion for the coarse-grained relevant
variables (see Eq. (52)). The initial step, towards deriving and interpreting
these equations of interest, is to obtain the coarse-grained projection
operators which shall be denoted by $\tilde{\mathcal{P}}$ and
$\tilde{\mathcal{Q}}$. For that, we need to revisit the thermodynamic
conjugate variables in Eq. (19) and Eq. (22). Then define $\delta a$, $\delta
b$ using the coarse-graining ansatz. Let us first establish the relations
between the coarse-grained thermodynamic conjugates $\delta a$ and $\delta
y_{\alpha}$ and their microscopic counterpart $\delta a_{{\bf g}}$. Drawing
analogy to Eq. (35), these quantities are

$\displaystyle\delta a({\bf q},t)=\dfrac{1}{n_{0}}\sum_{{\bf g}}n^{*}_{{\bf
g}}\delta a_{{\bf g}}({\bf q},t)$ (37a) $\displaystyle\delta y_{\alpha}({\bf
q},t)=i\sum_{{\bf g}}g_{\alpha}n^{*}_{{\bf g}}\delta a_{{\bf g}}({\bf q},t).$
(37b)

Next we substitute the $\delta n_{{\bf g}}$ with the fluctuating fields
$\delta n$ and $\delta u_{\alpha}$ using the ansatz in Eq. (34) in the
definitions of $\delta a_{{\bf g}}$ and $\delta b$ given in Eq. (22) and
repeated here with the substitutions,

$\displaystyle\delta a_{\mathbf{g}}({\bf q},t)$
$\displaystyle=\dfrac{1}{n_{0}}\sum_{{\bf g}^{\prime}}J^{*}_{{\bf g}{\bf
g}^{\prime}}n_{{\bf g}^{\prime}}\delta n({\bf q},t)-i\sum_{{\bf
g}^{\prime}}J^{*}_{{\bf g}{\bf g}^{\prime}}n_{{\bf
g}^{\prime}}g^{\prime}_{\alpha}\delta u_{\alpha}({\bf q},t)$ $\displaystyle-
U_{\bf g}({\bf q})\delta b({\bf q},t)$ (38a) $\displaystyle\delta b({\bf
q},t)$ $\displaystyle=-\dfrac{L^{-1}}{n_{0}}\sum_{{\bf g}^{\prime}}U^{*}_{{\bf
g}^{\prime}}({\bf q})n_{{\bf g}^{\prime}}\delta n({\bf q},t)$
$\displaystyle+iL^{-1}\sum_{{\bf g}^{\prime}}U^{*}_{{\bf g}^{\prime}}({\bf
q})n_{{\bf g}^{\prime}}g^{\prime}_{\alpha}\delta u_{\alpha}({\bf
q},t)+L^{-1}({\bf q})\delta e({\bf q},t).$ (38b)

The relations in Eq. (22) between the pairs of conjugate variables shown in
Eq. (19) requires the inversion of the static correlation matrix. Appendix B,
Eq. (87), explains how this involves the terms defined in Eq. (24). Here, for
deriving the relations given in Eq. (39), similar mathematical manipulations
[26] are used in conjunction with the coarse-graining ansatz Eq. (34).
Plugging in the expression for $\delta a_{{\bf g}}$ from Eq. (38)a to the
equations in Eq. (37) leads to expressions for $\delta a$ (Eq. (39)a) and
$\delta y_{\alpha}$ (Eq. (39)c) in terms of the generalised elastic
coefficients $\nu,\mu_{\alpha},\lambda_{{\alpha\beta}}$ given in Appendix C.1.
This allows us to write $\delta a,\delta b,\delta y_{\alpha}$ (Eq. (39)) in
terms of the fields of elasticity $\delta n,\delta u_{\alpha}$ and quantities
like $\nu,\mu_{\alpha},\lambda_{{\alpha\beta}},\theta,\tau_{\alpha}$ which
characterise the generalised material response of the solid. The relation
between $\delta j_{\alpha}$ and its conjugate $\delta v_{\alpha}$ is repeated
in this list to complete the set of relevant variables.

$\displaystyle\delta a({\bf q},t)=$ $\displaystyle\dfrac{\nu}{n^{2}_{0}}\delta
n({\bf q},t)-\dfrac{\mu_{\beta}({\bf q})}{n_{0}}\delta u_{\beta}({\bf q},t)$
$\displaystyle-$ $\displaystyle\frac{\theta^{*}({\bf q})}{n_{0}}\delta b({\bf
q},t),$ (39a) $\displaystyle\delta b({\bf q},t)=$ $\displaystyle L^{-1}({\bf
q})\left(-\frac{\theta({\bf q})}{n_{0}}\delta n({\bf q},t)+\tau_{\alpha}({\bf
q})\delta u_{\alpha}({\bf q},t)\right)$ $\displaystyle+$ $\displaystyle
L^{-1}({\bf q})\delta e({\bf q},t),$ (39b) $\displaystyle\delta
y_{\alpha}({\bf q},t)=$ $\displaystyle-\dfrac{\mu^{*}_{\alpha}({\bf
q})}{n_{0}}\delta n({\bf q},t)+\lambda_{\alpha\beta}({\bf q})\delta
u_{\beta}({\bf q},t)$ $\displaystyle+$ $\displaystyle\tau^{*}_{\alpha}({\bf
q})\delta b({\bf q},t)$ (39c) $\displaystyle\delta v_{\alpha}({\bf q},t)=$
$\displaystyle n^{-1}_{0}\delta j_{\alpha}({\bf q},t)$ (39d)

As already mentioned in section III.1, the variables $\delta y_{\alpha}({\bf
q})$ are components of the three dimensional vector ${\bf\delta y}$ which is
conjugate to the three dimensional vector ${\bf\delta u}$ constituted of the
displacement fields $\delta u_{\alpha}({\bf q})$. The relation (Eq. (35)) of
$\delta n$ and $\delta u_{\alpha}$ to the microscopic fields are reflected in
their respective coarse-grained conjugates $\delta a$ and $\delta y_{\alpha}$.
It would be pertinent here to note that the closed form representation (Eq.
(39)) of the coarse-grained conjugate fields $\delta a({\bf q},t)$, $\delta
b({\bf q},t)$ and $\delta y_{\alpha}({\bf q},t)$ require the use of a set of
generalised elastic coefficients $\lambda_{{\alpha\beta}}({\bf
q}),\mu_{\alpha}({\bf q})$, and $\nu({\bf q})$, which were derived and
validated in previous papers. Appendix C summarises their microscopic
definitions for the sake of completeness. References [14, 11, 10] present
these important mechanical properties and by detailing their relation to the
inverse density correlation matrix (Eq. (23)), show ways of determining them
in crystalline solids with local defects. Unlike this contribution, the
previous works formulate the theory keeping only the reversible contributions
(only the first term on the right hand side of Eq. (II.1)) to the time
evolution equations of the relevant variables and ignores all contributions of
energy fluctuations or transport. They enter via two new coefficients

$\displaystyle\tau_{\alpha}({\bf q})$ $\displaystyle=\mathrm{i}\sum_{\bf
g}U^{*}_{\bf g}({\bf q})n_{\mathbf{g}}g_{\alpha}$ (40a)
$\displaystyle\theta({\bf q})$ $\displaystyle=\sum_{\bf g}U^{*}_{\bf g}({\bf
q})n_{\mathbf{g}}$ (40b)

defined using quantities given in Eq. (24). Appendix C.2 discusses
$\tau_{\alpha}({\bf q})$, $\theta({\bf q})$ to explain their microscopic
origin and symmetries in the thermodynamic limit. Given the definition of
$\tau_{\alpha}({\bf q})$ in Eq. (40), Eq. (38)b transforms to Eq. (39)b.

As material properties at long wavelengths are of interest to us, here we
present the ${{\bf q}}\rightarrow 0$ limit of the variables defined in Eq.
(24) and Eq. (40). The mathematical manipulations, starting with the initial
definitions and leading to the final expressions, are detailed in Appendix
C.2. Here, we recapitulate the final expressions for all the quantities
$\tau_{\alpha}({\bf q}),\theta({\bf q}),L({\bf
q}),\lambda_{{\alpha\beta}}({\bf q}),\mu_{{\alpha\beta}}({\bf q}),\nu({\bf
q})$ at ${\bf q}\rightarrow 0$

$\displaystyle\lim_{{\bf q}\rightarrow 0}\theta({\bf q})$
$\displaystyle=\theta,$ (41a) $\displaystyle\lim_{{\bf q}\rightarrow
0}\tau_{\alpha}({\bf q})$
$\displaystyle=\mathrm{i}q_{\beta}\tau_{\alpha\beta},$ (41b)
$\displaystyle\lim_{{\bf q}\rightarrow 0}L({\bf q})$ $\displaystyle=L.$ (41c)

In the limit of macroscopic lengths, the leading order terms in $\theta({\bf
q})$ and $L({\bf q})$ become ${\bf q}$ independent real constants.
$\tau_{\alpha}({\bf q})$ at ${\bf q}\rightarrow 0$, on the other hand, has
symmetries similar to the generalised elastic constant $\mu_{\alpha}({\bf q})$
which is associated [11] with the coupling between the coarse-grained density
$\delta n$ and displacement fields $\delta u_{\alpha}$ in a defect rich
crystal; the symmetery holds $\tau_{\alpha\beta}=\tau_{\beta\alpha}$, as can
be shown with the rotational LMBW equation [11, 27, 28]. As we show in the
Appendix C.1, the elastic coefficients $\lambda_{{\alpha\beta}}$,
$\mu_{\alpha}$ and $\nu$ are the same as the ones in [14, 11, 10] and in the
hydrodynamic limit, they become

$\displaystyle\lim_{{\bf q}\to 0}\lambda_{\alpha\beta}({\bf
q})=\lambda_{\alpha\beta\gamma\delta}q_{\gamma}q_{\delta}$ (42a)
$\displaystyle\lim_{{\bf q}\to 0}\mu_{\alpha}({\bf
q})=i\mu_{\alpha\beta}q_{\beta}$ (42b) $\displaystyle\lim_{{\bf q}\to
0}\mu^{*}_{\alpha}({\bf q})=-i\mu_{\alpha\beta}q_{\beta}$ (42c)
$\displaystyle\lim_{{\bf q}\to 0}\nu({\bf q})=\nu$ (42d)

Therefore, the conjugate variables in Eq. (39), in the ${\bf q}\rightarrow 0$,
are as follows

$\displaystyle\delta a(t)=\dfrac{\nu}{n^{2}_{0}}\delta
n(t)-\dfrac{\mu_{\alpha\beta}}{n_{0}}\delta
u_{\alpha\beta}(t)-\frac{\theta}{n_{0}}\delta b(t)$ (43a) $\displaystyle\delta
b(t)=L^{-1}\left(-\frac{\theta}{n_{0}}\delta n(t)+\tau_{\alpha\beta}\delta
u_{{\alpha\beta}}(t)+\delta e(t)\right)$ (43b) $\displaystyle\delta
y_{{\alpha\beta}}(t)=-\dfrac{\mu_{\alpha\beta}}{n_{0}}\delta
n(t)+\tau_{\alpha\beta}\delta b(t)$ $\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \
+\left(\lambda_{\alpha\gamma\beta\delta}+\lambda_{\alpha\delta\gamma\beta}-\lambda_{\alpha\beta\gamma\delta}\right)\delta
u_{\gamma\delta}(t)$ (43c)

Here,

$\displaystyle\delta y_{\alpha}=-i\delta y_{{\alpha\beta}}q_{\beta}$ (44)

has been used as the definition of $\delta y_{{\alpha\beta}}$ to obtain Eq.
(43)c from $\delta y_{\alpha}$ of Eq. (39)c. Given the symmetries [11, 14]
$\mu_{{\alpha\beta}}=\mu_{\beta\alpha}$,
$\tau_{{\alpha\beta}}=\tau_{\beta\alpha}$ and
$\lambda_{{\alpha\beta}\gamma\delta}=\lambda_{\beta\alpha\gamma\delta}=\lambda_{{\alpha\beta}\delta\gamma}=\lambda_{\gamma\delta{\alpha\beta}}$,
the second rank tensor $\delta y_{\alpha\beta}=\delta y_{\beta\alpha}$ is also
symmetric. The definition of $\delta y_{\alpha\beta}$ in Eq. (44) is also a
reference to the fact that in our calculations, we have chosen to restrict the
displacement fluctuations to linear orders in strain by choosing
$iq_{\beta}\delta u_{\alpha}({\bf q}\rightarrow 0)=\delta u_{{\alpha\beta}}$.
Moreover, the orgin of the third term in Eq. (43c) can be tracked [26, 11] to
the definition of the symmetric strain field (Eq. (45)). Because of the
symmetry of $\mu_{{\alpha\beta}},\tau_{\alpha\beta}$, as well as of
$\lambda_{\alpha\beta\gamma\delta}$, only the symmetric part of the strain
tensor enters. We therefore can redefine :

$\displaystyle\delta u_{{\alpha\beta}}\to\frac{1}{2}(\delta
u_{{\alpha\beta}}+\delta u_{{\beta\alpha}})=\dfrac{1}{2}i(q_{\beta}\delta
u_{\alpha}+q_{\alpha}\delta u_{\beta}).$ (45)

As conjugates to strains $\delta u_{{\alpha\beta}}$, the quantity $\delta
y_{{\alpha\beta}}$ can be interpreted as linear stress variables. The relation
between $\delta u_{{\alpha\beta}}$ and $\delta y_{{\alpha\beta}}$ as a pair of
thermodynamic conjugates, is discussed in the context of thermodynamic free
energy, in section IV.2 (see Eq. (36) and Eq. (56)).

Having derived the coarse-grained hydrodynamic variables and their
thermodynamic conjugates, now we can define the coarse-grained projection
operators that will allow us to obtain the macroscopic equilibrium and
transport properties of the system. The conservation law (see Eq. (9)) implies

$\displaystyle\mathcal{Q}\mathcal{L}\;\delta\hat{\rho}_{{\bf g}}({\bf
q})=-(g+q)_{\alpha}\mathcal{Q}\;\delta\hat{j}_{\alpha}({\bf g}+{\bf q})$ (46)

where the microscopically defined projector $\mathcal{Q}$ is acting on the
mass conservation equation. Let us then take the time derivative of the
coarse-graining ansatz in Eq. (34) and look at the dynamics projected by
$\mathcal{Q}$

$\displaystyle i\mathcal{Q}\mathcal{L}\delta n_{{\bf g}}({\bf
q})=-i^{2}n_{{\bf g}}g_{\alpha}\mathcal{Q}\mathcal{L}\delta u_{\alpha}({\bf
q})+\dfrac{n_{{\bf g}}}{n_{0}}i\mathcal{Q}\mathcal{L}\delta n({\bf q})$ (47a)
$\displaystyle\implies$ $\displaystyle\mathcal{Q}\mathcal{L}\delta
n_{\mathbf{g}}({\bf
q})=-n_{\mathbf{g}}g_{\alpha}\tilde{\mathcal{Q}}\delta\dot{u}_{\alpha}({\bf
q})$ (47b)

In Eq. (47)a, the second term on the right, related to mass current
$\mathcal{L}\delta n$ vanishes if one postulates the conservation of number
density $\partial_{t}\delta n=iq_{\alpha}\delta j_{\alpha}({\bf q},t)$ and
realises that the linear momentum densities belonging to the eigen-space of
the projector do not contribute to the dynamics projected by
$\tilde{\mathcal{Q}}$. Arriving at Eq. (47)b through this
$\tilde{\mathcal{Q}}\mathcal{L}\delta n({\bf
q})=-q_{\alpha}\tilde{\mathcal{Q}}\delta j_{\alpha}({\bf q})=0$, uses the fact
that the components of linear momentum $j_{\alpha}({\bf q})$ in the 1st BZ are
conserved variables. However, that does not automatically imply conservation
of $j_{\alpha}({\bf g}+{\bf q})$ for ${\bf g}\neq 0$.

In Eq. (47) we arrive at the relation between the microscopic projector
$\mathcal{Q}$ and $\tilde{\mathcal{Q}}$ obtained after the coarse-graining
ansatz Eq. (34). The coarse-grained space of a considerably smaller number of
slow relevant variables, contains the number density $\delta n$, three
components of linear momentum $\delta j_{\alpha}$, three components of
displacement fields $\delta u_{\alpha}$ and energy density $\delta e$. These
eight variables are compatible with the eight hydrodynamic modes expected [3]
in an equilibrium crystalline solid at finite temperatures. This new smaller
set of coarse-grained relevant variables require appropriate projectors
written using them.

$\displaystyle\tilde{\mathcal{P}}=1-\tilde{\mathcal{Q}}$ $\displaystyle=\beta
V^{-1}\sum_{{\bf q}}\left[\delta a({\bf q})\rangle\langle\delta n^{*}({\bf
q})+\delta y_{\alpha}({\bf q})\rangle\langle\delta u^{*}_{\alpha}({\bf
q})+\delta b({\bf q})\rangle\langle\delta e^{*}({\bf q})+\delta
v_{\alpha}({\bf q})\rangle\langle\delta j^{*}_{\alpha}({\bf q})\right]$ (48)

A comparison of this coarse-grained projector to the microscopic projector in
Eq. (26) shows how all the terms involving fluctuations of Bragg peak
amplitudes at the reciprocal lattice vectors, get included in the first two
terms in the coarse-grained projector (Eq. (48)) through the use of the ansatz
Eq. (34). Before we can go on to obtain the dynamical equations of the coarse-
grained variables, we need to determine the implication of the ansatz Eq. (34)
for the dissipative terms in the equations of motion (Eq. (33)).

###  III.3 The dissipative terms or the memory terms derived for the coarse-
grained fields

Equipped with the coarse-grained projectors $\tilde{\mathcal{P}}$ and
$\tilde{\mathcal{Q}}$, we try to make sense of the transport coefficients in
Eq. (32). We rewrite them with the coarse-grained projector and using the
relation in Eq. (47)b.

$\displaystyle\frac{\Gamma^{*\rho\rho}_{{\bf g}{\bf g}^{\prime}}({\bf q})}{V}$
$\displaystyle=\frac{\beta}{V}n_{\mathbf{g}}n_{\mathbf{g}^{\prime}}^{*}g_{\alpha}g_{\beta}^{\prime}\int_{0}^{\infty}\mathrm{d}t\,\langle\delta\dot{u}_{\alpha}^{*}({\bf
q})\tilde{\mathcal{Q}}\mathrm{e}^{-\mathrm{i}\mathcal{L}t}\tilde{\mathcal{Q}}\delta\dot{u}_{\beta}({\bf
q})\rangle^{*}$
$\displaystyle=n_{\mathbf{g}}n_{\mathbf{g}^{\prime}}^{*}g_{\alpha}g_{\beta}^{\prime}\zeta_{\alpha\beta},$
(49a) $\displaystyle\frac{\Gamma^{*\rho e}_{{\bf g}}({\bf q})}{V}$
$\displaystyle=\frac{\beta}{V}n_{\mathbf{g}}g_{\alpha}q_{\beta}\int_{0}^{\infty}\mathrm{d}t\,\langle\delta\dot{u}_{\alpha}^{*}({\bf
q})\tilde{\mathcal{Q}}\mathrm{e}^{-\mathrm{i}\mathcal{L}t}\tilde{\mathcal{Q}}j_{\beta}^{e}({\bf
q})\rangle^{*}$
$\displaystyle=n_{\mathbf{g}}g_{\alpha}q_{\beta}\xi^{\top}_{\alpha\beta},$
(49b) $\displaystyle\frac{\Gamma^{*e\rho}_{{\bf g}}({\bf q})}{V}$
$\displaystyle=\frac{\beta}{V}n_{\mathbf{g}}^{*}g_{\beta}q_{\alpha}\int_{0}^{\infty}\mathrm{d}t\,\langle
j_{\alpha}^{e*}({\bf
q})\tilde{\mathcal{Q}}\mathrm{e}^{-\mathrm{i}\mathcal{L}t}\tilde{\mathcal{Q}}\delta\dot{u}_{\beta}({\bf
q})\rangle^{*}$
$\displaystyle=n_{\mathbf{g}}^{*}g_{\beta}q_{\alpha}\xi_{\alpha\beta}.$ (49c)

Note how the $\mathcal{Q}\mathcal{L}\delta n_{{\bf g}}$ replaced by $-n_{{\bf
g}}g_{\alpha}\tilde{\mathcal{Q}}\delta\dot{u}_{\alpha}({\bf q})$ in the Eq.
(III.3) is related to the fluctuations of linear momentum at ${\bf g}\neq 0$
through the conservation law in Eq. (9) and hence Eq. (46). In the limit of
long time and small wave vectors, the integrals in Eq. (49) represent the
transport coefficients $\zeta_{{\alpha\beta}}$, $\xi_{{\alpha\beta}}$ and
$\xi^{T}_{{\alpha\beta}}$. Their symmetries are consistent [26] with the
Onsager reciprocal relations [29, 30] dictated by the symmetry of the
Hamiltonian under time-reversal. The physical significance of these transport
coefficients will become evident when discussed, in section III.4 and section
IV.3, in the context of the dynamical equations of the hydrodynamic variables.

Next consider the Onsager transport coefficient
$\Gamma^{jj}_{{\alpha\beta}}({\bf q})$ and then, in it, substitute the term
$\mathcal{L}\delta\hat{j}_{\alpha}({\bf q})$ using the conservation of the
linear momentum (Eq. (13)). This leads to the viscosity tensor
$\eta_{{\alpha\beta}{\gamma\delta}}$ given by the integral in the Eq. (III.3).

$\displaystyle\frac{\Gamma^{jj}_{{\alpha\beta}}({\bf q})}{V}$
$\displaystyle=q_{\beta}q_{\gamma}\eta_{{\alpha\beta}{\gamma\delta}}$
$\displaystyle=\frac{\beta}{V}q_{\beta}q_{\gamma}\int_{0}^{\infty}\mathrm{d}t\,\langle\sigma_{\alpha\beta}^{*}({\bf
q})\tilde{\mathcal{Q}}\mathrm{e}^{-\mathrm{i}\mathcal{L}t}\tilde{\mathcal{Q}}\sigma_{\gamma\delta}({\bf
q})\rangle^{*}.$ (50)

Another transport coefficient, arises from the memory term associated with the
energy conservation (Eq. (14))

$\displaystyle\frac{\Gamma^{ee}({\bf q})}{V}$
$\displaystyle=q_{\alpha}q_{\beta}\alpha_{\alpha\beta}T$
$\displaystyle=\frac{\beta}{V}q_{\alpha}q_{\beta}\int_{0}^{\infty}\mathrm{d}t\,\langle
j_{\alpha}^{e*}({\bf
q})\tilde{\mathcal{Q}}\mathrm{e}^{-\mathrm{i}\mathcal{L}t}\tilde{\mathcal{Q}}j_{\beta}^{e}({\bf
q})\rangle^{*}.$ (51)

Note that both projectors $\mathcal{Q}$ and $\tilde{\mathcal{Q}}$ could be
used in Eqs. (III.3) and (III.3) as they act identically in both kernels. The
reason is time-parity.

The components of $\bm{\Gamma}^{\rho j}_{(N\times 3)}$,
$\bm{\Gamma}^{ej}_{(1\times 3)}$ and their conjugate transposes in the
dynamical equations Eq. (33) have been neglected in our calculations. It can
be shown [26] that the leading ${\bf q}$ dependent term for these components
arises from $\mathcal{O}\left({\bf q}({\bf g}+{\bf q})\right)$ while all the
other components of $\bm{\Gamma}$ has $\mathcal{O}({\bf g}+{\bf q})$ (see Eq.
(49), Eq. (III.3) and Eq. (III.3)) leading order terms.

After writing the complete equations of motions for the coarse-grained
relevant variables in section III.4, physical significance of these transport
coefficients will be discussed in section III.4 as well as in section IV.3.

###  III.4 Equations of motion of the coarse-grained fields of elasticity

Incorporating all the coarse-grained variables in the microscopic equations of
motion in Eq. (33) one obtains the full time evolution equations for the
fields of elasticity. Appendix E explains the steps involved in the derivation
of these final equations (Eq. (52)). This derivation of the hydrodynamic
equations Eq. (52) in conjunction with the definitions of the material
constants in section III.2 and the Onsager transport coefficients in section
III.3 provides a microscopic basis for all the reactive and dissipative
coefficients that govern the macroscopic static and dynamic properties of
crystalline solids with local-defects.

$\displaystyle\partial_{t}\delta n({\bf q},t)$
$\displaystyle=-\mathrm{i}n_{0}q_{\alpha}\delta v_{\alpha}({\bf q},t)$ (52a)
$\displaystyle\partial_{t}\delta e({\bf q},t)$
$\displaystyle=-\mathrm{i}\left(e_{0}+p_{0}\right)q_{\alpha}\delta
v_{\alpha}({\bf q},t)$ $\displaystyle-q_{\alpha}\xi_{\alpha\beta}\delta
y_{\beta\gamma}({\bf
q},t)q_{\gamma}-q_{\alpha}q_{\beta}\alpha_{\alpha\beta}T\delta b({\bf q},t)$
(52b) $\displaystyle\partial_{t}\delta j_{\alpha}({\bf q},t)$
$\displaystyle=-in_{0}q_{\alpha}\delta a({\bf q},t)+i\delta
y_{\alpha\beta}({\bf q},t)q_{\beta}$
$\displaystyle-\mathrm{i}\left(e_{0}+p_{0}\right)q_{\alpha}\delta b({\bf
q},t)-q_{\beta}q_{\gamma}\eta_{{\alpha\beta}{\gamma\delta}}\delta
v_{\delta}({\bf q},t)$ (52c) $\displaystyle\partial_{t}\delta{u_{\alpha}({\bf
q},t)}$ $\displaystyle=\delta v_{\alpha}({\bf q},t)$
$\displaystyle+i\zeta_{\alpha\beta}\delta y_{\beta\gamma}({\bf
q},t)q_{\gamma}-\mathrm{i}q_{\beta}\xi_{\alpha\beta}^{\top}\delta b({\bf
q},t).$ (52d)

At this point, we have the time evolution equation of the eight slow variables
in the system. Comparing the Eq. (52) to the analogous equations in case of a
fluid [6, 31], one immediately concludes that the terms with $\delta
y_{{\alpha\beta}}$ would be absent in the fluid. In case of a fluid, there are
only five hydrodynamic modes and a displacement field is ill-defined. So only
five equations of motion specifying the time evolution of number density,
energy density and linear-momentum density, would be relevant. This analogy to
fluid equations will allow us to identify the viscosity $\eta_{\alpha\beta}$
and transport coefficient for heat conductivity $\alpha_{\alpha\beta}$. We are
left with two more dissipative constants $\zeta_{{\alpha\beta}}$ and
$\xi_{{\alpha\beta}}$ each of which arises in an ordered solid where the
displacement fields need to be treated as separate variables justified by the
spontaneous breaking of translational invariance. This enters the possibility
of defect motion into Eq. (52), because the divergence of the displacement
field is not tied to the density change as would hold in an ideal solid
without defects.

##  IV Connection to thermodynamics

The original Mori-Zwanzig projection operator formalism, that eventually leads
to these equations, does not provide a recipe for choosing an optimal set of
slow variables that best captures the macroscopic properties of a given
system. Conventionally, the variables associated with conservation laws or
spontaneously broken continuous symmetries are found to be best suited for
this purpose. Validating our chosen set of relevant variables, requires
examination of how our equations of motion, in the small wave vector limit,
relates to the thermodynamic parameters determining the equilibrium ensemble
of the system. In section IV, we present a detailed derivation of the
connections between the dynamics of the relevant variables and thermodynamics
of the crystal with finite concentration of local defects.

###  IV.1 Free energy expansion and thermodynamic identities

Any thermodynamic system at equilibrium, is characterised by a minimum in the
free energy and a maximum in the entropy. Therefore, small changes in the
thermodynamic parameters characterising the system, result in small increments
in the free energy. Relying on this conceptual premise, the free energy
density of a point-defect rich crystalline solid [23], can be expanded around
the free energy density of a reference un-deformed
$(\nabla_{\beta}u_{\alpha}=0)$ crystal phase at density and temperature
$n_{0}$ and $T_{0}$ respectively.

$\displaystyle f(n_{0}+\delta n,0+\delta u_{\alpha\beta},T_{0}+\delta T)$
$\displaystyle=$ $\displaystyle f^{0}+\mu^{0}\delta
n+\frac{1}{2}A\left(\frac{\delta n}{n_{0}}\right)^{2}+B_{\alpha\beta}\delta
u_{{\alpha\beta}}\frac{\delta n}{n_{0}}$ $\displaystyle+$
$\displaystyle\frac{1}{2}C^{n}_{{\alpha\beta}{\gamma\delta}}\delta
u_{{\alpha\beta}}\delta u_{\gamma\delta}-s^{0}\delta T+D_{s}\frac{\delta
n}{n}\delta T$ $\displaystyle+$ $\displaystyle\frac{1}{2}E\delta
T^{2}+F_{\alpha\beta}\delta T\cdot\delta u_{{\alpha\beta}},$ (53)

The linear strain field $\delta
u_{\alpha\beta}=\dfrac{1}{2}(\nabla_{\beta}u_{\alpha}+\nabla_{\alpha}u_{\beta})$
is symmetrised in its definition because an asymmetry in the strain field in
the thermodynamic limit amounts to a free energy conserving rotation of the
entire system. The strain multiplied with the system volume $V$ is an
extensive thermodynamic variable and we define $h_{{\alpha\beta}}$ as its
conjugate intensive field. The coefficients of the quadratic couplings between
the thermodynamic parameters are related to the curvature of the free energy
along specific directions of the thermodynamic parameter space. Thus, the
leading order terms in this expansion, provides important information related
to the mechanical and thermal constants in the system.

The thermodynamic conjugates of the three independent variables $n$,
$u_{\alpha\beta}$ and $T$, are the chemical potential $\mu$, a tensor
$h_{\alpha\beta}$ and entropy density $s$. Here, we are working with the free
energy density $f=F/V$ and hence we use densities of the extensive variables
like number of particles or entropy per volume. From the free energy expansion
Eq. (IV.1) and knowledge [23] about the pair of thermodynamic conjugates, we
get the following relations,

$\displaystyle\delta\mu$ $\displaystyle=\frac{A}{n_{0}^{2}}\delta
n+\frac{B_{\alpha\beta}}{n_{0}}\delta
u_{{\alpha\beta}}+\frac{D_{s}}{n_{0}}\delta T$ (54a) $\displaystyle\delta
h_{\alpha\beta}$ $\displaystyle=\frac{B_{\alpha\beta}}{n_{0}}\delta
n+C^{n}_{{\alpha\beta}{\gamma\delta}}\delta
u_{\gamma\delta}+F_{\alpha\beta}\delta T$ (54b) $\displaystyle-\delta s$
$\displaystyle=\frac{D_{s}}{n_{0}}\delta n+F_{\alpha\beta}\delta
u_{{\alpha\beta}}+E\delta T.$ (54c)

for the deviations of the conjugate fields around their equilibrium values,
which are $\mu_{0}$ for the chemical potential, $s^{0}$ for the entropy
density, and $h^{0}_{{\alpha\beta}}=0$ for the stress field conjugate to
strain in this ensemble. From these relations, one can proceed to obtain the
following thermodynamic derivatives and their Maxwell relations

$\displaystyle\dfrac{\partial^{2}f}{\partial
n^{2}}\Big{|}_{u_{{\alpha\beta}},T}=\dfrac{\partial\mu}{\partial
n}\Big{|}_{u_{{\alpha\beta}},T}=\dfrac{A}{n^{2}_{0}}$ (55a)
$\displaystyle\dfrac{\partial^{2}f}{\partial n\partial
u_{\alpha\beta}}=\dfrac{\partial\mu}{\partial
u_{{\alpha\beta}}}\Big{|}_{n,T}=\dfrac{\partial h_{{\alpha\beta}}}{\partial
n}\Big{|}_{u_{\gamma\delta},T}=\dfrac{B_{{\alpha\beta}}}{n_{0}}$ (55b)
$\displaystyle\dfrac{\partial^{2}f}{\partial u_{{\alpha\beta}}\partial
u_{\gamma\delta}}\Big{|}_{n,T}=\dfrac{\partial h_{{\alpha\beta}}}{\partial
u_{\gamma\delta}}\Big{|}_{n,T}=C^{n}_{\alpha\beta\gamma\delta}$ (55c)
$\displaystyle\dfrac{\partial^{2}f}{\partial n\partial
T}=\dfrac{\partial\mu}{\partial
T}\Big{|}_{u_{{\alpha\beta}},n}=-\dfrac{\partial s}{\partial
n}\Big{|}_{u_{\gamma\delta},T}=\dfrac{D_{s}}{n_{0}}$ (55d)
$\displaystyle\dfrac{\partial^{2}f}{\partial
T^{2}}\Big{|}_{u_{\alpha\beta},n}=-\dfrac{\partial s}{\partial
T}\Big{|}_{u_{{\alpha\beta}},n}=E$ (55e)
$\displaystyle\dfrac{\partial^{2}f}{\partial T\partial
u_{\alpha\beta}}=-\dfrac{\partial s}{\partial
u_{{\alpha\beta}}}\Big{|}_{n,T}=\dfrac{\partial h_{{\alpha\beta}}}{\partial
T}\Big{|}_{u_{\gamma\delta},n}=F_{{\alpha\beta}}$ (55f)

The second derivatives of the free energy density provide coefficients of
elasticity and coupling constants between different thermodynamic variables
which can be used in obtaining certain useful measurable quantities like
isothermal compressibility, heat capacity per unit volume at constant volume
and thermal expansion coefficient (see e.g. [26]). One of the aims of this
paper is to provide the microscopic basis to these material constants. In the
subsequent section IV.2, we present the connections between the constants that
we identified in the relations Eq. (55) and the quantities that arise in the
reactive and dissipative terms of the hydrodynamic equations we derived (Eq.
(52)) using the projector Eq. (48) after identifying the relevant dynamical
variables and their respective thermodynamic conjugates(see Eq. (36)).

###  IV.2 Revisiting the Conjugate variables

In section II, when we first introduced the thermodynamic conjugate variables
to the microscopic relevant variables in Eq. (19) and then later re-framed
them in Eq. (43) using the coarse-graining ansatz, we never mentioned the
physical significance of these quantities in the broader context of
thermodynamic parameters. In this section and appendix D we address this
point.

In the appendix D, we connect $(i)$ the change in the density of an extensive
variable (like $n$) induced by the change in its intensive thermodynamic
conjugate (like $\beta\mu$) (see Eq. (113d)) to $(ii)$ correlations between
the densities of the extensive thermodynamic variables (see Eq. (113h)). The
matrix block ${\bm{\chi}_{nue}}$, corresponding to the correlations between
fluctuations in density ($n$), displacement fields ${\bf u}$ with components
$u_{\alpha}$ and energy density $e$, of the coarse-grained susceptibility
matrix $\bm{\chi}$ is identified as matrix of thermodynamic derivatives; see
Eq. (106). These relations follow from evaluating the variances in the
generalised grand-canonical ensemble [23, 6, 11]. The final outcome of
appendix D is the identification of $\delta a$, $\delta{\bf y}$ and $\delta b$
in terms of the intensive thermodynamic quantities showing up in the expansion
of the free energy in Eq. (IV.1). This we summarise here.

$\displaystyle\delta
a=\beta^{-1}\delta\left(\beta\mu\right)=T\delta\left(\frac{\mu}{T}\right)=\delta\mu-\frac{\mu^{0}}{T}\delta
T$ (56a) $\displaystyle\delta y_{{\alpha\beta}}=\beta^{-1}\delta(\beta
h_{{\alpha\beta}})=T\delta\left(\dfrac{h_{{\alpha\beta}}}{T}\right)=\delta
h_{{\alpha\beta}},(\because h^{0}_{{\alpha\beta}}=0)$ (56b)
$\displaystyle\delta
b=-\beta^{-1}\delta\beta=T\delta\left(-\dfrac{1}{T}\right)=\frac{1}{T}\delta
T.$ (56c)

Using the information in Eq. (54a) and Eq. (54b), we rewrite the expressions
for $\delta a$, $\delta y_{{\alpha\beta}}$ and $\delta b$

$\displaystyle\delta a=\delta\mu-\dfrac{\mu_{0}}{T}\delta
T=\dfrac{A}{n^{2}_{0}}\delta n+\dfrac{B_{\alpha\beta}}{n_{0}}\delta
u_{\alpha\beta}+\frac{D_{s}}{n_{0}}\delta T-\mu_{0}\dfrac{\delta T}{T},$ (57a)
$\displaystyle\delta y_{{\alpha\beta}}=\delta
h_{{\alpha\beta}}=\dfrac{B_{\alpha\beta}}{n_{0}}\delta
n+C^{n}_{{\alpha\beta}\gamma\delta}\delta
u_{\gamma\delta}+F_{{\alpha\beta}}\delta T.$ (57b)

and then compare them to Eq. (43) to obtain

$\displaystyle\dfrac{\partial\mu}{\partial
n}\Big{|}_{u_{{\alpha\beta}},T}=\dfrac{A}{n^{2}_{0}}=\dfrac{\nu}{n^{2}_{0}}$
(58a) $\displaystyle\dfrac{\partial\mu}{\partial
u_{{\alpha\beta}}}\Big{|}_{n,T}=\dfrac{B_{{\alpha\beta}}}{n_{0}}=-\dfrac{\mu_{{\alpha\beta}}}{n_{0}}$
(58b) $\displaystyle\dfrac{\partial h_{{\alpha\beta}}}{\partial
u_{\gamma\delta}}\Big{|}_{n,T}=C^{n}_{{\alpha\beta}}=\left(\lambda_{\alpha\gamma\beta\delta}+\lambda_{\alpha\delta\gamma\beta}-\lambda_{\alpha\beta\gamma\delta}\right)$
(58c)

Here we have invoked the set of Eq. (54) and Eq. (55) for the relations in Eq.
(58). First we focus on the coefficients of $\delta n$ and $\delta
u_{{\alpha\beta}}$ and take up the identification of the coefficients to
$\delta T$ separately. The thermodynamic relation given in Eq. (55b) shows
$B_{{\alpha\beta}}$ to be a coupling between the displacement and the density
fields. Previous theoretical perspectives derived by Szamel _et al._ [4, 5]
neglect the contribution of this coupling to the linear elastic response of
crystals with point-defects. However, following the definition given by Walz
_et al._ [10], later studies [11, 15, 14] evaluate this quantity in isothermal
crystalline systems [12, 32, 33, 15] with known direct correlation functions
from classical density functional theory. In this paper, with a more general
treatment of the thermodynamics, we recover (see Eq. (58b)) the definition of
this coupling term while identifying similar cross-correlations of density and
displacements with temperature (see Eq. (55d) and Eq. (55f)). We explain (see
Eq. (61)) how these quantities can be defined in terms of coefficients derived
from microscopic fluctuations.

At this stage, before proceeding to re-consider the interpretations of the
equations of motion (Eq. (52)), we need thermodynamic basis for the terms
$\theta,\tau_{{\alpha\beta}}$ and $L$, that we introduced in our derivations
of the dynamical equations. So we examine the connections between
$\theta,\tau_{{\alpha\beta}},L$ and $D_{s},E,F_{{\alpha\beta}}$ (Eq. (55)).
Following from the first and the second laws of thermodynamics, the relation
between the density of entropy and internal energy for a system at constant
volume and no external strain, is given by [6]

$\displaystyle\frac{1}{T}\delta e=\delta s+\frac{\mu}{T}\delta n.$ (59)

This relation is identical to the ones employed in the context of the
hydrodynamic description of a simple one-component fluid [6]. Now, let us
reconsider the variable $\delta b$ (Eq. (43b)), which now has been established
(Eq. (56c)) as an intensive thermodynamic field. In Eq. (43b), substitute
$\delta b$ and $\delta e$ using Eq. (56c) and Eq. (59) respectively.
Rearranging the substituted equation gives Eq. (60). Compare this expression
for $\delta s$ with Eq. (54c) to obtain relations between microscopically
defined variables like $\theta$, $\tau_{\alpha\beta}$, $L$ and the
thermodynamic coefficients in the expansion of the free energy density (Eq.
(IV.1)).

$\displaystyle-$ $\displaystyle\delta
s=\left(\dfrac{\mu}{T}-\dfrac{\theta}{n_{0}T}\right)\delta
n-\dfrac{L}{T^{2}}\delta T+\dfrac{\tau_{{\alpha\beta}}}{T}\delta
u_{{\alpha\beta}}$ (60)

Equating the coefficients of $\delta n,\delta T$ and $\delta
u_{{\alpha\beta}}$ in the equations Eq. (60) and Eq. (54c) leads to the
following relations

$\displaystyle\theta=\mu n_{0}-TD_{s}$ (61a) $\displaystyle L=-T^{2}E$ (61b)
$\displaystyle\tau_{\alpha\beta}=TF_{\alpha\beta}.$ (61c)

This concludes the thermodynamic interpretation of all the microscopically
defined variables derived and used in the equations of motion in section
III.4.

###  IV.3 Equations of motion

In this section we want to re-examine the hydrodynamic equations (Eq. (52))
with the knowledge of the thermodynamic relations derived in section IV.2.
Inserting the conjugate variables $\delta a$, $\delta b$ and $\delta
y_{\alpha}$ from Eq. (56) into the time evolution equations (Eq. (52)) of the
relevant variables, one obtains

$\displaystyle\partial_{t}\delta n({\bf q},t)$
$\displaystyle=-in_{0}q_{\alpha}\delta v_{\alpha}({\bf q},t)$ (62a)
$\displaystyle\partial_{t}\delta e({\bf q},t)$
$\displaystyle=-\mathrm{i}\left(e_{0}+p_{0}\right)q_{\alpha}\delta
v_{\alpha}({\bf q},t)$
$\displaystyle+q_{\alpha}q_{\gamma}\xi_{\alpha\beta}\delta
h_{\gamma\beta}({\bf q},t)-q_{\alpha}q_{\beta}\alpha_{\alpha\beta}\delta
T({\bf q},t)$ (62b) $\displaystyle\partial_{t}{\delta j_{\alpha}({\bf q},t)}$
$\displaystyle=\mathrm{i}q_{\beta}\delta h_{\alpha\beta}({\bf
q},t)-\mathrm{i}q_{\beta}\delta p({\bf q},t)\delta_{\alpha\beta}$
$\displaystyle-q_{\beta}q_{\gamma}\eta_{{\alpha\beta}{\gamma\delta}}\delta
v_{\delta}({\bf q},t)$ (62c) $\displaystyle\partial_{t}{\delta u_{\alpha}({\bf
q},t)}$ $\displaystyle=\delta v_{\alpha}({\bf q},t)$
$\displaystyle+iq_{\gamma}\zeta_{\alpha\beta}\delta h_{\gamma\beta}({\bf
q},t)-\mathrm{i}q_{\beta}\xi_{\alpha\beta}^{\top}\dfrac{\delta T({\bf
q},t)}{T}$ (62d)

Acquiring this version of the equation for the density of the linear momentum
(Eq. 62d) required further simplification. First, using $p=-e+\mu n+Ts$ and
the Gibbs–Duhem relation, $\delta p=n\delta\mu+s\delta T$, one arrives at

$\displaystyle n_{0}\delta a+(e_{0}+p_{0})\delta b$
$\displaystyle=n_{0}\delta\mu+\frac{e_{0}+p_{0}-\mu n_{0}}{T}\delta T$ (63a)
$\displaystyle=\delta p.$ (63b)

which simplifies the equation of motion of the momentum density given in Eq.
(52)c to the following

$\displaystyle\partial_{t}{\delta j_{\alpha}({\bf q},t)}=$
$\displaystyle-\mathrm{i}q_{\alpha}\delta p({\bf q},t)-\delta y_{\alpha}({\bf
q},t)$ $\displaystyle-
q_{\beta}q_{\gamma}\eta_{{\alpha\beta}{\gamma\delta}}v_{\delta}({\bf q},t).$
(64)

As we additionally identified all the intensive thermodynamic variables
conjugate to their extensive counterparts (see Eqs. (56)), the final
hydrodynamic equations of motion of a real crystal have been obtained starting
from microscopic variables. Now, we can consider their reactive and
dissipative contributions. One of the obvious outcomes of the reactive
couplings between momentum density and number density as well as the
displacement fields, are the coefficients of elasticity. These coefficients,
derived from the microscopic perspective of Mori-Zwanzig formalism has been
discussed in great detail in references [14, 10, 11]. The focus of these
previous studies were the derivation of the isothermal reversible reactive
couplings between the relevant hydrodynamic variables and their proper
validation through implementation in appropriate local-defect rich crystalline
solids. In this contribution, in addition to deriving the microscopic basis
for the dissipative terms, we also obtain the equation of motion of the energy
density (Eq. (62)b).

###  IV.4 Implications of our results

The phenomenon of heat transport in a solid medium at finite temperature is
observed to be diffusive. Consistent with existing macroscopic theory
perspectives [3, 20], our equations yield diffusive heat transport with two
separate transport coefficients $\xi_{\alpha\beta},\alpha_{{\alpha\beta}}$
(see Eq. (66)c) and their respective microscopic definitions (Eq. (49)c and
Eq. (III.3)). For an equilibrated fluid or a crystalline solid with constant
point-defect concentration and in absence of any external deforming fields,
one can define a variable $\hat{q}({\bf r},t)=\hat{e}({\bf
r},t)-\dfrac{e_{0}+p_{0}}{n_{0}}\hat{n}({\bf r},t)$ whose fluctuation, in the
thermodynamic limit, is

$\displaystyle\delta q=\delta e-\dfrac{e_{0}+p_{0}}{n_{0}}\delta n.$ (65)

It can be shown [6, 31] that $\delta q$ is related to the entropy $S$ of the
system through the relation $\delta q=Tn_{0}\delta(S/N)$. Thus, following the
second law of thermodynamics, $q$ can be interpreted as the heat density in
the system. Rearranging the terms in the Eq. (62)b for the energy density

$\displaystyle\partial_{t}\delta e({\bf
q},t)+\mathrm{i}\left(e_{0}+p_{0}\right)q_{\alpha}\delta v_{\alpha}({\bf
q},t)=q_{\alpha}q_{\beta}\xi_{\alpha\beta}\delta h_{\alpha\beta}({\bf
q},t)-q_{\alpha}q_{\beta}\alpha_{\alpha\beta}\delta T({\bf q},t)$ (66a)
$\displaystyle\implies$ $\displaystyle\partial_{t}\delta e({\bf
q},t)-\dfrac{e_{0}+p_{0}}{n_{0}}\partial_{t}\delta n({\bf
q},t)=q_{\alpha}q_{\beta}\xi_{\alpha\beta}\delta h_{\alpha\beta}({\bf
q},t)-q_{\alpha}q_{\beta}\alpha_{\alpha\beta}\delta T({\bf q},t)$ (66b)
$\displaystyle\implies$ $\displaystyle\partial_{t}\delta
q=q_{\alpha}q_{\beta}\xi_{\alpha\beta}\delta h_{\alpha\beta}({\bf
q},t)-q_{\alpha}q_{\beta}\alpha_{\alpha\beta}\delta T({\bf q},t)$ (66c)

where $\delta v_{\alpha}({\bf q})$ was substituted using the mass conservation
equation Eq. (62)a, one ends up with the two dissipative contributions (Eq.
(66))c to the diffusion of heat in a crystal with point defects. Here we point
out the distinctive aspects of the dissipative terms which distinguish a
crystalline solid from a fluid. The analogue to the dissipative coefficient
tensor $\alpha_{{\alpha\beta}}$ in a fluid is a scalar representing heat
conductivity. In case of the crystalline solid $\alpha_{{\alpha\beta}}$ is a
second rank tensor whose symmetry is determined by the symmetry of the
concerned crystalline structure. Moreover, the additional dissipative coupling
term $\xi_{\alpha\beta}$, arising due to the coupling between energy and
displacement fields, appears in the time evolution equation for energy (Eq.
(62)b). Because of Onsager symmetry, its transpose $\xi_{\alpha\beta}^{T}$,
appears in the time evolution equation for the displacement fields Eq. (62)d.
The tensor $\xi_{\alpha\beta}$ is characteristic for a system with long range
order.

In a crystalline solid, the two diffusive hydrodynamic modes are heat
transport and diffusion of point defects [2, 3]. Remaining consistent with our
previous contributions [10, 11, 14], we define the field of fluctuating point-
defect concentration as

$\displaystyle\delta c({\bf q},t)=-\delta n({\bf
q},t)-n_{0}\mathrm{i}q_{\alpha}\delta u_{\alpha}({\bf q},t).$ (67)

Taking the time derivative of Eq. (67), then substituting $\partial_{t}\delta
n({\bf q},t)$ and $\partial_{t}\delta u_{\alpha}({\bf q},t)$ using Eq. (62)a
and Eq. (62)d immediately shows that the ansatz Eq. (67) leads to the
following

$\displaystyle\partial_{t}\delta c({\bf
q},t)=n_{0}q_{\alpha}q_{\gamma}\zeta_{\alpha\beta}\delta h_{\gamma\beta}({\bf
q},t)-n_{0}q_{\alpha}q_{\beta}\xi_{\alpha\beta}^{\top}\dfrac{\delta T({\bf
q},t)}{T}$ (68)

equation for the time evolution of the point-defect concentration. Eq. (66)c
and Eq. (68) show that the leading order contributions to the time evolution
equations for the heat density and the point-defect density are of the order
${\bf q}^{2}$. This is consistent with the phenomenological theories [3, 20]
and predicts the transport processes for heat and point-defects in a
crystalline solid to be diffusive.

The definition of the fluctuation of the point-defect concentration (Eq. (67))
will help to rationalize how the present, more general, perspective reduces to
our previous microscopic dissipation-less description of isothermal elastic
properties of a point-defect rich crystal [11]. Substitution of the density
fluctuation $\delta n({\bf q},t)$ in terms of the defect density concentration
$\delta c({\bf q},t)$ and the displacement fields $\delta u_{\alpha}({\bf
q},t)$ in Eq. (121)a, followed by some mathematical manipulations [26],
results in the following equation for the density of linear momentum

$\displaystyle\partial_{t}\delta j_{\alpha}({\bf q},t)=$
$\displaystyle-\Lambda_{{\alpha\beta}}({\bf q})\delta u_{\beta}({\bf q},t)$
$\displaystyle-V_{\alpha}({\bf q})\delta c({\bf
q},t)-\dfrac{Z^{*}_{\alpha}({\bf q})}{T}\delta T({\bf q},t)$ (69)

where

$\displaystyle\Lambda_{a\beta}({\bf q})$
$\displaystyle=\lambda_{{\alpha\beta}}({\bf q})-iq_{\alpha}\mu_{\beta}({\bf
q})+iq_{\beta}\mu^{*}_{\alpha}({\bf q})+q_{\alpha}\nu({\bf q})q_{\beta}$
$\displaystyle\approx\Lambda_{{\alpha\beta}\gamma\delta}q_{\beta}q_{\delta}$
(70a) $\displaystyle\Lambda_{{\alpha\beta}\gamma\delta}$
$\displaystyle=\lambda_{\alpha\gamma\beta\delta}+\lambda_{\alpha\delta\gamma\beta}-\lambda_{\alpha\beta\gamma\delta}$
$\displaystyle+\delta_{\alpha\beta}\mu_{\gamma\delta}+\mu_{{\alpha\beta}}\delta_{\gamma\delta}+\nu\delta_{{\alpha\beta}}\delta_{\gamma\delta}$
(70b) $\displaystyle V_{\alpha}({\bf q})$
$\displaystyle=\dfrac{1}{n_{0}}\left[\mu^{*}_{\alpha}({\bf
q})-iq_{\alpha}\nu({\bf q})\right]$
$\displaystyle\approx\dfrac{1}{n_{0}}\left[-i\mu_{{\alpha\beta}}q_{\beta}-i\nu
q_{\alpha}\right]$ (70c) $\displaystyle Z^{*}_{\alpha}({\bf q})$
$\displaystyle=\tau^{*}_{\alpha}({\bf
q})+iq_{\alpha}\left[e_{0}+p_{0}-\theta^{*}({\bf q})\right]$
$\displaystyle\approx-i\tau_{{\alpha\beta}}q_{\beta}+iq_{\alpha}[e_{0}+p_{0}-\theta].$
(70d)

The second lines in the expressions for $\Lambda_{{\alpha\beta}}({\bf q})$
(Eq. (70a)), $V_{\alpha}({\bf q})$ (Eq. (70c)) and $Z^{*}_{\alpha}({\bf q})$
(Eq. (70d)) represent the leading terms in the ${\bf q}\rightarrow 0$ limit of
the respective quantities. Since the aim is to show how to recover our
previous results from a more general description, we have chosen to keep only
the reversible parts in Eq. (IV.4). In this limit, Eq. (62d) reduces to
$\partial_{t}\delta\mathbf{u}=\mathbf{v}$ in which $\mathbf{v}$ is related to
${\bf j}$ through the relation in Eq. (39d). If the system is isothermal,
setting $\delta T$ to zero, tive-derivative of Eq. (IV.4) leads to the
isothermal wave equation for the linear momentum, first derived in reference
[10]. Moreover, the quantity $\Lambda_{{\alpha\beta}}({\bf q})$ is identical
to the dynamical matrix associated with the isothermal elastic properties of
the crystal with constant defect concentration. Thus the definitions of the
quantities $\Lambda_{{\alpha\beta}}({\bf q})$ and $V_{\alpha}({\bf q})$ given
in Eq. (70) have the same interpretations as in references [10, 11, 14].

If one chooses to write the reversible part of the density of the linear
momentum while considering $\delta q({\bf q},t)$ and $\delta c({\bf q},t)$ as
the relevant variables, the following equation is obtained

$\displaystyle\partial_{t}\delta j_{\alpha}({\bf q},t)$
$\displaystyle=-\Lambda^{ad}_{{\alpha\beta}}({\bf q})\delta u_{\beta}({\bf
q},t)$ $\displaystyle-V^{ad}_{\alpha}({\bf q})\delta c({\bf
q},t)-L^{-1}Z^{*}_{\alpha}({\bf q})\delta q$ (71)

where

$\displaystyle\Lambda^{ad}_{a\beta}({\bf q})$
$\displaystyle=\Lambda_{{\alpha\beta}}({\bf q})+L^{-1}({\bf
q})Z^{*}_{\alpha}({\bf q})Z_{\beta}({\bf q})$
$\displaystyle\approx\Lambda_{{\alpha\beta}\gamma\delta}q_{\beta}q_{\delta}$
$\displaystyle+L^{-1}(\tau_{\alpha\gamma}q_{\gamma}-q_{\alpha}[e_{0}+p_{0}-\theta])(\tau_{\beta\delta}q_{\delta}-q_{\beta}[e_{0}+p_{0}-\theta])$
(72a) $\displaystyle V^{ad}_{\alpha}({\bf q})$ $\displaystyle=V_{\alpha}({\bf
q})-L^{-1}({\bf q})Z^{*}_{\alpha}({\bf
q})\left(\dfrac{e_{0}+p_{0}-\theta^{*}({\bf q})}{n_{0}}\right)$
$\displaystyle\approx\dfrac{-i}{n_{0}}(\mu_{{\alpha\beta}}q_{\beta}-\nu
q_{\alpha})$
$\displaystyle+\dfrac{i}{n_{0}}L^{-1}\left(\tau_{{\alpha\beta}}q_{\beta}-q_{\alpha}[e_{0}+p_{0}-\theta]\right)(e_{0}+p_{0}-\theta).$
(72b)

Examining the ${\bf q}\rightarrow 0$ limit of $V^{ad}_{\alpha}({\bf q})$ and
$L^{-1}Z^{*}_{\alpha}({\bf q})$ reveals that for both of them, the leading
order term is $\mathcal{O}({\bf q})$. This becomes clear from the ${\bf
q}\rightarrow 0$ limit, given in Appendix C, of the constituent variables
$\mu_{\alpha}({\bf q}),\nu({\bf q}),\tau_{\alpha}({\bf q}),\theta({\bf q})$ of
these two terms. From Eq. (66)c and Eq. (68) we know that time derivatives of
$\delta q$ and $\delta c$ has leading order $\mathcal{O}({\bf q}^{2})$
contributions. Therefore, time derivative of the momentum density equation Eq.
(IV.4) yields $\mathcal{O}({\bf q}^{2})$ contributions from the first term
while the ${\bf q}-$dependence is of the order $\mathcal{O}({\bf q}^{3})$ for
the terms associated with $\partial_{t}\delta q({\bf q},t)$ and
$\partial_{t}\delta c({\bf q},t)$. Ignoring the terms with higher order ${\bf
q}-$dependence, after substituting $\partial_{t}\delta u_{\alpha}({\bf q},t)$
with Eq. (62)d, we arrive at the wave equation

$\displaystyle\partial^{2}_{t}\delta j_{\alpha}({\bf
q},t)=-n^{-1}_{0}\Lambda^{ad}_{{\alpha\beta}}({\bf q})\delta j_{\beta}({\bf
q},t).$ (73)

The matrix $\Lambda^{ad}_{{\alpha\beta}}({\bf q})$, defined in Eq. (72)a, is
the adiabatic dynamical matrix [34] for the crystalline solid. It is
associated with all the adiabatic elastic coefficients and governs the
longitudinal and transverse speeds of sound in the crystalline solid under
adiabatic conditions where heat and local defects can adjust freely.

In defining the viscosity tensor $\eta_{{\alpha\beta}\gamma\delta}$ (Eq.III.3)
we have not decomposed the longitudinal and transverse contributions of the
stress fluctuations. In case of a simple fluid, the attenuation of the
propagative sound waves associated with longitudinal components of the current
correlation functions [35], is shown [6] to be related to the viscosity
tensor. Components of the viscosity tensor are also shown to govern the
diffusion of the transverse shear waves associated with the transverse
components of the current correlation functions. There, the reactive part of
the time evolution equation of the momentum density is given by the pressure
gradient alone, the second term in Eq. (62)c. Our calculations reveal that, in
case of a crystalline solid, an additional reactive contribution arises in the
equation of motion of the momentum density. From the literature [3, 20] on the
macroscopic hydrodynamics of crystalline solids it is well known that this
indicates the emergence of propagating shear waves. This marks another
important difference between a system of simple fluid and a solid with long
range order. A completely microscopic definition of the viscosity tensor
proves to be another important result of this paper.

##  V Conclusions and Outlook

Mechanical properties and transport coefficients, in a crystalline phase with
a finite concentration of point-defects, are derived in this paper from a
completely microscopic classical description. The hydrodynamic equations of
the three components of the displacement vectors and the variables governing
the local conservation of mass, energy, and linear momentum are derived. Our
results are based on the Zwanzig-Mori formalism, which requires a choice of
the considered variables, and on an ansatz for the microscopic density
fluctuations in terms of the coarse-grained fields, Eq. (34). The choice of
the relevant variables is dictated by conservation laws and spontaneous
symmetry breaking. The ansatz for the displacement field (Eq. (35b)), first
suggested by Szamel and Ernst [4] has been tested in hard spheres [15] and
cluster crystals [14]. How this hydrodynamic description relates to the
equilibrium thermodynamic properties like free energy and entropy was
rigorously worked out too. Our approach allows us to derive the transport
coefficients associated with the Green-Kubo relations. The rank of the
transport coefficient tensors reflect the crystal symmetries and their
explicit microscopic origin provides the possibility of evaluating them from
inputs obtained from atomistic simulations of appropriate systems. Theoretical
and simulation studies in reference [14] show how theoretical frameworks
derived in references [10, 11] can be successfully implemented to obtain
quantitative insights regarding mechanical response of solids. Though the
previous studies [11, 14] were restricted to reversible mechanical response in
isothermal crystals, the theories success in predicting elastic properties of
a model known to mimic crystalline phases in DNA-based dendritic
nanostructures [13], illustrates the scope of applicability of this
formulation.

This paper provides a complete microscopic formulation of all the constants of
linear elasticity and coefficients governing diffusive and wave transport
processes for crystals. It has been established [11, 14, 15], albeit in a much
more simplified purely reversible version of this theory, that the predictions
can be validated using inputs from atomistic simulations of crystalline phases
of particles interacting via ultra-soft potentials [36, 13]. Recent advances
in experimental techniques allow the tracking of single particle dynamics in a
thermodynamic ensemble of soft materials [37, 38]. This offers the exiting
possibility of understanding the emergent macroscopic properties from
experimentally observable microscopic dynamics through the lens of a
statistical mechanics theory derived from first principles.

From a macroscopic phenomenology perspective, heat transport is an
irreversible process related to the production of entropy in the system.
Reference [39], a study of meta-stable polycrystals by some of the authors,
presents a completely different point of view centered around the relation of
entropy to statistical microstates. In this case our extensive atomistic
simulations give us insights about the relative stability of topological-
defect rich polycrystalline configurations and allows us to discern the role
of entropy in it. Therefore, combining the theoretical phenomenology of this
paper and its predecessors [10, 11, 14] with simulations similar to [39], can
lead to fundamentally new ways of interpreting reversible or irreversible
thermo-mechanical response of crystalline solids and understanding transport
processes driven by gradients in intensive thermodynamic fields like
temperature.

Defining the displacement fields unambiguously, in a system with long range
periodicity and diffusing defects, proves to be one of the conceptual
challengers to phenomenological or atomistic theories attempting to explain
macroscopic mechanical response [40]. Our approach of representing the system
in reciprocal space, resolves this problem by defining displacements, particle
number density or defect density in terms of density fluctuations close to
Bragg-peaks. A natural extension of this would be to examine the possibility
of employing similar principles to study systems subjected to large
deformations near the onset of plasticity.

Some recent research employs spatial projection operators [41, 42, 43] to
segregate microscopic displacement fluctuations associated with macroscopic
elastic or plastic response in defect-free crystals. This interpretation, of
displacement fields, helps in explaining [44] the origin of rigidity and the
shear rate dependence of the yield-point [45] in an ideal crystal. In this
paper, the projection operators derive the dynamics of variables chosen
because of their slow relaxation time-scales. The propagating longitudinal or
transverse sound modes, described here, exhibit dispersion relations that
vanish linearly with decreasing wave-vector. The slopes of their dispersion
relations in the small wave-number limit also give the moduli of linear
elastic response encoded in the dynamical matrix appearing in the wave
equation of the displacement fields [10]. The spatial projection of reference
[41] may present a way to separate the contribution of affine macroscopic
deformations and local particle motions to the total microscopic displacement
fluctuations in a defect rich crystal. Therefore, one of the future avenues
for investigation will be an attempt to understand the onset of plasticity
through the convergence of these perspectives.

Identifying connections, between the current theoretical framework and
reversible or irreversible isothermal and adiabatic processes in the crystal
with local defects, will pave the path for future endeavours to evaluate
macroscopic mechanical constants and transport coefficients of materials of
theoretical and practical interest.

###### Acknowledgements.

Discussions with Gerhard Kahl, Martin Oettel, Thomas Franosch, and Johannes
Häring are gratefully acknowledged. This work is supported by Deutsche
Forschungsgemeinschaft through grant FU 309/11-1.

##  Appendix A Some conventions and definitions

###  A.1 The microscopic variables and their Fourier transforms

Assuming an interaction potential $V(|{\bf r}_{i}-{\bf r}_{j}|)=V(r_{ij})$
dependent on the distance ${\bf r}_{ij}$ between particles $i$ and $j$, in
Fourier space, the microscopic stress tensor and the microscopic energy
current reads,

$\hat{\boldsymbol{\sigma}}(\mathbf{q},t)=-\sum_{i}\frac{\hat{\bf p}({\bf
r}_{i},t)\hat{\bf p}({\bf
r}_{i},t)}{m}e^{-i\mathbf{q}\cdot\mathbf{r}_{i}(t)}+\frac{1}{2}\sum_{i\neq
j}\frac{\mathbf{r}_{ij}(t)\mathbf{r}_{ij}(t)}{r_{ij}(t)}V^{\prime}(r_{ij}(t))\frac{e^{-i\mathbf{q}\cdot\mathbf{r}_{j}(t)}-e^{-i\mathbf{q}\cdot\mathbf{r}_{i}(t)}}{i\mathbf{q}\cdot\mathbf{r}_{ij}(t)},$
(74)

$\displaystyle\hat{\mathbf{j}}^{e}(\mathbf{q},t)$
$\displaystyle=\sum_{i}E({\bf r}_{i},t)\frac{\hat{\mathbf{p}}({\bf
r}_{i},t)}{m}e^{-i\mathbf{q}\cdot\mathbf{r}_{i}(t)}-\frac{1}{4}\sum_{i\neq
j}\frac{\hat{\mathbf{p}}({\bf r}_{i},t)+\hat{\mathbf{p}}({\bf
r}_{j},t)}{m}\cdot\frac{\mathbf{r}_{ij}(t)\mathbf{r}_{ij}(t)}{r_{ij}(t)}V^{\prime}(r_{ij}(t))\frac{e^{-i\mathbf{q}\cdot\mathbf{r}_{j}(t)}-e^{-i\mathbf{q}\cdot\mathbf{r}_{i}(t)}}{i\mathbf{q}\cdot\mathbf{r}_{ij}(t)}$
(75)

And one can indeed verify that [46]

$\int
d\mathbf{r}e^{-i\mathbf{q}\cdot\mathbf{r}}\int_{0}^{1}ds\,\delta(\mathbf{r}-\mathbf{r}_{i}+s\mathbf{r}_{ij})=\frac{e^{-i\mathbf{q}\cdot\mathbf{r}_{j}}-e^{-i\mathbf{q}\cdot\mathbf{r}_{i}}}{i\mathbf{q}\cdot\mathbf{r}_{ij}}.$
(76)

In the microscopic definitions for the energy density (Eq. (15)), stress (Eq.
(74)) and energy current (Eq. (75)), the per particle mass $m$ appears. This
however, is set to one without any loss of generality and therefore it does
not appear in any other equations in the rest of the paper.

###  A.2 Equilibrium values

We can define an equilibrium pressure, using the virial equation for the
average pressure

$p_{0}=n_{0}k_{\text{B}}T-\frac{1}{6V}\langle\sum_{i\neq
j}r_{ij}V^{\prime}(r_{ij})\rangle.$ (77)

The last quantity that is not equally zero at equilibrium is the energy
density,

$\displaystyle\langle\hat{e}(\mathbf{r})\rangle$
$\displaystyle=\langle\sum_{i}E_{i}\delta(\mathbf{r}-\mathbf{r}_{i})\rangle$
$\displaystyle=\langle\sum_{i}\frac{\hat{\mathbf{p}}_{i}^{2}}{2m}\delta(\mathbf{r}-\mathbf{r}_{i})\rangle+\frac{1}{2}\langle\sum_{i\neq
j}V(r_{ij})\delta(\mathbf{r}-\mathbf{r}_{i})\rangle$
$\displaystyle=\frac{3}{2}k_{B}Tn(\mathbf{r})+\frac{1}{2}\langle\sum_{i\neq
j}V(r_{ij})\delta(\mathbf{r}-\mathbf{r}_{i})\rangle.$ (78)

Non-locally, taking an ensemble and volume average, we get

$e_{0}=\frac{1}{V}\int d\mathbf{r}\,\langle e(\mathbf{r})\rangle\\\
=\frac{3}{2}n_{0}k_{B}T+\frac{1}{2V}\langle\sum_{i\neq j}V(r_{ij})\rangle.$
(79)

In reciprocal space, they translate to

$e_{0}=\frac{1}{V}\langle\hat{e}(\mathbf{q}=0)\rangle\\\
=\frac{3}{2}n_{0}k_{B}T+\frac{1}{2V}\langle\sum_{i\neq j}V(r_{ij})\rangle.$
(80)

###  A.3 Matrix components

The momentum-momentum density correlation is straightforward with the
equipartition principle,

$\langle\delta\hat{\mathbf{j}}^{*}(\mathbf{q})\delta\hat{\mathbf{j}}(\mathbf{q})\rangle=n_{0}k_{B}TV\mathbb{I}$
(81)

With $\mathcal{L}=-i\frac{d}{dt}$ and Eq.9, we get the first element of the
frequency matrix

$\displaystyle\bm{\omega}^{j\rho}_{{\bf
g}\alpha}=\beta\langle\delta\hat{j}^{*}_{\alpha}(\mathbf{q})\mathcal{L}\delta\hat{\rho}_{{\bf
g}}(\mathbf{q})\rangle=-V(g+q)_{\alpha}n_{{\bf g}}$ (82a)
$\displaystyle\bm{\omega}^{\rho j}_{{\bf
g}\alpha}=\beta\langle\delta\hat{j}^{*}_{\alpha}(\mathbf{q})\mathcal{L}\delta\hat{\rho}_{{\bf
g}}(\mathbf{q})\rangle^{*}=-V(g+q)_{\alpha}n^{*}_{{\bf g}}$ (82b)

The second one reads

$\displaystyle\beta^{-1}\bm{\omega}^{je}=\beta^{-1}\bm{\omega}^{ej}=\langle\delta\mathbf{j}^{*}(\mathbf{q})\mathcal{L}\delta
e(\mathbf{q})\rangle=-{\bf
q}\cdot\langle\delta\hat{\mathbf{j}}^{*}(\mathbf{q})\hat{\mathbf{j}}^{e}(\mathbf{q})\rangle$
(83a) $\displaystyle=-\mathbf{q}\cdot\langle\sum_{i,k}E({\bf
r}_{i})\frac{\hat{\mathbf{p}}_{i}}{m}\hat{\mathbf{p}}_{k}e^{-i\mathbf{q}\cdot\mathbf{r}_{ik}}\rangle+\mathbf{q}\cdot\langle\frac{1}{4}\sum_{k}\sum_{i\neq
j}\frac{\hat{\mathbf{p}}_{i}+\hat{\mathbf{p}}_{j}}{m}\cdot\frac{\mathbf{r}_{ij}\mathbf{r}_{ij}}{r_{ij}}\hat{\mathbf{p}}_{k}V^{\prime}(r_{ij})\frac{e^{-i\mathbf{q}\cdot\mathbf{r}_{j}}-e^{-i\mathbf{q}\cdot\mathbf{r}_{i}}}{i\mathbf{q}\cdot\mathbf{r}_{ij}}e^{i\mathbf{q}\cdot\mathbf{r}_{k}}\rangle$
(83b)
$\displaystyle=-\mathbf{q}\cdot\langle\sum_{i,k}\frac{\hat{\mathbf{p}}_{i}^{2}}{2m}\frac{\hat{\mathbf{p}}_{i}}{m}\hat{\mathbf{p}}_{k}e^{-i\mathbf{q}\cdot\mathbf{r}_{ik}}\rangle-\mathbf{q}\cdot\langle\sum_{k}\sum_{i\neq
j}V(r_{ij})\frac{\hat{\mathbf{p}}_{i}}{2m}\hat{\mathbf{p}}_{k}e^{-i\mathbf{q}\cdot\mathbf{r}_{ik}}\rangle+\mathbf{q}\cdot\langle\frac{k_{B}T}{2}\sum_{i\neq
j}\frac{\mathbf{r}_{ij}\mathbf{r}_{ij}}{r_{ij}}V^{\prime}(r_{ij})\frac{\sin(\mathbf{q}\cdot\mathbf{r}_{ij})}{\mathbf{q}\cdot\mathbf{r}_{ij}}\rangle$
(83c)
$\displaystyle=-\mathbf{q}\frac{5}{2}N\left(k_{B}T\right)^{2}-\mathbf{q}\frac{1}{2}k_{B}T\langle\sum_{i\neq
j}V(r_{ij})\rangle+\mathbf{q}\frac{k_{B}T}{6}\langle\sum_{i\neq
j}r_{ij}V^{\prime}(r_{ij})\rangle+\mathcal{O}(\mathbf{q}^{3})$ (83d)
$\displaystyle=-\mathbf{q}k_{B}TV\left(e_{0}+p_{0}\right)+\mathcal{O}(\mathbf{q}^{3})$
(83e)

with at equilibrium
$r_{\alpha}r_{\beta}=\dfrac{1}{3}r^{2}\delta_{\alpha\beta}$, and

$\displaystyle\langle\hat{\mathbf{p}}_{i}\hat{\mathbf{p}}_{k}\rangle=mk_{B}T\delta_{ik}\mathbb{I},\
\ \langle\hat{\mathbf{p}}_{i}^{2}\rangle=3mk_{B}T,$
$\displaystyle\langle\hat{\mathbf{p}}_{i}^{2}\hat{\mathbf{p}}_{i}\hat{\mathbf{p}}_{k}\rangle=5\left(mk_{B}T\right)^{2}\delta_{ik}\mathbb{I},$
$\displaystyle\langle\hat{\mathbf{p}}_{i}^{2}\hat{\mathbf{p}}_{k}\hat{\mathbf{p}}_{k}\rangle=\left(3+2\delta_{ik}\right)\left(mk_{B}T\right)^{2}\mathbb{I}.$
(84)

##  Appendix B Inversion of the static susceptibility matrix

In the following matrix identities, the matrix blocks $\bm{A,B,C,D}$ are
square matrices. We use these identities (Eq. (85)) to perform the inversion
of the static susceptibility matrix $\bm{\chi}$ (Eq. (86)).

$\displaystyle\begin{bmatrix}\bm{A}&0\\\
0&\bm{B}\end{bmatrix}^{-1}=\begin{bmatrix}\bm{A}^{-1}&0\\\
0&\bm{B}^{-1}\end{bmatrix},\ \ \begin{bmatrix}\bm{A}&\bm{B}\\\
\bm{C}&\bm{D}\end{bmatrix}^{-1}=\begin{bmatrix}\bm{A}^{-1}+\bm{A}^{-1}\bm{B}\left(\bm{D}-\bm{C}\bm{A}^{-1}\bm{B}\right)^{-1}\bm{C}\bm{A}^{-1}&-\bm{A}^{-1}\bm{B}\left(\bm{D}-\bm{C}\bm{A}^{-1}\bm{B}\right)^{-1}\\\
-\left(\bm{D}-\bm{C}\bm{A}^{-1}\bm{B}\right)^{-1}\bm{C}\bm{A}^{-1}&\left(\bm{D}-\bm{C}\bm{A}^{-1}\bm{B}\right)^{-1}\end{bmatrix}$
(85) $\displaystyle\bm{\chi}({\bf q})$
$\displaystyle=\begin{bmatrix}\bm{\chi}^{\rho\rho}_{(N\times
N)}&\bm{\chi}^{\rho e}_{(N\times 1)}&\bm{\chi}^{\rho j}_{(N\times 3)}\\\
\bm{\chi}^{e\rho}_{(1\times N)}&\bm{\chi}^{ee}_{(1\times
1)}&\bm{\chi}^{ej}_{(1\times 3)}\\\ \bm{\chi}^{j\rho}_{(3\times
N)}&\bm{\chi}^{je}_{(3\times 1)}&\bm{\chi}^{jj}_{(3\times
3)}\end{bmatrix}=\begin{bmatrix}\bm{\chi}^{\rho\rho}_{(N\times
N)}&\bm{\chi}^{\rho e}_{(N\times 1)}&0\\\ \bm{\chi}^{e\rho}_{(1\times
N)}&\bm{\chi}^{ee}_{(1\times 1)}&0\\\ 0&0&\bm{\chi}^{jj}_{(3\times
3)}\end{bmatrix}$ (86)

Reiterating Eq. (20) and noting the block diagonal structure of $\bm{\chi}$
(Eq. (86)) first the individual diagonal blocks are inverted using identity
Eq. (85)a. The first diagonal block, comprising of correlation matrices
$\bm{\chi}^{\rho\rho}$, $\bm{\chi}^{\rho e}$ and $\bm{\chi}^{ee}$ is inverted
using identity Eq. (85)b. The second diagonal block, with correlation between
components of momentum densities, is easier to invert using the classical
equipartition (Eq. (21)). Thus the explicit expressions for the inverse matrix

$\displaystyle\chi^{-1}({\bf q})=\begin{bmatrix}({\bf J}^{\rho\rho}+{\bf
UL^{-1}U^{*}})&-{\bf L^{-1}U}&0\\\ -({\bf L^{-1}U)^{*}}&{\bf L^{-1}}&0\\\
0&0&n_{0}^{-1}\mathbb{I}.\end{bmatrix}$ (87)

is acquired where ${\bf J}^{\rho\rho}=(\bm{\chi}^{\rho\rho})^{-1}$ with
components $J_{{\bf g}{\bf g}^{\prime}}$ (Eq. (23)). The scalar $L(q)$ and the
components $U_{{\bf g}}$ of the vector $U$ are given in Eq. (24).

##  Appendix C Small wavelength limit of coefficients related to the inverse
density correlation function and hence the direct correlation function

###  C.1 The generalised elastic coefficients

The elastic coefficients, which are the same as the ones in [10, 14], are
summarised

$\displaystyle\lambda_{\alpha\beta}({\bf q})=\sum_{{\bf g}^{\prime}{\bf
g}^{\prime\prime}}g^{\prime}_{\alpha}n^{*}_{{\bf
g}^{\prime}}J^{*\rho\rho}_{{\bf g}^{\prime}{\bf g}^{\prime\prime}}n_{{\bf
g}^{\prime\prime}}g^{\prime\prime}_{\beta}=\lambda_{\alpha\beta\gamma\delta}q_{\gamma}q_{\delta}+\dots$
(88a) $\displaystyle\mu_{\alpha}({\bf q})=\sum_{{\bf g}^{\prime}{\bf
g}^{\prime\prime}}n^{*}_{{\bf g}^{\prime}}J^{*\rho\rho}_{{\bf g}^{\prime}{\bf
g}^{\prime\prime}}n_{{\bf
g}^{\prime\prime}}ig^{\prime\prime}_{\alpha}=i\mu_{\alpha\beta}q_{\beta}+\dots$
(88b) $\displaystyle\mu^{*}_{\alpha}({\bf q})=\sum_{{\bf g}^{\prime}{\bf
g}^{\prime\prime}}-ig^{\prime}_{\alpha}n^{*}_{{\bf
g}^{\prime}}J^{*\rho\rho}_{{\bf g}^{\prime}{\bf g}^{\prime\prime}}n_{{\bf
g}^{\prime\prime}}=-i\mu_{\alpha\beta}q_{\beta}+\dots$ (88c)
$\displaystyle\nu({\bf q})=\sum_{{\bf g}^{\prime}{\bf
g}^{\prime\prime}}n^{*}_{{\bf g}^{\prime}}J^{*\rho\rho}_{{\bf g}^{\prime}{\bf
g}^{\prime\prime}}n_{{\bf g}^{\prime\prime}}=\nu+\dots$ (88d)

and then derived here for the sake of completeness.

Substituting $J^{*}_{{\bf g}^{\prime}{\bf g}}$ using Eq. (23) in the
expression for $\lambda_{\alpha\beta}(\bf q)$ in Eq. (88)a and utilising the
expansion of the gradient of the average density distribution

$\nabla_{\alpha}n({\bf r})=\sum_{{\bf g}}ig_{\alpha}n_{{\bf g}}e^{i{\bf
g}.{\bf r}}$ (89)

in terms of the Bragg peak amplitudes $n_{{\bf g}}$, one obtains

$\displaystyle\lambda_{\alpha\beta}({\bf q})=$
$\displaystyle\frac{k_{B}T}{V}\int d^{3}r_{1}\int
d^{3}r_{2}\nabla_{\alpha}n({\bf r}_{1})\nabla_{\beta}n({\bf r}_{2})e^{-i{\bf
q}.({\bf r}_{1}-{\bf r}_{2})}$ $\displaystyle\left[\frac{\delta({\bf
r}_{1}-{\bf r}_{2})}{n({\bf r}_{1})}-c({\bf r}_{1},{\bf r}_{2})\right].$ (90)

Upon using an equation derived by Lovett, Mou, Buff, Wertheim (LMB [27]W
[28]),

$\frac{\nabla_{\alpha}(n({\bf r}))}{n({\bf r})}=\int d^{3}r^{\prime}c({\bf
r},{\bf r}^{\prime})\nabla_{\alpha}n({\bf r}^{\prime})$ (91)

and realising that the gradient of the equilibrium density
$\nabla_{\alpha}n({\bf r})$ is real i.e., $\sum_{{\bf g}}ig_{\alpha}n_{{\bf
g}}e^{i{\bf g}.{\bf r}}=\sum_{{\bf g}}-ig_{\alpha}n^{*}_{{\bf g}}e^{-i{\bf
g}.{\bf r}}$, one gets

$\displaystyle\lambda_{\alpha\beta}({\bf q})=$
$\displaystyle\frac{k_{B}T}{V}\int d^{3}r_{1}\int
d^{3}r_{2}\nabla_{\alpha}n({\bf r}_{1})\nabla_{\beta}n({\bf r}_{2})c({\bf
r}_{1},{\bf r}_{2})$ $\displaystyle\left(1-e^{-i{\bf q}\cdot({\bf r}_{1}-{\bf
r}_{2})}\right)$ (92a) $\displaystyle\approx$
$\displaystyle\lambda_{\alpha\beta\gamma\delta}q_{\gamma}q_{\delta}+\mathcal{O}(q^{4})$
(92b)

Similar arguments lead to the expression for

$\displaystyle\mu_{\alpha}({\bf q})=$ $\displaystyle\frac{k_{B}T}{V}\int
d^{3}r_{1}\int d^{3}r_{2}n({\bf r}_{1})\nabla_{\alpha}n({\bf r}_{2})c({\bf
r}_{1},{\bf r}_{2})$ $\displaystyle\left(1-e^{-i{\bf q}\cdot({\bf r}_{1}-{\bf
r}_{2})}\right)$ (93a) $\displaystyle\approx$ $\displaystyle
i\mu_{\alpha\beta}q_{\beta}+\mathcal{O}(q^{2})$ (93b)

It can, however, be shown [14] that for a crystal with inversion symmetry, the
correction in the small ${\bf q}$ expansion of $\mu_{\alpha}({\bf q})$ is
$\mathcal{O}(q^{3})$. Finally, the generalised elastic coefficient $\nu({\bf
q})$, whose leading order contribution comes from the homogeneous constant
$\nu$, is given by

$\displaystyle\nu({\bf q})=$ $\displaystyle\frac{k_{B}T}{V}\int d^{3}r_{1}\int
d^{3}r_{2}n({\bf r}_{1})n({\bf r}_{2})e^{-i{\bf q}.({\bf r}_{1}-{\bf r}_{2})}$
$\displaystyle\left[\frac{\delta({\bf r}_{1}-{\bf r}_{2})}{n({\bf
r}_{1})}-c({\bf r}_{1},{\bf r}_{2})\right]$ (94a) $\displaystyle\approx$
$\displaystyle\nu+\mathcal{O}(q^{2})$ (94b)

It can be shown [10] that, as a consequence of the ${\bf
r}_{1}\leftrightarrow{\bf r}_{2}$ symmetry, $\nu({\bf q})$ is real and has
contributions from even powers in a long wavelength expansion in ${\bf q}$.

###  C.2 The coefficients coupling to energy

In this section we derive the small wave vector limit for the coefficients
$\tau_{\alpha},\theta$ and $L$. The initial definition of $\tau_{\alpha}$ (Eq.
(40))

$\displaystyle\tau_{\alpha}({\bf q})=i\sum_{{\bf g}}U^{*}_{{\bf g}}({\bf
q})n_{{\bf q}}g_{\alpha}$ (95)

uses the abbreviated notation of $U_{{\bf g}}$ defined in terms of the inverse
density correlation function $J_{{\bf g}{\bf g}^{\prime}}$ (see Eq. (23)) and
the correlation between fluctuation in energy and Bragg peak amplitudes
$K_{{\bf g}}$ (see Eq. (24)). Therefore, plugging in these definitions,
$\tau_{\alpha}$ can be written in terms of the direct correlation function
$c({\bf r}_{1},{\bf r}_{2})$ and symmetry of these function can be exploited
to derive an expression for $\tau_{\alpha}$ in the small ${\bf q}$ limit.

$\displaystyle\tau_{\alpha}({\bf q})$ $\displaystyle=(\beta V)^{-1}\sum_{{\bf
g},{\bf g}^{\prime}}\int\mathrm{d}^{3}r_{1}\int\mathrm{d}^{3}r_{2}in_{{\bf
g}}g_{\alpha}e^{i{\bf g}\cdot{\bf r}_{1}}K^{*}_{{\bf g}^{\prime}}({\bf
q})e^{-i{\bf g}^{\prime}\cdot{\bf r}_{2}}e^{-i{\bf q}\cdot({\bf r}_{1}-{\bf
r}_{2})}\left[\dfrac{\delta({\bf r}_{1}-{\bf r}_{2})}{n({\bf r}_{1})}-c({\bf
r}_{1},{\bf r}_{2})\right]$ (96a) $\displaystyle=(\beta V)^{-1}\sum_{{\bf
g}^{\prime}}\int\mathrm{d}^{3}r_{2}\nabla_{\alpha}n({\bf r}_{1})K^{*}_{{\bf
g}^{\prime}}({\bf q})e^{-i{\bf g}^{\prime}\cdot{\bf r}_{2}}e^{-i{\bf
q}\cdot({\bf r}_{1}-{\bf r}_{2})}\left[\dfrac{\delta({\bf r}_{1}-{\bf
r}_{2})}{n({\bf r}_{1})}-c({\bf r}_{1},{\bf r}_{2})\right]$ (96b)
$\displaystyle=(\beta V)^{-1}\sum_{{\bf
g}^{\prime}}\int\mathrm{d}^{3}r_{2}K^{*}_{{\bf g}^{\prime}}({\bf q})e^{-i{\bf
g}^{\prime}{\bf r}_{2}}\int\mathrm{d}^{3}r_{1}\nabla_{\alpha}n({\bf
r}_{1})c({\bf r}_{1},{\bf r}_{2})\left(1-e^{-i{\bf q}\cdot({\bf r}_{1}-{\bf
r}_{2})}\right)$ (96c) $\displaystyle=V^{-2}\sum_{{\bf
g}^{\prime}}\int\mathrm{d}^{3}r_{2}\langle\sum_{jk}E_{j}e^{-i{\bf q}\cdot({\bf
r}_{j}-{\bf r}_{k})}e^{i{\bf g}^{\prime}\cdot{\bf r}_{k}}\rangle e^{-i{\bf
g}^{\prime}{\bf r}_{2}}\int\mathrm{d}^{3}r_{1}\nabla_{\alpha}n({\bf
r}_{1})c({\bf r}_{1},{\bf r}_{2})[1-i{\bf q}.{\bf r}_{12}+\dots]$ (96d)
$\displaystyle=iq_{\beta}\tau_{{\alpha\beta}}+\mathcal{O}(q^{2})$ (96e)

where the second rank tensor $\tau_{\alpha\beta}$ is

$\displaystyle\tau_{{\alpha\beta}}$ $\displaystyle=(\beta V)^{-1}\sum_{{\bf
g}^{\prime}}\int\mathrm{d}^{3}r_{2}\bar{K}^{*}_{{\bf g}^{\prime}}e^{-i{\bf
g}^{\prime}{\bf r}_{2}}\int\mathrm{d}^{3}r_{1}\nabla_{\alpha}n({\bf
r}_{1})c({\bf r}_{1},{\bf r}_{2})r_{12,\beta}$ (97)

Note the analogy between the definitions of $\mu_{\alpha}$ and
$\tau_{\alpha}$. From this analogy, symmetry arguments applicable for
$\mu_{\alpha\beta}$ (Eq. (88)) (also see reference [10]), holds for
$\tau_{{\alpha\beta}}$ as well. The symmetry $c({\bf r}_{1},{\bf
r}_{2})=c({\bf r}_{2},{\bf r}_{1})$ and the LMBW equation indicates
$\tau_{{\alpha\beta}}=\tau_{\beta\alpha}$. The term $K^{*}_{{\bf
g}^{\prime}}({\bf q})$, in Eq. (96), is a ${\bf q}$ dependent correlation
function (see Eq. (24)a) between the energy density and Bragg diffraction
amplitudes. In the ${\bf q}\rightarrow 0$ limit the leading order term, in the
expansion of the ${\bf q}$ dependent exponential in its expression, is a
constant $\bar{K}^{*}_{{\bf g}^{\prime}}$ and ${\bf q}$ independent
macroscopic property of the system. Similar to the correlations in Eq. (82),
this quantity is expected to have the periodicity of the lattice structure.

Next we take up the coefficient $\theta$. While deriving the explicit
expressions for $\theta$ in terms of $c({\bf r}_{1},{\bf r}_{2})$ we draw
attention to the analogy between the definitions of $\theta$ and $\nu$.

$\displaystyle\theta({\bf q})$ $\displaystyle=\sum_{{\bf g}}U^{*}_{{\bf
g}}({\bf q})n_{{\bf g}}=\sum_{{\bf g},{\bf g}^{\prime}}K^{*}_{{\bf
g}^{\prime}}({\bf q})J^{*}_{{\bf g}^{\prime}{\bf g}}({\bf q})n_{{\bf g}}$
(98a) $\displaystyle=(\beta V)^{-1}\sum_{{\bf g},{\bf
g}^{\prime}}\int\int\mathrm{d}^{3}r_{1}\mathrm{d}^{3}r_{2}n_{{\bf g}}e^{i{\bf
g}\cdot{\bf r}_{1}}\bar{K}^{*}_{{\bf g}^{\prime}}e^{-i{\bf
g}^{\prime}\cdot{\bf r}_{2}}e^{-i{\bf q}\cdot({\bf r}_{1}-{\bf
r}_{2})}\left[\dfrac{\delta({\bf r}_{1}-{\bf r}_{2})}{n({\bf r}_{1})}-c({\bf
r}_{1},{\bf r}_{2})\right]$ (98b) $\displaystyle=(\beta V)^{-1}\sum_{{\bf
g}^{\prime}}\int\int\mathrm{d}^{3}r_{1}\mathrm{d}^{3}r_{2}n({\bf
r}_{1})\bar{K}^{*}_{{\bf g}^{\prime}}e^{-i{\bf g}^{\prime}\cdot{\bf
r}_{2}}e^{-i{\bf q}\cdot({\bf r}_{1}-{\bf r}_{2})}\left[\dfrac{\delta({\bf
r}_{1}-{\bf r}_{2})}{n({\bf r}_{1})}-c({\bf r}_{1},{\bf r}_{2})\right],\ \
\because n({\bf r}_{1})=\sum_{{\bf g}}n_{{\bf g}}e^{i{\bf g}\cdot{\bf r}_{1}}$
(98c) $\displaystyle=(\beta V)^{-1}\sum_{{\bf
g}^{\prime}}\int\mathrm{d}^{3}r_{2}\bar{K}^{*}_{{\bf g}^{\prime}}e^{-i{\bf
g}^{\prime}\cdot{\bf r}_{2}}-(\beta V)^{-1}\sum_{{\bf
g}^{\prime}}\int\int\mathrm{d}^{3}r_{1}\mathrm{d}^{3}r_{2}n({\bf
r}_{1})\bar{K}^{*}_{{\bf g}^{\prime}}e^{-i{\bf g}^{\prime}\cdot{\bf
r}_{2}}c({\bf r}_{1},{\bf r}_{2})[1-i{\bf q}.{\bf r}_{12}+\dots]$ (98d)

The above equation shows how in the long wavelength limit,

$\displaystyle\theta$ $\displaystyle=(\beta V)^{-1}\sum_{{\bf
g}^{\prime}}\int\mathrm{d}^{3}r_{2}\bar{K}_{{\bf g}^{\prime}}e^{-i{\bf
g}^{\prime}\cdot{\bf r}_{2}}-(\beta V)^{-1}\sum_{{\bf
g}^{\prime}}\int\int\mathrm{d}^{3}r_{1}\mathrm{d}^{3}r_{2}n({\bf
r}_{1})\bar{K}^{*}_{{\bf g}^{\prime}}e^{-i{\bf g}^{\prime}\cdot{\bf
r}_{2}}c({\bf r}_{1},{\bf r}_{2})$ (99)

is a ${\bf q}$ independent real constant similar to $\nu$. Here, similar to
$\tau_{{\alpha\beta}}$, we have used $\bar{K}^{*}_{{\bf g}^{\prime}}$ as the
long wavelength expectation value for $K^{*}_{{\bf g}^{\prime}}({\bf q})$.

Utilising the microscopic expressions for the energy fluctuations (see Eq.
(15) and Eq. (16)) in the definition of $L$ (Eq. (24)c) and taking Taylor
expansion of the exponential functions of ${\bf q}$ leads to

$\displaystyle L({\bf q})$
$\displaystyle=\beta\langle\delta\hat{e}^{*}\delta\hat{e}\rangle-\sum_{{\bf
g},{\bf g}^{\prime}}K_{{\bf g}}({\bf q})J_{{\bf g}{\bf g}^{\prime}}({\bf
q})K^{*}_{{\bf g}^{\prime}}({\bf q})$ (100a)
$\displaystyle=\beta\langle\delta\hat{e}^{*}\delta\hat{e}\rangle-(\beta
V)^{-1}\sum_{{\bf g},{\bf
g}^{\prime}}\int\int\mathrm{d}^{3}r_{1}\mathrm{d}^{3}r_{2}\bar{K}_{{\bf
g}}e^{i{\bf g}\cdot{\bf r}_{1}}\bar{K}^{*}_{{\bf g}^{\prime}}e^{-i{\bf
g}^{\prime}\cdot{\bf r}_{2}}\left[\dfrac{\delta({\bf r}_{1}-{\bf
r}_{2})}{n({\bf r}_{1})}-c({\bf r}_{1},{\bf r}_{2})\right]e^{-i{\bf
q}\cdot({\bf r}_{1}-{\bf r}_{2})}$ (100b)
$\displaystyle=\beta\langle\delta\hat{e}^{*}\delta\hat{e}\rangle-(\beta
V)^{-1}\sum_{{\bf g},{\bf
g}^{\prime}}\int\int\mathrm{d}^{3}r_{1}\mathrm{d}^{3}r_{2}\bar{K}_{{\bf
g}}e^{i{\bf g}\cdot{\bf r}_{1}}\bar{K}^{*}_{{\bf g}^{\prime}}e^{-i{\bf
g}^{\prime}\cdot{\bf r}_{2}}\left[\dfrac{\delta({\bf r}_{1}-{\bf
r}_{2})}{n({\bf r}_{1})}-c({\bf r}_{1},{\bf r}_{2})\right][1-i{\bf q}.{\bf
r}_{12}+\dots]$ (100c)

with the small wavelength limit of $L({\bf q})$, given by a constant

$\displaystyle L$
$\displaystyle=\beta\langle\delta\hat{e}^{*}\delta\hat{e}\rangle-(\beta
V)^{-1}\sum_{{\bf g},{\bf
g}^{\prime}}\int\int\mathrm{d}^{3}r_{1}\mathrm{d}^{3}r_{2}\bar{K}_{{\bf
g}}e^{i{\bf g}\cdot{\bf r}_{1}}\bar{K}^{*}_{{\bf g}^{\prime}}e^{-i{\bf
g}^{\prime}\cdot{\bf r}_{2}}\left[\dfrac{\delta({\bf r}_{1}-{\bf
r}_{2})}{n({\bf r}_{1})}-c({\bf r}_{1},{\bf r}_{2})\right]$ (101)

Here, once again, in the ${\bf q}\rightarrow 0$ limit, the leading order
contributions from the terms $K_{{\bf g}}({\bf q})$ are ${\bf q}$ independent
constants like $\bar{K}_{{\bf g}}$. The arguments related to the term $J_{{\bf
g}{\bf g}^{\prime}}({\bf q})$, in the small ${\bf q}$ limit, are identical to
the ones given for the calculation of $\nu({\bf q})$ (see Eq. (94) and
reference [10, 14, 26]).

##  Appendix D coarse-grained static susceptibility

This appendix aims to derive relations between (i) the intensive thermodynamic
conjugate fields defined within the Mori-Zwanzig projection formalism and
appearing on the left hand side of the Eq. (36), and (ii) the thermodynamic
fields introduced in the expansion of the free energy in Eq. (IV.1) and now
used to represent the partition function in Eq. (107) (also appearing on the
left hand side of Eq. (112)).

The relation between the coarse-grained thermodynamic variables in Eq. (36) is
given by the static susceptibility matrix. The block diagonal structure of the
coarse-grained susceptibility matrix $\bm{\chi}$, similar to its microscopic
higher dimensional analogue in Eq. (20), decouples a $3\times 3$ diagonal
block of correlations between the linear momentum densities, $\bm{\chi^{jj}}$,
from the $5\times 5$ matrix block $\bm{\chi_{nue}}$, allowing one to
separately consider

$\displaystyle V^{-1}\bm{\chi_{nue}}\begin{bmatrix}\delta a\\\ \delta\bm{y}\\\
\delta b\end{bmatrix}=\begin{bmatrix}\delta n\\\ \delta\bm{u}\\\ \delta
e\end{bmatrix}.$ (102)

decoupled from the correlations to linear momentum. The matrix
$\bm{\chi_{nue}}$ represents the static correlations between fluctuations in
number density, displacement fields and energy density (see Eq. (106)). In
this appendix, we focus on interpreting the coarse-grained thermodynamic
conjugate variables after defining the static correlation matrix
$\bm{\chi_{nue}}$

$\displaystyle\lim_{{\bf q}\rightarrow 0}{\bm{\chi}_{nue}}({\bf q})$
$\displaystyle=\lim_{{\bf q}\rightarrow
0}\beta\left(\begin{array}[]{ccc}\langle\delta n^{*}({\bf q})\delta n({\bf
q})\rangle&\langle\delta n^{*}({\bf q})\delta u_{\gamma}({\bf
q})\rangle&\langle\delta n^{*}({\bf q})\delta e({\bf q})\rangle\\\
\langle\delta u_{\alpha}^{*}({\bf q})\delta n({\bf q})\rangle&\langle\delta
u_{\alpha}^{*}({\bf q})\delta u_{\gamma}({\bf q})\rangle&\langle\delta
u_{\alpha}^{*}({\bf q})\delta e({\bf q})\rangle\\\ \langle\delta e^{*}({\bf
q})\delta n({\bf q})\rangle&\langle\delta e^{*}({\bf q})\delta u_{\gamma}({\bf
q})\rangle&\langle\delta e^{*}({\bf q})\delta e({\bf
q})\rangle\end{array}\right)$ (106)

in terms of thermodynamic derivatives (see Eq. (113)). These relations follow
from evaluating the static correlations between the thermodynamic density
fields (see Eq. (108)) in the generalized grand-canonical ensemble [23, 6,
11]. But for that, first we will have to define the partition function (see
Eq. (107)). In case of three dimensional systems, ${\bf u}$ is a three
dimensional vector with components $u_{\alpha}$ corresponding to the three
Cartesian coordinates. For ease of representation, in Eq. (106), we choose to
show the correlations and the thermodynamic derivatives corresponding to one
of the components of ${\bf u}$. These expressions are representative of more
general susceptibility matrices with dimensions appropriate for the systems
concerned. Reference [11] explicitly derives these relations for an isothermal
crystalline solid without considering heat transport associated with the
fluctuations in energy density $\delta e$. For the isothermal case the coarse-
grained susceptibility matrix, representing the correlations between the
density $\delta n$ and displacement $\delta u_{\alpha}$ fluctuations, is a
matrix of dimensions $4\times 4$. In this paper, with the additional energy
fluctuations $\delta e$, the size of the static correlation matrix
$\bm{\chi_{nue}}$ increases to $5\times 5$ to account for the additional
thermodynamic correlations. Apart from this increase in the number of
correlations involved, the main thermodynamic arguments remain identical.

Now we examine how the correlations in Eq. (106) can be represented as
thermodynamic derivatives obtained starting from the partition function.
Drawing analogy to the statistical mechanics of ordered ferromagnetic states
in a system with the Heisenberg Hamiltonian [23], we introduce the elastic
energy ${\bf h}\cdot{\bf u}$ in the Hamiltonian $\mathcal{H}$. Here we follow
the definitions of intensive stress fields $h_{{\alpha\beta}}$ introduced as
the thermodynamic conjugates to the symmetric linear strain fields
$u_{{\alpha\beta}}$ in the free energy expansion in Eq. (IV.1). The partition
function $\mathcal{Z}$ corresponds to the equilibrium crystalline solid in the
limit of $h_{\alpha\beta}\rightarrow 0$.

$\mathcal{Z}=\int\mathrm{d}\Gamma e^{-\beta\mathcal{H}+\beta\mu N-\beta V{\bf
h}\cdot{\bf u}}$ (107)

Within the linear response picture [1], the vector ${\bf u}$ has six
components representing the Voigt symmetric strain fields $u_{{\alpha\beta}}$
: the independent components being for ${\alpha\beta}=11,22,33,23\text{ or
}32,13\text{ or }31,12\text{ or }21$. From a thermodynamics perspective, the
components of the conjugate stress fields $h_{\alpha\beta}$ can be defined as
in Eq. (55)c where the coefficients of the stiffness tensor
$C^{n}_{{\alpha\beta}\gamma\delta}$ are the constants connecting the conjugate
pairs. For a thermodynamic variable density $w$ measured at ${\bf
q}\rightarrow 0$, the following relations for various density correlations can
be derived from Eq. (107) [23]

$\displaystyle\frac{\partial w}{\partial(\beta\mu)}\big{|}_{\beta,\beta{\bf
h}}=\langle\delta w^{*}\delta n\rangle$ (108a) $\displaystyle-$
$\displaystyle\frac{\partial w}{\partial(\beta
h_{{\alpha\beta}})}\big{|}_{\beta,\beta\mu,\beta
h_{\gamma\delta}}=\langle\delta w^{*}\delta u_{{\alpha\beta}}\rangle$ (108b)
$\displaystyle-$ $\displaystyle\frac{\partial
w}{\partial(\beta)}\big{|}_{\beta\mu,\beta{\bf h}}=\langle\delta w^{*}\delta
e\rangle.$ (108c)

For our system of interest, the thermodynamic variable density $w$ denotes
number density $n$, energy density $e$ and linear symmetric strain fields
$u_{{\alpha\beta}}$, which can be equivalently represented in terms of the
displacement fields as given in Eq. (45). Consider representing fluctuations
in these three quantities in terms of the partial derivatives of the three
intensive fields $\beta\mu,\beta h_{{\alpha\beta}}$ and $\beta$ now introduced
in Eq. (108) through the definition of the partition function in Eq. (107).
With $\sum_{i=1}^{n}\dfrac{\partial f}{\partial x_{i}}\delta x_{i}=\delta f$
being a general form for the total derivative of a function $f$ of variables
$x_{1},x_{2},\dots,x_{n}$, the quantities $\delta n,\delta u_{\gamma\delta}$
and $\delta e$ can be written as follows

$\displaystyle\frac{\partial n}{\partial(\beta\mu)}\big{|}_{\beta,\beta{\bf