Phonon-limited transport coefficients in extrinsic graphene

# Phonon-limited transport coefficients in extrinsic graphene

Enrique Muñoz Instituto de Física, Pontificia Universidad Católica de Valparaíso
PO Box 4059, Valparaíso, Chile.
September 14, 2019
###### Abstract

The effect of electron-phonon scattering processes over the thermoelectric properties of extrinsic graphene was studied. Electrical and thermal resistivity, as well as the thermopower, were calculated within the Bloch theory approximations. Analytical expressions for the different transport coefficients were obtained from a variational solution of the Boltzmann equation. The phonon-limited electrical resistivity shows a linear dependence at high temperatures, and follows at low temperatures, in agreement with experiments and theory previously reported in the literature. The phonon-limited thermal resistivity at low temperatures exhibits a dependence, and achieves a nearly constant value at high temperatures. The predicted Seebeck coefficient at very low temperatures is , which shows a dependence with the density of carriers, in agreement with experimental evidence. Our results suggest that thermoelectric properties can be controlled by adjusting the Bloch-Grüneisen temperature through its dependence on the extrinsic carrier density in graphene.

###### pacs:
65.80.Ck, 72.80.Vp, 63.22.Rc

Recently, Efetov and Kim Efetov and Kim (2010) reported experimental measurements of the phonon-limited electrical resistivity of graphene samples, by achieving extremely high carrier densities by means of an electrolytic gate, in order to minimize the effect of other scattering mechanisms and to be away from the neutrality (Dirac) point. They observed that, at low temperatures, the electrical resistivity displays a dependence, whereas at higher temperatures it follows a linear trend, in agreement with an earlier theory by Hwang and Das Sarma Hwang and Das Sarma (2008). The later based their analysis on the Bloch theory Ziman (1960); Madelung (1978), which neglects phonon drag effects and umklapp processes, and modeled the interaction between electrons and longitudinal acoustic phonons by the deformation potential approximation. Moreover, they calculated the electrical resistivity Hwang and Das Sarma (2008) within the relaxation time approximation, which is valid for elastic scattering processes, a condition not strictly satisfied in electron-phonon scattering.

In this paper, phonon-limited transport coefficients in graphene are studied, by considering the presence of temperature and voltage gradients. Under similar assumptions as in Ref. (Hwang and Das Sarma (2008)), that is, Bloch theory and deformation potential approximation for the interaction between electrons and longitudinal-acoustic phonons, analytical expressions for the transport coefficients are obtained from a variational solution of the Boltzmann equation Ziman (1960); Madelung (1978). This method, at the lowest-order approximation, provides equivalent results to the relaxation time approximation and, at higher orders, converges to solutions of the Boltzmann equation which are exact to an arbitrary precision, within the limitations of the set of variational functions chosen. As has been pointed out before Peres and Lopes (2007), the semiclassical approach based on the Boltzmann equation cannot capture the subtleties of electronic transport in graphene right at the neutrality point. This limitation can, at least in part, be attributed to the inhomogeneity in the charge distribution experimentally observed in graphene samples at the neutrality pointPeres (2010). However, at high carrier densities the Fermi level is well above the Dirac point, and the Boltzmann equation then provides a correct description.

## I Electron-phonon interaction

As is well known from tight-binding calculations, electronic transport in graphene displays relativistic features. This is because pristine graphene, free of impurities and defects, behaves as a metal only in the so-called Dirac points where the conduction and valence bands touch. There exists two non-equivalent such points within the first Brillouin zone, defined by the vectors , with the direct lattice parameter. In the vicinity of each of these points, or valleys, it is shown that the dispersion relation is linear, , for . In each valley , the envelope electronic functions related to this dispersion relation can thus be obtained as the pseudo-spinor solutions of an effective two-dimensional Dirac’s equation,

 ^HKψ(±)k=vF^α⋅^pψ(±)k=E(±)kψ(±)k. (1)

Here are the Pauli matrices, we have re-defined for each valley, and the energy eigenvalues for positive (negative) helicity are . The corresponding eigenvectors, normalized by the total area of the sample, are given by the pseudo-spinors

 ψ(±)k(r)=1√2A(e−iφk/2±eiφk/2)eik⋅r≡u(±)keik⋅r√A (2)

with .

In the language of second quantization, the conduction electrons field for in the vicinity of the valley is

 ^ψη(r)=∑k,σψ(+)k(r)^ckση (3)

with Fermionic operators , that destroy (create) a conduction electron with momentum in the vicinity of the valley, with spin component . We model the electron-phonon interaction by the deformation potential approximation Hwang and Das Sarma (2008), which is a reasonable one when the Fermi surface possesses spherical symmetry Ziman (1960); Madelung (1978). The two-dimensional version of this is a Fermi ”circle” in graphene, which corresponds to the intersection of the Dirac cone and the constant energy plane defined by the Fermi energy at . In the deformation potential approximation, the operator representing lattice dilation in the second quantization language is

 ^Δ(r)=i∑q(2Aρmωq)−1/2(q⋅^eL)(^a†qeiq⋅r−^aqe−iq⋅r) (4)

Here, is the frequency for longitudinal acoustic phonons in graphene, is the creation operator for longitudinal acoustic phonon modes, and is the corresponding polarization vector for these modes. The electron-phonon interaction Hamiltonian, in the deformation potential approximation, is obtained from

 ^Hinte−ph=i∫d2r∑η=1,2^ψ†η(r)D^Δ(r)^ψη(r) (5) =i∑k,q,η,σ√ℏq2D22ρmωqcos(θ/2)(^aq−^a†−q)^c†k+q,ση^ckση

Notice the presence of the factor which arises from the inner product of the two pseudo-spinor functions defined in Eq. (2). This factor, which would be absent from the inner product of non-relativistic scalar Schrödinger wave-functions, is a fingerprint of the relativistic features of electrons in graphene. The more obvious consequence of its presence the suppression of backscattering (), a phenomenon observed in graphene and directly related to so called Klein tunneling Beenakker (2008); Peres (2010); Castro, Guinea and Peres (2009) in the context of relativistic quantum mechanics.

We shall neglect umklapp processes, and thus we have for normal processes. In this case, scattering events satisfy energy and ”crystal momentum” conservation, and . From the latter equation, one obtains the condition , with the scattering angle. Since the relevant scattering events occur within a thin layer near the Fermi surface, then , which then implies . In a Debye model approximation, this condition restricts the range of lattice momenta involved in electron-phonon scattering, , with the radius of the Debye circle (for a two-dimensional system). In geometrical terms, the most restrictive condition depends on whether the radius of the Fermi surface (a circle centered on each Dirac point for graphene) is larger or smaller than the Debye circle. When written in terms of the frequency of the longitudinal-acoustic phonon modes involved , the former condition reads , with the Debye temperature, and the Bloch-Grüneisen temperature. In the experimental setup presented in Ref.(Efetov and Kim (2010)), the carrier densities involved are such that the Fermi surface is contained within the Debye circle, and therefore the Bloch-Grüneisen rather than the Debye temperature imposes the characteristic energy scale for phonon scatterers. Notice that the analysis of the scattering angle above also yields the relation . Therefore, after Eq. (I) the matrix element for the electron-phonon interaction is

 Mq=−i(ℏ/2ρmωq)1/2Dq(1−(q/2kF)2)1/2. (6)

Here, is the deformation potentialZiman (1960); Madelung (1978); Efetov and Kim (2010) coupling constant for graphene, and kg/ the surface mass density.

## Ii Boltzmann Equation

We introduce, in the usual form, the deviation of the distribution function from equilibrium by the expression . Let us consider the electron-phonon interaction as the sole scattering mechanism. Therefore, after Eqs.(5,6) we should consider processes of the form . Under Bloch’s approximation, which is to assume that the phonon system is in quasi-equilibrium, the linearized form of the Boltzmann equation when both an external electric field and a temperature gradient are imposed can be expressed as

 −vk⋅∇T∂f0k∂T−vk⋅E(−e)∂f0k∂ϵk= gηgσkBT∑k′,q[{χk−χk′}Pk′kq−{χk′−χk}Pkqk′] (7)

Here, are the spin and valley degeneracies, respectively, while the transition rate for electron-phonon scattering processes is

 Pk′kq =(2π/ℏ)δG,k′−k−q|Mq|2δ(ϵk−ϵk′+ℏωq) (8) ×n0qf0k(1−f0k′)

Here is the Fermi-Dirac distribution for electron states with momentum , and the Bose distribution for phonon states with momentum . All vectors are two-dimensional, with a vector of the reciprocal lattice.

For extrinsic graphene, the Fermi level is located above the Dirac point, and depends on the density of carriers by the relation . In Ref.(Efetov and Kim (2010)), an experimental method is presented which allows to control the carrier density, and hence the Bloch-Grüneisen temperature. As discussed in Refs.(Efetov and Kim (2010); Fuhrer (2010)), this in turn becomes a way to control the phonon-limited electrical resistivity. In this work, we will show that also the electronic component of the thermal conductivity, as well as the thermopower (Seebeck coefficient) which determines thermoelectric effects, can be controlled in a similar way. We shall use the variational method to calculate the phonon-limited transport coefficients which, in contrast with the relaxation time approximation, it has the advantage of avoiding the assumption of quasi-elasticity in the scattering processes.

## Iii Coupled charge and heat transport

For a system with thermal and electrostatic potential gradients, the macroscopic entropy production rate is given in terms of the heat and charge fluxes by Ziman (1960); Madelung (1978)

 ∂tS=J⋅E/T+U⋅∇(1/T). (9)

The Onsager’s theorem Onsager (1931) indicates that the fluxes and generalized potentials are then linearly coupled,

 J=LEEE+LET∇T U=LTEE+LTT∇T (10)

The thermopower can be expressed in terms of the coefficients by considering the case in Eq. (10), to obtain . Hence, one concludes that . The thermal conductivity is obtained from the second equation, , which implies . If, on the other hand, one considers the case , then the relation is obtained. Thus, it is concluded that the electrical resistivity is . The variational method Ziman (1960); Madelung (1978) is based on the principle of maximal entropy production Onsager (1931), and provides an iterative but virtually exact procedure to solve the Boltzmann equation, and therefore to calculate the transport coefficients, within the limitations of the set of variational functions chosen. The principle is based on equating the macroscopic entropy production rate, as expressed in terms of the macroscopic currents Eq.(10), with the entropy production rate due to microscopic scattering eventsOnsager (1931), , as defined from the linearized Boltzmann equation scattering terms,

 ˙Sscatt=gηgσkBT2∑k,q,k′{χk−χk′}2Pk′kq (11)

The electric current and the heat current are given by the expressions

 J = gηgσ∑k(−e)vkχk∂f0k∂ϵk U = gηgσ∑kvk(ϵk−EF)χk∂f0k∂ϵk (12)

The technique is based on expanding the distribution function deviation from equilibrium in terms of a set of variational functions of the form . For a minimal set of trial functions and , we obtain the charge and heat currents

 Ji =(−e)gηgσ∑kvkϕi(k)∂f0∂ϵk Ui =gηgσ∑kvk(ϵk−EF)ϕi(k)∂f0∂ϵk (13)

Here, are the valley and spin degeneracies in graphene, respectively. The coefficients are defined as follows

 Pij=(gηgσ)2kBT∑k,k′,q{ϕi(k)−ϕi(k′)}{ϕj(k)−ϕj(k′)}Pk′kq

In terms of the variational solution, after Eq.(12) the electric and heat currents are written as

 J = ∑iαiJi U = ∑iαiUi (15)

whereas the macroscopic and microscopic entropy generation rates in Eqs.(9,11) are expressed by

 ∂tS = ∑iαi[T−1Ji⋅E−T−2Ui⋅∇T], ˙Sscatt = T−1∑i,jαiαjPij. (16)

The variational problem is to find the set of coefficients that maximize , subject to the constraint , which is enforced by a Lagrange multiplier ,

 δδαn[∑i,jαiαjPijT−λ∑iαi(Ji⋅ET−Ui⋅∇TT2)]=0

The solution to this constrained variational problem is given by , and

 αi=∑j[P−1]ij(Jj⋅E−Uj⋅∇TT) (18)

Here, is the inverse of the matrix whose elements are the coefficients defined by Eq.(LABEL:eq14). After substituting this solution into the expressions for the electric and heat currents, we obtain explicit analytical formulas for the transport coefficient tensors defined in Eq.(10),

 LEE = ∑ijJi[P−1]ijJj LET = −T−1∑ijJi[P−1]ijUj LTE = ∑ijUi[P−1]ijJj LTT = −T−1∑ijUi[P−1]ijUj (19)

Therefore, we proceed to calculate the currents and coefficients. Let us first consider the current . The group velocity in Eq.(13) is given by , with . The discrete sums are treated in the usual quasi-continuum limit as two-dimensional integrals, , and the integral is calculated by using the change of variables ,

 J1 =A(−e)vFπ2(ℏvF)3∫+π−πdφ^k^k⋅^u∫dϵϵ2∂f0ϵ∂ϵ (20) =^uAeπℏ(EFℏvF)2

For the current , we consider the quasi-continuum limit for the sum over wave-vectors in Eq.(13), with directions . Using the integration variables , we have

 U1 = AvFπ2(ℏvF)3∫π−πdφ^k^k⋅^u∫dϵϵ2(ϵ−EF)(∂f0ϵ∂ϵ) (21) = −^u2πA3(kBT)2ℏEF(ℏvF)2

Similarly, for the current we obtain

 U2 = AvFπ2(ℏvF)3∫π−πdφ^k^k⋅^u∫dϵϵ2(ϵ−EF)2(∂f0ϵ∂ϵ) (22) = −^uπA3(kBT)2ℏ(EFℏvF)2

Finally, it is straightforward to verify by inspection that .

Let us now consider the coefficient . From Eq.(8), we first notice that the transition probability rate imposes the condition . We define the directions and , with the scattering angle. As before, in the quasi-continuum limit for the sums over wave-vectors, we choose the integration variables , to obtain

 P11 =(kBT)−1A2(πℏvF)4∫dϵϵ∫dϵ′ϵ′ (23) ×∫π−πdθ∫π−πdφ′(q⋅^u)2Pk′kq

In Appendix A it is shown that the angular integral , with . Hence, substituting the explicit expression for the transition probability rate, Eq.(23) reduces to

 2(kBT)−1A2π2ℏ(ℏvF)4∫π−πdθq2n0q|Mq|2 ×∫dϵϵ(ϵ+ℏωq)f0ϵ(1−f0ϵ+ℏωq) (24)

The energy integral is evaluated using the technique presented in Appendix B. Using the change of variables , we obtain

 P11=(kBT)−1A2kFℏD2π2ρm(ℏvF)2∫2kF0dq(q4+4(vsvF)2 ×(TΘBG)2[π2q43+q62][1−(q2kF)2]1/2) ×(1−e−ℏωq/(kBT))−1(eℏωq/(kBT)−1)−1 (25)

Eq.(25) is written in terms of the functions defined in Appendix C,

 P11=ℏD2A2EF(kBT)4π2ρm(ℏvs)5(ℏvF)3[J4(T/ΘBG)+4(vsvF)2 ×(TΘBG)2{π23J4(T/ΘBG)+12J6(T/ΘBG)}] (26)

Let us now consider the coefficient , which after similar manipulations as before, can be expressed by the integral form

 P22 =(kBT)−1A2(πℏvF)4∫π−πdθ∫π−πdφ′∫dϵϵ∫dϵ′ϵ′ (27) ×[(q⋅^u)(ϵ−EF)+ℏωq(k′⋅^u)]2Pk′kq

After expanding the square in Eq.(27), it is shown in Appendix A that the angular integrals give, , , and . Then, after similar algebraic procedures as described previously, Eq.(27) reduces to

 P22=ℏD2A2E3F(kBT)4π2ρm(ℏvs)3(ℏvF)5[{1+4π23(1+12 ×(vsvF)2)}(TΘBG)2J4(T/ΘBG)−23{1−6 ×(vsvF)2}(TΘBG)2J6(T/ΘBG)+16π23(vsvF)2 ×(TΘBG)4{(vsvF)2J6(T/ΘBG)+110π2{1 −2(vsvF)2}J8(T/ΘBG)}] (28)

At last, we consider the coefficient . After our definition Eq.(LABEL:eq14), we have

 P12 =(kBT)−1A2(πℏvF)4∫π−πdθ∫π−πdφ′∫dϵϵ∫dϵ′ϵ′ (29) ×[(q⋅^u)2(ϵ−EF)+ℏωq(k′⋅^u)(q⋅^u)]Pk′kq

As shown in Appendix A, the angular integrals are given by , and . Therefore, using the procedure described previously, the integral in Eq.(29) is calculated and expressed in terms of the functions defined in Appendix C,

 P12=ℏD2A2E2F(kBT)4π2ρm(ℏvs)3(ℏvF)5[J4(T/ΘBG)−vsvF(TΘBG) ×J5(T/ΘBG)+8π23(TΘBG)2{J4(T/ΘBG)+32 ×(vsvF)2J4(T/ΘBG)+14π2J6(T/ΘBG)}+23 (30) ×(vsvF)3(TΘBG)3{J7(T/ΘBG)−2π2J5(T/ΘBG)}]

## Iv Electrical Resistivity

As discussed in Section III, the electrical resistivity is obtained by setting in Eq. (10), to obtain . By direct substitution of the coefficients and currents obtained in Section III, we obtain after Eqs.(19) that the leading term which defines the electrical resistivity due to electron-phonon scattering processes is

 ρe−ph(T)=ρ0[(TΘBG)4J4(ΘBG/T)+4π23(vsvF)2 ×(TΘBG)6{J4(ΘBG/T)+32π2J6(ΘBG/T)}] (31)

Here, is a coefficient with dimensions of a two-dimensional resistivity (), and is the Bloch-Grüneisen temperature. The functions and their asymptotic properties are defined in Appendix C. Eq. (31) depends on the dimensionless parameter , a feature previously obtainedHwang and Das Sarma (2008) within the relaxation time approximation. The functions , defined in Appendix C, have the property , with the Riemann zeta function. Therefore, the low-temperature () behavior of the electron-phonon contribution to the electrical resistivity is

 ρe−ph(T) =ρ04!ζ(4)(T/ΘBG)4+o(T/ΘBG)6 (32)

At high temperature (), from Eq.(31) we find the asymptotic limit . This behavior is in agreement with experiments Efetov and Kim (2010) and earlier theoretical results based on the relaxation time approximation Hwang and Das Sarma (2008), thus supporting the application of the variational method in this case. In Fig.(1) we compare the theoretical prediction with the experimental data reported in Ref. (Efetov and Kim (2010)). The experimental data are fitted to the expression

 ρ(T)=ρ(T=0)+ρe−ph(T), (33)

with defined by our theory Eq.(31), and being the value of the resistivity at zero temperature. This parameter can be interpreted as the (temperature independent) contribution of impurity scattering to the total electrical resistivity. The impurity contribution to the electrical resistivity has been discussed Stauber, Peres and Guinea (2007); Peres and Lopes (2007), and it is a sample-dependent property, since it is proportional to the surface concentration of impurities . It has been shown Stauber, Peres and Guinea (2007); Peres and Lopes (2007) that, for a short range scattering potential , the in-plane resistivity is , whereas for a long range Coulomb potential Stauber, Peres and Guinea (2007) , this last case showing an inverse dependency on the carrier concentration as well. Therefore, it is expected that the zero temperature resistivity exhibits a dependence on the carrier concentration of the form , with the value of the coefficients and depending on details of the sample, particularly on the concentration and distribution of impurities and defects. We verify that this equation is in excellent agreement with the data reported in Ref. (Efetov and Kim (2010)), as shown in the inset of Fig. (1). We extracted values () and (), with linear regression coefficient .

As seen in Fig. (1), Eq.(33) provides an excellent fit to the experimental data. The temperature-dependent electron-phonon contribution to the electrical resistivity is fitted with just two free parameters: The energy parameter eV in the deformation potential, and the speed of sound for longitudinal phonons in graphene Km/s. It is remarkable that both values are in excellent agreement with independent estimations reported in the literature, particularly the speed of sound Munoz and Lu (2010); Saito Dresselhaus (1998). It is also relevant to notice that the ratio is negligibly small, and hence terms proportional to this factor can in practice be neglected in the analytical expressions obtained for the electrical resistivity, as well as in other transport coefficients discussed along this work. This also explains the quantitative success when applying the relaxation time approximation to electron-phonon scattering as reported in previous theoretical work Hwang and Das Sarma (2008), even though this scattering mechanism is not strictly quasi-elastic.

## V Thermal Resistivity

The thermal conductivity is obtained from the general expression Eq. (10) by setting , which as discussed in Section III leads to the expression . By direct substitution of the coefficients and currents calculated in Section III, and neglecting terms proportional to , we obtain that the leading contribution to the thermal resistivity due to electron-phonon scattering is given by

 ×J4(ΘBG/T)−12π2(T/ΘBG)4J6(ΘBG/T)} (34)

Here, is the Lorenz number for the free electron gas. At very low temperatures (), Eq. (34) predicts the asymptotic limit

 [κe−ph(T)]−1=ρ0L0ΘBG34π24!ζ(4)TΘBG+o(TΘBG)3 (35)

It is remarkable that, according to Eq. (35), the phonon-limited thermal resistivity is linear at very low temperatures, in contrast with the typical behavior theoretically predicted and experimentally observed in normal, three-dimensional metals Ziman (1960); Madelung (1978). The thermal resistivity contribution due to elastic scattering with Coulomb impurities has been calculatedPeres and Lopes (2007), and corresponds to , with . The total thermal resistivity can be estimated from the expression . In Fig.(2) inset, we compare the value of the total three-dimensional thermal conductivity at room temperature, normalized by a nominal packing ”thickness” of Å for the graphene layer. For numerical evaluation, we assumed for the impurities, and from the experimental values of the zero temperature electrical resistivity, we extracted an average impurity concentration . For the experimental system reported in Ref.(Efetov and Kim (2010)), we estimated a dielectric constant representing the average between the substrate and the PEO polymer electrolyte. As seen in the inset of Fig.(2), the total normalized thermal conductivity at room temperature is on the order of . We can compare this with the thermal conductivity due to the phonon system, where the experimental valueGhosh and Calizo (2008); Balandin and Ghosh (2008), in excellent agreement with a theory previously reported by the present authorMunoz and Lu (2010), is aboutMunoz and Lu (2010); Ghosh and Calizo (2008); Balandin and Ghosh (2008) at room temperature. These results suggest that, in most of the temperature range, the phonon contribution to the thermal conductivity in graphene dominates over the electronic contribution, as the present author has pointed out elsewhere Munoz and Lu (2010).

## Vi Thermopower

As discussed in Section III, by setting in Eq.(10), we obtain that the thermopower (Seebeck coefficient) is given by . By direct substitution of Eq.(19), we obtain that the leading contribution to the thermopower is given by the expression

 Q(T)=−π23ek2BTEF[1+4π23(TΘBG)2]J4(ΘBG/T)+vsvF(TΘBG)J5(T/ΘBG)−2(TΘBG)2J6(ΘBG/T)[1+4π23(TΘBG)2]J4(ΘBG/T)−23(TΘBG)2J6(T/ΘBG) (36)

Therefore, we find that the thermopower at very low temperatures () becomes

 Q(T)∼−π23ek2BTEF[1+vsvF5!ζ(5)4!ζ(4)(TΘBG)]+o(T3)