Conductivity tensor of graphene dominated by spin-orbit coupling scatterers: A comparison between the results from Kubo and Boltzmann transport theories

The diagonal and Hall conductivities of graphene arising from the spin-orbit coupling impurity scattering are theoretically studied. Based on the continuous model, i.e. the massless Dirac equation, we derive analytical expressions of the conductivity tensor from both the Kubo and Boltzmann transport theories. By performing numerical calculations, we find that the Kubo quantum transport result of the diagonal conductivity within the self-consistent Born approximation exhibits an insulating gap around the Dirac point. And in this gap a well-defined quantized spin Hall plateau occurs. This indicates the realization of the quantum spin Hall state of graphene driven by the spin-orbit coupling impurities. In contrast, the semi-classical Boltzmann theory fails to predict such a topological insulating phase. The Boltzmann diagonal conductivity is nonzero even in the insulating gap, in which the Boltzmann spin Hall conductivity does not exhibit any quantized plateau.


Results
The Model Hamiltonian and the Green Function. In pristine graphene, the low-energy electron in K or K′ valley can be well described by a massless Dirac Hamiltonian 2,3 , which takes a form as , where A/B represents the sublattice and ↑ /↓ represents the electronic spin. According to relevant literature, the electronic SOC interaction in graphene can be incorporated into such a low-energy theoretical framework by adding a mass term to the above Dirac Hamiltonian 26 . Therefore, to model N i randomly distributed SOC scatterers in graphene, we can assume that the scattering potential takes a form as It gives a Dirac mass in the electronic spin and sublattice space and takes the form of delta function in the coordinate space. λ stands for the SOC potential strength. R j denotes the central position of an individual SOC impurity. It captures the SOC effect of the randomly distributed clusters of heavy adatoms on graphene sheet in an actual experimental situation. Noting that the delta function potential is only appropriate to the limiting case of the cluster size being far smaller than the electronic Fermi wavelength. It becomes invalid as the cluster size gets very large. In a recent literature 13 , it was found that large cluster size can destroy the QSH state of graphene. To study such an effect, we need to choose other function forms to mimic the finite range scattering potential in real space, e.g. the Gaussian function. It is an interesting issue, and could be studied in a future work.
Next we will introduce the Green function and the corresponding self-energies, in terms of which the DOS and the conductivity tensor of graphene in the presence of SOC scatterers are formulated. Prior to proceed, one can notice from the above Hamiltonian that the different valley and spin states are actually decoupled. Therefore, we can restrict our theoretical treatment in the K valley and spin-up subspace. The results in other subspaces can be derived similarly. For a spin-up electron in K valley, the unperturbed Hamiltonian is x y x y 0 Scientific RepoRts | 6:23762 | DOI: 10.1038/srep23762 and the scattering potential reduces to The eigenstate of the unperturbed Hamiltonian eq. (3) can be easily obtained, it is given by k r i i where s = + 1(s = − 1) denotes the conduction (valence) band, L is the linear size of the graphene sheet and ϕ is the angle between k and positive x axis, i.e. tan ϕ = k y /k x . The corresponding eigen energy is ε γ In the eigen representation of  H 0 , the matrix element of the Green function of the system is defined as For evaluating the macroscopically observable physical quantities, such as the conductivity tensor, we need to perform an average of the Green functions over all possible SOC impurity configurations, i.e.
Such an averaged Green function is connected to the proper self-energy Σ sk,s′k' (ε) by Dyson equation where G sk 0 is the unperturbed Green function Within SCBA, the self-energy in eq. (7) obeys an equation as follows  which can also be depicted by the Feynman diagrams shown in Fig. 1(a). In eq. (9),  means the impurity configuration averaging to the quantity in it, in the same way as explained above to the Green function. As a result of the impurity configuration averaging, the Green function 〈 G sk,s′k′ (ε)〉 is diagonal and isotropic about k, which can be seen by expanding the second diagram in Fig. 1(a) in terms of diagrams including only unperturbed Green functions, as shown in Fig. 1(b). Thus, we can use the shorthand G ss′k (ε) ≡ 〈 G sk,s′k′ (ε)〉 . Similarly, the self-energy is also k-diagonal. Moreover, it is independent of k. Up to now, the self-energy reduces to where n i = N i /L 2 is the concentration of the SOC impurities. To be more specific, we can define the even and odd self-energies as where the argument ε of the self-energies is dropped for simplicity. With the help of Dyson eq. (7), we can express the Green functions above in terms of the self-energies. Consequently, we get the self-consistent equations about the self-energies. They are More details about the derivation of eqs (13) and (14) are given in the Methods Section. We can take advantage of eqs (13) and (14) to evaluate the self-energies by numerical iteration method. When evaluating the integral I 0 numerically, we need to choose a cut-off upper limit k c . The choice of this cut-off wavevector is somewhat arbitrary, but the results depend only weakly on k c when it is large enough. When formulating the Hall conductivity given below, the derivatives of the self-energies with respect to the energy argument ε are involved. Differentiating eqs (13) and (14), we can get the equations about the derivatives of the self-energies where the shorthands ∂ Σ E/O ≡ (∂ Σ E/O /∂ ε) are adopted. I 1 is another k-integral similar to I 0 , which is given by In general, the self-energies are complex functions which depend on a complex energy argument. In particular, we can denote the retarded and advanced self-energies as In the weak scattering limit, namely n i λ → 0, one can readily find that Σ E → 0 and Σ O → − n i λ. Accordingly, the DOS has the following approximation This implies that within SCBA, the system has an energy gap which is proportional to the concentration and the strength of the SOC scatterers to the weak scattering limit. This gap is expected to be topologically nontrivial because it is opened by the SOC impurity. A more detailed derivation of eq. (21) is also given in the Methods Section.
The Conductivity Tensor from Kubo Formalism. To work out the expressions of the conductivity tensor, i.e. the diagonal and Hall conductivities, we start from the Kubo-Bastin formula for noninteracting electrons 32 where α or β denotes an arbitrary Cartesian coordinate x or y, and G ± = G(ε ± i0) is the retarded/advanced Green function. The delta functions can be expressed in terms of Green functions, i.e. δ(ε − H) = − π −1 Im G(ε + i0). Following the steps in ref. 33, we can simplify the expressions of diagonal and Hall conductivities at zero temperature. As a result, the diagonal conductivity depends only on the Fermi energy ε F and can be written as which is equivalent to Kubo-Greenwood formula 34 at zero temperature. However, the Hall conductivity is divided into two terms. One term depends only on the Fermi energy, whereas the other term is an integral up to the Fermi energy. Therefore, the Hall conductivity has the form which is equivalent to Kubo-Středa formula 35 . In eqs (23) and (24), the correlation function J αβ (ε, ε′ ) is defined as The velocity operator v α can be represented by Heisenberg In general, J αβ is a complex function of complex energy arguments ε and ε′ . Within SCBA, the correlation function J αβ (ε, ε′ ) can be obtained by summing the Feynman diagrams including the vertex corrections, as shown in Fig. 2. The summation turns out to be the sum of a geometrical series, and the results for J xx and J xy are where A = n i λ 2 /4L 2 , and the function φ is Substituting eqs (26)(27)(28)(29) into eqs (23) and (24) and taking advantage of the self-consistent eqs (13)(14)(15)(16)(17), we can express the diagonal and Hall conductivities in terms of the self-energies and the derivatives of them. Therefore, the conductivity tensor can be calculated directly from the self-energies and the derivatives of them, which can be calculated by self-consistent iteration method. For more details about the derivation of the conductivity tensor, see the Methods Section. We only give the result of the diagonal conductivity to the weak scattering limit here This result can also be obtained via the semi-classical Boltzmann transport theory, as can be seen in the following section.
The Conductivity Tensor from Boltzmann Theory. In this section, we derive the expression of diagonal and Hall conductivities from Boltzmann transport theory. We start from the semi-classical conductivity formula where E is an external electric field. W s′k′,sk is the scattering rate from |s, k〉 to |s′ , k′ 〉 and can be written in terms of the scattering matrix element as follows The summation of the unperturbed Green functions yields where ε c = γk c . Substituting eqs (34) and (35) into eq. (33), and dropping off the terms which are of the second or higher order of λ, we get the scattering rate which is the same as the weak scattering limit of the Kubo diagonal conductivity given by eq. (30). However, to get a nonzero Hall conductivity, we must include the second term of the inverse relaxation time, which corresponds to skew scattering 20 So much for the analytical results. We will discuss the numerical results in the next section.

Discussion
With the formulation of the conductivity tensor developed in the preceding section, we are in a position to perform a numerical calculation of the conductivity spectrums, i.e. the diagonal and the Hall conductivities as functions of the Fermi energy. And then a comparison between the results from Kubo and Boltzmann theories can be made, based on which the validity of the semi-classical Boltzmann transport theory to describe the conductivity tensor of graphene arising from SOC impurities can be clarified. However, before doing this, we would like to present numerical results about the DOS spectrums. From the approximate result of DOS to the weak scattering limit eq. (21), we know that the SOC impurity can open a band gap which is proportional to the product of the scatterer strength and concentration. The numerical results of DOS spectrums are shown in Fig. 3 for different strengths and concentrations of the SOC impurities. We can clearly see that a gap occurs around the Dirac point (zero-energy point) for all the cases. And at relatively weak strength and concentration of the SOC impurities, the approximate expression of DOS to the weak scattering limit, i.e. eq. (21), can give satisfactory results, well agreeing with the numerical results within SCBA. In the insights of Fig. 3, we give the dependence of the band gap on the strength and concentration of the SOC impurities. As given by eq. (21), the gap is 2n i λ in the weak scattering limit. The numerical results shown in the insights of Fig. 3 indicate that such a simple relation holds true only in the case of relatively weak scattering. However, some previous works reported that the simple linear dependence of the band gap on the product of the SOC scatterer strength and concentration, i.e. Δ g = 2n i λ, is still a good approximation even when the system is far away from the weak scattering limit 13,22,36 . The deviation of the band gap from the linear relation as shown in the insights of Fig. 3 is due to that the SCBA employed in this work excludes some high-order scattering processes. To illustrate this, we can go one step further to calculate the DOS under t-matrix approximation instead of SCBA. Within t-matrix approximation, we add up all the terms represented by the Feynman diagrams shown in Fig. 1(c) to obtain the self-energy self-consistently. The results are shown in Fig. 4(a). We can see that when the impurity strength is relatively weak (λ < 0.4), the t-matrix approximation and the SCBA almost give the same result. The band gap as a function of the SOC impurity strength obtained by t-matrix approximation is closer to the linear rule than that from SCBA, as seen in the insight of Fig. 4(a). Another issue we would like to study is the effect of usual scalar scatterers on the topologically nontrivial gap. To derive the SCBA self-energy in this case, we need to evaluate the Feynman diagrams shown in Fig. 1(d).
We assume that the strength of the scalar scatterers has a zero value on average. Thus, we characterize the disorder by the root mean square strength µ µ = rms 2 . The concentration of the scalar scatterers is denoted by n s . The DOS results are shown in Fig. 4(b). It can be readily seen that the band gap opened by the SOC scatterers survive a weak scalar disorder. However, a strong enough scalar disorder will close the gap, and hence destroy the QSH state. The derivation of the self-energy within t-matrix approximation and the SCBA self-energy in the presence of usual scalar scatterers are given in the Methods section. Now, we turn to discuss the numerical results about the conductivity tensor. The diagonal conductivity spectrums are shown in Fig. 5 for different SOC scatterer strengths and concentrations. From this figure, we can see that the result from Boltzmann theory is almost a nonzero constant even in the gap. It comes to the conclusion that the semi-classical Boltzmann theory fails to describe the electronic transport properties of graphene arising from the SOC scattering not only in the vicinity of the Dirac point, but also around the band edges. Notice that when the band gap is large, the band edges are far away from the Dirac point. As shown in Fig. 5(a,b), σ xx from Kubo formula vanishes in the SOC band gap. It increases gradually and tends to the Boltzmann result as the Fermi energy goes away from the gap. This result tells us that the Boltzmann theory is only valid in the short electronic wavelength regions where the quantum interference effect becomes weak. Before ending the discussion about the diagonal conductivity spectrum, we would like to point out that the electronic contributions to the diagonal conductivity from other subband spaces, e.g. the spin-down electrons or K′ valley electrons, are exactly the same as shown in Fig. 5.
In Fig. 6, the Hall conductivity spectrums obtained by Kubo formula are shown for different SOC scatterer strengths and concentrations. From this figure, we can see that for the K valley spin-up electron, the Hall conductivity shows a plateau with a height of e 2 /2h within the gap, regardless of the scattering strength and concentration. This result indicates that the gap opened by SOC impurities is topologically nontrivial. Considering that the spin-down electron contributes an opposite Hall conductivity due to the time reversal symmetry, the total Hall conductivity vanishes. However, the system is on a QSH state with a quantized spin Hall conductivity The results shown in Fig. 6 support our previous theoretical prediction that the randomly distributed SOC impurities can drive the graphene into a QSH state 22 . This quantized spin Hall conductivity has been obtained in ref. 23 where the SOC interaction does not act as a scatterer. Now the interesting thing is that the QSH state emerges even if the SOC is in the scatterer but not in the band. From Fig. 6, we can also see that outside the gap, the Hall conductivity is nonzero, though it decreases rapidly as the Fermi energy goes away from the gap. This indicates that in this region, the system is in a spin Hall regime even though the spin Hall conductivity is not quantized. The nonzero spin Hall conductivities outside the gap are also discussed in ref. 23, in which the intrinsic Hall conductivity agrees with our results. A typical numerical result of σ xy from Boltzmann transport theory is shown in Fig. 7. We can see that the Hall conductivity here is completely different from that calculated from Kubo theory. We cannot observe a spin Hall plateau in the gap region and the spin Hall conductivity tends to zero when ε F → 0. Although the spin Hall conductivity from Boltzmann theory is not quantized, it  has nonzero values especially when the Fermi energy is far away from the Dirac point. This nonzero spin Hall conductivity is the result of skew scattering, which is detailed studied in ref. 21.
In summary, a theoretical study on the conductivity tensor of graphene arising from the SOC impurities is presented in this work. The calculated results of the conductivity tensor treated by the semi-classical Boltzmann and the quantum transport approaches are compared. In the quantum transport approach realized by the Kubo-Středa formula within SCBA, the diagonal conductivity shows an insulating gap around the Dirac point. Meanwhile, in such a gap the spin Hall conductivity shows a well-defined quantized plateau. These features indicate unambiguously the realization of the topologically nontrivial state of graphene driven by the randomly distributed SOC impurities. In contrast, within the Boltzmann theory, the diagonal conductivity does not vanish and the spin Hall conductivity does not show any quantized plateau in the vicinity of the Dirac point. Thus, Boltzmann theory cannot well describe the low-energy transports in graphene arising from the SOC scatterers. Finally, we have to point out that recent experimental works on the transport properties of graphene dominated by heavy adatoms found no evidence of a SOC band gap 27,28 . The possible reason might be the clustering effects and the long range Coulomb scattering. Therefore, it is useful and interesting to perform a further theoretical study on the effects of the cluster size and the Coulomb scattering potential on the QSH state.

Methods
The Derivation of the Self-Energies. In this subsection, we give the details for the self-consistent equations of the self-energies. We start from the self-energy within SCBA , (1) , , which can be represented by the Feynman diagrams in Fig. 1(a). The first term is Because the average Green function is diagonal and isotropic about k, the shorthand notation G ss′k is adopted. Therefore, the second term of the self-energy yields We can see from eqs (41) and (43) that the self-energy is diagonal about k and does not depend on k. Therefore, the self-energy has the form of eq. (10). Moreover, there are only two different self-energies, the even self-energy Σ E and the odd self-energy Σ O , which are defined in eqs (11) and (12). By Dyson eq. (7), we can express the Green functions in terms of the self-energies, for example G −+k can also be evaluated by Dyson equation Combining eqs (44) and (45), we can solve G ++k and G −+k . By the same token, we can also solve the other two Green functions. The results are Substituting eqs (46) and (47) into eqs (11) and (12) and transforming the summation into integration, we obtain the self-consistent equations for the self-energies, i.e. eqs (13) and (14).
In the weak scattering limit, we first substitute Σ E = 0 and Σ O = − n i λ into the expression of I 0 but keep the imaginary part of Σ E in the denominator Making a variable substitution t = (γk) 2 and t c = (γk c ) 2 , and noticing Γ E is a positive infinitesimal quantity, we have The first term is a principle-value integral, and whether the second term vanishes depends on the sign of ε λ − n i 2 2 2 . Substituting eq. (49) into the self-consistent equations of the self-energies, we obtain the self-energies to the weak scattering limit The DOS in the weak scattering limit ρ (0) in the main text is obtained directly from Γ E . Notice that the four quantities above are of different orders of n i and λ, λ ε . Thus, we take the limits Σ E /ε → 0, Δ O → − n i λ and Γ O → 0 when evaluating the diagonal conductivity to the weak scattering limit, i.e. eq. (30).
Within t-matrix approximation, the self-energy can be expressed by a geometric series which can be represented by the Feynman diagrams shown in Fig. 1(c). A general n-order term is Taking advantage of eqs (46) and (47), the self-consistent equations of the self-energies within t-matrix approximation can be obtained In order to study the effect of a usual scalar disorder on the topologically nontrivial gap, we add the following scalar potential to the Hamiltonian where {μ j } is a set of random numbers to simulate the scalar disorder. We assume 〈 μ〉 = 0. In this case, we need to evaluate the Feynman diagrams shown in Fig. 1(d). Because of the zero value of the average impurity strength, the second term in Fig. 1(d) is zero. And the fourth term has a similar form as the third term. Therefore, the self-consistent equations of the self-energies within SCBA in this case are π λ µ ε ε Σ = + −Σ n n I 1 2 ( ) ( ) ( ) When k c is large, the integral above converges to a constant and hence φ(ε, ε) = − L 2 /(πγ 2 ). Thus, For the special case when the two energy arguments take the values ε + i0 and ε − i0, the correlation functions can be expressed by real functions where B = γ 2 /n i λ 2 . Next, we derive the formula for the derivative of the correlation function. Taking partial derivative with respect to ε′ of eq. (27), we get