First-principles study of the electrical and lattice thermal transport in monolayer and bilayer graphene

TL;DR

This study uses first-principles calculations and Boltzmann transport theory to analyze electrical and thermal transport in monolayer and bilayer graphene, revealing effects of doping, electric fields, and impurities on their properties.

Contribution

It provides new insights into how doping, electric fields, and impurities influence electrical and thermal transport in graphene from first-principles calculations.

Findings

Bloch-Grüneisen behavior observed in monolayer graphene resistivity.

Boron nitride substitution doubles the Seebeck coefficient.

Calculated lattice thermal conductivities agree with experimental data.

Abstract

We report the transport properties of monolayer and bilayer graphene from first principles calculations and Boltzmann transport theory (BTE). Our resistivity studies on monolayer graphene show Bloch-Gr${\rm \ddot{u}}$neisen behavior in a certain range of chemical potentials. By substituting boron nitride in place of a carbon dimer of graphene, we predict a twofold increase in the Seebeck coefficient. A similar increase in the Seebeck coefficient for bilayer graphene under the influence of a small electric field $\sim 0.3$ eV has been observed in our calculations. Graphene with impurities shows a systematic decrease of electrical conductivity and mobility. We have also calculated the lattice thermal conductivities of monolayer graphene and bilayer graphene using phonon BTE which show excellent agreement with experimental data available in the temperature range 300-700\ K.

Tables

Table 1. Table 1: Calculated phonon frequency ω 𝜔 \omega (in cm -1 ) and the Grüneisen parameter γ 𝛾 \gamma of MLG and BLG at the high-symmetric q 𝑞 q points in the hexagonal BZ for different vibrational modes.

System	$q$		ZA	TA	LA	ZO	TO	LO	ZO’
MLG		$\omega$ $\gamma$	0 -100	0 0.779	0 1.848	907 -0.086	1580 1.850	1580 1.605	- -
	$\Gamma$
BLG		$\omega$ $\gamma$	0 -50	0 0.936	0 1.640	915 0.110	1540 1.892	1544 1.878	108 -0.498
MLG		$\omega$ $\gamma$	545 -1.245	1000 0.510	1230 1.713	545 -1.245	1370 2.584	1230 1.713	- -
	K
BLG		$\omega$ $\gamma$	554 -1.06	1046 0.581	1205 1.753	557 -1.025	1287 2.779	1206 1.749	554 -1.049
MLG		$\omega$ $\gamma$	578 -1.057	631 0.139	1350 1.953	645 -1.004	1430 2.281	1372 1.184	- -
	M
BLG		$\omega$ $\gamma$	478 -0.982	672 0.267	1327 2.040	660 -0.762	1348 2.432	1363 1.195	485 -0.956

Equations

$$v_{\alpha }(i,\text{k})=\frac{1}{\hslash }\frac{\partial {\epsilon }_{i,k}}{\partial k_{\alpha }}$$

$$\frac{{\sigma }_{\alpha \beta }(T,\mu )}{\tau }$$

$$S_{\alpha \beta }(T,\mu )=\frac{1}{eT}\frac{\int v_{\alpha }(i,\mathbf{k})v_{\beta }(i,\mathbf{k})(ϵ-\mu )[\frac{-\partial f_{\mu }(T,ϵ)}{\partial ϵ}]dϵ}{\int v_{\alpha }(i,\mathbf{k})v_{\beta }(i,\mathbf{k})[\frac{-\partial f_{\mu }(T,ϵ)}{\partial ϵ}]dϵ}$$

$$\displaystyle S=-\frac{\pi^{2}k_{B}^{2}T}{3e}{d\over dE}\,[ln\,\sigma(E)]{\Bigr{|}}_{E=E_{F}}$$

$${\kappa }_L^{\alpha \beta }=\frac{1}{k_B{T}^2\Omega N}\sum _{\lambda }{f}_0(f_0+1)(\hslash {\omega }_{\lambda }{)}^2{v}_{\lambda }^{\alpha }{\tau }_{\lambda }^0(v_{\lambda }^{\beta }+{\Delta }_{\lambda }^{\beta })$$

$$n=\frac{k_F^2}{\pi }=\frac{{ϵ}_F^2}{\pi (\hslash v_F{)}^2}$$

$$D({ϵ}_F)=\frac{dn}{dϵ}=\frac{2\sqrt{n}}{\hslash v_F\sqrt{\pi }}$$

$$\displaystyle\rho(T)=\rho(0)+A\Big{(}\frac{T}{\Theta_{BG}}\Big{)}^{m}\int^{\frac{\Theta_{BG}}{T}}_{0}\frac{x^{m}}{(e^{x}-1)(1-e^{-x})}dx$$

$$H=H_0+H_{imp}$$

$$H_0=\sum _iϵ|i⟩⟨i|-\sum _{i,j}{t}_{ij}|i⟩⟨j|$$

$$H_{imp}=\Delta \sum _i|i⟩⟨i|$$

$${\mu }_{FE}=\frac{1}{e}\frac{d\sigma }{dn}$$

$$\displaystyle\begin{split}\gamma_{s}(q)&=\frac{-a_{0}}{2\,\omega_{s}(q)}\frac{\delta\omega_{s}(q)}{\delta a}\\ &\approx\frac{-a_{0}}{2\,\omega_{s}(q)}\Big{[}\frac{w_{+}-w_{-}}{da}\Big{]}\end{split}$$

Full text

I Introduction

Transport properties in monolayer graphene (MLG) and bilayer graphene (BLG) have drawn a great deal of theoretical and experimental attention in recent years due to interest in both fundamental physics and potential electronic applications. An interesting transport property is the thermoelectric effect (TE) in which a temperature difference $\Delta T$ across a conductor creates an electric potential $\Delta V=S\Delta T$ where $S$ is the thermoelectric power, also called Seebeck coefficient, and characterizes the thermoelectric sensitivity of the conductor. Thus a material with a high value of $S$ is a necessary condition for it to be used in thermoelectric applications. The performance of the thermoelectric material in such TE devices is determined by a dimensionless parameter $ZT=\frac{S^{2}\sigma}{\kappa}T$ , where $\sigma$ is the electrical conductivity, $\kappa=\kappa_{e}+\kappa_{p}$ is the thermal conductivity including both electron and phonon contributions, and $T$ is the temperature. Hence such materials should have a high power factor $S^{2}\sigma$ and a low $\kappa$ Rowe (2006); Tritt (2001); Mahan (1998).

Extensive experiments have been carried out by Kim $et\ al$ Zuev et al. (2009); Bolotin et al. (2008) on transport properties in graphene. For ultraclean suspended graphene, Bolotin $et\ al.$Bolotin et al. (2008) have shown that resistivity is strongly temperature dependent. At low temperatures ($T\sim 5K$) the electron transport is near-ballistic having a mobility of $\sim$ 170000 ${\rm cm^{2}}/{\rm Vs}$. For large carrier densities, resistivity increases with $T$ and is linear above 50K. This suggests that the carrier scattering is from acoustic phonons. They also reported the highest measured value of mobility ($\sim$ 120000 ${\rm cm^{2}}/{\rm Vs}$) at $T$=240K. Using a micro-fabricated heating system, electrical conductivity and Seebeck coefficients were measured for the first time by Zuev $et\ al$ Zuev et al. (2009). The gate-dependent and hence carrier-density-dependent conductance and Seebeck coefficient showed good agreement with the semiclassical Mott relation. It was shown that the sign of Seebeck coefficients changes as the majority carrier density shifts from the electron to that of the hole. The gate-dependent Seebeck coefficient and Nernst signal on graphene in a magnetic field was studied by Ong $et\ al.$ Checkelsky and Ong (2009). Quantities closely related to $\sigma$ are mobility ($\mu_{FE}$) and resistivity ($\rho$). Since a field effect transistor can measure the mobility, we use $\mu_{FE}$ to define mobility so as to not confuse it with chemical potential ($\mu$). Values of $\mu_{FE}$ as high as 20000 ${\rm cm}^{2}/{\rm V\ s}$, have been reported for MLG at low temperatures $T$ Novoselov et al. (2005); Zhang et al. (2005); Schedin et al. (2007); Tan et al. (2007); Chen et al. (2008). It is well known that the resistivity $\rho$ strongly depends on the electron-phonon (e-ph) scattering in the material and varies linearly with temperature, $\rho\sim T$, at temperatures comparable to or higher than the Debye temperature $\Theta_{D}$ , reflecting a classical behavior of the phonons. The bosonic nature of the phonons becomes evident below $\Theta_{D}$ resulting in a rapid decrease of resistivity which for a normal conductor $\rho\sim T^{5}$, known as the Bloch-Gr${\rm\ddot{u}}$neisen (BG) regime. This defines a new characteristic temperature scale for the low-density e-ph scattering, the BG temperature $\Theta_{BG}=\frac{2v_{ph}k_{F}}{k_{B}}<\Theta_{D}$, where $v_{ph}$ and $k_{F}$ are the phonon group velocity and Fermi momentum, respectively. It was theoretically predicted Hwang and Sarma (2008) that in a two-dimensional systems such as graphene, there is a smooth transition in the resistivity from a linear dependence $\rho\sim T$ to $\rho\sim T^{4}$. This two-dimensional (2D) BG was observed experimentally by Efetov and Kim Efetov and Kim (2010) in graphene.

Temperature dependence of electron transport in graphene and its bilayer have also been studied by Morozov $et\ al.$ Morozov et al. (2008). They have shown that if extrinsic disorders are eliminated, mobilities higher than $2\times 10^{5}{\rm cm^{2}/Vs}$ are attainable. Ponomarenko $et\ al.$ Ponomarenko et al. (2009) studied the scattering mechanism on graphene on various substrates and found no significant changes in carrier mobility on different substrates questioning the largely believed assumption that the dominant source of scattering in graphene is the charge impurities in the substrate.

The experimental works carried out by various groups such as Kim $et\ al$ Zuev et al. (2009), Geim $et\ al$Ponomarenko et al. (2009), Ong $et\ al$ Checkelsky and Ong (2009), Nam $et\ al$ Nam et al. (2010) and Wang $et\ al$ Wang et al. (2011) on gate-dependent electron transport properties of MLG and BLG have motivated us to carry out computational studies based on first-principles calculations of MLG-doped hexagonal boron nitride (h-BN) and pure BLG under the influence of an electric field, since there is a direct correspondence between gate voltage and chemical potential. It is observed that $S$ is strongly dependent on the amount of doping in the case of MLG. A marked increase of $S$ is observed when an electric field is applied perpendicularly to the plane of the bilayer graphene sheets.

A further motivation to study the effects of impurity scattering on graphene sheets was provided by a recent paper by Ghahari $et\ al$ Ghahari et al. (2016), which reported an enhanced Seebeck coefficient in graphene due to scattering. We have used a tight-binding model (TBM) for graphene with impurities and solved the Boltzmann transport equations for electrons to study the effect of impurities on electrical conductance. Our study used a large $k$ mesh to capture correctly the enhancement of the Seebeck coefficient and a constant decrease in electrical conductivity due to doping and decrease in mobility by an order of magnitude as observed experimentally Ghahari et al. (2016); Chen et al. (2008); Tan et al. (2007).

Finally, we have reexamined the lattice thermal conductivity in MLG and BLG using the lattice Boltzmann transport method (Li et al., 2014a) employing the phonon bandstructure calculated from first-principles DFT Giannozzi and et al. (2009). Our calculated phonon dispersion agrees with previous results Yan et al. (2008); Nika et al. (2009); Kong et al. (2009); Yao et al. (2000); Maultzsch et al. (2004), whereas the lattice thermal conductivity shows quantitative agreement with recent experimental data Li et al. (2014b).

The purpose of this paper is to illustrate that all essential thermoelectric parameters of MLG and BLG, both for electrical and lattice thermal transport, can be accurately calculated using first-principles electronic and phonon band structure together with the Boltzmann transport theory. Our results are in excellent agreement with experimental data.

II Method of Calculation

II.1 Density functional theory based electronic structure calculation

All electronic structure calculations were carried out using the density functional theory (DFT) based plane-wave method as implemented in the Quantum Espresso code Giannozzi and et al. (2009). A hexagonal unit cell has been used in all these calculations. For the exchange-correlation potential, we use the generalized gradient approximation (GGA) Perdew et al. (1996) and the ultrasoft pseudopotential Vanderbilt (1990) to describe the core electrons. Monkhorst–Pack $k$-point grids Monkhorst and Pack (1976) of $120\times 120\times 1$ and $70\times 70\times 4$ were implemented for monolayer and bilayer graphene in the self-consistent calculations with a 40Ry kinetic energy cutoff and a 160Ry charge density energy cutoff, respectively. The Kohn-Sham equations are solved self-consistently to achieve an accuracy of $10^{-9}$ Ry in the total energy. Plane-wave methods incorporate periodicity and hence in order to avoid interaction between sheets we have a vacuum spacing of 22Å along the $z$ direction. The van der Waals interaction was included for multilayered systems.

II.2 Boltzmann transport theory for band electrons

The semiclassical Boltzmann transport theory (BTE) applied to band electrons, as implemented in the Boltztrap code Madsen and Singh (2006), was used to calculate the transport properties. The BTE allows us to calculate the thermoelectric parameters along two orthogonal principal axes in the plane of the two-dimensional graphene layers. Thus the calculated thermoelectric parameters are taken as average over those along the principal directions. The energies and their corresponding k-points were taken from the electronic structure calculations to deduce various transport parameters. For example, the group velocity of the electrons $v_{\alpha}$ is given by,

[TABLE]

In the relaxation-time approximation electrical conductivity is then expressed as,

[TABLE]

Detailed studies Schulz et al. (1992); Allen et al. (1993) have shown that to a good approximation the relaxation time $\tau$ is isotropic.

The Seebeck coefficient is then calculated using,

[TABLE]

where, $f_{\mu}=\frac{1}{1+e^{\beta(\epsilon-\mu)}}$ is the Fermi-Dirac distribution, $\mu$ is the chemical potential and $k_{B}$ is the Boltzmann constant. MLG and BLG are known to be semi-metallic where the transport occurs only near the Fermi level. We can therefore use the Sommerfeld expansion of Eqs. 2 and 3 to obtain,

[TABLE]

Eq. 4 is known as the Mott formula Cutler and Mott (1969).

II.3 Transport properties from tight-binding method

In order to study the effect of dilute impurities in graphene on its electrical conductivity, we have used a tight-binding (TB) model of graphene using the method described in Saito et al. (1998) for MLG and McCann and Koshino (2013) for BLG. Thermoelectrical parameters are then calculated from the Boltzmann transport equations described above. The TB model allows us to calculate using a very large unit cell containing hundred atoms or more with $k$ pont grids of $120\times 120\times 1$ for MLG and $70\times 70\times 4$ for BLG, not accessible by the DFT-based plane wave methods.

II.4 Lattice thermal conductivity from BTE

We have used the phonon Boltzmann transport method to solve the Boltzmann transport equations for phonons starting from a set of interatomic force constants (IFCs) obtained from the phonon dispersion obtained using ab initio calculations Giannozzi and et al. (2009), as implemented in the ShengBTE code Li et al. (2014a). Lattice thermal conductivity calculations require the harmonic second-order interatomic force constants (IFCs) as well as the anharmonic third-order IFCs. The harmonic IFCs were calculated using density functional perturbation theory. To attain the anharmonic IFCs, we use a real-space supercell approach. The phonon Boltzmann transport calculates the converged set of phonon scattering rates and uses them to obtain the lattice thermal conductivity $\kappa_{L}$ using the expression,

[TABLE]

where, $k_{B}$ is the Boltzmann constant, $\Omega$ is the volume of the unit cell, $N$ is total number of $q$ points in the Brillouin zone sampling, $\tau_{\lambda}^{0}$ is the relaxation time of mode $\lambda$ which is obtained from perturbation theory and used within the relaxation time approximation (RTA) and $f_{0}$ is the Bose-Einstein statistic. $\omega_{\lambda}$ and $v_{\lambda}$ are the angular frequency and group velocity of the phonon mode $\lambda$, where $\lambda$ comprises of a phonon branch index and wave vector. $\Delta_{\lambda}$ is the measure of how much the population of a specific phonon mode (and hence associated heat current) deviates from the RTA prediction. Detailed information on the work-flow can be found here Li et al. (2014a).

For both cases, MLG and BLG, a $4\times 4\times 1$ supercell along with a $5\times 5\times 2$ Monkhorst–Pack $k$ mesh was chosen to evaluate the forces. All calculations involved in the anharmonic third-order IFCs accounted for interactions up to the third-nearest neighbors. The phonon calculations were performed using an ultrasoft pseudopotential Rappe et al. (1990), which reproduces correctly the phonon dispersion of both MLG and BLG. In the present study the van der Waals interaction was included for BLG as prescribed by Grimme Grimme et al. (2010).

III Results and Discussion

III.1 Electrical conductivity and Seebeck coefficient of MLG

The calculated lattice constants for both MLG and BLG were 2.47 Å using GGA and 2.46 Å using LDA, respectively. The interlayer separation for BLG was found to be 3.32 Å.

As the precise numerical value of the electron relaxation time $\tau$ for graphene is not known and can only be estimated from experiments Tan et al. (2007), we show in Figs. 1(a,b) the scaled electrical conductivity $\sigma$/$\tau$, calculated from the Bolzmann transport theory (Eq. 2), as a function of energy and carrier concentration in the temperature range 40K$-300$K, respectively. It is observed that close to the Fermi energy $\epsilon_{F}$, $\sigma$ increases with increasing temperature but is almost independent of $T$ away from Fermi energy. This behavior has been reported by Morozov $et\ al$. Morozov et al. (2008). In Fig. 1(b) we observe $\sigma$ of MLG to be proportional to $\sqrt{n}$ where $n$ is the charge carrier density.

This $\sqrt{n}$ behavior can be explained using a tight-binding model with an energy expansion around the $k$ point to get $\epsilon_{F}=\hbar v_{F}k_{F}$, where $\hbar$ is the reduced Planck constant, $v_{F}$ is the Fermi velocity and $k_{F}$ is the Fermi wavenumber. For a 2-dimensional sample,

[TABLE]

Therefore the MLG’s density of states $D(\epsilon_{F})$ is given by,

[TABLE]

Substituting this in the Einstein relation $\sigma\propto D(\epsilon_{F})$, we obtain $\sigma\propto\sqrt{n}$. This behavior has been experimentally confirmed for pristine graphene Bolotin et al. (2008).

The Seebeck coefficient was calculated using the Mott relation Eq. 4 numerically and is plotted in Fig. 1(c) as a function of energy. Our calculated form of $S$ is in very good agreement with the experimental results reported by Zuev $et\ al$ Zuev et al. (2009). Using a back-gated field effect transistor, the Fermi energy can be tuned by adjusting the gate voltage. Zuev $et\ al$ Zuev et al. (2009) have used this method for a mesoscopic graphene sample using the formula $S=-\frac{\pi^{2}k_{B}^{2}T}{3e}\frac{1}{\sigma}\frac{d\sigma}{dV_{g}}\frac{dV_{g}}{dE}|_{E=E_{f}}$. It can be easily seen that this is identical to Eq. 4. Therefore there is a direct correspondence between the gate voltage and the chemical potential. In Fig. 1(d) we plot $S$ against $T$ for two different values of chemical potential. The linear dependence of $S$ on $T$, which has also been reported by experimental measurements by Zuev $et\ al.$ Zuev et al. (2009), suggests that the mechanism of thermoelectric transport is diffusive Scarola and Mahan (2002).

III.2 Resistivity and fit to the Bloch-Grüneisen form

In Fig. 2 we plot $\rho$ as a function of $T$ calculated from the Boltzmann transport equations for electron transport. The electron relaxation time $\tau$ of graphene depends on the degree of doping and carrier concentration as indicated by experimental measurements Tan et al. (2007) and can vary in the range 10  fs to 1  ps. For simplicity we have assumed \tau\sim 1\times 10^{-14}\s; however, our temperature dependent behaviour of $\rho$ should not depend on the choice of $\tau$. The behaviour of $\rho$ at certain values of $\mu$ having the Bloch-Gr${\rm\ddot{u}}$neisen form has been observed experimentally Efetov and Kim (2010); Bolotin et al. (2008). In order to understand the behaviour, we fit our calculated values of $\rho$ to the Bloch-Gr${\rm\ddot{u}}$neisen formula,

[TABLE]

with $m$ and $T$ as fitting parameters. Our best fit resulted in $m=4$ and $\Theta_{BG}$ as shown in Fig. 2. This suggests that in the low-$T$ regime, resistivity scales as $T^{4}$ and smoothly scales to a linear $T$ behavior at higher $T$ regimes. The $T^{4}$ behavior of resistivity reflects the two-dimensional nature of electrons in graphene. At high temperatures the quantization of lattice waves is irrelevant; therefore the scattering is proportional to the square of the amplitude of the fluctuations about their equilibrium position that is proportional to $\sqrt{T}$, which leads to the linear behavior of resistivity at higher temperatures. The reduction from the $T^{5}$ behavior of the resistivity for a typical bulk metal to the $T^{4}$ one, as seen in graphene, is due to reduced electron Fermi surface of graphene, smaller than the size of its phonon Brillouin zone. Therefore, only a small fraction of acoustic phonons will scatter off the electrons. An excellent illustration and explanation of the high, low and Bloch-Gr${\rm\ddot{u}}$neisen temperature regimes, by taking into account different sizes of the Fermi surface of graphene and a conventional metal at different temperatures, was given by Fuhrer Fuhrer (2010).

As mentioned previously, $\Theta_{BG}$ is defined as $\Theta_{BG}=\frac{2v_{ph}k_{F}}{k_{B}}$, where $v_{ph}$ is the phonon velocity. However when expressed as a function of carrier concentration, it can be easily shown Sarma et al. (2011) to have the form $\Theta_{BG}=A_{0}\sqrt{n}\,$K, with carrier concentration measured in $10^{12}$cm*-2*. In our calculation, when fitted to the above form, we obtain $A_{0}$=45.5, which is close to the earlier estimated value of 54 Sarma et al. (2011).

In the inset of Fig. 2 we plot $\Delta\rho$ ($=\rho(T)-\rho(0)$) against temperature. This parameter will now show us the increase in resistivity with respect to temperature and when plotted in logarithmic scale, we can see that the increase in resistivity in the lower regime has a $\sim T^{4}$ feature and $\sim T$ feature at higher temperatures. The red and green lines in the inset are $\propto T^{4}$ and $\propto T$ equations where the constant of proportionality was found by a best fit method. $\Theta_{BG}$ is in the same order of magnitude range as reported by Efetov and Kim (2010). We must note that this behavior of $\rho$ is found only for a certain range of chemical potential. If we compare our Seebeck coefficient result to that of the experimental results by Zuev $et\ al$ Zuev et al. (2009) we find that as we vary our chemical potential from 1 eV to -1 eV, the same behaviour is found when they vary their gate voltage from 40V to -40V. This implies that a small change in the chemical potential is equivalent to a large change in gate voltage. This is the reason why we choose only a small range of chemical potentials to demonstrate the Bloch-Gr${\rm\ddot{u}}$neisen behavior as observed experimentally. If one looks at the study by Kim $et\ al$ Bolotin et al. (2008), the linear part of resistivity (higher $T$ regime) increases quickly with an increase of gate voltage of only 1eV. As observed experimentally Efetov and Kim (2010), we found that the slope of the $\rho$ vs $T$ curves in the linear-$T$ regime increases with $\mu$. The present study of the temperature dependent resistivity, which has been theoretically predicted Hwang and Sarma (2008) and experimental observed Efetov and Kim (2010), is here done using the Boltzmann transport theory applied to band electrons.

III.3 Enchancement of Seebeck coefficient upon doping and in presence of an electric field

As mentioned previously, the performance of a thermoelectric material is measured by a dimensionless parameter, the figure of merit($ZT=\frac{S^{2}\sigma}{\kappa}$). Hence a technique to increase the Seebeck coefficient and simultaneously decrease the thermal conductivity is highly desired. Very recent studies Lee and Lee (2015); Goharshadia and Mahdizadeha (2015) have shown that doping of and impurities in graphene sheets will decrease the thermal conductivity. Pop $et\ al.$ Pop et al. (2012) have mentioned that any surplus residue from sample fabrication or any form of disorder will reduce the thermal conductivity further. There have been many experimental Kan et al. (2008); Ding et al. (2009); Dutta et al. (2009); Pruneda (2010); Bhowmick et al. (2011); Liu et al. (2011) and theoretical Bernardi et al. (2012); Kumar et al. (2013); D’Souza and Mukherjee (2015) reports to show the formation of a band gap by doping graphene with boron nitride. Similarly there are reports to show the formation of band gaps in bilayer graphene Min et al. (2007); Castro et al. (2007); Mak et al. (2009) when under the influence of an electric field. We have successfully used the Boltzmann transport equation for electron transport to calculate the thermoelectric parameter of Cx(BN)1-x D’Souza and Mukherjee (2015) and graphene/h-BN/graphene heterostructures (D’Souza and Mukherjee, 2016a). Therefore in this section, using the Boltzmann transport equations, we study the behavior of these transport parameters focusing mainly on enhancing the Seebeck coefficient.

In Fig. 3 we show the calculated electrical conductivity ($\sigma$) and the Seebeck coefficient ($S$) for MLG upon doping by BN. To study the behaviour of doped graphene, we have substituted one and two dimers of boron and nitrogen in our graphene unit cell as shown in Fig. 3.

Since boron is an acceptor and nitrogen is a donor, the total number of charge carriers remains unchanged, leading to a gap at the Fermi energy Mukherjee and Kaloni (2012); D’Souza and Mukherjee (2015). This band gap transforms the metal to a semiconductor, thereby decreasing the electrical conductivity as shown in Fig. 3 (a) and (c). The $\frac{1}{\sigma}$ term in Eq. 4 increases the Seebeck coefficient. This is evident in Fig. 3(b). We can therefore predict that doping graphene with boron and nitrogen will increase the Seebeck coefficient.

Since a band gap decreases the electrical conductivity thus increasing the Seebeck coefficient, we apply an electric field perpendicularly to the monolayer graphene sheets. BLG in an electric field has been shown to have a band gap Min et al. (2007). We have thus considered the effect of an electric field for three different values of external potential ($U$), i.e., $U$=0.2, 0.3, 0.5 eV. The averaged Coulomb potential plotted as a function of its perpendicular length ($z$), is shown in Fig. 4(a).

In Fig. 4(b) we plot the Seebeck coefficient of BLG under the influence of an electric field, showing that an increase in the external potential results in an increase in $S$.

Experimentally the effect of electric field on $S$ for BLG has been studied by Wang et al. Wang et al. (2011), showing enhancement of $S$ with increasing electric field as our calculations indicate.

III.4 Impurity scattering in graphene

Increasing the Seebeck ($S$) coefficient of low-dimensional materials such as graphene has always been a pursuit in thermoelectric applications. Enhancement of $S$ by inelastic scattering has been reported recently by Ghahari $et\ al$ Ghahari et al. (2016). Experimental electrical conductivity studies of graphene doped with potassium as a function of charge carrier has been reported by Chen $et\ al$ Chen et al. (2008). From Mott’s formula Eq. 4, it can be readily seen that a decrease in electrical conductivity will enhance the Seebeck coefficient. As pointed in the previous section doping of boron nitride decreases the electrical conductivity and hence increases $S$. The unit cells used in those calculations are relatively small to mimic the experimental behavior of doping since there is no long-range nature of charge-impurity scattering. In order to understand the behavior of electrical conductivity as reported by Chen $et\ al$ Chen et al. (2008) one would require a long range nature of impurity scattering and hence a large simulation cell making first-principles (DFT) calculations extremely hard.

As only the $\pi$ states are responsible for transport in MLG and BLG D’Souza and Mukherjee (2016b), a tight-binding band calculation would be more useful since it would allow incorporation of a very large unit cell to account for the long-range nature of the impurity scattering. In this section we report the calculation of electrical conductivity of graphene with impurities with the help of our tight-binding method using the simple orthogonal nearest-neighbor tight-binding model which has the potential of modeling several impurity properties as described by Pedersen $et\ al$ Pedersen and Pedersen (2013).

In this model, we consider an infinite graphene sheet but with one or more atoms replaced by an impurity representing the graphitic impurities. Since we are interested only in the transport properties, our Hamiltonian states would consists only of the $p_{z}$ $\pi$ states. The Hamiltonian would then be written as

[TABLE]

with $H_{0}$ defined as

[TABLE]

where $\epsilon$ is the on-site energy of carbon, and $t_{ij}$ is the hopping integral which is described by $H_{0}$ between the $p_{z}$ $\pi$ states on site $i$ and $j$. The added $H_{imp}$ is the impurity Hamiltonian. It depends on which site an atom has been replaced by an impurity and takes the form

[TABLE]

where $\Delta$ is the increase or decrease of the on site energy on site $i$. It must be noted that the difference between $H_{imp}$ and the first term of $H_{0}$ is that the summation in $H_{imp}$ runs over those sites where an atom has been replaced by an impurity whereas in the first term of $H_{0}$ it runs over all available sites.

From Fig. 5 it can be seen that the behavior of electrical conductivity when plotted against charge carrier is that which is observed experimentally Chen et al. (2008) for doped systems. Larger values of $\Delta$ for the same number of atoms replaced by an impurity decreases the electrical conductivity and hence from Eq. 4 would increase the Seebeck coefficient. One can also see the behavior for electrical conductivity on pristine graphene ($\sigma\sim\sqrt{n}$) tending to the linear behaviour of electrical conductivity ($\sigma\sim n$) for graphene sheets with impurities, a fact that has be observed experimentally Bolotin et al. (2008).

III.5 Mobility ($\mu_{FE}$) of doped and undoped MLG

In Fig.6(a,b) we show the scaled mobility ($\mu_{FE}/\tau$) of MLG as function of energy and carrier concentration $n$, calculated from the Boltzmann transport equations, showing similar trends as to those seen in experiments by tuning the gate voltage and carrier concentration Ponomarenko et al. (2009); Ghahari et al. (2016), respectively. We have calculated $\mu_{FE}$ using,

[TABLE]

It is easy to understand the behavior of $\mu_{EF}$ since the derivative of electrical conductivity with respect to $n$ should be proportional to $\frac{1}{\sqrt{n}}$. This behavior has been reported by Ponomarenko $et\ al$. Ponomarenko et al. (2009). The method used to introduce impurities is discussed in detail in the previous section. Experimental data Tan et al. (2007); Chen et al. (2008) for graphene samples with increasing doping concentrations have been shown to reduce its mobility by an order of magnitude. This effect can be seen in Fig. 6(b), where variation of $\mu_{FE}/\tau$ is shown as a function of carrier concentration $n$ plotted in logarithmic scale. Our results on mobility show similar behavior to that observed experimentally Ponomarenko et al. (2009); Ghahari et al. (2016).

Expressing mobility in units of $\tilde{n}=\frac{n}{10^{10}}$ cm*-2* and $\tilde{\sigma}=\frac{\sigma}{h/e^{2}}$ k$\Omega^{-1}$, and assuming the electron relaxation time $\tau\sim 1\times 10^{-14}$s, our calculations result in a value of $\mu_{FE}\approx 1.6\times 10^{4}\frac{\tilde{\sigma}}{\tilde{n}}$ [cm2/Vs] which is close to the earlier estimate of $\mu_{FE}\approx 2.42\times 10^{4}\frac{\tilde{\sigma}}{\tilde{n}}$ [cm2/Vs] Sarma et al. (2011).

III.6 Phonon dispersion, Grüneisen parameter, and lattice thermal conductivity

III.6.1 Phonon dispersion

In Fig. 7 we show our calculated phonon band structure along the high-symmetric points in the irreducible hexagonal Brillouin zone (BZ) for the monolayer and bilayer graphene. Accurate calculation of phonon dispersion of MLG and BLG is necessary to understand the thermal conduction in these materials. Based on the harmonic second-order IFCs, we calculate the phonon dispersion of MLG and BLG along high-symmetric $q$ points obtained within the linear response framework by employing density functional perturbation theory (DFPT) Baroni et al. (1987), as implemented in the Quantum Espresso code Giannozzi and et al. (2009) described earlier.

The out-of-plane (ZA), in-plane longitudinal (LA), and in-plane transverse (TA) modes, which arise from the $\Gamma$ point of the BZ of MLG, correspond to the acoustic mode while the remaining branches correspond to the optical modes (ZO, LO and TO) de Juan et al. (2015). The TA and LA modes show linear $q$ dependence at low $q$, as is usually seen for acoustic modes. The out-of-plane ZA mode shows a quadratic ($q^{2}$) dependence, which is a distinctive feature of layered crystals as observed experimentally Mounet and Marzari (2005); Zabel (2001). An explanation of this quadratic dependence could be due to the two-dimensional out-of-plane phonon mode and threefold rotational symmetry for BLG (sixfold for MLG) Saito et al. (1998). The LO and TO modes are degenerate at $\Gamma$ having a frequency of 1580 cm*-1*. Our calculated value of the degenerate frequency is in good agreement with the result using inelastic x-ray scattering measurements by Maultzsch $et\ al.$ Maultzsch et al. (2004) having a value of 1587 cm*-1*. High-voltage transport measurements by Yao $et\ al$ Yao et al. (2000) estimated that for graphite the frequency of zone-boundary phonons should be around 1300 cm*-1*. Our calculations show that at $K$, the BZ corner, the phonon energy of the in-plane transverse optical (TO) mode has a frequency of 1370 cm*-1* for MLG and 1287 cm*-1* for BLG. This suggests that our calculations agree well with the experiment. The phonon dispersion of BLG is very similar to that of MLG except for a characteristic feature of an additional low-frequency optical mode with energy nearly about 108 cm*-1* at $\Gamma$. This layer breathing mode arises due to the interlayer movements. The phonon dispersions shown in Fig. 7, calculated using harmonic IFCs, are consistent with both experimental and previous theoretical studies Yan et al. (2008); Nika et al. (2009); Kong et al. (2009); Yao et al. (2000); Maultzsch et al. (2004).

III.6.2 Grüneisen parameter ($\gamma$)

To carry out a precise calculation of the lattice thermal conductivity, effects from the harmonic and anharmonic lattice displacements should be taken into account to include contributions of higher order phonon-phonon scattering processes Broido et al. (2005). Since the Grüniesen parameter ($\gamma$) provides useful information on the phonon relaxation time and the anharmonic interactions between lattice waves and the degree of phonon scattering, we have therefore calculated the mode-dependent Grüneisen parameter ($\gamma$) for MLG and BLG. We employ the method as developed previously Kong et al. (2009); Cai et al. (2014); Mounet and Marzari (2005) to calculate the degree of phonon scattering. It is carried out by dilating the lattice applying a biaxial strain of $\pm$ 0.5%. In Fig. 8 we show the calculated mode-dependent Grüneisen parameter along the high-symmetric $q$ points. $\gamma$ is expressed as

[TABLE]

where $a_{0}$ is the relaxed lattice constant without strain, $\omega$ is the phonon frequency, $\omega_{+}$ and $\omega_{-}$ are the phonon frequencies under positive and negative biaxial strain, respectively, and $da$ is the difference in the lattice constant when the system is under positive and negative biaxial strain.

MLG and BLG both have negative values for $\gamma$ along the high-symmetric $q$ points for the out-of-plane acoustic (ZA) and optical (ZO and ZO’) modes while $\gamma$ has only positive values for the in-plane longitudinal and transverse modes. Negative (positive) $\gamma$ implies an increase (decrease) in phonon frequency when the lattice constant is increased. The slight difference in the Grüneisen parameters in MLG and BLG is, in the case of BLG near the long-wavelength limit ($\Gamma$ point), $\gamma_{ZO}$ corresponding to the out-of-plane optical mode changing sign unlike in MLG. This suggests that near the long-wavelength limit the atom vibrations perpendicular to the plane of the sheet between the two layers lose their coherence hence decreasing the phonon frequencies of $\omega_{ZO}$ when under a biaxial strain. Since the TO mode in graphene is Raman active at the $\Gamma$ point, the Grüneisen parameter can be measured experimentally using Raman spectroscopy. Mohiuddin $et\ al$ Mohiuddin et al. (2009) have measured $\gamma_{E_{2g}}$ to be 1.99, which is in excellent agreement with our calculated value of $\gamma_{TO}=1.85$ ($\gamma_{TO}=1.89$) for MLG (BLG). In Table 1 we show the phonon frequencies and the Grüneisen parameters of MLG and BLG at high-symmetric $q$ points $\Gamma$, K, and M for different vibrational modes.

Our calculated $\omega$ for MLG are in very good agreement with previous calculations Mounet and Marzari (2005), and also $\gamma$ show excellent agreement with previous calculations Kong et al. (2009); Mounet and Marzari (2005) and experiment Mohiuddin et al. (2009). Our results for $\gamma$ of BLG are similar to that of graphite Mounet and Marzari (2005) but differs from BLG results by Kong et al Kong et al. (2009) in the low-$q$ region for the ZO’ and ZA modes. The phonon dispersion of graphite Mounet and Marzari (2005) is very similar to that of BLG, so a similar behavior in $\gamma$ can be expected. We feel that the discrepancy in the BLG results by Kong et al Kong et al. (2009) can be attributed to their unstable phonon frequencies for ZO’ and ZA modes near $\Gamma$ point, which may not be correct. It should be noted that the thermal conductivity can be calculated with the help of the mode-dependent Grüneisen parameter using the Callaway-Klemens approach Klemens (1951); Callaway (1959); Nika et al. (2009). However, we have not used this approach but a real-space supercell approach as implemented in the ShengBTE method Li et al. (2014a).

III.6.3 Lattice thermal conductivity ($\kappa_{L}$)

In the ShengBTE method Li et al. (2014a), the third-order anharmonic interatomic force constants (IFCs) were also taken into account apart from the usual second-order harmonic IFCs which produced the phonon dispersion, in the calculations of thermal conductivity ($\kappa_{L}$). The third-order anharmonic IFCs were calculated using a finite-difference supercell approach with a set of displaced supercell configurations depending on the size of the system. We have used a $4\times 4\times 1$ supercell for both MLG and BLG, which generated 72 and 156 configurations for MLG and BLG, respectively. The three-phonon scattering amplitudes are then computed from a set of third-order derivatives of energy, calculated from these configurations using the Quantum Espresso code Giannozzi and et al. (2009).

In Fig. 9 we show our calculated lattice thermal conductivity ($\kappa_{L}$) of MLG and BLG using Eq. 5 as implemented in the ShengBTE code Li et al. (2014a). In the inset we compare our results to experimental data of Li et al Li et al. (2014b) available in the temperature range 300K to 700K. Our results are in very good agreement with experimental measurements. Graphene at room temperature (RT) has one of the highest know $\kappa_{L}$.

The experimental results of $\kappa_{L}$ for MLG Chen et al. (2011, 2012); Balandin et al. (2008); Ghosh et al. (2008); Cai et al. (2010); Jauregui et al. (2010) have shown that for freely suspended samples $\kappa_{L}$ lies between 2000-5000 Wm*-1K-1*. This wide variation in experimental estimate of $\kappa_{L}$ is presumably due to any disorder or residue from fabrication leading to an increase in the phonon scattering. We have, therefore, taken the most recent data by Li $et\ al.$ Li et al. (2014b) to compare with our calculations. Our calculated $\kappa_{L}$ at RT of MLG and BLG were found to be 2870 Wm*-1K-1* and 1730 Wm*-1K-1*, respectively, which is within the range seen experimentally, and are in good agreement with the previous literature Nika et al. (2009); Li et al. (2014b); Pop et al. (2012); Chen et al. (2011, 2012); Ghosh et al. (2010). Our calculations also show that at higher temperatures, $\kappa_{L}$ does not change significantly by addition of another layer which is consistent with the report by Koh $et\ al.$ Koh et al. (2010), suggesting that $\kappa_{L}$ between graphene and its environment has a much larger influence than that of individual graphene sheets. We find that $\kappa_{L}$ increases initially from 20 K to 170 K for MLG and to 230 K for BLG, before decreasing. For MLG, if we compare our results to the experimental data by Chen $et\ al$ Chen et al. (2011), we find that the maximum values of $\kappa_{L}$ seen experimentally occuring at a temperature between 150 K and 200 K are in agreement to our calculations. Therefore, calculations involving both harmonic and anharmonic IFCs, solving the BTE for phonons as done in the ShengBTE method, provides an accurate method for the calculation of the lattice thermal conductivity.

IV Summary

We have reported various transport properties such as electrical conductivity, resistivity, the Seebeck coefficient, mobility and lattice thermal conductivity of MLG and BLG graphene using first-principles DFT calculations and Boltzmann transport equations. We were able to capture many essential features seen in MLG and BLG, for example $\textendash$ the $\sqrt{n}$ behavior of electrical conductivity and its temperature dependence, the increase of the Seebeck coefficient with temperature, and a linear dependence of the Seebeck coefficient on temperature for a constant chemical potential, as observed in experiments. For a particular range of chemical potentials we obtained the Bloch-Gr${\rm\ddot{u}}$neisen behavior of resistivity in MLG, where the resistivity increased linearly at higher temperatures whereas it showed a $\sim T^{4}$ behavior at lower temperatures, as observed experimentally. We have also observed an order of magnitude decrease in mobility when the energy on impurity sites is decreased, a fact that has been verified experimentally. The Seebeck coefficient was found to increase almost twofold upon doping by boron nitride. Our results for graphene with impurities show a systematic decrease in electrical conductivity (and hence mobility) when we decrease the on-site terms of particular atoms in the sheet. We also observe that for a high concentration of impurities, the electrical conductivity was found to change from a $\propto\sqrt{n}$ behavior to a $\propto n$ one.

Our calculated phonon dispersion and Grüneisen parameters with harmonic and anharmonic IFCs, respectively, for both MLG and BLG show good agreement with available experimental data and previously published calculations. We finally show the result of the lattice thermal conductivity, calculated using phonon Boltzmann transport theory and first-principles phonon bandstructure including both harmonic and anharmonic interactions, showing excellent agreement with recent experimental data Li et al. (2014b) available in the temperature range 300-700 K. Further experimental measurements are needed to verify the occurrence of a peak in $\kappa_{L}$ near $T\sim$150-200 K for both MLG and BLG.

V Acknowledgment

We thank Dr. Jesús Carrete for his insightful correspondence on the ShengBTE code for the calculation of thermal conductivity. All calculations were performed in the High Performance Cluster platform at the S.N. Bose National Centre for Basic Sciences. RD acknowledges support through a Senior Research Fellowship of the S.N. Bose National Centre for Basic Sciences.