Enhancement of thermoelectric figure-of-merit of Graphene upon BN-doping and sample length reduction

Enhancement of thermoelectric figure-of-merit of Graphene upon BN-doping and sample length reduction

Ransell D’Souza ransell.d@gmail.com; ransell.dsouza@bose.res.in    Sugata Mukherjee sugata@bose.res.in; sugatamukh@gmail.com Department of Condensed Matter Physics and Materials Science, S.N. Bose National Centre for Basic Sciences, Block JD, Sector III, Salt Lake, Kolkata 700106, India

Using first-principles density functional perturbation theory based calculations of length-dependent lattice thermal conductivity () and using our previously calculated results (Phys Rev B 95 085435 (2017)) of electrical transport, we report results of thermoelectric figure-of-merit () of monolayer and bilayer Graphene. We find nearly ten-fold increase in for the graphene sample doped with boron nitride and reduced sample length. We also compare calculated using the iterative real space method with conventional analytical Callaway-Klemens method and obtain the flexural (ZA) phonon modes to be dominant in thermal transport unlike in the latter method. Our calculations are in good agreement with available experimental data.

I Introduction

Graphene is single-atom thick, hybridized carbon atoms arranged in a two-dimensional (2D) honeycomb crystal structure having two atoms in its unit cell Geim and Novoselov (2007). Bilayer graphene consists of two monolayers, with four atoms in the unit cell arranged in an AB-type stacking known as Bernal-stacked. Monolayer graphene (MLG), has a linear band dispersion and is a semi-metal with a high electrical mobility Bolotin et al. (2008a); Zuev, Chang, and Kim (2009a); Novoselov et al. (2005); Zhang et al. (2005); Schedin et al. (2007); Tan et al. (2007); Chen et al. (2008). On the other hand bilayer graphene (BLG) can be used as a tunable band gap semiconductor Min et al. (2007). These intriguing properties along with the fact that MLG and BLG are extremely atomically stable make them ideal candidates for a good testing ground since they can be supported between two leads. Being a planar 2D structure, one of the main advantage of MLG and BLG over many other materials is that they can be readily used in circuit designs with standard lithography methods and hence can be used to fabricate transistors at true nanoscale limit. In the interest of thermoelectric devices and fundamental physics, MLG and BLG have therefore drawn a large extent of experimental, computational and theoretical attention for their transport properties.

Graphene and other related 2D nanomaterials exhibit crucial important properties which are useful for their application in renewable energy Repp et al. (2018); Pham et al. (2016). Graphene nanoribbons are found to be promising candidates for power generation in thermoelectric devices. Experimentally it has been shown that the thermal transport of graphene based devices can be tailored by defects, isotope engineering, edge roughness and by techniques to introduce nano-holes Bonaccorso et al. (2015). From a theoretical point of view, study of the lattice thermal conductivity of MLG and BLG is of great significance Lindsay et al. (2014); Saito, Mizuno, and Dresselhaus (2018).

The quality for being able to efficiently generate thermoelectric power in transport devices is related to the dimensionless quantity figure-of-merit (), expressed by . Where, is the Seebeck coefficient, is the electrical conductivity, is the temperature and is total (contribution due to electrons and the lattice) thermal conductivity. An efficient material will hence require a high power factor () and low thermal conductivity. We have recently D’Souza and Mukherjee (2017a) showed that BN-doped MLG and BLG exhibit increase in the power factor, as such doping induces a small band gap, thereby increasing . Therefore calculations on thermal conductivity of these materials are highly desirable because if the thermal conductivity decreases, doping or inducing graphene with impurities would be a useful technique to increase the power performance of graphene devices.

The relative contributions to heat conduction by the acoustic in-plane and out-of-plane phonons are still debatable. Several studies Kuang et al. (2016); Lindsay and Broido (2011); Lindsay, Broido, and Mingo (2011); Seol et al. (2010); Lindsay, Broido, and Mingo (2010) indicate that out-of-plane ZA phonon modes to be the most dominant while few others Shen et al. (2014); Alofi and Srivastava (2013); Kong et al. (2009); Aksamija and Knezevic (2011); Nika et al. (2009a); Nika, Pokatilov, and Balandin (2011); Nika and Balandin (2012); Wei et al. (2014) report the opposite. Therefore an analytical expression for the acoustic modes would play a very important role in solving these discrepancies. In this paper, we attempt to resolve such issues.

Allen has shown that the Callaway method Callaway (1959) underestimates the suppression of the normal processes and has proposed an improved method Allen (2013) which has been compared with the iterative method before for various materials Ma and Luo (2014). However, this has been done only for three dimensional materials. There have also been earlier studies on the length dependence thermal conductivity of single layer graphene using either a Monte Carlo simulation on a quadratic and linear fit to the acoustic phonon modes Mei et al. (2014) or the improved Callaway model Majee and Aksamija (2016); Guo and Wang (2017).

In all of these calculations, none of the relaxation times were calculated beyond the relaxation time approximation (RTA) using an iterative method. All length dependence calculations were done at room temperature. Moreover, those calculations do not take the symmetry of the sample into account.

Previous studies using density functional perturbation theory (DFPT)to calculate the thermal conductivity by solving the BTE exactly have been reported for graphene Lindsay et al. (2014), bilayer graphene Fugallo et al. (2014), and, B and N doped graphene Polanco and Lindsay (2018). For example, Fugallo et al.Fugallo et al. (2014) have solved the BTE exactly for MLG and BLG using phonon-phonon scattering rates derived from DFPT. Their calculations concentrate on the collective phonon excitations in comparison to the single phonon excitations. Our calculations, on the other hand, deal with solving the BTE beyond the RTA for each acoustic mode at various lengths and temperatures. The temperature dependent behavior of each acoustic mode using the iterative method has been demonstrated by Lindsay et. al Lindsay et al. (2014) only for MLG.

Previous theoretical calculations on mode dependent lattice thermal conductivity for BLG have used the Tersoff potential Lindsay, Broido, and Mingo (2011) and hence a first-principle calculation showing the length, temperature dependence of the lattice thermal conductivity is of great importance. Concentration dependent lattice thermal conductivity for B and N doped graphene calculations suggest a decrease in lattice thermal conductivity Polanco and Lindsay (2018), a feature similar to what is seen in our calculations. Reproducing these previously reported results for graphene justifies our predictions for the temperature dependent lattice thermal conductivity for BLG. Merging the lattice thermal conductivity calculations, which are in good agreement with previous theoretical and experimental data, with the calculations of the electrical transport parameters ensure that our calculated figure of merit is accurate for MLG, BLG, and BN doped MLG.

In this paper we provide an analytical solution using the Callaway-Klemens method using a quadratic fit to the out-of-plane ZA and linear fit to the in-plane LA and TA acoustic phonon dispersion, respectively, and show that the results are in perfect agreement with experiments for MLG at room temperature only. The Callaway-Klemen method overestimates at lower temperatures for MLG and overestimates at all temperatures and lengths for BLG when compared to experimental measurements. This motivated the present study on the thermal conductivity by calculating the scattering rates for each acoustic mode beyond the RTA using the first-principles iterative ShengBTE Li et al. (2014a) method for MLG and BLG at various lengths as well as temperatures and compare our results with various experimental data. Moreover, using the scattering rates from the ShengBTE method, we also calculate the length dependent for doped MLG treating BN dimers as point defects.

Using the above accurate calculations of and our earlier calculated data D’Souza and Mukherjee (2017a) on and for these materials, we calculate the thermoelectric figure-of-merit () of MLG and BLG and study the effect of sample length and BN-doping on . Purpose of this paper is to demonstrate that the iterative ShengBTE method gives accurate results of and its dependence on sample length for both MLG and BLG without any fitting parameters in quantitative agreement with experimental data. We have also found to decrease () for MLG upon BN-doping. Our calculated decrease in is consistent with such decrease observed in oxygen defects in Graphene using Raman spectroscopy Anno et al. (2017). We have presented a comparison of our results with experimental data using the local activation model Jorio et al. (2011). This decrease in together with increase in and leads to nearly ten fold increase in the thermoelectric figure-of-merit () upon BN-doping and decrease in sample length.

In the next section our calculational method is presented together with the theoretical framework used for the calculation of lattice thermal conductivity . Results on the phonon dispersion, Grüneisen parameter, and eventually are given in subsequent sections followed by a summary.

Ii Method of Calculation

ii.1 Electronic and phonon band dispersion

Geometry optimizations were carried out on a hexagonal unit cell for both monolayer (MLG) and bilayer graphene (BLG) using the first-principles Density functional theory (DFT) as implemented in the QUANTUM ESPRESSO code Giannozzi and et al. (2009). The unit cell consists of 2 and 4 carbon atoms in the -plane for MLG and BLG, respectively. Ultrasoft pseudopotential was used to describe the exchange-correlation potential kernel in the local density approximation (LDA) Rappe et al. (1990). A vacuum spacing of 20 Å was introduced in the -direction for the periodic supercell, which avoid interactions between atoms in different planes for MLG. For BLG two such layers were taken in such supercell where Van der Waals interaction Grimme (2006) between the planes was included. For the -point sampling, we have chosen a Monkhorst-Pack Monkhorst and Pack (1976) grid of and for MLG and BLG, respectively. A 160 Ry charge density energy cut-off and 40 Ry kinetic energy cut-off were used in solving the Kohn-Sham equation self consistently with an accuracy of 10 Ry.

The phonon dispersion along the high-symmetry points in the two-dimensional (2D) hexagonal Brillouin zone (BZ) (), and phonon density of states (PDOS) were calculated on the geometrically optimized structures using the density functional perturbation theory (DFPT) Baroni, Giannozzi, and Testa (1987). A Monkhorst-Pack -grid of for MLG and for BLG was used in the self-consistent calculations with a phonon frequency threshold of 10 cm.

ii.2 Theoretical methods for calculation of

ii.2.1 Analytical Callaway-Klemens method

The lattice thermal conductivity () by Callaway-Klemens Klemens (1958); Callaway (1959) method and modified by Nika et. al. Nika et al. (2009a) for an isotropic two-dimensional system is expressed as


where, is the reduced Planck constant, is the hight between two consecutive layers, number of layers, is the Boltzmann constant. The mode dependent phonon frequency and velocity at wave vector and corresponding to branch is denoted by and . is wave vector corresponding to the Debye frequency while is the wave-vector that corresponds to the sample length () dependent low cut-off frequency. We have used this formulation to calculate for 2D single- and multilayered BN and its length dependence previously D’Souza and Mukherjee (2017b).

The total phonon relaxation time comprises of the contributions from, (i) the phonon-phonon Umklapp scattering, (ii) the boundary scattering, and (iii) scattering due to point defects. The phonon-phonon Umklapp scattering rate () for a given mode is given by Klemens (1958); Callaway (1959); Nika et al. (2009a)


where, is the total mass of the atoms in the unit cell and is the mode dependent Grüneisen parameter at wave vector . The Grüneisen parameter () was calculated by us D’Souza and Mukherjee (2017a) recently by applying a 0.5% biaxial strain to MLG and BLG using .

The rough boundary scattering rate are shown to be given by Nika et al. (2009a),


where, is the width of the sample. The specularity parameter () depends on the roughness of the edges. For example an ideal smooth sample would have a specularity parameter . The scattering rate due to point defects is written asNika et al. (2009a),


where, is a dimensionless parameter to determine the strength of the point-defect scattering, given by , where is the average atomic mass, is the fractional concentration of the impurity atoms with mass . The cross-sectional area per atom of the lattice is denoted by . Each of the mentioned scattering rates can be combined to calculate the total phonon relaxation time which is given by the Matthiessen’s rule Ashcroft and Mermin (1976),


The integral in Eq.1 over the entire BZ would require setting the lower cut-off to zero which would lead to an infinite thermal conductivity. Klemens Klemens and Pedraza (1994) offered a physical reason for selecting the cut-off frequency corresponding to the ZO mode at the point for bulk graphite. This method would work with BLG but would not for MLG due to the absence of the ZO mode.

Another approach to avoid the divergence is by curbing the phonon mean free path (MFP) on the boundaries of the sheets Nika, Pokatilov, and Balandin (2011). This is achieved by making a condition that the MFP cannot exceed the physical size of the sample. This has been incorporated in both the approaches in the present paper.

By changing the integral variable from wave-vector () to phonon frequency (), the upper and lower limits of the integral in Eq. 1 would then correspond to the Debye () and cut-off ( frequency, respectively, which has been described in a recent paper D’Souza and Mukherjee (2017b).

A simple analytical expression for the mode-dependent using the Callaway-Klemens method with point defects and boundary scattering is given in the appendix.

ii.2.2 Iterative real-space method (ShengBTE)

We have also used an alternative first-principles iterative method based on DFT and DFPT utilizing third-order anharmonic inter-atomic force constants for calculating of MLG and BLG, as implemented in ShengBTE code Li et al. (2014a). In this method the lattice thermal conductivity tensor is calculated by solving the phonon Boltzmann transport equation from the converged set of phonon scattering rates using the expression,


Here, is the volume of the unit-cell, denotes the number of -points in the BZ sampling. Here, , and are the total phonon relaxation time, frequency and the phonon group velocity, respectively. denotes the Bose-Einstein distribution function and accounts for the correction to the phonon group velocity due to the deviation from the relaxation time approximation Li et al. (2014a); D’Souza and Mukherjee (2017b). It must be noted that in the iterative ab-initio method, the scattering rate processes can not be classified into Umklapp or Normal scattering processes as done in the Callaway-Klemens method (Eq. 2). The solutions to the self consistent equations are the combined three-phonon scattering rates Li et al. (2014a).

Lattice thermal conductivity calculations for MLG and BLG using ShengBTE method, having two and four atoms in its unit cell, using third nearest neighbor interactions for a supercell yields 72 and 156 displaced supercell configurations, respectively. From a set of third order derivatives of energy, calculated by implementing the plane-wave method on these displaced supercell configurations, the third-order anharmonics IFCs are constructed. The number of such configurations increase exponentially with increase in number of atoms in the unit cell. For example, five layered boron nitride, having 20 atoms in the unit cell, the same number of nearest neighbor interaction for the calculation of required 828 configurations as used earlier D’Souza and Mukherjee (2017b).

Finally, from our previous results D’Souza and Mukherjee (2017a) of the electrical conductivity () and the Seebeck coefficient (), we obtain the figure-of-merit () of undoped MLG and BLG and BN-doped MLG, using , where is the electronic thermal conductivity, found to be times smaller than .

Iii Results and discussions

iii.1 Phonon dispersion, density of states and the Grüneisen parameters

In Fig. 1 we show the calculated phonon dispersion along the high symmetry -points in the 2D irreducible hexagonal BZ for MLG and BLG D’Souza and Mukherjee (2017a) and the phonon densities of states (PDOS). Using standard group theoretical methods, it can be shown D’Souza and Mukherjee (2017b) that for MLG the six allowed phonon modes are represented by the point-group (PG) symmetry + + + , whereas for BLG the twelve allowed modes are represnted by the PG symmetry 2 + 2 + 2 + 2. Transitions corresponding to the basis are Infrared active while transitions corresponding to product of those basis () are Raman active. The momentum conservation requires that the first-order Raman scattering processes are limited to the phonons at the center of BZ (). Calculated allowed Raman and Infrared phonon frequencies for MLG and BLG and available experimental data are shown in table 1.

Figure 1: Calculated phonon dispersion and phonon density of states of MLG (above) and BLG (below) along the high symmetry points of the 2D hexagonal Brillouin zone. The magenta dashed lines are the best linear and quadratic fit to the in-plane and out-of-plane wave dependent fit to the phonon dispersion.

We have used linear fit to the in-plane acoustic phonon modes (LA, TA) and a quadratic fit to the out-of-plane acoustic mode (ZA), for obtaining an analytical form of mode-dependent for these materials as discussed in the Appendix.

(cm) Expt. (Sys.) MLG (P.G. Sym.) BLG (P.G. Sym.)
867.9 (MLG) 907 () 915 ()
1579.7(MLG) 1580.0 () 1544 ()
1579.7 (MLG) 1580.0 () 1540 ()
99 (BLG) - 108 ()
32 (BLG) - 22.16 ()
32 (BLG) - 22.16 ()
Expt (Sys.) MLG BLG
1.99 (MLG) 1.85 1.89
  • Experimental data, Ref. Malard et al. (2009); Yanagisawa et al. (2005).
    Experimental data derived from overtone Raman peaks, Ref. Lui and Heinz (2013).
    Experimental Raman data, Ref. Tan et al. (2012).
    Experimental Raman data, Ref. Mohiuddin et al. (2009).

Table 1: Phonon frequencies at the point derived from our earlier calculations in comparison with experimentally measured Raman frequencies.

As in the case of the first and second order Raman and Infrared spectroscopy Lui and Heinz (2013); Tan et al. (2012), where our calculations based on harmonic IFCs are justified experimentally, the measurements based on the pressure dependence of Raman lines Mohiuddin et al. (2009) sheds light on the anharmonic IFCs. The mode dependent Grüneisen parameters () can be obtained experimentally at the high-symmetry point using the mentioned techniqueMohiuddin et al. (2009).

Grüneisen parameters along the to K direction for MLG and BLG are plotted in Fig. 2. Alongside the first-principles Grüneisen parameter calculations, we plot the fitted constant values that we use in our study to calculate the contributions to from the in-plane acoustic phonon modes (LA, TA) and the fitted inverse square dependence on the wave-vector () for the out-of-plane acoustic mode (ZA).

The inverse wave vector squared () dependence of for the ZA mode can be easily understood from the definition of Grüneisen parameter for two-dimensional crystal D’Souza and Mukherjee (2017b), , and from the -dependence of the phonon dispersion of the ZA mode. Here, and denote the lattice constant and its change under strain, respectively. Under small positive and negative strain , since the second term in will not depend on .

Figure 2: The calculated Grüneisen () parameters of all the modes along the to K direction of the 2D Brillouin zone of the hexagonal unit cell for MLG(left) and BLG(right). The maroon dashed lines are the best constant and inverse squared wave dependent fits to the in-plane (LA,TA) and out-of-plane (ZA) parameters, respectively.

It is clear from the figure that the parameters for the in-plane modes (LA,TA) do not deviate much from their average value justifying a constant approximation made in our calculations. Similarly, the fit to the parameters using the fitted inverse square wave dependence for the out-of-plane ZA mode is in good agreement with our first-principle calculations. The parameters provide information on the degree phonon scattering and anharmonic interactions between lattice waves Mounet and Marzari (2005). From Eq. 2 it is clear that the relaxation time is highly dependent on which in turn is dependent on .

In the next section we calculate, using the phonon dispersion and the Grünesien parameters, the lattice thermal conductivity with and without point defects.

iii.2 Lattice thermal conductivity using the Callaway-Klemens method

We first obtain the analytical solution of for an ideal sheet of MLG and BLG and then calculate numerically for them with defects and specularity parameter () with values other than one.

Figure 3: Temperature dependence of using the analytical solutions of the Callaway-Klemens method for each of the acoustic modes at a constant length for (a) MLG and (b) BLG. The blue dots are experimental values of at room temperature for sample lengths 1.4m and 5m for MLG and BLG, respectively. Inset: Length dependence of at constant temperatures, =120K and =300K. The maroon dotted lines are the length dependence with point defects with parameters used to fit the experimental data Pettes et al. (2011); Xu et al. (2014).

Fig. 3 shows the acoustic mode dependent as a function of temperature for (a) MLG and (b) BLG at a constant length of 1.4 m. The insets in Fig. 3 (a) and (b) are the length dependent at two constant temperatures 120K and 300K. Length is defined as the direction along which the heat propagates. The fitting parameters, discussed in the previous section, used in this study are shown in table 2. The length dependent analytical form for for MLG at RT is in excellent agreement with a recent experiment Xu et al. (2014). However, our results of overestimate the experimental data at =120K. One explanation for this could be that in the lower temperature range (0-100 K), increases rapidly and hence a small change in the temperature in this range would result into a large change in the lattice thermal conductivity, making comparison with experimental data difficult. Another possible reason could be that we have considered a sample with an ideal sheet without any form of defects or impurity. We find that, using a large specularity parameter of and an extremely small value of =0.001 of Eq. 4, our length dependent calculations with point defects for MLG agree with experimental measurements at =120K.

The major difference between the thermal conductivity of MLG and BLG is due the out-of-plane ZA phonon mode. MLG has a total of twelve process involving the flexural phonons (ZA). Seol et. al. Seol et al. (2010) obtained a selection rule for the three-phonon scattering rates stating that only an even number of ZA phonons is attributed to each process. The four allowed processes involving flexural phonon-modes have been listed by Shen et. al. Shen et al. (2014). Therefore, the scattering rate given by Eq. 2 needs to be multiplied by three for the case of MLG. Our calculations suggest that the LA and TA modes contribute maximum to the total .

MLG 18021.5 12968.4 1.70 0.65 5.64 -7.7
BLG 18014.3 12624.9 1.75 0.72 5.89 -7.47
Table 2: Parameters used in the analytical solutions of the Callaway-Klemens method.

iii.3 Lattice thermal conductivity using the Iterative method

Figure 4: The calculated mode-dependent () and total plotted as a function of sample-length in logarithmic scale(a,b) and linear scale (c,d) at two temperatures, K and K. Solid(dashed) curves refer to calculation on MLG(BLG). The orange diamond points are the values of D’Souza and Mukherjee (2017a) at the thermodynamic limit (). (d) Zoomed box in (c) comparing our calculations with available experimental data Pettes et al. (2011); Xu et al. (2014).

Figs. 4(a) and 4(b) show, in logarithmic scale, the length dependence of the contribution from each of the acoustic modes (LA,TA,ZA) to the total lattice thermal conductivity () at two fixed temperatures, K and K, respectively, for MLG and BLG, calculated using the iterative ShengBTE method Li et al. (2014a). The iterative method clearly shows that the out-of-plane acoustic ZA mode contributes the most to the total lattice thermal conductivity. At room temperature (RT) and at the thermodynamic limit, the contribution for MLG (BLG) are 79% (70%), 19% (26%) and 2% (4%) from the ZA, TA and LA modes, respectively.

We find for BLG, a 9% drop in in comparison to that of MLG due to the ZA mode. The major difference between the phonon dispersions between MLG and BLG is the additional out-of-plane optical mode ZO. Due to this additional low-frequency mode, more phase space states are now available for phonon scattering and is one of the reason for the decrease in in BLG. As evident from Fig. 4(b), we find at small sample lengths and at lower temperatures, the mode dependent contributions to are identical for both MLG and BLG indicating that the phonon transport is ballistic and independent to the number of layers.

In Fig. 4(c), we plot the total as a function of sample length. Since most of the lattice thermal conductivity measurements were carried out at small sample lengths, in order to compare our calculations to experimental data, we show in Fig. 4(d) the zoomed data in the thin rectangular box of Fig. 4(c) where experimental measurements are available for the given sample length range.

The orange diamond points shown in Figs. 4(a), 4(b) and 4(c) are values of at the thermodynamic limit, reported previously D’Souza and Mukherjee (2017a), at the corresponding temperatures. The thermodynamic value of for MLG at 120K is higher than that its value at RT while reverse is case for BLG. This suggests that the temperature dependence of has a peak closer to 120K for MLG, whereas this peak shifts to a higher value for BLG. Lindsay el. al. Lindsay et al. (2014) have shown that the mode dependence of for MLG depends only slightly on strain. The length dependent calculations of using first principles calculation based on DFPT Lindsay et al. (2014) at RT, referring to each of the acoustic modes, for the unstrained MLG are in very good agreement with our calculations shown in Fig. 4(a). Our length dependent calculations are in excellent agreement with earlier theoretical calculations Fugallo et al. (2014) shown in Fig. 4(d).

Figure 5: The calculated mode-dependent contributions to at three different sample lengths (a,b), and total (c,d) as a function of temperature. Black circular dots are available experimental data Pettes et al. (2011); Xu et al. (2014).

Fig. 5 shows the temperature dependence of each of the acoustic modes and the total lattice thermal conductivity at three constant lengths, =1.4 m, 5 m and 9 m, calculated using the ShengBTE method Li et al. (2014a), along with the available experimental data. The ZA out-of-plane mode is shown to be the most sensitive to length as compared to the in-plane, LA and TA modes. This suggests that the ZA phonons travel ballistically in the sheets while the TA and LA modes travel diffusively. Measurements of graphene Seol et al. (2010) on a SiO substrate show a reduction in which has been explained with a scattering model where the contributions from the out-of-plane are the most dominant, in line with our calculations.

Experimental results at the thermodynamic limit () of at room temperature for graphite show a value of 2000 WmKFugallo et al. (2014). Our calculated thermodynamic limit of for BLG is 1700 WmK. This proximity of between BLG and graphite implies that the interlayer interactions are short ranged.

In Figs. 5(a) and 5(b), it can be seen that at low temperatures, the ZA mode is always larger than the in-plane acoustic modes (LA,TA). This behavior can be understood by considering the phonon density of states (PDOS) which is proportional the number of phonon-modes per frequency interval Lindsay et al. (2014). Using the definition of the 2D density of states, , one can measure the contributions from each phonon modes to the total thermal conductivity. Denoting and as the PDOS for the out-of-plane and in-plane modes, it can be easily shown that, assuming a quadratic () and linear () fit to the out-of-plane and in-plane phonon modes, respectively, . Where, , (=LA,TA), are the fitting parameters to the phonon velocity and phonon frequency shown in Table 2 and is plotted in Fig. 1. Substituting the values from Table 2, it is evident that at the long wavelength limit (), .

iii.4 Comparison between Callaway-Klemens and Iterative method

Though the calculations of the total for MLG using the Callaway-Klemens and the iterative ShengBTE methods exhibit excellent agreement, the results on mode-dependent differ sigficantly. Our analytical solutions (see Appendix) using the Callaway-Klemens method suggest that, due to the quadratic nature of phonon dispersion of the out-of-plane ZA phonon-mode and the large negative values of Grüneisen parameters, contribution from the ZA mode should contribute the least to the total lattice thermal conductivity. On the other hand the LA and TA phonon-modes exhibit smaller Grüneisen parameters and linear phonon dispersion, making the phonon group velocities almost constant and large along the boundary of the Brouillon zone.

The iterative ShengBTE method yields, due to the PDOS and symmetry of MLG, the ZA phonon-modes to contribute the most in the total as discussed in the previous subsection. Therefore, calculating the thermal conductivity of BLG using the two mentioned methods should end this disparity since the selection rules for the ZA modes is broken for BLG as compared to that in MLG. It should be noted that, in both methods, the difference between the total in MLG and BLG is due to the contributions of the ZA modes. Note, that a multiplicative factor of three in the scattering rate (Eq.2) in Callaway-Klemens method for MLG was introduced heuristically which is not an artefact of the theory, as discussed in Section III(B). Absence of the multiplicative factor would imply that the thermal conductivity of MLG would be similar to that of BLG since the magnitude of phonon group velocity and Grüneisen parameters are similar for both the materials. Calculations of by Kong et. al. Kong et al. (2009) have reported that , in line with our calculations using the analytical form without the multiplicative factor D’Souza and Mukherjee (2017b). However, lattice thermal conductivity experimental data at RT yield 0.68 Li et al. (2014b). The iterative ShengBTE method at RT yields 0.60 . As seen from Fig.3, at room temperatures calculations using the Callaway-Klemens method overestimates for BLG.

Moreover, length dependence calculations of of MLG at 120K using the iterative method are in decent agreement with experimental data at the same temperature, unlike in Callaway-Klemens method where point defects with additional fitting parameters were required to fit the calculations with the experiment (Fig. 3). Since our calculations using the iterative method are in excellent agreement with the experimental measurements at different lengths and temperatures Li et al. (2014b); Xu et al. (2014); Pettes et al. (2011) along with other available theoretical calculations based on Tersoff empirical inter-atomic potentials Lindsay, Broido, and Mingo (2010, 2011) and first-principles calculations Lindsay et al. (2014), the relaxation times calculated by iterative method are accurate as compared to those calculated by the Callaway-Klemens method. Calculations of implementing the iterative method on layered hexagonal boron-nitride have shown to be in excellent agreement using the Tersoff potentials Lindsay and Broido (2011, 2012) as well as first-principles methods D’Souza and Mukherjee (2017b). We therefore use only derived from the iterative method for the calculations of the figure of merit () in the next section.

iii.5 Figure of Merit of undoped MLG and BLG

Figure 6: The calculated Figure of merit () for undoped MLG (left) and BLG (right) at three different lengths, =1.4m, =5m and =9m together with its thermodynamic limit (). The black dashed line refers to the experimental data Anno et al. (2017). Inset: Calculated as a function of temperature for fixed chemical potential.

Fig. 6 shows the figure of merit () of undoped MLG and BLG at three different lengths together with thermodynamic limit (). Our calculated at in thermodynamic limit appears to be in good agreement with recent experimental data available for pristine graphene Anno et al. (2017), =0.55 (shown as black dashed line). The electrical Boltzmann transport equations using the RTA yields an electrical relaxation time () scaled electrical conductivity (). Berger et. al. Berger et al. (2006); Fang et al. (2015) have experimentally measured the resistivity to be cm which results in an electrical conductivity 7.1 . Adopting a Drude model, Tan et. al. Tan et al. (2007) have estimated as a function of charge density having values in the range 10fs-100ps. Using the lower bound for the relaxation time we obtain , as calculated by us recently D’Souza and Mukherjee (2017a), in the same range as seen experimentally Berger et al. (2006). Therefore, we use = 10 fs in all our calculations for the estimate of the figure of merit for MLG.

It should be noted that in general is a function of temperature and the electron momentum and therefore, depends on the direction in the Brillouin zone. Durczewski et al. Durczewski and Ausloos (2000) have devised a formalism to calculate the electron relaxation time and Zahedifar et al. Zahedifar and Kratzer (2018) have used this model to calculate the figure of merit of half-Heusler semiconductors. A realistic model of is of great importance but difficult to be taken into account. The implementation of in the Boltzmann transport equations in the BoltzTrap code Madsen and Singh (2006) are treated to be isotropic based on the formulation of Schulz et al. Schulz, Allen, and Trivedi (1992). Moreover, the behavior of the electrical resistivity, conductivity, mobility and Seebeck coefficients calculated using an isotropic D’Souza and Mukherjee (2017a) are in excellent agreement to various experimental measurements Zuev, Chang, and Kim (2009b); Bolotin et al. (2008b); Novoselov et al. (2005); Morozov et al. (2008), which indicates that an isotropic can be a good approximation for graphene and related materials.

In the inset of Fig. 6 we see that the for a fixed chemical potential, in the temperature range 50-300K, is larger for smaller sample lengths. Our calculated for both, MLG and BLG, are symmetric along the chemical potential. Due to the linear and parabolic electronic bandstructure of MLG and BLG, respectively; MLG has one peak while BLG has two in their as a function of chemical potential. Both, MLG and BLG are semi-metals and hence transport would occur only near the Fermi energy because for electrons away from the Fermi energy, there are no available states within a small energy window. At the Fermi energy, the is zero because the electronic density of states corresponding to chemical potential at the Fermi energy is zero.

iii.6 Decrement of and Enhancement of in BN-doped MLG

Defects are commonly considered to be destructive to the properties of a material used in solid states devices. Nonetheless, defects can occasionally be useful in supplying dopants to control their carrier concentration depending on the carriers either being -type or -type. Araujo, Terrones, and Dresselhaus (2012). Systems such as graphene have defects introduced in them for technological applications. Point defects arise within the planes of graphene mostly in the form of impurity atoms and lattice vacancies. Foreign impurities such as boron and nitrogen are common -type and -type dopants for graphene.

Micro-Raman spectroscopy is a method to characterize in-plane defects in graphene-like systems Genc et al. (2017); Araujo, Terrones, and Dresselhaus (2012); Jorio et al. (2011). The disorder-induced band, also known as the D-band, is a characteristic Raman feature in graphene-like systems. The D-band has no intensity in the absence of any defects and any given impurity that breaks the translation symmetry of the lattice introduces a D-band intensity in the Raman spectrum. Along with the D-band, the G-band in Raman spectrum also gives information in understanding defects in graphene-like materials predominantly when the impurity atoms dopes the material to change the bonding strength of the foreign species in the host carbon atom. Therefore, the ratio of the intensity of the D band to the G-band ()in the Raman spectrum plays an vital role in understanding the defects due to impurity scattering in graphene-like systems.

Study of the disorder due to defects in graphene caused by low energy Ar ion bombardment was done by Lucchese et al.Lucchese et al. (2010) using Raman scattering. This was carried out by varying the densities of the defects induced with different doses in the ion bombardment. The results of the experiment were modelled by inferring that a single impact of an ion on the graphene sheet would modify the sheet on two length scales Jorio et al. (2011). The model is known as the local activation model. The two length scales are referred to as and which are the radii of two circular areas measured from the impact point as shown in Fig. 7. The shorter radius, , is the structural disorder from the impact point and is know as the structurally-disordered region or the S-region. At distance for radii greater but smaller than causes a mixing of Bloch states near the K point and hence enhances the intensity in the D-band in the Raman spectrum. This region is termed as the activated or A-region beyond which the lattice structure is preserved and absent from any defect or impurity Jorio et al. (2011).

Figure 7: Unit cell containing 50 atoms with one BN-dimer (left) and two BN-dimers (right) embedded in graphene used in our calculation. The red and black arrow correspond to the radius of the activated region () and the structurally defective () region, respectively.

The local activation model for the ratio is a function of the average distance between two defects, and is expressed as Jorio et al. (2011); Araujo, Terrones, and Dresselhaus (2012); Anno et al. (2017),


where and are adjustable dimensionless parameters. For graphene-like materials the value of and are found to be 4.2 and 0.87 respectively Jorio et al. (2011); Araujo, Terrones, and Dresselhaus (2012); Anno et al. (2017) and are the values used in our calculations. The remaining parameters used in our paper are as follows: The average length between the defects is the same as the length of the unit cell, Å. The radii used for the A-region and S-region for one BN-dimer are Å and Å, respectively. Similarly, the radii used for the A-region and S-region for two BN-dimers are Å and Å  respectively (See Fig.7). The resulting ratio of one and two BN-dimers are calculated to be 0.253 and 0.451, respectively. The method to calculate the lattice thermal conductivity of BN-doped graphene will be discussed shortly and the results are plotted in the inset of Fig. 8 as a function of the ratio.

We have used our previous results D’Souza and Mukherjee (2017a) on electrical transport of BN-doped MLG obatined using first-principles DFT based electronic band structure and Boltzmann transport equations for the band electrons for obtaining and , which are then used to evaluate .

For calculating for BN-doped MLG, we have used the iterative ShengBTE method taking BN dimer as point defects in graphene sheets. Calculation of the thermal conductivity of doped MLG was performed by extracting the phonon frequency () dependent phonon relaxation time from the iterative method for MLG, adding the -dependent point defects (Eqs. 14,15) with calculated parameters, and solving Eqs. 8 and 9 with the new calculated phonon relaxation time derived from the Matthiessen’s rule (Eq. 5). The required parameters for one and two BN-doping were calculated to be , which enter in Eqs. 14,15.

Polanco et al. Polanco and Lindsay (2018) have calculated the scattering rates due to point defects by various atoms including boron and nitrogen in graphene. The point defect formula used in their paper is very similar to the Eqs. 14 and 15, with a linear fit (for the LA,TA modes) and a quadratic fit (for the ZA mode) to the phonon dispersion, as done by Lindsay et. al. Lindsay et al. (2014).

Figure 8: The calculated Figure of merit () of one and two BN dimers doped Graphene at three different lengths, =1.4m, =5m and =9m along with its thermodynamic limit. Inset: The lattice thermal conductivity plotted as a function of the ratio of intensity of the D and G band. The black circular points refer to the experimental data of Anno et al.Anno et al. (2017). The red circular points refer to the present calculations.

In Fig. 8 we show the figure of Merit for MLG doped with one and two BN dimers at three different sample lengths, =1.4m, =5m and =9m along with its thermodynamic limit. The two fold increase in , thereby increasing , for MLG upon doping D’Souza and Mukherjee (2017a) is attributed to the occurance of a small band gap. Further increase in for smaller sample-lengths is attributable to decrease in as shown in earlier sections. In the inset of Fig. 8 we plot the lattice thermal conductivity as a function of the ratio. We find that the behavior as a function of ratio is in decent agreement with experiments introducing defects in graphene using oxygen plasma treatment Anno et al. (2017).

Our calculations predict that is almost symmetric around the Fermi energy showing an increase with gate voltage for both n(p)-type doping. As one goes to higher values of energy (or gate voltage), there are additional peaks in separated by minima at around 0.7eV above and below the Fermi energy. The higher values of found at various energy range may lead to increased thermoelectric performance of doped graphene based devices.

Graphene is semi-metallic and gapless which leads to extremely small thermoelectric power factor (). However, a band gap is created at the Fermi level when graphene is doped simultaneously with boron and nitrogen Mukherjee and Kaloni (2012), which leads to enhancement in its thermoelectric power factor. Elaborate work have been carried out by various groups D’Souza and Mukherjee (2015); D’Souza, Mukherjee, and Saha-Dasgupta (2017); Bernardi, Palummo, and Grossman (2012) on band gap engineering of boron-nitride doped graphene by varying their constituent concentration. Therefore, graphene doped with two BN dimers have negligible around the Fermi level for a larger chemical potential range as compared to MLG doped with one BN dimer which in turn has negligible for a larger chemical potential range as compared to pristine MLG. This is the product of two BN dimers doped MLG having band gaps greater than one BN dimer doped MLG and that pristine MLG has a no band gap. In order to have a better understanding of all the peak seen when is plotted versus the chemical potential, we performed a model calculation described in Appendix B.

Iv Summary

In summary, using first-principles DFT based electronic and phonon band structure methods together with Boltzmann transport equations for electron and phonon, we have calculated the phonon dispersion and Grüneisen parameters for MLG and BLG and find our results to be in good agreement with experimental data (Raman spectroscopy and HREELS). Making a linear and quadratic fit to the in-plane and out-of-plane acoustic phonon dispersion along with constant in-plane and an out-of-plane inverse square wave vector dependent Grüneisen parameters, we find an analytical solution to the mode, length and temperature dependent lattice thermal conductivity for the Callaway-Klemens method. The Callaway-Klemens method suggests that the out-of-plane ZA modes contribute the least to the total lattice thermal conductivity due to the large negative Grüneisen parameters and vanishing velocities at the long wavelength limit and that the major contribution to are due to the in-plane modes, LA and TA, due to their large velocities and small Grüneisen parameters. The lattice thermal conductivity was also calculated beyond the RTA using an iterative method implemented in the ShengBTE code. The iterative method suggests that, in direct contrast to the Callaway-Klemens method, that the ZA modes contribute the most to the total while the in-plane modes contribute the least.

The Callaway-Klemens and iterative method both yield excellent agreement to total for MLG at RT and is the reason as to which mode contributes the most to . In order to solve this discrepancy, we calculate the mode, length, and temperature dependent of BLG since the selective rule is broken in the ZA modes for BLG. We find that the Callaway-Klemens method overestimates the thermal conductivity and additional point defects parameters are required to make the theory fit with the experiments. However, using the iterative method, we observe that all our calculations are in excellent agreement with many available experiments without the use of any fitting parameters. We therefore conclude that the thermal conductivity has its major contribution from the ZA mode and is also the most sensitive to the sample length. We also conclude that the mode dependent relaxation time calculated from the Callaway-Klemens method are not accurate and one must go beyond the RTA to solve the relaxation times especially for 2D materials like MLG and BLG.

Along with the electrical transport parameters like electrical conductivity, Seebeck coefficient and hence power factor calculated by us earlier, we have calculated the figure of merit of MLG and BLG. The lattice thermal conductivity used in our calculations were only taken from the iterative method. Our calculation for pristine graphene at the thermodynamic limit are in excellent agreement with available experimental data. We also find an enhancement of the figure of merit when the sample lengths are in order of  m as compared to that of the thermodynamic limit. Implementing the activation model, our estimate of the ratio for graphene doped with one and two BN-dimers treated as point defects are in excellent agreement with an experiment where defects were introduced by oxygen plasma treatment. Finally, we show that when pristine graphene is doped with one or two boron nitride dimers, the figure of merit is found to be enhanced over a wide range in chemical potential. We have therefore found a new route to enhance the figure of merit of graphene and hence improve graphene based devices over a wide range in gate voltage.

V Acknowledgments

We would like to thank J. Carrete, D.L. Nika and A. Balandin for their helpful correspondance. The calculations were performed in the High Performance Cluster platform of the S.N. Bose National Centre for Basic Sciences (SNBNCBS). RD acknowledges support from SNBNCBS through a Senior Research Fellowship.

Vi Appendix

vi.1 Simple analytic expression for

In order to have an analytical expression for , we make the approximations - (i) a linear phonon dispersion for the in-plane acoustic modes, (ii) a quadratic dispersion for the out-of-plane mode, (iii) constant Grüneisen parameters corresponding to the in-plane acoustic modes, and (iv) an inverse square wave-vector dependent Grüneisen parameters corresponding to the out-of-plane acoustic mode.

For the case of an ideal 2D material i.e., a material without any point defects having a specularity parameter of 1, with these approximations, substituting equation 2 in 1 it can be shown that the contribution to the total from LA, TA and ZA modes will all have a closed form Nika et al. (2009b); D’Souza and Mukherjee (2017a).

We now extend our calculations with point defects and specularity parameters for values of less than 1. With these approximations, Eq. 1 for the LA,TA and ZA mode can be easily shown to be,


where is given by and the separate relaxation times for the in-plane and out-of-plane modes with these approximation become,


vi.2 Model calculation of and of BN-doped Graphene

The behavior of the can be apprehended by studying the Seebeck coefficient. For semiconductors with small energy band gaps, the Seebeck coefficient, at a constant temperature, can be shown to be D’Souza and Mukherjee (2017a); Cutler and Mott (1969). Where, is the electrical conductivity. The electrical conductivity as a function of wave vector was derived from the wave vector dependent velocity (), . is the energy dispersion derived from the electronic bandstructure. The energy dependent electrical conductivity and velocity are then calculated using, and , respectively. The dummy variable corresponds to the band index. In our model calculation, runs from to , two bands below and above the Fermi energy. This sections aims to understand the behavior of and hence all of the constants in our calculations are set to 1.

Figure 9: (a) Bandstructure of one BN-doped graphene, (b) group velocity of electrons belonging to the two closest bands to the Fermi energy, shown in red and blue, (c) electrical conductivity of these electron system, (d) their Seebeck coefficient. The blue and red curves in (b,c,d) refer to the bands of the same colour as in (a). The green circles in (d) are the first-principles calculations of taking contributions from all bands D’Souza and Mukherjee (2017a).

Fig. 9 (a) shows the the bandstructure of one BN dimer doped graphene. The red curves are the two bands closest to the Fermi energy and the blue curves are the next closest. The energy dependent velocity, electrical conductivity and Seebeck coefficient are plotted in Fig. 9 (b,c,d) respectively. The colour conventions for these curves correspond to the colour of the bands in Fig.9 (a). It is evident from our calculations that the zeros in the Figure of Merit are due to the vanishing electron velocities and hence electrical conductivities. Our results using this model calculation show that the features at 1 eV, which are absent in pristine graphene, are due to the bands which are second to the closest bands to the Fermi energy. The form of the Seebeck coefficient shown in Fig. 9 (d) is in decent agreement to the form calculated using the Boltzmann equations implemented in the BOLTZTRAP code, shown in green circlesD’Souza and Mukherjee (2017a); Madsen and Singh (2006).


  • Geim and Novoselov (2007) A. K. Geim and K. S. Novoselov, Nature materials 6, 183 (2007).
  • Bolotin et al. (2008a) K. I. Bolotin, K. J. Sikes, J. Hone, H. L. Stormer,  and P. Kim, Physical Review Letters 101, 096802 (2008a).
  • Zuev, Chang, and Kim (2009a) Y. Zuev, W. Chang,  and P. Kim, Physical Review Letters 102, 096807 (2009a).
  • Novoselov et al. (2005) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos,  and A. A. Firsov, Nature 438, 197 (2005).
  • Zhang et al. (2005) Y. Zhang, Y. W. Tan, H. L. Stormer,  and P. Kim, Nature 438, 201 (2005).
  • Schedin et al. (2007) F. Schedin, A. K. Geim, S. V. Morozov, E. W. Hill, P. Blake, M. I. Katsnelson,  and K. S. Novoselov, Nature Mater. 6, 652 (2007).
  • Tan et al. (2007) Y. W. Tan, Y. Zhang, K. Bolotin, Y. Zhao, S. Adam, E. H. Hwang, S. D. Sarma, H. L. Stormer,  and P. Kim, Phys. Rev. Lett. 99, 246803 (2007).
  • Chen et al. (2008) J. H. Chen, C. Jang, S. Adam, M. S. Fuhrer, E. D. Williams,  and M. Ishigami, Nature Physics 4, 377 (2008).
  • Min et al. (2007) H. Min, B. Sahu, S. K. Banerjee,  and A. H. MacDonald, Phys. Rev. B 75, 155115 (2007).
  • Repp et al. (2018) S. Repp, E. Harputlu, S. Gurgen, M. Castellano, N. Kremer, N. Pompe, J. Wörner, A. Hoffmann, R. Thomann, F. M. Emen, S. Weber, K. Ocakoglubf,  and E. Erdem, Nanoscale 10, 1877 (2018).
  • Pham et al. (2016) C. V. Pham, S. Repp, R. Thomann, M. Krueger, S. Weber,  and E. Erdem, Nanoscale 8, 9682 (2016).
  • Bonaccorso et al. (2015) F. Bonaccorso, L. Colombo, G. Yu, M. Stoller, V. Tozzini, A. C. Ferrari, R. S. Ruoff,  and V. Pellegrini, Science 347, 1246501 (2015).
  • Lindsay et al. (2014) L. Lindsay, W. Li, J. Carrete, N. Mingo, D. A. Broido,  and T. L. Reinecke, Phys. Rev. B 89, 155426 (2014).
  • Saito, Mizuno, and Dresselhaus (2018) R. Saito, M. Mizuno,  and M. S. Dresselhaus, Phys. Rev. Appl. 9, 024017 (2018).
  • D’Souza and Mukherjee (2017a) R. D’Souza and S. Mukherjee, Phys. Rev. B 95, 085435 (2017a).
  • Kuang et al. (2016) Y. Kuang, L. Lindsay, S. Shi, X. Wang,  and B. Huang, International Journal of Heat and Mass Transfer 101, 772 (2016).
  • Lindsay and Broido (2011) L. Lindsay and D. A. Broido, Phys. Rev. B 84, 155421 (2011).
  • Lindsay, Broido, and Mingo (2011) L. Lindsay, D. A. Broido,  and N. Mingo, Phys. Rev. B 83, 235428 (2011).
  • Seol et al. (2010) J. H. Seol, I. Jo, A. L. Moore, L. Lindsay, Z. H. Aitken, M. T. Pettes, X. Li, Z. Yao, R. Huang, D. Broido, N. Mingo, R. S. Ruoff,  and L. Shi, Science 328, 213 (2010).
  • Lindsay, Broido, and Mingo (2010) L. Lindsay, D. A. Broido,  and N. Mingo, Phys. Rev. B 82, 115427 (2010).
  • Shen et al. (2014) Y. Shen, G. Xie, X. Wei, K. Zhang, M. Tang, J. Zhong, G. Zhang,  and Y. W. Zhang, Appl. Phys. Lett. 115, 063507 (2014).
  • Alofi and Srivastava (2013) A. Alofi and G. P. Srivastava, Phys. Rev. B 87, 115421 (2013).
  • Kong et al. (2009) B. D. Kong, S. Paul, M. B. Nardelli,  and K. W. Kim, Phys. Rev. B 80, 033406 (2009).
  • Aksamija and Knezevic (2011) Z. Aksamija and I. Knezevic, Applied Physics Letters 98, 141919 (2011).
  • Nika et al. (2009a) D. L. Nika, E. P. Pokatilov, A. S. Askerov,  and A. A. Balandin, Phys. Rev. B 79, 155413 (2009a).
  • Nika, Pokatilov, and Balandin (2011) D. L. Nika, E. P. Pokatilov,  and A. A. Balandin, Physica Status Solidi B 248, 2609 (2011).
  • Nika and Balandin (2012) D. L. Nika and A. A. Balandin, J. Phys.: Condens. Matter 24, 233203 (2012).
  • Wei et al. (2014) Z. Wei, J. Yang, K. Bi,  and Y. Chen, Journal of Applied Physics 116, 153503 (2014).
  • Callaway (1959) J. Callaway, Phys. Rev. 113, 1046 (1959).
  • Allen (2013) P. B. Allen, Phys. Rev. B 88, 144302 (2013).
  • Ma and Luo (2014) J. Ma and W. L. X. Luo, Physical Review B 90, 035203 (2014).
  • Mei et al. (2014) S. Mei, L. N. Maurer, Z. Aksamija,  and I. Knezevic, Journ. of Appl. Phys. 116, 164307 (2014).
  • Majee and Aksamija (2016) A. K. Majee and Z. Aksamija, Phys. Rev. B 93, 235423 (2016).
  • Guo and Wang (2017) Y. Guo and M. Wang, Phys. Rev. B 96, 134312 (2017).
  • Fugallo et al. (2014) G. Fugallo, A. Cepellotti, L. Paulatto, M. Lazzeri, N. Marzari,  and F. Mauri, Nano Lett. 14, 6109 (2014).
  • Polanco and Lindsay (2018) C. A. Polanco and L. Lindsay, Phys. Rev. B 97, 014303 (2018).
  • Li et al. (2014a) W. Li, J. Carrete, N. A. Katcho,  and N. Mingo, Comp. Phys. Commun. 185, 1747 (2014a).
  • Anno et al. (2017) Y. Anno, Y. Imakita, K. Takei, S. Akita,  and T. Arie, 2D Materials 4, 025019 (2017).
  • Jorio et al. (2011) A. Jorio, M. S. Dresselhaus, R. Saito,  and G. Dresselhaus, Raman Spectroscopy in Graphene Related Systems (Wiley, Weinheim, 2011).
  • Giannozzi and et al. (2009) P. Giannozzi and et al., J. Phys. Condens. Matter 21, 395502 (2009).
  • Rappe et al. (1990) A. M. Rappe, K. M. Rabe, E. Kaxiras,  and J. D. Joannopoulos, Phys. Rev. B 41, 1227 (1990).
  • Grimme (2006) S. Grimme, J. Comp. Chem. 27, 1787 (2006).
  • Monkhorst and Pack (1976) H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  • Baroni, Giannozzi, and Testa (1987) S. Baroni, P. Giannozzi,  and A. Testa, Physical Review Letters 58, 1861 (1987).
  • Klemens (1958) P. G. Klemens, in Solid State Physics: Advances in Research and Applications, Vol. 7, edited by F. Seitz and D. Turnbull (Academic, New York, 1958) p. 1.
  • D’Souza and Mukherjee (2017b) R. D’Souza and S. Mukherjee, Phys. Rev. B 96, 205422 (2017b).
  • Ashcroft and Mermin (1976) N. W. Ashcroft and N. D. Mermin, Solid State Physics (Holt, Reinehart and Winston, New York, 1976).
  • Klemens and Pedraza (1994) P. G. Klemens and D. F. Pedraza, Carbon 32, 735 (1994).
  • Malard et al. (2009) L. Malard, M. Pimenta, G. Dresselhaus,  and M. Dresselhaus, Physics Reports 473, 51 (2009).
  • Yanagisawa et al. (2005) H. Yanagisawa, T. Tanaka, Y. Ishida, M. Matsue, E. Rokuta, S. Otani,  and C. Oshima, Surf. Interface Anal. 37, 133 (2005).
  • Lui and Heinz (2013) C. H. Lui and T. F. Heinz, Phys. Rev. B 87, 121404 (2013).
  • Tan et al. (2012) P. H. Tan, W. P. Han, W. J. Zhao, Z. H. Wu, K. Chang, H. Wang, Y. F. Wang, N. Bonini, N. Marzari, N. Pugno, G. Savini, A. Lombardo,  and A. C. Ferrari, Nature materials 11, 294 (2012).
  • Mohiuddin et al. (2009) T. M. G. Mohiuddin, A. Lombardo, R. R. Nair, A. Bonetti, G. Savini, R. Jalil, N. Bonini, D. M. Basko, C. Galiotis, N. Marzari, K. S. Novoselov, A. K. Geim,  and A. C. Ferrari, Physical Review B 79, 205433 (2009).
  • Mounet and Marzari (2005) N. Mounet and N. Marzari, Physical Review B 71, 205214 (2005).
  • Pettes et al. (2011) M. T. Pettes, I. Jo, Z. Yao,  and L. Shi, Nano Lett. 11, 1195 (2011).
  • Xu et al. (2014) X. Xu, L. F. Pereira, Y. Wang, J. Wu, K. Zhang, X. Zhao, S. Bae, C. T. Bui, R. Xie, J. T. Thong, B. H. Hong, K. P. Loh, D. Donadio, B. Li,  and B. Özyilmaz, Nature Communications 5, 3689 (2014).
  • Li et al. (2014b) H. Li, H. Ying, X. Chen, D. L. Nika, A. I. Cocemasov, W. Cai, A. A. Balandin,  and S. Chen, Nanoscale 6, 13402 (2014b).
  • Lindsay and Broido (2012) L. Lindsay and D. A. Broido, Phys. Rev. B 85, 035436 (2012).
  • Berger et al. (2006) C. Berger, Z. Song, X. Li, X. Wu, N. Brown, C. Naud, D. Mayou, T. Li, J. Hass, A. N. Marchenkov, E. H. Conrad, P. N. First,  and W. A. de Heer, Science 312, 1191 (2006).
  • Fang et al. (2015) X.-Y. Fang, X.-X. Yu, H.-M. Zheng, H.-B. Jin, L. Wang,  and M.-S. Cao, Phys. Lett. A 379, 2245 (2015).
  • Durczewski and Ausloos (2000) K. Durczewski and M. Ausloos, Phys. Rev. B 97, 5303 (2000).
  • Zahedifar and Kratzer (2018) M. Zahedifar and P. Kratzer, Phys. Rev. B 97, 035204 (2018).
  • Madsen and Singh (2006) G. K. H. Madsen and D. J. Singh, Computer Physics Communication 175, 67 (2006).
  • Schulz, Allen, and Trivedi (1992) W. Schulz, P. Allen,  and N. Trivedi, Phys. Rev. B 45, 10886 (1992).
  • Zuev, Chang, and Kim (2009b) Y. Zuev, W. Chang,  and P. Kim, Physical Review Letters 102, 096807 (2009b).
  • Bolotin et al. (2008b) K. I. Bolotin, K. J. Sikes, J. Hone, H. L. Stormer,  and P. Kim, Physical Review Letters 101, 096802 (2008b).
  • Morozov et al. (2008) S. V. Morozov, K. S. Novoselov, M. I. Katsnelson, F. Schedin, D. C. Elias, J. A. Jaszczak,  and A. K. Geim, Phys. Rev. Lett. 100, 016602 (2008).
  • Araujo, Terrones, and Dresselhaus (2012) P. T. Araujo, M. Terrones,  and M. S. Dresselhaus, Mater. Today 15, 98 (2012).
  • Genc et al. (2017) R. Genc, M. O. Alas, E. Harputlu, S. Repp, N. Kremer, M. Castellano, S. G. Colak, K. Ocakoglu,  and E. Erdem, Scientific Reports 7, 1122 (2017).
  • Lucchese et al. (2010) M. Lucchese, F. Stavale, E. M. Ferreira, C. Vilani, M. Moutinho, R. B.Capaz, C. Achete,  and A.Jorio, Carbon 48, 1592 (2010).
  • Mukherjee and Kaloni (2012) S. Mukherjee and T. P. Kaloni, Journal of Nanoparticles Research 14, 1059 (2012).
  • D’Souza and Mukherjee (2015) R. D’Souza and S. Mukherjee, Physica E 69, 138 (2015).
  • D’Souza, Mukherjee, and Saha-Dasgupta (2017) R. D’Souza, S. Mukherjee,  and T. Saha-Dasgupta, Journal of Alloys and Compounds 708, 437 (2017).
  • Bernardi, Palummo, and Grossman (2012) M. Bernardi, M. Palummo,  and J. Grossman, Phys. Rev. Lett 108, 226805 (2012).
  • Nika et al. (2009b) D. L. Nika, S. Ghosh, E. P. Pokatilov,  and A. A. Balandin, Appl. Phys. Lett. 94, 203103 (2009b).
  • Cutler and Mott (1969) M. Cutler and N. F. Mott, Physical Review 181, 1336 (1969).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description