This is the accepted version of the article: Dinh Van Tuan, Shaffique Adam, Stephan Roche. Spin dynamics in bilayer graphene: Role of electron-hole puddles and Dyakonov-Perel mechanism. Physical Review B, 94(4):2016, article 041405(R) Available at: https://doi.org/10.1103/PhysRevB.94.041405 All rights reserved Spin Dynamics in Bilayer Graphene : Role of Electron-Hole Puddles and the Dyakonov-Perel Mechanism Dinh Van Tuan,1 Shaffique Adam,2,3 and Stephan Roche1,4 1 Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC and The Barcelona Institute of Science and Technology, Campus UAB, Bellaterra, 08193 Barcelona, Spain 2 Centre for Advanced 2D Materials and Physics Department, National University of Singapore 3 Yale-NUS College, 138614, Singapore 4 ICREA, Institució Catalana de Recerca i Estudis Avançats, 08070 Barcelona, Spain (Dated: April 29, 2016) Abstract We report on spin transport features which are unique to high quality bilayer graphene, in absence of magnetic contaminants and strong intervalley mixing. The time-dependent spin polarization of propagating wavepacket is computed using an efficient quantum transport method. In the limit of vanishing effects of substrate and disorder, the energy-dependence of spin lifetime is similar to monolayer graphene with a M-shape profile and minimum value at the charge neutrality point, but with an electron-hole asymmetry fingerprint. In sharp contrast, the in- corporation of substrate-induced electron-hole puddles (characteristics of supported graphene either on SiO2 or hBN) surprisingly results in a large enhancement of the low-energy spin lifetime and a lowering of its high-energy values. Such feature, unique to bilayer, is explained in terms of a reinforced Dyakonov-Perel mechanism at the Dirac point, whereas spin relaxation at higher energies is driven by pure dephasing effects. This suggests further electrostatic control of the spin transport length scales in graphene devices. PACS numbers: 72.80.Vp, 73.63.-b, 73.22.Pr, 72.15.Lh, 61.48.Gh 1 Introduction.- Owing to its long spin diffusion length at room-temperature (exceeding several tens of micrometers), single layer graphene (SLG) stands as an unquestionable candidate for the realization of practical devices harvesting the spin degree of freedom, and for more innovation in spintronic appli- cations [1–5]. Such ability to propagate spins over very long distances is due to an intrinsically small spin-orbit coupling (SOC) and hyperfine interaction [6]. However, it has been recently demonstrated that the unavoidable coupling between spin and pseudospin degrees of freedom [7], and resulting interwined quantum dynamics, produce a minimum value of spin lifetime (τs ) at the Dirac point, followed by an enhancement of τs with energy [8, 9], as commonly observed experimentally whatever the substrate and material quality [10, 11]. Bilayer graphene (BLG) differs from SLG by a parabolic band dispersion, however preserving the chiral nature of low-energy electronic excitations. Besides, in contrast to SLG, an electronic bandgap can be induced and tuned in BLG under external electric fields [12]. Transport measurements show critical differences between SLG and BLG [13–15], which are attributed to varying bandstructure and electron- hole puddles characteristics (such as screening strength [16]). Finally, spin lifetimes are generally found to differ substantially from the SLG case [17]. These differences include variations in the absolute values as well as an opposite scaling of τs versus charge density, and a dominating Dyakonov-Perel mechanism in BLG [17]. Kochan and coworkers [18] have recently proposed a phenomenological interpretation of the energy profile of spin lifetime driven by local magnetic moments and thermal or disorder averaging, without capturing any difference in the nature of the spin relaxation mechanism between SLG and BLG. In this Letter, thanks to a fully quantum treatment of spin dynamics in real space, we report on spin transport features which are unique to high quality BLG. The critical role of the substrate to reproduce the typical energy-dependent profile of spin lifetimes (observed experimentally) is revealed, with a large enhancement of τs near the charge neutrality point driven by electron-hole puddles and local electric field effects. By comparison with the unsupported (pristine) graphene case, our findings identify that a reinforcement of the Dyakonov-Perel mechanism occurs at low energy, whereas higher energy spin lifetimes are dictated by quantum dephasing effects. Our study points towards the uniqueness of the BLG bandstructure in presence of spin-orbit interaction and non-uniform energy-dependent spin precession frequency to capture spin transport fingerprints, and suggest the possibility to electrostatically monitor the spin diffusion length scales. BLG Hamiltonian in presence of spin-orbit interaction and electron-hole puddles.- BLG can be con- 2 sidered as two coupled SLGs with the top layer shifted a carbon bond from the bottom layer (Fig.1(a)). Consequently, BLG consists of four carbon atoms in its unit cell, two carbons A1 , B1 from the unit cell of the bottom SLG and A2 , B2 from the top layer where B2 places on the top of A1 , namely dimer sites and B1 , A2 are called non-dimer sites. In the tight-binding model the full Hamiltonian for BLG reads : H = HSLG T + HSLG B + H0Inter + HSOC Inter (1) where the first and second terms (intralayer parts HSLG l ) are the Hamiltonians for each single layer involving the SOCs effect (intrinsic λI and Rashba types λlR ), the different potential energies ∆ of the top (l = 1) and bottom layer (l = 2) as well as the long range potential simulating the electron-hole puddles V (r) [16]. ∑ 2i ∑ l + HSLG l = − γ0 a+ l,i bl,j + λR al,i⃗z · (⃗s × d⃗ij )bl,j 3 ⟨ij⟩ ⟨ij⟩ 2i ∑ + λI a+ s · (d⃗kj × d⃗ik )bl,j + h.c. l,i⃗ 9 ⟨⟨ij⟩⟩ ∆∑ { } − (−1)l a+ + l,i al,i + bl,i bl,i + V (r) (2) 2 i Where al,i (bl,i ) is the annihilation operators acting on Ai (Bi ) in layer l. d⃗ij is the unit vector pointing from j to i, k is the common nearest neighbor of i and j. In this paper we use intralayer-intrinsic SOC λI = 12µeV [19] (Fig.1(a)), and intralayer-Rashba SOC λlR = λ0 − (−1)l 2λBR (Fig. 1(c)) which involves two contributions, one from the bulk-inversion-asymmetry induced by the adjacent layer with λ0 = 5µeV [19] and the other 2λBR = 10 × E[V/nm]µeV which is field dependent. In this calculation we choose 2λBR = 2.5µeV corresponding to an electric field E = 0.25[V/nm] independent of the charge density [19], which is a reasonable approximation from experimental considerations [20]. The third term in Eq.(1) is the non-spin-orbit coupling part of the interlayer Hamiltonian ∑ ∑ H0Inter = γ1 a+ 1,i b2,i + γ3 b+ 1,i a2,j i ⟨ij⟩ ∑{ } − γ4 a+ + 1,i a2,j + b1,i b2,j + h.c. ⟨ij⟩ where the first term in above Hamiltonian describes the interlayer hopping (γ1 = 340 meV) between dimer sites {A1 , B2 } [19]. The second term denotes the interlayer coupling between B1 and its adjacent A2 with γ3 = 280 meV. The third term corresponds to hopping integral from A1 to its adjacent A2 , and from B1 to its adjacent B2 , with γ4 = 145 meV (Fig. 1(b)). All these parameters have been derived from 3 the ab-initio calculations [19, 21]. Finally, the SOC part of the interlayer interaction is described by the final term in Eq.(1), HSOC Inter which reads : 2iλ4 ∑ { + ⃗∥ ⃗∥ } a1,i⃗z · (⃗s × dij )a2,j + b1,i⃗z · (⃗s × dij )b2,j + h.c. + 3 i,j ∥ where d⃗ij is the unit vector of the projection of vector d⃗ij on the horizontal plane. This part is the Rashba-type spin-orbit interaction (λ4 = −12µeV [19]) between sites with interlayer hopping term described by γ4 (Fig. 1(c)). Additional much smaller terms as described in Ref. [19] are here neglected. ∑ The electron-hole puddles in BLG are simulated by a long-range potential V (r) = N j=1 ϵj exp[−(r− Rj )2 /(2ξ 2 )], with ξ = 3.5 nm, a value extracted from self-consistent calculations [16]. The magnitude of the onsite potential ϵj is randomly chosen within [−W, W ] with W = 35 meV for SiO2 substrate and W = 11 meV for hBN substrate [16]. The impurity concentration is ni = 1012 cm−2 (0.04%) for a SiO2 substrate and ni = 1011 cm−2 (0.004%) for a hBN substrate. A2 γ0 B2 Dimer Non-dimer b) a) λI γ3 γ4 γ4 A2 B2 γ1 A1 B1 γ0 λI Top layer A1 B1 A2 λ +2λ B 0 BR 2 c) Bottom layer λ4 λ4 λ0-2λBR A1 B1 FIG. 1: (color online): (a) Sketch of BLG composed of a top (in red) and a bottom (in black) layers. The intralayer intrinsic SOC is shown in blue. (b) Schema of non-spin-orbit interaction between one type of carbon atom to the nearest carbons of the other kind. (c) The same in (b) but for Rashba-type SOC. Spin dynamics methodology.- The spin dynamics of electron in BLG is investigated using the time- dependent evolution of the spin polarization Pz (E, t) of propagating wavepackets [8], which is computed through ⟨Ψ(t)|sz δ(E − H) + δ(E − H)sz |Ψ(t)⟩ Pz (E, t) = (3) 2⟨Ψ(t)|δ(E − H)|Ψ(t)⟩ where sz is the z component of the Pauli matrices and δ(E − H) is the spectral measure operator. The evolution of the wavepackets |Ψ(t)⟩ is obtained by solving the time-dependent Schrödinger equation 4 [22], starting from a wavepacket |Ψ(t = 0)⟩ in an out-of-plane (z direction) polarization. An energy broadening parameter η = 13.5 meV is introduced for expanding δ(E − H) through a continued fraction expansion of the Green’s function [22]. This method has been previously used to investigate spin relax- ation in gold-decorated graphene [8], hydrogenated graphene [23], or SOC coupled graphene under the effect of electron-hole puddles [9]. 1400 1 1200 0.5 Pz (t) 1000 0 E=0meV -0.5 E=-100meV Time (ps) 800 E=+100meV 0 100 200 300 400 500 600 t (ps) 400 200 TΩ τs 0 -300 -200 -100 0 100 200 300 E (meV) FIG. 2: (color online) Main frame: Spin lifetime τs (red line) and spin precession time TΩ (black line) in spin- orbit coupled BLG, in absence of electron-hole puddles. In this limit, τs has the same characteristic M-shape as for monolayer graphene, but with some electron-hole asymmetry due to skew interlayer hopping. Inset: spin polarization Pz (t) at some energies (solid lines) together with fits to the function Pz (t) = cos(2πt/TΩ )e−t/τs (dashed lines) from which TΩ and τs are extracted (main frame). Spin dynamics and dephasing in the ultraclean limit- We first consider the situation of pristine BLG in absence of microscopic disorder (V (r) = 0; ∆ = 0), and in which only the uniform SOC and the small energy broadening dictate the spin lifetime characteristics. Fig. 2 (inset) shows Pz (E, t) (solid lines), which exhibits an oscillatory pattern typical for spin precession together with an exponential decay which dictates the loss of spin information. A significant electron-hole asymmetry is observed for two chosen energies E = ±100 meV, which correspond to a charge density of ±5 × 1012 cm−2 . A faster oscillation of the spin signal is seen at E = 100 meV (green line) in contrast to the slower oscillation observed at E = −100 meV (red line). The spin relaxation at the charge neutrality point is even faster and likely driven by the same interwoven dynamics of spin and pseudospin degrees of feeedom, as unveiled for SLG [8]. 5 From the fits of the numerical data using Pz (t) = cos(2πt/TΩ ) exp(−t/τs ) (dashed lines), the spin precession time TΩ and τs are extracted (see Fig. 2 (inset), dashed lines). The electron-hole asymmetry is observed in both TΩ and τs with larger values for the hole side compared to the electron side. This phenomenon has been noticed by Diez and Burkard [24], who proved that the contribution of skew interlayer hopping term γ3 leads to the reduction (increase) of spin splitting energy ∆E in the hole (electron) side (note that TΩ ∼ 1/∆E). The largest variation of TΩ (E) occurs in the vicinity of the charge neutrality point, and this non-uniformity results in stronger dephasing effects and shortest spin lifetime, as discussed for the SLG case [8, 9]. On the other hand, the obtained value for τs varies from 100 ps to 1.4 ns, which is the typical range of experimental data [17]. Finally, τs (E) for BLG shares another similar additional feature with SLG, that is a downturn at higher energies, which has been related to the contribution of trigonal warping [9]. 1400 1 E=0meV 1200 E=-100meV 0.5 Pz (t) E=+100meV 1000 0 800 0 200 400 600 800 1000 τs (ps) t (ps) 600 τs τs: pristine BLG 400 200 0 -300 -200 -100 0 100 200 300 E (meV) FIG. 3: (color online) Main frame: spin lifetime τs (solid line) for graphene on SiO2 substrate in the comparison with τs of pristine spin-orbit coupled BLG (dashed line). Inset: spin polarization Pz (t) at some energies (solid lines) and their fits to the function Pz (t) = cos(2πt/TΩ )e−t/τs (dashed lines) To study the substrate effect, we introduce electron-hole puddles through the long range potential V (r), which well reproduce the measured charge density fluctuations for graphene either supported on SiO2 or hBN [16]. Fig. 3 (inset) shows the evolution of spin polarization Pz (t) (solid lines) for graphene on SiO2 substrate and the corresponding fits (dashed lines) from which τs is extracted (main frame, green solid line). Remarkably, the electron-hole puddles associated to SiO2 substrate provokes an inversion in 6 the energy-dependent profile of τs , with a peak at CNP, as reported in experiments in the low temperature regime [17]. The absolute values, energy dependence as well as the electron-hole asymmetry of the extracted τs provide a consistent support to the analysis of state-of-the-art experimental data without the need to introduce any magnetism in the problem [17]. It is worth mentioning that the obtained value for τs (green solid line) close to CNP is larger than in the case of pristine spin-orbit coupled BLG with no puddles (red dashed line). Such enhancement of τs close to CNP, and driven by electron- hole puddles, can be rationalized as a reinforcement of the Dyakonov-Perel mechanism [24]. Indeed, the presence of electron-hole puddles generates elastic scattering which act on the spin precession and produce a motional narrowing phenomenon. This is particularly strong at the charge neutrality owing to the specific bandstructure of BLG where parabolic bands and higher density of states are obtained and favour an enhancement of the scattering probability, when compared to SLG [25]. This enhanced scattering more efficiently impedes spin precession, which in absence of disorder, is the mechanism driving to relaxation. Such effect is inactive in SLG [9], bringing an essential difference between both types of structures. 600 1 E=0meV E=-100meV 500 0.5 E=+100meV Pz(t) 0 400 Times (ps) -0.5 0 100 200 300 400 500 300 t (ps) 200 τs: hBN, ∆=2meV 100 TΩ: pristine BLG 0 -300 -200 -100 0 100 200 300 E (meV) FIG. 4: (color online): Main frame: Spin relaxation time τs (blue lines) for graphene on hBN substrate. TΩ for the pristine spin-orbit coupled BLG is shown in black dashed line. Inset: spin polarization Pz (t) at some energies (solid lines) together with fits to Pz (t) = cos(2πt/TΩ )e−t/τs (dashed lines). For states close to CNP, the fit is made with Pz (t) = Pz (t0 )e−(t−t0 )/τs . For the case of hBN substrate, the situation becomes more complicated. Scattering due to electron- hole puddle alone (W = 11 meV and ni = 1012 cm−2 ) is too weak to modify the results when compared 7 to the unsupported case (Fig.1). However, the scattering strength close to CNP due to electron-hole puddles is likely enhanced by the formation of a pseudogap induced by the electric field, so far neglected. Indeed in presence of an external electric field, the weak interaction between hBN and graphene layers breaks the symmetry between top and the bottom layers, an effect which can be modelled by adding an energy difference between layers [12, 26]. To account for it, we thus introduce a small energy difference between the top and the bottom layers, which would open a real gap of ∆ = 2 meV for the pristine BLG. Here however, the disorder potential stemming from the electron-hole puddles is strong enough to wipe out the gap (W ≫ ∆), while maintaining spatially uncorrelated local energy fluctuations between layers. Fig. 4 shows the resulting spin lifetime for BLG on hBN. The presence of this extra energy asymmetry between top and bottom layers turns out to play an essential role on scattering effects close to the CNP. Similarly to the case of SiO2 substrate, the Dyakonov-Perel mechanism is reinforced close to the CNP and spin relaxation is reduced. The fits of the spin polarization (solid lines in the inset) with the function Pz (t) = cos(2πt/TΩ )e−t/τs (dashed lines) allow the extraction of the spin lifetime τs (blue line in the main frame) exhibiting a sharp peak close to CNP. Since the spin polarization shows an initial transient fast decay [23], the long time spin relax- ation physics is better captured by fitting the data close to CNP using from the function Pz (t) = Pz (t0 )e−(t−t0 )/τs (fit starts from t0 = 60ps, see blue dashed line in inset). The extracted τs is seen to be very sharp close to CNP (as observed in a recent experiment [20]). The enhancement of τs close to CNP originates from the more efficient motional narrowing driven by the Dyakonov-Perel mechanism. Indeed, the strong variation of spin precession time TΩ (black dashed line) close to CNP gives rise to spin dephas- ing and short spin lifetime. Here, for Gaussian correlated disorder with the chosen parameters, a semi- ( 2 ) ~ ~ classical transport calculation gives a minimum momentum scattering time τp = 4π 1 mξ 2 W W ≈ 100 fs for a hBN substrate, and about 10 fs for bilayer graphene on a SiO2 substrate. We note that experimen- tally τp for BLG on SiO2 [20, 27] and on hBN [26] are consistent with these estimates, and in all cases τp ≪ TΩ which satisfies the criterion to enter into the DP regime and supports our interpretation of an enhanced DP mechanism at the Dirac point. Overall, scattering events induced by electron-hole puddle together with the pseudogap act against the spin dephasing of cleaner samples and consequenly leads to the enhancement of spin relaxation time close to CNP, indicating that the Dyakonov-Perel mechanism govers the low-energy spin lifetime in BLG. Finally, It is worth mentioning that the pseudogap of graphene on hBN is not only induced by electric 8 field but also by the staggered potential which captures the interaction between the graphene lattice and hBN [28]. This can explain the fact that we observe a sharp peak for spin relaxation time on hBN whereas a broaden peak is experimentally observed for graphene on SiO2 . Discussion and conclusion.- Recently, D. Kochan and coworkers [18] suggested that the origin for spin relaxation in BLG is the same with SLG, namely resonant scattering by magnetic impurities. By varying the strength of electron-hole puddles and broadening factors, the experimental data could roughly reproduce τs with the upturn at CNP in BLG or the downturn in SLG. However, this scenario predicts an Elliot-Yafet type of relaxation [23], whereas experiments on BLG clearly evidence a scaling behavior as τs ∼ 1/τp , indicating the predominance of the Dyakonov-Perel mechanism [17]. Here we have found that even for long mean free paths, the impact of electron-hole puddles due to silicon oxide or hBN substrates make the Dyakonov-Perel mechanism predominating over spin dephasing, a fact unique to BLG. A characteristic peak at the CNP will be seen for the SiO2 or will be very sharp around the CNP for the hBN substrate. Those findings provide a consistent interpretation of all reported experiments on BLG, without the need of introducing additional relaxation mechanism driven by magnetic impurities [17, 20]. This work has received funding from the European Union Seventh Framework Programme under grant agreement 604391 Graphene Flagship. S.R. acknowledges the Spanish Ministry of Economy and Competitiveness for funding (MAT2012-33911), the Secretaria de Universidades e Investigacion del De- partamento de Economia y Conocimiento de la Generalidad de Catalua and the Severo Ochoa Program (MINECO SEV-2013-0295). S.A. acknowledges support from the National Research Foundation of Singapore under its Fellowship programme (NRF-NRFF2012-01) and computational resources from the Centre for Advanced 2D Materials (an NRF mid-sized centre). [1] N. Tombros, C. Jozsa, M. Popinciuc, H.T. Jonkman, B.J. Van Wees, Nature (London), 448 (2007), 571 [2] P. Seneor, B. Dlubak B, M.B. Martin, A. Anane, H. Jaffres, and A. Fert, MRS Bulletin 37, 1245 (2012). B. Dlubak, M.-B. Martin, C. Deranlot, B. Servet, S. Xavier, R. Mattana, M. Sprinkle, C. Berger, W.A. de Heer, F. Petroff, A. Anane, P. Seneor and A. Fert, Nature Phys. 8, 557 (2012). [3] H. Dery, H. Wu, B. Ciftcioglu, M. Huang, Y. Song, R. Kawakami, J. Shi, I. Krivorotov, I. Zutic, L.J. Sham, IEEE Trans. Electron Devices 59, 259262 (2012). [4] W. Han, R.K. Kawakami, M. Gmitra and J. Fabian, Nature Nanotech. 9, 794 (2014) . 9 [5] S. Roche, J. Åkerman, B. Beschoten, J.-Ch. Charlier, M. Chshiev, S. Prasad Dash, B. Dlubak, J. Fabian, A. Fert, M. Guimares, F. Guinea, I. Grigorieva, C. Schönenberger, P. Seneor, C. Stampfer, S. O. Valenzuela, X. Waintal, and B. van Wees, 2D Materials 2, 030202 (2015). [6] D. Huertas-Hernando, F. Guinea and A. Brataas. Phys. Rev. B 74, 155426 (2006). M. Gmitra, S. Konschuh, C. Ertler, C. Ambrosch-Draxl and J. Fabian, Phys. Rev. B 80, 235431 (2009). A. H. Castro Neto and F. Guinea, Phys. Rev. Lett. 103, 026804 (2009). V. K. Dugaev, E. Ya. Sherman, and J. Barnas, Phys. Rev. B 83, 085306 (2011). [7] E. I. Rashba, Phys. Rev. B 79, 161409 (2009). [8] D. Van Tuan, F. Ortmann, D. Soriano, S. O. Valenzuela and S. Roche, Nature Phys. 10, 857 (2014). [9] D. Van Tuan, F. Ortmann, A.W. Cummings, D. Soriano, and S. Roche, Scientific Reports 6, 21046 (2016). A. W. Cummings and S. Roche, Phys. Rev. Lett. 116, 086602 (2016). [10] K. Pi, W. Han, K.M. McCreary, A.G. Swartz, Y. Li, and R.K. Kawakami, Phys. Rev. Lett. 104, 187201 (2010). S. Jo, D.-K. Ki, D. Jeong, H.-J. Lee, and S. Kettemann, Phys. Rev. B 84, 075453 (2011). P.J. Zomer, M. H. D. Guimaraes, N. Tombros, and B.J. van Wees, Phys. Rev. B 86, 161416(R) (2012). T. Maassen, J.J. van den Berg, N. IJbema, F. Fromm, Th. Seyller, R. Yakimova, and B.J. van Wees, Nano Lett. 12, 1498 (2012). W. Han, J.-R. Chen, D. Wang, K.M. McCreary, H. Wen, A. G. Swartz, J. Shi, and R.K. Kawakami, Nano Lett. 12 3443 (2012). M. H. D. Guimarães, A. Veligura, P.J. Zomer, T. Maassen, I.J. Vera-Marun, N. Tombros, and B.J. van Wees, Nano Lett. 12, 3512-3517 (2012). [11] M.H.D. Guimarães, P.J. Zomer, J. Ingla-Aynes, J.C. Brant, N. Tombros and B.J. van Wees, Phys. Rev. Lett. 113, 086602 (2014). M. Drögeler, F. Volmer, M. Wolter, B. Terres, K. Watanabe, T. Taniguchi, G. Güntherodt, C. Stampfer and B. Beschoten, Nano Lett. 14, 6050 (2014). M. Venkata Kamalakar, C. Groen- veld, A. Dankert, S.P. Dash, Nature Comm. 6, 6766 (2015). [12] E. McCann, K. Kechedzhi, V. I. Falko, H. Suzuura, T. Ando, and B. L. Altshuler, Phys. Rev. Lett. 97, 146805 (2006). E. V. Castro, K.S. Novoselov, S.V. Morozov, N. M. R. Peres, J.M.B. Lopes Dos Santos, J. Nilsson, F. Guinea, A.K. Geim and A.H. Castro Neto, Phys. Rev. Lett. 99, 216802 (2007). [13] C. Ojeda-Aristizabal, M. Monteverde, R. Weil, M. Ferrier, S. Gueron, and H. Bouchiat, Phys. Rev. Lett. 104, 186802 (2010). J. B. Oostinga, H. B. Heersche, X. Liu, A. F. Morpurgo and L. M. K. Vandersypen, Nature Materials 7, 151 (2007). [14] K. Kechedzhi, V.I. Falko, E. McCann, B.L. Altshuler, Phys. Rev. Lett. 98 176806 (2007). [15] M. I. Katsnelson, Phys. Rev. B 76, 073411 (2007). 10 [16] S. Adam, S. Jung, N. N. Klimov, N. B. Zhitenev, J. A. Stroscio and M. D. Stiles, Phys. Rev. B 84, 235421 (2011). [17] T.- Y. Yang, J. Balakrishnan, F. Volmer, A. Avsar, M. Jaiswal, J. Samm, S. R. Ali, A. Pachoud, M. Zeng, M. Popinciuc, G. Güntherodt, B. Beschoten, and B. Ozyilmaz, Phys. Rev. Lett. 107, 047206 (2011). W. Han and R. K. Kawakami, Phys. Rev. Lett. 107, 047207 (2011). A. Avsar, T.-Y. Yang, S. Bae, J. Balakrishnan, F. Volmer, M. Jaiswal, Z. Yi, S.R. Ali, G. Güntherodt, B.H. Hong, B. Beschoten, B. Ozyilmaz, Nano Lett. 11 2363 (2011). J. Ingla-Aynes, M. H. D. Guimarães, R. J. Meijerink, P. J. Zomer, and B. J. van Wees, Phys. Rev. B 92, 201410 (R) (2015). [18] D. Kochan, M. Gmitra, and J. Fabian, Phys. Rev. Lett. 112, 116602, (2014). D. Kochan, S. Irmer, M. Gmitra, and J. Fabian, Phys. Rev. Lett. 115, 196601 (2015). [19] S. Konschuh, M. Gmitra, D. Kochan, and J. Fabian, Phys. Rev. B 85, 115423 (2012). [20] A. Avsar, I.J. Vera-Marun, J.Y. Tan, K. Watanabe, T. Taniguchi, and B. Ozyilmaz, arXiv:1602.07817. An estimation of the experimental electric field needs to electrostatically dope the BLG to 8 × 1012 cm−2 is ∼ 0.7 V/nm. [21] J. Jung and A. H. MacDonald, Phys. Rev. B 89, 035405 (2014). [22] L. E. F. Foa Torres, S. Roche, and J. C. Charlier in Introduction to Graphene-Based Nanomaterials: From Electronic Structure to Quantum Transport (Cambridge University Press, Cambridge, 2014). [23] D. Soriano et al., 2D Materials, 2, 022002 (2015). [24] M. Diez and G. Burkard, Phys. Rev. B 85, 195412 (2012). [25] J. Nilsson, A. H. Castro Neto, F. Guinea, and N. M. R. Peres, Phys. Rev. B 78, 045405 (2008). [26] A. Ramasubramaniam, D. Naveh and E. Towe, Nano Lett., 11 (3), 1070-1075 ( 2011); P.J. Zomer, M.H. D. Guimaraes, J.C. Brant, N. Tombros and B.J. van Wees, Appl. Phys. Lett. 105, 013101 (2014). [27] M. Monteverde, C. Ojeda-Aristizabal, R. Weil, K. Bennaceur, M. Ferrier, S. Gueron, C. Glattli, H. Bouchiat, J.N. Fuchs, and D.L. Maslov, Phys. Rev. Lett. 104, 126801 (2010). [28] J. Jung, A.M. DaSilva, A. H. MacDonald, and S. Adam, Nature Comm. 6, 6308 (2015). 11