Anharmonic properties from a generalized third order ab initio approach: theory and applications to graphite and graphene

# Anharmonic properties from a generalized third order ab initio approach: theory and applications to graphite and graphene

Lorenzo Paulatto    Francesco Mauri    Michele Lazzeri IMPMC, Université Pierre et Marie Curie, CNRS, 4 place Jussieu, F-75005 Paris, France
July 19, 2019
###### Abstract

We have implemented a generic method, based on the 2n+1 theorem within density functional perturbation theory, to calculate the anharmonic scattering coefficients among three phonons with arbitrary wavevectors. The method is used to study the phonon broadening in graphite and graphene mono- and bi-layer. The broadening of the high-energy optical branches is highly nonuniform and presents a series of sudden steps and spikes. At finite temperature, the two linearly dispersive acoustic branches TA and LA of graphene have nonzero broadening for small wavevectors. The broadening in graphite and bi-layer graphene is, overall, very similar to the graphene one, the most remarkable feature being the broadening of the quasi acoustical ZO’ branch. Finally, we study the intrinsic anharmonic contribution to the thermal conductivity of the three systems, within the single mode relaxation time approximation. We find the conductance to be in good agreement with experimental data for the out-of-plane direction but to underestimate it by a factor 2 in-plane.

###### pacs:
63.20.kg,63.20.dk,63.22.Rc,65.80.Ck

## I Introduction

Thermal transport is currently attracting much attention; the main applications of interest are materials for thermoelectric energy conversion Mahan et al. (1997); *disalvo99 and materials used for thermal dissipation in microelectronics Cahill et al. (2003). While, in the first case, the goal is to engineer the smallest possible thermal conduction, in the second case a good thermal conduction is required. In general, our understanding of thermal properties is heavily based on theoretical modeling and the use of precise and reliable approaches, such as the ab initio computational methods, is highly desirable.

The presence of a temperature gradient in a solid induces a heat flux. The heat carriers can be lattice vibrations (phonons) or electronic excitations. In general, lattice conduction is the dominant mechanism in the presence of an electronic gap (semiconductors and insulators) or when the gap is zero but the density of electronic states at the Fermi level is small (semi-metals). Lattice thermal resistance is then dictated by phonon scattering, which can be induced by extrinsic mechanisms (isotopic disorder, structural defects, finite-size of the crystals, etc.) or by intrinsic ones (anharmonic phonon-phonon scattering). Determining the intrinsic anharmonic scattering is in itself a very complex task; it has been attempted ab initio, within density functional theory (DFT), using finite difference derivation Narasimhan and Vanderbilt (1991); Tang et al. (2010); *tang11 or molecular dynamics techniques Donadio and Galli (2007) or from linear response theory.Broido et al. (2007); Ward et al. (2009); *ward10; Garg et al. (2011); Bonini et al. (2012)

In a crystal, the intrinsic lattice thermal conduction can be obtained by knowing harmonic phonon energies and anharmonic phonon-phonon scattering coefficients. Harmonic phonon energies are determined by the second order derivative of the system total energy, with respect to atomic displacements. This second derivative can be efficiently calculated ab initio by using density functional perturbation theory (DFPT) Baroni et al. (2001), which allows the determination of the phonon dynamical matrix for an arbitrary q wavevector in the Brillouin zone. DFPT is implemented in the Quantum ESPRESSO package www (); *qe, within the plane-waves and pseudopotential approaches. The anharmonic scattering coefficients can be determined by the third order derivative of the energy with respect to three phonon perturbations, characterized by the wavevectors q, , . For the thermal transport problem, it is necessary to know these derivatives with respect to three arbitrary wavevectors (possibly , , ), with the only condition , where G is a reciprocal lattice vector. In principle, these coefficients can be obtained within DFPT by using the, so called, “2n+1” theorem as formulated by Ref. Gonze and Vigneron, 1989. This theorem allows us to access the 3rd derivative of the total energy by using only the 1st derivative of ground state density and wavefunctions; contrary to the finite differences approach, we do not have to perform expensive supercell calculations.

The first implementations of the “2n+1” approach were limited to the scattering of one zero-momentum phonon towards two phonons with opposite arbitrary momenta (0, -q, q). Ref. Debernardi et al., 1995 implemented this approach for insulating and semiconducting materials. Later, Ref. Lazzeri and de Gironcoli, 2002 generalized the approach to metals and zero-gap materials. The (0, -q, q) anharmonic coefficients can be computed within Quantum ESPRESSO www (); *qe, using the d3 code which was developed in Ref. Lazzeri and de Gironcoli, 2002. These coefficients can be easily used to compute the anharmonic broadening of a phonon (see e.g. Ref. Lazzeri et al., 2003). By using a super-cell approach one can also compute the broadening of phonons having q commensurate with the super cell. However, the super-cell approach (which was used in Ref. Garg et al., 2011; Bonini et al., 2012; Lazzeri et al., 2003) can be computationally very demanding. Recently, Ref. Deinzer et al., 2003 further extend the method to three arbitrary phonons, (), although only for insulators/semiconductors. This was done within a proprietary non-publicly-available software.

We have developed and implemented an extension of the d3 code of the Quantum ESPRESSO software to compute, within the DFPT “2n+1” approach, the three-phonons anharmonic coefficients for three arbitrary wavevectors for insulators/semiconductors and also for metallic or zero gap systems. The coefficients thus obtained can be used in a straightforward way to compute the anharmonic broadening of a phonon with an arbitrary wavevector q and the intrinsic thermal conductivity within the single-mode relaxation time approximation (SMA).  Srivastava (1974); *asenpalmer97; *cao04; *khitun01 The first applications of the method are devoted to graphite, graphene and graphene bi-layer. Indeed, the thermal properties of these systems have attracted significant attention, Balandin (2011); Seol et al. (2010) being carbon systems excellent thermal conductors.

Sec. II describes the method. Sec. III reports the results and the discussion. Conclusions are summarized in Sec. IV.

## Ii Method

In Sec. II.1 we provide the expressions for the phonon anharmonic broadening and for the phonon thermal conduction. In Sec. II.2 the method is briefly described and we provide the relevant computational details. A more detailed description of the method is reported in the Appendix.

### ii.1 Anharmonic decay and thermal transport

Let us consider the total energy for a crystal , where is the displacement from the equilibrium position of the atom along the Cartesian coordinate in the unit cell identified by the lattice vector . We define

 uq,s,α=1N∑Re−iq⋅RvR,s,α, (1)

where the sum is performed on the lattice vectors and is the number of cells involved in the summation. We define the dynamical matrix

 D2(qss′αα′)=1N∂2Etot∂u−q,s,α∂uq,s′,α′, (2)

the angular frequency of a phonon with wavevector q and branch index is obtained by solving

 ∑s′,α′1√msms′D2(qss′αα′)zq,js′,α′=ω2q,jzq,js,α, (3)

where are the orthogonal phonon eigenmodes normalized in the unit cell and is the atom mass. We define the three-phonon scattering coefficients as

 V(3)qj,q′j′,q′′j′′=1N∂3Etot∂Xq,j∂Xq′,j′∂Xq′′,j′′, (4)

where

 ∂∂Xq,j=∑s,α√ℏ2msωq,jzq,js,α∂∂uq,s,α. (5)

has the dimension of an energy and does not depend on , while is adimensional. Because of the translational symmetry of the crystal, the coefficients from Eq. 4 are only when , where G is any reciprocal lattice vector.

With these definitions, the lifetime due to anharmonic phonon–phonon interaction, , and the corresponding broadening (full width at half maximum) of the phonon are Lazzeri et al. (2003):

 1τqj(T)=γqj(T)=πℏ2Nq∑q′,j′,j′′∣∣V(3)qj,q′j′,q′′j′′∣∣2×[(1+nq′j′+nq′′j′′))δ(ωqj−ωq′j′−ωq′′j′′)×[+2(nq′j′−nq′′j′′)δ(ωqj+ωq′j′−ωq′′j′′)]. (6)

Where is the temperature, is the Bose-Einstein statistics occupation of phonon , and is the Dirac distribution. The sum is performed over a sufficiently fine grid of q-points in the Brillouin zone (BZ) and . and depend on only through the phonon occupations .

The r.h.s. of Eq. 6 is usually interpreted as the sum of scattering processes in which a phonon of wavevector q decays into two phonons , , (third line of Eq. 6) or in which the phonon q coalesces with and emits (fourth line of Eq. 6). The energy conservation of the processes are guaranteed by the Dirac delta. One can also distinguish between Normal and Umklapp processes: by choosing q and such that they belong to the first BZ, the scattering is Normal when also belongs to the first BZ; on the contrary, when does not belong to the first BZ, the scattering is Umklapp.

By knowing the anharmonic scattering coefficients, Eq. 4, one can determine the lattice thermal conductivity within the framework of the Boltzmann transport equation (BTE) for phonons Ziman (2001). In general, an exact solution of the BTE is a difficult task; a commonly used approximation to the problem is the so-called single mode relaxation time approximation (SMA)Garg et al. (2011); Srivastava (1974); Asen-Palmer et al. (1997); Cao et al. (2004); Khitun and Wang (2001). Within the SMA, the lattice thermal conductivity tensor becomes:

 κα,βL=ℏ2NqΩKBT2∑qjcαqjcβqjω2qjnqj(nqj+1)τqj. (7)

Here, is the volume of the unit cell, is the Boltzmann constant and is the phonon group velocity of mode along Cartesian direction : . The SMA conductivity from Eq. 7 can be obtained in a straightforward way once the anharmonic lifetimes have been computed from Eq. 6. is a tensor which takes into account the possible anisotropies and transversal conductance. However, in high-symmetry crystals, as graphene and graphite, the off-diagonal elements are zero, if two axes lie in the graphene plane. Moreover, in both graphene and graphite the two in-plane and components are identical. The out-of-plane component is not defined in the bidimensional graphene systems, but it is meaningful in graphite.

The validity limits of the SMA are discussed in literature Garg et al. (2011); Srivastava (1974); Asen-Palmer et al. (1997); Cao et al. (2004); Khitun and Wang (2001). Here, we wish to remind that, for a generic material, the SMA is expected to be valid (that is, to provide the correct solution to the BTE) at room conditions and to break down only at very small temperatures (see, e.g. Ref. Ward et al., 2009; Ward and Broido, 2010).

### ii.2 Density-functional theory calculation

Calculations of the phonon properties are done within density functional perturbation theory Baroni et al. (2001) as implemented in Ref. Giannozzi et al., 2009. The third order coefficients defined in Eq. 4 are computed using a code which has been developed for the present work. This code has been written on the top of a previous less general implementation available within the Quantum ESPRESSO package: the d3 code, which was implemented in Ref. Lazzeri and de Gironcoli, 2002. The method is described in detail in Appendix A, and  B.

We use local-density approximation and the carbon atom is described by a norm-conserving pseudopotential which includes four electrons in valence. Plane waves kinetic energy cutoff is  Ry. For all the systems, the in-plane lattice parameter is  Å, which is the theoretical equilibrium value for graphite. For graphite, we use . This value, which is only slightly different from the experimental value , is chosen phenomenological to accurately reproduce the low frequency phonon dispersion along the -A direction (see the discussion below). For graphene, the periodic replicas of the planes are spaced along the direction with  Å of vacuum. The two layers of the graphene bilayer are spaced with the inter-planar distance of bulk graphite; periodic images are then spaced with  Å of vacuum.

The computational parameters are listed in Appendix D. We remind here that the electronic integration has to be done with a small value of smearing (and a consequent fine-grained k-point grid) due to the presence of a Kohn anomaly for the highest optical branch near K Piscanec et al. (2004) (usually called TO). The phonon frequencies and the third-order coefficients , used in Eqs. 6 and 7, are calculated in a slightly different way. On one hand, phonon energies are corrected using an ad hoc procedure (based on DFT+GW renormalization of the electron-phonon interaction as in Ref. Lazzeri et al., 2008, see Appendix D). This correction affects only the TO branch, it does not touch the other branches, and it provides better agreement with measurements, Fig. 1. On the other hand, the third-order coefficients are computed within standard (less precise) DFT. The using of these two different procedures for the and calculations is not consistent. However, this should not affect the results in a major way: phonon broadening results from a sum over different processes which are selected by energy conservation enforced by the two Dirac in Eq. 6. The intensity of the processes is then proportional to the square of the coefficients. Consequently, the computational accuracy of and that of the coefficients affect the result in a very different way. An error in the phonon dispersion can affect the lifetime in a not predictable way and, thus, a special care should be taken into finding the best possible description of the phonon dispersion. The same care is not strictly necessary for the third order calculations.

Fig. 1 compares measured with calculated phonon dispersions for graphite. Notice that plain DFT calculations do not provide a satisfactory description of the highest optical TO branch near K, while DFT+GW ones do much better. The lower panel of Fig. 1 shows in detail the low frequency dispersion. This region is characterized by the splitting of the acoustic phonon branches of the two graphene planes in the graphite unit cell. These branches are particularly sensitive to the actual value of . In particular, in that region, by changing the lattice parameters from (which is the value used throughout the paper) to the value of the phonon branches change by almost 14%.

Actual DFT calculations are done on a relatively coarse grid of q wavevectors, described in Appendix D. The dynamical matrices and the third order coefficients, necessary to compute the broadening and the thermal conductivity (Eqs. 6 and  7), are then obtained on a finer grid via the Fourier interpolation technique described in Appendix C. Eqs. 6 and  7 are evaluated by performing the sum over a discrete grid of points and by substituting the with a Gaussian function characterized by an artificial smearing . This approximation is valid as long as is smaller than the thermodynamic fluctuation, which is of order . The grids and the values are specified in Appendix D. Here we just remark that the results shown in Sects. III.1 and  III.2 are obtained using a particularly fine-grained sampling. This is only necessary to produce the very sharp features which are present in the broadening of the higher optical bands or to produce the correct behavior of the broadening of the acoustic branches in the vicinity of . Indeed, a much coarser grid is sufficient for most applications, included those described in Sec. III.3.

## Iii Results and discussion

This section reports and describes the results. Section III.1 analyzes the anharmonic phonon broadening in the graphene monolayer. Special relevance is given to the three acoustic branches which are the most important for the thermal transport. Sec. III.2 analyzes the broadening in graphene bilayer and graphite. Sec. III.3 is dedicated to thermal transport.

Fig. 2 shows the calculated graphene phonon dispersion, the respective anharmonic broadening and the vibrational density of states (VDOS). The branches are labeled in the usual way.Mounet and Marzari (2005) There are three acoustic branches (ZA, TA, LA) and three optical branches (ZO, TO, LO). ZA and ZO correspond to an atomic motion perpendicular to the graphene plane ( direction), all the other branches are polarized parallel to the plane. In the vicinity of , TA and TO are quasi transverse, while LA and LO are quasi longitudinal. In the following, these labels will be used to classify the branches all along the high symmetry lines (as in Fig. 2), although this distinction is not meaningful for an arbitrary wavevector in the Brillouin zone (BZ). Because of symmetry, the modes perpendicular polarized (ZA and ZO) are separated from the others all over the BZ. Moreover, the TA branch is always well separated from the other parallel polarized branches (labeled as ). In Fig. 2, we can, thus, separate the VDOS in three distinct components labeled as Z, TA, and . The two dimensional character of the phonon dispersion is associated with some specific features. The ZA branch is quadratic near and, thus, in the limit the VDOS does not go to zero (Fig. 2). The presence of a local maximum in the phonon dispersion (as the one at 1008 cm for the TA branch near K or the one at 904 cm for the ZO one near ) is associated with a step in the VDOS. The presence of a saddle point in the dispersion (as those at 477 cm, 631 cm, 643 cm, and 1432 cm at the M point) is associated with a sharp peak in the VDOS.

Fig. 3 reports in more detail the calculated anharmonic phonon broadening, along high symmetry lines, and its decomposition into the different allowed decay channels. For symmetry reasons, the -polarized branches can only decay toward one Z and one non-Z phonons.Bonini et al. (2007); Lindsay et al. (2010) The other bands can only decay towards two phonons which are either both or neither -polarized. The two most striking features in Figs. 23 are the small behavior of the acoustic branches and the highly non uniform behavior of the broadening.

The existence of a finite broadening at small for the TA and LA acoustic branches is problematic. Indeed, the concept itself of phonon is meaningful only when , being the broadening (i.e. the inverse of the phonon lifetime). From the present calculations, the condition is satisfied for both the TA and LA branches for , with , being the in-plane lattice spacing. Thus, for , the TA or LA frequency can become smaller that the broadening. In this region the present treatment is, obviously, not valid (see the discussion in Ref. Bonini et al., 2012) and a proper treatment of the phenomenon is beyond the present scope. In practice, however, represents a tiny portion of the Brillouin zone (the corresponding region in Figs. 23 has width of the order of the thickness of the vertical line passing through ). As a consequence, we can assume that the properties obtained as a sum over the Brillouin zone (such as the thermal conductivity of Eq. 7) are not affected by a major error.

Concerning the global appearance of Figs. 23, the many sharp peaks in the broadening can be ascribed to different mechanisms. Those in the high energy part of the spectrum are, in general, associated with peaks in the VDOS: when one or both of the final states (i.e. of the states that meet the energy and momentum conservation requirements in Eq. 6) produce a peak in the VDOS, there the broadening will typically exhibit a peak. For example, the large scattering probability, predicted for the M point on the LO branch ( cm), corresponds to a decay toward a ZO phonon close to ( cm) and a ZA phonon close to M ( cm). As the VDOS, Fig. 2, has maximum in both region, this transition is particularly favored. On the other hand, for q0.66M (along the M direction) or for q0.62K (along K), the LO broadening displays a sudden increase. This is because, for these wavevectors, the LO energy becomes small enough to activate the decay channel towards the ZO branch (At q0.66M and q0.62K the LO phonon decays into a ZO with the same wavevector and a ZO with q).

We remark that the presence of sharp peaks which are essentially determined by energy and momentum conservation in the decay process implies that even a small change in the phonon dispersion used in the calculation could induce significant differences in the calculated broadening.

Concerning the three acoustic branches the broadening near is almost entirely due to Normal scattering, Fig. 4. The peaks which are observed at q0.44 M and q0.41 K for the LA branch, and at q0.65 M and q0.54 K for the TA one, are associated with the activation of Umklapp scattering towards the ZA phonons. To have a more comprehensive view, Fig. 5 reports the broadening in the entire Brillouin zone. The TA and LA branches exhibit a feature-rich behavior in a wide region, far from , which roughly starts at about halfway to the first BZ edge. In this region, the anharmonic decay presents a component of Umklapp processes, which is absent in the vicinity of , where the scattering is almost entirely Normal. On the other hand, the ZA broadening is relatively feature-less and isotropic; it is quadratic in q in the center of the first BZ then it saturates and remains roughly constant.

#### iii.1.1 Phonon mean free path

An alternative way to represent the effect of the phonon broadening is to plot the single-phonon mean free path (MFP):

 λqj=τqj∣∣cqj∣∣, (8)

where is the phonon group velocity and is the phonon lifetime, from equation 6. In figure 6 we have plotted the MFP at  K for the three acoustic branches. The MFP for the TO and LO bands is of order  nm or smaller. We have verified that it does not get substantially higher at lower temperatures, except in the vicinity of where it diverges at  K. At room temperature, the MFP of the ZA bands is one order of magnitude larger, i.e. of order  m, in the center region of the Brillouin zone. Furthermore, it increases as when temperature decreases. Obviously, especially at small temperatures, when the intrinsic anharmonic MFP is too big, other effects (typically the scattering on the borders of the sample) become important and limit the value of the MFP. We remark that the MFP of acoustical phonons is only one order of magnitude smaller than the typical dimensions of high-quality graphene samples. It is, also, definitely larger than the transverse dimension of graphene nano-ribbons. This result suggests that ballistic phonon-driven conductance could play a relevant role in this kind of systems.

#### iii.1.2 Temperature dependence

The intrinsic anharmonic broadening of a specific phonon (qj) has, in general, a typical dependence of the temperature : It is almost constant below a certain characteristic temperature , then it rapidly becomes linear in . Such a behavior is reproduced by Eq. 6. A quadratic dependence on can be observed experimentally only at relatively high and it is due to terms of order higher than those included in Eq. 6.  Balkanski et al. (1983)

From Eq. 6, one can check that

 limT→∞γqj(T)=αqjT+O(1/T), (9)

where

 αqj= πKBℏ3Nq∑q′,j′,j′′∣∣V(3)qj,q′j′,q′′j′′∣∣2 ×[(1ωq′j′+1ωq′′j′′)δ(ωqj−ωq′j′−ωq′′j′′) ×+2(1ωq′j′−1ωq′′j′′)δ(ωqj+ωq′j′−ωq′′j′′)] (10)

does not depend on . One can thus be tempted to approximate the overall dependence on of the broadening by

 γqj(T)≃~γqj(T)=αqjθSqj coth⎛⎝θSqjT⎞⎠, (11)

where is the hyperbolic cotangent function, , and is the broadening from Eq. 6. Indeed, from Eq. 11 is almost constant for and tends to for . Moreover, for . To check the validity of this approximation (Eq. 11), we systematically computed the graphene broadening for different phonons in the temperature range between 0 and 1500 K, using Eq. 6. These results are reasonably well reproduced by Eq. 11 with an error less than 5%.

As a consequence, for a given phonon mode (), the knowledge of the two corresponding parameters and is enough to determine the overall temperature behavior of the broadening, by using Eq. 11. The two parameters can be extracted from Fig. 7, for the graphene acoustic branches, along high symmetry lines.

Finally, the broadening of the LA and TA branches can be fitted with an isotropic function of and of of the form:

 γ(q,T)=qBcoth(qAT). (12)

By defining , where is the cell parameter, we have: for the LA band,  cm and  K; for the TA band,  cm and  K. These fitted parameters reproduce the computed linewidth with an error of generally less than % for (smaller for higher temperature and lower ). For the ZA band, can be assumed, resulting in this simple fitting function:

 γ(q,T)=Bq2T. (13)

Taking  cm we reproduce the value of linewidth to % accuracy in the range , except for systematically underestimating it in the very small region where it is negligible ( cm).

### iii.2 Graphite and bilayer graphene

We now discuss the anharmonic broadening in graphite and graphene bilayer. Each of the six phonon branches of the graphene monolayer splits into two branches for both graphite and graphene bilayer. The three acoustic branches of graphene (ZA, TA, LA) split into three acoustic (ZA, TA, LA) and three quasi-acoustic branches (ZO, TO, LO). The quasi-acoustic branches are almost degenerate with their respective acoustic ones, except in the vicinity of the -A line. The remaining six optical branches are pair-wise quasi-degenerate in the entire Brillouin zone, and they will be referred in pairs simply as ZO, TO and LO or, in some cases as ZO, ZO, etc. This notation does not hold along the line in graphite, as degeneracy changes. It is however still possible to name the branches by continuity below  cm.

Fig. 8 shows a general view of the calculated phonon dispersion and broadening in bulk graphite and graphene bilayer. Fig. 9 and Fig. 10, compare the broadening of the acoustic and quasi-acoustic branches of, respectively, graphite and bilayer graphene with those of the single layer graphene. Fig. 11, shows in more detail the low frequency region. In the high energy part of the spectrum, graphite and graphene monolayer and bilayer are almost indistinguishable, Fig. 8, meaning that, also in graphite, the physics is ruled by the two dimensional character of the phonon dispersion. Some of the graphene sharp features are a slightly broader in graphite due to the out-of-plane phonon dispersion acting like an effective smearing. The most striking differences between three-dimensional bulk graphite and two-dimensional graphene monolayer and bilayer are associated with the acoustic and quasi-acoustic branches.

First, we remind that, in graphene, the vibrational density of states (VDOS in Fig. 2) has a finite constant value for energies approaching zero. This is because the ZA branch has a quadratic dispersion (not linear as usual) and the VDOS is calculated in a two dimensional Brillouin zone. On the contrary, in graphite, the VDOS goes to zero almost linearly for energies going to zero, Fig. 8. This happens in spite of the fact that in graphite, when the out-of-plane component , the ZA branch is not very different from the graphene one. In particular, on the scale of Fig. 8, the graphite ZA branch (for ) appears quadratic as the graphene one from Fig. 2. However, for graphite, the VDOS is calculated in a three dimensional BZ and the phonons have an infinitesimal weight. Also, notice in the graphite VDOS, the presence of a peak at 132 cm corresponding to the ZO frequency at . This peak and this phonon do not have a correspondence in graphene. Finally, from Fig. 8, the bilayer VDOS has a finite constant value for zero energy (as for the monolayer) and it also shows the ZO peak at 94 cm (as in graphite).

Let us consider the broadening of the LA and TA branches in Fig. 9. As already said, in graphene, these broadening do not vanish for and, for small , they are relatively constant over a wide range of values (e.g. for the LA mode we are considering the region with and in Fig. 9). This, “plateau” is due to Normal scattering towards the ZA phonons and its characteristics stem for the fact that the ZA dispersions is quadratic and that the integration (the sum in Eq. 6) is done on a two dimensional Brillouin zone. On the contrary, for three dimensional graphite, the broadening of both LA and TA modes vanishes for . It is remarkable, however, that the graphite broadening still presents a Normal scattering plateau, similar to the one of graphene, for sufficiently large . This is particularly evident for the LA mode for and in Fig. 9. The LA and TA broadening in bilayer graphene, Fig. 10, are rather more similar to the graphite one than to the graphene ones, indicating that, for a higher number of layers, the broadening should rapidly converge to the bulk graphite one. Note that, graphene bilayer presents nonvanishing broadening of the TA and LA bands at , its magnitude being about half for bilayer than for graphene.

Finally, let us consider the polarized branches. The ZA broadening in graphite is similar to the graphene one. On the other hand, the ZO broadening of graphite is much larger than the ZA one, in spite of the fact that the ZO and ZA branches are strictly related. Particularly striking is the sudden increase in the ZO broadening for certain values of q in the vicinity of (e.g for along the direction). This peak of the broadening is due to the decay of a ZO phonon, having a finite wavevector q into a ZA phonon near and a LO (or TO) phonon near . This kind of decay is possible only when the energy difference between the ZO and the ZA is equal or smaller than the energy of the LO and TO at (the LO and TO are degenerate at ), see Fig. 11. This condition is verified only for q sufficiently far from . Thus, by increasing , the sudden availability of this new decay channel produces the peak in the broadening. Finally, for the bilayer, the ZO broadening presents a structure which is not fundamentally different from the one in graphite.

### iii.3 Thermal transport

The intrinsic anharmonic thermal conductivity has been computed within the single mode time relaxation approximation using Eq. 7 . For the two dimensional materials (graphene monolayer and bilayer) we have used the convention that the volume in Eq. 7 is the surface planar unit cell multiplied by the inter-layer distance of graphite,  Å. Fig. 12 reports the thermal conductivity and its decomposition into different branch contributions. In the temperature range considered, the conductivity is almost entirely due to the acoustic and the quasi-acoustic branches. Calculations of Fig. 12 are done for  K. Below that temperature, converged results can be obtained only by using a much finer q-points grid than those presently used.

Let us, first, consider graphene. From Fig. 12, the ZA contribution increases by decreasing the temperature (actually, it diverges for ), while the LA and TA contributions are non monotonic and reach a maximum near 300 K. The difference in the two behaviors can be understood by considering that, for small , the phonons mostly occupied have small , and that, for , the anharmonic broadening (the inverse of the appearing on the r.h.s of Eq. 7) of the ZA mode goes to zero at any . This is not the case for the broadening of the TA and LA branches, Fig. 3 . Now, let us compare in Fig. 12 graphene with graphite. The ZA contribution in graphene corresponds, in graphite, to the two separate contributions ZA and ZO. These two are quite different already at room temperature. Below room temperature, the ZA increases and diverges for (as for the ZA in graphene), while the ZO does not. The TA contribution in graphene corresponds, in graphite, to the TA and TO ones. Above 200 K, the TA and TO contributions are almost indistinguishable, in spite of the fact that only the TA branch is actually acoustical. Important differences between TA and TO contributions appear only below 50 K (not shown). The same considerations hold for the LA LO couple. By comparing in Fig. 12 graphite with the bilayer, above 200 K, the overall behavior of the two systems is relatively similar, in spite of the different dimensionality. Indeed, in the same temperature range, the total conductivity of two dimensional graphene mono- and bi-layer and that of three dimensional graphite are relatively very similar, Fig. 13, with a difference between graphene and graphite of less than 10%.

Fig. 12 also reports the measured in-plane thermal conductivity. This is done only for bulk graphite, because of the abundance of experimental data. At present, experimental measurements on graphene exist only for small samples, where border effects are important, consequently the range of measured values is large and the number of samples in temperature limitedBalandin (2011). From Fig. 12, the calculated graphite conductivity (which is obtained within the SMA) grossly underestimates the measured one, by about a factor two. It is unlikely that this disagreement is due to density functional theory. Indeed, DFT reproduces accurately the measured graphite phonon dispersion, Fig. 1, suggesting that the most important quantities used in Eq. 7 are correct. On the other hand, as already explained at the end of Sec. II.1, the thermal conductivity calculated according to Eq. 7 derives from the single mode relaxation time approximation and not from an exact solution of the transport equation. Indeed, according to Refs. Lindsay et al., 2011, 2010, the SMA cannot be used to properly describe the in-plane thermal conductivity in graphitic materials and the exact solution of the Boltzmann transport equation is required. However, the results of Refs. Lindsay et al., 2011, 2010 are obtained by using a semi-empirical interatomic potential and a direct comparison with the present results is not meaningful. Further investigation on this point is required.

We now consider the thermal conductivity along the axis, perpendicular to the graphene planes. Fig. 14 shows calculations and compares them with measurements. The quasi totality of the conduction is due to the acoustic and quasi-acoustic phonons polarized along and, as expected, the conduction is much smaller than the in-plane one. The transport in-plane and the one along are quite different. The phonons relevant for the conduction have much smaller velocity and are localized around in reciprocal space. Indeed, only phonons with a nonzero dispersion along the -axis can conduct since they have a nonzero velocity and, thus, can give a contribution to the sum in Eq. 7. These phonons belong a small region surrounding the -A line. If we only integrate the contribution from a z-oriented q-space cilinder centered around , we estimate that, at 300 ,K the central 20% of the WS cell contributes about 85% of the transverse transport, but only aout 40% of the in-plane transport.

We also remark that the the conduction is extremely sensitive to the value of the lattice parameter. Indeed, a small change in results in a systematic increase or decrease of the frequencies of all the phonons relevant for the transport along , see Fig. 1. Moreover, a systematic rescale of phonon frequencies by a certain factor results in a rescale of the conduction by a factor which can be much bigger than the initial (see Eq. 7). The agreement with measurements from Fig. 14 is satisfactory; the theoretical calculations slightly overrestimates the experimental value, as we expect since we are omitting isotopic effects. We judge it compatible with the assumption that the SMA correctly describes the thermal transport along the axis.

## Iv Conclusions

We have developed and implemented a generic method for the calculation of anharmonic three-phonon scattering coefficients within density functional theory. The three phonons can have three arbitrary wavevectors . The approach works for materials with an electronic gap and also for metallic or zero gap materials. The method has been implemented in the Quantum ESPRESSO package www (); *qe and generalizes a previously existing code developed in Ref. Lazzeri and de Gironcoli, 2002. The anharmonic coefficients which can be obtained can be used in a straightforward way to compute the anharmonic broadening of a phonon with an arbitrary wavevector q and the intrinsic thermal conductivity within the single-mode relaxation time approximation (SMA). The first applications have been devoted to study graphite, graphene mono- and bi-layer.

We have reported a detailed analysis of the anharmonic phonon broadening in graphene. Interesting, the broadening of the high-energy optical branches is highly nonuniform and presents a series of sudden steps and spikes of various origin. At finite temperature, the two linearly dispersive acoustic branches TA and LA have nonzero broadening for , as noticed in Ref. Bonini et al., 2012. This anomalous behavior is due to Normal scattering towards two ZA phonons (which are quadratically dispersive), for small . The activation of the Umklapp scattering for sufficiently large is associated with a sudden increase of the broadening at nearly half the Wigner Seitz cell. We provide a set of expressions which ca be used to fit the anharmonic scattering time and broadening for the acoustic phonon branches, which are the most relevant in thermal transport.

The broadening of graphite and bi-layer is, overall, very similar to the graphene one. The most remarkable feature is the broadening of the quasi acoustical ZO branch, which is much larger and very different from the one of the strictly related ZA acoustic branch. On the other hand, the broadening of the TA and LA branches of graphite, displays a certain number of similarities with that of graphene mono and bi-layer, in spite of the different dimensionality of the systems.

Finally, we have calculated the intrinsic anharmonic thermal conductivity within the SMA. The in-plane conductivity in graphite, graphene mono- and bi-layer are very similar in spite of the differences among these systems. The calculated SMA conductivity heavily underestimates the measured one (for graphite) by almost a factor two, in the temperature range from 200 to 2000 K. This finding is compatible with the conclusions of Ref. Lindsay et al., 2011, 2010, which state that the SMA cannot be used to properly describe the in-plane thermal conductivity in graphitic materials. On the other hand, the calculated SMA conductivity for graphite along the direction perpendicular to the plane is in good agreement with measurements.

## Acknowledgments

This work was financed by ANR project accattone. Calculations were done at IDRIS (Orsay, France), Project No. 096128 and CINES Project imp6128. We thank G. Fugallo for discussions.

## Appendix A Third order calculation

This section describes the method used to calculate the third order anharmonic scattering coefficients. In order to fix the notation, Sec. A.1, and  A.2 resume DFT and linear response to DFT (alias DFPT). Third order calculations are described in Sec. A.3. Sec. A.4, and  A.5 gives the explicit expressions for certain terms. Sec. A.6 describes the implementation of the nonlinear core corrections.

### a.1 Kohn-Sham equations

Within DFT the total energy a system can be determined from the ground state electronic charge density . In turn, can be obtained by solving self-consistently the Kohn Sham (KS) equationsMartin (2004), which, in a periodic crystal are:

 [Tkin+VKS]|ψk,i⟩=ϵk,i |ψk,i⟩ (14) VKS(r)=vion(r)+δEI[n]δn(r) (15) n(r)=∑k,i~θk,i|⟨ψk,i|r⟩|2 ; ∫n(r)dr=Nel (16)

In Eq. 14, is the single-particle kinetic energy operator, is the KS potential, are the Bloch eigenstates with wavevector k, band index , and energy . , being r the position and R a lattice vector. is the external potential due to the ions, is the interaction functional (Hartree energy plus exchange-correlation contribution). The sum in Eq. 16 is done on a sufficiently fine grid of k-points. is the occupation of an electronic state: for valence band electrons and for conduction ones. Here and in the following is the integral over all the space. is the total number of valence electrons (we consider ). The total energy of the system is:

 Etot=Eion+∑k,iϵk,i~θk,i+EI[n]−∫δEI[n]δn(r)n(r)dr, (17)

where is the ionic contribution.

In this form, the KS equations are suitable for semiconductor or insulators, where the electronic gap is different from zero. if the electronic gap vanishes (metal and semi-metal case) it is customary Methfessel and Paxton (1989) to introduce a smearing function , which is characterized by a smearing width , and which becomes a step function in the limit . The KS equation are still written as in Eqs. 14-16, but now , where the Fermi energy has to be determined self consistently from

 ∑k,i~θk,i=Nel. (18)

Furthermore, in the metallic case, a proper definition of the energy requires de Gironcoli (1995) Eq 17 to include the term

 ∑k,i∫ϵF−ϵk,i−∞xδσ(x) dx, (19)

where .

### a.2 Linear response

The derivative of the electronic charge distribution with respect to the q periodic displacement, as defined in Eq. 1 (for simplicity from now on we will drop the indexes and ), can be obtained from first order perturbation theory Baroni et al. (2001). For the metallic case, linear response can be implemented following Ref. de Gironcoli, 1995; we will follow the equivalent, but slightly different approach of Ref. Lazzeri and de Gironcoli, 2002.

Let us consider a uniform grid of electronic k points and a phonon wavevector q which belongs to this grid. First, one has to solve the KS equations and obtain the ground state change density and the corresponding KS wavefunctions . Then, one has to define a “cutoff” energy which separates the electronic states which are completely empty from those which are occupied or partially occupied. In the semiconductor/insulator case, can be any energy within the gap; in the metallic case any is a reasonable choice. We define as the projector on the manifold spanned by the empty states and as the projector onto the occupied and partially occupied states.

The derivative of the charge and the derivative of the KS wavefunctions projected onto the conduction bands can be obtained by solving self-consistently the equations:

 [Tkin+VKS+αPv−ϵk,i]|ϕqki⟩=−Pc∂VKS∂uq|ψki⟩ (20)
 ∂VKS(r)∂uq=∂vion(r)∂uq+∫δ2EI[n]δn(r)δn(r′)∂n(r′)∂uqdr′ (21)
 ∂n(r)∂uq =∑k⟨r|{∑i~δkiϵqF|ψki⟩⟨ψki| +v∑i,j~θki−~θk+q,jϵki−ϵk+q,j|ψk+q,j⟩⟨ψk+q,j|Vq|ψki⟩⟨ψki| +∑i~θki[|ϕqki⟩⟨ψki|+|ψki⟩⟨ϕ−qki|]}|r⟩ (22)

is a constant chosen in such a way that the linear system of Eq. 20 is not singular. indicates that the sum is to be performed only on the partially occupied states.

The first two terms in the right hand side of Eq. 22 are different from zero only in the metallic case and are written using the notation and . for q0; it has to be determined self-consistently from

 ϵqF=∑k,i⟨ψk,i|Vq|ψk,i⟩∑k,i~δk,i. (23)

In the second line of Eq. 22 we have used a compact notation which reads: when the denominator is equal to zero one has to substitute the ratio with its limit as the denominator approaches zero. The same substitution can be used when the denominator is very small in order to gain numerical stability. Thus, when one can substitute with . We, finally, note that the present approach is different from the the one described in Ref. de Gironcoli, 1995 and that the wavefunctions presently defined are different from the of Ref. de Gironcoli, 1995.

### a.3 Third order

Let us consider three phonon displacements , , such that the sum of their wavevectors is . By solving the linear response equations one can obtain and the corresponding to the three phonons. The third derivative of the energy can then be obtained from

 ∂3Etot∂uq∂uq′∂uq′′=16[~Eq,q′,q′′+~Eq′′,q,q′+~Eq′,q′′,q+~Eq,q′′,q′+~Eq′′,q′,q+~Eq′,q,q′′] (24)
 ~Eq,q′,q′′=Zq,q′,q′′+∂3Eion∂uq∂uq′∂uq′′+∫n(r)∂3vion(r)∂uq∂uq′∂uq′′dr+3∫∂n(r)∂uq∂2vion(r)∂uq′∂uq′′dr+∭δ3EI[n]δn(r)δn(r′)δn(r′′)∂n(r)∂uq∂n(r′)∂uq′∂n(r′′)∂uq′′drdr′dr′′ (25)

For the semiconductor/insulator case, we can follow Ref. Debernardi et al., 1995 and write:

 Zq,q′,q′′=6∑k[∑i~θki⟨ϕ−qki|Vq′|ϕq′′ki⟩−v∑i,j~θki⟨ϕ−qk+q,i|ϕq′k−q′,j⟩⟨ψk−q′,j|Vq′′|ψk+q,i⟩] (26)

For the metallic case, we follow Ref. Lazzeri and de Gironcoli, 2002:

 Zq,q′,q′′=∑k{6∑i~θki⟨ϕ−qki|Vq′|ϕq′′ki⟩+6v∑i,j[~θk+q,i⟨ϕ−qk+q,i|Vq′|ψk−q′,j⟩−~θk−q′,j⟨ψk+q,i|Vq|ϕq′k−q′,j⟩]⟨ψk−q′,j|Vq′′|ψk+q,i⟩ϵk+q,i−ϵk−q′,j+2v∑i,j,l[⟨ψk,i|Vq|ψk−q,j⟩⟨ψk−q,j|Vq′|ψk+q′′,l⟩⟨ψk+q′′,l|Vq′′|ψk,i⟩+2v∑i,j,l[×~θk,i(ϵk−q,j−ϵk+q′′,l)+~θk−q,j(ϵk+q′′,l−ϵk,i)+~θk+q′′,l(ϵk,i−ϵk−q,j)(ϵk,i−ϵk−q,j)(ϵk−q,j−ϵk+q′′,l)(ϵk+q′′,l−ϵk,i)]+3ϵqF⎡⎣v∑i,j~δk,i−~δk+q′′,jϵk,i−ϵk+q′′,j⟨ψk,i|Vq′|ψk+q′′,j⟩⟨ψk+q′′,j|Vq′′|ψk,i⟩+2∑i~δk,i⟨ψk,i|Vq′|ϕq′′k,j⟩⎤⎦+3ϵqFϵq′F(∑i~δ(1)k,i⟨ψk,i|Vq′′|ψk,i⟩)−ϵqFϵq′Fϵq′′F(∑i~δ(1)k,i)}, (27)

where . Eq. 27 is written with the same compact notation of 22: when one of the denominators in Eq. 27 vanishes the corresponding term is replaced with its limit. In particular, when in the second line of Eq. 27, the argument of the sum can be written as

 −[~θk+q,i⟨ϕ−qk+q,i|ϕq′k−q′,j⟩+~δk+q,i⟨ψk+q,i|Vq|ϕq′k−q′,j⟩]⟨ψk−q′,j|Vq′′|ψk+q,i⟩.

The limits of the factor in the fourth line of Eq. 27 are

 for ϵk,i∼ϵk−q,j≠ϵk+q′′,l: ⎡⎣~θk,i−~θk+q′′,lϵk,i−ϵk+q′′,l+~δk,i⎤⎦1ϵk+q′′,l−ϵk,i for ϵk−q,j∼ϵk+q′′,l≠ϵk,i: [~θk−q,j−~θk,iϵk−q,j−ϵk,i+~δk−q,j]1ϵk,i−ϵk−q,j for ϵk+q′′,l∼ϵk,i≠ϵk−q,j: ⎡⎣~θk+q′′,l−~θk−q,jϵk+q′′,l−ϵk−q,j+~δk+q′′,l⎤⎦1ϵk−q,j−ϵk+q′′,l for ϵk,i∼ϵk−q,j∼ϵk+q′′,l: −12~δ(1)k,i. (28)

Finally, in the fifth line of Eq. 27, when one can substitute with .

Once the the derivative in Eq. 24 has been determined, one can obtain the phonon scattering coefficients by combining Eq. 4 and Eq. 5. Provided that G is vector of the reciprocal lattice, we also remark that

 ∂3Etot∂uq∂