Spin–orbit coupling in LaAlO₃/SrTiO₃ interfaces: magnetism and orbital ordering

Rashba spin–orbit-coupling (SOC) effects have been mostly studied in weakly interacting systems such as semiconductor heterostructures designed for spintronics applications [1], and occur in two-dimensional (2D) systems without mirror symmetry [2]. However, the effects of Rashba SOC in a 2D system with strongly interacting electrons, found for example in interfaces, are emerging as a new frontier. There is thus a pursuit for new emergent phases of matter in this regime, both theoretically [3] and experimentally [4], with the electron gas at the interface between the two non-magnetic insulators LaAlO₃ and SrTiO₃ (LAO/STO) the widest-studied such example. The observed ferromagnetic instability at this interface [5] could well be of Stoner type, which suggests that electronic correlations may be enhanced due to low dimensionality and poor screening at low densities. Hence, combined with its demonstrated tunability [6, 7], the LAO/STO interface is an ideal testbed for physics of Rashba SOC in correlated electron systems.

Recent observation of magnetic anisotropy may signal further richness in the phase diagram of the interface. Specifically, Bert et al [8] and Li et al [9] observed strong in-plane preference for magnetization (see figure 1(a)). Bert et al attributed this observed anisotropy to the shape anisotropy of the interface: energetic bias toward a certain magnetization direction e.g. along the longest axis of an ellipsoid, driven by an anisotropic demagnetization field. However, for ultra thin films consisting of only a few atomic layers this effect is subdominant next to microscopic effects [10, 11].

Zoom In Zoom Out Reset image size

If the shape anisotropy is negligible in the interfaces, the dominant source of magnetic anisotropy would be SOC. However, the typical alignment of the spin with the largest angular momentum in itinerant d-electron systems [12] would predict an out-of-plane magnetization upon occupation of d_xz and d_yz orbitals, in disagreement with experiments [8, 9, 13]. In this work, we show that Rashba SOC leads to the unusual circumstance of an anisotropic (spin) susceptibility. Assuming Stoner ferromagnetism close to a van Hove singularity near band edges due to SOC, we argue that this in turn leads to a magnetization anisotropy.

We further investigate the possibility of nematic order in the form of orbital ordering between the d_xz and d_yz orbitals, since the large density of states near the band edge also leads to an enhanced tendency toward such ordering. This order can couple to an in-plane magnetization anisotropy and we show here that a consequence of the Rashba SOC is a strong such coupling in the low-density regime near the band edge. This implies that orbital ordering accompanies the magnetic phase and leads to an additional magnetization anisotropy within the plane.

We use a three-band model for the Ti t_2g orbitals in the xy plane to describe the electronic structure of the interface [14]. For now, we ignore the atomic SOC in order to gain more analytic insight. In the presence of an external magnetic field $\vec {H}$, the Hamiltonian reads

$$\begin{equation} \mathcal{H} = \mathcal{H}_0 + \mathcal{H}_{\mathrm{soc}}^{\mathrm{R}} -\mu_{\mathrm{B}}\vec{H}\cdot\vec{S}. \end{equation} \tag{ 1 }$$

Here, $\mathcal {H}_0$ is the hopping Hamiltonian

$$\begin{equation} \mathcal{H}_0=\sum_{l,{\mathbf{k}},s} \xi_{l{\mathbf{k}}}^{(0)} c_{l{\mathbf{k}} s}^\dagger c^{\phantom{\dag}}_{l{\mathbf{k}} s} \end{equation} \tag{ 2 }$$

with bare dispersions ξ⁽⁰⁾_lk = k²_x/2m^l_x + k²_y/2m^l_y − μ_l and c^†_lks creates an electron in band l = (1,2,3) ≡ (d_xz, d_yz, d_xy) with momentum k, and spin s. We use mass parameters from [15]: the light masses m = m¹_x = m²_y = m³_x = m³_y = 0.7m_e and the heavy masses M = m¹_y = m²_x = 15m_e with m_e the electron mass. The chemical potentials are related by μ = μ₁ = μ₂ = μ₃ + Δ and we use in the following Δ = 50 meV. In the Zeeman term (the last term in equation (1)), $\vec {S}=\sum _{l,{\mathbf {k}},s,s'}c^\dagger _{l{\mathbf {k}} s}\vec {\sigma }_{ss'}c^{\phantom {\dag }}_{l{\mathbf {k}} s'}$ is the total spin with $\vec {\sigma }$ being the Pauli matrices. Finally, $\mathcal {H}^{\mathrm {R}}_{\mathrm { soc}}$ is the Rashba SOC at the interface due to the absence of the in-plane mirror symmetry, which can phenomenologically be introduced as a relativistic effect due to an electric field $\vec {E}$ in z direction: the spin of an electron moving with velocity $\vec {v}$ couples to an effective magnetic field $(\vec {E}/c \times \vec {v})$. Since the velocity v_l = ∂ξ⁽⁰⁾_lk/∂k for an electron in band l is k-dependent, the Rashba SOC is

$$\begin{equation} \mathcal{H}_{\mathrm{soc}}^{\mathrm{R}} = \alpha\! \sum_{l,{\mathbf{k}},s, s'} \vec{g}_{l{\mathbf{k}}}\cdot(c^{\dagger}_{l{\mathbf{k}} s}\vec{\sigma}_{ss'}c^{\phantom{\dag}}_{l{\mathbf{k}} s'}), \end{equation} \tag{ 3 }$$

where $\vec {g}_{l{\mathbf {k}}} = (v_{l,y}, -v_{l,x}, 0)$ and the overall scale α ≈ 10⁻¹¹ eV m [16]. The Rashba term $\mathcal {H}^{\mathrm {R}}_{\mathrm { soc}}$ changes the zero-field bandstructure by splitting the spin degeneracy of the individual bands, as shown in figure 1(b). Notice how different bands have different band-edge configurations: a ring of lowest energy momenta for the 2D d_xy band and a saddle point for the two quasi-one-dimensional (quasi-1D) bands stemming from the d_xz and d_yz orbitals. This difference is due to the isotropic (anisotropic) momentum dependence of the Rashba coupling $\vec {g}_{l{\mathbf {k}}}$ for the 2D (quasi-1D) bands and leads to different types of divergences in the density of states ρ(ε) at the respective band edges: A $1/\sqrt {\varepsilon }$ divergence at the bottom of the d_xy band and a logarithmic divergence near the bottom of the d_xz, d_yz bands. In the regime of low electron densities, this system can thus have instabilities to broken-symmetry phases even for weak interactions.

Now we turn to the impact of $\mathcal {H}^{\mathrm {R}}_{\mathrm { soc}}$ on the in-field ($\vec {H}\neq 0$) bandstructure, which leads to one of our main results: the anisotropy in the bare uniform susceptibility near band edges. In order to calculate the bare susceptibility, we first diagonalize the Hamiltonian (1) for each momentum k to obtain the in-field spectrum $\xi _{{\mathbf {k}}\nu }(\vec {H})$, where ν = 1,...,6. Then, the (diagonal) susceptibility is given through $\chi _{i} = \left .\frac {\partial ^2\omega }{\partial H_i^2}\right |_{\vec {H}=0}$, with ω = Ω/N the grand potential per lattice site,

$$\begin{eqnarray} &&\fl \chi_{i} = \frac1N\sum_{\nu,{\mathbf{k}}}\left\{\frac{1}{4T\cosh[\xi_{{\mathbf{k}}\nu}(\vec{H})/(2T)]^2}\left[\frac{\partial \xi_{{\mathbf{k}}\nu}(\vec{H})}{\partial H_i}\right]^2 - n_{\mathrm{F}}[\xi_{{\mathbf{k}}\nu}(\vec{H})]\left[\frac{\partial^2\xi_{{\mathbf{k}}\nu}(\vec{H})}{(\partial H_i)^2}\right]\right\}|_{\vec{H} =0}, \end{eqnarray}\noindent \tag{ 4 }$$

where T is the temperature and n_F(ξ) is the Fermi distribution function. In the absence of the Rashba term, the Hamiltonian (1) is diagonal for the spin-quantization direction parallel to $\vec {H}$ and $\xi _{{\mathbf {k}}\nu }(\vec {H})$ is linear in $\vec {H}$. Hence, only the first term in equation (4) contributes to the susceptibility: the usual Pauli susceptibility. Since the Rashba term $\mathcal {H}^{\mathrm {R}}_{\mathrm { soc}}$ does not commute with the Zeeman term, $\xi _{{\mathbf {k}}\nu }(\vec {H})$ is not linear in $\vec {H}$ once $\mathcal {H}^{\mathrm {R}}_{\mathrm { soc}}$ is present and the second term of equation (4), the so-called van Vleck susceptibility, becomes non-zero. The direction-dependent balance between the two contributions determines a possible anisotropy χ_z ≠ χ_x.

As the Hamiltonian in equation (1) is block diagonal in the orbital basis, the total bare susceptibility is the sum of contributions χ_i(l) from each orbital l. We start with the contribution from the quasi-1D orbitals. The two d_xz bands have dispersions

$$\begin{equation} \xi^{xz}_{{\mathbf{k}} \pm}(\vec{H}) = \frac{k_x^2}{2m} + \frac{k_y^2}{2M} - \mu \pm|(\alpha \vec{g}_{xz,{\mathbf{k}}} - \mu_{\mathrm{B}}\vec{H})|. \end{equation} \tag{ 5 }$$

They contribute to the total susceptibility through

$$\begin{equation} \chi^{\mathrm{P}}_{i}(xz) = \mu_{\mathrm{B}}^2\sum_{{\mathbf{k}},\pm} (\hat{g}_{xz,{\mathbf{k}}}^i)^2\frac{1}{4T\cosh[\xi^{xz}_{{\mathbf{k}}\pm}(0)/(2T)]^2} \end{equation} \tag{ 6 }$$

and

$$\begin{equation} \chi^{\mathrm{vV}}_{i}(xz) = \mu_{\mathrm{B}}^{2}\sum_{{\mathbf{k}}}[1-(\hat{g}_{xz,{\mathbf{k}}}^i)^2] \frac{n_{\mathrm{F}}[\xi^{xz}_{{\mathbf{k}}-}(0)]-n_{\mathrm{F}}[\xi^{xz}_{{\mathbf{k}}+}(0)]}{|\vec{g}_{xz,{\mathbf{k}}}|} \end{equation} \tag{ 7 }$$

with $\hat {g}_{xz,{\mathbf {k}}}$ the unit vector along $\vec {g}_{xz,{\mathbf {k}}}$. For T → 0, we substitute $\tilde {k}_y = (M/m) k_y$ and change the sums in equations (6) and (7) into (cylindrical) integrals. For μ > 0, we obtain

$$\begin{equation} \chi^{\mathrm{P}}_{i}(xz) =\frac{\mu_{\mathrm{B}}^2 M}{2\pi^2}\int \mathrm{d}\phi \frac{(\tilde{g}_{xz, \phi}^i)^2}{\cos^2\phi + \frac M m \sin^2\phi} \end{equation} \tag{ 8 }$$

and

$$\begin{equation} \chi^{\mathrm{vV}}_{i}(xz) =\frac{\mu_{\mathrm{B}}^2 M}{2\pi^2} \int \mathrm{d}\phi\frac{1-(\tilde{g}_{xz, \phi}^i)^2}{\cos^2\phi + \frac M m \sin^2\phi} \end{equation} \tag{ 9 }$$

with $\tilde {g}_{xz, \phi }^{i} = (\sin \phi, -\!\cos \phi , 0)$ and ϕ is the angle relative to the crystalline x-axis. Clearly, the total d_xz contribution to the susceptibility $\chi _i(xz)= \mu _{\mathrm {B}}^2\sqrt {mM}/\pi$ is independent of the field direction i or the Rashba SOC strength α.

On the other hand, near the band edge of the 1D bands, i.e. −α²/2m < μ < 0, the susceptibilities in equations (6) and (7) yield

$$\begin{eqnarray} &&\fl \chi^{\mathrm{P}}_{i}(xz) =\frac{\mu_{\mathrm{B}}^2 M}{2\pi^2}\int \mathrm{d}\phi \frac{(\tilde{g}_{xz, \phi}^i)^2}{\cos^2\phi + \frac M m \sin^2\phi} \left[1 + 2\,\mu m\left(\cos^2\phi + \frac M m \sin^2\phi\right)\!\Bigg/\alpha^2\right]^{-1/2} \end{eqnarray}\noindent \tag{ 10 }$$

and

$$\begin{eqnarray} &&\fl \chi^{\mathrm{vV}}_{i}(xz) =\frac{\mu_{\mathrm{B}}^2 M}{2\pi^2} \int \mathrm{d}\phi\frac{1-(\tilde{g}_{xz, \phi}^i)^2}{\cos^2\phi + \frac M m \sin^2\phi}\left[1 + 2\,\mu m\left(\cos^2\phi + \frac M m \sin^2\phi\right)\Bigg/\alpha^2\right]^{1/2}. \end{eqnarray}\noindent \tag{ 11 }$$

The total contribution to the susceptibility is thus anisotropic at the band edge.

For the total bare susceptibility, we also need to consider the 2D orbital d_xy contribution. We can read off χ_i(xy) from the quasi-1D contribution χ_i(xz) discussed above by setting M = m and shifting the band edge by Δ. Therefore, χ_i(xy) = χ^xy₀ is isotropic for μ > −Δ. When contributions from all the orbitals are combined, the total susceptibility shows two regions in the chemical potential where χ_x = χ_y > χ_z: near the band edge of the 2D bands and that of the quasi-1D bands. Figure 2(a) shows the total susceptibility for the three bands near the quasi-1D band edge as a function of chemical potential, where we shaded the anisotropic region.

Zoom In Zoom Out Reset image size

Figure 2(b) summarizes two ways atomic SOC impacts the phase diagram: (i) it widens the region of anisotropic susceptibility (gray region), (ii) it causes anisotropy in the spin direction by aligning the spin with the largest angular momentum when the susceptibility is isotropic (white region). The first effect results from adding the atomic SOC,

$$\begin{equation} \mathcal{H}_{\mathrm{SOC}}=\mathrm{i} \frac{\alpha^{\mathrm{at}}}{2}\sum_{lmn}\epsilon_{lmn}\sum_{{\mathbf{k}}, s, s'}c^{\dagger}_{l{\mathbf{k}} s}c^{\phantom{\dag}}_{m{\mathbf{k}} s'}\sigma^n_{ss'} \end{equation} \tag{ 12 }$$

to the Hamiltonian (1) and evaluating the susceptibility (4) numerically. For the second effect, the atomic SOC is treated as a perturbation in the magnetic phase [12]. This again predicts in-plane magnetization below the gray region where only the 2D band is occupied, but a switch to out-of-plane magnetization above this region.

In order to compare our phase diagram with experiments, we need to translate the chemical potential to a gate voltage. While such a translation is non-trivial, Hall measurements under gate-voltage sweeps can offer hints as to where the as-grown samples lie. Joshua et al [13] showed that the interface acquires heavy carriers when under a gate voltage, which indicates that the as-grown samples are near the quasi-1D band edge. Though the gray region is narrow in figure 2(b), the density changes by a factor of ≈2 in this range. This pushes the density upper bound for in-plane magnetization substantially, consistent with experiments [17].

In the regime of low electron densities, the van Hove singularity can promote broken symmetry phases in the presence of suitable interactions. While we focused so far on a magnetic instability, which could for instance be driven by repulsive intra-orbital interactions, we now turn to orbital-ordering possibilities. Even though the three-fold degeneracy of the t_2g orbitals is already broken due to the interface symmetry, the d_xz and d_yz orbitals remain degenerate. A spontaneous orbital symmetry breaking described by a non-vanishing order parameter η ≡ n₁ − n₂, with n_1/2 the occupation of the d_xz/d_yz orbital, can be driven by inter-orbital repulsive interactions [18]. We first analyze the tendency toward such an instability by studying the (bare) nematic susceptibility. For this purpose, we introduce the field H_η conjugate to η, which enters the Hamiltonian (1) through μ_1/2 = μ ± H_η, i.e. it acts through an opposite shift of the chemical potential for the two orbitals. The nematic susceptibility then yields

$$\begin{equation} \chi_{\eta} = \left.\frac{\partial^2\omega}{\partial H_{\eta}^2}\right|_{H_\eta = 0}=\frac{\partial^2 \omega}{\partial\mu_1^2} + \frac{\partial^2\omega}{\partial \mu_2^2} \end{equation} \tag{ 13 }$$

and is given by the contribution of the d_xz and d_yz orbitals to the density of states. Hence, the van Hove singularity near the band edge of the quasi-1D bands (see figure 1(b)) allows for a nematic order for sufficiently strong inter-orbital interaction [18].

Next, we consider the coexistence of nematic and magnetic order. While a system with C₄ symmetry does not allow for coupling of η and $|\vec {M}|$ directly, there is an allowed tri-linear coupling between η and an in-plane anisotropy M²_x − M²_y, as both acquire a factor of −1 under C₄ rotation. Specifically, a coupling of the form λ(M²_x − M²_y)η enters the free energy. Given the Hamiltonian (1), we can explicitly calculate the coupling constant λ. It is given by the generalized susceptibility

$$\begin{equation} \lambda \propto \chi_{\eta xx} = \left(\frac{\partial^3 \omega }{\partial H_\eta\partial H_x^2}\right)_{|\vec{H}|=0}, \end{equation}\noindent \tag{ 14 }$$

which measures the change in the spin susceptibility χ_x upon shifting the chemical potential of the two quasi-1D bands against each other. Figure 3 shows the result with and without atomic SOC in the presence of Rashba SOC. (We find the coupling to be negligible without Rashba SOC.) The coupling becomes substantial in the above-identified density range with the in-plane preference in the bare spin susceptibility. This is due to the inequivalence between the two quasi-1D band contributions to χ_x shown as dotted and dash/dotted lines in figure 2(a), more specifically the difference in the slope of χ_x(xz) and χ_x(yz). The sign of λ determines the relative sign between η and (M²_x − M²_y) and whether the majority orbitals will be along or perpendicular to the magnetization axis. Notice how our result shown in figure 3 predicts a change of the sign of λ over the chemical potential range exhibiting the magnetization anisotropy.

Zoom In Zoom Out Reset image size

Given the observed magnetism, we finally investigate the effects of coupling between magnetic order and nematic fluctuations associated with the nearby nematic phase. In particular, we may ask on the one hand how the proximity to a nematic instability influences the in-plane magnetization, and on the other hand, how a magnetization effects the orbital ordering. Assuming an XY ferromagnet in the absence of a coupling to nematic order, the Landau free energy becomes

$$\begin{equation} \fl f(M_x, M_y, \eta) = f_0 + \frac{a(T)}2 |\vec{M}|^2 + \frac b4 |\vec{M}|^4 + \lambda(M_x^2 - M_y^2)\eta + \frac{a_{\eta}}2 \eta^2+\cdots, \end{equation} \tag{ 15 }$$

where M_x (M_y) is the magnetization along the crystalline x (y)-axis. Since we assume no independent instability towards an orbital-ordered nematic, a_η > 0⁴. Integrating out η by minimizing the free energy with

$$\begin{equation} \eta = -\frac{\lambda}{a_{\eta}}(M_x^2 - M_y^2), \end{equation} \tag{ 16 }$$

the free energy for the magnetization becomes

$$\begin{equation} f(M_x, M_y) = f_0 + \frac{a(T)}2 |\vec{M}|^2 + \frac b4|\vec{M}|^4 -\frac{\lambda^2}{a_{\eta}}(M_x^2-M_y^2)^2. \end{equation} \tag{ 17 }$$

For finite λ, the coupling to orbital order therefore locks the magnetization along one of the crystal axes, leading to an additional in-plane anisotropy.

For the orbital ordering, the magnetization anisotropy acts as a driving field. For a_η very small, i.e. the system close to an instability, we should include additional terms to the free energy for η,

$$\begin{equation} f(\eta;\vec{M}) = \frac{a_{\eta}}2\eta^2 + \frac{b_{\eta}}4\eta^4 + \frac{c_{\eta}}6 \eta^6 +\lambda(M^2_x - M_y^2)\eta+\cdots. \end{equation} \tag{ 18 }$$

For b_η < 0, the system undergoes either a metanematic crossover, or a first-order transition at a critical magnetization. This is in analogy to metamagnetic transitions observed in systems close to ferromagnetism [19] and could be driven indirectly here by an applied magnetic field.

We have shown that the combination of Rashba SOC and atomic SOC leads to an electron-density dependent magnetization anisotropy in LAO/STO interfaces. While experiments so far appear to lie in the in-plane-magnetization region, we predict a switch to out-of-plane magnetization at sufficiently high gate voltages. We have identified a regime near the band edge with anisotropic susceptibility as a non-trivial effect of the Rashba SOC in a low carrier density system. The high density of states near the band edge in principle also allows for a spontaneous orbital ordering and we predict in this regime an enhanced coupling between the magnetization direction and this kind of nematic order. This coupling locks the in-plane magnetization direction to be along one of the crystal axes and promotes Ising nematicity.

Next, we comment on the issue of heterogeneity detected in [8]. The observed heterogeneity is likely driven by both extrinsic and intrinsic effects. For instance, a recent study showed that strong Rashba SOC can promote phase separation [20]. The proposed magnetization–nematicity coupling has important consequences in both extrinsic and intrinsic fronts. On the one hand, oxygen vacancies and other spatial inhomogeneities act as a random field for the Ising nematic and in turn cause a distribution of moment directions. On the other hand, the reduction of the magnetic order parameter symmetry due to the coupling changes the type of magnetic textures and their energetics. Furthermore, we expect the coupling to cause non-trivial in-plane anisotropy in the magneto transport [21]. Moreover, the sign change in the coupling λ could be observed through the rotation in the dominant direction upon gate voltage sweep. Finally, we note that the range in the density with magnetization–nematicity near the band edge of the quasi-1D bands has also been shown to exhibit critical scaling in recent Hall measurement [13].

We are grateful to J Bert, H Hwang, B Kalisky and K Moler for useful discussions. MHF and E-AK acknowledge support from NSF grant no. DMR-0955822 and from NSF grant no. DMR-1120296 to the Cornell Center for Materials Research. SR acknowledges support from the LDRD program at SLAC and the Alfred P Sloan Research fellowship.