Electronic properties of graphene and graphene nanoribbons with"pseudo-Rashba"spin-orbit coupling

We discuss the electronic properties of graphene and graphene nanoribbons including"pseudo-Rashba"spin-orbit coupling. After summarizing the bulk properties, we first analyze the scattering behavior close to an infinite mass and zigzag boundary. For low energies, we observe strong deviations from the usual spin-conserving behavior at high energies such as reflection acting as spin polarizer or switch. This results in a spin polarization along the direction of the boundary due to the appearance of evanescent modes in the case of non-equilibrium or when there is no coherence between the two one-particle branches. We then discuss the spin and density distribution of graphene nanoribbons.


Introduction
Graphene, the single-layer allotrope of carbon, is undoubtedly one of the most active fields in today's both experimental and theoretical condensed matter physics [1,2,3]. Among an entire plethora of phenomena and proposals, the issue of spin-orbit coupling has generated particular interest [4,5,6,7,8,9]. A detailed understanding of spin-orbit interaction in graphene is crucial for the interpretation of ongoing experiments on spin transport performed by various groups [10,11,12,13,14,15,16,17,18]. Other issues include various device proposals [19,20] and theoretical predictions [21,22,23] related to spins and spin-orbit coupling in graphene.
In the present paper we investigate a single layer of graphene in the presence of spinorbit interaction of the "pseudo-Rashba" type, coupling the sublattice or pseudo spin to the physical electron spin [4,5,6,7,8,9,24]. Our interest is based on the fact that for graphene on Ni with intercalation of Au a 100-fold enhancement of the "pseudo-Rashba" spin-orbit coupling has been reported [25]. Furthermore, impurities which induce a sp3distortion will lead to a "pseudo-Rashba" spin-orbit coupling with a value comparable to the one found in diamond and other zinc-blende semiconductors [26]. The latter result indicates that the "pseudo-Rashba" spin-orbit coupling can be controlled via the impurity coverage.
In this paper, we will concentrate on the scattering behavior of spin densities near boundaries created either by an infinite mass or a zigzag edge. Our presentation is organized as follows: In section 2 we introduce the basic Hamiltonian and discuss its general bulk solution in the absence of a mass term; the technically more complicated case of a nonzero mass is deferred to the appendices. In the following section 3 we investigate in detail the scattering properties and spin dephasing at hard boundaries for various types of incoming spinors and energy ranges. This discussion is extended in section 4 to averaged spin polarizations obtained from continuous distributions of incoming directions. In section 5, we analyze the spin and density distribution of graphene nanoribbons. We close with a summary in section 6. Throughout this manuscript, we use parameters of Ref. [25].
2 Dirac fermions with "pseudo-Rashba" spin-orbit coupling The single-particle Hamiltonian of monolayer graphene with "pseudo-Rashba" spin-orbit interaction can be formulated as [4,5,6,24] where, among standard notation, λ is the spin-orbit coupling parameter, and the Pauli matrices τ , σ describe the sublattice and the electron spin degree of freedom, respectively. For a given wave vector k this Hamiltonian reads explicitly: From experience with the "classic" Dirac equation of relativistic quantum mechanics, it is occasionally of use not to study just a given Hamiltonian but also its square. Here we find where the positive sign corresponds to the eigenvectors while for the negative sign we have where ϑ ∈ [0, π] and cos ϑ = |λ| In the basis (|α 1 , |β 1 , |α 2 , |β 2 ) the Hamiltonian reads and Now it is straightforward to obtain the full eigensystem: We find a gaped pair of eigenvalues with eigenspinors (type I) and With g V = 2 being the valley degeneracy, the corresponding density of states reads The other pair of dispersion branches does not exhibit a gap, and has eigenspinors (type II) The corresponding density of states reads Let us now consider expectation values within the eigenstates with wave functions 1 −k 2 1 y −k y k y ϕ ϕ σ Figure 1: A plane wave of type I with spin perpendicular to the momentum k = (k x , k 1 y ) (ϕ = arctan(k 1 y /k x ), ϕ σ = ϕ + π/2) is reflected at the boundary into a plane wave with k ′ = (k x , −k 1 y ) and k ′′ = (k x , −k 2 y ) with perpendicular spin, but anti-parallel with respect to each other (see Eq. (24) for the definition of k 1/2 y ).
µ ∈ {1, 2}, and A being the area of the system. Here we find and Here, ϕ is the usual azimuthal angle of the wave vector, k = k(cos ϕ, sin ϕ). Note that as usual for Rashba spin-orbit coupling, and where for sin ϑ < 1 sublattice and electron spin degree of freedom are entangled which each other.

Spin dephasing due to reflection on a hard wall
In this section, we will study the scattering behavior from a hard wall which will lead to spin dephasing as depicted in Fig. 1. For that, a general plane wave with fixed momentum k x and energy E ≥ 2|λ| is written as µ ∈ {1, 2} and the normalization constant N k . For energies E < 2|λ|, some modification in Eq. (23) have to be made which shall be discussed in more detail below. In the following, we will discuss the reflection at a hard wall at y = 0 for the two types of plane waves, i.e., we will first set A 1 = 1, A 2 = 0 (type I) and then A 1 = 0, A 2 = 1 (type II). The discussion is based on the reflected spin direction which shall be denoted by ϕ ′ σ . It is obtained from the expectation value of the spin-density operator at the boundary ρ = σδ(ˆ r), ρ ≡ ψ E,kx | ρ|ψ E,kx via Due to translational invariance in x-direction, ρ will only depend on the y-coordinate. For the following discussion, we will also discuss the at r = 0 normalized expectation value σ = ρ / n with n ≡ ψ E,kx |δ(ˆ r)|ψ E,kx . This shall not be confused with the bulk expectation of σ as it appears in the Hamiltonian. We will distinguish the two different cases of the half-plane y ≥ 0 (scattering from the lower or bottom boundary) and y ≤ 0 (scattering from the upper or top boundary). We shall further assume a plane wave with k x > 0 moving in positive x-direction. The results for k x < 0 are then obtained by changing bottom to top boundary and vice versa. The results for the K ′ -point can also be deduced from the following discussion (see appendix A). The sign of λ determines the sign of the expectation value of τ and σ. In the following, we set λ = |λ|, but in some of the following expression we explicitly use |λ| for sake of clarity.
We will discuss two different types of confinement. First, we use the fact that Dirac fermions can be confined by an infinite mass boundary, first discussed by Berry and Mondragon [27]. We then also study the reflection from a zigzag boundary first addressed in Ref. [28].

Infinite mass boundary
With ψ E,kx = (ψ 1 , ψ 2 , ψ 3 , ψ 4 ) T , the boundary conditions at the infinite mass boundary read (see appendix B and C) Note that there are different boundary conditions depending on whether one approaches the boundary from below or above.
Let us first discuss the scattering behavior from the lower boundary. For k x = k cos ϕ, the incident spin direction is given by ϕ σ = π/2 − |ϕ|. On the left hand side of Fig. 2  At large energies with ǫ = λ/(hv F k) ≪ 1 and ǫ ≪ sin 2 ϕ, we have R 1 = (1 − ǫ) cos ϕ and R 2 = i sin ϕ − 2ǫ cos ϕ and the spin polarization is approximately conserved. The expansion of Eq. (25) yields For energies close to the band gap energy of the type I-spinors, E → 2λ, scattering from the boundary acts as a spin polarizer since This is a surprising result since R 1 → −1 and incoming and reflected wave seem to compensate. But even though R 2 → − √ 6ǫe −iϕ sin ϕ tends to zero, its admixture has a dominating effect.
For the upper boundary, we obtain the expansion Note that the different sign compared to Eq. (30) results in a different asymptotic behavior for large energies since ϕ ′ σ (E = 2λ) is larger than the maximal incident spin direction ϕ σ = π. This different behavior is illustrated on the right hand side of Fig. 2.
The sign is determined to yield an exponential decay in the reflected region. In Eqs. (32) and (33), where the upper (lower) sign holds for reflections from the upper (lower) boundary, and s 1 by Let us first discuss the scattering behavior from the lower boundary. On the left hand side of Fig. 3, the reflected spin direction is plotted against the incident spin direction rotated by π. For large energies and normal incident ϕ ≈ π/2, we again obtain ϕ ′ σ = ϕ σ . But for nearly parallel incident such that (E − 2|λ|)/(E + 2|λ|) < (cos ϕ) 2 , we obtain ϕ ′ σ = ±π/2. For energies close to the band-gap E → 2λ, all reflected modes of type I are evanescent and scattering from the wall acts as a switch which leads to either ϕ ′ σ = π/2 or ϕ ′ σ = −π/2. Let us understand the appearance of the two extreme values of ϕ ′ σ = ±π/2 in the regime where k 1 y is imaginary. Since z 1 is real and the incident and reflected wave of type |χ 2,+ compensate, the expectation value in x-direction σ x = 0. For the incident wave, |σ y | incident is negative and for small incident angle, we thus have ϕ ′ σ = −π/2. But if |R 1 | is large, the admixture of |χ 1,+ can lead to ϕ ′ σ = π/2. Additionally, the spin in z-direction σ z assumes a non-zero value to guarantee | σ | = 1. On the left hand side of Fig. 4, this general behavior is shown whether the reflected spin angle (rotated by π), the expectation values σ i (i = x, y, z) and the absolute value of the reflection amplitudes |R 1 | and |R 2 | is plotted versus the incident spin direction at y = 0 at energy E = 3λ.
The scattering behavior from the upper boundary is considerably simpler. There, only two regimes appear with are marked by whether k 1 y is real or imaginary. This can be seen on the right hand side of Figures 3 and 4.

Scattering behavior for plane waves of type II with E < 2λ
For energies with E < 2λ, one of the reflected modes becomes evanescent which leads to σ x = 0. For a more detailed analysis, we have to distinguish the two cases E > λ and   Figure 4: The reflected spin angle (rotated by π), the expectation values σ i (i = x, y, z) and the absolute value of the reflection amplitudes |R 1 | and |R 2 | versus the incident spin direction at y = 0 and A 2 = 0 for energies E = 3λ. We usehv F = 5.6eVÅ and λ = 6meV. Left: Reflection from the lower boundary. Right: Reflection from the upper boundary.
For λ < E < 2λ, the reflected momentum k 1 y = ±iq is imaginary with the same expression as in Eq. (34). The sign is determined to yield an exponential decay in the reflected region. With the ansatz we obtain the same expressions for R 1 → R 1 and R 2 as in Eqs. (32) and (33) with the replacement c 1 → (1 + cos ϑ 1 )/2, s 1 → i (cos ϑ 1 − 1)/2, and z 1 → −i(k x ∓ q)/ q 2 − k 2 x , where the upper (lower) sign holds for reflections from the upper (lower) boundary.
For energies with 0 < E < λ, there is no reflected wave of type I, |χ 1,+ , but one of the reflected momenta of |χ 2,+ is imaginary, k 2 y = ±iq with the same definition as in Eq. (34). With we have  Figure 5: The reflected versus the incident spin direction at y = 0 with A 2 = 0 (left hand side) and A 1 = 0 (rotated by π) (right hand side) for various energies E ≥ 2λ in the case of a zigzag boundary. We usehv F = 5.6eVÅ and λ = 6meV.
We obtain σ y = −1 for the upper and σ y = 1 for the lower boundary, respectively which is independent of the incident direction nor of the energy. Here, we assumed that the bottom boundary is terminated by sublattice A and the top boundary by sublattice B.

Zigzag boundary
For a general plane wave Eq. (23) scattered at y = 0 with energy E ≥ 2λ, the boundary conditions for the bottom boundary (sublattice A) Eq. (39) yield the following expressions for R 1 , R 2 : The boundary conditions for the upper boundary (sublattice B) yield the following expressions for R 1 , R 2 :  The abbreviations are the same as for the infinite mass boundary. Since the reflected angle is symmetric around normal incident, we will only discuss the reflection from the bottom boundary for k x > 0. In Fig. 5, the reflected versus the incident spin direction at y = 0 is shown for the two types of incident plane waves. As in the case of the infinite mass boundary, σ x = 0 for incident plane waves of type II with cos 2 ϕ > (E − 2|λ|)/(E + 2|λ|). But in contrary to the infinity mass boundary, the spin-polarization in out-of-plane direction assumes a non-zero value even when the reflected wave of type I is extended. For this case, i.e., k 1 y ∈ R, we obtain The K ′ -point yields the opposite sign such that there is no net-polarization in z-direction.
For energies E < 2λ a similar discussion as in the case of infinite mass boundary applies.

Spin polarization close to the boundary
So far we have only discussed the polarization properties at the boundary y = 0. For finite y, we expect an oscillatory behavior of the reflected spin polarization. For E → 2λ and plane wave scattering of type I, k 1 y → 0 and the period will thus be solely determined by k 2 y → √ 2(2λ/hv F ). This oscillatory behavior is again independent of the incident spin polarization and results in a striped phase for the reflected spin-polarization. For E > 2λ, two periods related to k 1/2 y contribute and a more complicated pattern emerges which also depends on the incident spin polarization and whether one deals with a reflection from the top or from the bottom. This hints to the fact that a Dirac particle in a box shows quasi-chaotic behavior [29].
In the following, we will study the spin polarization averaged over the incident direction for fixed A 1 , A 2 and including the two K-points as function of the y-direction. We will further average over positive and negative k x -momenta. With an incident wave of type µ We only discuss the spin polarization at the lower boundary which depends on the sign of λ (here we choose λ = |λ|). The spin polarization on the upper boundary is obtained by reversing the sign. In Fig. 6, the angle-averaged spin density A ρ x µ ( r) is shown as function of y for various energies E ≥ 2λ where A denotes the area of the sample. We show the results for an incident plane wave of type I (left hand side) and type II (right hand side) with an infinite mass boundary. There is a clear difference between the two types for low energies which is due to the appearance of imaginary momenta k 1 y = ±iq for type II-reflections. For low energies, most incident angles of the initial plane wave of type II lead to evanescent modes and thus to σ x = 0. For large energies E ≥ 10 3 λ, the spin polarization of the two types have approximately the same absolute value, but differ in sign.
Obviously, the above ensemble average breaks time-reversal symmetry since there is one incident plane wave with fixed k y -direction and two reflected plane waves. But if there is no coherence between the incident plane waves of type I and II, e.g., due to temperature, then time-reversal symmetry is effectively broken and we find a net polarization in xdirection by adding the two contributions ρ I and ρ II (and possibly weighting them with the corresponding density of states). This is demonstrated in Fig. 7, where the sum of the two contributions A µ ρ µ is shown for a infinite mass boundary (left) and for a zigzag boundary (right). Moreover, we expect spin polarization in x-direction for various non-equilibrium situations.
In the other two directions, we find no net spin polarization if the two inequivalent K-points are included. We note, however, that ρ y and ρ z assume a finite value for one K-point, only. This opens up the possibility of spin polarization in these directions in the presence of ripples or a magnetic field. Especially surface states due to, e.g., zigzag boundaries which effectively break the sublattice symmetry and which are not included in our continuous model should give rise to a finite spin polarization.

Dirac electrons with "pseudo-Rashba" spin-orbit coupling in nanoribbons
In this section, we will consider graphene nanoribbons and the quantization properties of the transverse momenta in the presence of "pseudo-Rashba" spin-orbit coupling. We will then discuss the density and spin distribution at various energies.

Quantization of the transverse momentum
Let us first consider infinite mass boundaries. For a general plane wave with fixed momentum k x and energy E, ψ E,kx (x, y) ≡ (ψ 1 , ψ 2 , ψ 3 , ψ 4 ) T , there are four conditions that have to be satisfied, i.e., ψ 1 = ±ψ 3 and ψ 2 = ±ψ 4 at y = 0, and ψ 1 = ∓ψ 3 and ψ 2 = ∓ψ 4 at y = W , where the upper (lower) sign stands for the K(K ′ )-point and W the width of the nanoribbon. For a zigzag nanoribbon which terminates on sublattice A at the bottom and on sublattice B at the top, the four conditions read ψ 1 = ψ 2 = 0 at y = 0 and ψ 3 = ψ 4 = 0 at y = W .
Let us first assume two propagating waves as in Eq. (23), see also Fig. 8. In order to have a non-trivial solution, a necessary condition is with the bar denoting the complex conjugate. For infinite mass boundaries, the above matrices read at the K-point and for zigzag boundaries, we have where we introduced w µ = e ik µ y W and used the definitions of section 3. det M in Eq. (46) is real and thus yields the quantization of the transverse momentum in y-direction.
For (hv F k 2 y ) 2 < 4Eλ, there is the appearance of evanescent modes since k 1 y = ±iq is imaginary. In this case, a general plane wave with fixed momentum k x and energy E ≥ 2|λ|, ψ E,kx (x, y) ≡ (ψ 1 , ψ 2 , ψ 3 , ψ 4 ) T , is written as ψ E,kx (x, y) = N k e ikxx A 1 e −q(W −y) |χ 1,+ (k x , −iq) + A 2 e ik 2 y y |χ 2,+ (k x , k 2 y ) + R 1 e −qy |χ 1,+ (k x , iq) + R 2 e −ik 2 y y |χ 2, with Again, in order to have a non-trivial solution, Eq. (46) must hold, but this time the matrices for infinite mass boundaries at the K-point read and for zigzag boundaries, we have with w 1 = e −qW , z ± 1 = (k x ±q)/ k 2 x − q 2 , c 1 → (1 + cos ϑ 1 )/2 and s 1 → i (cos ϑ 1 − 1)/2. The definitions for the plane wave of type II remain unchanged. Since the wave function of the evanecent mode is now real, the matricesĀ,B are not the complex conjugates of A, B, but given bȳ for infinite mass boundaries, and for zigzag boundaries they read It is now preferable to write Eq. (46) in powers of w 1 . For zigzag boundaries, this yields which is purely imaginary and thus again yields a quantization of the transverse momentum in y-direction. For infinite mass boundaries, we obtain a similar expression.

Spin and density distribution
The particle density n at energy E is now obtained by summing over all transverse modes n that obey the above boundary conditions or the corresponding boundary conditions for the K ′ -point. Denoting the n-th transverse momentum of type II by k 2 y,n , we have In Fig. 9, the density distribution of a graphene nanoribbon of width W = 100nm for various low energies with infinite mass (left) and zigzag (right) boundaries is shown. In general, the number of modes is the same with and without "pseudo-Rashba" spinorbit coupling and the resulting density distributions only differ slightly. But for zigzag boundaries at E = 4λ, we observe strong deviations due to the fact that there are 8 modes in the case with spin-orbit coupling in contrast to 12 modes in the case without spin-orbit coupling. Also note that whereas for the case without spin-orbit coupling, all modes are extended, some modes for the case with spin-orbit coupling are evanescent for the type I-branch. For zigzag boundaries, e.g., we have no extended and 4 evanescent (type-I) modes at E = 2λ, 6 extended and 2 evanescent modes at E = 4λ and 6 extended and 6 evanescent modes at E = 6λ. The spin polarization at the boundaries is in all cases zero; in x-direction it is zero also for only one K-point, in y-and z-direction it is non-zero for one K-point, but averages to zero when two K-points are included. This is an immediate consequence of time-reversal symmetry. In Ref. [30], spin polarization in z-direction is reported for various k x -values within a lattice model of a zigzag nanoribbon. At equilibrium, this can only be attributed to edge-states which effectively break the sublattice symmetry and which are not included in our continuous model.

Summary
In this paper, we have investigated the spin dephasing of Dirac fermions with "pseudo-Rashba" spin-orbit coupling due to the reflection from a hard wall. In order to confine the Dirac electrons, we used infinite mass and zigzag boundaries. For large energies compared to the spin-orbit coupling, we obtained the expected result that there is hardly spindephasing due to the scattering process. But for energies close to the band gap for plane waves of type I, E ≈ 2λ, strong spin dephasing is observed. If the incident plane wave is of type II (gapless branch), even stronger effects are seen like the appearance of evanescent modes. We also observe the rotation of the spin in out-of-plane direction away from the boundary and for incident plane waves of type II also at the boundary. This effect will be canceled by averaging over the two inequivalent K-points.
We also discussed the spin polarization averaged over the incident direction and including the two K-points. We find that for energies E ≥ 2λ, there is a finite spin polarization in x-direction when there is no coherence between the two branches. This polarization differs in sign for the upper and lower boundary, respectively. Also for non-equilibrium situations, there will be a spin polarization in this direction.
We finally analyzed the spin and density distribution of graphene nanoribbons. At certain energies, the number of transverse modes does not match the one of a corresponding nanoribbon without "pseudo-Rashba" spin-orbit coupling. This results in significant changes in the density distribution. But generally, the "pseudo-Rashba" spin-orbit coupling leads to marginal differences, only. Further, there is no spin polarization if both K-points are included, but we find a finite spin polarization in y-and z-direction for one K-point, only. Surface states due to, e.g., zigzag boundaries which only live on one sublattice and thus break the valley-symmetry should therefore yield a finite spin polarization.

Acknowledgments
T.S. wants to thank Nuno Peres and João Lopes dos Santos for illuminating discussions and support. We further thank E. Rashba

Appendix A: The full model including the two Kpoints
The full model including the two K-points reads where κ z = ±1 denotes the two inequivalent K-points. For a given wave vector k the Hamiltonian around the K ′ -point (κ z = −1 ) reads The Hamiltonian around the K ′ -point can thus be obtained from the Hamiltonian around the K-point by interchanging the pseudo-spin index and reversing the sign. All previous results without the mass term can thus be used. The results involving the mass term are obtained by M → −M . This leads to a change in the boundary conditions, i.e., 9 Appendix B: Massive Dirac fermions with "pseudo-Rashba" spin-orbit coupling Massive Dirac fermions with "pseudo-Rashba" spin-orbit interaction can be described by where, among standard notation, λ is the spin-orbit coupling parameter, and the Pauli matrices τ , σ describe the sublattice and the electron spin degree of freedom, respectively.
Here we have assumed a positive spin-orbit coupling parameter, λ = |λ|, and ϕ is the usual azimuthal angle of the wave vector, k = k(cos ϕ, sin ϕ). Note that massive Dirac fermions assume a non-zero expectation value for the pseudo-spin and spin in z-direction.
10 Appendix C: Scattering from infinite mass boundary Dirac fermions can be confined by an infinite mass boundary, first discussed by Berry and Mondragon [27]. In the following, we will study the scattering behavior from a boundary located at y = 0 and y = W . Within the strip 0 < y < W , the mass of the Dirac fermions shall be zero; outside the strip, the mass shall be infinite. A general plane wave within the strip with fixed momentum k x and energy E > 0 can be written as ψ E,kx (x, y) = e ikxx A 1 e ik 1 y y |χ 1,+ (k x , k 1 y ) + A 2 e ik 2 y y |χ 2,+ (k x , k 2 y ) + R 1 e −ik 1 y y |χ 1,+ (k x , −k 1 y ) + R 2 e −ik 2 y y |χ 2,+ (k x , −k 2 y ) , withh µ ∈ {1, 2}. The wave function of the transmitted electron is also decomposed by the two eigenfunctions |χ µ,+ , ψ E,kx (x, y) = e ikxx T 1 e ik 1 y y |χ 1,+ (k x , k 1 y ) + T 2 e ik 2 y y |χ 2,+ (k x , k 2 y ) , withh In the limit M → ∞, the transmitted plane wave simplifies and ψ E,kx (x, W ) = e ikxx T 1 with s λ = λ/|λ|. The different expressions at y = 0 and y = W originate from the different sign ofhv F k y → ±iM v 2 that has to be chosen to yield an exponential decay in the infinite mass region. It therefor only depends on whether one deals with the upper or lower boundary. At the boundaries y = 0 and y = W , the four components have to be continuous to guarantee a continuous current which leads to the following two sets of equations: With ψ E,kx = (ψ 1 , ψ 2 , ψ 3 , ψ 4 ) T , this translates to the familiar boundary condition from Ref. [27] for the two spin channels, respectively: