Energy spectrum and Landau levels in bilayer graphene with spin-orbit interaction

We present a theoretical study of the bandstructure and Landau levels in bilayer graphene at low energies in the presence of a transverse magnetic field and Rashba spin-orbit interaction in the regime of negligible trigonal distortion. Within an effective low energy approach (L\"owdin partitioning theory) we derive an effective Hamiltonian for bilayer graphene that incorporates the influence of the Zeeman effect, the Rashba spin-orbit interaction, and inclusively, the role of the intrinsic spin-orbit interaction on the same footing. Particular attention is spent to the energy spectrum and Landau levels. Our modeling unveil the strong influence of the Rashba coupling $\lambda_R $ in the spin-splitting of the electron and hole bands. Graphene bilayers with weak Rashba spin-orbit interaction show a spin-splitting linear in momentum and proportional to $\lambda_R $, but scales inversely proportional to the interlayer hopping energy $\gamma_1$. However, at robust spin-orbit coupling $\lambda_R $ the energy spectrum shows a strong warping behavior near the Dirac points. We find the bias-induced gap in bilayer graphene to be decreasing with increasing Rashba coupling, a behavior resembling a topological insulator transition. We further predict an unexpected assymetric spin-splitting and crossings of the Landau levels due to the interplay between the Rashba interaction and the external bias voltage. Our results are of relevance for interpreting magnetotransport and infrared cyclotron resonance measurements, including also situations of comparatively weak spin-orbit coupling.


Introduction
Graphene and its bilayer (BLG) posses very distinctive physical properties. [1,2] Neglecting interactions, the low energy quasiparticles in pristine single and double layer graphene obey, respectively, linear and quadratic dispersion laws at the K(K ′ ) symmetry points. [4,3,5] In the presence of a quantizing magnetic field B > 0, the relativistic massless Dirac fermions of graphene exhibits Landau levels (LLs) non-equidistant in energy, E n ∝ |n|B, with n = 0, ±1, ±2, ..., the LL index. [6]The latter gives rise to the half integer quantum Hall effect at room temperature, [7] and to the fractional quantum Hall effect at high magnetic fields in the presence of many body effects. [8,9] On the other hand, LLs for BLG show a rather intricate index sequence instead; with a roughly linear B-field dependence for low LLs, and a √ B for high LLs. [10,11,12]At low energies, the LLs in unbiased BLG follows the sequence E n ∝ |n|(|n| + 1) for n ≥ 1 with a double degenerate zero-energy level, E o = 0 for n = 0. [13,14,15] Experimentally, the LLs dipole-allowed transition energies in single layer and BLG have been studied in detail by cyclotron resonance. [16,17,18] Most recently, phonon-mediated inter LL excitations have been explored by magneto-Raman scattering experiments. [19] The particularly high interest in graphene spin physics is strongly motivated by its expected prospects in nanoelectronics and spin-based devices for spintronics. [20] In this realm, the role of the spin-orbit interaction (SOI) effects in graphene sheets is a phenomenon under intense scrutiny. Two types of SOI in graphene have been identified; (i) the induced by carbon intra-atomic SOI (intrinsic-SOI), and (ii), the SOI coming from the breaking of the space inversion symmetry of the hexagonal lattice (extrinsic-SOI), this is the so-called Rashba-SOI. The latter can have different origins, among those is the presence of a substrate, buckling, ad-atoms, or external electric fields.
The magnitude of the excitation gap η of the intrinsic-SOI in single layer is predicted to be very small (0.86-50 µeV). [21,22,23,24,25,26,27] Likewise, estimates of the Rashba coupling λ R , leads to small values (few tens of µeV) at typical electric fields (∼ 0.16 mV/nm). [23,24,25] However recent spin-resolved photoemission experiments in graphene/Au/Ni(111) showed enhancements of the Rashba coupling as large as 13 meV. [28] Induced distortions by neutral impurities (ad-atoms) [29,30] and the interplay of buckling with Rashba-SOI also yield significant enhancement of the spin-splittings (up to 40 meV). [31] The role of the impurities and lattice deformations seems to be crucial for the observed long spin relaxation times (up to 2 ns in BLG [32]) linked to SOI effects [33,34], as well as for a non-monotonic dependence on bias voltage of the spin relaxation times in BLG due to the D'yakonov-Perel' spin-precession mechanism. [35] Large spin-splittings (∼ 0.22 eV) of the graphene π−states attributed to Rashba-SOI has been reported also in epitaxial graphene on a Ni substrate. [36] In theory, it has been shown that the Rashba-SOI can induce significant changes in the bandstructure of graphene [37,38] as well as in BLG. [39] In graphene monolayers with Rashba coupling, the Landau levels are described by E (1) n,µ± = µ E (1) n,± , (in units of λ R ), with [37] E (1) n,± = (2n − 1)Γ 2 + 1 2 1 ± 1 + 4(2n − 1)Γ 2 + 4Γ 4 (1) for n ≥ 2, whereΓ = v F /l B |λ R |, v F is the Fermi velocity (∼ 10 6 m s −1 ), l B = c /eB is the magnetic length, the Plank's constant over 2π, c is the light velocity in vacuum, and −e is the electron charge. Here µ = ± stands for the electron/hole LL branch. The lowest n = 0 is given by E 0 = 0 whiles the n = 1 level gives rise to three modes. A zero mode E (0) 1 = 0, in addition to its two nondegenerate modes at E 1µ = µ 1 + 2Γ 2 . Exact solutions for the LLs in monolayer graphene under the influence of a Zeeman field and spin-orbit interactions has been also reported recently by De Martino et al.. [40] In this paper we show within low energy effective theory that for biased BLG in which the Rashba effect is the dominant SOI, its LLs must follow E (2) n,µ± = µ E (2) n,± , with for n ≥ 2, with U the interlayer bias energy, Γ = 2 √ 2λ R v F /γ 1 l B , and ω = 2v 2 F 2 /γ 1 l 2 B , being γ 1 inter-layer hopping energy. As it occurs in graphene, in BLG with Rashba-SOI we obtain three modes for n = 1, namely, the nondegenerate E 1− = −U, whereas for n = 0 its eigenvalue also vanishes in the absence of gating, E 0− = U. Eq.(2) comprises one of the main results of this contribution.
The aim of this study is to investigate the energy spectrum and the Landau levels in bilayer graphene under the influence of sizable spin-orbit interactions (SOI) of both, intrinsic and Rashba types. An effective low energy Hamiltonian for bilayer graphene that includes both types of SOI and Zeeman effect is derived within Löwdin partitioning theory. Whiles the Rashba SOI in single layer graphene is known to modify its otherwise conic spectrum, to a spectrum that includes two zero gap bands and two gapped branches of width 2λ R (with parabolic shape, similar to unbiased BLG); [37] here in contrast, the (unbiased) BLG with Rashba-SOI shows a spin-splitting which is linear in momentum and proportional to λ R , but inversely proportional to the interlayer hopping parameter γ 1 . We predict also a strong influence of the Rashba spin-orbit interaction in the warping of the low energy bandstructure of biased bilayer at comparatively weak spin-orbit coupling λ R . It is found that the bias-induced gap in bilayer graphene decreases as the Rashba strength coupling is increased. It is further predicted the rise of an unexpected asymmetric spin-splitting of the Landau levels due to the interplay among the Rashba coupling and the external bias voltage.
The remainder of the paper is organized as follows. In Sec. 2 we outline the model and derivation for the low energy BLG effective Hamiltonian. The Landau level spectrum in the presence of Rashba-SOI is discussed in Sec. 3. The band spectrum properties and the Landau levels of BLG with Rashba-SOI are analyzed in detail in Sec. 4 and Sec. 5. A summary of our results is given in Sec. 6 . Additionally, we provide three appendixes. In Appendix A we outline the derivation (Löwdin partitioning theory) of the low energy Hamiltonian in the presence of intrinsic and Rashba types of SOI, as well as the Zeeman effect. In Appendix B a detailed description of the eigenvalues of the low energy Hamiltonian is given, and finally in Appendix C the Landau levels for BLG with Rashba-SOI are derived.

Low energy effective Hamiltonian for bilayer graphene
Here we focus in the derivation of the low energy effective Hamiltonian for BLG with SOI in the presence of magnetic field that eventually leads to Eq.(2). We start by considering a pile of two graphene layers (BLG) in which the sites A 2 of the upper layer 2 lies directly on top of the B 1 sites of the bottom layer 1 (AB-Bernal stacking). At the vicinity of the Dirac K symmetry point, the effective noninteracting bilayer graphene Hamiltonian H o , written in terms of the spin-dependent basis where π = π x + iπ y , with π = p + eA/c = (π x , π y ) is the canonical momentum, and A is the vector potential. Here γ ≡ v F = γ o a √ 3/2 , with γ o ∼ 3.1 eV (intralayer hopping energy) and a = 0.246 nm is the lattice parameter. [1] The electrostatic potential ±U of the bottom/upper layer is gate voltage adjustable and opens a gap of 2U in the spectrum. [2] The dominant interlayer interaction in BLG is described to first approximation by the term where γ 1 is the nearest neighbor (interlayer) hopping energy (∼ 0.1γ o ). The terms proportional to the velocities v 3 = γ 3 a √ 3/2 and v 4 = γ 4 a √ 3/2 arise due to second nearest neighbor (interlayer) hopping processes associated with γ 3 and γ 4 tightbinding parameters, respectively. [15,41] The former produces a trigonal warping whose characteristic energy is E 3 = 1 2 m * v 2 3 = γ 1 (γ 3 /2γ o ) 2 ∼ 1 meV, being m * = γ 1 /2v 2 F the electron effective mass [15], whiles the latter has yet a smaller characteristic energy, Since typically E 3,4 ≪ γ 1 in BLG, the contributions in Eq.(4) coming from the terms with velocities v 3 and v 4 are negligible. However we should notice that, when considering Rashba spin-orbit interaction, there might be situations in which v 3 and v 4 may be important, particularly at very weak Rashba-strengths of the order of E 3,4 . We would like to mention nevertheless, that there is experimental (and theoretical) evidence of Rashba-SOI spin-splittings of one order of magnitude [28,29,31] and even larger [36] than the typical distortion energies E 3,4 . In any case, its inclusion in the model can be readily incorporated in the general derivation of the reduced effective Hamiltonian via the Löwdin perturbation theory described in Appendix A. The regime in which there is a possible interplay among the v 3,4 terms and Rashba SOI is out of the scope of the present analytical study and for sake of clarity and simplicity these trigonal terms will be disregarded hereafter.
When taking into account explicitly the presence of intrinsic-SOI, Rashba-SOI and Zeeman splitting, the total eight-band effective Hamiltonian will read The second term to the right in Eq.(5) arises due the influence of an effective electric field perpendicular to the BLG plane producing a Rashba type of SOI. The the leading contribution to the Rashba-SOI Hamiltonian, H R , is modeled as follows: here λ R parameterize the strength of the intra-layer Rashba-SOI, as in monolayer graphene, with σ ± = 1 2 (σ x ± iσ y ), being (σ x , σ y ) the usual 2 × 2 Pauli spin matrices. The intensity of the Rashba-SOI can be sizable (λ R ∼ 10 meV) due for instance to the presence of a metallic substrate. [28] Within tight-binding theory, it is understood that the Rashba-SOI arises because of the effective nearest-neighbor hopping of two p z orbitals with opposite spins under the presence of an applied transverse electric field. [26] Recently it has been predicted that λ R can be even larger (of a few tens of meV) due to buckling effects in conjunction with external electric fields. [31] Furthermore by varying the electric field the Rashba parameter can be tuned. A possible inter-layer Rashba spin-orbit coupling of strength λ ⊥ R can in principle be present in BLG as well, however such contributions will be ignored here because of its predicted slight influence on the energy bands for λ ⊥ R /γ o 0.3. [39] Additionally, in the same basis set above, the intrinsic-SOI Hamiltonian for BLG (third term to the r.h.s. in Eq.(5) ) shall follows the 8 × 8 matrix form with η the intrinsic SOI constant and s z is the spin operator along the z-axis, perpendicular to the BLG plane. The intrinsic SOI is a second order tunneling process (within tight-binding theory), which involves next-nearest-neighbor hopping events of p z electrons of a given spin. As mentioned in the introduction, the value of the excitation energy η has been inferred to be very small in monolayer graphene in both K(K ′ ) symmetry points ( 50 µeV), even if one consider interactions up to first order by including the unoccupied d and higher orbitals. [25,26] Interestingly, in BLG, taking into account the interlayer overlapping of the π and σ bonds yields enhanced values of the intrinsic-SOI; about one order of magnitude larger than in single layer graphene ( ∼ 0.1 meV). [42] We would like to emphasize here that such values still somewhat weak, compared with those relatively large strengths, of which reportedly, the Rashba parameter λ R can acquire (similar to the values attained in III-V semiconductors). Nevertheless, for generality, the intrinsic SOI has been incorporated in the present derivation of the low energy effective Hamiltonian. This will be helpful when considering BLG in the extreme limit, i.e. when the intrinsic-SOI η is much stronger than the Rashba-SOI λ R parameter (η ≫ λ R ). However, we shall concentrate here our discussion on the bandstructure and Landau levels of BLG for the case λ R ≫ η, the opposite limit will be treated elsewhere. The last term in Eq.(5) arises if an external magnetic field B is present, affecting the energetic of the quasiparticles in the form of a Zeeman interaction, H Z , for a field perpendicular to the BLG plane it reads where I is a 4 × 4 unit matrix, 2∆ = gµ B B is the Zeeman splitting energy, g is the electron Landé factor, and µ B is the Bohr magneton. Note that even at relatively high magnetic fields (B = 10 T), the Zeeman splitting it is still somewhat small (∆ ∼ 1.1 meV), whiles it is practically negligible at low fields. Note that the condition, typically holds at finite fields (B 0.1 T). This condition will allow us to work safely within the low energy theory and derive an effective Hamiltonian for BLG, including the extrinsic (Rashba)-SOI, intrinsic-SOI, and Zeeman effect on the same footing. We finally should remark that the total Hamiltonian (5) is valid near K symmetry point only. For the K ′ point of the Brillouin zone, the Hamiltonian H K ′ = Σ y H K Σ −1 y , with Σ y = σ y ⊗ I, should be used instead.

Low energy bilayer Hamiltonian
Using Löwdin partitioning theory [44,45] the full 8×8 Hamiltonian H K can be projected through a canonical transformation [46] into a 4 × 4 low energy effective Hamiltonian H in an appropriate basis (see Appendix A). The projected low energy Hamiltonian can be further expressed in terms of Kronecker products of 2 × 2 matrices and conveniently separated into the sum of the Hamiltonians (keeping terms up to 1/γ 2 1 ), in which the term independent of the interlayer hopping parameter γ 1 reads, where σ z is the z−component of the Pauli matrices and σ o is the 2 × 2 unit matrix. The dominant contribution to H is described by the Hamiltonian H (1) , given by where we have defined the operator s ± = 1 2 (σ o ± σ z ). Without Rashba-SOI (λ R = 0), Eq.(12) decouples to the usual effective BLG Hamiltonian obtained within low energy theory in the absence of trigonal warping effects. [13] Such term gives rise to well known parabolic spectrum of the massive Dirac quasiparticles in BLG. The second term in H (1) is linear in momentum and can be viewed as a renormalization of the Rashba coefficient due to the presence of the higher bands. Notice that it scales as the inverse of the interlayer hopping energy γ 1 .
The remaining terms proportional to 1/γ 2 1 in Eq.(10) are compacted into the sum The low energy effective Hamiltonian H described in Eq. (10) is valid within the energy range ǫ γ 1 . Notice that it will be fairly sufficient to keep only the leading order contribution h (2) 1 in H (2) given the typical smallness of the ratios λ 2 R /γ 2 1 , λ R ∆/γ 2 1 , and λ R η/γ 2 1 appearing in H (2) together with the assumption U < γ 1 . Hence the description for the low energy (and momentum) effective Hamiltonian will be given by If we further neglect the Zeeman and the intrinsic SOI (∆ = η = 0) the effective bilayer Hamiltonian with Rashba-SOI written in the atomic basis where we have defined the parameters ξ = 2Uγ 2 , α = 2γ λ R and β = γ 1γ 2 , with γ = γ/γ 1 .

Bilayer graphene spectrum with Rashba effect at zero field
Without magnetic field (B = 0), π † = k − = (k x − ik y ) and π = k + = (k x + ik y ) and the eigenvalues of Eq.(13) are readily determined by (Appendix B) with k = k 2 x + k 2 y . Here µ = ± describe the electron/hole branch, whiles s = ± characterizes its spin chirality. Therefore, the low-energy spectrum consist of four spinsplitted bands, two conduction and two valence bands. For the unbias voltage case (U = 0) the spectrum reduces simply to We note that in contrast with single layer graphene, [47] in BLG the Rashba-SOI it is induced a linear spin-splitting in momentum of the conduction and valence bands, in close analogy with the Rashba interaction arisen in two-dimensional electron gases in semiconductor heterostructures. In addition we observe that the cyclotron effective mass at the Fermi energy m (s) µs /∂ k) turns out to be spin-dependent as long as λ R = 0 following the relation, with λ β = β/α = γ/2λ R , and we have expressed the Fermi wave number in terms of the 2D carrier density via k F = √ πn. Notice that in the limit case λ β → ∞, the cyclotron mass m (s) c → µm * = µγ 1 /2v 2 F , as one expects for unbiased BLG in the absence of SOI. However as the carrier density n → 0 the cyclotron effective mass does not diverges as predicted by tight-binding models.
The normalized eigenvectors |ψ of H corresponding to the electron and hole bands µ = ±, respectively, for the case of spin up (s = +) are written in the four-vector form (Appendix B) as, whiles the normalized eigenvectors for the spin down (s = −) electron/hole bands are specified by in which we have defined the dimensionless parameter with θ = tan −1 (2βk/α), and φ is the azimuthal angle of the in-plane wave vector, k = k(cos φ, sin φ). The denominator of Eq. (19) is explictly, R µs = µ ε o µs (k) = |ε o µs (k)|, which implies R +s = R −s , and therefore the relation χ (+) σ χ (−) σ = −1 it is always satisfied. Without external bias voltage (U = 0), the parameter χ (±) σ reduces to ±1 for all k.
τ y µs = 2χ whereas the components of the spin polarization (in units of /2) satisfy S y µs = + s sin θ cos(φ), As it occurs with the standard Rashba SOI in semiconductors, in BLG the orientation of the spin-polarization in the plane is always perpendicular to the direction of the momentum, S · k = 0. We notice also that, in contrast with single layer graphene, in BLG the dot product S · τ = 0 in general. Interestingly that is, S µs is forced to lay on the BLG plane. Clearly, as the Rashba SOI coefficient λ R → 0, i.e. α → 0, the magnitude of the spin-polarization reaches it maximum value, | S µs | → 1 as the electron/hole spin is conserved. Finally, we notice that Eq. (26) for the spin orientation in unbiased BLG is formally identical to that obtained in single layer graphene with Rashba SOI. [37] 3. Landau levels in bilayer graphene with Rashba SOI For a magnetic field B = 0 perpendicular to the BLG plane, the operators π and π † do not commute any more since its components fails to do so, and of course, care has to be exercised in their ordering. By making the substitution in Eq. (13) of the momentum operators in terms of the Bose operators, π † = √ 2 a † /l B and π = √ 2 a/l B with a, a † = 1, the effective Hamiltonian in the limit of low bias (U ≪ γ 1 ) takes the form with the notation Γ = √ 2 α/l B and ω = 2 2 β/l 2 B . The eigenfunctions ofĤ L can now be written in the form, |ψ n = (c and a|n = √ n|n − 1 . Consequently one can write the expectation value ψ n |Ĥ L |ψ n = φ n |H n |φ n , with and |φ n = (c ) T satisfying the normalization condition φ n |φ n = 1, whiles The eigenvalues of (28) leads to the Landau spectrum given by Eq. (2), which, in the absence of a bias gate voltage across the layers (U = 0), reads (Appendix C) for n ≥ 2. In order to gain further physical insight of the behavior of the Landau levels in BLG, it is illustrative to analyze the limit cases at zero, weak(large) Rashba-SOI relative to the magnetic field strength, with and without bias voltage. (ii) Weak Rashba-SOI (Γ/ω ≪ 1). At very weak Rashba-SOI strength values (large fields), the LL's still evolve approximately linear with B, but shifted by a small energy proportional to λ 2 R , described by with Γ o = 2ω o /m * (λ R /v F ) and n ≥ 2. Since Γ 2 o /4ω o = λ 2 R /γ 1 , only for rather large Rashba SOI coefficient (λ R ≃ γ 1 ) strengths gives rise to significant broken degeneracies of the electron/hole LLs as described next in the opposite regime.
In addition, the energy spectrum with LL level index n = 1 and n = 0 are special cases, giving rise to three levels, one at zero energy (ǫ 0− = 0) and two nondegenerate levels for n = 1. In the high field limit, in the weak field regime.
(i) Zero Rashba-SOI in the ω ≪ U limit. In this case the quantum states should follows clearly, in addition to the aperture of an energy gap of 2U between the spin-degenerate electrons/hole LLs, the presence of the gate voltage induces a deviation from the linear dependence in B occurring at U = 0, to a parabolic behavior with B instead. This is also a known result in the literature. [11,12].
(ii) Weak Rashba-SOI in the regime U ≫ ω ≫ Γ. Here the LLs still behave quadratic in B but with a spin-dependent shift linear in B given as the term proportional to Γ 2 o is responsible for the anticrossings of the fan spectrum of the ν = ± states.
(iii) Strong Rashba-SOI in the regime U ≫ Γ ≫ ω. In contrast with the case of U = 0, at very large Rashba-SOI strengths the LL spectrum follows (to leading order) a linear behavior with B for the positive chirality states (ν = +), whereas for the negative (ν = −) quantum states develops a cubic dependence instead. Explicitly they are given by, such drastic change on the field dependence of the LLs with different spin chirality states (ν = ±) will enhance dramatically the degree of their spin-splitting and on the multiplicity of the level crossing as studied in the next section. All of this discussed above holds for n > 2. The cases n = 0, 1 are treated separately as in the condition for U = 0. In the limit U ≫ ω ≫ Γ (small Rashba SOI or large fields) ε 1µ+ ≃ µ(U + ω 2 o B 2 /U), and ε 1µ+ ≃ µ(U + Γ 2 o B/U) for the large Rashba-SOI (or small Rashba-SOI) case (U ≫ Γ ≫ ω), whiles ε 1− = −U for both regimes. The case for the zero Landau level index is ε 0− = U, which also holds in both regimes.

Spin-polarization of the Landau levels
The eigenfunctions of Eq. (28) for the n−th Landau level of a given electron(hole) band µ in biased BLG are given in Appendix C. From Eqs. (C.18) and the orthogonality . We notice that the valley polarization in the perpendicular direction turns to be k−independent, and that for the unbiased case, it vanishes for all LL. Furthermore, the z-th component of the spin polarization of the n−th LL with U = 0, reduces to where, as before ν = ± denotes the plus/minus n−LL of a given µ = ± electron/hole branch. In the limit of high field S (n) z µ± → ∓1, a full polarization is reached, and the states ν = ± coincides with the spin-magnetization signs (∓) of the LLs. If the limit B → 0 is taken, then the spin-polarization S (n)

Band structure properties in BLG with Rashba SOI
The low quasiparticle energy band structure for unbiased bilayer graphene with Rashba-SOI, as predicted by Eq.(15) at zero field, is illustrated in Fig.1 for different values of λ R strength. In the absence of Rashba-SOI (λ R = 0) the well recognized parabolic spin-degenerated conduction and valence bands touching at its extrema at k = 0 are depicted in Fig.1(a). For non-zero λ R the spin-degeneracy of the bands is broken inducing a k-linear spin-splitting of width ∆ s (k) = |E ±,∓ − E ±,± | = αk. For relatively weak Rashba-SOI (4β 2 k 2 ≫ α 2 ) the band dispersion follows a parabolic behavior, ε o µs (k) ≃ µ(β 2 k 2 − 1 2 sαk), as shown in Fig.1(b). On the other hand, if the condition 4β 2 k 2 ≪ α 2 holds, then the relation (15) evolves to a linear spectrum for the inner bands (ε o µs (k) ≃ µαk for s = −), and to a k-cubic spectrum for the outermost bands (ε o µs (k) ≃ µβ 2 k 3 /α for s = +) as plotted in Fig.1(d). Similar drastic changes in the bandstructure due to large Rashba-SOI strengths were reported earlier numerically by van Gelderen et al. [39] within a tight-binding framework; see for instance the low energy bands near the Dirac points of Fig.1(d) of that reference.
At the intermediate regime ( Fig.1(c)), the innermost bands interpolates from a k-linear behavior for electron/hole momentum very close to the Dirac point, to a kcubic dependence for high momentum. In contrast, the s = + bands seems to be well described by the cubic spectrum for all values of k. Such remarkable asymmetry behavior of the spectrum of BLG with Rashba-SOI is certainly unique, since it is not seen in monolayer graphene neither in semiconductors with the Rashba-type of SOI. These peculiar characteristics of the spectrum of BLG would have interesting consequences on the electronic and spin-transport properties.
In Fig.2 we show the bilayer graphene low energy spectrum for finite bias voltage (U = 0.050 eV) at various Rashba coupling strengths (λ R = 0, λ o , 4λ o and 12λ o ). As in the unbias case, without Rashba-SOI, the bands are spin degenerated (Fig.2(a)). A  gap of 2|U| at k = 0 is opened between the conduction and valence bands turning BLG into a semiconductor. Moreover a band-bending appears at small wave-numbers (low momentum) due interplay with the bias gate voltage U. This is the so-called "Mexicanhat-like" shape of the lowest energy bands well reported in the literature. From Eq.(B.8), in this regime (λ R = 0) and ka 1 the bands are reasonably well described by ε µs (k) ≃ µ(U − ξk 2 + (β 2 /2U)k 4 ). For nonzero λ R (Figs.2(b)-(d)) the spin-degeneracy is lifted. As the Rashba-SOI parameter increase the lowest/highest conduction/valence bands becomes more warped and the gap tend to decrease as the lowest conduction band E ++ evolves from a Mexican-hat like shape to an inverted one, and vice versa for the highest valence band E −+ , see Fig.3(d) calculated using Eq.(A.14). The behavior of the gaps ∆ g+ = E ++ − E −+ (for the innermost bands) and ∆ g− = E +− − E −− (for the outmost bands) as a function of the ratio λ R /γ 1 is plotted in Fig.3 for different bias voltages U. The gap ∆ g+ in BLG closes as the Rashba parameter increases, reaching its minimum (zero gap) at λ R = γ 1 / √ 2, to then gradually and linearly opens again as λ R /γ 1 is increase up to 1. For λ R /γ 1 > 1 the gap ∆ g+ remains constant. Inset (a) of Fig.3 shows the bandstructure for biased BLG with U = 0.025 eV illustrating the zero gap condition. An analogous behavior can be seen in the numerical plot depicted in Fig.7(b) of Ref. [39]. Notice that our analytical low effective modeling for the bilayer graphene bandstructure allow us to capture also these anomalous behavior of the lowest bands (near the Dirac points) under finite bias and relatively large Rashba coupling. In addition it predicts that the closing of the gap occurs provided that λ R = γ 1 / √ 2, regardless of the magnitude of the bias voltage applied.

Landau level spectrum in BLG with Rashba-SOI
The Landau level energy spectrum as a function of magnetic field (up to 12 T) for BLG with Rashba-SOI according to formula Eq.(2) is plotted (from n = 0 to n = 41) in Fig.4 and Fig.5 for the unbiased and biased case, respectively. The zero gate voltage (U = 0) and without Rashba-SOI case shows a LL fan diagram which is linear with B and degenerate in spin ( Fig.4(a)), as expected, because of the parabolic behavior of the energy bands and since there is no spin-dependent mechanism here to break the spin-symmetry. When the Rashba-SOI is present the spin-degeneracy is lifted inducing multiple crossings of the LLs at the Fermi energy, similar as it occurs in semiconductor 2DEGs with Rashba SOI, and stems because for sizable λ R strengths the LLs with high index and with spin chirality s = + have lower energies than those with spin chirality s = − as the field is increased. This is basically an intermediate regime between those LLs discussed in section 3.1(ii) and 3.1(iii). Notice that for a relatively weak intensity of the Rashba parameter (λ R = λ o ) the LLs behave roughly linear with the field (Fig.4(b)), no matter its spin chirality. However, for large values (λ R = 4λ o ) a drastic and unusual change in the LLs spectrum arises (Fig.4(c)); the LLs with spin ν = + evolves as B 1/2 whiles for those with ν = − develops a B 3/2 dependence as described in section 3.1(iii). Such difference in the field dependence effectively squeezes the LLs with ν = − to lower energies as λ R is increased promoting multiple crossings between the LLs. This surprising result suggests a strong spin polarization of the Landau levels induced by a significant increase of the Rashba parameter.
The LLs for the biased U = 0 and without Rashba coupling case shows instead a parabolic behavior with B, and the fan diagram is split by a gap of 2U (Fig.5(a)) as also described by Eq. (32). The linear behavior observed for U = 0 and λ R = λ o of the LLs transforms as well to a B 2 dependence for U = 0, as predicted also by Eq. (33). The presence of the Rashba-SOI manifest as well as crossings of the LLs as seen in Fig.5(b). The behavior illustrated in Fig.5(c) at relatively large λ R is similar to that seen in Fig.4(c) but with an open gap and with almost a linear behavior with field for ν = +, and roughly a cubic dependence in B for the ν = − LLs, se also Eq. (34). Such asymmetric response at large Rashba-SOI in BLG contrast radically with the linear behavior with field that occurs in single layer graphene [37] at the same regime.

Summary
We have studied the problem of the influence of the Rashba spin-orbit interaction on the bandstructure of biased and unbiased bilayer graphene. Using low energy effective theory we have derived a low energy Hamiltonian for bilayer graphene in the presence of an external magnetic field and spin-orbit interactions. Analytical formulae for the energy spectrum of a graphene bilayer with Rashba spin-orbit interaction are obtained. We show that for relatively weak Rashba coupling the spin-degeneracy of the electron and hole bands is broken inducing a k-linear spin-splitting, very similar to that found in semiconductors heterostructures. At the intermediate strengths of the Rashba effect, the innermost bands interpolates from a k-linear behavior at small momentum, to a k-cubic dependence at high momentum. In contrast, the outermost bands seems to be well described by the cubic spectrum for all values of k. For large values of the Rashba coupling there is a remarkable warping behavior of the spectrum near the Dirac point. Such behavior is unique in biased bilayer graphene. It is found that the bias-induced gap in bilayer graphene decreases as the Rashba is increased, showing a behavior resembling a topological insulator transition phenomena. These peculiar characteristics of the spectrum of bilayer graphene with Rashba spin-orbit interaction may have important consequences on its electronic and spin-transport properties.
We also obtained an analytical expression for the Landau levels and spinpolarization in biased bilayer graphene with Rashba effect valid in the low bias regime. It is further predicted the appearance of an unexpected asymmetric spin-splitting and crossings of the Landau levels due the combined effect between the Rashba interaction and the bias voltage. These results suggests significant consequences on the Shubnikovde Hass oscillations and magnetotransport in bilayer graphene as quantum and spin Hall effects, under the presence of sizable Rashba spin-orbit interaction in the range of few meVs.

Appendix A. Derivation of the low energy Hamiltonian
In this appendix we derive the low-energy Hamiltonian for bilayer graphene in the presence of both, intrinsic and Rashba spin-orbit interaction, as well as Zeeman effect. The low energy Hamiltonian is obtained via Löwdin partitioning, also known as van Vleck's perturbation theory in the context of atomic physics. [44,45,46] First, it is convenient to express the Hamiltonian H k of Eq.(5) in the new spin-dependent atomic basis which now can be written as the sum with δ + = 0 , and δ − = 1. On the other hand Notice that H s is nothing but the single-layer graphene Hamiltonian with Rashba-SOI in the basis {ψ A 2(1)↑ , ψ A 2(1)↓ , ψ B 1(2)↑ , ψ B 1(2)↓ }.

Appendix B. Eigenvalues of the low energy Hamiltonian
By squaring the low energy Hamiltonian (13), a straightforward diagonalizable system is obtained at zero magnetic field (B = 0).