Magnetic and nematic orders of the 2DEG at oxide (111) surfaces and interfaces

TL;DR

This paper investigates the magnetic and nematic orders in two-dimensional electron gases at oxide (111) surfaces, highlighting how multiorbital interactions can lead to various symmetry-breaking phenomena relevant to observed nematicity and superconductivity.

Contribution

It introduces an effective field theory capturing the interplay of magnetism and nematic charge order in (111) 2DEGs, considering multiorbital interactions and their effects.

Findings

Multiorbital interactions induce diverse magnetic and orbital orders.

Broken symmetries may explain observed nematicity in oxide 2DEGs.

Implications for transport and superconductivity are discussed.

Abstract

Recent experiments have explored two-dimensional electron gases (2DEGs) at oxide (111) surfaces and interfaces, finding evidence for hexagonal symmetry breaking in SrTiO$_3$ at low temperature. We discuss many-body instabilities of such (111) 2DEGs, incorporating multiorbital interactions in the $t_{2g}$ manifold which can induce diverse magnetic and orbital orders. Such broken symmetries may partly account for the observed nematicity, cooperating or competing with phonon mechanisms. We present an effective field theory for the interplay of magnetism and nematic charge order, and discuss implications of the nematicity for transport and superconductivity in (111) 2DEGs.

Equations

$$H_0=\sum _{\mathbf{k}\ell {\ell }^{\prime }\sigma }{c}_{\ell \sigma }^{†}(\mathbf{k})h_{\ell {\ell }^{\prime }}^{}(\mathbf{k})c_{{\ell }^{\prime }\sigma }^{}(\mathbf{k})$$

$$\underline{h}(\mathbf{k})=\left(\begin{array}{ccc}{ϵ}_{\mathbf{k}}^c+{\eta }_{\mathbf{k}}^{ab}& {\gamma }_{\mathbf{k}}^a& {\gamma }_{\mathbf{k}}^b\\ {\gamma }_{\mathbf{k}}^a& {ϵ}_{\mathbf{k}}^b+{\eta }_{\mathbf{k}}^{ca}& {\gamma }_{\mathbf{k}}^c\\ {\gamma }_{\mathbf{k}}^b& {\gamma }_{\mathbf{k}}^c& {ϵ}_{\mathbf{k}}^a+{\eta }_{\mathbf{k}}^{bc}\end{array}\right),$$

$$H_{\mathrm{int}}$$

$${\left({\underline{\chi }}_0(\mathbf{q},\Omega )\right)}_{{\ell }_1{\ell }_2;{\ell }_3{\ell }_4}=\frac{1}{N}\sum _{ij\sigma {\sigma }^{\prime }}{\int }_0^{\beta }\mathrm{d}\tau {\mathrm{e}}^{\mathrm{i}\mathbf{q}\cdot ({\mathbf{r}}_i-{\mathbf{r}}_j)-\mathrm{i}\Omega \tau }$$

$$\times ⟨c_{i{\ell }_1\sigma }^{†}(\tau )c_{i{\ell }_2\sigma }(\tau )c_{j{\ell }_3{\sigma }^{\prime }}^{†}(0)c_{j{\ell }_4{\sigma }^{\prime }}(0)⟩.$$

$$\displaystyle\begin{array}[]{lcl}(U^{c})_{\ell\ell;\ell\ell}=-U,&&(U^{c})_{\ell\ell^{\prime};\ell\ell^{\prime}}=-2V,\\ (U^{c})_{\ell\ell;\ell^{\prime}\ell^{\prime}}=V-\frac{3}{2}J,&&(U^{c})_{\ell\ell^{\prime};\ell^{\prime}\ell}=-2J^{\prime},\end{array}$$

$$\displaystyle\begin{array}[]{lcl}(U^{s})_{\ell\ell;\ell\ell}=U,&&(U^{s})_{\ell\ell^{\prime};\ell\ell^{\prime}}=J,\\ (U^{s})_{\ell\ell;\ell^{\prime}\ell^{\prime}}=\frac{1}{2}J-V,&&(U^{s})_{\ell\ell^{\prime};\ell^{\prime}\ell}=2J^{\prime},\end{array}$$

$${\mathcal{L}}_{\psi }$$

$${\mathcal{L}}_{\phi }$$

$${\mathcal{L}}_{\psi \phi }$$

$${\left({\chi }_0(\mathbf{q},\Omega )\right)}_{{\ell }_1{\ell }_2;{\ell }_3{\ell }_4}=\frac{1}{N\beta }\sum _{\mathbf{k},{\omega }_n}{G}_{{\ell }_2{\ell }_3}(\mathbf{k},{\omega }_n)G_{{\ell }_4{\ell }_1}(\mathbf{k}+\mathbf{q},{\omega }_n+\Omega )$$

$$G_{\ell {\ell }^{\prime }}(\mathbf{k},{\omega }_n)=\sum _{m=1}^3\frac{u_{\ell m}(\mathbf{k})u_{{\ell }^{\prime }m}^{\ast }(\mathbf{k})}{\mathrm{i}{\omega }_n-{\epsilon }_m(\mathbf{k})+\mu }.$$

$$h^R(\mathbf{k})=2\mathrm{i}{\lambda }_{OR}\left(\begin{array}{ccc}0& \sin (k_b)-\sin (k_c)& \sin (k_a)+\sin (k_c)\\ \sin (k_c)-\sin (k_b)& 0& \sin (k_b)-\sin (k_a)\\ -\sin (k_a)-\sin (k_c)& \sin (k_a)-\sin (k_b)& 0\end{array}\right)$$

$$\begin{array}{rl}{\mathcal{F}}_{\mathrm{var}}=& \sum _{i\ell {\ell }^{\prime }}\left\{⟨n_{i,\ell }⟩⟨n_{i,{\ell }^{\prime }}⟩\left[\frac{U}{4}{\delta }_{\ell {\ell }^{\prime }}+\frac{1}{2}\left(U-\frac{5}{2}J\right)(1-{\delta }_{\ell {\ell }^{\prime }})\right]-⟨{\mathbf{S}}_{i,\ell }⟩\cdot ⟨{\mathbf{S}}_{i,{\ell }^{\prime }}⟩\left[U{\delta }_{\ell {\ell }^{\prime }}+J(1-{\delta }_{\ell {\ell }^{\prime }})\right]\right\}\\ & +{\int }_{FBZ}({\mathrm{d}}^2\mathbf{k})\mathrm{Tr}({\underline{u}}^{†}(\mathbf{k})\underline{h}(\mathbf{k})\underline{u}(\mathbf{k})\underline{n_F}(\mathbf{k}))-TS\end{array}$$

$$H_{\mathrm{spiral}}=\sum _{\mathbf{k}}{\Psi }^{†}(\mathbf{k})\mathcal{H}(\mathbf{k})\Psi (\mathbf{k})$$

$$\mathcal{H}(\mathbf{k})=\left(\begin{array}{cccccc}{\epsilon }_{\mathbf{k}}^{yz}& {\gamma }_{\mathbf{k}}^{yz,zx}& {\gamma }_{\mathbf{k}}^{yz,xy}& {\zeta }_{yz}& 0& 0\\ {\gamma }_{\mathbf{k}}^{yz,zx\ast }& {\epsilon }_{\mathbf{k}}^{zx}& {\gamma }_{\mathbf{k}}^{zx,xy}& 0& {\zeta }_{zx}& 0\\ {\gamma }_{\mathbf{k}}^{yz,xy\ast }& {\gamma }_{\mathbf{k}}^{zx,xy\ast }& {\epsilon }_{\mathbf{k}}^{xy}& 0& 0& {\zeta }_{xy}\\ {\zeta }_{yz}& 0& 0& {\epsilon }_{\mathbf{k}+\mathbf{q}}^{yz}& {\gamma }_{\mathbf{k}+\mathbf{q}}^{yz,zx}& {\gamma }_{\mathbf{k}+\mathbf{q}}^{yz,xy}\\ 0& {\zeta }_{zx}& 0& {\gamma }_{\mathbf{k}+\mathbf{q}}^{yz,zx\ast }& {\epsilon }_{\mathbf{k}+\mathbf{q}}^{zx}& {\gamma }_{\mathbf{k}+\mathbf{q}}^{zx,xy}\\ 0& 0& {\zeta }_{xy}& {\gamma }_{\mathbf{k}+\mathbf{q}}^{yz,xy\ast }& {\gamma }_{\mathbf{k}+\mathbf{q}}^{zx,xy\ast }& {\epsilon }_{\mathbf{k}+\mathbf{q}}^{xy}\end{array}\right)$$

$${\epsilon }_{\mathbf{k}}^{\ell }$$

$${\zeta }_{\ell }$$

$${\mathbf{b}}_{i\ell }={\mathfrak{b}}_0+{\mathfrak{b}}_1(-1{)}^{m+n}{\delta }_{\ell ,yz}+{\mathfrak{b}}_2(-1{)}^m{\delta }_{\ell ,zx}+{\mathfrak{b}}_3(-1{)}^n{\delta }_{\ell ,xy}$$

$$\hat{\sigma}_{a}=\left(\begin{array}[]{cc}\sigma_{1}&0\\ 0&\sigma_{2}\end{array}\right),\quad\hat{\sigma}_{b}=\left(\begin{array}[]{cc}\frac{1}{4}\sigma_{1}+\frac{3}{4}\sigma_{2}&\frac{\sqrt{3}}{4}(\sigma_{1}-\sigma_{2})\\ \frac{\sqrt{3}}{4}(\sigma_{1}-\sigma_{2})&\frac{3}{4}\sigma_{1}+\frac{1}{4}\sigma_{2}\end{array}\right),\quad\hat{\sigma}_{c}=\left(\begin{array}[]{cc}\frac{1}{4}\sigma_{1}+\frac{3}{4}\sigma_{2}&-\frac{\sqrt{3}}{4}(\sigma_{1}-\sigma_{2})\\ -\frac{\sqrt{3}}{4}(\sigma_{1}-\sigma_{2})&\frac{3}{4}\sigma_{1}+\frac{1}{4}\sigma_{2}\end{array}\right),$$

$$\hat{\sigma}\approx\frac{3}{2}(\sigma_{1}+\sigma_{2})\left(\begin{array}[]{cc}1&0\\ 0&1\end{array}\right)+\frac{1}{2}(\sigma_{1}-\sigma_{2})\left(\begin{array}[]{cc}{\rm Re}\,\psi_{n}&{\rm Im}\,\psi_{n}\\ {\rm Im}\,\psi_{n}&-{\rm Re}\,\psi_{n}\end{array}\right).$$

Full text

Introduction.— Transition metal oxide heterostructures and interfaces can realize exotic low-dimensional electronic phases and allow for engineering oxide-based devices Chakhalian et al. (2014). Extensive research Ohtomo and Hwang (2004); Thiel et al. (2006); Nakagawa et al. (2006); Takizawa et al. (2006); Reyren et al. (2007); Caviglia et al. (2008, 2010); Ariando et al. (2011); Bert et al. (2011); Li et al. (2011); Khalsa and MacDonald (2012); Mehta et al. (2012); Zhong et al. (2013); Fischer et al. (2013a); Kim et al. (2013); Banerjee et al. (2013); Chen and Balents (2013); Park and Millis (2013); Ruhman et al. (2014); Diez et al. (2015); Nandy et al. (2016); Tolsma et al. (2016); Atkinson et al. (2017) has focused on the two-dimensional electron gas (2DEG) at the (001) LaAlO3-SrTiO3 (LAO-STO) interface induced by a combination of the polar catastrophe and oxygen vacancies. This 2DEG shows evidence of correlated magnetism in torque magnetometry and scanning SQUID measurements Bert et al. (2011); Li et al. (2011). In addition, it exhibits superconductivity (SC) Reyren et al. (2007); Dikin et al. (2011) which may be tied to that of doped bulk STO, though the interface might harbor modulated SC pairing Michaeli et al. (2012) or Majorana modes Fidkowski et al. (2013).

Recently, various groups have started to probe 2DEGs at oxide (111) surfaces and interfaces, for instance induced by photon Rödel et al. (2014) or ion McKeown Walker et al. (2014); Miao et al. (2016) irradiation at the (111) STO surface, as well as that at the (111) LAO-STO interface Rout et al. (2017); Davis et al. (2017a, b); Monteiro et al. (2017); Davis et al. (2017c). Part of this interest stems from proposals for realizing topologically nontrivial phases along this growth direction Xiao et al. (2011); Rüegg and Fiete (2011); Cook and Paramekanti (2014); Okamoto et al. (2014); Hu et al. (2015); Baidya et al. (2016); Si et al. (2016); Kim and Kee (2017). The [111] growth direction is polar for STO due to alternating Ti4+ and (SrO${}_{3})^{4-}$ layers, and the internal electric fields could lead to stronger confinement Baidya et al. (2015) of the (111) 2DEG, potentially enhancing correlation effects relative to (001) 2DEGs. Angle resolved photoemission spectroscopy (ARPES) on the (111) STO surface reveals a Fermi surface (FS) composed of all three $t_{2g}$ orbitals, which appears to preserve the expected hexagonal symmetry Rödel et al. (2014); McKeown Walker et al. (2014). However, very recent experiments have discovered, via measurements of magnetotransport Miao et al. (2016); Rout et al. (2017); Davis et al. (2017a, b) and the resistive transition into the SC state Davis et al. (2017c), that this (111) 2DEG exhibits an anisotropy which sets in at low temperatures, spontaneously breaking the hexagonal symmetry. While a weak resistive anisotropy may arise from the ${\sim 100}\;\rm K$ pseudo-cubic to pseudo-tetragonal transition of bulk STO Aharony et al. (1977), the onset temperature seen in these experiments is much lower, ${T_{\rm ani}\sim 4\text{-}30\;\rm K}$ depending on the sample and the electron density. For bulk STO, it is known that the transition into the pseudo-tetragonal phase is sensitive to stress along the [111] direction Müller et al. (1991), and proceeds via an intermediate trigonal phase; it remains to be tested if the lower symmetry at the (111) surface or interface leads to a low temperature surface phonon instability.

In light of these developments, it is in any case also important to consider the impact of electron-electron interactions on (111) 2DEGs, in order to (i) study possible interaction induced many-body instabilities, and (ii) ask if there are electronic mechanisms for the observed anisotropies of the (111) 2DEG which may cooperate or compete with phonon instabilities. Such an interplay has been actively investigated in the iron pnictide superconductors (see Ref. Fernandes et al., 2014 for a review).

Motivated by these questions, we examine a model for $t_{2g}$ electronic states of the (111) STO surface 2DEG, which is consistent with the ARPES measurements, and study its instabilities driven by multiorbital interactions. Our main findings, summarized in Figs. 1 and 2, based on a combination of random phase approximation (RPA) calculations supplemented by mean field theory, is that there is a range of densities over which this 2DEG is unstable to ferromagnetic (FM) or antiferromagnetic (AF) order, accompanied by ferro-orbital order. Even if thermal fluctuations melt such magnetic orders in 2D, the orbital order and the fluctuating magnetism are expected to survive to higher temperatures, leading to a nematic fluid Oganesyan et al. (2001); Khavkine et al. (2004); Kee and Kim (2005); Fradkin et al. (2007); Fang et al. (2008); Fradkin et al. (2010) which breaks hexagonal symmetry. We present a Landau theory of this nematic, and discuss implications for transport measurements and superconductivity. Such nematicity induced by orbital or spin order has been previously considered for the (001) 2DEG Zhou et al. (2015); Fischer et al. (2013b); Li et al. (2014); Tolsma et al. (2016). Our results should be broadly applicable to a wide class of oxide (111) 2DEGs.

Model. — We begin with a tight-binding model of Ti $t_{2g}$-orbitals on a 2D triangular lattice which captures the FS seen in ARPES Rödel et al. (2014); McKeown Walker et al. (2014) for the (111) 2DEG at the STO surface:

[TABLE]

with $\ell\equiv yz,zx,xy$, and

[TABLE]

where $\epsilon^{\alpha}_{{\bf{k}}}\!=\!-2t\cos k_{\alpha}$ and $\eta^{\alpha\beta}_{{\bf{k}}}\!=\!-2t_{\perp}(\cos k_{\alpha}\!+\!\cos k_{\beta})$ determine the intraorbital dispersion which leads to elliptical FSs, while $\gamma^{\alpha}_{{\bf{k}}}=-2t^{\prime}\cos k_{\alpha}$ captures weak interorbital hopping. Here, we have defined $k_{\alpha}={{\bf{k}}}\cdot\hat{\alpha}$ ($\alpha=a,b,c$), with $\hat{a}=\hat{x}$, $\hat{b}=\hat{x}/2+\hat{y}\sqrt{3}/2$, and $\hat{c}=-\hat{x}/2+\hat{y}\sqrt{3}/2$. We work in units where the triangular lattice constant $d\approx 5.66\;\rm\AA$ is set to unity. The ARPES data Rödel et al. (2014); McKeown Walker et al. (2014) can be reasonably fit by choosing ${t=320\;\text{meV}}$ and $t_{\perp}=0.04t$, and an electron density of $\bar{n}=0.3$ electrons per site, corresponding to $10^{14}$ cm*-2*; we therefore study a range of densities around this value. The interorbital terms appear to be small; for concreteness, we set $t^{\prime}=-0.04t$. The resulting FSs are shown overlaid on the paramagnetic phases in Fig. 2. The real 2DEG wave functions will be spread over a few layers, so $H_{0}$ should only be viewed as the simplest 2D tight-binding parameterization of the ARPES data. We omit spin-orbit coupling (SOC), but comment on its effects later. The local multiorbital interactions are

[TABLE]

where $i$ denotes the site and $\ell$ the orbital. Below, we fix $V\!=\!(U-5J/2)$ and $J^{\prime}\!=\!J$ as appropriate for $t_{2g}$ orbitals, and explore broken symmetry states driven by varying the interactions $J/t,U/t$. These interactions should be scaled down compared to atomic values by the number of layers over which the 2DEG is spread.

RPA analysis. — To identify the leading weak-coupling instabilities we use an unbiased multi-orbital RPA approach Graser et al. (2009), with the matrix response ${\underline{\chi}^{(c,s)}_{RPA}({\bf{q}},\Omega)=\underline{\chi}_{0}({\bf{q}},\Omega)(1-\underline{U}^{(c,s)}\underline{\chi}_{0}({\bf{q}},\Omega))^{-1}}$, where

[TABLE]

is the bare response function (see Supplemental Material (SM)). Here, $N$ is the number of sites, $c$ and $s$ respectively denote charge and spin responses. The non-zero interaction vertex matrices are,

[TABLE]

where $\ell\neq\ell^{\prime}$. When the largest eigenvalue of $\underline{U}^{(c,s)}\underline{\chi}_{0}({\bf{q}},\Omega=0)$ is $\lambda^{(c,s)}({\bf{q}})=1$, the response function diverges, indicating an instability towards an ordered state, with corresponding eigenvectors, $f^{(c,s)}_{\ell\ell^{\prime}}({\bf{q}})$. Figure 1 shows the largest eigenvalues, $\lambda^{c}({\bf{q}})$ and $\lambda^{s}({\bf{q}})$ along high symmetry lines in the Brillouin zone (BZ) for ${\bar{n}=0.3}$, $U/t=2$ and temperature $T/t=0.02$, demonstrating the emergence of various instabilities as we vary Hund’s coupling. For our choice of experimentally motivated parameters, the leading instabilities are nearly orbital diagonal, ${f^{(c,s)}_{\ell\ell^{\prime}}({\bf{q}})\approx f^{(c,s)}_{\ell}({\bf{q}})\delta_{\ell\ell^{\prime}}}$.

When $J=0$, the Hubbard interaction $U$ drives a leading instability in the spin channel, with ${\bf{q}}$ generically incommensurate (close to the $M$ points of the BZ for the range of densities investigated), with all the weight on a single orbital. This instability indicates a tendency towards incommensurate spiral or commensurate stripe antiferromagnetic (AF) order – for each orbital in a different direction. In the charge channel, we find a subleading instability, with two degenerate eigenvalues at ${\bf{q}}=0$, indicating a tendency towards ferro-orbital order which will lead to nematicity associated with broken lattice rotational symmetry. With increasing $J$, the AF instability gives way to a ferromagnetic (FM) instability seen in the spin response at the $\Gamma$-point. At the same time, the ferro-orbital response is strongly suppressed. Below, we use mean field theory (MFT) in order to further characterize the broken symmetry phases.

Mean field theory. — We study the phase diagram of our model, Eqns. (1)-(3), using a momentum space MFT within a single-${{\bf{Q}}}$ spiral ansatz with a spatially uniform but orbital-dependent density. This is captured by a mean field Hamiltonian, ${H_{\rm var}=H_{0}-\sum_{i\ell}(\phi_{\ell}+\mu)n_{i\ell}-\sum_{i\ell}{\bf b}_{\ell}\cdot{\bf S}_{i\ell}{\rm e}^{{\mathrm{i}}{{\bf{Q}}}\cdot{{\bf{r}}}_{i}}}$, which we use to generate variational ground states $\ket{\psi_{\rm var}}$ at the desired charge density by tuning $\mu$. The fields $\phi_{\ell},{\bf{b}}_{\ell}$, and the wavevector ${{\bf{Q}}}$ are selected to minimize $\bra{\psi_{\rm var}}H_{0}+H_{\rm int}\ket{\psi_{\rm var}}$ (see SM).

Fig. 2 (a)-(c) shows the MFT phase diagrams for densities $\bar{n}=0.15,0.30,0.45$. Broadly, we find three phases consistent with RPA: (i) a paramagnetic (PM) metal where $\phi_{\ell}\!=\!0,{\bf b}_{\ell}\!=\!0$, (ii) a commensurate stripe-AF (AF-1) or incommensurate spiral (ICS) metal driven by $U$, where $\phi_{\ell}\neq 0$ and ${\bf b}_{\ell}\neq 0$, which has higher density in one of the three orbitals, and (iii) a FM metal driven by Hund’s coupling where ${\bf b}_{\ell}\!=\!{\bf b}$, and ${{\bf{Q}}}\!=\!0$, either having the same density in all orbitals (FM), or with one orbital having a lower density than the other two (i.e., a nematic FM: ${\cal N}$FM). Interestingly, we find a direct (generically first-order) transition from the PM into the AF-1 phase for all three densities, in contrast to the RPA which finds AF-1 only at a fine-tuned density.

To go beyond the single-${{\bf{Q}}}$ ansatz, we have also studied the commensurate AF by minimizing the free energy at $T\!=\!0.02t$ assuming a $2\!\times\!2$ unit cell. The result for $\bar{n}\!=\!0.3$ is shown in Fig. 2(d) (see also SM); we find reasonably good agreement with the single-${{\bf{Q}}}$ MFT, but discover small regimes where the single-${{\bf{Q}}}$ order gives way to multi-${{\bf{Q}}}$ condensates where two (AF-$2$) or three (AF-$3$) wavevectors are simultaneously present. A robust feature is the presence of simultaneous AF and ferro-orbital order in the AF-$1$ and ${\text{AF-}2}$ phases. A similar competition between single-${{\bf{Q}}}$ and multi-${{\bf{Q}}}$ phases also appears in single-orbital honeycomb and triangular lattice Hubbard models Chern et al. (2012); Nandkishore et al. (2012).

The ${\cal N}$FM, AF-$1$, AF-$2$, and ICS phases feature discrete ferro-orbital order which breaks the hexagonal lattice symmetry. Thus, even if fluctuations melt the magnetic order itself, there may be large regimes in the phase diagram where the electronic nematicity survives. Below we use Landau theory to further understand this interplay of magnetism and nematicity.

Effective field theory. — Landau theory is a powerful tool to analyze magnetic orders Banerjee et al. (2013) and study spin textures such as skyrmions which might arise at the (001) LAO-STO interface Li et al. (2014). For the (111) 2DEG, our RPA analysis suggests that the soft electronic modes include a complex nematic charge mode ${\psi_{n}=\delta\rho_{xy}+\omega\delta\rho_{yz}+\omega^{2}\delta\rho_{zx}}$ (with $\omega={\rm e}^{{\mathrm{i}}2\pi/3}$) constructed from the slowly varying orbital densities $\delta\rho_{\ell}$ at the $\Gamma$-point, and the spin modes, $\vec{\varphi}_{0}$ at the $\Gamma$-point, and $\vec{\varphi}_{\alpha}$ at the magnetic wavevectors ${{\bf{Q}}}_{\alpha}$ ($\alpha=1,2,3$), which can describe both FM and AF orders. Since the interorbital hopping is small, $\vec{\varphi}_{1,2,3}\sim\vec{\varphi}_{xy,yz,zx}$ but with weak orbital admixture. The spin modes $\vec{\varphi}_{\alpha}$ are complex for incommensurate ${{\bf{Q}}}_{\alpha}$, and real if ${{\bf{Q}}}_{\alpha}$ correspond to the commensurate $M$-points. The nematic order parameter $\psi_{n}$ transforms under anticlockwise lattice ${\pi/3\text{-rotations}}$ as $\psi_{n}\to\omega^{2}\psi_{n}$, and under reflections about the ${\hat{x}\text{-axis}}$ as $\psi^{\phantom{\dagger}}_{n}\to\psi^{*}_{n}$. Turning to the spin modes, $\vec{\varphi}_{0}$ is invariant under lattice symmetries, anticlockwise $\pi/3$-rotations leads to ${\vec{\varphi}_{1}\to\vec{\varphi}_{2},\vec{\varphi}_{2}\to\vec{\varphi}_{3},\vec{\varphi}_{3}\to\vec{\varphi}^{*}_{1}}$. Under spin rotations, $\psi_{n}$ is invariant but all spin modes undergo $SO(3)$ rotations, $\vec{\varphi}_{0,\alpha}\to\mathfrak{R}\vec{\varphi}_{0,\alpha}$. Time-reversal sends $\vec{\varphi}_{0,\alpha}\to-\vec{\varphi}_{0,\alpha}$. Armed with this, the mean field Landau free energy is ${\cal F}=\int\mathrm{d}^{2}x\;({\cal L}_{\psi}+{\cal L}_{\varphi}+{\cal L}_{\psi\varphi})$, with

[TABLE]

where we have defined a complex magnetic nematic order ${\cal S}_{n}=|\vec{\varphi}_{1}|^{2}+\omega|\vec{\varphi}_{2}|^{2}+\omega^{2}|\vec{\varphi}_{3}|^{2}$ which transforms analogous to $\psi_{n}$. In this effective field theory, the PM phase corresponds to $(r_{\psi},r_{0},r_{Q}>0)$, the FM phase corresponds to $(r_{0}<0,r_{\psi},r_{Q}>0)$, and the various AF phases correspond to $(r_{Q}<0,r_{0}>0)$. Note that we never find a ground state nematic charge order unaccompanied by spin order in the MFT. The various types of AF orders will be dictated by higher order (quartic and sixth order) terms denoted above by ellipsis. In turn, this can lead to a ‘pinning field’ for the charge nematic order via the cubic interaction $\lambda_{1}$ in ${\cal L}_{\psi\varphi}$; our mean field results indicate $w_{\psi},\lambda_{1},\lambda_{2}>0$. Below, we discuss some implications of this Landau theory, deferring its microscopic derivation to a future publication.

(a) Incommensurate spiral (ICS): For the generic incommensurate instability, the simplest spin order (which will not lead to any charge modulation) is a single-mode coplanar spiral at one ${{\bf{Q}}}_{\alpha}$, with the complex ${\vec{\varphi}_{\alpha}=\hat{\Omega}_{1}+{\mathrm{i}}\hat{\Omega}_{2}}$ and $\hat{\Omega}_{1}\cdot\hat{\Omega}_{2}=0$. This leads to ${{\cal S}_{n}\sim\omega^{\alpha-1}}$, so the cubic interaction $\lambda_{1}>0$ will pin $\psi_{n}\sim\omega^{\alpha-1}$, causing density enhancement in orbital $\alpha$.

(b) Commensurate AF: For commensurate stripe AF, ${{\bf{Q}}}_{\alpha}\equiv M_{\alpha}$, in which case $\vec{\varphi}_{\alpha}$ are real fields. This case, for which we have also carried out a real space MFT, leads to three orders. (i) AF-1 has condensation at a single $M_{\alpha}$, which corresponds to a collinear stripe order with $\vec{\varphi}_{\alpha}\sim\hat{\Omega}_{1}$. This has ${\cal S}_{n}\sim\omega^{\alpha-1}$, which pins ${\psi_{n}\sim\omega^{\alpha-1}}$, leading to a charge nematic order similar to the ICS state. (ii) AF-2 features condensation at a pair of wavevectors $M_{\alpha},M_{\beta}$, with $\langle\vec{\varphi}_{\alpha}\rangle\sim\hat{\Omega}_{1}$ and $\langle\vec{\varphi}_{\beta}\rangle\sim\hat{\Omega}_{2}$. This state can be either collinear, $\hat{\Omega}_{1}=\hat{\Omega}_{2}$, or coplanar, $\hat{\Omega}_{1}\cdot\hat{\Omega}_{2}=0$. In both case, however, ${{\cal S}_{n},\psi_{n}\sim\omega^{\alpha-1}+\omega^{\beta-1}}$, displaying charge nematic order. (iii) Finally, AF-3 is a triple-${{\bf{Q}}}$ spin crystal, similarly featuring either collinear or non-coplanar tetrahedral order of the spins. Both cases are obtained by condensation at all three $M$-points, with $\langle\vec{\varphi}_{1,2,3}\rangle\sim\hat{\Omega}_{1,2,3}$, and ${\cal S}_{n}=0$, so no charge nematic is induced. In the collinear case $\hat{\Omega}_{1}=\hat{\Omega}_{2}=\hat{\Omega}_{3}$, while in the non-coplanar spin order $\hat{\Omega}_{1}\cdot\hat{\Omega}_{2}=0$, and $\hat{\Omega}_{3}=\pm\hat{\Omega}_{1}\times\hat{\Omega}_{2}$. The latter case, in the presence of interorbital hopping, will feature an anomalous Hall effect Martin and Batista (2008); Kato et al. (2010). The collinear AF-2 and AF-3 will also break translational symmetries with associated charge modulation driven by terms proportional to $\vec{\varphi}_{\alpha}\cdot\vec{\varphi}_{\beta}\neq 0$; such orders may be favored by repulsive interactions between neighboring sites.

(c) Nematic FM (${\cal N}$FM): Starting with a uniform FM state, the quartic coupling $\lambda_{2}>0$, if sufficiently large, can drive charge nematicity, since it can change the ‘mass’ of the nematic field to $(r_{\psi}-\lambda_{2}|\vec{\varphi}_{0}|^{2})$. This coupling between FM and nematic orders is not linear in $\psi_{n}$, unlike the above AFM/ICS cases. Thus, the nematicity in this case is not symmetry-enforced. For $w_{\psi}>0$, the ${\cal N}$FM will have depletion of the density in one orbital, as we find from the MFT.

Fluctuation/disorder effects. — In 2D, without SOC, thermal fluctuations will destroy long-range magnetic order at any $T>0$. In this case, ICS, the ${\cal N}$FM, the collinear ${\text{AF-}1}$, and orthogonal AF-2, will melt into a charge nematic, reflected in a nonzero $\psi_{n},{\cal S}_{n}$, which will undergo symmetry restoration via a $Z_{3}$ clock (or $3$-state Potts) transition. Within MFT this is a first order transition, but thermal fluctuations render it a continuous transition Wu (1982). The non-coplanar AF-3 state will lead to a magnetically melted state with only chiral order (linked to the $\pm$ choice of $\hat{\Omega}_{3}$), featuring a nonzero anomalous Hall effect that vanishes above an Ising transition at which time-reversal symmetry is restored Kato et al. (2010).

We speculate that disorder might also weakly suppress long-range magnetic order, even with SOC, leaving vestigial nematic order Nie et al. (2014) down to $T=0$. This suggestion is motivated by Sr3Ru2O7, where the observation of nematic transport Borzi et al. (2007) near the metamagnetic critical point was recently attributed, via neutron scattering, to arise from nearly ordered SDW phases Lester et al. (2015).

Discussion. — The electronic nematic phases we have proposed in (111) 2DEGs will lead to transport anisotropies. On symmetry grounds, the scaled resistive anisotropy $(\rho_{xx}\!-\!\rho_{yy})/(\rho_{xx}\!+\!\rho_{yy})$ will track the nematic order parameter Fradkin et al. (2000); Basak and Carlson (2017). A simple Drude picture (see SM) shows that, for the coordinates used above, $\rho_{xx}\!-\!\rho_{yy}\!\sim\!{\rm Re}\,\psi$ while $\rho_{xy}\!=\!\rho_{yx}\!\sim\!{\rm Im}\,\psi$. Further signatures of nematic order may be observed in Friedel oscillations which can be probed using scanning tunneling spectroscopy. Even in a conventional phonon-induced SC state, the presence of such background nematic order would lead to an anisotropy of the vortex shape and the mobility as well as the critical current, explaining the anisotropy observed in the resistive transition into the SC state Davis et al. (2017c). If such orbital order is weak, it will be less evident in ARPES Rödel et al. (2014); McKeown Walker et al. (2014) than in transport probes.

Rashba SOC Zhong et al. (2013); Park et al. (2013, 2012, 2011) does not significantly impact the (111) FS for relevant densities $\bar{n}\sim 0.3$, or lead to a significant spin-splitting near the tips of the elliptical FSs where we find the magnetic instability (see SM) since orbital mixing is negligible at those momenta. Thus, we expect SOC will not significantly modify the phase diagram at these densities; however, it can pin the magnetic order or convert the uniform FM into a long wavelength spiral Banerjee et al. (2013); Li et al. (2014). SOC will have a more significant impact on low density 2DEGs, and transport properties which average over the entire FS of all bands.

Our work has not taken into account random oxygen vacancies - these can locally pin the nematic order but cannot induce anisotropies in macroscopically averaged resistivity measurements. However, such ‘nematogen’ defects could amplify weak resistive anisotropies, both of the bulk tetragonal phase in dilute 2DEGs, as well as of the higher density nematic phases with orbital order. This interplay, which has been studied in the pnictides Gastiasoro et al. (2014), would be worth exploring in the oxide 2DEGs.

Finally, tetragonal lattice distortions, described by a 3-state Potts theory Aharony et al. (1977) similar to ${\cal L}_{\psi}$, will couple linearly to the nematic order parameter $\psi_{n}$ itself, affecting the SDW degrees of freedom as well. For instance, a tetragonal distortion with elongation of the ${c\text{-axis}}$ will favor the nematic order associated with the AF-1 state. The interplay of electronic nematicity explored here, with anisotropies induced by surface phonon mechanisms, is an important topic for future research.

We thank F. Y. Bruno, Q. Li, L. Miao, A. Chubukov, and E. Fradkin for useful discussions and feedback. This research was supported by NSERC of Canada. AP acknowledges the support and hospitality of the International Center for Theoretical Sciences (Bangalore) during completion of this manuscript.

I Supplemental Material

I.1 Details of RPA calculation

The bare susceptibility, which is the building block of the RPA approximations, can be written in terms of the single particle Green’s functions as

[TABLE]

Given the unitary matrix $\underline{u}({{\bf{k}}})$ which diagonalizes $\underline{h}({{\bf{k}}})$, and its eigenvalues $\varepsilon_{m}({{\bf{k}}})$, the Green function is

[TABLE]

Using this we compute the RPA response ${\underline{\chi}^{(c,s)}_{RPA}({\bf{q}},\Omega)=\underline{\chi}_{0}({\bf{q}},\Omega)(1-\underline{U}^{(c,s)}\underline{\chi}_{0}({\bf{q}},\Omega))^{-1}}$. In order to understand the instabilities, it suffices to focus on the matrix product $\underline{U}^{(c,s)}\underline{\chi}_{0}({\bf{q}},\Omega=0)$; if its eigenvalues exceed unity, it signals an RPA instability. The RPA results for different electron densities are presented in figure 3. In the small $J$ regime and for small densities, the leading instability is a $2k_{F}$ instability in the spin channel in the $\Gamma\to$ M direction, which is commensurate for specific densities and incommensurate in general. The $\bar{n}=0.3$ electrons per site is an example of a commensurate order, which we study in the main text using a $2\times 2$ unit cell. For larger densities, when $2k_{F}$ exceeds the length of the reciprocal lattice vectors, the ordering wavevector moves to the M $\to$ K direction. When Hund’s coupling is large, we get a FM transition corresponding to a leading instability at the $\Gamma$ point.

I.2 Spin-orbit coupling

Surface and interface 2DEGs are subject to Rashba spin-orbit coupling (SOC). This arises from a combination of inversion symmetry breaking at the interface, combined with atomic SOC.Park et al. (2011, 2012, 2013) The inversion breaking term is odd under ${{\bf{k}}}\rightarrow-{{\bf{k}}}$, and can be phenomenologically incorporated by allowing for inter-orbital hoppings which are ordinarily forbidden by the presence of inversion symmetry. Such terms have been labelled ‘orbital Rashba’ terms.Park et al. (2011, 2012, 2013) For example, an electron in the $xy$ orbital could not move to a $yz$ orbital along the $\hat{a}$ direction if inversion symmetry is not broken. However, the broken symmetry at the interface will distort the orbital cloud, leading to an inter-orbital overlap which will be either stronger or weaker depending on the hopping direction (e.g., an $xy\rightarrow yz$ hopping is strong in the $-\hat{a}$ direction but weak in the $+\hat{a}$ direction), giving rise to inter-orbital hopping terms $\propto\sin({{\bf{k}}}\cdot\textbf{r})$. When combined with atomic SOC, this leads to the Rashba effect which has been shown to play a role in (001) oxide 2DEGs. Here, we generalize this idea to the (111) 2DEG.

Explicitly, for the (111) 2DEG in momentum space, we get for the inversion breaking ‘orbital Rashba’ term $H_{OR}=\sum_{{{\bf{k}}},\ell\ell^{\prime},\sigma}c^{\dagger}_{\ell\sigma}({{\bf{k}}})h^{R}_{\ell\ell^{\prime}}({{\bf{k}}})c_{\ell^{\prime}\sigma}({{\bf{k}}})$, where

[TABLE]

and $\ell\equiv(yz,zx,xy)$. Atomic spin-orbit coupling takes the form $H_{A}=-\lambda_{A}\;\vec{L}\cdot\vec{S}$ with $\lambda_{A}\approx 20$ meV (see e.g., Ref.Zhong et al., 2013). In the $(c_{yz,\uparrow},c_{zx,\uparrow},...,c_{xy,\downarrow})^{T}$ basis, we can write this as $H_{A}={\mathrm{i}}\frac{\lambda_{A}}{2}\sum_{{{\bf{k}}}}\varepsilon^{\phantom{\dagger}}_{\ell mn}\tau^{n}_{\sigma\sigma^{\prime}}c^{\dagger}_{\ell\sigma}({{\bf{k}}})c^{\phantom{\dagger}}_{m\sigma^{\prime}}({{\bf{k}}})$, where $\varepsilon_{\ell mn}$ is the Levi-Civita symbol and $\tau^{n}$ the Pauli matrices.

What is the impact of such atomic SOC and inversion breaking on the Fermi surface? Comparing figure 4(a) ($\lambda_{A}=0,\lambda_{OR}=0$) to figure 4(b) ($\lambda_{A}=20$ meV, $\lambda_{OR}=0$) we see that the main effect of atomic spin-orbit coupling on the Fermi surface is the disappearance of the small band near the $\Gamma$ point, which gets pushed below the Fermi level. The leading magnetic instability uncovered in our RPA (reported in the main text and above) occurs at the wavevector connecting the tips of the elliptical Fermi surfaces. Near these tips, we find that the Fermi surface is nearly unaffected which we can understand since the Fermi surface near these regions is almost entirely composed of a single orbital. Using $\lambda_{A}=\lambda_{OR}=20$ meV (see Ref.Zhong et al., 2013), we obtain the Fermi surface shown in figure 4(c). We can see that the Rashba spin-orbit coupling breaks the spin degeneracy, previously protected by inversion and time-reversal symmetry. The tips of ellipses are once again nearly unaffected by this term, the main changes to Fermi surfaces being in the $\Gamma\rightarrow\text{K}$ direction where the different orbitals tend to strongly hybridize. However, even here the changes are small. We include figures 4(d) and 4(e) to show the effect of even larger inversion breaking on the Fermi surfaces; however we point out that these Fermi surfaces are not consistent with ARPES data on the STO surface 2DEG, suggesting that they are too large to be of relevance.

I.3 Details of the MFT calculation

The mean field theory can be cast in the form of a variational minimization where we pick the variational state to be the ground state of a variational Hamiltonian ${H_{\rm var}=H_{0}-\sum_{i\ell}(\phi_{i\ell}+\mu)n_{i\ell}-\sum_{i\ell}{\bf b}_{i\ell}\cdot{\bf S}_{i\ell}}$, with $\{\phi_{i\ell},{\bf b}_{i\ell}\}$ being variational parameters which allow in general for site and orbital dependent order parameters. The expression for the variational free energy is given by:

[TABLE]

where $\underline{u}(\mathbf{k})$ is the matrix which diagonalizes $H_{\rm var}$, $\underline{n_{F}}$ is the Fermi-Dirac distribution, and the thermal expectation values $\langle.\rangle$ are taken with respect to the equilibrium density matrix of the variational Hamiltonian, i.e., ${\rm e}^{-H_{\rm var}/T}$, and $S$ is the fermion entropy for given field configuration $\phi_{i\ell}$ and ${\bf b}_{i\ell}$.

We discuss two cases in the main text, a single-${{\bf{Q}}}$ spiral mean field theory which can be formulated directly in momentum space with a single site unit cell, and a case where we study broken symmetry states with $2\times 2$ real space unit cell which allows us to also explore multi-${{\bf{Q}}}$ order parameters for the commensurate AF.

I.3.1 Spiral MFT (zero temperature)

The $T=0$ spiral mean field theory is performed with only intra-orbital orders, uniform density and spiral magnetization with fixed wavevector ${\bf{q}}$, so $\langle n_{i\ell\sigma}\rangle=\frac{1}{2}\rho_{\ell}$, and $\langle S^{\pm}_{i\ell}\rangle=m_{\ell}{\rm e}^{\pm{\mathrm{i}}{\bf{q}}\cdot{{\bf{r}}}_{i}}$. This leads to the effective Hamiltonian

[TABLE]

where $\Psi^{\dagger}({{\bf{k}}})\equiv(c^{\dagger}_{yz,\uparrow}({{\bf{k}}}),c^{\dagger}_{zx,\uparrow}({{\bf{k}}}),c^{\dagger}_{xy,\uparrow}({{\bf{k}}}),c^{\dagger}_{yz,\downarrow}({{\bf{k}}}+{\bf{q}}),c^{\dagger}_{zx,\downarrow}({{\bf{k}}}+{\bf{q}}),c^{\dagger}_{xy,\downarrow}({{\bf{k}}}+{\bf{q}}))^{T}$, and

[TABLE]

Here

[TABLE]

with $\varepsilon^{\ell(0)}_{{\bf{k}}}$ being the bare intra-orbital dispersion for orbital $\ell$, and $\gamma_{{\bf{k}}}$ being inter-orbital hybridization. Using expression (17), we compute the energy for various ${\bf{q}}$ and pick the wavevector which minimizes the energy. This leads us to the phase diagrams shown in Figs. 2(a-c) of the main text.

I.3.2 MFT with $2\times 2$ real space unit cell

For the $2\times 2$ unit cell, the most general $\bf b$ field can be expressed in terms of 4 vector parameters ${\mathfrak{b}}_{\alpha}$ ($\alpha=0,1,2,3$),

[TABLE]

where $(m,n)$ correspond to lattice coordinates of site $i$ along the $\hat{a}$ and $\hat{c}$ directions respectively. We thus formulate this mean field theory partly in momentum space where we work in the corresponding reduced BZ. The paramagnetic state corresponds to $\mathfrak{b}_{\alpha}=0$ for all $\alpha$, and the ferromagnetic state to only $\mathfrak{b}_{0}\neq 0$. The different AF states all correspond to $\mathfrak{b}_{0}=0$ and different nonzero $\mathfrak{b}_{\alpha}$ (for $\alpha\neq 0$). The collinear AF-1 (or AF-2 or AF-3) state are obtained when respectively only a single component of one (or two or three) $\mathfrak{b}_{\alpha}$ is nonzero. The orthogonal AF-2 (or AF-3) corresponds to two (three) orthogonal choices of nonzero $\mathfrak{b}_{\alpha}$ fields. Minimizing the free energy, we are led to the phase diagram shown in Fig. 2(d) of the main text.

As seen from the figure, there are large regimes where we find robust AF-1 and FM phases. We also find smaller regimes where AF-2 and AF-3 phases are stabilized which involved the superposition of two or three ${{\bf{Q}}}$ orders, but the free energy per site for these states is lower than the AF-1 state by a very small amount (of order $\sim 10^{-5}t$). Furthermore, the free energy difference per site between different superpositions (collinear, coplanar, orthogonal) in the multi-${{\bf{Q}}}$ states (AF-2, AF-3) are even smaller, $\sim 10^{-7}t$.

We have studied the temperature dependence of a few points in the AF phases and found that as the temperature is raised, the AF-2 and AF-3 become an AF-1 phase, consistent with Nandkishore et al. (2012). The mean field critical temperature into the PM starting significantly deep in the AF and FM phases is $\simeq 0.18t$; without SOC, this will be suppressed to zero, leaving only vestigial orbital and nematic orders. A numerical study of the renormalized $T_{c}$ for such vestigial orders taking full thermal fluctuations into account will be discussed elsewhere.

I.4 Nematicity in the Drude conductivity

Simple Drude-like considerations are enough to determine how the electric conductivity depends on the nematic order parameter $\psi_{n}$. We begin by assuming that only $\psi_{n}$ has condensed and that there is no magnetic order. For simplicity, we further assume that the inter-orbital hopping is zero. Thus, each orbital, denoted by its preferred hopping direction, $\ell=a(xy),b(xz),c(yz)$, has a separate contribution to the conductivity tensor, $\hat{\sigma}=\sum_{\ell}\hat{\sigma}_{\ell}$. In terms of Cartesian coordinates, $x,y$, used in the main text, with $\hat{a}=\hat{x}$, $\hat{b}=\hat{x}/2+\hat{y}\sqrt{3}/2$, and $\hat{c}=-\hat{x}/2+\hat{y}\sqrt{3}/2$, we get

[TABLE]

where, in general $\sigma_{1}\neq\sigma_{2}$. A six-fold rotation operation, $C_{6}$, relates the tensors above: $\hat{\sigma}_{b}=C_{6}\hat{\sigma}_{a}C_{6}^{-1}$ and $\hat{\sigma}_{c}=C_{6}^{-1}\hat{\sigma}_{a}C_{6}$. Within a Drude picture, the matrix elements are proportional to the density of electrons, in this case, for each orbital separately. Therefore, we can approximate $\hat{\sigma}\approx\sum_{\ell}\rho_{\ell}\hat{\sigma}_{\ell}/\rho_{0}$, where, $\rho_{\ell}$ is the electron density in each orbital, while the average density per orbital is denoted by $\rho_{0}$. The nematic order parameter captures the imbalance in orbital densities, $\psi_{n}\sim((\rho_{a}-\rho_{0})+\omega(\rho_{b}-\rho_{0})+\omega^{*}(\rho_{c}-\rho_{0}))/\rho_{0}$. Putting everything together, we find

[TABLE]

Finally, assuming $\psi_{n}$ is small, one can obtain a qualitatively similar dependence of the resistivity tensor $\hat{\rho}$ on $\psi_{n}$.