Surface orbitronics: new twists from orbital Rashba physics

When the inversion symmetry is broken at a surface, spin-orbit interaction gives rise to spin-dependent energy shifts - a phenomenon which is known as the spin Rashba effect. Recently, it has been recognized that an orbital counterpart of the spin Rashba effect - the orbital Rashba effect - can be realized at surfaces even without spin- orbit coupling. Here, we propose a mechanism for the orbital Rashba effect based on sp orbital hybridization, which ultimately leads to the electric polarization of surface states. As a proof of principle, we show from first principles that this effect leads to chiral orbital textures in $\mathbf{k}$-space of the BiAg$_2$ monolayer. In predicting the magnitude of the orbital moment arising from the orbital Rashba effect, we demonstrate the crucial role that the Berry phase theory plays for the magnitude and variation of the orbital textures. As a result, we predict a pronounced manifestation of various orbital effects at surfaces, and proclaim the orbital Rashba effect to be a key platform for surface orbitronics.

T He spin Rashba effect as a fundamental manifestation of the spin-orbit coupling (SOC) at surfaces has revolutionized the spintronics research, and served as a foundation for a new field of spintronics rooted in relativistic effects − spinorbitronics [1][2][3][4]. As a consequence of the spin Rashba effect, the surface states created as a result of SOC and surface potential gradient exhibit finite spin polarization, which forms chiral textures in reciprocal k-space. This gives rise to a multitude of prominent phenomena such as Dzyaloshinskii-Moriya interaction [5][6][7][8][9], spin Hall effect [10], direct and inverse Edelstein effects [11,12], quantum anomalous Hall effect [13], and current-induced spin-orbit torques [14][15][16].
Only recently it has been realized that many of the spinorbitronic effects can be rethought in terms of their "parent" orbital analogs brought to the "spin" level by SOC [17][18][19][20][21][22][23]. This leads to the notion of orbitronics as a promising branch of electronics which operates with the orbital degree of electrons rather than their spin. The orbital moment (OM) of electrons in a solid, being an evasive quantity until very recently [24,25], offers higher flexibility with respect to its magnitude and internal structure as compared to spin, and it is a key quantity for novel and promising effects in orbitronics such as gyrotropic magnetic effect [26] and orbital Edelstein effect [27]. The orbitronic analogue of the spin Rashba effect is the orbital Rashba effect (ORE), which results in the emergence of chiral OM textures of the states in k-space even without SOC, and manifests in the "conventional" Rashba effect in the spin channel upon including SOC into the picture [17][18][19][20][21]. Indeed, circular dichroism measurements of Au(111) [20] and Bi 2 Se 3 [21] surfaces are consistent with the prediction of the ORE forming non-trivial chiral OM textures in k-space. On the other hand, quasi-particle interference experiments performed on another Rashba systems, the surface of Pb/Ag(111) [28] and Bi/Ag(111) [29], demonstrate clear fingerprints of spin-orbital-flip scattering, which implies the * godw2718@postech.ac.kr vital role of the orbital degree of freedom of electrons for surface scattering. As the spin Rashba effect played a key role in the development of spintronics, our understanding of the ORE and the ability to control its properties is pivotal for bringing orbitronics to the surface realm.
In this work, we introduce a mechanism for the ORE at surfaces, which originates in the orbital sp hybridization and gives rise to pronounced ORE even without SOC in sp surface alloys such as BiAg 2 monolayer. We elucidate that the sp hybridization gives rise to the orbital magnetoelectric coupling and intrinsically links the ORE with electric polarization of states at the surface. Moreover, by referring to explicit first principles calculations, we demonstrate that including nonlocal effects [30][31][32][33][34][35][36] can drastically enhance the magnitude of the ORE-driven orbital polarization of the surface states as compared to the k-dependent OM obtained from simple atomic arguments. We speculate that such an enhancement, which is particularly prominent in the vicinity of band crossings in the surface electronic structure, can result in a giant magnitude of various orbitronics phenomena such as orbital Hall effect [22,23], gyrotropic magnetic effect [26], and current-induced orbital magnetization [27]. This marks the ORE as the most promising platform for realizing orbitronics at surfaces.

Results
Tight-binding model. We start the discussion of the ORE with tight-binding considerations. While previously only pure pand d-orbital systems have been studied [18,37], here we study monolayer of sp-derived surface alloys such as BiAg 2 [28,38,39], which is different in that sp hybridization is important. The tight-binding Hamiltonian can be generally written as:  . The combination of a strong surface potential gradient and a pronounced hybridization between Bi 6pand Ag 5s-states gives rise to the orbital Rashba effect. (b) The hexagonal Brillouin zone and high-symmetry points. The tight-binding (c) and the first principles (d) band structures were calculated with and without taking into account spin-orbit coupling (SOC). States of distinct orbital symmetry are denoted as p−1, p0, and p+1 following the convention of the main text. The tight-binding parameters were chosen so as to closely reproduce the first principles bandstructure.
three p orbitals and one s orbital at each site, we choose the basis states as |ϕ nk = (1/ √ N ) R e ik·R |φ nR . Here, φ nR denotes the n-th (n = p x , p y , p z , s) atomic orbital centered at the Bravais lattice R, and N is the number of lattice sites. In this basis, the hybridization assumes the form: where γ sp is the nearest-neighbor hopping amplitude between p x(y) and s orbitals, a is the lattice constant, V z (k) = eE z ϕ pzk | z |ϕ sk describes the effect of sp z hybridization by a surface potential gradient E z , and e > 0 is the elementary positive charge. The effect of SOC is included into the Hamiltonian as H soc = λ socL ·Ŝ, where λ soc is the spin-orbit strength,Ŝ is the operator of spin, and L = RL R , is the representation of the atomic contribution of the orbital angular momentum operator [41]. Here, r and p denote the canonical position and momentum operators, respectively. Although the general conclusions that we draw from the tightbinding model do not depend on the exact choice of the hopping parameters, they can be easily tuned such that the tightbinding band structure closely resembles the first principles one of the BiAg 2 monolayer (Fig. 1).
Orbital Rashba effect. Neglecting for the moment the effect of SOC and assuming |H p (k) − H s (k)| |h(k)| as is the case for BiAg 2 in the long-wavelength limit of k → 0, we can perturbatively downfold the h(k) term to arrive at an effective Hamiltonian (see Supplementary Information): where H p,eff (k) = H p (k) + H OR (k). The expression is known as the orbital Rashba Hamiltonian since it resembles the conventional spin Rashba Hamiltonian with the orbital angular momentum operator replacing that of spin. Remarkably, the combined effect of sp hybridization and surface potential gradient is concisely described by the orbital Rashba Hamiltonian within the subspace spanned by p orbitals. Thus, the orbital Rashba physics arises not from the electron's spin but from the orbital degrees of freedom even without SOC. In analogy to the Rashba constant of the conventional spin Rashba model, the parameter is called the orbital Rashba constant, with ∆E sp (k) denoting the energy gap between sand p-derived bands, and η ∼ 1 being a parameter dependent on the lattice structure. Using the specific tight-binding model parameters for BiAg 2 (see Supplementary Information), we estimate the orbital Rashba constant to be about 1 eV·Å. From Eq. (6) it is clear that the ORE roots in the sp orbital hybridization, and that the strength of the orbital Rashba effect is directly proportional to the value of the surface potential gradient associated with the buckling of Bi atoms in BiAg 2 . This implies that the desired properties of the ORE can be designed by controlling the band hybridization via chemical and structural engineering.
Crystal field splitting. In contrast to the conventional Rashba effect, the ORE is very sensitive to the crystal field splitting (CFS), which quenches the OM. In the absence of the CFS, that is, when , E pz as corresponding energy eigenvalues of p x(y) -and p z -derived states at k = 0, the expectation value of the OM is given by within the so-called atom-centered approximation for the OM, which takes into account only intra-atomic contributions.
Here,ẑ andk are unit vectors along the z-axis and and vector k, respectively, ψ p l k is the p l -derived eigenstate of H p,eff (k), and l = {−1, 0, +1} is the angular momentum quantum number with respect to the quantization axisẑ ×k. In the vicinity of the Fermi energy, the p-derived bands can thus be denoted in terms of their dominant orbital character as p −1 , p 0 , and p +1 in the order of increasing energy ( Fig. 1c-d). The corresponding k-dependent eigenenergies are E −1 (k), E 0 (k), and E +1 (k). The associated k-linear orbital-dependent energy splitting arising due to the orbital Rashba term in H p,eff (k) amounts to ∆E OR = lα OR |k| in the long-wavelength limit. However, in the presence of the CFS, although its direction remains intact, the expectation value of the OM is reduced by a factor of 2α OR |k|/∆ CFS , when assuming |∆ CFS | |α OR k|. For this reason, orbital-dependent energy splitting appears from the second order in k if ∆ CFS = 0, which can be seen from Eq. (5).
Relation to electric polarization. Consider an orbitalcoherent state ψ p l k , that is, an eigenstate of H hop (k) with ∆ CFS = 0 which exhibits a quantized value of OM (e.g., due to an applied in-plane magnetic field or the ORE): where l = {−1, 0, +1} is the orbital angular momentum quantum number with respect to the direction of some inplane unit vectorn k . While in previous studies the electric polarization within the ORE was introduced phenomenologically [17,19], one can show from perturbation theory arguments that such an orbital-coherent state naturally exhibits electric polarization perpendicular to the surface plane: where the parameter is what we call the orbital Rashba susceptibility. The structure of χ OR reflects that the electric polarization arises from the sp hybridization (Fig. 2a). From Eq. (6) it is clear that α OR = E z χ OR , and that the orbital Rashba Hamiltonian as given by Eq. (5) can be understood as a dipole coupling between an electric field at the surface with the operator of electric polarization, 2b). Physically, it means that within the ORE the nontrivial orbital texture of states in k-space arises so as to gain maximal energy by the interaction of the states' polarization with the surface electric field. Thus, the ORE can be seen as a consequence of orbital magnetoelectric coupling at surfaces.
Atom-centered approximation from first principles. In order to verify our predictions in a realistic system, we evaluate the kand band-resolved value of the OM in BiAg 2 from first principles, which is given by m ACA n (k) = −(µ B /h) µ ψ nk | r µ × p |ψ nk µ in ACA. Here, ψ nk is an eigenstate, r µ is the position operator relative to atom µ, the summation is performed over all atoms in the lattice, and the real-space integration is restricted locally to atom-centered muffin-tin spheres. Computed in such a way OM is a direct generalization to the first principles framework of the OM computed in ACA from tight-binding. In Fig. 3a-c, we show the distribution of the OM in k-space for the p-derived bands p −1 , p 0 , and p +1 in absence of SOC, and observe clockwise (−1), zero, and counter-clockwise (+1) type of chiral behavior of this distribution around the Γ-point, respectively. The magnitude of the OM without SOC along different highsymmetry lines is shown in Fig. 3d. While the OM reaches as much as 0.27µ B for the p +1 band in the middle of the Brillouin zone, we can use its behaviour with k in the vicinity of the Γ-point to estimate the magnitude of the orbital Rashba constant α OR . First, we note that a small but finite value of the slope in the OM of p 0 as well as visible differences in the behavior of the OM of p −1 and p +1 bands can be observed, in contrast to our prediction Eq. (8). The reason for the discrepancy lies in a non-zero crystal field splitting, which we can estimate from Fig. 1d to be ∆ CFS ≈ 0.76 eV. Taking this into account, the orbital Rashba constant of p −1 and p +1 bands amounts to 0.96 eV·Å and 1.38 eV·Å, respectively. Similarly, we also find orbital chiralities near K-point. Investigating further the influence of SOC, we found no qualitative changes in the distribution of the OM in the Brillouin zone. However, we clearly see that the resulting "Rashba" spin texture in reciprocal space, emerging upon including SOC, is strictly bound to the local direction of the OM. Details on these results are provided in the Supplementary Information. Finally, we display in Figs. 4f and 4g the chirality of the OM at the Fermi surface both with and without SOC taken into consideration, respectively.
Berry phase theory in films. The primary manifestation of the ORE in solids is the kand band-dependent generation of finite OM. For understanding this fundamental effect and predicting its magnitude, it is crucial to evaluate the magnitude of the OM properly without assuming any approximations such as the ACA. Recently, it was shown that a rigorous treatment of OM in solids within the complete Berry phase description [30][31][32][33] naturally accounts for non-local effects [35]. Thereby, theoretical estimations of OM in magnetic materials of various nature have significantly improved [24,25,36].
Since the ORE manifests in the generation of in-plane OM, following the procedure of Ref. [33], we extended the previous formulation to the case of in-plane components of OM of a film finite along the z-axis. The expression for the in-plane components of the OM is given by where u nk is an eigenstate of the lattice-periodic Hamiltonian H(k) with the eigenvalue E n (k), and E F is the Fermi energy. Equation (11) can be further decomposed into the so-called local circulation term The latter expression is connected to the projected Berry curvature which is closely related to the well-known Berry curvature in bulk systems by formally replacing z with i∂ kz . Inserting the model orbital Rashba Hamiltonian, Eq. (5), directly into Eq. (11), we find that in the long-wavelength limit, where l is the angularmomentum quantum number defined in Eq. (7). Assuming the typical values E z ∼ 0.1 eV/Å and (E p x(y) +E pz −2E F ) ∼ 1 eV, we find not only that ACA and modern OM exhibit the same chirality of the distribution, but also that m mod l (k) ∼ m ACA l (k) in the vicinity of the Γ-point within the model analysis.
Berry phase theory from first principles. To investigate whether significant differences between the ACA and modern treatment of OM arise in a realistic situation, we evaluate from first principles m ACA n (k) and m mod n (k). In Fig. 4, we show the OM distribution evaluated from the modern theory (Figs. 4a, 4c) and ACA (Figs. 4b, 4d), where the individual contributions of all occupied bands were summed up for each k point. In this figure, the distribution of the OM within the modern approach and ACA is similar around the Γ-point both in magnitude and overall distribution, in accordance with our model considerations. Near the K-point, however, the distribution of the modern OM without SOC deviates significantly from the ACA result. If SOC is taken into account, the difference between the two approaches is even more drastic as it amounts to one order of magnitude. In Fig. 4, the visible discontinuity of the OM distribution occurs along the Fermi surface lines (Figs. 4e-4h). Directly at the Fermi surface, the itinerant contribution to the OM in the modern theory vanishes and we may restrict ourselves to visualizing in Figs. 4e and 4g the local circulation from the modern theory, without and with SOC, respectively. In contrast, within ACA there is no such decomposition of the Fermi-surface contribution (Figs. 4f, 4h). The OM of Fermi-surface states plays a crucial role in various orbital magnetoelectric phenomena [26,27].
Within the Berry phase theory, one of the most remarkable features of the OM in bulk is its correlation with the Berry curvature in k-space, which often exhibits a spiky behavior in the vicinity of band crossings [34]. In our formalism for the inplane OM in thin films as expressed by Eq. (11), the projected Berry curvature, Eq. (12), is a key ingredient. It behaves similarly to the conventional Berry curvature in that it can exhibit singular behavior at band crossings as a consequence of the rapid variation of the wave functions with k. Following this spirit, we also seek for such a spiky behavior in the local circulation m LC (k) by studying the OM in the vicinity of the band crossing in the electronic structure of BiAg 2 which is close to the M-point and about 0.2 eV below the Fermi level. Setting E F in Eq. (11) to the energy of the crossing and treating all bands below this energy as occupied, we observe a singular behavior of the OM magnitude within the modern theory in the vicinity of the crossing point (Fig. 5a). This behavior can be directly correlated with sizable OM contributions of those states which constitute the band crossing. At the point of singularity, the modern-theory OM is purely due to the local circulation and it reaches as much as 1.01 µ B in magnitude, while the ACA predicts a tiny value of 0.03 µ B . Remarkably, both ACA and modern theory agree qualitatively in their prediction of the behavior of the OM away from the band crossing (see Supplementary Information). Two-band model near the band crossing. To study the aforementioned behavior of the OM near the band crossing from the model point-of-view, we consider the following two-band Hamiltonian: where σ is the vector of Pauli matrices reflecting the orbital degree of freedom of two chosen basis states: |ϕ s,k , and |ϕ p,k = (1/ √ 2) |ϕ px,k − (i/ √ 2) |ϕ pz,k . We assume that d x (k) = 0, d y (k) = γ sp k x a + V z , and d z (k) = δ, where a is the lattice constant, δ is the energy gap of the massive Dirac-like dispersion along k x , γ sp is the nearestneighbor hopping amplitude between s and p x orbitals, and V z = −(eE z / √ 2) ϕ s,k | z |ϕ pz,k expresses the breaking of inversion symmetry. The band structure of this Hamiltonian is shown in Fig. 5a, where we assumed a = 5Å, V z = 0.1 eV, γ sp = 1 eV, and δ = 5 meV. We represent the position operator z in the specified basis as where p x = p z = 0 and p y = −(1/ √ 2) ϕ s,k | z |ϕ pz,k . Taking into account the representation (15) and applying Eq. (11), we obtain an expression for the in-plane OM from the modern theory for the general two-band Hamiltonian (14): whered(k) is the direction of d(k) and the "+" ("−") sign stands for upper (lower) band. In the second line, owing to the fact that in our model (ẑ × ∂ k )p(k) ≈ 0, we related the modern OM to the derivative of P ± z (k) = ±(−e)[d(k) · p(k)], that is, the z-component of the electric polarization of the upper (lower) band as defined in Eq. (9).
From this generic expression we clearly observe that the origin of the non-vanishing OM within the modern theory lies in the gradient of the electric polarization in reciprocal space. Using the parameters stated above and setting E F = 3 meV and p y = 2Å, we compute the OM of all occupied states in the vicinity of the band crossing from both modern theory, Eq. (16), and ACA, Eq. (3). By replacing further E F with |d(k)| in Eq. (16), we obtain the Fermi-surface contribution due to the self-rotation of the wave packet [30], which determines the orbital magnetoelectric response [26,27]. Based on the results shown in Fig. 5b, we conclude that while the values of m ACA ± (k) are strictly bound by µ B everywhere as can be confirmed analytically, the pronounced singular behavior of m mod ± (k) with large values within the gap can be attributed solely to the rapid variation of the electric polarization in the vicinity of the crossing.

Discussion
We have shown that sp orbital hybridization is the main mechanism for the ORE at surfaces of sp-alloys, which is manifest already without SOC. Just like the spin Rashba effect follows from the ORE via SOC, we can expect other orbital-dependent phenomena to be formulated and discovered, from which SOC recovers their spin analogues. One remarkable example is orbital analogue of the quantum anomalous Hall effect, where the quantized edge state is orbital-polarized [42]. Another ex-ample is the recent formulation of the orbital version of the Dzyaloshinskii-Moriya interaction governing the formation of chiral structures such as domain walls and skyrmions [43].
Our simulations reveal the complexity of the orbital textures driven by ORE, and their sensitivity to the electronic structure of realistic materials. This means that the desired properties of the ORE, and phenomena it gives rise to, can be designed by proper electronic-structure engineering. Moreover, by making use of both spin and orbital degrees of freedom, one can generate and manipulate arising spin and orbital textures in k-space. While for the case of BiAg 2 considered here the spin aligns collinearly to the OM in the presence of SOC, a Rashba effect with respect to the total angular momentum emerges in each of the j = 3/2 and j = 1/2 branches in the regime where SOC is dominant over the ORE [17,19]. Based on a similar idea, an orbital version of the Chern insulator state in a situation of very large SOC has been recently proposed [44].
We predict that within the ORE, in contrast to the spin Rashba effect, the magnitude of the OM in the vicinity of band crossings can reach gigantic values. This observation has very far-reaching consequences for the magnitude of effects which are directly associated with the OM at the Fermi surface. This particularly concerns the orbital magnetoelectric effect, as recently discussed in the context of orbital Edelstein effect [27] and gyrotropic magnetic effect [26]. Within the orbital Edelstein effect, a finite OM at the Fermi surface is generated by an asymmetric change in the distribution function created by an electric field. In the gyrotropic magnetic effect, discussed intensively these days with respect to topological metals [45,46], an external magnetic field gives rise to an electrical current [26]. Both phenomena rely crucially on the magnitude of the local orbital moments at the Fermi surface of materials, and we predict here that they can be drastically enhanced by tuning the electronic structure such that the singularities in the OM, which we disclose in our work, are positioned at the Fermi level. The generally enhanced predicted magnitude of the Berry-phase OM of the occupied states, as compared to the OM computed from the commonly used approximation, can also manifest in the reevaluation of the magnitude of other effects such as the orbital Hall effect [22,23].
An additional flavor to the ORE is the intrinsic valley degree of freedom inherent to systems of the type studied here: because of time-reversal symmetry the points of singularity in the OM always come in pairs. Since the OM is opposite for opposite valleys, the overall OM integrated over the Brillouin zone is zero in the ground state, and orbital magnetoelectric response becomes strongly valley-dependent. Exploiting the ORE for the purpose of generating sizable ground state net OM at surfaces has to be done in combination with generating a non-vanishing exchange field and magnetization in the system. This can give rise to a plethora of effects relying on intertwined spin and orbital degrees of freedom in complex magnetic materials.

Methods
First principles calculation. We performed self-consistent density-functional theory calculations of the electronic structure of BiAg 2 using the film mode of the FLEUR code [47], which implements the FLAPW method [48,49]. Exchange and correlation effects were treated within the generalized gradient approximation [50]. We assumed a ( √ 3 × √ 3)R30 • unit cell (see Fig. 2a) with the in-plane lattice constant a = 9.47 a 0 , where a 0 is Bohr's radius. The surface relaxation of Bi was set to d = 1.61 a 0 , and the muffin-tin radii of Bi and Ag were chosen as 2.80 a 0 and 2.59 a 0 , respectively. We used K max = 4.0 a −1 0 as plane-wave cutoff and sampled the irreducible Brillouin zone using 110 points. Spin-orbit coupling was included self-consistently within the second-variation scheme [51].
Based on the converged charge density, maximallylocalized Wannier functions (MLWFs) were obtained in a post-processing step employing an equidistant 16 × 16 kmesh. Starting from sp 2 and p z trial orbitals on Bi as well as s, p, and d trial functions on Ag, we constructed 44 MLWFs out of 120 energy bands using the WANNIER90 program [52]. The frozen window was set 2.78 eV above the Fermi energy. Subsequently, we calculated the OM in (i) the ACA, and (ii) the Berry phase theory according to the scheme proposed by Lopez et al. [25].