Enhancement of thermoelectric figureofmerit of Graphene upon BNdoping and sample length reduction
Abstract
Using firstprinciples density functional perturbation theory based calculations of lengthdependent lattice thermal conductivity () and using our previously calculated results (Phys Rev B 95 085435 (2017)) of electrical transport, we report results of thermoelectric figureofmerit () of monolayer and bilayer Graphene. We find nearly tenfold 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 CallawayKlemens 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 singleatom thick, hybridized carbon atoms arranged in a twodimensional (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 ABtype stacking known as Bernalstacked. Monolayer graphene (MLG), has a linear band dispersion and is a semimetal 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 nanoholes 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 figureofmerit (), 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 BNdoped 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 inplane and outofplane 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 outofplane 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 phononphonon 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 firstprinciple 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 CallawayKlemens method using a quadratic fit to the outofplane ZA and linear fit to the inplane 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 CallawayKlemen 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 firstprinciples 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 figureofmerit () of MLG and BLG and study the effect of sample length and BNdoping 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 BNdoping. 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 figureofmerit () upon BNdoping 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 firstprinciples 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 exchangecorrelation 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 MonkhorstPack Monkhorst and Pack (1976) grid of and for MLG and BLG, respectively. A 160 Ry charge density energy cutoff and 40 Ry kinetic energy cutoff were used in solving the KohnSham equation self consistently with an accuracy of 10 Ry.
The phonon dispersion along the highsymmetry points in the twodimensional (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 MonkhorstPack grid of for MLG and for BLG was used in the selfconsistent calculations with a phonon frequency threshold of 10 cm.
ii.2 Theoretical methods for calculation of
ii.2.1 Analytical CallawayKlemens method
The lattice thermal conductivity () by CallawayKlemens Klemens (1958); Callaway (1959) method and modified by Nika et. al. Nika et al. (2009a) for an isotropic twodimensional system is expressed as
(1)  
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 wavevector that corresponds to the sample length () dependent low cutoff 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 phononphonon Umklapp scattering, (ii) the boundary scattering, and (iii) scattering due to point defects. The phononphonon Umklapp scattering rate () for a given mode is given by Klemens (1958); Callaway (1959); Nika et al. (2009a)
(2) 
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),
(3) 
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),
(4) 
where, is a dimensionless parameter to determine the strength of the pointdefect scattering, given by , where is the average atomic mass, is the fractional concentration of the impurity atoms with mass . The crosssectional 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),
(5) 
The integral in Eq.1 over the entire BZ would require setting the lower cutoff to zero which would lead to an infinite thermal conductivity. Klemens Klemens and Pedraza (1994) offered a physical reason for selecting the cutoff 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 wavevector () to phonon frequency (), the upper and lower limits of the integral in Eq. 1 would then correspond to the Debye () and cutoff ( frequency, respectively, which has been described in a recent paper D’Souza and Mukherjee (2017b).
A simple analytical expression for the modedependent using the CallawayKlemens method with point defects and boundary scattering is given in the appendix.
ii.2.2 Iterative realspace method (ShengBTE)
We have also used an alternative firstprinciples iterative method based on DFT and DFPT utilizing thirdorder anharmonic interatomic 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,
(6) 
Here, is the volume of the unitcell, 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 BoseEinstein 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 abinitio method, the scattering rate processes can not be classified into Umklapp or Normal scattering processes as done in the CallawayKlemens method (Eq. 2). The solutions to the self consistent equations are the combined threephonon 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 planewave method on these displaced supercell configurations, the thirdorder 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 figureofmerit () of undoped MLG and BLG and BNdoped 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 pointgroup (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 firstorder 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.
We have used linear fit to the inplane acoustic phonon modes (LA, TA) and a quadratic fit to the outofplane acoustic mode (ZA), for obtaining an analytical form of modedependent 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 
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 highsymmetry 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 firstprinciples Grüneisen parameter calculations, we plot the fitted constant values that we use in our study to calculate the contributions to from the inplane acoustic phonon modes (LA, TA) and the fitted inverse square dependence on the wavevector () for the outofplane 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 twodimensional 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 .
It is clear from the figure that the parameters for the inplane 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 outofplane ZA mode is in good agreement with our firstprinciple 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 CallawayKlemens 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.
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 (0100 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 outofplane 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 threephonon scattering rates stating that only an even number of ZA phonons is attributed to each process. The four allowed processes involving flexural phononmodes 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 .
System 






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 
iii.3 Lattice thermal conductivity using the Iterative method
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 outofplane 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 outofplane optical mode ZO. Due to this additional lowfrequency 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).
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 outofplane mode is shown to be the most sensitive to length as compared to the inplane, 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 outofplane 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 WmK Fugallo 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 inplane acoustic modes (LA,TA). This behavior can be understood by considering the phonon density of states (PDOS) which is proportional the number of phononmodes 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 outofplane and inplane modes, it can be easily shown that, assuming a quadratic () and linear () fit to the outofplane and inplane 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 CallawayKlemens and Iterative method
Though the calculations of the total for MLG using the CallawayKlemens and the iterative ShengBTE methods exhibit excellent agreement, the results on modedependent differ sigficantly. Our analytical solutions (see Appendix) using the CallawayKlemens method suggest that, due to the quadratic nature of phonon dispersion of the outofplane ZA phononmode 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 phononmodes 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 phononmodes 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 CallawayKlemens 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 CallawayKlemens 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 CallawayKlemens 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 interatomic potentials Lindsay, Broido, and Mingo (2010, 2011) and firstprinciples calculations Lindsay et al. (2014), the relaxation times calculated by iterative method are accurate as compared to those calculated by the CallawayKlemens method. Calculations of implementing the iterative method on layered hexagonal boronnitride have shown to be in excellent agreement using the Tersoff potentials Lindsay and Broido (2011, 2012) as well as firstprinciples 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
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 10fs100ps. 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 halfHeusler 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 50300K, 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 semimetals 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 BNdoped 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.
MicroRaman spectroscopy is a method to characterize inplane defects in graphenelike systems Genc et al. (2017); Araujo, Terrones, and Dresselhaus (2012); Jorio et al. (2011). The disorderinduced band, also known as the Dband, is a characteristic Raman feature in graphenelike systems. The Dband has no intensity in the absence of any defects and any given impurity that breaks the translation symmetry of the lattice introduces a Dband intensity in the Raman spectrum. Along with the Dband, the Gband in Raman spectrum also gives information in understanding defects in graphenelike 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 Gband ()in the Raman spectrum plays an vital role in understanding the defects due to impurity scattering in graphenelike 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 structurallydisordered region or the Sregion. 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 Dband in the Raman spectrum. This region is termed as the activated or Aregion beyond which the lattice structure is preserved and absent from any defect or impurity Jorio et al. (2011).
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),
(7)  
where and are adjustable dimensionless parameters. For graphenelike 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 Aregion and Sregion for one BNdimer are Å and Å, respectively. Similarly, the radii used for the Aregion and Sregion for two BNdimers are Å and Å respectively (See Fig.7). The resulting ratio of one and two BNdimers are calculated to be 0.253 and 0.451, respectively. The method to calculate the lattice thermal conductivity of BNdoped 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 BNdoped MLG obatined using firstprinciples 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 BNdoped 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 BNdoping 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).
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 samplelengths 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 semimetallic 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 SahaDasgupta (2017); Bernardi, Palummo, and Grossman (2012) on band gap engineering of boronnitride 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 firstprinciples 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 inplane and outofplane acoustic phonon dispersion along with constant inplane and an outofplane 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 CallawayKlemens method. The CallawayKlemens method suggests that the outofplane 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 inplane 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 CallawayKlemens method, that the ZA modes contribute the most to the total while the inplane modes contribute the least.
The CallawayKlemens 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 CallawayKlemens 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 CallawayKlemens 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 BNdimers 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 inplane acoustic modes, (ii) a quadratic dispersion for the outofplane mode, (iii) constant Grüneisen parameters corresponding to the inplane acoustic modes, and (iv) an inverse square wavevector dependent Grüneisen parameters corresponding to the outofplane 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,
(8)  
(9) 
where is given by and the separate relaxation times for the inplane and outofplane modes with these approximation become,
(10)  
(11)  
(12)  
(13)  
(14)  
(15) 
vi.2 Model calculation of and of BNdoped 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.
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).
References
 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 SahaDasgupta (2017) R. D’Souza, S. Mukherjee, and T. SahaDasgupta, 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).