Unsteady ballistic heat transport in two-dimensional harmonic graphene lattice

Graphene has a regular honeycomb lattice with a unit cell containing two carbon atoms (see figures 1(a) and (b)). Each atom has two independent translational degrees of freedom; as a result, the unit cell has four degrees of freedom. Unit cells are numbered by pairs of indices {n, m}. Position vectors of unit cells {n, m} and {s, p} are related as ⁵

$$\begin{equation}{\mathbf{r}}_{n,m}={\mathbf{r}}_{s,p}+(n-s){\mathbf{b}}_{1}+(m-p){\mathbf{b}}_{2}.\end{equation} \tag{ 1 }$$

Here primitive vectors of the lattice, b₁, b₂, are represented in Cartesian basis i, j as

$$\begin{equation}{\mathbf{b}}_{1}=a\frac{\sqrt{3}}{2}\left(\mathbf{i}+\sqrt{3}\mathbf{j}\right),\qquad {\mathbf{b}}_{2}=a\frac{\sqrt{3}}{2}\left(-\mathbf{i}+\sqrt{3}\mathbf{j}\right),\end{equation} \tag{ 2 }$$

where a is an equilibrium bond length in graphene; in figure 1(b) vector i is horizontal and vector j is vertical.

Zoom In Zoom Out Reset image size

Vectors a₁, ..., a₆ connect two atoms of the unit cell with their nearest neighbors (figure 1(b)). The vectors are introduced such that a₄ = −a₁, a₅ = −a₂, a₆ = −a₃.

We consider the equations of motion for two atoms from the unit cell with indices {n, m}. Each atom has two degrees of freedom. Displacements of atoms from the unit cell {n, m} are denoted as 2D vectors U_n,m , V_n,m . Following [27], we use the harmonic approximation for the potential energy of the lattice. Each atom is connected with three nearest neighbors by linear springs (bonds) with stiffness c. Additionally, the nearest bonds between particles are connected by angular springs with stiffness g. Equations of motion for this system were obtained using the Euler–Lagrange formalism in paper [27]. For the sake of completeness, we derive these equations in the appendix A by entirely different means. The resulting equations are

$$\begin{align}\hfill M{\ddot{U}}_{n,m}^{x}& =\frac{c\sqrt{3}}{4}\left[\sqrt{3}{V}_{n+1,m}^{x}+{V}_{n+1,m}^{y}-2\sqrt{3}{U}_{n,m}^{x}\right.\hfill \\ \hfill & \quad \left.+\sqrt{3}{V}_{n,m+1}^{x}-{V}_{n,m+1}^{y}\right]\hfill \\ \hfill & \quad +\frac{g}{4}\left[2\left(-{U}_{n-1,m}^{x}+\sqrt{3}{U}_{n-1,m}^{y}\right.\right.\hfill \\ \hfill & \quad \left.\left.-{U}_{n,m-1}^{x}-\sqrt{3}{U}_{n,m-1}^{y}\right)+{U}_{n+1,m-1}^{x}+\right.\hfill \\ \hfill & \quad \left.+\sqrt{3}{U}_{n+1,m-1}^{y}+{U}_{n-1,m+1}^{x}-\sqrt{3}{U}_{n-1,m+1}^{y}\right.\hfill \\ \hfill & \quad \left.-2\left({U}_{n+1,m}^{x}+{U}_{n,m+1}^{x}\right)-\right.\hfill \\ \hfill & \quad \left.+6\left({V}_{n+1,m}^{x}-\sqrt{3}{V}_{n+1,m}^{y}-5{U}_{n,m}^{x}\right.\right.\hfill \\ \hfill & \quad \left.\left.+4{V}_{n,m}^{x}+{V}_{n,m+1}^{x}+\sqrt{3}{V}_{n,m+1}^{y}\right)\right],\hfill \\ \hfill M{\ddot{U}}_{n,m}^{y}& =\frac{c}{4}\left[\sqrt{3}{V}_{n+1,m}^{x}+{V}_{n+1,m}^{y}-6{U}_{n,m}^{y}+4{V}_{n,m}^{y}\right.\hfill \\ \hfill & \quad \left.-\sqrt{3}{V}_{n,m+1}^{x}+{V}_{n,m+1}^{y}\right]+\hfill \\ \hfill & \quad +\frac{\sqrt{3}g}{4}\left[{U}_{n-1,m+1}^{x}-\sqrt{3}{U}_{n-1,m+1}^{y}\right.\hfill \\ \hfill & \quad \left.-{U}_{n+1,m-1}^{x}-\sqrt{3}{U}_{n+1,m-1}^{y}+2{U}_{n+1,m}^{x}-\right.\hfill \\ \hfill & \quad \left.-2{U}_{n,m+1}^{x}-10\sqrt{3}{U}_{n,m}^{y}\right.\hfill \\ \hfill & \quad \left.+6\left({V}_{n,m+1}^{x}+\sqrt{3}{V}_{n,m+1}^{y}-{V}_{n+1,m}^{x}\right.\right.\hfill \\ \hfill & \quad \left.\left.+\sqrt{3}{V}_{n+1,m}^{y}\right)\right].\hfill \end{align} \tag{ 3 }$$

$$\begin{align}\hfill M{\ddot{V}}_{n,m}^{x}& =\frac{\sqrt{3}c}{4}\left[\sqrt{3}{U}_{n-1,m}^{x}+{U}_{n-1,m}^{y}-2\sqrt{3}{V}_{n,m}^{x}\right.\hfill \\ \hfill & \quad \left.+\sqrt{3}{U}_{n,m-1}^{x}-{U}_{n,m-1}^{y}\right]+\hfill \\ \hfill & \quad +\frac{g}{4}\left[6\left({U}_{n-1,m}^{x}-\sqrt{3}{U}_{n-1,m}^{y}\right.\right.\hfill \\ \hfill & \quad \left.\left.+{U}_{n,m-1}^{x}+\sqrt{3}{U}_{n,m-1}^{y}+4{U}_{n,m}^{x}\right)-\right.\hfill \\ \hfill & \quad \left.-2\left({V}_{n+1,m}^{x}-\sqrt{3}{V}_{n+1,m}^{y}+{V}_{n,m+1}^{x}\right.\right.\hfill \\ \hfill & \quad \left.\left.+\sqrt{3}{V}_{n,m+1}^{y}\right)\right.\hfill \\ \hfill & \quad \left.+{V}_{n+1,m-1}^{x}-\sqrt{3}{V}_{n+1,m-1}^{y}\right.\hfill \\ \hfill & \quad \left.+{V}_{n-1,m+1}^{x}+\sqrt{3}{V}_{n-1,m+1}^{y}-30{V}_{n,m}^{x}\right.\hfill \\ \hfill & \quad \left.-2({V}_{n,m-1}^{x}+{V}_{n-1,m}^{x})\right],\hfill \\ \hfill M{\ddot{V}}_{n,m}^{y}& =\frac{c}{4}\left[\sqrt{3}{U}_{n-1,m}^{x}+{U}_{n-1,m}^{y}\right.\hfill \\ \hfill & \quad \left.+4{U}_{n,m}^{y}-6{V}_{n,m}^{y}-\sqrt{3}{U}_{n,m-1}^{x}+{U}_{n,m-1}^{y}\right]\hfill \\ \hfill & \quad \times \frac{\sqrt{3}g}{4}\left[2\sqrt{3}\left(-\sqrt{3}{U}_{n-1,m}^{x}+3{U}_{n-1,m}^{y}\right.\right.\hfill \\ \hfill & \quad \left.\left.-5{V}_{n,m}^{y}+\sqrt{3}{U}_{n,m-1}^{x}+3{U}_{n,m-1}^{y}\right)\right.\hfill \\ \hfill & \quad \left.+{V}_{n+1,m-1}^{x}-\sqrt{3}{V}_{n+1,m-1}^{y}-{V}_{n-1,m+1}^{x}\right.\hfill \\ \hfill & \quad \left.-\sqrt{3}{V}_{n-1,m+1}^{y}-2\left({V}_{n,m-1}^{x}-{V}_{n-1,m}^{x}\right)\right],\hfill \end{align} \tag{ 4 }$$

where M is the mass of the carbon atom; c, g are stiffnesses of linear and angular springs respectively. In further calculations, we use the following values of parameters for graphene:

$$\begin{equation}c=730.2\enspace \mathrm{N}\enspace {\mathrm{m}}^{-1}\qquad g=66.9\enspace \mathrm{N}\enspace {\mathrm{m}}^{-1},\qquad a=0.142\enspace \mathrm{n}\mathrm{m}.\end{equation} \tag{ 5 }$$

In paper [28] it is shown that the given values of stiffnesses yield correct in-plane elastic moduli of graphene.

We consider random initial conditions corresponding to a given initial kinetic temperature profile T₀(r_n,m ) (see definition (10)):

$$\begin{align}\hfill {U}_{n,m}^{x}& ={U}_{n,m}^{y}=0,\qquad {V}_{n,m}^{x}={V}_{n,m}^{y}=0,\hfill \\ \hfill {\dot{U}}_{n,m}^{x}& ={\beta }_{n,m}^{x}\sqrt{\frac{{k}_{\mathrm{B}}}{M}{T}_{0}({\mathbf{r}}_{n,m})},\qquad {\dot{U}}_{n,m}^{y}={\beta }_{n,m}^{y}\sqrt{\frac{{k}_{\mathrm{B}}}{M}{T}_{0}({\mathbf{r}}_{n,m})},\qquad \hfill \\ \hfill {\dot{V}}_{n,m}^{x}& ={\gamma }_{n,m}^{x}\sqrt{\frac{{k}_{\mathrm{B}}}{M}{T}_{0}({\mathbf{r}}_{n,m})},\qquad {\dot{V}}_{n,m}^{y}={\gamma }_{n,m}^{y}\sqrt{\frac{{k}_{\mathrm{B}}}{M}{T}_{0}({\mathbf{r}}_{n,m})}.\hfill \end{align} \tag{ 6 }$$

where k_B is the Boltzmann constant. It is assumed that the function T₀ is defined for all points of the graphene plane and slowly changes at distances of order of the interatomic distance. In (6), parameters ${\beta }_{n,m}^{x},{\beta }_{n,m}^{y}$ and ${\gamma }_{n,m}^{x},{\gamma }_{n,m}^{y}$ are uncorrelated random values with zero mathematical expectation and unit variance, i.e.

$$\begin{align}\hfill \left\langle {\beta }_{n,m}^{x}\right\rangle & =\left\langle {\beta }_{n,m}^{y}\right\rangle =\left\langle {\gamma }_{n,m}^{x}\right\rangle =\left\langle {\gamma }_{n,m}^{y}\right\rangle =0,\hfill \\ \hfill \left\langle {({\beta }_{n,m}^{x})}^{2}\right\rangle & =\left\langle {({\beta }_{n,m}^{y})}^{2}\right\rangle =\left\langle {({\gamma }_{n,m}^{x})}^{2}\right\rangle =\left\langle {({\gamma }_{n,m}^{y})}^{2}\right\rangle =1,\hfill \\ \hfill \left\langle {\beta }_{n,m}^{x}{\beta }_{n,m}^{y}\right\rangle & =\left\langle {\gamma }_{n,m}^{x}{\gamma }_{n,m}^{y}\right\rangle =\left\langle {\beta }_{n,m}^{x}{\gamma }_{n,m}^{y}\right\rangle =\left\langle {\gamma }_{n,m}^{x}{\beta }_{n,m}^{y}\right\rangle \hfill \\ \hfill & =0.\hfill \end{align} \tag{ 7 }$$

Here and below 〈‥〉 stands for the mathematical expectation. We note that the initial conditions (6) are such that mathematical expectations of kinetic energies (and temperatures (9)), corresponding to four degrees of freedom of the unit cell are equal.

To define thermodynamic properties of the system, we consider an infinite number of realizations with random initial conditions (6). Initially, kinetic temperatures corresponding to four degrees of freedom of the unit cell are equal. Further it is shown that during ballistic heat transfer these temperatures are generally different. Therefore we introduce the temperate matrix [19], T , such that

$$\begin{equation}{k}_{\mathrm{B}}\boldsymbol{T}=M\left(\begin{matrix}\hfill \left\langle {\dot{U}}^{x2}\right\rangle \hfill & \hfill \left\langle {\dot{U}}^{x}{\dot{U}}^{y}\right\rangle \hfill & \hfill \left\langle {\dot{U}}^{x}{\dot{V}}^{x}\right\rangle \hfill & \hfill \left\langle {\dot{U}}^{x}{\dot{V}}^{y}\right\rangle \hfill \\ \hfill \left\langle {\dot{U}}^{x}{\dot{U}}^{y}\right\rangle \hfill & \hfill \left\langle {\dot{U}}^{y2}\right\rangle \hfill & \hfill \left\langle {\dot{U}}^{y}{\dot{V}}^{x}\right\rangle \hfill & \hfill \left\langle {\dot{U}}^{y}{\dot{V}}^{y}\right\rangle \hfill \\ \hfill \left\langle {\dot{U}}^{x}{\dot{V}}^{x}\right\rangle \hfill & \hfill \left\langle {\dot{U}}^{y}{\dot{V}}^{x}\right\rangle \hfill & \hfill \left\langle {\dot{V}}^{x2}\right\rangle \hfill & \hfill \left\langle {\dot{V}}^{x}{\dot{V}}^{y}\right\rangle \hfill \\ \hfill \left\langle {\dot{U}}^{x}{\dot{V}}^{y}\right\rangle \hfill & \hfill \left\langle {\dot{U}}^{x}{\dot{V}}^{y}\right\rangle \hfill & \hfill \left\langle {\dot{V}}^{x}{\dot{V}}^{y}\right\rangle \hfill & \hfill \left\langle {\dot{V}}^{y2}\right\rangle \hfill \end{matrix}\right).\end{equation} \tag{ 8 }$$

In this subsection, arguments r_n,m , t of temperatures as well as indices n, m of velocities are omitted for brevity.

The diagonal elements of the temperature matrix define the kinetic temperatures, corresponding to four degrees of freedom of the unit cell:

$$\begin{align}\hfill {k}_{\mathrm{B}}{T}_{11}=M\left\langle {\dot{U}}^{x2}\right\rangle ,\qquad {k}_{\mathrm{B}}{T}_{22}=M\left\langle {\dot{U}}^{y2}\right\rangle ,\\ \hfill {k}_{\mathrm{B}}{T}_{33}=M\left\langle {\dot{V}}^{x2}\right\rangle ,\qquad {k}_{\mathrm{B}}{T}_{44}=M\left\langle {\dot{V}}^{y2}\right\rangle .\end{align} \tag{ 9 }$$

We also introduce the conventional (average) kinetic temperature as

$$\begin{equation}T=\frac{1}{4}\left({T}_{11}+{T}_{22}+{T}_{33}+{T}_{44}\right).\end{equation} \tag{ 10 }$$

Under initial conditions (6), all temperatures are equal at t = 0, i.e. T_ii = T = T₀.

The main goals of the present paper are to investigate evolution of the temperature field in the harmonic graphene lattice and to show that during ballistic heat transport the temperatures T_ii are different.

Ballistic heat transfer is carried out by elastic waves traveling in the lattice. Therefore description of the heat transfer requires knowledge of the dispersion relation for the lattice and corresponding group velocities (see e.g. papers [10–12, 14, 19]). To obtain the dispersion relation we seek the solution of equations of motion (3) and (4) in the form of plane waves

$$\begin{align}\hfill {\mathbf{U}}_{n,m}& =\mathbf{U}(\mathbf{k}){\mathrm{e}}^{\mathrm{i}(\omega t-\mathbf{k}\cdot {\mathbf{r}}_{n,m})},\qquad {\mathbf{V}}_{n,m}=\mathbf{V}(\mathbf{k}){\mathrm{e}}^{\mathrm{i}(\omega t-\mathbf{k}\cdot {\mathbf{r}}_{n,m})},\hfill \end{align} \tag{ 11 }$$

where i² = −1; U(k), V(k) are the amplitudes of the waves; ω is a frequency; r_n,m stands for position vector of the unit cell; k is a wave vector represented in a reciprocal basis ${\tilde{\mathbf{b}}}_{i}$ such that ${\tilde{\mathbf{b}}}_{i}\cdot {\tilde{\mathbf{b}}}_{j}=2\pi {\delta }_{ij}$:

$$\begin{align}\hfill \mathbf{k}& ={k}_{1}{\tilde{\mathbf{b}}}_{1}+{k}_{2}{\tilde{\mathbf{b}}}_{2},{\tilde{\mathbf{b}}}_{1} = \frac{2\pi }{3a}\left(\sqrt{3}\mathbf{i}+\mathbf{j}\right),{\tilde{\mathbf{b}}}_{2} = \frac{2\pi }{3a}\left(-\sqrt{3}\mathbf{i}+\mathbf{j}\right).\hfill \end{align} \tag{ 12 }$$

Here and below k₁, k₂ ∈ [0; 1] are dimensionless components of the wave vector. Taking formulae (12) into account, we substitute (11) into (3) and (4). The substitution leads to a homogeneous system of four linear equations with respect to components of the amplitudes. The system has nontrivial solutions only if its determinant is equal to zero. This condition yields quartic equation with respect to ω². Solving the equation numerically, we obtain the dispersion surfaces ω_j (k₁, k₂), j = 1, 2, 3, 4. Plots of the surfaces are presented in paper [27].

Given known the dispersion relation, the group velocities ${\mathbf{v}}_{g}^{j},j=1,2,3,4$ are calculated as

$$\begin{equation}{\mathbf{v}}_{g}^{j}({k}_{1},{k}_{2})=\frac{\partial {\omega }_{j}}{\partial {k}_{1}}{\mathbf{b}}_{1}+\frac{\partial {\omega }_{j}}{\partial {k}_{2}}{\mathbf{b}}_{2},\qquad j=1,\dots ,4.\end{equation} \tag{ 13 }$$

Further, the frequencies ω_j and group velocities ${\mathbf{v}}_{g}^{j}$ are employed for description of the heat transfer.

Analytical description of the ballistic heat transfer, presented e.g. in paper [19], is based on the assumption that the initial temperature field T₀ slowly changes in space on the distances of order of the interatomic distance. Using this assumption the following formulae, describing evolution of the temperature field, are derived:

$$\begin{align}\hfill T(\mathbf{r},t)& ={T}_{\mathrm{F}}(\mathbf{r},t)+{T}_{\mathrm{S}}(\mathbf{r},t),\hfill \\ \hfill & {T}_{\mathrm{F}}(\mathbf{r},t)=\frac{{T}_{0}(\mathbf{r})}{8}\sum\limits _{j=1}^{4}{\int }_{\mathbf{k}}\mathrm{cos}(2{\omega }_{j}t)\mathrm{d}\mathbf{k},\hfill \\ \hfill {T}_{\mathrm{S}}(\mathbf{r},t)& =\frac{1}{8}\sum\limits _{j=1}^{4}{\int }_{\mathbf{k}}{T}_{0}(\mathbf{r}+{\mathbf{v}}_{g}^{j}t)\mathrm{d}\mathbf{k},\hfill \\ \hfill & {\int }_{\mathbf{k}}\dots \mathrm{d}\mathbf{k}\enspace \stackrel{\;\text{def}}{=}\enspace {\int }_{0}^{1}{\int }_{0}^{1}\dots \mathrm{d}{k}_{1}\enspace \mathrm{d}{k}_{2},\hfill \end{align} \tag{ 14 }$$

where k₁, k₂ are defined by formula (12). Detailed derivation of formula (14) is given in paper [19].

In formula (14), the first term, T_F, describes high-frequency oscillations of the temperature caused by local transition to thermal equilibrium at short times. This process in the graphene lattice is considered in detail in paper [27]. In particular, it is shown that T_F practically vanishes at times of order of ten periods of atomic vibrations. Here we focus on behavior of the second term, T_S, describing slow changes in the temperature profile due to ballistic heat transport. Characteristic time scale of the heat transfer is much larger than the time scale of the transition to thermal equilibrium (see section 6.1 for details). Therefore in further analysis T_F is neglected.

In the following section, formula (14) is employed for description of evolution of sinusoidal and circular initial temperature profiles.

We consider sinusoidal initial temperature profile in the direction given by unit vector e (e.g. for zigzag e = i, while for armchair e = j):

$$\begin{equation}{T}_{0}(\mathbf{r})={T}_{\mathrm{b}}+{\Delta}T\enspace \mathrm{sin}\frac{2\pi \mathbf{r}\cdot \mathbf{e}}{L},\end{equation} \tag{ 15 }$$

where T_b, ΔT are constants such that T_b > ΔT; L is the wavelength of sine. In this case, the general solution (14) is simplified. Additionally, we obtain an approximate solution in the closed form.

Substituting formula (15) into (14) after some transformations yields

$$\begin{align}\hfill T& ={T}_{\mathrm{F}}+{T}_{\mathrm{S}},\hfill \\ \hfill & {T}_{\mathrm{F}}=\frac{1}{8}\left({T}_{\mathrm{b}}+{\Delta}T\enspace \mathrm{sin}\frac{2\pi \mathbf{r}\cdot \mathbf{e}}{L}\right)\sum\limits _{j=1}^{4}{\int }_{\mathbf{k}}\mathrm{cos}(2{\omega }_{j}t)\mathrm{d}\mathbf{k},\hfill \\ \hfill {T}_{\mathrm{S}}& =\frac{{T}_{\mathrm{b}}}{2}+\frac{{\Delta}T}{8}\mathrm{sin}\frac{2\pi \mathbf{r}\cdot \mathbf{e}}{L}\sum\limits _{j=1}^{4}{\int }_{\mathbf{k}}\mathrm{cos}\frac{2\pi \mathbf{e}\cdot {\mathbf{v}}_{g}^{j}t}{L}\mathrm{d}\mathbf{k}.\hfill \end{align} \tag{ 16 }$$

Formula (16) shows that the temperature profile remains sinusoidal at any moment in time. Therefore instead of calculating the entire temperature field, we only compute amplitude, A, of the sine defined as

$$\begin{equation}A(t)=\frac{2}{L}{\int }_{0}^{L}T(\mathbf{r},t)\mathrm{sin}\frac{2\pi z}{L}\mathrm{d}z,\quad z=\mathbf{r}\cdot \mathbf{e}.\end{equation} \tag{ 17 }$$

Substitution of (16) into (17) yields

$$\begin{equation}A(t)=\frac{{\Delta}T}{8}\sum\limits _{j=1}^{4}{\int }_{\mathbf{k}}\left(\mathrm{cos}(2{\omega }_{j}t)+\mathrm{cos}\frac{2\pi \mathbf{e}\cdot {\mathbf{v}}_{g}^{j}t}{L}\right)\mathrm{d}\mathbf{k}.\end{equation} \tag{ 18 }$$

According to formula (16), T_F and T_S have different time scales, proportional to 1/ω_* and L/c_* respectively, where ${\omega }_{\ast }=\sqrt{c/M}$, c_* = ω_* a. The first time scale is determined by frequencies of vibrations of individual atoms. The second time scale is determined by a time required for a sound wave to pass distance L. The ratio of these time scales, L/a, is usually a large parameter. Therefore the time scales are well separated. In further calculations, we focus on behavior of sine at large times of order of L/c_*. Neglecting the first term in brackets in formula (18), we obtain

$$\begin{equation}A(t)\approx \frac{{\Delta}T}{8}\sum\limits _{j=1}^{4}{\int }_{\mathbf{k}}\mathrm{cos}\frac{2\pi \mathbf{e}\cdot {\mathbf{v}}_{g}^{j}t}{L}\mathrm{d}\mathbf{k}.\end{equation} \tag{ 19 }$$

Formula (18) shows that time evolution of amplitude A depends on direction, e, of the initial thermal perturbation. Therefore, in general, the ballistic heat transport in graphene is anisotropic, although graphene's in-plane elastic properties are isotropic for small strains [36]. The accuracy of formula (19) is examined below.

How strong is the anisotropy of ballistic heat transfer in graphene? To address this question, we derive an approximate expression for A(t) at short times.

Consider behavior of A(t) at 1/ω_* ≪ t ≪ L/c_*. Series expansion of the cosine in (19) then yields

$$\begin{equation}A(t)\approx \frac{{\Delta}T}{2}\left(1-\beta {t}^{2}\right),\quad \beta =\frac{{\pi }^{2}}{2{L}^{2}}\sum\limits _{j=1}^{4}\mathbf{e}\cdot {\int }_{\mathbf{k}}{\mathbf{v}}_{g}^{j}{\mathbf{v}}_{g}^{j}\mathrm{d}\mathbf{k}\cdot \mathbf{e}.\end{equation} \tag{ 20 }$$

It is seen that the dependence of A on the direction e is determined by the second rank tensors ${\int }_{\mathbf{k}}{\mathbf{v}}_{g}^{j}{\mathbf{v}}_{g}^{j}\mathrm{d}\mathbf{k}$. Since the graphene lattice has the third order symmetry then all these tensors are spherical and then

$$\begin{equation}\beta =\frac{{\pi }^{2}}{4{L}^{2}}\sum\limits _{j=1}^{4}{\int }_{\mathbf{k}}{{\mathbf{v}}_{g}^{j}}^{2}\mathrm{d}\mathbf{k}.\end{equation} \tag{ 21 }$$

Numerical evaluation of the integrals in (21) yields

$$\begin{equation}\beta {L}^{2}/{c}_{\ast }^{2}\approx 1.88,\end{equation} \tag{ 22 }$$

where c_* = ω_* a, and ${\omega }_{\ast }=\sqrt{c/M}$.

The complete solution (19) and approximate short-time solution (20) are shown in figure 2. It is seen that at times t ∼ L/c_* sinusoidal profiles in zigzag (e = i) and armchair (e = j) directions decay almost identically, as predicted by formula (20).

Zoom In Zoom Out Reset image size

Thus for t < L/c_* the decay of amplitude, A, of the sinusoidal temperature profile is described by the parabola (20). For symmetry reasons, the coefficient of the parabola is independent on the direction, e, of the initial temperature perturbation. Therefore at short times the anisotropy of heat transfer is negligible. The behavior of the amplitude at larger times is considered below.

In this subsection, we compare analytical predictions with results of numerical solution of equations of motion (3) and (4) with initial conditions (6).

Sinusoidal temperature profiles in zigzag (x = r ⋅ i) and armchair (y = r ⋅ j) directions are considered, i.e.

$$\begin{equation}{T}_{0}={T}_{b}+{\Delta}T\enspace \mathrm{sin}\frac{2\pi x}{L}\quad \mathrm{o}\mathrm{r}\quad {T}_{0}={T}_{b}+{\Delta}T\enspace \mathrm{sin}\frac{2\pi y}{L}.\end{equation} \tag{ 23 }$$

Analytical solutions of heat transfer problems with initial condition (23) are given by formula (16) with e = i and e = j. Since the profile remains sinusoidal then we only compute its amplitude A. Analytical expression for A is given by (19). Integral in this formula is evaluated numerically. The first term in this integral, describing high frequency oscillations of temperature, is neglected.

When solving the equations of motion (3) and (4) numerically, we obtain the temperature field T(r, t) and calculate the amplitude, A, as

$$\begin{equation}A(t)=\frac{2}{{L}^{2}}{\int }_{0}^{L}{\int }_{0}^{L}T(\mathbf{r},t)\mathrm{sin}\frac{2\pi x}{L}\mathrm{d}x\enspace \mathrm{d}y\end{equation} \tag{ 24 }$$

for the temperature profile in zigzag (x) direction (for the armchair direction the formula is similar). Note that in contrast to (17), formula (24) involves additional averaging in y direction. The averaging is added to reduce the influence of randomness of the initial conditions. The temperature field in formula (24) is calculated by averaging kinetic energies of particles with respect to large number $\left(1{0}^{4}\right)$ of realizations. In simulations ΔT/T_b = 0.9901, L = 182a. We note that since the lattice is harmonic then neither absolute values of ΔT, T_b nor their ratio influence the results shown below.

Comparison of predictions of formula (19) with results of numerical solution of equations of motion are shown in figure 3. It is seen that analytical and numerical results are in a good quantitative agreement. Figure 3 demonstrates several specific features of ballistic heat transport in graphene. Firstly, the curves corresponding to zigzag and armchair directions practically coincide for t < L/c_*. Therefore at short times the heat transfer is nearly isotropic, as predicted by formula (20). For t ≫ L/c_*, the anisotropy is significant. The anisotropy is caused by the fact that propagation of short waves in graphene is directional dependent. Secondly, the decay of amplitude at large times is non-monotonic. This is due to ballistic (wave) nature of heat transfer. Thirdly, in both directions the amplitude decays inversely proportional to time (see the subplot in figure 3), while in systems with the Fourier heat transport the decay is exponential.

Zoom In Zoom Out Reset image size

Note that, strictly speaking, the rate of decay of the amplitude should be estimated using asymptotic methods, e.g. the stationary phase method [29]. However, in two dimensions application of this method is not straightforward and requires a separate study, which is beyond the scope of the present paper. Therefore in this paper, we limit our-self to less rigorous numerical arguments. The inset plot in figure 3 shows the oscillations of the function Ac_* t/(ΔTL) over the time scale corresponding to the main plot. The amplitude of these oscillations neither decrease nor increase in time. We have checked numerically that similar behavior is also observed at large times (at least up to c_* t/L = 150). Therefore we conclude that the amplitude decays approximately as 1/t in the most interesting time interval in which the changes in temperature are most noticable.

In this subsection, we compare results of numerical solution of equations of motion with the analytical solution (14).

In numerical simulations, we consider two square samples of length L = 182a and L = 1818a. We compute the temperature profiles in these samples at t = 15.91τ_* and t = 176.8τ_* respectively. These moments in time are chosen such that the temperature front reaches the boundary of the corresponding sample.

To construct the analytical solution, we substitute the initial temperature distribution (25) into the general solution (14). Integral in this formula is evaluated numerically. Temperature profiles at t = 15.91τ_* and t = 176.8τ_* calculated using formula (14) are shown in figure 4.

Zoom In Zoom Out Reset image size

To compare analytical solution (14) with numerical simulation results, we consider convergence of the temperature field with respect to number of realizations. In each realization, kinetic energies of all atoms are calculated. To compute the temperature field, the kinetic energies are averaged over 10, 10², 10⁴ and 10⁵ realizations (see figure 5). Figure 5 shows that the temperature field converges with increasing number of realizations. The convergence is practically achieved for 10⁴ realizations.

Zoom In Zoom Out Reset image size

We also check that numerical results converge to the analytical solution (14). For this purpose we plot the distribution of temperature T(x, 0) along the x axis (see figure 6) at t = 15.91τ_* and at t = 176.8τ_*. In this calculation, atoms that are at a distance less than or equal to 0.5a from the x axis are taken into account. If two atoms belong to the same unit cell then the average of their temperatures is plotted. To compute the temperatures, kinetic energies of the atoms are averaged over 10⁵ realizations for a small sample and over 10⁴ realizations for a large one. Figure 6 shows that the analytical solution (14) describes results of numerical simulations with acceptable accuracy.

Zoom In Zoom Out Reset image size

In this subsection, we discuss several peculiarities of the solution, shown in figures 5 and 6.

To start, we consider similar problem in a material with the Fourier heat transfer. In this case, the temperature profile has Gaussian-like shape, i.e. it has one local maximum at the center and it monotonically decrease with increasing distance from the center. In other words, the hottest point is always at the center of the sample.

In the case of ballistic heat transfer, the shape of the temperature field is qualitatively different. Firstly, the temperature profile has well-defined circular front, moving with the maximum group velocity. Secondly, in addition to the local maximum at the center, there are multiple local maxima, moving with different speeds. Thirdly, the position of the global maximum of the temperature changes in time, i.e. the hottest point is not always at the center. For example, at t = 15.91τ_* the hottest point is at the center, while at t = 176.8τ_* there are six identical hottest points located at distance about 600a from the center (see figures 5 and 6). At larger times, the hottest points move along the main lattice directions and form a regular hexagon.

To explain these peculiarities, we analyze the analytical solution (14). The solution may be interpreted in terms of the kinetic theory of heat transfer ⁶ . In the kinetic theory, the temperature is carried out by quasi-particles moving with group velocities. These quasi-particles are sometimes erroneously associated with phonons. We believe that they should be associated with the wave-packets (see discussion in paper [37]). Then the kinetic temperature at some point is proportional to density of quasi-particles at this point. Formula (14) shows that there are four types of quasi-particles, corresponding to four branches of the dispersion relation. Initially, the quasi-particles are uniformly distributed in a circle. Since the quasi-particles move freely then their motion is completely determined by the distribution of group velocities.

We compute the distribution of group velocities as follows. We calculate absolute values of the group velocities $\vert {\mathbf{v}}_{g}^{j}({k}_{1}^{s},{k}_{2}^{p})\vert$ at vertices of the square mesh ${k}_{1}^{s}=s{\Delta}k$, ${k}_{2}^{p}=p{\Delta}k$, s, p = 0, ..., N_k − 1, where Δk = 1/N_k . The values of group velocities form four arrays of length ${N}_{k}^{2}$. These arrays are sorted in ascending order and maximum group velocities $\mathrm{max}\vert {\mathbf{v}}_{g}^{j}\vert$ are calculated. Each interval from 0 to $\mathrm{max}\vert {\mathbf{v}}_{g}^{j}\vert$ is divided into equal segments of length Δv. Number of array elements in all these intervals is calculated. Then the distribution function is calculated as

$$\begin{equation}{\phi }_{j}(v)=\frac{{n}_{j}(v-{\Delta}v/2,v+{\Delta}v/2)}{{\Delta}v{N}_{k}^{2}},\end{equation} \tag{ 26 }$$

where n_j (v − Δv/2, v + Δv/2) is a number of elements from array j in the interval from v − Δv/2 to v + Δv/2. Function ϕ_j (v) shows the relative number of quasi-particles with velocities close to v. Further we assume that each quasi-particle carries the same amount of thermal energy [37] (see formula (14)).

The functions (26) are shown in figure 7. It can be noted that the maximum group velocities $\mathrm{max}\vert {\mathbf{v}}_{g}^{j}\vert ,j=1,2,3,4$ are equal to 0.583c_*, 0.795c_*, 0.535c_*, and 0.277c_* respectively. It is seen that the distributions are non-uniform. In particular, three out of four distributions have multiple sharp local maxima. The most pronounced maxima are at points 0.068c_*, 0.277c_*, 0.424c_*, 0.517c_*, 0.525c_*, 0.795c_*. Since the temperature is proportional to the density of quasi-particles then the presence of local maxima on distributions causes the maxima of the temperature field shown in figures 4–6. In particular, the highest maxima in figure 6(b) (at x/a ∼ 580 and x/a ∼ 620) are formed by the quasi-particles moving with velocities close to 0.517c_* and 0.525c_*. Since the difference between the velocities is small then these maxima merge at short times and form a single maxima at x/a ∼ 50 (see figure 6(a)).

Zoom In Zoom Out Reset image size

Thus interpretation of the analytical solution (14) in terms of the kinetic theory yields simple explanation of the presence of multiple local maxima of the temperature field (hot points) shown in figures 4–6. The hot points are formed by the quasi-particles moving at speeds corresponding to maxima of the distribution function ϕ_j , defined by formula (26).

The theory suggests that in heat conducting harmonic crystals the kinetic temperatures, corresponding to different degrees of freedom of the unit cell, are generally different (see paper [19] for discussion). In particular, in papers [19, 38] this phenomenon have been observed in the harmonic one-dimensional chain with alternating masses. In this subsection, we show that this phenomenon is also present in the graphene lattice.

The presence of several distinct temperatures is demonstrated by numerical solution of the problem with circular initial temperature profile (25). Simulation results show that during heat transfer each atom has two distinct temperatures, corresponding to motions in zigzag (x) and armchair (y) directions even though initially these temperatures are equal. For each unit cell, the temperatures (9) satisfy the relations (except for a small number of cells, where all four temperatures are equal):

$$\begin{equation}{T}_{11}\approx {T}_{33},\qquad {T}_{22}\approx {T}_{44},\qquad {T}_{11}\ne {T}_{22}.\end{equation} \tag{ 27 }$$

To demonstrate this phenomenon, we plot temperatures (9) of individual atoms at t = 15.91τ_*, obtained in computer simulation (see figures 8(a) and (b)). Figure 8(a) shows temperatures T₁₁/T₁ and T₃₃/T₁ corresponding to motion along the x axis, while figure 8(b) shows temperatures T₂₂/T₁ and T₄₄/T₁ corresponding to motion along the y axis. It is seen that kinetic energies of thermal motion (and corresponding kinetic temperatures) are significantly different.

Zoom In Zoom Out Reset image size