Tunable spin and orbital polarization in SrTiO₃-based heterostructures

In the search for a novel platform for spintronics devices, oxides heterostructures has emerged as one of the most promising candidates due to their excellent electrical and magnetic properties [1]. The formation of two-dimensional electron gas [2–5] and magnetism [6–10] at the interface of transition metal oxides such as ${\mathrm{LaAlO}}_{3}/{\mathrm{SrTiO}}_{3}$ (LAO/STO) gives rise to a unique opportunity to study and implement spintronics in a single platform [11].

A key requirement to spintronics application is the presence of spin–orbit coupling (SOC). As the inversion symmetry at the interface of the two transition metal oxides is inherently broken, there exists SOC of Rashba type [12–14]. Magnetotransport measurements have detected a substantial spin–orbit splitting of 2–10 meV at the LAO/STO interface [14–17], while a much stronger spin-splitting of 90 meV has been observed on the surface of STO [18]. Strikingly, the spin splitting is not just of k-linear [14] but also k-cubed Rashba types [19]. Furthermore, the interfacial Rashba coupling is versatile, and can be tunable via electrical gate [12–14], and via crystal orientation selection [20–22]. For comparison, the above splitting is higher than that of 1–5 meV in typical semiconductor heterostructures [23]. The strong Rashba coupling in LAO/STO has led to the realization of various spintronics effects such as spin-to-charge conversions [24–28] and spin–orbit torques [29]. The fact that these phenomena are interfacial in nature allows us a much better degree of external control of its electrical and magnetic properties compared to conventional spintronic materials such as ferromagnetic metals [30, 31]. More importantly from the application standpoint, the conducting layer is protected against external perturbations as it is sandwiched between conventional oxides insulators [30, 31], which makes it more robust in comparison to the exposed surface states of topological insulators.

At the same time, theoretical models have been proposed to elucidate the rich SOC properties of the STO-based heterostructures, which are not captured by the standard linear Rashba theory [32]. Combined first-principles calculations and tight-binding approach [13, 33, 34] revealed a multi-orbital Rashba effects, of the linear and cubic types, at the LAO/STO interface layer. Nevertheless, the crystal field, one of the ingredients of the SOC which accounts for the properties of the heterostructure, had been deliberately introduced, and the contribution of the bulk STO as well as the inter-orbital hopping was ignored. Significantly, the large spin splitting at the anti-crossing region, a signature of the multi-orbital effects, has not been fully elucidated. Recently, a theory based on k.p model [35, 36], including the contribution of the bulk STO, predicted a spin splitting induced by the inter-orbital hopping and interfacial band-bending. However, a complete picture of the multi-orbital Rashba SOC is still lacking.

In this work, we aim to construct the effective Hamiltonian of Rashba spin–orbit coupling system in STO-based structures, taking into account the underlying factors such as the confinement effect, as well as interface effects including inversion symmetry breaking (ISB) and band-bending. We derive analytical expressions of properties, e.g. Rashba parameter, effective mass, band edge energy and orbital occupancy, as functions of material and tunable heterostructure parameters. We show that the type and strength of spin–orbit coupling as well as the orbital selectivity are tunable by controlling geometric factor such as STO thickness and by applying electric field via gate voltage. The thickness and gate control of SOC provides a possible avenue for optimization of Rashba SOC for spintronics applications and devices.

We consider a 2DEG system in STO-based heterostructure. The 2DEG can be generated by interfacing STO with either LAO [37] or vacuum with δ-doping [22]. In the low energy limit, the electron in bulk STO is described by following Hamiltonian

$$\begin{eqnarray}&&{H}_{\mathrm{full}}={H}_{{t}_{2g}}+{H}_{\mathrm{ISB}}+{H}_{\mathrm{SO}}+V(z).\end{eqnarray} \tag{ 1 }$$

In the above ${H}_{{t}_{2g}}$ is the kinetic energy of the d-electron in the ${{\rm{\Gamma }}}_{25}^{+}$ band in cubic crystal, ${H}_{\mathrm{ISB}}$ is the inter-orbital hopping, ${H}_{\mathrm{SO}}$ is the atomic spin–orbit interaction, and the last term is the electrical confinement potential in the z-direction $V(z)={{ \mathcal E }}_{z}z$.

In the t_2g orbital-spin basis ${{\rm{\Phi }}}_{{t}_{2g}}=\left\{\left|{yz},\uparrow \right\rangle ,\left|{zx},\uparrow \right\rangle ,\left|{xy},\uparrow \right\rangle ,\left|{yz},\downarrow \right\rangle ,\left|{zx},\downarrow \right\rangle ,\left|{xy},\downarrow \right\rangle \right\}$, with yz, zx and xy representing the d_yz, d_zx, and d_xy orbitals, respectively, the electron in the bulk Ti t_2g band can be described by the effective ${\bf{k}}\cdot {\bf{p}}$ Hamiltonian [35, 38]

$$\begin{eqnarray}&&{H}_{{t}_{2g}}={{\mathbb{I}}}_{2}\otimes {h}_{{{\rm{\Gamma }}}_{25}^{+}},\end{eqnarray} \tag{ 2 }$$

where ${{\mathbb{I}}}_{2}$ is 2 × 2 unit matrix. The diagonal terms of ${h}_{{{\rm{\Gamma }}}_{25}^{+}}$ is given by $\langle {yz},\sigma | {h}_{{{\rm{\Gamma }}}_{25}^{+}}\left|{yz},\sigma \right\rangle ={{Lk}}_{x}^{2}+M({k}_{y}^{2}+{k}_{z}^{2}),$ and the off-diagonal terms are $\langle {yz},\sigma | {h}_{{{\rm{\Gamma }}}_{25}^{+}}\left|{zx},\sigma \right\rangle ={{Nk}}_{x}{k}_{y}$, with effective mass parameters L, M, N, and $\sigma =\uparrow ,\downarrow$. The other elements follow by exchanging the x, y, and z indices. In the above, the off-diagonal elements can be considered as the inter-orbital hopping in the bulk STO. The fit of the ${\bf{k}}\cdot {\bf{p}}$ model to the energy dispersion from our DFT calculations yield the mass parameters as $L=0.68\ \mathrm{eV}\,{\mathring{\rm A} }^{2},M=9.23\ \mathrm{eV}\,{\mathring{\rm A} }^{2}$, and $N=1.54\ \mathrm{eV}\,{\mathring{\rm A} }^{2}$(see appendix A for the details).

At the same time, in the presence of interface, the inversion symmetry is broken which leads to an additional Hamiltonian [34, 39–42]

$$\begin{eqnarray}{H}_{{\rm{ISB}}}={{\mathbb{I}}}_{2}\otimes \gamma \left[\begin{array}{ccc}0 & 0 & 2{{\rm{i}}{k}}_{x}\\ 0 & 0 & 2{{\rm{i}}{k}}_{y}\\ -2{{\rm{i}}{k}}_{x} & -2{{\rm{i}}{k}}_{y} & 0\end{array}\right],\end{eqnarray} \tag{ 3 }$$

up to the first order in k, where γ is the ISB parameter. In general, the ISB parameter is position-dependent, having a maximum value at the interface and rapidly decreases across at the deeper layers [41, 42].

The above Hamiltonian is symmetrical with respect to the xy, yz, and zx orbitals. However, this orbital degeneracy will be lifted in the presence of the atomic SOC. The SOC splits the six-fold degeneracy into a four-fold J = 3/2 multiplet of symmetry ${{\rm{\Gamma }}}_{8}^{+}$ and a twofold J = 1/2 multiplet of symmetry ${{\rm{\Gamma }}}_{7}^{+}$. In the original t_2g basis, the SOC is expressed as

$$\begin{eqnarray}{H}_{\mathrm{SO}}={\xi }_{\mathrm{SO}}\left[\begin{array}{cccccc}0 & {\rm{i}} & 0 & 0 & 0 & -1\\ -{\rm{i}} & 0 & 0 & 0 & 0 & {\rm{i}}\\ 0 & 0 & 0 & 1 & -{\rm{i}} & 0\\ 0 & 0 & 1 & 0 & -{\rm{i}} & 0\\ 0 & 0 & {\rm{i}} & {\rm{i}} & 0 & 0\\ -1 & -{\rm{i}} & 0 & 0 & 0 & 0\end{array}\right],\end{eqnarray} \tag{ 4 }$$

where the SO splitting is calculated via DFT to be ${{\rm{\Delta }}}_{\mathrm{SO}}=3{\xi }_{\mathrm{SO}}=29.85\,\mathrm{meV}$ (see appendix A).

We now turn to derive the Rashba split band structure of a STO-based heterostructure with 2DEG. For a STO film with finite thickness, we first find the eigen-states of (B2) at the Γ point ${H}_{0}={H}_{\mathrm{full}}({k}_{x}={k}_{y}=0)$. Subsequently, an effective Hamiltonian at other k values can be obtained by projecting the full Hamiltonian into the space of these eigenstates.

To simplify the analytical calculation, we transform the Hamiltonian into the basis $\left|J,{m}_{j}\right\rangle$ of ${{\rm{\Gamma }}}_{8}^{+}$ and ${{\rm{\Gamma }}}_{7}^{+}$ multiplets , ${{\rm{\Phi }}}_{J}=\left\{\left|\tfrac{1}{2},\tfrac{1}{2}\right\rangle ,\left|\tfrac{3}{2},\tfrac{1}{2}\right\rangle ,\left|\tfrac{3}{2},\tfrac{3}{2}\right\rangle ,\left|\tfrac{1}{2},-\tfrac{1}{2}\right\rangle ,\left|\tfrac{3}{2},-\tfrac{1}{2}\right\rangle ,\left|\tfrac{3}{2},-\tfrac{3}{2}\right\rangle \right\}$, in which the atomic SOC matrix in equation (4) is diagonalized (see appendix B). In the above basis, the Γ-point Hamiltonian is block-diagonalized and it can be decomposed into three subsets as ${H}_{0}={h}_{+}\oplus {h}_{0}\oplus {h}_{-}$, where

$$\begin{eqnarray}{h}_{\pm }=\left[\begin{array}{cc}a & \pm q\\ \pm q & b\end{array}\right],\ {h}_{0}=\left[\begin{array}{cc}c & 0\\ 0 & c\end{array}\right],\end{eqnarray} \tag{ 5 }$$

in which $a=\tfrac{1}{3}(L+2M){k}_{z}^{2}\,+\,2{\xi }_{\mathrm{SO}},b=\tfrac{1}{3}(2L+M){k}_{z}^{2}-{\xi }_{\mathrm{SO}}$, $c={{Mk}}_{z}^{2}-{\xi }_{\mathrm{SO}}$, and $q=\tfrac{\sqrt{2}}{3}(M-L){k}_{z}^{2}$.

For STO slab with finite thickness d, the electron states are quantized along the $z-$ direction, and we can make a substitution ${k}_{z}\to -{\rm{i}}{\partial }_{z}$. The eigenvectors and subband energies of electron can be found by solving the Schrödinger equation $({H}_{0}+V)\psi (z)=E\psi (z),$ in which ψ(z) is six-component eigenvector.

In solving the above Schrödinger equations in a STO slab with finite thickness d, we assume that the electron is confined in a quantum well with open boundaries, i.e. ψ(z = 0, d) = 0. This assumption may be understood by considering the character of the boundaries. On the one hand, at the LAO/STO interface, the La t_2g states in LAO layer lie far above the Ti t_2g states in the STO layer [5]. Therefore, the penetration of the Ti t_2g states into the LAO can be neglected [5] and one can assume an infinite barrier at the interface for the simplicity. On the other hand, the other surface of the STO slab is open to vacuum, and thus an infinite barrier potential is imposed. Similar boundary conditions are also applied for the δ-doped STO slab, where both STO surfaces are open to vacuum.

In the framework of variational approach, the solutions are given by $\left|{\psi }_{0}^{\sigma }\right\rangle ={f}_{0}(z)\left|\tfrac{3}{2},\sigma \tfrac{3}{2}\right\rangle ,$ and $\left|{\psi }_{\pm }^{\sigma }\right\rangle ={f}_{\pm }(z)\left(\cos {\chi }_{\pm }\left|\tfrac{1}{2},\sigma \tfrac{1}{2}\right\rangle +\sigma \sin {\chi }_{\pm }\left|\tfrac{3}{2},\sigma \tfrac{1}{2}\right\rangle \right),$ where the spatial envelope functions read as

$$\begin{eqnarray}&&{f}_{\tau }(z)={C}_{\tau }{{\rm{e}}}^{-{\beta }_{\tau }z}\sin \lambda z.\end{eqnarray} \tag{ 6 }$$

Here, $\lambda =n\tfrac{\pi }{d}$, C_τ is the normalization constant, and β_τ is the variational parameter for the minimization of the corresponding eigen-energy ${E}_{\tau }={E}_{\tau }^{(0)}+{{ \mathcal E }}_{z}F({\beta }_{\tau })$, where ${E}_{\tau }^{(0)}$ is the eigenvalue energy solution when V = 0 given by ${E}_{0}^{(0)}=M{\lambda }^{2}-{\xi }_{\mathrm{SO}},{E}_{\pm }^{(0)}=\tfrac{1}{2}\left[{\xi }_{\mathrm{SO}}+(M+L){\lambda }^{2}\right]\pm \tfrac{1}{2}{{\rm{\Delta }}}_{I}$, in which

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{I}={E}_{+}^{(0)}-{E}_{-}^{(0)}=\sqrt{{\eta }^{2}+8{\xi }_{\mathrm{SO}}^{2}},\end{eqnarray} \tag{ 7 }$$

with $\eta ={\xi }_{\mathrm{SO}}+(M-L){\lambda }^{2}$, and F(β) is the correction when $V={{ \mathcal E }}_{z}z$ (see appendix B for details).

Rearranging the sequence of the six eigenvectors as $\left\{{\psi }_{+}^{\uparrow },{\psi }_{+}^{\downarrow },{\psi }_{0}^{\uparrow },{\psi }_{0}^{\downarrow },{\psi }_{-}^{\uparrow },{\psi }_{-}^{\downarrow }\right\}$, and by projecting the full Hamiltonian to the spin-orbital space of these vectors, we obtain an effective Hamiltonian beyond the Γ-point in block form as (see appendix C for the full expression)

$$\begin{eqnarray}{ \mathcal H }=\left[\begin{array}{ccc}{{ \mathcal H }}_{+} & {{ \mathcal G }}_{3} & {{ \mathcal G }}_{1}\\ {{ \mathcal G }}_{3}^{\dagger } & {{ \mathcal H }}_{0} & {{ \mathcal G }}_{2}\\ {{ \mathcal G }}_{1}^{\dagger } & {{ \mathcal G }}_{2}^{\dagger } & {{ \mathcal H }}_{-}\end{array}\right].\end{eqnarray} \tag{ 8 }$$

In the above, the diagonal blocks are given by (for τ = 0, ±)

$$\begin{eqnarray}&&{{ \mathcal H }}_{\tau }={E}_{\tau }+\displaystyle \frac{{{\hslash }}^{2}{{\boldsymbol{k}}}^{2}}{2{m}_{\tau }^{* }}+{\alpha }_{\tau }({\boldsymbol{k}}\times {\boldsymbol{z}})\cdot {\boldsymbol{\sigma }},\end{eqnarray} \tag{ 9 }$$

which are characterised by bottom band edges E_τ, effective masses ${m}_{\tau }^{* }$, and Rashba parameters α_τ. Here, ${\boldsymbol{k}}=({k}_{x},{k}_{y})$ is the electron momentum in the x–y plane, ${\boldsymbol{\sigma }}$ is the vector of the Pauli matrices, and ${\boldsymbol{z}}$ is the unit vector perpendicular to the STO slab.

First, we look at the effective masses given by $\tfrac{1}{{m}_{0}^{* }}=\tfrac{(L+M)}{{{\hslash }}^{2}},\tfrac{1}{{m}_{\pm }^{* }}=\tfrac{1}{{m}_{0}^{* }}+\tfrac{(M-L)}{2{{\hslash }}^{2}}\left(1\pm \tfrac{\eta }{{{\rm{\Delta }}}_{I}}\right).$ Recall that M > L, and η/Δ_I < 1, we have ${m}_{\pm }^{* }\lt {m}_{0}^{* }$, which indicates that ${{ \mathcal H }}_{\pm }$ describe light-electron (denoted as ${{le}}^{\pm }$) bands, whereas ${{ \mathcal H }}_{0}$ corresponds to the heavy-electron (he) band. Furthermore, it can be verified that E₊ > E₀ > E₋, so that the two le-bands are sandwiched by the he-band. This arrangement leads to the anti-crossing between the bottom (le⁻)and the middle (he) bands (see figure 1(a)).

Zoom In Zoom Out Reset image size

3.1. Linear Rashba SOC

Now let us consider the expressions of linear RSOC given by

$$\begin{eqnarray}&&{\alpha }_{0}=0,\ {\alpha }_{\pm }=\pm \displaystyle \frac{4{\xi }_{\mathrm{SO}}{\bar{\gamma }}_{\pm }}{{{\rm{\Delta }}}_{I}},\end{eqnarray} \tag{ 10 }$$

where Δ_I is given in equation (7) regarded as effective crystal field, and

$$\begin{eqnarray}&&{\bar{\gamma }}_{\pm }=\int {\rm{d}}{z}\gamma (z)| {f}_{\pm }(z){| }^{2},\end{eqnarray} \tag{ 11 }$$

is considered as the effective inversion-symmetry breaking, with the integration taken over the STO slab. In the above, $| f(z){| }^{2}$ is the spatial distribution function, thus $\bar{\gamma }$ is the average ISB parameter across the STO slab. As the ISB constitutes an interface effect, it is only finite near the interface layer, and becomes vanishingly small in the bulk [41, 42]. Thus, equation (11) indicates that the localization of electron at the interface is significant for inducing large effective ISB, which can be achieved by having a large confinement effect induced by external electric field and/or charge density [43]. To obtain a simple analytical expression of the average ISB parameter, we consider the case where the ISB is finite only at the interface layer and zero elsewhere in the bulk, i.e. $\gamma (z)={\gamma }_{0}\theta (a-z)$, where θ(x) is the Heaviside step function, and a = 3.870 Å is the lattice constant of STO crystal. In this case, the average ISB parameter evaluates as

$$\begin{eqnarray}&&{\bar{\gamma }}_{\pm }{/\gamma }_{0}=\displaystyle \frac{{a}^{3}}{{M}_{\pm }}{{ \mathcal E }}_{z}+2{\pi }^{2}{\left(\displaystyle \frac{4}{3{M}_{\pm }^{2}}\right)}^{2/3}\displaystyle \frac{{a}^{3}}{{d}^{2}}{{ \mathcal E }}_{z}^{1/3},\end{eqnarray} \tag{ 12 }$$

in which ${M}_{\pm }=M/2\left(1\pm \tfrac{\eta }{{{\rm{\Delta }}}_{I}}\right)$. Equation (12) shows that the average ISB, and thus the RSOC, is scaled with the interface confining electric field ${{ \mathcal E }}_{z}$ to the first order as depicted in figure 2(c).

Zoom In Zoom Out Reset image size

Furthermore, it can be seen that reducing the QW width d (STO thickness) will increase the effective ISB due to the wavefunction localization near the interface. However, the effective crystal field in the thin film limit reads as

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{I}\approx \displaystyle \frac{(M-L){\pi }^{2}}{{d}^{2}},\end{eqnarray} \tag{ 13 }$$

which also increases with decreasing d. This trend is ascribed to the quantum confinement effect in the STO quantum well, which causes the subbands further split in the thin STO limit. From equation (10), these competing trends above between the effective ISB and effective crystal field lead to non-monotonic thickness dependence of the RSOC (see figure 2(a)). In general, the RSOC diminishes with reducing STO thickness, which is one of our main results.

On the other hand, at large STO thickness, the effective crystal field reduces to ${{\rm{\Delta }}}_{I}\approx 3{\xi }_{\mathrm{SO}}$, and ${\bar{\gamma }}_{\pm }\approx \tfrac{{a}^{3}{{ \mathcal E }}_{z}}{{M}_{\pm }}$, which result in a saturation of the linear Rashba to ${\alpha }_{\pm }\approx \tfrac{4{\gamma }_{0}{a}^{3}{{ \mathcal E }}_{z}}{3{M}_{\pm }}$ as shown in figure 2(a). It is interesting to see that in the large thickness limit, the saturation value of RSOC has a weak dependence on the atomic SOC ξ_SO. However, this does not deny the vital role of the atomic SOC, as the RSOC expression in equation (10) would vanish if ξ_SO is set to zero in the first place. Assuming typical parameter values, the linear Rashba coefficient is found to be in the range of 1–10 meV Å, corresponding to spin-splitting of 1–10 meV, which is consistent with previous computational [13, 33, 34] and experimental works [14–17].

We note that in other transition metal oxides such as ${\mathrm{SrIrO}}_{3}$ [44], the atomic SOC is much stronger than that in the ${\mathrm{SrTiO}}_{3}$. However, the crystal field in such oxides is also much larger at the same time [44]. Therefore, following equation (10), the strong atomic SOC oxides does not always translate into a strong Rashba SOC. Experimentally, it has been shown that the Rashba SOC in ${\mathrm{SrIrO}}_{3}$ is about 5 meV Å [45], a strength which is comparable to that in the ${\mathrm{SrTiO}}_{3}$.

3.2. Cubic Rashba SOC

Now we turn to discuss role of the off-diagonal elements ${ \mathcal G }$ in equation (8), which describe the coupling between the m_j = ±3/2 and m_j = ±1/2 states. It has been shown that the above coupling is the prerequisite for the presence of cubic SOC in the m_j = 3/2 states [23], which in our case is the he bands. Indeed, by reducing the 6 × 6 Hamiltonian (8) to three 2 × 2 matrices by means of quasi-degenerate partitioning [23], we obtain cubic corrections as

$$\begin{eqnarray}\begin{array}{rcl}\delta {{ \mathcal H }}_{{{le}}^{-}}^{\mathrm{cubic}} & = & {\beta }_{3}({k}_{x}^{2}-{k}_{y}^{2})({k}_{y}{\sigma }_{x}-{k}_{x}{\sigma }_{y})-{\eta }_{3}{k}_{x}{k}_{y}({k}_{x}{\sigma }_{x}-{k}_{y}{\sigma }_{y}),\end{array}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\delta {{ \mathcal H }}_{he}^{{\rm{c}}{\rm{u}}{\rm{b}}{\rm{i}}{\rm{c}}}={\beta }_{3}\left.\left({k}_{x}^{2}-\right.{k}_{y}^{2}\right)\,\left({k}_{y}{\sigma }_{x}-{k}_{x}{\sigma }_{y}\right)+{\eta }_{3}{k}_{x}{k}_{y}({k}_{x}{\sigma }_{x}-{k}_{y}{\sigma }_{y}),\end{eqnarray} \tag{ 15 }$$

in which the k-cubed SOC parameters are given by

$$\begin{eqnarray}&&{\beta }_{3}=\displaystyle \frac{(M-L){\xi }_{\mathrm{SO}}R}{{{\rm{\Delta }}}_{I}{{\rm{\Delta }}}_{-}},\ {\eta }_{3}=-\displaystyle \frac{2N{\xi }_{\mathrm{SO}}R}{{{\rm{\Delta }}}_{I}{{\rm{\Delta }}}_{-}},\end{eqnarray} \tag{ 16 }$$

in which $R=(2{\gamma }_{-0}-N{\kappa }_{-0}){\rho }_{-0},{{\rm{\Delta }}}_{-}={E}_{0}-{E}_{-}$, where we defined ${\gamma }_{{ij}}=\int \gamma (z){f}_{i}{f}_{j},{\kappa }_{{ij}}=\int {f}_{i}{\partial }_{z}{f}_{j}$, and ${\rho }_{{ij}}=\int {f}_{i}{f}_{j}$ (see appendix D). In equation (16), there are two cubic RSOC terms with different anisotropy in k-space. The first term depends on the mass difference, i.e. M – L, which is responsible for the mass anisotropy in the bulk STO. The second term stems from the the bulk inter-orbital hopping N, which is responsible for the asymmetry in the bulk. The cubic SOC parameters have a similar thickness dependence as in the linear case, in that they saturate at large thickness limit as shown in figure 2 (b). However, the cubic SOC diminishes with increasing electric field, in contrast to the opposite trend shown by the linear RSOC plotted in figures 2(c) and (d). With typical parameter values,, the cubic RSOC is estimated to be in the range 1−4 eV ų, which is in agreement with previous experimental and computational works [13, 19]. Furthermore, the strength of cubic SOC has been shown to decrease at high electron density [19], a trend which is consistent with our theory taking into consideration the linear correlation between the interface electric field and electron density [43].

3.3. Orbital polarization

In the above, we obtain three pairs of bands, each of which is a mix of d_xy, d_yz, and d_zx orbitals. To analyze the contribution of these orbitals in each band, we introduce the orbital occupancy which is obtained by projecting the eigenstates of the effective Hamiltonian (8) onto the t_2g basis states (see appendix E for details). Denote ${{ \mathcal N }}_{i}^{\tau }$ as the occupancy of orbital i = xy, yz, zx in band $\tau ={{le}}^{\pm },{he}$, we can define the orbital polarization as

$$\begin{eqnarray}&&{{ \mathcal P }}_{\tau }=\displaystyle \frac{{{ \mathcal N }}_{{xy}}^{\tau }-{{ \mathcal N }}_{{yz}+{zx}}^{\tau }}{{{ \mathcal N }}_{{xy}}^{\tau }+{{ \mathcal N }}_{{yz}+{zx}}^{\tau }},\end{eqnarray} \tag{ 17 }$$

which characterizes the contribution of the preferable d_xy in the interfacial 2DEG channel. In figure 1(e), the orbital polarization is depicted, which shows drastically different contribution of the orbitals in each band at different energy region. Around the Γ-point, the orbitals are well polarized with d_xy and d_yz/zx predominant in the two lowest bands, respectively. In this case, the orbital polarization can be obtained as ${{ \mathcal P }}_{{he}}=-1,{{ \mathcal P }}_{{{le}}^{\pm }}=\mp \tfrac{\eta }{{{\rm{\Delta }}}_{I}}.$ Close to the crossing region, stronger mixing between the d_xy-like and d_yz/zx-like bands will neutralize the orbital polarization. This trend is consistent with previous computational works [46, 47], where it has been shown that the d_xy orbital filling is less at higher energy region. At the same time, the relative contribution of the different orbitals can be controlled by tuning the heterostructure's parameters such as STO thickness. For example, the d_xy orbital contribution becomes more dominant in the lowest le⁻ band when the STO thickness is reduced, and approaches ${{ \mathcal P }}_{{{le}}^{-}}\approx 1$ in the thin limit, whereas its weight reduces as

$$\begin{eqnarray}&&{{ \mathcal P }}_{{{le}}^{-}}\approx \displaystyle \frac{1}{3}\left[1+\displaystyle \frac{8(M-L){\pi }^{2}}{9{\xi }_{\mathrm{SO}}{d}^{2}}\right],\end{eqnarray} \tag{ 18 }$$

in the thick limit.

The thickness dependence of the d_xy occupancy is in line with the diminishing of the RSOC in the thin limit. The RSOC originates from the ISB, which represents the hopping between the d_xy and ${d}_{{yz}/{zx}}$ at the interface as shown in equation (3). This means that only bands that occupied by both d_xy and d_yz/zx orbitals will show the linear RSOC. Therefore, in thin STO slabs, the d_yz/zx disappears in the le⁻ bands resulting in the vanishing of the hopping as well as the RSOC. Similarly, since there is no d_xy in the he bands, the linear RSOC is absent in these bands at all STO thickness. Instead, the he band shows cubic SOC as discussed in the previous section.

In this work, we formulate the analytical Hamiltonian of Rashba splitting in the t_2g bands of STO heterostructure based on ${\bf{k}}\cdot {\bf{p}}$ formalism. First, we express the Rashba parameter as functions of the heterostructure parameters and show that it is thickness dependent. Strikingly, the RSOC diminishes as the STO thickness reduces, and it saturates in thick STO layer. The hybridization between the orbitals results in orbital-dependent Rashba spin splittings. Explicitly, the linear Rashba SOC is associated with the d_xy orbital which preferably occupies the lowest band at low energy region. On the other hand, at higher energy crossing region, the cubic Rashba SOC accompanies the enriched d_yz/zx orbitals. These results reveal the capability to control the spin–orbit coupling and orbital selection in ${\mathrm{LaAlO}}_{3}/{\mathrm{SrTiO}}_{3}$ heterostructures for spintronics applications and devices.

MBAJ would like to acknowledge MOE Tier I (NUS Grant No. R-263-000-D66-114), MOE Tier II MOE2018-T2-2-117 (NUS Grant No. R-398-000-092-112) grants and NRF-CRP12-2013-01 (NUS Grant No. R-263-000-B30-281) for financial support. This work is also supported by the Ministry of Educations (MOE2015-T2-2-065, MOE2015-T2-1-099, MOE2015-T2-2-147, and MOE2017-T2-1-135), Singapore National Research Foundation under its Competitive Research Funding (NRF-CRP 8-2011-06 and NRF-CRP15-2015-01) and under its Medium Sized Centre Programme (Centre for Advanced 2D Materials and Graphene Research Centre), and FRC (R-144-000-379-114 and R-144-000-368-112).

To determine the ${\bf{k}}\cdot {\bf{p}}$ parameters, first-principles calculations were performed using density-functional theory (DFT) based Vienna ab initio simulation package (VASP5.4.4.18) with local density approximation (LDA) for the exchange-correlation interaction. Projector augmented wave potentials were selected to account for the interactions between electrons and ions. The LDA+U method was used with on-site effective U = 1.2 eV for the Ti d orbitals. The kinetic cut off energy for the expansion of plane-wavefunction was set to 500 eV. Monkhorst–Pack-based k-point grids for sampling the first Brillouin-zone were set to 8 × 8 × 8 and 8 × 8 × 1 for bulk STO and LAO/STO respectively. The Hellman–Feynman forces on each atom were minimized with a tolerance value of 0.01 eV Å^–1 for bulk STO, bulk LAO, and ${(\mathrm{LAO})}_{5.5}/{(\mathrm{STO})}_{30.5}$ heterostructure. Based on these settings, the equilibrium lattice constant of the bulk STO is calculated a = 3.870 Å, which is consistent with the experimental result (3.905 Å).

Without atomic SOC—Along k_[100], the dispersions are given by ${{Lk}}_{[100]}^{2},{{Mk}}_{[100]}^{2}$. By fitting the calculated energies along k_[100] direction to these ${\boldsymbol{k\cdot p}}$ formula (as shown in figure A1(a)), we can obtain the values of $L=0.68\ \mathrm{eV}\,{\mathring{\rm A} }^{2},M=9.23\ \mathrm{eV}\,{\mathring{\rm A} }^{2}$, respectively. Similarly, along k_[111] the dispersions are read as $1/3(L+2M-N){k}_{[111]}^{2},1/3(L+2M+2N){k}_{[111]}^{2}$. The different between these two bands are ${{Nk}}_{[111]}^{2}$, and thus by fitting the dataset to this expression, we obtain $N=1.54\ \mathrm{eV}\,{\mathring{\rm A} }^{2}$.

Zoom In Zoom Out Reset image size

With SOC—Similarly, in the presence of the atomic SOC, the three t_2g bands are no longer degenerate at the Γ-point, and the gap is determined by the SOC splitting Δ_SO = 29.85 meV (as shown in figure A1(b)).

We introduce $\left|J,{m}_{j}\right\rangle$ basis given by

$$\begin{eqnarray}\begin{array}{rcl}{u}_{1} & = & \left|\displaystyle \frac{1}{2},\displaystyle \frac{1}{2}\right\rangle =\displaystyle \frac{1}{\sqrt{3}}\left(-{\rm{i}}\left|{yz},\downarrow \right\rangle +\left|{zx},\downarrow \right\rangle -{\rm{i}}\left|{xy},\uparrow \right\rangle \right),\\ {u}_{2} & = & \left|\displaystyle \frac{3}{2},\displaystyle \frac{1}{2}\right\rangle =\displaystyle \frac{1}{\sqrt{6}}\left(-{\rm{i}}\left|{yz},\downarrow \right\rangle +\left|{zx},\downarrow \right\rangle +2{\rm{i}}\left|{xy},\uparrow \right\rangle \right),\\ {u}_{3} & = & \left|\displaystyle \frac{3}{2},-\displaystyle \frac{1}{2}\right\rangle =\displaystyle \frac{1}{\sqrt{2}}\left(-{\rm{i}}\left|{yz},\downarrow \right\rangle -\left|{zx},\downarrow \right\rangle \right),\\ {u}_{4} & = & \left|\displaystyle \frac{3}{2},\displaystyle \frac{3}{2}\right\rangle =\displaystyle \frac{1}{\sqrt{3}}\left(-{\rm{i}}\left|{yz},\uparrow \right\rangle -\left|{zx},\uparrow \right\rangle +{\rm{i}}\left|{xy},\downarrow \right\rangle \right),\\ {u}_{5} & = & \left|\displaystyle \frac{1}{2},-\displaystyle \frac{1}{2}\right\rangle =\displaystyle \frac{1}{\sqrt{6}}\left({\rm{i}}\left|{yz},\uparrow \right\rangle +\left|{zx},\uparrow \right\rangle +2{\rm{i}}\left|{xy},\downarrow \right\rangle \right),\\ {u}_{6} & = & \left|\displaystyle \frac{3}{2},-\displaystyle \frac{3}{2}\right\rangle =\displaystyle \frac{1}{\sqrt{2}}\left({\rm{i}}\left|{yz},\uparrow \right\rangle -\left|{zx},\uparrow \right\rangle \right),\end{array}\end{eqnarray} \tag{ B1 }$$

in which the atomic SOC Hamiltonian is diagonalized ${H}_{{\rm{SO}}}={\xi }_{{\rm{SO}}}\mathrm{diag}\,\left[2,-1,-1,2,-1,-1\right]$. The full Hamiltonian in this basis is given by

$$\begin{eqnarray}{H}_{\mathrm{full}}=\left[\begin{array}{cc}h(q) & {\rm{\Delta }}\\ {{\rm{\Delta }}}^{\dagger } & h{\left(-q\right)}^{* }\end{array}\right],\end{eqnarray} \tag{ B2 }$$

where

$$\begin{eqnarray}h(q)=\left[\begin{array}{ccc}a & q & r\\ q & b & r\\ {r}^{* } & {r}^{* } & c\end{array}\right],\ {\rm{\Delta }}=\gamma \left[\begin{array}{ccc}\tfrac{4}{3}{k}_{-} & \tfrac{{\rm{i}}\sqrt{2}}{3}{k}_{-} & \tfrac{{\rm{i}}\sqrt{2}}{\sqrt{3}}{k}_{+}\\ -\tfrac{{\rm{i}}\sqrt{2}}{3}{k}_{-} & \tfrac{4{\rm{i}}}{3}{k}_{-} & \tfrac{2{\rm{i}}}{\sqrt{3}}{k}_{+}\\ \tfrac{{\rm{i}}\sqrt{2}}{\sqrt{3}}{k}_{+} & \tfrac{2{\rm{i}}}{\sqrt{3}}{k}_{+} & 0\end{array}\right]+{{Nk}}_{z}\left[\begin{array}{ccc}0 & -\tfrac{{k}_{-}}{\sqrt{2}} & -\tfrac{{k}_{+}}{\sqrt{6}}\\ -\tfrac{{k}_{-}}{\sqrt{2}} & 0 & \tfrac{{k}_{+}}{\sqrt{3}}\\ -\tfrac{{k}_{+}}{\sqrt{6}} & -\tfrac{{k}_{+}}{\sqrt{3}} & 0\end{array}\right]\end{eqnarray} \tag{ B3 }$$

in which ${k}_{\pm }={k}_{x}\pm {{\rm{i}}{k}}_{y}$, and $a=\tfrac{1}{3}(L+2M)({k}_{x}^{2}+{k}_{y}^{2}+{k}_{z}^{2})+2{\xi }_{{\rm{SO}}}+\tfrac{{\xi }_{T}}{3}$, $b=\tfrac{1}{6}(L+5M)({k}_{x}^{2}+{k}_{y}^{2})+\tfrac{1}{3}(2L+M){k}_{z}^{2}-{\xi }_{{\rm{SO}}}+\tfrac{2{\xi }_{T}}{3}$, $c=\tfrac{1}{2}(L+M)({k}_{x}^{2}+{k}_{y}^{2})+{{Mk}}_{z}^{2}-{\xi }_{{\rm{SO}}}$, $q=\tfrac{1}{3\sqrt{2}}(L-M)({k}_{x}^{2}+{k}_{y}^{2}-2{k}_{z}^{2})-\tfrac{\sqrt{2}}{3}{\xi }_{T}$, and $r=\tfrac{1}{\sqrt{6}}(L-M)({k}_{x}^{2}-{k}_{y}^{2})-\tfrac{2{\rm{i}}}{\sqrt{6}}{{Nk}}_{x}{k}_{y}.$

Schrödinger equations—The Γ−point Hamiltonian is given by

$$\begin{eqnarray}&&{H}_{0}={H}_{\mathrm{full}}({k}_{x}={k}_{y}=0).\end{eqnarray} \tag{ B4 }$$

Following the decomposition of the above Hamiltonian, the six eigenstates can be written in general forms as

$$\begin{eqnarray}&&\left|{\psi }_{\pm }^{\sigma }\right\rangle ={f}_{1}^{\pm }(z)\left|\displaystyle \frac{1}{2},\sigma \displaystyle \frac{1}{2}\right\rangle +\sigma {f}_{2}^{\pm }(z)\left|\displaystyle \frac{3}{2},\sigma \displaystyle \frac{1}{2}\right\rangle ,\end{eqnarray} \tag{ B5 }$$

$$\begin{eqnarray}&&\left|{\psi }_{0}^{\sigma }\right\rangle ={f}_{0}(z)\left|\displaystyle \frac{3}{2},\sigma \displaystyle \frac{3}{2}\right\rangle ,\end{eqnarray} \tag{ B6 }$$

corresponding to eigen-energies E_± , E₀, respectively, in which σ = ± stand for $\uparrow ,\downarrow$, respectively. In the above, the spatial functions ${\left({f}_{1}^{-},{f}_{2}^{-}\right)}^{T}$ and ${\left({f}_{1}^{+},{f}_{2}^{+}\right)}^{T}$ are determined by solving the following coupled differential equations

$$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{1}{3}\left[\begin{array}{cc}-(L+2M){\partial }_{z}^{2}+(6{\xi }_{{\rm{SO}}}+{\xi }_{T}) & \sqrt{2}((L-M){\partial }_{z}^{2}-{\xi }_{T})\\ \sqrt{2}((L-M){\partial }_{z}^{2}-{\xi }_{T}) & -(2L+M){\partial }_{z}^{2}-(3{\xi }_{{\rm{SO}}}-2{\xi }_{T})\end{array}\right]\left(\begin{array}{c}{f}_{1}^{\pm }(z)\\ {f}_{2}^{\pm }(z)\end{array}\right)\\ \quad +V(z)\left(\begin{array}{c}{f}_{1}^{\pm }(z)\\ {f}_{2}^{\pm }(z)\end{array}\right)={E}_{\pm }\left(\begin{array}{c}{f}_{1}^{\pm }(z)\\ {f}_{2}^{\pm }(z)\end{array}\right),\end{array}\end{eqnarray} \tag{ B7 }$$

and f₀(z) is the solution of following differential equation

$$\begin{eqnarray}&&-{{Mf}}_{0}^{{\prime\prime} }(z)+V(z){f}_{0}(z)=({E}_{0}-{\xi }_{{SO}}){f}_{0}(z).\end{eqnarray} \tag{ B8 }$$

Boundary conditions—In solving the above Schrödinger equations in a STO slab with finite thickness d, we assume that the electron is confined in a quantum well with open boundaries, i.e. $\psi (z=0,d)=0$. This assumption may be understood by considering the character of the boundaries. On the one hand, at the LAO/STO interface, the La t_2g states in LAO layer lie far above the Ti t_2g states in the STO layer [5]. Therefore, the penetration of the Ti t_2g states into the LAO can be neglected [5] and one can assume an infinite barrier at the interface for the simplicity. On the other hand, the other surface of the STO slab is open to vacuum, and thus an infinite barrier potential is imposed. Similar boundary conditions are also applied for the δ-doped STO slab, where both STO surfaces are open to vacuum.

Solutions in zero field—First, we find the the solutions of (B7) and (B8) in the absence of the electric field, which can be found as

$$\begin{eqnarray}&&\left|{\psi }_{\pm }^{(0)},\sigma \right\rangle ={f}_{\pm }^{(0)}(z)\left(\cos {\chi }_{\pm }\left|\displaystyle \frac{1}{2},\sigma \displaystyle \frac{1}{2}\right\rangle +\sigma \sin {\chi }_{\pm }\left|\displaystyle \frac{3}{2},\sigma \displaystyle \frac{1}{2}\right\rangle \right),\end{eqnarray} \tag{ B9 }$$

$$\begin{eqnarray}&&\left|{\psi }_{0}^{(0)},\sigma \right\rangle ={f}_{0}^{(0)}(z)\left|\displaystyle \frac{3}{2},\sigma \displaystyle \frac{3}{2}\right\rangle ,\end{eqnarray} \tag{ B10 }$$

corresponding to eigen-energies

$$\begin{eqnarray}&&{E}_{\pm }^{(0)}=\displaystyle \frac{1}{2}\left({\xi }_{{\rm{SO}}}+{\xi }_{T}+(M+L){\lambda }^{2}\right)\pm \displaystyle \frac{1}{2}{\left[9{\xi }_{{\rm{SO}}}^{2}-2{\xi }_{{\rm{SO}}}({\xi }_{T}+(L-M){\lambda }^{2})+{\left({\xi }_{T}+(L-M){\lambda }^{2}\right)}^{2}\right]}^{1/2},\end{eqnarray} \tag{ B11 }$$

$$\begin{eqnarray}&&{E}_{0}^{(0)}=M{\lambda }^{2}-{\xi }_{{\rm{SO}}},\end{eqnarray} \tag{ B12 }$$

respectively. In the above, the spatial coefficients are simple sinusoidal functions as

$$\begin{eqnarray}&&{f}_{\pm }^{(0)}(z)={f}_{0}^{(0)}(z)=\sqrt{\displaystyle \frac{2}{d}}\sin \lambda z,\ \lambda =n\displaystyle \frac{\pi }{d},\end{eqnarray} \tag{ B13 }$$

and the spinor angles in equation (B9) are read as

$$\begin{eqnarray}&&{\chi }_{\pm }={\tan }^{-1}\left[\displaystyle \frac{2{\xi }_{{SO}}+{\xi }_{T}+{\lambda }^{2}(L+2M)-3{E}_{\pm }^{(0)}}{\sqrt{2}({\xi }_{T}+(L-M){\lambda }^{2})}\right].\end{eqnarray} \tag{ B14 }$$

It is easy to verify that $\chi ={\chi }_{+}={\chi }_{-}+\tfrac{\pi }{2}$, which secures the orthogonality of the eigen-wavefunctions $\langle {\psi }_{\tau }^{(0)},\sigma \left|{\psi }_{\tau ^{\prime} }^{(0)},\sigma ^{\prime} \right\rangle ={\delta }_{\tau \tau ^{\prime} }{\delta }_{\sigma \sigma ^{\prime} }$. Furthermore, this angle is in the range $0\lt \chi \lt \arctan \tfrac{1}{\sqrt{2}}$.

Solutions under applied field—In the presence of the confinement potential $V(z)={{ \mathcal E }}_{z}z$, with ${{ \mathcal E }}_{z}$ being the electric field, the solutions can be found via variational method [48], in which the eigen wavefunctions are expressed as

$$\begin{eqnarray}&&\left|{\psi }_{\pm }^{\sigma }\right\rangle ={f}_{\pm }(z)\left(\cos {\chi }_{\pm }\left|\displaystyle \frac{1}{2},\sigma \displaystyle \frac{1}{2}\right\rangle +\sigma \sin {\chi }_{\pm }\left|\displaystyle \frac{3}{2},\sigma \displaystyle \frac{1}{2}\right\rangle \right),\end{eqnarray} \tag{ B15 }$$

$$\begin{eqnarray}&&\left|{\psi }_{0}^{\sigma }\right\rangle ={f}_{0}(z)\left|\displaystyle \frac{3}{2},\sigma \displaystyle \frac{3}{2}\right\rangle ,\end{eqnarray} \tag{ B16 }$$

in which

$$\begin{eqnarray}&&{f}_{\tau }(z)={C}_{\tau }{{\rm{e}}}^{-{\beta }_{\tau }z/d}{f}_{\tau }^{(0)}(z),\end{eqnarray} \tag{ B17 }$$

where C_τ are the normalization constants given by

$$\begin{eqnarray}&&{C}_{\tau }={\left[\displaystyle \frac{2{\beta }_{\tau }({\pi }^{2}+{\beta }_{\tau }^{2})}{{\pi }^{2}(1-{{\rm{e}}}^{-2{\beta }_{\tau }})}\right]}^{1/2}\end{eqnarray} \tag{ B18 }$$

and β_τ are the variational parameters for minimization of the corresponding energies

$$\begin{eqnarray}&&{E}_{\tau }={E}_{\tau }^{(0)}+{{ \mathcal E }}_{z}{\rm{d}}{F}({\beta }_{\tau }),\end{eqnarray} \tag{ B19 }$$

in which

$$\begin{eqnarray}&&F({\beta }_{\tau })=\displaystyle \frac{{M}_{\tau }{\beta }_{\tau }^{2}}{{{ \mathcal E }}_{z}{d}^{3}}+\displaystyle \frac{1}{2}\left(1+\displaystyle \frac{1}{{\beta }_{\tau }}+\displaystyle \frac{2{\beta }_{\tau }}{{\pi }^{2}+{\beta }_{\tau }^{2}}-\coth {\beta }_{\tau }\right),\end{eqnarray} \tag{ B20 }$$

where M₀ = M, and ${M}_{\pm }=\tfrac{1}{3}(L+2M){\cos }^{2}{\chi }_{\pm }+\tfrac{1}{3}(2L+M){\sin }^{2}{\chi }_{\pm }+\tfrac{\sqrt{2}}{3}(M-L)\sin 2{\chi }_{\pm }$.

The variational parameters can be obtained by minimizing the above energies, which increase as the electric field and the well width increase, as depicted in figure B1. In the limit of weak electrostatic energy ${\delta }_{\tau }=\tfrac{{\pi }^{2}{{ \mathcal E }}_{z}d}{{M}_{\tau }{\lambda }^{2}}\ll 1$, the variational parameter and corresponding energy can be approximated as

$$\begin{eqnarray}&&{\beta }_{\tau }=\displaystyle \frac{{\pi }^{2}-6}{12{\pi }^{2}}{\delta }_{\tau },\ {E}_{\tau }={E}_{\tau }^{(0)}+\displaystyle \frac{{{ \mathcal E }}_{z}d}{2}-({M}_{\tau }{\lambda }^{2})\displaystyle \frac{{\left({\pi }^{2}-6\right)}^{2}}{144{\pi }^{2}}{\delta }_{\tau }^{2}.\end{eqnarray} \tag{ B21 }$$

On the other hand, in the strong electrostatic energy limit δ ≫ 1, we have

$$\begin{eqnarray}&&{\beta }_{\tau }=\displaystyle \frac{{3}^{1/3}{\delta }_{\tau }^{1/3}}{{2}^{2/3}},\ {E}_{\tau }={E}_{\tau }^{(0)}+({M}_{\tau }{\lambda }^{2})\displaystyle \frac{3}{4}{\left(6{\delta }_{\tau }\right)}^{2/3}.\end{eqnarray} \tag{ B22 }$$

Zoom In Zoom Out Reset image size

Having derived the general expression of the eigenvectors at the Γ-point, we rearrange their sequence as

$$\begin{eqnarray}&&\left\{{\psi }_{+}^{\uparrow },{\psi }_{+}^{\downarrow },{\psi }_{0}^{\uparrow },{\psi }_{0}^{\downarrow },{\psi }_{-}^{\uparrow },{\psi }_{-}^{\downarrow }\right\}.\end{eqnarray} \tag{ C1 }$$

Map the full Hamiltonian (1) to the spin-orbital space of these vectors as ${ \mathcal H }={\int }_{0}^{d}{\rm{d}}{z}{\left({\psi }_{\tau }^{\sigma }\right)}^{* }H{\psi }_{\tau ^{\prime} }^{\sigma ^{\prime} }$, we obtain an effective Hamiltonian beyond the Γ-point, which can be written in block form as

$$\begin{eqnarray}{ \mathcal H }=\left[\begin{array}{ccc}{{ \mathcal H }}_{+} & {{ \mathcal G }}_{3} & {{ \mathcal G }}_{1}\\ {{ \mathcal G }}_{3}^{\dagger } & {{ \mathcal H }}_{0} & {{ \mathcal G }}_{2}\\ {{ \mathcal G }}_{1}^{\dagger } & {{ \mathcal G }}_{2}^{\dagger } & {{ \mathcal H }}_{-}\end{array}\right],\end{eqnarray} \tag{ C2 }$$

C.1. Diagonal blocks

The diagonal blocks of (C2) are given by (for τ = 0, ±)

$$\begin{eqnarray}&&{{ \mathcal H }}_{\tau }={E}_{\tau }+\displaystyle \frac{{{\hslash }}^{2}{k}^{2}}{2{m}_{\tau }^{* }}+{\alpha }_{\tau }({k}_{y}{\sigma }_{x}-{k}_{x}{\sigma }_{y}).\end{eqnarray} \tag{ C3 }$$

Rashba spin–orbit coupling

$$\begin{eqnarray}&&{\alpha }_{0}=0,\end{eqnarray} \tag{ C4 }$$

$$\begin{eqnarray}&&{\alpha }_{\pm }=\pm \displaystyle \frac{{\bar{\gamma }}_{\pm }}{3}(4\cos 2\chi -\sqrt{2}\sin 2\chi )=\pm \displaystyle \frac{4{\xi }_{{\rm{SO}}}{\bar{\gamma }}_{\pm }}{{{\rm{\Delta }}}_{I}},\end{eqnarray} \tag{ C5 }$$

where ${\bar{\gamma }}_{\pm }=\int {\rm{d}}{z}\gamma (z)| {f}_{\pm }{| }^{2}$, and have substituted χ = χ₊ so that $(4\cos 2\chi -\sqrt{2}\sin 2\chi )=4{\xi }_{{\rm{SO}}}/{{\rm{\Delta }}}_{I}$, with

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{I}={E}_{+}^{(0)}-{E}_{-}^{(0)}=\sqrt{{\eta }^{2}+8{\xi }_{{\rm{SO}}}^{2}},\end{eqnarray} \tag{ C6 }$$

with $\eta ={\xi }_{{\rm{SO}}}+(M-L){\lambda }^{2}$.

Effective mass

$$\begin{eqnarray}&&\displaystyle \frac{1}{{m}_{0}^{* }}=\displaystyle \frac{1}{{{\hslash }}^{2}}(L+M),\end{eqnarray} \tag{ C7 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{1}{{m}_{\pm }^{* }} & = & \displaystyle \frac{1}{2{{\hslash }}^{2}}(L+3M)\pm \displaystyle \frac{1}{6{{\hslash }}^{2}}(M-L)\left(2\sqrt{2}\sin (2\chi )+\cos (2\chi )\right),\\ & = & \displaystyle \frac{1}{{m}_{0}^{* }}+\displaystyle \frac{(M-L)}{2{{\hslash }}^{2}}\left(1\pm \displaystyle \frac{\eta }{{{\rm{\Delta }}}_{I}}\right).\end{array}\end{eqnarray} \tag{ C8 }$$

Effective ISB parameter—As the ISB value rapidly decays in the second and deeper layers, we assume that it is only finite in the interface layer and zero in the bulk, i.e. $\gamma (z)={\gamma }_{0}\theta (a-z)$, with a being the lattice constant. In this case, the effective ISB parameter is read as

$$\begin{eqnarray}&&{\bar{\gamma }}_{\pm }={\gamma }_{0}\displaystyle \frac{\left(\coth \left({\beta }_{\pm }\right)+1\right)}{2{\pi }^{2}}\left({\pi }^{2}-{{\rm{e}}}^{-2x{\beta }_{\pm }}\left(2{\beta }_{\pm }^{2}{\sin }^{2}(\pi x)+\pi {\beta }_{\pm }\sin (2\pi x)+{\pi }^{2}\right)\right),\ x=\displaystyle \frac{a}{d}.\end{eqnarray} \tag{ C9 }$$

which can be obtained once the variational parameters are known. In the the weak field regime, by substituting equation (B21) the above can be explicitly expressed as

$$\begin{eqnarray}&&{\bar{\gamma }}_{\pm }{/\gamma }_{0}=\displaystyle \frac{2{\pi }^{2}{a}^{3}}{3{d}^{3}}+\displaystyle \frac{\left({\pi }^{2}-6\right){a}^{3}{{ \mathcal E }}_{z}}{18{M}_{\pm }},\end{eqnarray} \tag{ C10 }$$

while in the strong field limit, we have

$$\begin{eqnarray}&&{\bar{\gamma }}_{\pm }{/\gamma }_{0}=\displaystyle \frac{{a}^{3}}{{M}_{\pm }}{{ \mathcal E }}_{z}+2{\pi }^{2}{\left(\displaystyle \frac{4}{3{M}_{\pm }^{2}}\right)}^{2/3}\displaystyle \frac{{a}^{3}}{{d}^{2}}{{ \mathcal E }}_{z}^{1/3}.\end{eqnarray} \tag{ C11 }$$

C.2. Off-diagonal blocks

The off-diagonal blocks of (C2) are given by

$$\begin{eqnarray}&&{{ \mathcal G }}_{1}={\epsilon }_{1}(k)+{\delta }_{1}\left({k}_{y}{\sigma }_{x}-{k}_{x}{\sigma }_{y}\right),\end{eqnarray} \tag{ C12 }$$

$$\begin{eqnarray}&&{{ \mathcal G }}_{2}={\epsilon }_{2}(k)+{\delta }_{2}\left({k}_{x}{\sigma }_{x}-{k}_{y}{\sigma }_{y}\right),\end{eqnarray} \tag{ C13 }$$

$$\begin{eqnarray}&&{{ \mathcal G }}_{3}={\epsilon }_{3}(k)+{\delta }_{3}\left({k}_{x}{\sigma }_{x}-{k}_{y}{\sigma }_{y}\right),\end{eqnarray} \tag{ C14 }$$

where the hybridization parameters are

$$\begin{eqnarray}&&{\delta }_{1}=\displaystyle \frac{1}{6}\left(3\sqrt{2}N{\kappa }_{+-}+2{\gamma }_{+-}\left(4\sin (2\chi )+\sqrt{2}\cos (2\chi )\right)\right),\end{eqnarray} \tag{ C15 }$$

$$\begin{eqnarray}&&{\delta }_{2}=-\displaystyle \frac{{\rm{i}}}{2\sqrt{3}}\left(\sqrt{2}\cos (\chi )+2\sin (\chi )\right)\left(2{\gamma }_{-0}-N{\kappa }_{-0}\right),\end{eqnarray} \tag{ C16 }$$

$$\begin{eqnarray}&&{\delta }_{3}=-\displaystyle \frac{{\rm{i}}}{2\sqrt{3}}\left(\sqrt{2}\cos (\chi )-2\sin (\chi )\right)\left(2{\gamma }_{+0}-N{\kappa }_{+0}\right),\end{eqnarray} \tag{ C17 }$$

and

$$\begin{eqnarray}&&{\epsilon }_{1}=\displaystyle \frac{{\rho }_{+-}}{12}(M-L)\left({k}_{x}^{2}+{k}_{y}^{2}\right)\left(2\sqrt{2}\cos (2\chi )-\sin (2\chi )\right),\end{eqnarray} \tag{ C18 }$$

$$\begin{eqnarray}&&{\epsilon }_{2}=\displaystyle \frac{{\rho }_{-0}}{2\sqrt{3}}\left(\cos (\chi )-\sqrt{2}\sin (\chi )\right)\left((M-L){\sigma }_{z}\left({k}_{x}^{2}-{k}_{y}^{2}\right)+2{{\rm{i}}{N}{k}}_{x}{k}_{y}\right),\end{eqnarray} \tag{ C19 }$$

$$\begin{eqnarray}&&{\epsilon }_{3}=-\displaystyle \frac{{\rho }_{+0}}{2\sqrt{3}}\left(\sin (\chi )+\sqrt{2}\cos (\chi )\right)\left((M-L){\sigma }_{z}\left({k}_{x}^{2}-{k}_{y}^{2}\right)+2{{\rm{i}}{N}{k}}_{x}{k}_{y}\right).\end{eqnarray} \tag{ C20 }$$

In the above, we have defined the integrals ${\rho }_{\tau \tau ^{\prime} }=\int {{\rm{d}}{z}{f}}_{\tau }{f}_{\tau ^{\prime} },{\kappa }_{\tau \tau ^{\prime} }=\int {{\rm{d}}{z}{f}}_{\tau }{\partial }_{z}{f}_{\tau ^{\prime} },{\gamma }_{\tau \tau ^{\prime} }=\int {\rm{d}}{z}\gamma (z){f}_{\tau }{f}_{\tau ^{\prime} }$. Explicitly, they are given as

$$\begin{eqnarray}&&{\rho }_{\tau \tau ^{\prime} }={C}_{\tau }{C}_{\tau ^{\prime} }\displaystyle \frac{4(1-{{\rm{e}}}^{-{\beta }_{\tau \tau ^{\prime} }}){\pi }^{2}}{{\beta }_{\tau \tau ^{\prime} }(4{\pi }^{2}+{\beta }_{\tau \tau ^{\prime} }^{2})},\ {\beta }_{\tau \tau ^{\prime} }={\beta }_{\tau }+{\beta }_{\tau ^{\prime} }\end{eqnarray} \tag{ C21 }$$

$$\begin{eqnarray}&&{\kappa }_{\tau \tau ^{\prime} }={\rho }_{\tau \tau ^{\prime} }\displaystyle \frac{{\beta }_{\tau }-{\beta }_{\tau ^{\prime} }}{2d},\end{eqnarray} \tag{ C22 }$$

$$\begin{eqnarray}&&{\gamma }_{\tau \tau ^{\prime} }\approx {\rho }_{\tau \tau ^{\prime} }\displaystyle \frac{{\beta }_{\tau \tau ^{\prime} }(4{\pi }^{2}+{\beta }_{\tau \tau ^{\prime} }^{2})}{5(1-{{\rm{e}}}^{-{\beta }_{\tau \tau ^{\prime} }})}{x}^{3},\ x=\displaystyle \frac{a}{d}\ll 1.\end{eqnarray} \tag{ C23 }$$

The Hamiltonian (C2) can be decomposed into three subspaces via the Lowdin partitioning [23], i.e.

$$\begin{eqnarray}{ \mathcal H }=\left[\begin{array}{ccc}{{ \mathcal H }}_{+} & {{ \mathcal G }}_{3} & {{ \mathcal G }}_{1}\\ {{ \mathcal G }}_{3}^{\dagger } & {{ \mathcal H }}_{0} & {{ \mathcal G }}_{2}\\ {{ \mathcal G }}_{1}^{\dagger } & {{ \mathcal G }}_{2}^{\dagger } & {{ \mathcal H }}_{-}\end{array}\right]\to \left[\begin{array}{ccc}{{ \mathcal H }}_{{{le}}^{+}} & 0 & 0\\ 0 & {{ \mathcal H }}_{{he}} & 0\\ 0 & 0 & {{ \mathcal H }}_{{{le}}^{-}}\end{array}\right],\end{eqnarray} \tag{ D1 }$$

where the effective Hamiltonian of the each subspace is derived as

$$\begin{eqnarray}&&{{ \mathcal H }}_{{he}}={{ \mathcal H }}_{0}+\displaystyle \frac{{{ \mathcal G }}_{2}\cdot {{ \mathcal G }}_{2}^{\dagger }}{{E}_{0}-{E}_{-}}+\displaystyle \frac{{{ \mathcal G }}_{3}^{\dagger }\cdot {{ \mathcal G }}_{3}}{{E}_{0}-{E}_{+}},\end{eqnarray} \tag{ D2 }$$

$$\begin{eqnarray}&&{{ \mathcal H }}_{{{le}}^{+}}={{ \mathcal H }}_{+}+\displaystyle \frac{{{ \mathcal G }}_{1}\cdot {{ \mathcal G }}_{1}^{\dagger }}{{E}_{+}-{E}_{-}}+\displaystyle \frac{{{ \mathcal G }}_{3}\cdot {{ \mathcal G }}_{3}^{\dagger }}{{E}_{+}-{E}_{0}},\end{eqnarray} \tag{ D3 }$$

$$\begin{eqnarray}&&{{ \mathcal H }}_{{{le}}^{-}}={{ \mathcal H }}_{-}+\displaystyle \frac{{{ \mathcal G }}_{1}^{\dagger }\cdot {{ \mathcal G }}_{1}}{{E}_{-}-{E}_{+}}+\displaystyle \frac{{{ \mathcal G }}_{2}^{\dagger }\cdot {{ \mathcal G }}_{2}}{{E}_{-}-{E}_{0}}.\end{eqnarray} \tag{ D4 }$$

Substituting the off-diagonal matrix ${{ \mathcal G }}_{i}$ into the above, and to the leading order in the ξ_SO, we obtain the k-cubed corrections as

$$\begin{eqnarray}&&\delta {{ \mathcal H }}_{\tau }={\beta }_{3}^{\tau }({k}_{x}^{2}-{k}_{y}^{2})({k}_{y}{\sigma }_{x}-{k}_{x}{\sigma }_{y})+{\eta }_{3}^{\tau }{k}_{x}{k}_{y}({k}_{x}{\sigma }_{x}-{k}_{y}{\sigma }_{y}),\end{eqnarray} \tag{ D5 }$$

for $\tau ={he},{{le}}^{\pm }$, in which the cubic SOC parameters are derived as

$$\begin{eqnarray}&&{\beta }_{3}^{{he}}={\beta }_{3}^{{{le}}^{+}}+{\beta }_{3}^{{{le}}^{-}},{\eta }_{2}^{{he}}=-({\eta }_{3}^{{{le}}^{+}}+{\eta }_{2}^{{{le}}^{-}}),\end{eqnarray} \tag{ D6 }$$

$$\begin{eqnarray}&&{\beta }_{3}^{{le}+}=\displaystyle \frac{(M-L){\xi }_{{\rm{SO}}}}{{{\rm{\Delta }}}_{I}}\displaystyle \frac{{R}_{+}}{{{\rm{\Delta }}}_{+}},\ {\eta }_{3}^{{le}+}=\displaystyle \frac{2N{\xi }_{{\rm{SO}}}}{{{\rm{\Delta }}}_{I}}\displaystyle \frac{{R}_{+}}{{{\rm{\Delta }}}_{+}},\end{eqnarray} \tag{ D7 }$$

$$\begin{eqnarray}&&{\beta }_{3}^{{le}-}=\displaystyle \frac{(M-L){\xi }_{{\rm{SO}}}}{{{\rm{\Delta }}}_{I}}\displaystyle \frac{{R}_{-}}{{{\rm{\Delta }}}_{-}},\ {\eta }_{3}^{{le}-}=\displaystyle \frac{2N{\xi }_{{\rm{SO}}}}{{{\rm{\Delta }}}_{I}}\displaystyle \frac{{R}_{-}}{{{\rm{\Delta }}}_{-}},\end{eqnarray} \tag{ D8 }$$

in which ${R}_{+}=(2{\gamma }_{+0}-N{\kappa }_{+0}){\rho }_{+0},\ {R}_{-}=(2{\gamma }_{-0}-N{\kappa }_{-0}){\rho }_{-0},\ {{\rm{\Delta }}}_{+}={E}_{+}-{E}_{0},\ {{\rm{\Delta }}}_{-}=({E}_{0}-{E}_{-})$. As the le⁺ band lie far above the he and le⁻ bands, its corresponding wavefunction is more spreading into the bulk. As a consequence, one can show that ${\gamma }_{+0}\ll {\gamma }_{-0},{\rho }_{+0}\ll {\rho }_{-0}$, and ${\kappa }_{+0}\ll {\kappa }_{-0}$, so that ${P}_{+}\ll {P}_{-}$. Therefore, in equations (D6) we retain only the second term,i.e., ${\beta }_{3}^{{he}}={\beta }_{3}^{{{le}}^{-}},{\eta }_{2}^{{he}}=-{{eta}}_{2}^{{{le}}^{-}}$.

In the weak electric field regime, we have following approximations to the leading order in the E-field

$$\begin{eqnarray}&&{R}_{\pm }\approx \displaystyle \frac{4{\pi }^{2}{\gamma }_{0}{a}^{3}}{3{d}^{3}}+{{ \mathcal E }}_{z}\left[\displaystyle \frac{{\gamma }_{0}{a}^{3}({\pi }^{2}-6)}{18}\left(\displaystyle \frac{1}{{M}_{0}}+\displaystyle \frac{1}{{M}_{\pm }}\right)+\displaystyle \frac{{{Nd}}^{2}({\pi }^{2}-6)}{24{\pi }^{2}}\left(\displaystyle \frac{1}{{M}_{\pm }}-\displaystyle \frac{1}{{M}_{0}}\right)\right],\end{eqnarray} \tag{ D9 }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\pm }\approx {{\rm{\Delta }}}_{\pm }^{(0)}\pm {{ \mathcal E }}_{z}^{2}\displaystyle \frac{{\left({\pi }^{2}-6\right)}^{2}{d}^{4}}{144{\pi }^{4}}\left(\displaystyle \frac{1}{{M}_{0}}-\displaystyle \frac{1}{{M}_{\pm }}\right),\end{eqnarray} \tag{ D10 }$$

where ${{\rm{\Delta }}}_{\pm }^{(0)}=\pm ({E}_{\pm }^{(0)}-{E}_{0}^{(0)})$. The above equations show the dependence of the cubic RSOC on the electric field as $\sim \tfrac{{R}_{\pm }}{{{\rm{\Delta }}}_{\pm }}\sim \tfrac{1}{{{ \mathcal E }}_{z}}$.

The occupancy of the orbitals in each band can be obtained by representing the eigenstates of the effective Hamiltonian in terms of the t_2g basis as

$$\begin{eqnarray}&&\left|{{\rm{\Psi }}}_{\tau }^{\sigma }\right\rangle =\displaystyle \sum _{i}{c}_{i}^{\tau }\left|i\right\rangle \end{eqnarray} \tag{ E1 }$$

for τ = 0, ±, and i = xy, yz, zx. Then ${{ \mathcal N }}_{i}^{\tau }=| {c}_{i}^{\tau }{| }^{2}$ is referred to as the occupancy of the i-orbital in the τ bands. Around the Γ-point, the occupancies can be analytically obtained as

$$\begin{eqnarray}&&{{ \mathcal N }}_{{xy}}^{{he}}=0,\ {{ \mathcal N }}_{{yz}}^{{he}}={{ \mathcal N }}_{{zx}}^{{he}}=\displaystyle \frac{1}{2},\end{eqnarray} \tag{ E2 }$$

$$\begin{eqnarray}&&{{ \mathcal N }}_{{xy}}^{{{le}}^{-}}=\displaystyle \frac{1}{3}{\left(\sqrt{2}\cos \chi +\sin \chi \right)}^{2},\ {{ \mathcal N }}_{{yz}}^{{{le}}^{-}}={{ \mathcal N }}_{{zx}}^{{{le}}^{-}}=\displaystyle \frac{1}{6}{\left(\cos \chi -\sqrt{2}\sin \chi \right)}^{2},\end{eqnarray} \tag{ E3 }$$

$$\begin{eqnarray}&&{{ \mathcal N }}_{{xy}}^{{{le}}^{+}}=\displaystyle \frac{1}{3}{\left(\cos \chi -\sqrt{2}\sin \chi \right)}^{2},\ {{ \mathcal N }}_{{yz}}^{{{le}}^{+}}={{ \mathcal N }}_{{zx}}^{{{le}}^{+}}=\displaystyle \frac{1}{6}{\left(\sqrt{2}\cos \chi +\sin \chi \right)}^{2},\end{eqnarray} \tag{ E4 }$$

where $\tau ={he},{{le}}^{\pm }$ are equivalent to 0, ± bands, respectively. It can be verified that the orbital occupancies satisfy ${\sum }_{i}{{ \mathcal N }}_{i}^{\tau }=1$ for each band.

Substitute χ = χ₊ given in (B14), the above expressions reduce to

$$\begin{eqnarray}&&{{ \mathcal N }}_{{xy}}^{{he}}=0,\ {{ \mathcal N }}_{{xy}}^{{{le}}^{\pm }}=\displaystyle \frac{1}{2}\left(1\mp \displaystyle \frac{\eta }{{{\rm{\Delta }}}_{I}}\right),\end{eqnarray} \tag{ E5 }$$

$$\begin{eqnarray}&&{{ \mathcal N }}_{{yz}/{zx}}^{\tau }=\displaystyle \frac{1}{2}(1-{{ \mathcal N }}_{{xy}}^{\tau }).\end{eqnarray} \tag{ E6 }$$