Phonon hydrodynamics, thermal conductivity and second sound in 2D   crystals

Pablo Scuracchio; Karl H. Michel; Francois M. Peeters

arXiv:1904.06327·cond-mat.stat-mech·April 15, 2019

Phonon hydrodynamics, thermal conductivity and second sound in 2D crystals

Pablo Scuracchio, Karl H. Michel, Francois M. Peeters

PDF

TL;DR

This paper derives hydrodynamic equations for 2D crystals to analyze phonon-related phenomena like second sound and thermal conductivity, providing temperature-dependent transport coefficients without relaxation time assumptions.

Contribution

It introduces a novel derivation of hydrodynamic equations for 2D crystals, linking thermal and mechanical phenomena and calculating transport coefficients directly from correlation functions.

Findings

01

Identification of temperature-dependent transport coefficients

02

Observation of coupled thermal and mechanical resonances

03

Validation of results for various 2D hexagonal crystals

Abstract

Starting from our previous work where we have obtained a system of coupled integro-differential equations for acoustic sound waves and phonon density fluctuations in 2D crystals, we derive here the corresponding hydrodynamic equations and study their consequences as function of temperature and frequency. These phenomena encompass propagation and damping of acoustic sound waves, diffusive heat conduction, second sound and Poiseuille heat flow, all of which are characterized by specific transport coefficients. We calculate these coefficients by means of correlation functions without using the concept of relaxation time. Numerical calculations are performed as well in order to show the temperature dependence of the transport coefficients and of the thermal conductivity. As a consequence of thermal tension mechanical and thermal phenomena are coupled. We calculate the dynamic…

Tables2

Table 1. Table 1: Parameters used for the calculation of the structure functions S s s subscript 𝑆 𝑠 𝑠 S_{ss} and S θ θ subscript 𝑆 𝜃 𝜃 S_{\theta\theta} in the diffusive regime.

$m$	$q_{1}$ (m^-1)	$c_{L}$ (ms^-1)	$β^{', 2}$	$q_{1} {\hat{η}}_{11}$	$q_{1} λ$	$T$ (K)
$2 m_{C}$	$4 \times 10^{5}$	$23 \times 10^{3}$	$c_{L}^{2} / 10$	$c_{L} / 10$	$c_{L} / 5$	$100$

Table 2. Table 2: Parameters used for the calculation of the structure functions S s s subscript 𝑆 𝑠 𝑠 S_{ss} and S θ θ subscript 𝑆 𝜃 𝜃 S_{\theta\theta} in the second sound regime.

$m$	$q_{1}$ (m^-1)	$c_{L}$ (ms^-1)	$β^{', 2}$	$q_{1} {\hat{η}}_{11}$	$q_{1} \tilde{λ}$	$V$	$T$ (K)
$2 m_{C}$	$4 \times 10^{5}$	$23 \times 10^{3}$	$c_{L}^{2} / 10$	$c_{L} / 10$	$c_{L} / 7$	$c_{L} / 2$	$100$

Equations279

\displaystyle H_{h}=\frac{1}{2}\sum_{\vec{q},\alpha}\Big{[}{P^{\alpha}_{\vec{q}}}^{\dagger}P^{\alpha}_{\vec{q}}+\omega^{2}(\vec{q},\alpha){Q^{\alpha}_{\vec{q}}}^{\dagger}Q^{\alpha}_{\vec{q}}\Big{]}\;,

\displaystyle H_{h}=\frac{1}{2}\sum_{\vec{q},\alpha}\Big{[}{P^{\alpha}_{\vec{q}}}^{\dagger}P^{\alpha}_{\vec{q}}+\omega^{2}(\vec{q},\alpha){Q^{\alpha}_{\vec{q}}}^{\dagger}Q^{\alpha}_{\vec{q}}\Big{]}\;,

\displaystyle\Phi^{(3)}=\frac{1}{2}\sum_{\vec{p},\vec{k},\vec{q},i}\Phi^{(3)}\left(\!\begin{array}[]{ccc}i&\zeta&\zeta\\ \vec{q}&\vec{k}&\vec{p}\end{array}\!\right)Q^{i}_{\vec{q}}Q^{\zeta}_{\vec{k}}Q^{\zeta}_{\vec{p}}\;.

s_{i} (q) = \frac{1}{m} Q_{q}^{α} δ_{α i},

s_{i} (q) = \frac{1}{m} Q_{q}^{α} δ_{α i},

p_{i} (q) = m P_{q}^{α} δ_{α i} .

p_{i} (q) = m P_{q}^{α} δ_{α i} .

\displaystyle\Phi^{(3)}\left(\!\begin{array}[]{ccc}i&\zeta&\zeta\\ \vec{q}&\vec{k}&\vec{p}\end{array}\!\right)=\frac{i8}{\sqrt{Nm^{3}}}\sum_{s}\phi^{(3)}_{izz}(A_{1};B_{s})

\displaystyle\Phi^{(3)}\left(\!\begin{array}[]{ccc}i&\zeta&\zeta\\ \vec{q}&\vec{k}&\vec{p}\end{array}\!\right)=\frac{i8}{\sqrt{Nm^{3}}}\sum_{s}\phi^{(3)}_{izz}(A_{1};B_{s})

\displaystyle\times\text{cos}\Big{[}\frac{(\vec{q}+\vec{k}+\vec{p})\cdot\vec{\rho}(B_{s})}{2}\Big{]}\text{sin}\Big{[}\frac{\vec{q}\cdot\vec{r}(B_{s})}{2}\Big{]}\text{sin}\Big{[}\frac{\vec{k}\cdot\vec{r}(B_{s})}{2}\Big{]}

\displaystyle\times\text{sin}\Big{[}\frac{\vec{p}\cdot\vec{r}(B_{s})}{2}\Big{]}\Delta(\vec{q}+\vec{k}+\vec{p})\;.

Q_{q}^{α}

Q_{q}^{α}

P_{q}^{α}

\displaystyle\Big{[}m\omega^{2}\delta_{ij}-v_{2D}q_{k}q_{l}C_{ik,jl}\Big{]}s_{j}(\vec{q},\omega)=-F_{i}(\vec{q},\omega)

\displaystyle\Big{[}m\omega^{2}\delta_{ij}-v_{2D}q_{k}q_{l}C_{ik,jl}\Big{]}s_{j}(\vec{q},\omega)=-F_{i}(\vec{q},\omega)

+ i q_{j} \frac{ℏ}{N} k, α \sum δ_{α ζ} h_{j i} (k, ζ) n^{'} (k, α) ν (k, α; q, ω) .

\displaystyle\Big{[}\omega-q_{i}v_{i}(\vec{k},\alpha)\Big{]}\nu(\vec{k},\alpha;\vec{q},\omega)+i\hbar\omega q_{i}\delta_{\alpha\zeta}h_{ij}(\vec{k},\zeta)s_{j}(\vec{q},\omega)

\displaystyle\Big{[}\omega-q_{i}v_{i}(\vec{k},\alpha)\Big{]}\nu(\vec{k},\alpha;\vec{q},\omega)+i\hbar\omega q_{i}\delta_{\alpha\zeta}h_{ij}(\vec{k},\zeta)s_{j}(\vec{q},\omega)

+ q_{i} v_{i} (k, α) ℏ ω (k, α) \frac{Θ ( q , ω )}{T} = - i \hat{C} ν (k, α; q, ω),

h_{ij} (k, ζ)

h_{ij} (k, ζ)

\times (k \cdot r (B_{s}))^{2} r_{j} (B_{s}) .

n^{'} (k, α) = \frac{n ( k , α ) ( 1 + n ( k , α ))}{k _{B} T} = - \frac{1}{ℏ} \frac{\partial n ( k , α )}{\partial ω ( k , α )},

n^{'} (k, α) = \frac{n ( k , α ) ( 1 + n ( k , α ))}{k _{B} T} = - \frac{1}{ℏ} \frac{\partial n ( k , α )}{\partial ω ( k , α )},

n (k, α; q, ω) = \tilde{n} (k, α; q, ω) + n^{'} (k, α) ν (k, α; q, ω),

n (k, α; q, ω) = \tilde{n} (k, α; q, ω) + n^{'} (k, α) ν (k, α; q, ω),

\tilde{n} (k, α; q, ω) = (exp (ℏ ω (k, α; q, ω) / k_{B} T) - 1)^{- 1},

\tilde{n} (k, α; q, ω) = (exp (ℏ ω (k, α; q, ω) / k_{B} T) - 1)^{- 1},

ω (k, α; q, ω) = ω (k, α) - i δ_{α ζ} q_{i} \frac{h _{ij} ( k , ζ )}{N} s_{j} (q, ω) .

ω (k, α; q, ω) = ω (k, α) - i δ_{α ζ} q_{i} \frac{h _{ij} ( k , ζ )}{N} s_{j} (q, ω) .

\hat{C} ν (k, α; q, ω)

\hat{C} ν (k, α; q, ω)

+ (1 - δ_{α ζ}) C^{(2)} ν (k, α; q, ω),

C^{(1)} ν (k, ζ) = \frac{1}{n ^{'} ( k , ζ )} \times

C^{(1)} ν (k, ζ) = \frac{1}{n ^{'} ( k , ζ )} \times

\displaystyle\sum_{\vec{p}\vec{h},i}\bigg{\{}W\left(\!\begin{array}[]{ccc}\zeta&\zeta&i\\ \vec{h}&\vec{k}&\vec{p}\end{array}\!\right)\Big{[}\nu(\vec{k},\zeta)-\nu(\vec{h},\zeta)+\nu(\vec{p},i)\Big{]}

\displaystyle+W\left(\!\begin{array}[]{ccc}i&\zeta&\zeta\\ \vec{p}&\vec{h}&\vec{k}\end{array}\!\right)\Big{[}\nu(\vec{k},\zeta)-\nu(\vec{p},i)+\nu(\vec{h},\zeta)\Big{]}

\displaystyle+W\left(\!\begin{array}[]{ccc}\zeta&i&\zeta\\ \vec{k}&\vec{p}&\vec{h}\end{array}\!\right)\Big{[}\nu(\vec{k},\zeta)-\nu(\vec{p},i)-\nu(\vec{h},\zeta)\Big{]}\bigg{\}}\;,

C^{(2)} ν (k, j) = \frac{1}{n ^{'} ( k , j )} \times

C^{(2)} ν (k, j) = \frac{1}{n ^{'} ( k , j )} \times

\displaystyle\sum_{\vec{p}\vec{h}}\bigg{\{}\frac{1}{2}W\left(\!\begin{array}[]{ccc}j&\zeta&\zeta\\ \vec{k}&\vec{p}&\vec{h}\end{array}\!\right)\big{[}\nu(\vec{k},j)-\nu(\vec{p},\zeta)-\nu(\vec{h},\zeta)\big{]}

\displaystyle+W\left(\!\begin{array}[]{ccc}\zeta&j&\zeta\\ \vec{p}&\vec{k}&\vec{h}\end{array}\!\right)\Big{[}\nu(\vec{k},j)-\nu(\vec{p},\zeta)+\nu(\vec{h},\zeta)\Big{]}\bigg{\}}\;.

\displaystyle W\left(\!\begin{array}[]{ccc}\zeta&j&\zeta\\ \vec{k}&\vec{p}&\vec{h}\end{array}\!\right)

\displaystyle W\left(\!\begin{array}[]{ccc}\zeta&j&\zeta\\ \vec{k}&\vec{p}&\vec{h}\end{array}\!\right)

\times k_{B} T n^{'} (k, ζ) n^{'} (p, j) n^{'} (h, ζ)

\times δ (ω (h, ζ) - ω (k, ζ) + ω (p, j)),

\displaystyle W\left(\!\begin{array}[]{ccc}j&\zeta&\zeta\\ \vec{p}&\vec{h}&\vec{k}\end{array}\!\right)

\times k_{B} T n^{'} (p, j) n^{'} (h, ζ) n^{'} (k, ζ)

\times δ (ω (h, ζ) - ω (p, j) + ω (k, ζ)),

\displaystyle\Psi^{(3)}\left(\!\begin{array}[]{ccc}j&\zeta&\zeta\\ -\vec{p}&\vec{h}&\vec{k}\end{array}\!\right)=\frac{\Phi^{(3)}\left(\!\begin{array}[]{ccc}j&\zeta&\zeta\\ -\vec{p}&\vec{h}&\vec{k}\end{array}\!\right)}{\sqrt{8\omega(\vec{k},\zeta)\omega(\vec{p},j)\omega(\vec{h},\zeta)}}\;.

\displaystyle\Psi^{(3)}\left(\!\begin{array}[]{ccc}j&\zeta&\zeta\\ -\vec{p}&\vec{h}&\vec{k}\end{array}\!\right)=\frac{\Phi^{(3)}\left(\!\begin{array}[]{ccc}j&\zeta&\zeta\\ -\vec{p}&\vec{h}&\vec{k}\end{array}\!\right)}{\sqrt{8\omega(\vec{k},\zeta)\omega(\vec{p},j)\omega(\vec{h},\zeta)}}\;.

⟨ a_{1} ∣ a_{2} ⟩ = \frac{1}{N} p, α \sum a_{1}^{*} (p, α) n^{'} (p, α) a_{2} (p, α),

⟨ a_{1} ∣ a_{2} ⟩ = \frac{1}{N} p, α \sum a_{1}^{*} (p, α) n^{'} (p, α) a_{2} (p, α),

⟨ ν_{1} ∣ \hat{C} ∣ ν_{2} ⟩ = ⟨ ν_{2} ∣ \hat{C} ∣ ν_{1} ⟩,

⟨ ν_{1} ∣ \hat{C} ∣ ν_{2} ⟩ = ⟨ ν_{2} ∣ \hat{C} ∣ ν_{1} ⟩,

\displaystyle\langle\nu_{1}|\hat{C}|\nu_{2}\rangle=\frac{1}{N}\sum_{\vec{p}}\Big{\{}\nu_{1}^{*}(\vec{p},\zeta)C^{(1)}\nu_{2}(\vec{p},\zeta)

\displaystyle\langle\nu_{1}|\hat{C}|\nu_{2}\rangle=\frac{1}{N}\sum_{\vec{p}}\Big{\{}\nu_{1}^{*}(\vec{p},\zeta)C^{(1)}\nu_{2}(\vec{p},\zeta)

\displaystyle+\sum_{j}\nu_{1}(\vec{p},j)^{*}C^{(2)}\nu_{2}(\vec{p},j)\Big{\}}

\displaystyle=\frac{1}{N}\sum_{\vec{k}\vec{p}\vec{h},j}\bigg{\{}W\left(\!\begin{array}[]{ccc}\zeta&j&\zeta\\ \vec{k}&\vec{p}&\vec{h}\end{array}\!\right)\big{[}\nu_{1}^{*}(\vec{k},\zeta)-\nu_{1}^{*}(\vec{h},\zeta)-\nu_{1}^{*}(\vec{p},j)\big{]}

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Phonon hydrodynamics, thermal conductivity and second sound in 2D crystals

P. Scuracchio

[email protected]

K. H. Michel

[email protected]

F. M. Peeters

[email protected]

Universiteit Antwerpen, Department of Physics, Groenenborgerlaan 171, Be-2020 Antwerpen, Belgium

Abstract

Starting from our previous work where we have obtained a system of coupled integro-differential equations for acoustic sound waves and phonon density fluctuations in 2D crystals, we derive here the corresponding hydrodynamic equations and study their consequences as function of temperature and frequency. These phenomena encompass propagation and damping of acoustic sound waves, diffusive heat conduction, second sound and Poiseuille heat flow, all of which are characterized by specific transport coefficients. We calculate these coefficients by means of correlation functions without using the concept of relaxation time. Numerical calculations are performed as well in order to show the temperature dependence of the transport coefficients and of the thermal conductivity. As a consequence of thermal tension mechanical and thermal phenomena are coupled. We calculate the dynamic susceptibilities for displacement and temperature fluctuations and study their resonances. Due to the thermo-mechanical coupling the thermal resonances such as Landau-Placzek peak and second sound doublet appear in the displacement susceptibility and conversely the acoustic sound wave doublet appears in the temperature susceptibility, Our analytical results do not only apply to graphene but are also valid for arbitrary 2D crystals with hexagonal symmetry like 2D h-BN, 2H-transition metal dichalcogenides and oxides.

I Introduction

Phonon hydrodynamic phenomena such as propagation and attenuation of sound waveslandau , thermal conductivitypeierlsfirst ; peierls ; casimir ; klemens ; leibfriedscho ; ziman , Poiseuille flowsussmann ; gurzhi ; mezhov ; guyer1 and second soundward1 ; ward2 ; dingle ; sussmann ; prohofsky ; mezhov ; ackerman ; guyer ; GM1 ; GM2 ; guyer1 in insulating 3D crystals have been investigated by experiment and theory since a long time. From a microscopic point of view these phenomena are closely related to anharmonic phonon interactionpeierls ; ziman . Recently the subject has known a revival due to the discovery of 2D crystalsnovo . Layered and 2D crystals exhibit a peculiar anharmonic interaction, linear in in-plane displacements and quadratic in out-of-plane displacements (flexural modes).This leads to unusual static effects such as negative thermal expansionliftshitz ; mounet and is also at the origin of unusual high thermal conductivity in graphenebalandin1 ; seol . In particular, it was foundseol ; lindsay that the ZA (flexural acoustic) modes carry most of the heat in suspended graphene and that the selection rule for anharmonic phonon scattering in 2D crystals restricts momentum dissipating umklapp scattering. Numerical calculations using the single-mode relaxation time approximation confirm that thermal transport is dominated by ZA modesbonini .

A study of phonon lifetimes as function of the wave-vector throughout the 2D Brillouin zone has shown that in a broad temperature range below room temperature the decay rate of flexural modes is much less affected by umklapp processes than the decay rate of in-plane modesseba . The distinctive role played by flexural modes in anharmonic scattering processes has stimulated recent theoretical studies of hydrodynamic phenomena such as second sound and Poiseuille flow in 2D crystalslee_naturecom ; cepe_naturecom ; cepe_prx . In these works the derivation of hydrodynamic phenomena is based on the solution of a linearized Peierls-Boltzmann transport equation (kinetic equation) for the space and time dependent phonon density distribution, either by iterative methodslee_naturecom or by a variational approachcepe_naturecom ; cepe_prx . It was concluded that in 2D crystals hydrodynamic phenomena should exist in a broad temperature range up to and above room temperature. For a most recent review of the literature and discussion on thermal transport in 2D materials, see Ref. (revchinos, ).

While in the aforementioned studies all emphasis has been concentrated on the kinetic equation, the present authors have recently derived a system of dynamic equations for in-plane displacement correlations coupled to phonon density fluctuationsSMP . A linear coupling between in-plane lattice displacements and ZA phonon density fluctuations is a direct consequence of the anharmonic interaction mentioned before. As a result of linear response theory and Green’s functions techniques we find that phonon density fluctuations enter as an additional driving term in the sound wave equation for in-plane lattice deformations while the kinetic equation contains an additional term due to coherent lattice displacements. This calls for a study of the ensuing thermo-mechanical effects on phonon hydrodynamics in 2D crystals. Such an investigation paves also the way for a broader use of experimental techniquesgriffin since thermal effects will appear in mechanical response and vice-versa.

The content of this paper is the following. In Sec. II we recall the basic concepts and results of our previous workSMP . In Sec. III we present solutions for the coupled system of dynamic equations, first by assuming only energy conservation in the collision term of the kinetic equation, and second by taking also into account approximate conservation of crystal momentum. In the former case temperature variations will be described by a diffusion equation coupled to the in-plane lattice deformations, in the latter case depending on the frequency of the perturbation, temperature fluctuations will be described either by a wave equation (second sound) coupled to lattice deformations (first sound) or by Poiseuille flow. In both cases the equation for in-plane sound waves will be coupled to temperature fluctuations. In the steady-state case phonon transport in a confined crystal is described by Poiseuille flow. In Sec. IV we derive thermal and mechanical dynamic response functions and study their resonances in the diffusive and second sound regime. In Sec. V we calculate transport coefficients such as lattice viscosity, kinematic phonon viscosity and thermal diffusion by using a correlation function method developed by Götze and one of the present authorsGM3 . In Sec. VI we present numerical evaluations of these transport coefficients as functions of temperature for the case of graphene. Concluding remarks in Sec. VII and a comment on recent experiments close the paper.

II Basic concepts

We recall some basic concepts and results from our previous workSMP on sound waves and flexural mode dynamics in 2D hexagonal crystals. The previous paper is hereafter referred to as I. The Hamiltonian operator which describes the coupling between acoustic in-plane and out-of-plane displacements was defined as $H=H_{h}+\Phi^{(3)}$ . In terms of normal coordinates $Q^{\alpha}_{\vec{q}}$ and their conjugate momenta $P^{\alpha}_{\vec{q}}$ one has

[TABLE]

The wave vectors $\vec{q}$ , $\vec{k}$ , and $\vec{p}$ all belong to the first 2D Brillouin zone (BZ) and $\alpha=\{1,2,\zeta\}$ refers to the 3 acoustic polarizations LA (longitudinal acoustic in-plane), TA (transverse acoustic in-plane) and ZA (flexural acoustic), respectively. It is understood that $\omega^{2}(\vec{q},\alpha)$ is the renormalized flexural phonon frequency and hence $\omega(\vec{q},\zeta)$ is linear in $\vec{q}$ at long wavelengthsseba . Along this paper, polarizations denoted with Latin indexes $i,j=\{1,2\}$ refer only to in-plane modes.

At long wavelengths, the Fourier-transformed center of mass displacements per unit cell $s_{i}(\vec{q})$ are related to the normal coordinates by

[TABLE]

where $m$ is the total mass per unit cell, i. e. $m=2m_{C}$ , $m_{B}+m_{N}$ or $m_{Mo}+2m_{S}$ in the case of graphene, h-BN or MoS2. The corresponding conjugate momenta are

[TABLE]

The third order anharmonic coupling, calculated for the case of a central force interaction potential between nearest neighbors in a 2D hexagonal crystal of $N$ unit cells is given by

[TABLE]

Here the interaction parameters $\phi^{(3)}_{izz}(A_{1};B_{s})$ and the vectors $\vec{r}(B_{s})$ , $\vec{\rho}(B_{s})$ have been defined in I. The crystal momentum conservation function $\Delta(\vec{q}+\vec{k}+\vec{p})=\sum_{\vec{G}}\>\delta_{\vec{q}+\vec{k}+\vec{p},\vec{G}}$ accounts for normal (N-processes) and umklapp processes (U-processes) when the reciprocal 2D lattice vector $\vec{G}=0$ or $\vec{G}\neq 0$ , respectively.

In terms of phonon creation and annihilation operators $b^{\alpha\dagger}_{\vec{q}}$ and $b^{\alpha}_{\vec{q}}$ , the normal coordinates are written as

[TABLE]

where $b^{\alpha}_{\vec{q}}$ and $b^{\alpha\dagger}_{\vec{q}}$ satisfy the usual commutation rules for Bose operators. Using linear response theory and thermal Green’s functions techniques we have derived in I the following system of coupled dynamic equations. The first one describes sound waves due to coherent in-plane lattice deformations $s_{i}(\vec{q},\omega)$ :

[TABLE]

The second one is a kinetic equation

[TABLE]

where the quantity $\nu(\vec{k},\alpha;\vec{q},\omega)$ describes the deviation from local equilibrium of the dynamic phonon density distribution (see Eq. (15) below). Both equations are coupled by the thermoelastic coupling

[TABLE]

In agreement with Onsager’s reciprocity principle the coupling is symmetric in Eqs. (11) and (12). The structure of the anharmonic coupling $\Phi^{(3)}$ , Eq. (3), implies that only the flexural phonon density deviation $n^{\prime}(\vec{k},\zeta)\nu(\vec{k},\zeta;\vec{q},\omega)$ is coupled to the in-plane deformations. The factor $n^{\prime}(\vec{k},\zeta)$ stands for

[TABLE]

where $k_{B}$ is the Boltzmann constant and $n(\vec{k},\alpha)=[\textrm{exp}(\hbar\omega/k_{B}T)-1]^{-1}$ is the equilibrium Bose-Einstein distribution for acoustic phonons of energy $\hbar\omega(\vec{k},\alpha)$ , at equilibrium temperature $T$ . The remaining quantities occurring in Eqs. (11) and (12) are the Fourier transformed external mechanical force $F_{i}(\vec{q},\omega)$ and temperature perturbation $\Theta(\vec{q},\omega)$ , the collision operator $\hat{C}$ (see Eqs. (18) - (33) below), the isothermal elastic constants $C_{ik,jl}$ , the phonon group velocity $v_{i}(\vec{k},\alpha)=\partial\omega(\vec{k},\alpha)/\partial k_{i}$ and the area of the crystal unit cell $v_{2D}$ . As shown in I, $\nu(\vec{k},\alpha;\vec{q},\omega)$ is related to the non-equilibrium phonon density $n(\vec{k},\alpha;\vec{q},\omega)$ by

[TABLE]

where

[TABLE]

is the local equilibrium phonon distribution with frequency

[TABLE]

As originally pointed out by Akhiezerakhiezer , a sound wave passing through a crystal causes a disturbance of the distribution of thermal phononsbomel ; woodruff ; GM1 . It was shown in I that as a consequence of collisions (quadratic effects in $\Phi^{(3)}$ ) the in-plane phonon density is driven out of equilibrium too. The collision term in the kinetic equation (12) then stands for

[TABLE]

where

[TABLE]

and

[TABLE]

Here and in the following we simplify the actual notation by not mentioning explicitly the $\vec{q},\omega$ dependence of the quantity $\nu(\vec{k},\alpha)$ . We recall that expressions (18)- (33) correspond to linearized collision terms of Peierls-Boltzmann type kinetic equations. The transition probabilities $W$ in (27) and (33) are given by

[TABLE]

with

[TABLE]

While the $\delta$ -functions account for phonon energy conservation in scattering processes involving three phonons, the $\Delta$ -function contained in $\Phi^{(3)}$ ensures that crystal momentum is only conserved up to a reciprocal lattice vector. Since U-processes involve large wave vectors close to the BZ’s boundaries, the corresponding phonon energies $\hbar\omega_{BZ}$ are relatively large (see e.g. acoustic phonon dispersions in graphene)mohr . The phonon occupation factors $n^{\prime}$ contained in the transition probabilities $W$ are only significant for temperatures $k_{B}T\approx\hbar\omega_{BZ}$ . As it is well known (see e.g. Zimanziman ), U-processes are frozen out as $\approx\textrm{e}^{-\Theta_{D}/T}$ , with $\Theta_{D}$ of the order of the Debye temperature. With decreasing $T$ only N-processes survive. We recall that the hydrodynamic regime corresponds to a situation where the system is close to full thermodynamic equilibrium. Such a local equilibrium situation is established by N-processes while the subsequent decay to total equilibrium is described by U-processes or other crystal momentum destroying processes such as impurity and boundary scattering.

While Eq. (11) describes long wavelength and low frequency sound waves, the reduction of the kinetic equation (12) to the hydrodynamic case, i.e. to an equation for local temperature fluctuations requires an additional procedure. The principle is based on the existence of conservation laws and the solution of the kinetic equation in terms of the corresponding local densitiesHuang . We observe that Eq. (12) is an inhomogeneous linear integral equation with an Hermitian kernel. Indeed defining a scalar product for any two functions $a_{1}(\vec{k},\alpha)$ and $a_{2}(\vec{k},\alpha)$ of the phonon variables as

[TABLE]

we obtain for the collision term

[TABLE]

where

[TABLE]

Consequently $\hat{C}$ is an Hermitian operator and admits a spectral representation

[TABLE]

where $\omega^{l}$ and $\chi^{l}$ are the corresponding eigenvalues and orthogonalized eigenfunctions. Since the transition probabilities satisfy $W\geq 0$ , the eigenvalues $\omega^{l}$ (also called relaxation frequencies) satisfy $\omega^{l}\geq 0$ . From the conservation of energy in the phonon scattering processes Eq. (44) it follows that

[TABLE]

is an eigenfunction of $\hat{C}$ with eigenvalue $\omega^{0}=0$ . Strictly speaking there are no further conservation laws for the scattering processes. However at low temperature where U-processes are frozen out, crystal momentum $\hbar\vec{k}$ can be treated as a conserved quantity.

III Hydrodynamic equations

The solution of an inhomogeneous linear integral equation with an Hermitian kernel can be obtained by means of a series expansion in terms of the eigenfunctionscourant . We will start from an expansion of the non equilibrium phonon density $\nu(\vec{k},\alpha;\vec{q},\omega)$ in terms of the eigenfunctions $\chi^{l}(\vec{k},\alpha)$ of the collision operator. We proceed in two steps, first (subsection A) we will assume that only phonon energy is conserved in collisions and second (subsection B) we will in addition treat crystal momentum as an almost conserved quantity.

III.1 Energy Conservation

The non equilibrium phonon density is written as

[TABLE]

Here we have separated off the coefficient $a^{0}(\vec{q},\omega)=\theta(\vec{q},\omega)$ which refers to phonon energy conservation. In the following we will calculate the resolvent of the integral equation (12) by using perturbation theory. Following Ref. (GM2, ) we define the hydrodynamic regime as that one for which the external frequency $\omega$ (and corresponding wave vector $\vec{q}$ ) is small in comparison to the relaxation frequencies $\omega^{l}$ . We consider the situation of high $T$ such that a distinction between N-processes and U-processes is irrelevant. Then $\chi^{l}(\vec{k},\omega)$ with $l>0$ exhausts the spectrum of crystal momentum conserving and destroying scattering processes. Recalling from Sec. II that $\nu$ has the dimension of energy, we see that the quantity $\theta(\vec{q},\omega)$ has the dimension of temperature. In order to perform a reduction from the $\vec{k},\alpha$ dependent phonon density $\nu(\vec{k},\alpha;\vec{q},\omega)$ to an equation for the temperature fluctuations $\theta(\vec{q},\omega)$ , we insert expression (61) into Eqs. (11) and (12). We then isolate the quantities $\theta(\vec{q},\omega)$ and $a^{l}(\vec{q},\omega)$ by multiplying Eq. (12) scalarly with $\chi^{0}(\vec{k},\alpha)$ and $\chi^{l}(\vec{k},\alpha)$ respectively. In evaluating the scalar products we take into account the parity of the relevant functions: while $\omega(\vec{k},\alpha)$ and $h_{ij}(\vec{k},\zeta)$ are even in $\vec{k}\to-\vec{k}$ , $v_{i}(\vec{k},\alpha)$ is uneven.

Defining the specific heat per unit cell as

[TABLE]

and the thermal tensionseba by

[TABLE]

we obtain from Eq. (11)

[TABLE]

and from Eq. (12)

[TABLE]

and

[TABLE]

Taking into account that $0\leq\omega<\omega^{l}$ we use the approximation

[TABLE]

and eliminate $a^{l}(\vec{q},\omega)$ in Eqs. (64) and (65) with the following results:

[TABLE]

and

[TABLE]

Here we have introduced the tensors

[TABLE]

As already mentioned, $\theta(\vec{q},\omega)$ has the dimension of temperature. Taking the static limit ( $\omega\to 0$ ) of Eq. (69), we see that $\theta(\vec{q},\omega=0)$ is equal to the external temperature source $\Theta(\vec{q},\omega=0)$ . We therefore infer that $\theta(\vec{q},\omega)$ is the local temperature in the hydrodynamic regime. Eq. (69) has the meaning of a temperature diffusion equation where the quantity $\lambda_{ij}$ is the thermal diffusion coefficient. On the right hand side of this equation, beside the external temperature source, the lattice deformations $s_{i}(\vec{q},\omega)$ act as an additional perturbation. The coupling between thermal and mechanical quantities, in casu local temperature and coherent lattice deformations respectively, is mediated by the thermal tension coefficient $\beta_{ij}$ . Similarly, in the equation of motion (68) for the in-plane lattice deformations, in addition to the elastic restoring forces there appears now a dissipation term with dynamic viscosity $\eta_{ij,kl}$ , while the local temperature variations, again mediated by the thermal tension, act as an additional driving force.

After inverse Fourier-transforming we rewrite Eqs. (68) and (69) as hydrodynamic equations where $s_{i}(\vec{r},t)$ and $\theta(\vec{r},t)$ are slowly space and time dependent variables describing the deformation of the crystalline lattice and the changes of the temperature respectively:

[TABLE]

and

[TABLE]

Here $\partial_{j}$ stands for $\partial/\partial r_{j}$ and $\partial_{t}=\partial/\partial t$ . From now on, this scenario will be addressed as the diffusion regime (DR).

In absence of an external heat source and of coupling to the lattice the last equation reduces to the Fourier temperature diffusion equation. In the static limit ( $\partial_{t}=0$ ) and in absence of mechanical forces, Eqs. (72) and (73) reduce to

[TABLE]

i.e. the equation of thermal expansion. We observe that so far we did not make use of hexagonal symmetry, hence Eqs. (68) and (69) or equivalently (72) and (73) are valid for 2D crystals of any symmetry.

III.2 Energy and crystal momentum conservation

In subsection III.1 we have shown that if only phonon energy is conserved in scattering processes, temperature changes satisfy a diffusion equation. Here we will assume a situation where U-processes and other crystal momentum destroying processes such as phonon scattering with static impurities are relatively unimportant, while N-processes are still efficient to establish local equilibrium. Treating crystal momentum as an almost conserved quantity we replace Eq. (61) by

[TABLE]

where the crystal momentum

[TABLE]

is treated as an eigenfunction with eigenfrequency $\omega^{i}$ and where $a^{i}(\vec{q},\omega)$ has the dimension of velocity. We assume that $\chi^{i}(\vec{k})$ with $i=1,2$ exhausts the spectrum of momentum destroying processes and that $\chi^{l}(\vec{k},\alpha)$ with $l>2$ refers to N-processes. Such a separation is meaningful at low $T$ where U-processes freeze out. If in addition other momentum destroying processes such as elastic impurity scattering are included in the collision operator, we assume that the concentration of impurities is sufficiently weak such that the combined relaxation frequency $\omega^{i}$ still satisfies $\omega^{i}<\omega^{l}$ , for $i=1,2$ and $l>2$ . The meaning of the first two terms on the right hand side of Eq. (75) can be traced back by linearizing the local displaced or drifting phonon distribution:

[TABLE]

where $\vec{a}(\vec{q},\omega)$ is the drift velocity, and $\omega(\vec{k},\alpha;\vec{q},\omega)$ , Eq. (17), is the acoustic phonon frequency in a deforming lattice. The last term in Eq. (75) accounts for N-processes with relaxation frequencies $\omega^{l}$ . It has been noticed that Eq. (77) corresponds to a distribution function towards which normal processes relaxwoodruff in a deformed lattice. In the steady-state case we recover the displaced phonon distribution with a constant drift velocityklemens ; callaway .

In extending the procedure of subsection III.1 we now calculate $\theta$ , $a^{l}$ and $a^{i}$ . By taking the scalar products we use the fact that $\chi^{i}(\vec{k})$ is an uneven function of $\vec{k}$ . Substituting $\nu(\vec{k},\alpha;\vec{q},\omega)$ , Eq. (75), into Eq. (12), multiplying scalarly by $\chi^{0}(\vec{k},\alpha)$ and eliminating $a^{l}$ by using now $0\leq\omega<\omega^{l}$ with $l\geq 2$ we obtain

[TABLE]

Here the quantities $\lambda^{\prime}_{ik}$ are defined by an expression similar to Eq. (71) except that the $l$ -sum does not include terms with $l=0,1,2$ that refer to conserved quantities. In the last term on the right hand side we have defined

[TABLE]

The quantities $a^{i}(\vec{q},\omega)$ are found to satisfy

[TABLE]

for fixed $i$ . Here we have defined the transport coefficient

[TABLE]

which is due to N-processes and which plays the role of a kinematic viscosity of the phonon gas (see second sound regime and Poiseuille flow below). In writing down Eq. (80) we have neglected an additional transport coefficient which couples the phonon drift velocity to the lattice deformations $s_{i}$ since there is already a coupling to $s_{i}$ by the thermal tension in Eq. (78). Defining the matrix

[TABLE]

we obtain

[TABLE]

Though Eq. (83) would be an appropriate tool to discuss the frequency dependence of the thermal conductivitymajee without resort to the Callaway modelcallaway , we will not pursue this path here but restrict ourselves to the study of second sound and phonon Poiseuille flow.

In order to simplify the algebra, we will only retain the diagonal elements of the matrix $\Lambda_{ik}$ and restrict ourselves to a 2D crystal with hexagonal symmetrychaikin . Then second rank tensors will reduce to scalars viz. $\lambda^{\prime}_{ij}=\lambda^{\prime}\delta_{ij}$ , $\beta_{ij}=\beta\delta_{ij}$ , $d_{ij}=d\delta_{ij}$ , and $\omega^{i}\equiv\omega_{V}$ , $i=\in\{1,2\}$ , and Eq. (83) becomes

[TABLE]

which shows that the phonon drift velocity is damped by crystal momentum dissipating processes $\omega_{V}$ and by N-processes ( $q^{2}\pi^{\prime}$ ). In the following we will study the system of coupled equations (78) and (84) for two distinct hydrodynamic frequency regimes: (i) $\omega_{V}<\omega<\omega^{l}$ , and (ii) $\omega<\omega_{V}<\omega^{l}$ .

Case (i) corresponds to the frequency window conditionGM2 , necessary for the existence of second sound. Substituting $a^{k}(\vec{q},\omega)$ into Eq. (78), we get for a 2D crystal with hexagonal symmetry the wave equation

[TABLE]

where

[TABLE]

has the dimension of velocity squared and $q^{2}=q_{1}^{2}+q_{2}^{2}$ . In obtaining this result we have neglected quadratic contributions in the dissipative terms such as $\omega_{V}\lambda^{\prime}$ , $\pi^{\prime}\lambda^{\prime}$ as well as $\omega_{V}\beta$ , $\pi^{\prime}\beta$ . In analogy with the corresponding phenomenon originally studied in superfluid Heliumkhalatnikov , a wave equation for temperature or equivalently phonon density fluctuations is called second sound. In Eq. (85) the squared second sound velocity $V_{\theta}$ accounts for the restoring forces while the damping is determined by the relaxation frequency $\omega_{V}$ and the dissipation coefficients $\lambda^{\prime}$ and $\pi^{\prime}$ that are inversely proportional to the relaxation frequencies $\omega^{l}$ for N-processes. Besides the external temperature source, dynamic lattice deformations act as periodic temperature perturbations. The converse mechanism is readily studied by inserting the ansatz (75) into Eq. (11) for the lattice displacements. Applying again the method of successive determination of $a^{i}$ and eliminating $a^{l}$ we obtain as result an equation similar to Eq. (68) where $\theta(\vec{q},\omega)$ is determined by Eq. (85). Transforming Eq. (85) to a space and time dependent equation we get

[TABLE]

The corresponding equation for the lattice deformation $s_{i}(\vec{r},t)$ is of the same form as Eq. (72) where now $\theta(\vec{r},t)$ satisfies Eq. (87). This scenario will from now on be addressed as the second sound regime (SSR). In Fig. 1 we plot the second sound velocity $V_{\theta}$ calculated for graphene as a function of temperature $T$ . It ranges from $\sim 550$ ms*-1* at $T=1$ K to $\sim 3.6$ kms*-1* at $T=430$ K and increases monotonically with temperature, in agreement with results of Lee et al.lee_naturecom .

Turning to the case (ii) we take $\omega=0$ and fix the external heat sources such that a constant temperature gradient is maintained in the crystalmezhov ; gurzhi . In real space Eq. (84) leads to

[TABLE]

This equation describes the steady-state flow of a phonon gas driven by a temperature gradient. In analogy with the flow of a fluid of material particles in a confined space under the influence of a pressure differencelandauhydro , Eq. (88) describes the Poiseuille flow of a phonon gassussmann ; gurzhi ; guyer1 . To be specific we consider a 2D hexagonal crystal with rectangular shape of length $L\to\infty$ and width $w$ . A constant temperature gradient is applied along $L$ taken as $\vec{x}$ direction. Then $\vec{a}$ is solely a function of $y$ and $\pi^{\prime}_{ki,ki}$ reduces to $\pi^{\prime}_{12,12}=\pi^{\prime}_{66}$ . We use Voigt’s notation with $(11)\equiv 1$ , $(22)\equiv 2$ and $(12)\equiv 6$ . We then solve the differential equation (88) with boundary conditions $\vec{a}(y=\pm w/2)=0$ . Finally we average over the width $w$ and obtain

[TABLE]

where $\alpha_{w}=\sqrt{\omega_{V}/\pi^{\prime}_{66}}$ has the dimension of inverse length. Using dimensional arguments we identify $c/\omega_{V}$ and $\pi^{\prime}_{66}/c$ , where $c$ is of the order of the sound velocity, with the mean free paths $l_{V}$ and $l_{N}$ for momentum dissipating processes and N-processes respectively. These definitions were previously introduced by Gurzhigurzhi . Noting that the corresponding heat current per unit cell reads

[TABLE]

we define the thermal conductivity $\kappa_{2D}$ per unit cell from comparison of Eqs. (89) and (90) as

[TABLE]

The thermal conductivity per unit volumebalandin1 ; seol $\kappa$ is obtained by division of $\kappa_{2D}$ by the effective volume $v_{3D}=v_{2D}\times h$ , where $h$ is the interlayer distance in the corresponding 3D material:

[TABLE]

This expression provides a unified description of the temperature variation of the thermal conductivity. We notice that $\alpha_{w}$ is determined by both momentum dissipating processes and by N-processes through $\omega_{V}$ and $\pi^{\prime}_{66}$ respectively. Anticipating the results of Secs. (V) and (VI) where we show that $\pi^{\prime}_{66}$ is a strongly increasing function with decreasing temperature, and assuming a sample of finite width such that $\alpha_{w}w<2$ , we see that $\alpha_{w}\to 0$ at low $T$ . Then expansion of $\textrm{tanh}(\frac{\alpha_{w}w}{2})$ in Eq. (92) leads to

[TABLE]

Here $\kappa$ is independent of U-processes and intrinsic impurity scattering but limited by the scattering of N-processes against the sample boundaries. At high $T$ where $\alpha_{w}w>>2$ , Eq. (92) gives

[TABLE]

A quantitative discussion of $\kappa$ as a function of $T$ requires the explicit knowledge of $\omega_{V}$ and $\pi^{\prime}_{66}$ , and will be given in Sec. VI.

We observe that the results of this subsection have been obtained for $\omega_{V}<\omega^{l}$ . At high $T$ or with large impurity scattering this condition breaks down. At low frequencies and long wavelengths such that $\omega$ and $q^{2}\pi^{\prime}$ can be neglected in comparison with $\omega_{V}$ in the denominator of Eq. (84), we find that Eq. (78) reduces to Eq. (69) with $\lambda$ given by Eq.(71), i.e. we recover the situation of subsection III.1.

We conclude by showing how to calculate $\omega_{V}$ without impurity scattering. Starting from Eq. (58) with $\nu(\vec{k},\alpha)=k_{i}$ we obtain for the relaxation frequency $\omega_{U}$ due to U-processes:

[TABLE]

Symmetrizing the first term on the right hand side with respect to $(\vec{p},j)$ and $(\vec{h},\zeta)$ and re ordering terms,

[TABLE]

Typical calculations of thermal conductivity by ab initio methodslee_naturecom ; cepe_naturecom ; xunhydro involve also all combinations of (only) in-plane scattering processes. For a quantitative comparison, the inclusion of in-plane scattering to $\omega_{U}$ can be done by adding to the right hand side of Eq. (105) the terms

[TABLE]

where

[TABLE]

With this addition the relaxation frequency $\omega_{U}$ will be referred as $\tilde{\omega}_{U}$ instead. The expression for $\Phi^{(3)}(\vec{k},i;\vec{p},j;\vec{h},k)$ is similar to Eq. (9), except for the interatomic potentials $\phi^{(3)}_{ijk}(A_{1};B_{s})$ , that refer to in-plane anharmonic interactions. For more details on the specific values and symmetries see Ref. (seba, ).

IV Temperature and Displacement Response

Having obtained the hydrodynamic solutions to Eqs. (11) and (12), we now investigate their implications on the structure of the correlations functions. Such a study is also of relevance in view of experimentsgriffin .

In the following we will restrict ourselves again to 2D crystals of hexagonal symmetry. Then there are only two independent elastic constants $C_{11,11}$ and $C_{11,22}$ , while $2C_{12,12}=C_{11,11}-C_{11,22}$ . Using Voigt’s notation and observing that in 2D crystals these quantities have the dimension of tension coefficients, we write $\gamma_{11}\equiv C_{11,11}$ , $\gamma_{12}\equiv C_{11,22}$ , and $\gamma_{66}\equiv C_{12,12}$ . The other fourth rank tensors $\eta_{ij,kl}$ and $\pi^{\prime}_{ij,kl}$ also have the same symmetries as the elastic constants. In order to calculate the dynamic structure functions $S_{\gamma}(\vec{q},\omega)$ with $\gamma=\{ss,\theta\theta\}$ , we first determine the dynamic susceptibilities $\chi_{\gamma}(\vec{q},\omega)$ , take the imaginary parts $\chi^{\prime\prime}_{\gamma}(\vec{q},\omega)$ and use the fluctuation dissipation theoremcallen ; kadanoff

[TABLE]

with $n(\omega)=[\textrm{exp}(\hbar\omega/k_{B}T)-1]^{-1}$ . In accordance with subsections III.1 and III.2 we will consider separately the diffusive regime and the second sound regime.

IV.1 Diffusive regime

In order to calculate the displacement-displacement response function $S_{ss}$ we consider Eqs. (68) and (69) in the absence of the external temperature source $\Theta(\vec{q},\omega)$ . Taking $\vec{q}=\{q_{1},0\}$ we solve Eq. (69) with respect to $\theta(\vec{q},\omega)$ and then substitute the result into Eq. (68). For the longitudinal case $\vec{s}=\{s_{1},0\}$ we obtain:

[TABLE]

Here we have defined

[TABLE]

where $c_{L}$ is the isothermal sound velocity and $\beta^{\prime 2}$ accounts for the coupling to local temperature fluctuations. Eq. (118) describes damped longitudinal sound waves. Beside the kinematic viscosity $\hat{\eta}_{11}$ , temperature diffusion acts as a mechanism of sound absorption.

Differentiating Eq. (118) on both sides with respect to the external force $F_{1}(\vec{q},\omega)$ we obtain the longitudinal displacement-displacement susceptibility

[TABLE]

The resonances can be analyzed by hand in important limiting cases. At large sound wave frequencies where temperature fluctuations can not follow, i. e. $\omega=c_{L}q_{1}>q_{1}^{2}\lambda$ , the heat diffusion pole at $\omega=-iq_{1}^{2}\lambda$ can be neglected and the resonances are damped sound waves centered at $\pm c_{L}^{ad}q_{1}$ , where

[TABLE]

is the adiabatic sound velocity. In the opposite regime, i. e. $\omega<q_{1}^{2}\lambda$ , the heat diffusion pole leads to an additional resonance at $\omega=0$ , the Landau-Placzek peakbls . The relative strength of these resonances is most conveniently studied by observing that the denominator of Eq. (120) is a polynomial of third order in $\omega$ . Finding the roots by performing a perturbation expansion in powers of small $q_{1}$ and making subsequently a partial fraction decomposition we obtain

[TABLE]

The first two terms on the right hand side of Eq. (122) reflect the sound wave resonances while the third term is due to temperature diffusion. We next determine the imaginary part $\chi^{\prime\prime}_{ss}(\vec{q},\omega)$ and calculate the sum ruleGM2

[TABLE]

The first term on the right hand side has its origin in the Brillouin doublet and the second term is due to the heat diffusion peak. We see that as a consequence of the coupling of lattice deformations to temperature diffusion, the elastic sum rule is not completely exhausted by the isothermal sound waves but there is in addition a contribution of weight $(\beta^{\prime}/c_{L}^{ad})^{2}$ due to the central Landau-Placzek peak. In the left panel of Fig. 2 we have plotted the dynamic structure function $S_{ss}(\vec{q},\omega)$ for a hypothetical 2D crystal with the parameters listed in Table 1.

As we can observe, the Brillouin doublet and the Landau-Placzek peak are clearly present.

In order to calculate the dynamic temperature-temperature susceptibility, we put $F_{i}(\vec{q},\omega)=0$ in Eq. (68) and eliminate $s_{j}$ from the coupled system of Eqs. (68) and (69). Differentiating the resulting equation for $\theta(\vec{q},\omega)$ with respect to the external temperature source $\Theta(\vec{q},\omega)$ we obtain the susceptibility

[TABLE]

with

[TABLE]

Repeating the same steps made above for $\chi_{ss}(\vec{q},\omega)$ , we obtain again three resonances: a sound wave doublet at $\pm c_{L}^{ad}q_{1}$ and a heat diffusion peak at $\omega=0$ . The difference with $\chi_{ss}(\vec{q},\omega)$ is that the strength (weight) of the resonances has been reversed: here the heat diffusion peak is the dominating feature, while the sound wave doublet has weight $(\beta^{\prime}/c_{L}^{ad})^{2}$ . The corresponding scattering law $S_{\theta\theta}(\vec{q},\omega)$ is again readily obtained by taking the imaginary part of Eq. (124) and applying relation (117). The result is plotted in the right panel of Fig. 2 where again we have used the parameters listed in Table 1.

IV.2 Second sound regime

Here we start from the wave equation (85) for local temperature fluctuations. As already noticed at the end of Sec. III, the wave equation for the lattice deformations is formally the same as Eq. (68). We follow the same steps as subsection IV.1, now applied to Eqs. (68) and (85). The longitudinal displacement-displacement susceptibility is readily obtained as

[TABLE]

where

[TABLE]

The resonances of $\chi_{ss}(\vec{q},\omega)$ are now given by the zeros of a polynomial of fourth order in $\omega$ . In the limit of small wave vectors we apply again perturbation theory and perform a partial fraction decomposition. The result reads

[TABLE]

with

[TABLE]

The symbol $\tilde{c}_{L}$ stands for the first sound velocity renormalized by thermal tension due to heat waves (second sound), while $\tilde{V}_{\theta}$ is the second sound velocity renormalized by the presence of coherent dynamic lattice deformations (first sound waves). From Eq. (127) we see that the resonances are given by a sound wave doublet centered at $\pm\tilde{c}_{L}q_{1}$ and a second sound doublet at $\pm\tilde{V}_{\theta}q_{1}$ . The weight of these resonances in the dynamic displacement-displacement susceptibility or equivalently in the dynamic structure function $S_{ss}(\vec{q},\omega)$ is given by $1-\beta^{\prime 2}V^{2}_{\theta}/(c_{L}^{2}-V^{2}_{\theta})^{2}$ and $\beta^{\prime 2}V^{2}_{\theta}/(c_{L}^{2}-V^{2}_{\theta})^{2}$ , respectively. In the right panel of Fig. 3 we have plotted the dynamic structure function $S_{ss}(\vec{q},\omega)$ for a 2D crystal in the second sound regime with the parameters listed in Table (2). According to the partial fraction decomposition in Eq. (127) there are four distinguishable peaks. The doublet with larger intensity appears at $\pm\tilde{c}_{L}q_{1}\approx\pm 1.06\>c_{L}q_{1}$ while the smaller one, closer to $\omega=0$ , appears at $\pm\tilde{V}_{\theta}q_{1}\approx\pm 0.93\>V_{\theta}q_{1}$ . As we can see, the renormalization barely changes the position of the first and second sound peaks (with the a priori chosen parameters) and so the main difference with the diffusive regime is the splitting of the central Landau-Placzek peak into a second sound doublet.

The dynamic temperature-temperature susceptibility is obtained by differentiating Eq. (85) for $\theta(\vec{q},\omega)$ with respect to the external temperature source $\Theta(\vec{q},\omega)$ , after having eliminated $s(\vec{q},\omega)$ by means of Eq. (68). The result reads

[TABLE]

where

[TABLE]

The resonances of $\chi_{\theta\theta}(\vec{q},\omega)$ are given by the zeros of the same polynomial of fourth order in $\omega$ as in the case of $\chi_{ss}(\vec{q},\omega)$ . Carrying out again a partial fraction decomposition we obtain

[TABLE]

Comparison with Eq. (127) for $\chi_{ss}$ shows that the (first) sound doublet has now the smaller weight $\beta^{\prime 2}c_{L}^{2}/(c_{L}^{2}-V^{2})^{2}$ while the second sound doublet has the larger weight $1-\beta^{\prime 2}c_{L}^{2}/(c_{L}^{2}-V^{2})^{2}$ . This trend can be observed clearly in the right panel of Fig. 3 where we plot the values of $\chi_{\theta\theta}(\vec{q},\omega)$ in the second sound regime for the same parameters listed in Table. 2.

Here a remark on the $T$ -dependence of the thermal tension $\beta$ is on order. At low $T$ where flexural modes are dominant, the thermal expansion $\alpha_{T}$ or equivalently $\beta$ is negative in 2D crystalsliftshitz ; mounet as has been found by experiments in grapheneyoon ; bao . At higher $T$ the excitation of in-plane phonons gives in addition a positive contribution to $\alpha_{T}$ and $\beta$ which even can lead to a change of sign as is suggested by atomistic Monte Carlo simulationszakhar . This effect is materialsevic and size dependentseba .

V Transport coefficients

In Sec. III we derived a series of hydrodynamic equations by applying a formal solution method to the kinetic equations. Thereby the non-equilibrium phonon density has been expanded in terms of eigenfunctions of the collision operator $\hat{C}$ . The corresponding eigenvalues are the relaxation frequencies of the crystal. Separating the spectrum into two classes corresponding to slowly varying secular variables (conserved and quasi conserved quantities) and fast varying variables respectively, we have used perturbation theory to formulate the hydrodynamic equations for the secular variables such as lattice deformations and local temperature. In addition to restoring forces that account for the oscillatory part (e.g. elastic constants and sound velocities), the hydrodynamic equations contain dissipative terms that account for the transformation of kinetic energy into heat. The dissipative terms are characterized by the so-called kinetic coefficients or transport coefficients such as viscosity and thermal conductivity.

The transport coefficients $\eta$ , $\lambda$ , $\lambda^{\prime}$ and $\pi^{\prime}$ have been formulated in terms of projections of the currents of secular variables onto the eigenfunctions $\chi^{l}$ belonging to the fast variables of the collision operator. While formally compact, expressions (70), (71) and (81) are not very practical, since neither $\chi^{l}$ nor the corresponding eigenfrequencies $\omega^{l}$ are explicitly known. Here we will follow a more systematic approach where transport coefficients are written as correlation functions of currents of secular variablesKubo ; mori in the limit of zero frequencykadanoff . Taking in turn the currents as secular variables, Götze and one of the present authorsGM3 have expressed the correlation functions in terms of memory kernels for which closed expressions in form of multiple integrals over the Brillouin zone are obtained.

As shown in Ref. GM3, , a transport coefficient $\Lambda$ that describes the decay of a secular variable $A$ (or a set of secular variables) can be cast into the form of a matrix equation

[TABLE]

where the factors on the right hand side can be calculated by perturbation theory. The tensor properties of the factors depend on the nature of the transport coefficient. Here $\chi^{JJ}(0)$ is the static current-current susceptibility, $J$ stands for the non-secular part of the current of $A$ , and $\chi^{AA}(0)$ is the static $A$ - $A$ susceptibility. The memory function $\hat{m}^{KK}(i0)$ is the zero frequency limit $\omega\to+i0$ of a higher hierarchy current-current correlation function, where the current $K$ is obtained from $J$ taken in turn as a secular variable. Closed expressions of the memory function are obtained by evaluating the K-K correlation functions by means of decoupling techniqueszubarev ; GM3 and by taking subsequently the zero frequency limitkadanoff . See Eqs. (140), (150) and (161) below.

In the following, we will present the main steps of implementation of Eq. (131) for the calculation of the transport coefficients $\hat{\eta}_{11}$ , $\lambda$ , $\lambda^{\prime}$ and $\pi^{\prime}_{66}$ . For additional details of the method, see Ref. (GM3, ). In order to calculate the in-plane kinematic viscosity $\hat{\eta}_{11}$ , we take the in-plane momentum $\vec{p}(\vec{q})$ , Eq. (5) conjugate to the in-plane deformation $\vec{s}(\vec{q})$ as secular variable. Then Eq. (131) leads to

[TABLE]

The indexes $p$ , $q$ , $n$ and $l$ here correspond to Voigt’s notation: $\hat{\eta}_{11,11}\equiv\hat{\eta}_{11}$ , $\hat{\eta}_{22,22}\equiv\hat{\eta}_{22}$ , $\hat{\eta}_{11,22}\equiv\hat{\eta}_{12}$ and $\hat{\eta}_{12,12}\equiv\hat{\eta}_{66}$ . These definitions apply to all fourth rank tensors, including current-current susceptibilities and memory kernels as well. We now specify all factors that occur in $\eta_{pq}$ . The momentum-momentum static susceptibility is given by $\chi_{1}^{pp}(0)=1$ . The non secular part of the current of $\vec{p}(\vec{q})$ reads

[TABLE]

where

[TABLE]

withSMP

[TABLE]

Here $h^{(3)}$ is the anharmonic coupling constantseba between in-plane and out-of-plane displacements, and $a$ is the lattice constant. The current-current susceptibility reads

[TABLE]

Symmetry implies $\chi_{11}^{JJ}(0)=\chi_{22}^{JJ}(0)$ and $\chi_{12}^{JJ}(0)=\chi_{21}^{JJ}(0)$ . The subtraction of the term proportional to $\hbar\omega(\vec{k},\zeta)$ on the right hand side of Eq. (134) corresponds to the absence of the $l=0$ term in Eq. (70). The elements of the matrix $\hat{m}^{KK}(i0)$ are

[TABLE]

where $\omega(\vec{k},i)$ with $i=\{1,2\}$ are the in-plane phonon frequencies and $\Psi(\vec{k},i;\vec{k}_{1},\zeta;\vec{k}_{2},\zeta)$ has been defined by Eq. (49). While a quantitative evaluation of $\hat{\eta}$ and the other kinetic coefficients is only possible by numerical calculations (see next section), the asymptotic behavior as a function of temperature has to be obtained by analytical means. Assuming a linear dispersion for the in-plane phonon frequencies and a quadratic dispersion for the flexural mode, we have studied the $T$ -dependence of the various susceptibilities and memory kernels. We then find by means of Eqs. (136) and (140) that in the limit $T\to 0$ the susceptibility $\chi^{JJ}(0)$ diverges while $m^{KK}(i0)$ vanishes $\propto T$ . Hence we conclude that $\hat{\eta}$ diverges with $T\to 0$ . At high $T$ we find that $\chi^{JJ}(0)\propto T$ and $m^{KK}(i0)\propto T^{2}$ , hence $\hat{\eta}$ tends to a constant value.

Next we consider the thermal diffusion coefficient $\lambda_{ij}=\lambda\delta_{ij}$ which has been introduced in Sec. III, Eq. (71) The corresponding secular variable in this case is the harmonic phonon energy density

[TABLE]

Expression (131) leads to the thermal diffusion coefficient

[TABLE]

where all factors on the right hand side are scalars. The energy-energy density susceptibility is given by

[TABLE]

and the corresponding current reads

[TABLE]

with

[TABLE]

The current-current susceptibility reads

[TABLE]

and the memory function is obtained as

[TABLE]

Proceeding as outlined above, we find by analytical reasoning that at low $T$ the susceptibility $\chi^{\epsilon\epsilon}(0)$ and $\chi^{JJ}(0)$ vary as $T^{2}$ and $T^{3}$ respectively and that $m^{KK}(i0)\propto T^{3}$ . Hence we conclude that the thermal diffusion $\lambda$ vanishes $\propto T$ for $T\to 0$ . At high $T$ we obtain $\chi^{\epsilon\epsilon}\propto T$ , $\chi^{JJ}(0)\propto T$ and $m^{KK}(i0)\propto T^{2}$ , which means that $\lambda$ vanishes as $T^{-1}$ .

Finally we consider the kinetic coefficients $\lambda^{\prime}$ and $\pi^{\prime}$ introduced in Sec. III. The calculation of $\lambda^{\prime}$ requires the subtraction of the crystal momentum as a secular variable from the energy current. Eq. (145) has to be replaced by

[TABLE]

where $g^{\lambda}=\langle k_{1}v_{1}|\omega\rangle/\langle k_{1}|k_{1}\rangle$ . The kinetic coefficient $\lambda^{\prime}$ is then obtained by replacing $\gamma_{1}$ by $\gamma^{\prime}_{1}$ in the expression for $\chi^{JJ}$ and $\hat{m}^{KK}$ that enter Eq. (142) for $\lambda$ . The asymptotic $T$ -behavior of $\lambda^{\prime}$ is the same as that for $\lambda$ .

The kinetic coefficient $\pi^{\prime}_{ij,kl}$ defined in Eq. (81) plays the role of kinematic viscosity of the phonon gas and the secular variable is the crystal momentum density

[TABLE]

Since $\pi^{\prime}$ is a fourth rank tensor it has the symmetry of the elastic constants. Using again Voigt’s notation we may write $\pi^{\prime}_{pq}$ as:

[TABLE]

The crystal momentum susceptibility reads

[TABLE]

The momentum currents are given by

[TABLE]

with

[TABLE]

where $g^{\pi}_{11}=g^{\pi}_{22}=\langle\omega|k_{1}v_{1}\rangle/\langle\omega|\omega\rangle$ . The static current-current susceptibilities are

[TABLE]

and the relevant elements of the memory function are obtained as

[TABLE]

At low $T$ we find by analytical means that in leading order of $T$ , $\chi^{\pi\pi}(0)\propto T$ , $\chi^{JJ}(0)\propto T$ and $m^{KK}(i0)\propto T^{3}$ . Hence $\pi^{\prime}$ diverges as $T^{-2}$ for $T\to 0$ . At high $T$ we get $\chi^{\pi\pi}(0)\propto T$ , $\chi^{JJ}(0)\propto T$ and $m^{KK}(i0)\propto T^{2}$ , and we conclude that $\pi^{\prime}$ vanishes as $T^{-1}$ . Here and before high $T$ means $T$ larger then $\hbar\omega_{D}/k_{B}$ where $\omega_{D}$ is a BZ boundary frequency.

Notice that in obtaining the results Eqs. (140), (150) and (140) we do not make use of the concept of a wave-vector dependent relaxation time. The multiple integral expressions are reminiscent of transport theory results obtained by means of variational methodskoler ; leibfriedscho ; ziman , however in the present approach we do not resort to any trial function. A numerical evaluation of the transport coefficients is then reduced to a calculation of the single and multiple $\vec{k}$ -integrals occurring in the factors $\chi^{AA}(0)$ , $\chi^{JJ}(0)$ and $m^{KK}(i0)$ . Such a program will be carried out in the next section where we study the temperature evolution of the transport coefficients.

VI Temperature dependence of transport coefficients

As follows from the analytical results of Secs. III and V the evaluation of the transport coefficients $\eta$ , $\lambda$ , $\pi^{\prime}$ and the thermal conductivity $\kappa$ can be achieved by the calculation of the corresponding susceptibilities and memory functions. For this purpose we need to calculate simple two-folded and triple two-folded integrals inside the BZ. These integrals involve the harmonic frequencies $\omega(\vec{k},\alpha)$ , the group velocities $\vec{v}(\vec{k},\alpha)$ , the anharmonic interactions $\Phi(\vec{k}_{1},\alpha_{1};\vec{k}_{2},\alpha_{2};\vec{k}_{3},\alpha_{3})$ , and the momentum and energy conservation functions $\Delta(\vec{k}_{1},\vec{k}_{2},\vec{k}_{3},\vec{G})$ and $\delta(\omega(\vec{k}_{1},\alpha_{1}),\omega(\vec{k}_{2},\alpha_{2}),\omega(\vec{k}_{3},\alpha_{3}))$ . We will use the corresponding input parameters for grapheneseba .

Following the steps of Ref. (seba, ) we first generate the hexagonal BZ for graphene and then for each of the $N$ $\vec{k}$ -points in the BZ we calculate the acoustic phonon dispersions by direct diagonalization of the dynamical matrixmohr . Because of band-crossing, we also make use of an auxiliary algorithm to keep the correct sorting of phonon modes after diagonalization.

The evaluation of simple two-folded integrals for the calculation of scalar products, e.g. $\beta$ and $d$ , or susceptibilities, e.g. $\chi^{JJ}(0)$ and $\chi^{\epsilon\epsilon}(0)$ , is performed in a straightforward way by summing over all $\vec{k}$ -points and the three acoustic polarizations. In Figs. 1, 4, 5 and 6 the number of $\vec{k}$ -points used is $N=25600$ . Triple two-folded integrals are required specifically for memory functions $\hat{m}^{KK}(i0)$ . In order to evaluate them we generate a set of all possible combinations $\{\vec{k}_{1},\vec{k}_{2},\vec{k}_{3}\}$ from the original list of $\vec{k}$ -points in the BZ. Among all these combinations, those which satisfy $\vec{k}_{1}+\vec{k}_{2}+\vec{k}_{3}=0$ belong to the N-processes group and those which satisfies $\vec{k}_{1}+\vec{k}_{2}+\vec{k}_{3}=\vec{G}$ , for any reciprocal vector $\vec{G}$ , belong to the U-processes group. The explicit distinction between normal and umklapp groups allows us to determine their contributions separately, as it will be shown below. Energy conservation is enforced by means of a Gaussian function $(\epsilon\sqrt{\pi})^{-1}\>\textrm{exp}[(-\omega(\vec{k}_{1},\alpha_{1})+\omega(\vec{k}_{2},\alpha_{2})+\omega(\vec{k}_{3},\alpha_{3}))/\epsilon]^{2}$ , where a broadening parameter $\epsilon=5$ cm*-1* is chosen. Memory functions are then evaluated by summing over normal and umklapp groups.

Left panel of Fig. 4 shows in solid blue line the temperature evolution of the dynamic viscosity $\eta_{11}$ calculated for graphene according to Eqs. (132 - 140). Because of the even $\vec{k}$ -dependence of the thermoelastic coupling $h_{ij}(\vec{k},\zeta)$ , the scalar products $\langle h_{ij}|\chi^{k}\rangle$ for $k=\{1,2\}$ in Eq. (70) are zero. This makes the summation condition $l>2$ equivalent to $l>0$ , resulting in $\eta_{ij}$ having the same value in both diffusive and second sound regimes. As observed the dynamic viscosity shows a steep increase with $T\to 0$ in accordance with the $T$ -dependence discussion in the previous section and then presents a local minimum at $T\sim 140$ K. The reason of this local minimum is related to the temperature dependence of the thermal tension $\beta$ , which remains negative according to Eq. (63) (the change in sign of $\beta$ at higher temperatures is a direct consequence of in-plane scatteringseba ). Right panel shows the values of the thermal diffusion coefficients $\lambda$ and $\lambda^{\prime}$ calculated according to Eqs. (141 - 151). In addition we plot the corresponding N-processes contributions to $\lambda$ and $\lambda^{\prime}$ . As expected the low $T$ behavior in both the diffusion and second sound regimes is dominated entirely by N-processes. The contribution of U-processes becomes only noticeable above $T\sim 60$ K.

The kinematic viscosity coefficient $\pi^{\prime}_{66}$ for graphene, that enters in the expression for the thermal conductivity $\kappa$ Eq. (92), is plotted in the left panel of Fig. 5. According to Eqs. (152 - 161), this coefficient presents a steep increase with $T\to 0$ and a monotonically decreasing behavior for increasing $T$ . Right panel of Fig. 5 shows the temperature evolution of the relaxation frequency $\omega_{U}$ calculated according to Eq. (105) and the values of $\tilde{\omega}_{U}$ , i.e. the relaxation frequency $\omega_{U}$ with the addition of the in-plane scattering terms Eqs. (111 - 116). At low $T$ both frequencies are essentially the same, meaning that in-plane scattering can be completely neglected. With $T$ increasing above $\sim 150$ K the relative difference between $\omega_{U}$ and $\tilde{\omega}_{U}$ becomes noticeable and in particular it reaches a value of $\sim 70\%$ at $T=300$ K.

Fig. 6 shows the thermal conductivity $\kappa$ calculated for graphene according to Eq. (92) for three different widths $w=100$ , $50$ , and $20$ $\mu$ m. At low $T$ , $\kappa$ behaves as predicted by Eq. (93) and tends to zero with $\pi^{\prime}_{66}\to\infty$ . With increasing $T$ , $\pi^{\prime}_{66}$ decreases and, as far as the relaxation frequency $\omega_{V}$ of crystal momentum destroying processes is small, $\kappa$ increases. At high $T$ where U-processes are dominant (we assume that impurity scattering is negligible and hence $\omega_{V}\equiv\omega_{U}$ ), $\kappa$ decreases with increasing $T$ , in accordance with Eq. (94). In the intermediate temperature regime $\kappa$ reaches a maximum that is determined by the subtle interplay of the $T$ -dependence of $\pi^{\prime}_{66}$ and $\omega_{U}$ . The maximum shifts to higher $T$ with lower sample width $w$ , in agreement with Refs. fugallo, ; xunhydro, .

A quantitative comparison with the results of other workscepe_naturecom ; lee_naturecom ; fugallo ; xunhydro shows that our values of $\kappa$ for $T$ above $150$ K are larger. Also the inclusion of in-plane scattering does not resolve this discrepancy. A further reduction of the values of $\kappa$ requires the inclusion of impurity scattering and possibly of fourth order scattering processes of flexural (ZA) phononsfourthdegradation . We expect that the latter lead to a reduction of $\pi^{\prime}_{66}$ or equivalently to a decrease of the mean free path $l_{N}$ for N-processes.

VII Conclusions

Starting from coupled dynamic equations for in-plane lattice displacements (11) and flexural phonon kinetics (11) we have used perturbative methods to obtain the hydrodynamic equations describing elastic sound waves and local temperature fluctuations. We have distinguished two scenarios: (i) if phonon energy is the only conserved quantity, temperature fluctuations are described by a diffusion equation; (ii) if in addition crystal momentum is taken as an almost conserved quantity and the frequency window condition is fulfilledprohofsky , temperature fluctuations are described by a wave equation, the so-called second sound. In the zero frequency limit and in presence of a constant temperature gradient, scenario (ii) leads to Poiseuille flow.

The present study differs from previous theoretical works on the thermal properties of 2D crystals (see references in Sec. I) in several aspects. We start from coupled dynamic equations derived from a microscopic HamiltonianSMP and obtain therefrom a unified description for elastic and thermal hydrodynamic phenomena. The coupling is mediated by the elastic tension coefficient which is an anharmonic effect due to the coupling between in-plane and flexural lattice displacements. By investigating the elastic and thermal response functions, we have shown how the thermal resonances (second sound or Landau-Placzek peak) appear in the displacement-displacement response function while conversely the sound wave doublet is present in the thermal response. The interplay of elastic and thermal phenomena in various scattering laws is a challenge for new experiments. Indeed most recently second sound has been observed in graphite at temperatures above 100 K by means of time-resolved optical measurementsssgraphite . The method is based on the transient thermal grating techniquettgrating . The spatial and temporal decay of the thermal grating by second sound heat transport is reflected, as a consequence of thermo-mechanical coupling, in the surface displacement field that acts as a transient diffraction grating.

A further distinct aspect of our work is the calculation of various transport coefficients such as the first sound viscosity $\eta$ , the thermal diffusion $\lambda$ and the phonon viscosity $\pi^{\prime}$ . Together with the crystal momentum relaxation frequency $\omega_{V}$ these coefficients characterize the broadening of elastic and thermal resonances and determine the strength of the thermal conductivity $\kappa$ . Using Kubo-Mori response theory we have written the transport coefficients as relaxation functions of generalized currents of secular variables and calculated the current-current relaxation functions by means of an equation of motion method due to Götze and one of the present authorsGM3 . Thereby we have obtained closed expressions for the transport coefficients in form of multiple integrals over the Brillouin zone that are reminiscent of variational theory results leibfriedscho ; ziman ; sevic . A quantitative evaluation of the transport coefficients and their $T$ -dependence has been performed by numerical calculations in Sec. VI. In addition we have study the asymptotic behavior at low and high temperatures by analytical means. In particular we find that the viscosities $\eta$ and $\pi^{\prime}$ diverge at low $T$ while the thermal diffusion $\lambda$ tends to zero. Special care has been devoted to the calculation of Poiseuille flowgurzhi and the static thermal conductivity $\kappa$ . We have solved the differential equation (88) that describes phonon drift in presence of a static temperature gradient for a rectangular 2D crystal with width w. The thermal conductivity, Eq. (92), then depends on $\omega_{V}$ and $\pi^{\prime}_{66}$ . At low $T$ and for a sample with finite width, it follows from Eq. (93) and from the already mentioned divergence of $\pi^{\prime}$ that $\kappa$ vanishes with $T\to 0$ . At high $T$ where $\omega_{V}$ is large and $\pi^{\prime}_{66}$ vanishes, Eq. (94) shows that $\kappa$ tends asymptotically to zero with increasing $T$ . The complete temperature evolution of $\kappa$ is presented in Fig. (6). We recall that $\pi^{\prime}_{66}$ is due to N-processes, see Eq. (81) and remark following Eq. (89). The vanishing of $\kappa$ is then the ultimate consequence of the disappearance of thermal phonons as heat carriers in the zero temperature limit.

VIII Acknowledgements

This work was supported by the Flemish Science Foundation (FWO-Vl).

Bibliography61

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] L. D. Landau and E. M. Lifshitz. Theory of elasticity . Pergamon Press, London, 1959.
2[2] R. E. Peierls. Zur kinetischen theorie der wärmeleitung in kristallen. Ann. der Physik , 3:1055, 1929.
3[3] R. E. Peierls. Quantum Theory of Solids . International series of monographs on physics. Clarendon Press, 1955.
4[4] H. B. Casimir. Note on the conduction of heat in crystals. Physica , 5(6):495 – 500, 1938.
5[5] P. G. Klemens. Encyclopedia of Physics . Springer Verlag, Berlin, 1956.
6[6] G. Leibfried and E. Schlömann. Wärmeleitung in elektrisch isolierenden kristallen. Nachrichten der Akad. d. Wiss. Göttingen, Math. Phys. Klasse , page 71, 1954.
7[7] J. M. Ziman. Electrons and Phonons: The Theory of Transport Phenomena in Solids . International series of monographs on physics. OUP Oxford, Oxford, 1960.
8[8] J. A. Sussmann and A. Thellung. Thermal conductivity of perfect dielectric crystals in the absence of umklapp processes. Proceedings of the Physical Society , 81(6):1122, 1963.