We have studied the phonon specific heat in single-layer, bilayer, and twisted bilayer graphene. The calculations were performed using the Born-von Karman model of lattice dynamics for intralayer atomic interactions and spherically symmetric interatomic potential for interlayer interactions. We found that at temperature *T* < 15 K, specific heat varies with temperature as *T ^{n}*, where

*n*= 1 for graphene,

*n*= 1.6 for bilayer graphene, and

*n*= 1.3 for the twisted bilayer graphene. The phonon specific heat reveals an intriguing dependence on the twist angle in bilayer graphene, which is particularly pronounced at low temperature. The results suggest a possibility of phonon engineering of thermal properties of layered materials by twisting the atomic planes.

Specific heat, *C*, is one of the key parameters that characterize the phonon and thermal properties of materials. It is defined as $C=\delta Q/\delta T$, where $\delta Q$ is the change in energy density of a material when temperature changes by $\delta T$.^{1} In the Debye model, the phonon specific heat at low temperature varies with temperature *T* as *T ^{3}* in three-dimensional (3D) case,

*T*in two-dimensional (2D) case and

^{2}*T*in one-dimensional case, provided that the linear dependence of the phonon frequency $\omega $ on the phonon wavenumber

*q*is assumed.

^{1}Graphene, as a strictly 2D system,

^{2,3}provides a unique possibility for investigation of 2D electron and phonon gases.

^{2–6}Due to its high electron mobility and high thermal conductivity, graphene and few-layer graphene (FLG) are considered promising for applications in field-effect transistors,

^{7}interconnects,

^{8,9}thermal interface materials,

^{10,11}and heat spreaders.

^{12–14}

Phonon energy spectra and thermal conductivity in graphene have been investigated both theoretically and experimentally. Experimental studies of optical phonon spectra were performed with Raman spectroscopy^{15,16} while the thermal conductivity measurements were carried out using optothermal.^{4–6} or electrical self-heating methods.^{17,18} Theoretical investigations employed different models of the lattice vibrations, molecular dynamics simulations, density functional method and Boltzmann transport equation approach.^{13,14,19–23} Recently the attention has been shifted to twisted few-layer graphene (T-FLG), which demonstrates intriguing electron^{24–27} and phonon properties.^{28–36} Several independent experimental Raman studies of twisted graphene found series of phonon peaks in different energy regions, which are absent in FLG without the twist.^{28–30,32–37} These experimental results motivate theoretical research to describe phonons in T-FLG.

While the electronic properties of T-FLG were extensively investigated theoretically,^{26,27} the theoretical studies of phonons in T-FLG are rare.^{31,32} In one reported study, the Brenner potential for the intralayer interaction and Lennard-Jones potential for the interlayer interactions have been used to investigate phonon modes in a circular twisted bilayer graphene (T-BLG) with a finite radius *R* = 0.5 nm – 3 nm.^{31} It was found that the rotational angle strongly influences ZO phonon modes. We have theoretically investigated the hybrid folded phonons in twisted bilayer graphene using the Born-von Karman model for the intralayer interactions and Lennard-Jones potential for the interlayer interactions.^{32} We have shown that the twist angle affects the low-energy phonons the most while the phonon density of states (PDOS) and phonon average group velocity depend only weakly on the rotational angle.^{32} While the thirst theoretical studies of the phonon spectra in T-FLG appeared, there have been no reported investigations of thermal properties such as specific heat or thermal conductivity. In this Letter, we report the results of the theoretical study of specific heat of T-BLG, which reveal its intriguing dependence on the twist angle.

We consider a bilayer graphene structure with a relative rotational angle $\theta $ between the parallel graphene sheets, and define the initial stacking configuration with $\theta =0\xb0$ as AB stacking. To obtain T-BLG, we rotate one atomic layer of carbon atoms relative to another layer by an angle $\theta $ as shown in Figure 1(a). The chosen rotation scheme limits our consideration to rotational angles between $\theta =0\xb0$ and $\theta =60\xb0$, i.e., between AB and AA stacking configuration. The commensurate structures, i.e., structures with translation symmetry, exist for certain rotational angles only determined by the following condition:^{38} $cos\u2009\theta (p,n)=$ $(3p2+3pn+n2/2)/(3p2+3pn+n2)$, where *p* and *n* are co-prime positive integer numbers. If *n* is not divisible by 3, the number of atoms in T-BLG unit cell is given by^{32}

The unit cells of T-BLG with larger indices (*p,n*) contain a larger number of carbon atoms. For instance, the unit cell of T-BLG with $\theta (1,1)=21.8\xb0$ contains the smallest possible number of atoms *N* = 28 while a rotation by $\theta (2,1)=13.2\xb0$ increases this number to *N* = 76. In order to construct the Brillouin zone (BZ) of the T-BLG with the rotation angle $\theta (p,n)$, one should determine the corresponding reciprocal vectors $g\u21921$ and $g\u21922$ of T-BLG.^{32} The size of BZ in T-BLG depends on the rotation angle. In Figure 1(b), we show BZ of the single layer graphene (blue curves) and T-BLG with $\theta (1,1)=21.8\xb0$ (red curve), determined using a general method described by us in Ref. 32.

In the harmonic approximation, the equation of motion of an atom in T-BLG is given by^{32}

Here, $\omega $ is the phonon frequency, $q\u2192=(qx,qy)$ is the 2D phonon wave vector, *s* denotes the nearest atomic spheres of the atom *i; j* denotes the atoms of the atomic sphere *s*; *k* and $k\u2032$ designate the Cartesian coordinate components, $m=19.9442\xd710\u221227$ kg is the mass of a carbon atom, $\Phi kk\u2032ij(s)$ are the interatomic force constants describing the interaction between the atom *j* and the atom *i*, located at the center of the atomic spheres, and $r\u2192ij=r\u2192j\u2212r\u2192i$, where $r\u2192i$[ $r\u2192j$] is the radius vector of the atom *i*[*j*]. The phonon energies and thermal conductivity in graphene depend strongly on the interatomic potentials and their parameters.^{19,32,39} Therefore, a proper choice of the force constant matrices $\Phi $ is crucial for obtaining the correct phonon frequencies and for accurate calculation of phonon specific heat. We describe the strongly anisotropic interatomic interactions in T-BLG using the Born–von Karman (BvK) model of lattice dynamics for intralayer interaction and a spherical potential for interlayer interactions.

The general form of a BvK force constant matrices describing the interaction of an atom with its *s*th-nearest neighbor in a graphene sheet is^{32,40,41}

In our model, we take into account 4 nearest atomic spheres and 12 independent force constants,^{42,43} which were fitted to the experimental in-plane phonon dispersion of bulk graphite.^{44} In case of the interlayer coupling, the atomic configuration and force constant matrices are directly dependent on the rotational angle $\theta $. Therefore, an interatomic potential, which explicitly takes into account the dependence of the interaction strength on the atomic positions, should be used. One possible candidate is a well-known Lennard-Jones potential,^{45} which has been used to model interlayer interaction in graphite, FLG, and T-BLG.^{21,31,32,46} However, it was pointed out^{32,46} that the Lennard-Jones potential underestimates the frequencies of the shear vibrations in graphite.

Here, we describe the interlayer interaction by a spherically symmetric interatomic potential with the following components of corresponding force constant matrices (FCM):

where $\delta $ is the force constant of the interlayer coupling, which depends only on the distance between interacting atoms. Since the interlayer coupling in graphite is very weak and the interaction strength between atoms strongly decreases with $rij$, we model the dependence of the force constant $\delta $ on $rij$ as a decaying exponent: $\delta (rij)=A\u2009exp(\u2212rij/B)$ with the two fitting parameters *A* = 573.76 N/m and *B* = 0.05 nm determined from comparison with the experimental phonon frequencies from Г-A direction in graphite.^{47} The cutoff radius of 1 nm was used for $\delta (rij)$ to simplify summation over *j* in Eq. (2). In Figure 2(a) we show the results of our calculations of the phonon energy dispersion in bulk graphite along Г-A direction calculated using FCMs from Eq. (4) (red curves) and Lennard-Jones potential (black curves) with parameters $\epsilon $ = 4.6 meV, $\sigma $ = 0.3276 nm, taken from Ref. 21. The blue triangles represent the experimental points from Ref. 47. The curves, calculated with FCMs, are in excellent agreement with the experimental data while the frequencies of the doubly degenerate lowest phonon branches, calculated with the Lennard-Jones potential, are underestimated. The presented data provide experimental validation for our theoretical approach.

In Figures 2(b) and 2(c), we show the phonon energy spectra in AB-stacked bilayer graphene (panel b) and twisted bilayer graphene (panel c) with $\theta =21.8\xb0$. The experimental phonon energies of graphite, reproduced from Ref. 44, are indicated by the blue triangles. In the single layer graphene (SLG), there are six phonon branches: in-plane transverse acoustic (TA) and optic (TO), in-plane longitudinal acoustic (LA) and optic (LO) and out-of-plane acoustic (ZA) and optic (ZO).^{19,20} Due to the weak van der Waals interactions in bilayer graphene (BLG), these branches split to the twelve nearly-degenerated branches TA_{1}/TA_{2}, TO_{1}/TO_{2}, LA_{1}/LA_{2}, LO_{1}/LO_{2}, ZA_{1}/ZA_{2}, ZO_{1}/ZO_{2}^{6,18,32} with energies close to those in SLG. From Figure 2(b), we can conclude that our model provides an excellent agreement between the theoretical and experimental phonon frequencies. An evolution of phonon energies in T-BLG with increasing the angle between graphene layers were recently investigated by us in Ref. 32. The augmentation of the twist angle leads to a decrease of the Brillouin zone size and increase in the number of phonon branches. The hybrid folded phonons appear in T-BLG due to mixing of phonons from different directions of BZ in regular BLG (see Figure 2(c)). The energies of these phonons depend on the twist angle. The rotation angle-dependent Raman peaks, associated with hybrid folded phonons, have been observed experimentally,^{28–30,32–37} confirming the theoretical predictions.

For calculation of the phonon specific heat in T-BLG, we use the following formula:^{1,48}

where $\omega $ is the phonon frequency, $\omega max$ is the maximum phonon frequency, *f* is the 2D normalized phonon DOS, *T* is the temperature, *N _{A}* is the Avogadro constant,

*k*is the Boltzmann constant, and $\u210f$ is the Planck constant. The normalized phonon DOS is given by

_{B}where $g(\omega )$ is the 2D phonon DOS given by the relation

Here, *s*′ numerates phonon branches. In order to calculate $g(\omega )$, we applied a 200 × 200 2D grid to a 1/4th part of BZ of T-BLG (shown as a green segment in Figure 1(b)), and then calculated phonon frequencies for every (*q _{x}*,

*q*) point in this grid.

_{y}The dependence of specific heat on temperature for SLG and AB-BLG is presented in Figure 3(a). The experimental values of graphite heat capacity reported in Ref. 49 are also shown for comparison by the blue triangles. The low-temperature heat capacity of SLG is higher than that in graphite due to inequality of 2D and 3D phonon DOS. The difference in the heat capacity of SLG, BLG, and graphite diminishes with temperature rise. The heat capacities become identical within 0.01% deviation for *T* > 2000 K. At temperatures *T* > 2500 K, all heat capacities approach the classical Dulong-Petit limit $cV$ = 24.94 J K^{−1 }mol^{−1}. At small frequencies, ZA phonons demonstrate a quadratic dispersion $\omega \u223cq2$, leading to $cV(ZA)\u223cT$, while TA and LA phonons possess linear dispersions $\omega \u223cq$, resulting in $cV(LA,TA)\u223cT2$. The total heat capacity reveal a complicated temperature dependence $Tn$, where $1\u2264n\u22642$ resulting from the specifics of the phonon dispersion and population of *ZA*, *LA*, and *TA* modes. Our results show (see Fig. 3(a)) that in SLG, *n* = 1 for *T* $\u2264$ 15 K; *n* = 1.1 for 15 K < *T* $\u2264$ 35 K; *n* = 1.3 for 35 K < *T* $\u2264$ 70 K, and *n* = 1.5 for 70 K < *T* $\u2264$ 240 K. The power factor *n* increases with temperature not only due to a greater contribution of *LA* and *TA* phonons but also due to the deviation of *ZA* phonon frequency dispersion from the quadratic law, leading to both nearly linear dependence of *ZA* phonon frequencies on *q* in the interval 6 nm^{−1 }< *q* *<** *12 nm^{−1} (in $\Gamma \u2013K$ direction) and almost linear dependence of *ZA* PDOS on frequency for $\omega $ > 100 cm^{−1}. Our results for heat capacity in SLG are qualitatively similar to those described in Refs. 50 and 51 with the exception of the power factor values *n* = 1 (Ref. 50) and *n* = 1.1 (Ref. 51) for *T* < 100 K. The discrepancy in exact numerical values is explained by the fact that in Refs. 50 and 51 the authors did not take into account the non-parabolicity of *ZA* dispersion for $\omega $ > 100 cm^{−1}. In AA-BLG or AB-BLG, $cV\u223cT1.6$ for wide temperature range *T* $\u2264$ 170 K due to the changes in low-frequency PDOS in comparison with those in SLG.

In Figure 3(b), we plot a difference between the specific heat in AB-BLG and T-BLG as a function of temperature: $\Delta cv(\theta )=cv(AB)\u2212cv(\theta )$ for $\theta =21.8\xb0$, $\theta =13.2\xb0$ and $\theta =9.4\xb0$. The change in the specific heat due to twisting is relatively weak in a wide temperature range of 20 K–2000 K. It attains its maximum value ∼0.028 J K^{−1 }mol^{−1} at *T* ∼ 250 K. At the same time, at low temperatures the relative difference between specific heat in AB-BLG and T-BLG $\eta =(1\u2212cV(\theta )/cV(AB))\xd7100%$ constitutes substantial 10%–15% at *T* = 1 K and ∼3%–6% at *T* = 5 K in dependence on $\theta $ (see blue, red and green curves from the inset to Figure 3(b)). The low temperatures specific heat depends stronger on the twist angle because twisting affects the low-frequency ZA phonon modes the most.^{32} A somewhat similar effect of modulation of the specific heat was recently reported for SLG under strain.^{52,53} The authors found that a stronger modification of the specific heat occurs at temperatures below 1 K. The temperature dependence of low-temperature specific heat in T-BLG with $\theta =21.8\xb0$ differs from SLG or BLG: $cV\u223cT1.3$ for T < 10 K and $\u223cT1.6$ for 10 K $\u2264$ *T* $\u2264$ 100 K.

In summary, we found that the phonon specific heat reveals an intriguing dependence on the twist angle in T-BLG, which is particularly pronounced at low temperature. The results suggest a possibility of phonon engineering of thermal properties of layered materials by twisting the atomic planes. Specifically, the twisting decreased the specific heat in T-BLG by 10%–15% in comparison with reference BLG at *T* ∼ 1 K. It changes the temperature dependence of the heat capacity from $cV(BLG)\u223cT1.6$ to $cV(TBLG)\u223cT1.3$ for *T* < 10 K. The maximum deviation of the specific heat in TBG from that in AB-bilayer graphene was found at *T* ∼ 250 K.

This work was supported by the National Science Foundation (NSF) Grant No. ECCS-1307671 Two Dimensional Performance with Three Dimensional Capacity: Engineering the Thermal Properties of Graphene. D.L.N. and A.I.C. acknowledges the financial support through the Moldova State projects 11.817.05.10F, 14.819.16.02F, and Moldova-STCU project 14.820.18.02.012 STCU.A/5937.