Quantum Hall Effect in Graphene with Interface-Induced Spin-Orbit Coupling Tarik P. Cysne,1 Jose H. Garcia,2 Alexandre R. Rocha,3 and Tatiana G. Rappoport1 1Instituto de Fı́sica, Universidade Federal do Rio de Janeiro, Caixa Postal 68528, Rio de Janeiro 21941-972, RJ, Brazil 2Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC and The Barcelona Institute of Science and Technology, Campus UAB, 08193 Barcelona, Spain. 3Instituto de Fı́sica Teórica, Universidade Estadual Paulista (UNESP), Rua Dr. Bento T. Ferraz, 271, So Paulo, SP 01140-070, Brazil. (Dated: November 22, 2017) We consider an effective model for graphene with interface-induced spin-orbit coupling and calculate the quantum Hall effect in the low-energy limit. We perform a systematic analysis of the contribution of the different terms of the effective Hamiltonian to the quantum Hall effect (QHE). By analysing the spin-splitting of the quantum Hall states as a function of magnetic field and gate-voltage, we obtain different scaling laws that can be used to characterise the spin-orbit coupling in experiments. Furthermore, we employ a real-space quantum transport approach to calculate the quantum Hall conductivity and investigate the robustness of the QHE to disorder introduced by hydrogen impurities. For that purpose, we combine first-principles calculations and a genetic algorithm strategy to obtain a graphene-only Hamiltonian that models the impurity. I. INTRODUCTION The interaction between a graphene sheet and different sub- strates has attracted great attention in recent years due to the appearance of interesting effects in the graphene layer2–6. This is particularly relevant in context of van der Walls heterostructures7, a vertical stack of bidimensional materials that promise to lead to novel electronics devices and applica- tions. One possibility in this direction is the use of graphene- based van der Walls heterostructures for spintronics. Graphene is a non-magnetic material with very weak spin- orbit coupling (SOC)8–10, due to lightness of carbon atoms. Since its discovery, there were several proposals to introduce and control spin-dependent properties in graphene by induc- ing spin-orbit coupling11–16. Recently, there has been a signif- icant progress in engineering those properties by proximity ef- fect, which allows introducing spin-dependent features while preserving graphene’s high electronic mobility. Graphene has been proximity-coupled with a magnetic thin layer of YIG for transport17 and spin-pumping18 measurements. SOC was also induced by proximity effect in graphene on top of different transition metal Dichalcogenites (TMDC)19–22 and graphene decorated with Gold23. More recently, spins were optically injected in graphene/TMDC systems24,25, also indicating the presence of SOC in the graphene layer. However none of these measurements give clear indications of the type of underlying spin-orbit coupling mechanism that generates the observed phenomena. In this article, we exploit the possibility of using quantum Hall measurements to extract the characteristics of the spin- orbit coupling in graphene. For that purpose, we employ the effective model for graphene with interface-induced spin-orbit coupling from Ref. 1. We then use Landau operators to cal- culate the quantum Hall effect in the low-energy limit. By performing a systematic analysis of the spin-splitting of the Landau levels (LL) in the quantum Hall regime as a function of magnetic field and gate-voltage, we obtain characteristic scaling laws for the splittings produced each type of spin-orbit coupling. The same type of analysis can be used to charac- terise and estimate the SOC in transport experiments. Fur- thermore, we employ a real-space quantum transport approach based on a Chebyshev polynomial expansion of disordered Green functions to calculate the quantum Hall conductivity40. We use this numerical approach to investigate the robustness of the QHE to disorder introduced by hydrogen impurities. The impurities are modelled by an ab initio-derived tight- binding model. For the extraction of the tight-binding param- eters, it is necessary to perform a multiparametric fit. Deter- ministic approaches for the fit can be quite difficult, because of the occurrence of large number of extrema. Therefore, here we proposed the use of an heuristic algorithm to perform this task efficiently. The article is organised as follows: in section II, we intro- duce the tight-binding and low-energy models for graphene with interface-induced SOC. We also present the Hamilto- nian of the system under an external magnetic field, written in terms of Landau operators. In section III, we discuss our analytical results for the energy spectra and Hall conductivity and introduce the scaling laws that can be used to discrimi- nate the different SOC. In section IV we present the heuristic algorithm that was used to extract the tight-binding parame- ters from density functional theory spectra and our numerical approach for the conductivity calculations. We then discuss the results for the effect of hydrogenation on the QHE for graphene with interface-induced SOC. Finally, we present our conclusions in section V. II. THEORY With a combination of density functional theory (DFT) and group theory analysis, Kochan et al. proposed a tight-binding Hamiltonian and its corresponding low-energy approxima- tion for graphene stacked on transition metal dichalcogenides (TMDC) 1. Here, we consider their Hamiltonian with the two most relevant spin-orbit terms for the electronic properties of graphene in the vicinity of the Dirac point: ar X iv :1 71 1. 04 81 1v 2 [ co nd -m at .m es -h al l] 2 1 N ov 2 01 7 2 H = ∑ 〈i,j〉 t c†iscjs + ∑ i ∆ ηci c † iscis + 2i 3 λR ∑ 〈i,j〉 c†iscjs′ [ (̂s× dij)z ] ss′ + i 3 ∑ 〈〈i,j〉〉 λci√ 3 c†iscjs′νij ŝz (1) where c†is = ( a†is, b † is ) and cis = (ais, bis) are the creation and annihilation operators for an electron on a lattice site i and spin s belonging to sublattice A or B, respectively. The first term represents the hopping with amplitude t between π or- bitals in a honeycomb lattice and the second term is an energy offset ∆ between sub-lattices A (ηai = 1) and B (ηbi = −1) due to the super-lattice effect originated by the incommensu- rability of the two lattices. The third contribution is the typical Rashba SOC with strength λR26,27. It arises because the inver- sion symmetry is broken when graphene is placed on top of a TMDC.The fourth contribution is a sublattice-dependent in- trisic SOC with coupling intensities λa and λb. Hξ, the valley- dependent low-energy limit of the Hamiltonian of Eq. 1, has four terms1, Hξ = ~vf (ξσxkx + σyky) + ∆σz + λR(ξσxsy − σysx) + 1 2 ξ(λa(σz + σ0) + λb(σz − σ0))sz. (2) The first two terms are the spinless contributions where vF is Fermi velocity given by vF = 3ta 2~ with lattice constant a, σ is a pseudospin Pauli matrix related to sublattices A and B and kx and ky are the components of the electronic moment relative to the Dirac points. ξ is related with the valley degree of freedom, ξ = + for valleyK and ξ = − for valleyK ′. The third and fourth terms are the Rashba and sublattice resolved intrinsic spin-orbit couplings respectively. The fourth term contains the well known Kane-mele term28,29 with strength λI = (λa + λb)/2, and a valley-Zeeman SOC with strength λVZ = (λa − λb)/2 that couples spin and valley degrees of freedom. In Figure 1 we show the energy spectrum of the Hamilto- nian of Eq. 2 for different combinations of λa, λb and λR to better understand the characteristics of the novel intrin- sic spin-orbit coupling. Differently from the usual case of λa = λb, if only the valley-Zeeman contribution λa = −λb (i. e. λVZ = λa) is present, the spectrum is gapless and the spin degeneracy is broken (see Figure 1(a)). If only Rashba SOC is present, the spectrum is also gapless27. However, any combination λR 6= 0 and λVZ 6= 0 opens a gap in the energy spectrum as shown in Figure 1, panels b and c. For compari- son, Figure 1 (d) presents the energy spectrum for λa = λb (i. e. λI = λa) and λR 6= 0. It is important to note that valley- Zeeman combined with Rashba coupling preserves particle- hole symmetry ( Figures1(b) and 1(c)), while, when Kane- Mele and Rashba are both present, the particle-hole symmetry is broken (Figure1(d)). If a uniform perpendicular magnetic field is applied, ~p → ~π = ~p + e ~A, where ~A is the vector potential in the Landau gauge ~A = B(−y, 0, 0), andB is the intensity of the magnetic field. The low-energy Hamiltonian can be written as -0.1 -0.05 0. 0.05 0.1 E (k ) (a) (b) -0.03 -0.01 0.01 0.03 -0.1 -0.05 0. 0.05 0.1 k (c) -0.03 -0.01 0.01 0.03 k (d) Figure 1. Energy spectrum of the effective Hamiltonian of Eq. 2 for ∆ = 0: (a) λR = 0 and λVZ 6= 0, (b) λVZ � λR, (c) λVZ � λR and (d) λI � λR. Hξ = ~ω(σξaξ + σ−ξa † ξ)− 2iξλR(σ−ξs+ − σξs−), + 1 2 ξ(λa(σz + σ0) + λb(σz − σ0))sz + ∆σz, (3) where σ± = (1/2)(σx ± iσy), s± = (1/2)(sx ± isy), the cyclotron frequency is given by ω = √ 2vF /lB , the magnetic length is lB = √ ~/eB. Landau operators aξ are creation and annihilation operators on valley K and K ′ for ξ = ± : aξ = 1√ 2 ξ ( lB ~ px − 1 lB y + lB i~ py ) ξ , (4) a†ξ = 1√ 2 ξ ( lB ~ px − 1 lB y − lB i~ py ) ξ . (5) The Hamiltonians Hξ=± are block diagonal with each block indexed by an occupation number n. The two lowest blocks are 1 × 1 and 3 × 3 matrices and higher blocks are 4 × 4 matrices for both valleys. The energies Eξn,i and eigenvectors |ψξn,i〉 are indexed by the valley ξ, the occupation number n and i, that labels the eigenvalues and eigenvectors of a given 3 block n. For more details, see Appendix A. The transverse Hall conductivity σxy can be calculated in the framework of the Kubo formula30 σxy = ie3Bv2F 2π ∑ ξ=± ∑ n,n′ ∑ i,i′ ( f(Eξn,i)− f(Eξn′,i′) ) 〈 ψξn,i ∣∣∣σx∣∣∣ψξn′,i′ 〉〈 ψξn′,i′ ∣∣∣σy∣∣∣ψξn,i〉 (Eξn,i − E ξ n′,i′)(E ξ n,i − E ξ n′,i′ + i0+) (6) where f(E) is the Fermi-Dirac distribution function. We will not consider thermal effects here and the Fermi-Dirac distri- bution has a Heaviside function profile, f(E) = Θ(Ef − E). Expressions for the eigenvectors are computed in appendix A. III. ANALYTICAL RESULTS A. Weak spin-orbit coupling Here, we show and discuss our results on the quantum Hall effect in under the effect of Rashba and valley-Zeeman SOCs, that are present in graphene-TMDC heterostructures. Exper- iments and first-principle calculations point out to values of the coupling constants in the range of 0.1 − 10 meV21,22. We express our results in terms of a gate voltage which con- trols the Fermi energy. This voltage is related to the elec- trons’ Fermi momentum in graphene by kf = Ef/~vf =√ απVg 31, where α depends on the substrate. Here, we use α = 7.2 × 1010V −1cm−2, which is an appropriate value for either silicon oxide or TMDC substrates. To be compatible with the analysis of experimental results, instead of showing our energy spectra, we present fan dia- grams of the density of states as a function of gate voltage VG and magnetic field B. We begin by presenting the well known spin-splitting generated by λR. In this case, the n = 0 Lan- dau level is spin-degenerate, while all other levels are spin- split, as illustrated in Fig. 2. It is also clear, from the Landau level splittings of Fig2(a) and the quantum Hall conductivity of Fig.2(b) that, as expected, the splitting increases with the level number n. We proceed to analyse the spin-splitting generated by λVZ: similarly to the previous case, the n = 0 Landau is spin- degenerate while all other levels are spin-split, as illustrated in Fig.3. Again, the Landau level splittings of Fig. 3(a) and the Quantum Hall conductivity of Fig. 3(b) show the increase of the splitting with the level number n and it is considerably larger than the spin-splitting produced by Rashba, even for weak couplings and small n. To compare the spin-splitting produced by Rashba with the one by the valley-Zeeman SOC, we select the difference in energy of the n = 3 levels ∆VG as a function of the spin- orbit strength for various values of the external magnetic field . From Fig.4, we see a very different dependence of ∆VG as a function of λR (panel (a)) and λVZ (panel (b)). ∆VG varies linearly with λVZ and quadratically with λR. For the case of valley-Zeeman SOC, ∆VG is 1-5 V even for very weak cou- plings of the order of λVZ = 1 − 5 meV and can be resolved Figure 2. (a) Landau Fan diagram and (b) Hall conductivity as function of gate voltage for λR = 10 meV, λVZ = 0. Figure 3. (a) Landau Fan diagram and (b) Hall conductivity as function of gate voltage for λR = 0, λVZ = 3 meV. experimentally. The combined effect of Rashba and valley-Zeeman in the weak coupling regime can be seen in Figures 5 and 6. In Fig- 4 Figure 4. (a) Splitting (∆VG) of n=3 Landau level as function of λR with λVZ = 0.4 meV fixed and different magnetic fields. (b) Voltage splitting (∆VG) of n=3 Landau level as function of λVZ with λR = 1 meV fixed, and different magnetic fields. ure 5, we consider the case where λVZ > λR while in Figure 6 we consider the case where λVZ < λR. Because of the differ- ent functional laws presented in Fig.4, if both couplings have similar strengths the splitting is dominated by λVZ. However, as can be seen in these figures, it is difficult to extract informa- tion on the nature and strength of the spin-orbit coupling in a graphene heterostructure from fan diagrams and quantum Hall measurements. The spin-splitting of the Landau Levels and the quantum Hall plateaus for graphene with valley-Zeeman SOC are very similar to the ones for graphene with Rashba SOC. To address this issue, we found two different scaling laws for the spin-split states that can be used to distinguish between systems with valley-Zeeman SOC and ones with Rashba SOC. For the case of valley-Zeeman SOC, the scaled quantity is ∆VG, the voltage difference between spin-split states with same quantum number n. If we consider ∆VG/n as a func- tion of the magnetic field, all curves for different values of n collapse perfectly into a single one, as can be seen in Figure 7(a). If λR 6= 0, the scaling starts to fail whenever the main contribution for the spin-splitting is λR (see Figure 7(b)-(c)), indicating that this scaling law is a characteristic of valley- Zeeman SOC. On the other hand, if we consider λR 6= 0 and λVZ = 0, we need to use a different scaling law to collapse all curves. In this case, we use the quantity 4(V 2 G,s− V 2 G,s′)/n 2 where VG,s and VG,s′ are the voltage of the two spin-split states (repre- sented by s and s′) for a given n. This simple analysis, that can be used in transport measurements, is able to determinate what is the main SOC in graphene and the deviations from Figure 5. (a) Landau Fan diagram and (b) Hall conductivity as func- tion of gate voltage for λR = 1 meV, λVZ = 5 meV. Figure 6. (a) Landau Fan diagram (b) Hall conductivity as function of gate voltage for λR = 6 meV, λVZ = 3 meV. the scaling laws can also estimate the relative contributions of valley-Zeeman and Rashba SOC. B. Strong Spin-Orbit coupling Regime Let us now discuss the case where the spin-orbit cou- pling is of the order of tens of meV, as estimated via weak anti-weak localisation measurements on graphene-TMDC heterostructures21. Rashba combined with valley-Zeeman SOCs produce particle-hole symmetric Landau levels spectra in both strong and weak coupling regimes. This is a conse- 5 Figure 7. 2∆V 2 G/n as a function of the magnetic field B for the Landau levels n = 2, 3, 4, 5 with λVZ = 3 meV (a) λR = 0, (b) λR = 3 meV and (c) λR = 6 meV. Figure 8. 4(V 2 G,s − V 2 G,s′)/n 2 as a function of the magnetic field B for the Landau levels n = 2, 3, 4, 5, 6 with λR = 10 meV and (a) λVZ = 0 , (b) λVZ = 2 meV and (c) λVZ = 4 meV. Figure 9. Landau Fan diagrams for (a) λR = 20 meV, λVZ = 30 meV and (b) λR = 40 meV, λVZ = 20 meV. quence of particle-hole symmetric spectra in absence of mag- netic field (see Figure 1 (b) and (c) ). However, as we can seen in Figure 9, the structure of spin-splittings and plateaus is more complex in the strong coupling regime, due to the presence of several level crossings, that make it difficult to perform quantitative analysis and estimations based solely on spectra and conductivity profiles. On the other hand, Rashba imprints a clear signature in the fan diagram: for strong λR, the fan diagram acquires two extra lateral fans and the separa- tion between the main fan and the satellite ones is proportional to λR. Furthermore, if both λR and λVZ are strong, there is a splitting of the n = 0 level that does not occur in the case of pure Rashba SOC. C. Possible applications to other systems Figure 10. Landau Fan diagram for for (a) λR = 20 meV, λI = 30 meV and (b) λR = 40 meV, λI = 20 meV. Instead of considering the interplay between Rashba and valley-Zeeman SOCs, we can look at the structure of the Lan- dau fan diagram of graphene with interface-induced Rashba and intrinsic SOCs (λR and λI,) as in the case of graphene intercalated with gold. An example of this intercalated struc- ture is a recent experiment on van der Waals heterostructures of graphene-gold-hBN23. We show in figure 7, the spectra for (a) λI > λR and (b) λI < λR. In both cases, the particle-hole symmetry is broken, a characteristic of the interplay between Kane-Mele SOC and Rashba SOC. For λI > λR, it presents a topological gap at V = 0 while the gap is closed if λI < λR. IV. NUMERICAL RESULTS A. Graphene-only Hamiltonian for hydrogen on graphene. Hydrogen adsorption is one of the main sources of contam- ination when manufacturing graphene. Hydrogen hybridizes 6 with pz orbitals in graphene, modifying the local symme- try from sp2 to sp3, and create midgap-states32,33. To anal- yse the effect of hydrogenation in the quantum Hall effect, we employ a real-space quantum transport approach40. For an optimized use of computational resources, that allow us to study large systems, we need to find a single-orbit tight- binding Hamiltonian for hydrogen on graphene that repro- duces a band-structure obtained by density functional theory (DFT). We propose the following graphene-only Hamiltonian for hydrogen adatoms incorporated in a graphene layer Hh = εIc † IcI + tI ∑ 〈I,j〉 c†Icj + h.c+ ∑ 〈〈I,j〉〉 t (2) I c†Icj (7) where c†I and cI are the creation and annihilation operators for an electron on the adsorption site I , and 〈I, j〉 and 〈〈I, j〉〉 represent the sum over nearest and next-nearest neighbours of the adsorption site respectively. tI and t(2)I are parameters that have to be adjusted to properly fit the DFT band-structure. In Fig. 12.a we show the band-structure of a 5×5 graphene supercell with one hydrogen adatom. The first-principles calculations are based on DFT34,35. as implemented in the SIESTA code36, using GGA functional approximation follow- ing the PBE approach37. The pseudopotential were obtained through Troullier-Martins scheme38 and a double-ζ polarised basis set was used to described the electronic orbitals. The self-consistent cycle was performed using 16 × 16 × 1 k- sampling of the Brillouin zone. The structural relaxation was performed using conjugate gradient minimization until the forces were smaller than 0.01 eV/Å. The gap in the band struc- ture artificially originates from the broken sub-lattice symme- try due to the arbitrary choice of an adsorbtion site, which in this case belongs to the sub-latticeA, and disappears when the adatoms are placed randomly, because in average the symme- try will be restored. Multiparametric fits such as the present one are usually dif- ficult to be performed by deterministic approaches due to the occurrence of large number of extremas. Therefore, we pro- pose the use of an heuristic algorithm to efficiently perform this task. A simple algorithm that is inspired in the natural evolutions to find an optimal solution for a given problem is called genetic Algorithm, the scheme is outlined in Fig. 11 and follows the logic 1. Primordial Generation: A population Nt individual with random genome vector xp is generated, within a hypercube of allowed genomes. 2. Breeding: 2Np < Nt pairs of individuals are choicen randomly and is combined to give raise to No < Np ofsprings. 3. Ranking: All individuals parents are offspings are avaliated and ranked through the fittness function. 4. Reinsetion: The least fitned No individuals are termi- nated while keeping the best Nt individuals. 5. Mutation: An stochastic alteration of the i-th genome vector occurs with probability pM . Then the breeding phase occurs again For the primordial generation phase we consider each com- ponent of the genome vector to lie within a range of x ∈ [−x0, x0], with x0=1eV. The selection of pairs is performed by following stochastic universal sampling with Np = Nt/2, and it is combined using an intermediated recombination xoff = pxp1 + (1− p)xp2 (8) with p = −0.25, producing a set of No = Nt/4 offspings, the fitness function is f = ∣∣Y DFT − Y TB ∣∣ (9) where Y is the band structure shown in Fig. 12(a) Figure 11. Schematics of a genetic algorithm. In Fig. 12, we show a comparison between (a) the band structure and (b) density of states of our fitted tight-binding Hamiltonian and the first-principle results. In the inset of Fig. 12 we also compare our results with a fit that considers an in- dependent orbital for hydrogen39. The fit obtained with a ge- netic algorithm for a graphene-only Hamiltonian is compara- ble with the one with independent orbitals and has the advan- tage (if compared with other multiparametric approaches) of being a semi-automatic procedure and using a reduced Hilbert space. B. Effects of hydrogenation in the quantum Hall effect. The transport coefficients were computed using a real-space O(N) method based on the Chebyshev expansion of a variant of the Kubo formula, the Kubo-Bastin formula40 σαβ(µ, T ) = i~ Ω ∫ ∞ −∞ dεf(T, µ, ε) (10) × Tr 〈 jαδ(ε−H)jβ dG+(H, ε) dε − jα dG−(H, ε) dε jβδ(ε−H) 〉 , 7 Γ K M K -1.5 -1 -0.5 0 0.5 1 1.5 E n er g y - E F ( eV ) Γ K M K -0.05 0 0.05 E n er g y - E F ( eV ) 0 5 10 -1.5 -1 -0.5 0 0.5 1 1.5 Figure 12. (a) Band structure showing the conduction (blue), va- lence (red) and midgap (black) bands and their corresponding (b) Density of states obtained through DFT (dashed line) and the fit- ted TB hamiltonian (thick line) for a hydrogen adatom on a 5× 5 graphene’s supercell. Inset: comparison between are graphene-only model, the DFT result and the model of ref. 39 (yellow line) for the midgap band. The parameters of the fit are t = −2.56 eV , t2 = 0.010580 eV, εI = −1.5694 eV, tI = 2.6538 eV and t (2) I = 0.29617 eV where δ(ε−H) the δ-function operator, jα the α-component of the current operators defined as jα ≡ (1/i~)[xα, H], G+(H, ε) and G−(H, ε) the advanced and retarded Green’s functions and f(T, µ, ε) the Fermi-Dirac distribution. In this method, the Green’s functions and the δ-functions are numer- ically calculated using the kernel polynomial method15,40–43. The magnetic field was incorporated by following Peierls’s substitution Hi,j = Hi,j(B = 0)eiφi,j , with φi,j =∫Rj Ri A(r) · dr the Peierl’s phase and A the vector potential, which was chosen using the Landau Gauge A = (B0y, 0, 0) We proceed to analise the effect of disorder in the QHE in graphene with SOC. Figure 13(a) shows de density of states for λR = 20 meV, λVZ = 30 meV - which are the same val- ues of the SOC of Figure 9(a)- for B = 13T and different concentrations of hydrogen. The spin-split LL states are still visible with 0.1% of hydrogen, that is the concentration that can be expected from contamination of the samples in exper- iments. Still, for 1% of hydrogen, the gap at E = 0 is closed and the LL spectrum is destroyed. 13(b) presents longitudinal and transverse conductivities for 0.1% of hydrogen. Although not all quantum Hall plateaus are visible in the presence of disorder, the Laudau levels, including the spin-split states, are still visible in the longitudinal conductivity. This indicates that our analysis, based on the scaling of spin-split states from longitudinal conductivity data, should be effective even in the presence of hydrogenation that produces both intra-valley and inter-valley scattering. 0 0.1 0.2 D O S -0.2 -0.1 0 0.1 0.2 Energy ( t) -8 -6 -4 -2 0 2 4 6 8 σ xx (e 2 /h ) -8 -6 -4 -2 0 2 4 6 8 σ xy (e 2 /4 h) Figure 13. (a) Density of states for graphene under the action of an external magnetic field of strength B = 13 T, with a proximity- induced SOC defined by the parameters λR = 20 meV, λa = −λb = 30 meV, for different concentration of hydrogen adatoms: xp=0 (dashed blue), 0.001 (solid black) and 0.01 (dotted green). In all cases, there is also a small Anderson disorder with width W = 5meV. (b) Dissipative conductivity (circle Black) and Hall conductivity (thick red) for graphene under the action of an external magnetic field of strength B = 13 T, with a proximity-induced SOC defined by the parameters λR = 20 meV, λVZ = 30 meV with 0.1% hydrogenated impurities V. CONCLUSIONS Inferring the type and strength of spin-orbit coupling in a graphene heterostructure is a difficult experimental task. Quantum Hall measurements provide a useful tool to under- stand the physics that govern the charge carriers in 2D ma- terials. The Landau spectrum offers hints about the low en- ergy physics of electrons (or holes) and the possible SOCs. Each Landau level in pristine graphene is eightfold degener- ate due to spin, pseudospin and isospin quantum numbers. In the presence of SOC, this degenerescence is lifted in a specific way depending on the type of coupling induced in graphene. Motivated by experimental results on magneto- transport in graphene heterostructures, we studied QHE on graphene with different SOCs and provided a simple way to characterize SOCs induced by proximity effect by analyzing the spin-splitting of the Landau Levels. By using two dif- ferent scaling laws, we were able to determinate what was the dominant contribution to the SOC in a graphene layer. Additionally, we used an efficient genetic algorithm strategy together with ab-initio calculations to obtain a realistic all- graphene tight-binding Hamiltonian that models hydrogena- tion in graphene. With this novel Hamiltonian and a quantum transport approach, we analysed the effect of hydrogenation on the QHE in graphene with interface-induced SOC. The nu- merical results indicate that the scaling laws can in principle 8 be applied even with 0.1% of hydrogenation. All results presented here are based on a free parti- cle approximation. This approximation is appropriate in most situations but with evolution of the manufacturing of graphene, cleaner samples have been obtained and the effects of Coulomb interactions may become important. Due to the Dirac-Weyl nature of electrons on graphene, the role of in- teractions are different from typical materials, and still under intensive discussion44. The signature of Coulomb interaction for the quantum Hall physics of graphene is the appearance of extra four-fold splittings in the energy spectra, implying in- termediate plateaus in the Hall conductivity. For low Landau levels, this splitting has already been observed experimentally with high magnetic fields45,46. The importance of Coulomb interaction increases exponentially with magnetic field47 and the symmetry breaking is more pronounced for lower Landau levels, high magnetic fields and high mobility48. The quality of sample is a important factor to observe Coulomb splitting as disorder decreases the mobility of electrons. Approaches based on a free particle picture cannot capture the physics of Coulomb interactions, what limits our calculations for sam- ples with moderate electron mobilities and moderate magnetic fields (typically smaller than 15 Tesla). The study of effect of Coulomb interaction on the Landau spectra is beyond scope of present work. However, it is important to mention that for experimental transport measurements, the analysis presented here shows that a two-fold splitting of LLs can be associated to the presence of spin-orbit coupling while a four-fold splitting is a signature of Coulomb interaction. VI. ACKNOWLEDGEMENTS Tatiana G. Rappoport acknowledges the support from the Newton Fund and the Royal Society (U.K.) through the Newton Advanced Fellowship scheme (ref. NA150043), Tarik Cysne, A. R. R and T. G. R. thank the Brazil- ian Agency CNPq for financial support. JHG acknowl- edge the Severo Ochoa Program (MINECO, Grant SEV- 2013-0295), the Spanish Ministry of Economy and Com- petitiveness (MAT2012- 33911), the Secretaria de Univer- sidades e Investigacion del Departamento de Economa y Conocimiento de la Generalidad de Catalua, and the European Union Seventh Framework Programme under grant agree- ment 604391 Graphene Flagship.This research used computa- tional resources from the Santos Dumont supercomputer at the National Laboratory for Scientific Computing (LNCC/MCTI, Brazil). Appendix A: Eigenvalues and Eigenstates The hamiltonians Hξ is block diagonal with each block indexed by valley occupation number nξ (eigenvalue of the number operator N̂ξ = a†ξaξ). The two lowest blocks are 1× 1 and 3× 3 matrices and higher blocks are 4× 4 matrices for both valleys. The eigenvectors of first blocks are ∣∣ψ+ 0,0 〉 = ∣∣0, B,−〉 (A1)∣∣ψ−0,0〉 = ∣∣0, A,−〉 (A2) with energies E+ 0,0 = λb −∆ (A3) E−0,0 = λa + ∆. (A4) The first excited block of each valley is given by H±1 =  ε±α 0 ~ω 0 ε±β ∓2iλR ~ω ±2iλR ε±γ  . (A5) in the basis of |1, B(A),−〉, |0, B(A),+〉 and |0, A(B),−〉 for ξ = +(−), where ε±α = λb(a) ∓∆, ε±β = −λb(a) ∓∆ and ε±γ = −λa(b) ±∆. The eigenvectors, indexed by i = 1, 2, 3, are given by∣∣ψ±1,i〉 = α±1,i ∣∣1, B(A),− 〉 + β±1,i ∣∣0, B(A),+ 〉 +γ±1,i ∣∣0, A(B),− 〉 (A6) (A7) with coefficients α±1,i = ( (2λR)2 − (ε±β − E ± 1,i)(ε ± γ − E±1,i) ) / √ D±,(A8) β±1,i = ±i(2λR)(~ω)/ √ D±, (A9) γ±1,i = (~ω)(ε±β − E ± 1,i)/ √ D±, (A10) whereD± = (~ω)2(ε±β −E ± 1,i) 2+((2λR)2−(ε±β −E ± 1,i)(ε ± γ − E±1,i)) 2 + (~ω)2(2λR)2. Finally, for nξ > 2 the the blocks, written in the basis |n±− 2, A(B),+〉, |n±, B(A),−〉, |n± − 1, B(A),+〉 and |n± − 1, A(B),−〉, are given by H±n =  E±a 0 ~ω √ n± − 1 0 0 E±b 0 ~ω√n± ~ω √ n± − 1 0 E±c ∓2iλR 0 ~ω√n± ±2iλR E±d  where E±a = λa(b)±∆, E±b = λb(a)∓∆, E±c = −λb(a)∓∆ and E±d = −λa(b) ±∆. The eigenstates are given by ∣∣ψ±n,i〉 = a±n,i ∣∣n± − 2, A(B),+ 〉 + b±n,i ∣∣n±, B(A),− 〉 +c±n,i ∣∣n± − 1, B(A),+ 〉 + d±n,i ∣∣n± − 1, A(B),− 〉 (A11) and their coefficients can be written in terms of four eigenval- ues of each block E±n,i : 9 b±n,i = ( 1 + ( E±b − E ± n,i )2 n± ( ~ω )2 + n± ( ~ω 2λR )2( 1 + (~ω)2(n± − 1) (E±a − E±n,i)2 )( 1− ( E±b − E ± n,i )( E±d − E ± n,i )( ~ω )2 n± )2)−1/2 , (A12) a±n,i = − i ( ~ω )2√ n± ( n± − 1 ) (±2λR)(E±a − E±n,i) ( 1− ( E±b − E ± n,i )( E±d − E ± n,i )( ~ω )2 n± ) b±n,i, (A13) c±n,i = i~ω√n± (±2λR) ( 1− ( E±b − E ± n,i )( E±d − E ± n,i )( ~ω )2 n± ) b±n,i, (A14) d±n,i = − (E±b − E ± n,i) ~ω√n± b±n,i. (A15) - 1 Martin Gmitra, Denis Kochan, Petra Hgl, and Jaroslav Fabian, Phys. Rev. B 93 , 155104 (2016) 2 M. Venkata Kamalakar, André Dankert, Johan Bergsten, Tommy Ive, and Saroj P. Dash, Appl. Phys. Lett. 105, 212405 (2014) 3 P. J. Zomer, M. H. D. Guimarães, N. Tombros, and B. J. van Wees, Phys. Rev. B 86, 161416 (2012) 4 C. R. Dean, L. Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao, J. Katoch, M. Ishigami, P. Moon, M. Koshino, T. Taniguchi, K. Watanabe, K. L. Shepard, J. Hone and P. Kim, Nature (London) 497, 598 (2013) 5 B. Hunt, J. D. Sanchez-Yamagishi, A. F. Young, M. Yankowitz, B. J. LeRoy, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, R. C. Ashoori, Science 340, 1427 (2013) 6 W. Han, R. K. Kawakami, M. Gmitra and J. Fabian, Nat. Nan- otechnol. 9, 794 (2014). 7 A. K. Geim and I. V. Grigorieva, Nature (London) 499, 419 (2013). 8 Sergej Konschuhm, Martin Gmitra and Jaroslav Fabian, Phys. Rev. B 82 24245412 (2010). 9 Hongki Min, J. E. Hill, N. A. Sinitsyn, B. R. Sahu, Leonard Klein- man, and A. H. MacDonald, Phys. Rev. B 74 165310 (2006). 10 Yugui Yao, Fei Ye, Xiao-Liang Qi, Shou-Cheng Zhang, and Zhong Fang, Phys. Rev. B (R) 75 041401 , (2007). 11 A. H. Castro Neto and F. Guinea, Phys. Rev. Lett. 103, 026804 (2009). 12 A. Ferreira, T. G. Rappoport, M. A. Cazalilla, and A. C. Neto, Phys. Rev. Lett. 112, 066601 (2014). 13 H.-Y. Yang, C. Huang, H. Ochoa, and M. A. Cazalilla, Phys. Rev. B 93, 085418 (2016). 14 Mirco Milletari and Aires Ferreira, Phys. Rev. B 94, 134202 (2016). 15 J. H. Garcia and T. G. Rappoport, 2D Materials 3, 024007 (2016). 16 Jose Hugo Garcia Aguilar, Aron W. Cummings, and Stephan Roche, Nano lett., 17, 5078 (2017). 17 Zhiyong Wang, Chi Tang, Raymond Sachs, Yafis Barlas, Jing Shi, Phys. Rev. Lett. 114, 016603 (2015). 18 J. B. S. Mendes, O. Alves Santos, L M. Meireles, R. G. Lacerda, L. H. Vilela-Leao, F. L. A. Machado, R. L. Rodriguez-Suarez, A. Azevedo, and S. M. Rezende, Phys. Rev. Lett. 115, 226601 (2015). 19 A. Avsar, J.Y. Tan, T. Taychatanapat, J. Balakrishnan, G.K.W. Koon, Y. Yeo, J. Lahiri, A. Carvalho, A.S. Rodin, E.C.T. O’Farrell, G. Eda, A.H.C. Neto, and B. Ozyilmaz, Nat. Commun. 5, 4875 (2014). 20 Z. Wang, D.-K. Ki, H. Chen, H. Berger, A.H. MacDonald, and A.F. Morpurgo, Nat. Comm. 6, 8339 (2015). 21 Zhe Wang, Dong-Keun Ki, Jun Yong Khoo, Diego Mauro, Hel- muth Berger, Leonid S. Levitov, and Alberto F. Morpurgo, Phys. Rev. X 6, 041020 (2016). 22 B. Yang, M.-F. Tu, J. Kim, Y. Wu, H. Wang, J. Alicea, R. Wu, M. Bockrath and J. Shi , 2D materials 3, 031012(2016). 23 E. C. T. O’Farrell, J. Y. Tan, Y. Yeo, G. K. W. Koon, B. Ozyil- maz, K. Watanabe, and T. Taniguchi, Phys. Rev. Lett. 117, 076603 (2016). 24 Yunqiu Kelly Lou, Jinsong Xu, Tiancong Zhu, Guanzhong Wu, Elizabeth J. McCormick, Wenbo Zhan, Mahesh R. Neupane, and Roland K. Kawakami, Nano Lett., 17, 3877 (2017). 25 A. Avsar, D. Unuchek, J. Liu, O. L. Sanchez, K. Watanabe, T. Taniguchi, B. Ozyilmaz, and A. Kis, arXiv:1705.10267 (2017). 26 Y. A. Bychkov and E. I. Rashba, J. Phys. C 17, 6039 (1984). 27 E. I. Rashba, Phys. Rev. B 79, 161409(R) (2009). 28 C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005). 29 C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005). 30 G. D. Mahan, Many-Particle Physics, 3rd ed. (Plenum,New York, 2010). 31 N. M. R. Peres, Rev. Mod. Phys. 82, 2673 (2010). 32 V. M. Pereira, J. M. B. Lopes dos Santos, and A. H. Castro Neto, Phys. Rev. B 77, 115109 (2008). 33 M. Gmitra, D. Kochan, and J. Fabian, Phys. Rev. Lett. 110, 246602 (2013). 34 P. Hohenberg and W. Kohn, Phys. Rev. 136 B864 (1964). 35 W. Kohn and L. J. Sham, Phys. Rev. 140 (4A) (1965). 36 J. M. Soler, E. Artacho, J. D. Gale, A. Garcı́a, J. Junquera, P. Ordejón, and D. Sánchez-Portal, J. Phys.: Cond. Matt. 14, 2745 (2002). 37 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996). 38 N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991). 39 M. Gmitra, D. Kochan, and J. Fabian, Phys. Rev. Lett. 110, 246602 (2013). 40 J. H. Garcia, L. Covaci, and T. G. Rappoport, Phys. Rev. Lett. 114, 116602 (2015). 41 A. Weisse, G. Wellein, A. Alvermann, and H. Fehske, Rev. Mod. Phys. 78, 275 (2016). http://arxiv.org/abs/1705.10267 10 42 L. Covaci, F.M. Peeters, M. Berciu, Phys. Rev. Lett. 105, 167006 (2010) 43 R. Silver and H. Rder, Int. J. Mod. Phys. C 5, 735 (1994). 44 V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea, and A. H. Castro Neto, Rev. Mod. Phys. 84, 1067 (2012) 45 Y. Zhang, Z. Jiang, J. P. Small, M. S. Purewal, Y.-W. Tan, M. Fazlollahi, J. D. Chudow, J. A. Jaszczak, H. L. Stormer, and P. Kim, Phys. Rev. Lett. 96, 136806 (2006). 46 A. F. Young, C. R. Dean, L. Wang, H. Ren, P. Cadden-Zimansky, K. Watanabe, T. Taniguchi, J. Hone, K. L. Shepard, and P. Kim, Nat. Phys. 8, 550 (2012). 47 M. O. Goerbig, Rev. Mod. Phys. 83, 1193 (2011). 48 K. Nomura and A. H. MacDonald, Phys. Rev. Lett. 96, 256602 (2006). Quantum Hall Effect in Graphene with Interface-Induced Spin-Orbit Coupling Abstract I Introduction II Theory III Analytical Results A Weak spin-orbit coupling B Strong Spin-Orbit coupling Regime C Possible applications to other systems IV Numerical Results A Graphene-only Hamiltonian for hydrogen on graphene. B Effects of hydrogenation in the quantum Hall effect. V Conclusions VI Acknowledgements A Eigenvalues and Eigenstates References