Topological states at the (001) surface of SrTiO3

TL;DR

This paper reveals that the (001) surface of SrTiO3 hosts topological edge states due to inverted parity of electronic wavefunctions, suggesting new topological phenomena in this material.

Contribution

Theoretical demonstration of topological edge states at SrTiO3 (001) surface arising from inverted wavefunction parity and Dirac cone dispersion.

Findings

Presence of four gapped Dirac cones with anisotropic dispersion

Existence of topological edge states at the boundary

Predicted Josephson radiation signatures of topology

Abstract

Defect-free SrTiO3 (STO) is a band insulator but Angle Resolved Photoemission Spectroscopy (ARPES) experiments have demonstrated the existence of a nanometer thin two-dimensional electron liquid (2DEG) at the (001) oriented surface of this compound. The bulk is a trivial insulator, but our theoretical study reveals that the parity of electronic wavefunctions in this 2DEG is inverted in the vicinity of special points in reciprocal space where the low-energy dispersion consists of four gapped Dirac cones with a tilted and anisotropic shape. This gives rise to linearly dispersing topological edge states at the one-dimensional boundary. We propose to probe these modes by measuring the Josephson radiation from gapless bound Andreev states in STO based junctions, as it is predicted that they display distinctive signatures of topology.

Equations

$$H_{\text{surf/2DEG}}=H_{\text{kin}}+H_{\text{so}}+H_{\text{om}}.$$

$${ϵ}_{d_{yz}}=$$

$${ϵ}_{d_{xy}}=$$

$${ϵ}_{d_{xz}}=$$

$$H_{\text{so}}=\lambda \mathbf{L}\cdot \mathbf{S}$$

$$⟨xy;\sigma ⟩H_{om}yz;\sigma =2ig\sin k_x$$

$$⟨xy;\sigma ⟩H_{om}xz;\sigma =2ig\sin k_y$$

$$\left[\begin{array}{cccccc}{ϵ}_{d_{xy}}& 2ig\sin k_x& 2ig\sin k_y& 0& {\lambda }^{{}^{\prime }}& -i{\lambda }^{{}^{\prime }}\\ -2ig\sin k_x& {ϵ}_{d_{yz}}& i{\lambda }^{{}^{\prime }}& -{\lambda }^{{}^{\prime }}& 0& 0\\ -2ig\sin k_y& -i{\lambda }^{{}^{\prime }}& {ϵ}_{d_{xz}}& i{\lambda }^{{}^{\prime }}& 0& 0\\ 0& -{\lambda }^{{}^{\prime }}& -i{\lambda }^{{}^{\prime }}& {ϵ}_{d_{xy}}& 2ig\sin k_x& 2ig\sin k_y\\ {\lambda }^{{}^{\prime }}& 0& 0& -2ig\sin k_x& {ϵ}_{d_{yz}}& -i{\lambda }^{{}^{\prime }}\\ i{\lambda }^{{}^{\prime }}& 0& 0& -2ig\sin k_y& i{\lambda }^{{}^{\prime }}& {ϵ}_{d_{xz}}\end{array}\right]$$

$$|B⟩$$

$$|C⟩$$

$$\left[\begin{array}{cccc}{E}_b& R_1^{\ast }& 0& 0\\ R_1& E_b& 0& 0\\ 0& 0& E_c& R_2^{\ast }\\ 0& 0& R_2& E_c\end{array}\right],$$

$$|B_1⟩$$

$$|C_1⟩$$

$${\text{BO}}_{\text{b}}=|B_1⟩+|B_1^{\ast }⟩$$

$${\text{ABO}}_{\text{b}}=|B_1^{\ast }⟩-|B_1⟩$$

$${\text{BO}}_{\text{c}}=|C_1⟩+|C_1^{\ast }⟩$$

$${\text{ABO}}_{\text{c}}=|C_1^{\ast }⟩-|C_1⟩$$

$$\left[\begin{array}{cccc}{E}_b-\left|R_1\right|& v\sin k_y& w1\sin k_x& 0\\ v\sin k_y& E_c+\left|R_2\right|& 0& -w1\sin k_x\\ w1\sin k_x& 0& E_c-\left|R_2\right|& v\sin k_y\\ 0& -w1\sin k_x& v\sin k_y& E_b+\left|R_1\right|\end{array}\right]$$

$$h(k)=\left[\begin{array}{cc}ϵ(k)+M(k)& i\alpha k_x-\gamma k_y\\ -i\alpha k_x-\gamma k_y& ϵ(k)-M(k)\end{array}\right]$$

$$E=-B_x{k}_x^2+M-B_y{\lambda }_{+}{\lambda }_{-}+\frac{\alpha }{\gamma }{B}_y{k}_x({\lambda }_{+}+{\lambda }_{-}),$$

$$\left[\begin{array}{cc}M-(D_x+B_x)k_x^2-(D_y+B_y)k_y^2& i\alpha k_x-\gamma k_y\\ -i\alpha k_x-\gamma k_y& -M+(-D_x+B_x)k_x^2+(-D_y+B_y)k_y^2\end{array}\right].$$

$$\begin{gathered} {\lambda }_{\pm }^2=\frac{B_x}{B_y}{k}_x^2+\frac{{\gamma }^2}{2B_y^2}-\frac{M}{B_y}\pm \\ \sqrt{{\left(\frac{{\gamma }^2}{2B_y^2}-\frac{M}{B_y}\right)}^2+\frac{E^2-M^2}{B_y^2}+\frac{k_x^2}{B_y^2}\left(\frac{B_x}{B_y}{\gamma }^2-{\alpha }^2\right)} \end{gathered}$$

$$|{\varphi }_B⟩=\left[\begin{array}{c}a\\ b\end{array}\right](e^{-{\lambda }_{+}y}-e^{-{\lambda }_{-}y}),$$

$$\frac{a}{b}=\frac{M-E-(B_x{k}_x^2-B_y{\lambda }_{+}^2)}{-i(\alpha k_x+\gamma {\lambda }_{+})}=\frac{M-E-(B_x{k}_x^2-B_y{\lambda }_{-}^2)}{-i(\alpha k_x+\gamma {\lambda }_{-})},$$

$$E=M-B_x{k}_x^2+\frac{\alpha }{\gamma }{B}_y{k}_x({\lambda }_{-}+{\lambda }_{+})-{\lambda }_{+}{\lambda }_{-}{B}_y,$$

$$\sqrt{M-E}=\text{sign}(B_y)\sqrt{M+E},$$

Full text

I Introduction

Material science has sparked a recent interest in the research community with the emergence of certain innovative compounds with desirable transport properties. Of equal interest are the compounds that are endowed with interfacial or surface states whose properties vastly differ from their bulk parent compounds. For example, one finds metallicity confined at the boundary between two wide band gap bulk insulators, LaAlO3 (LAO) and SrTiO3 (STO) Ohtomo and Hwang (2004). Such unusual behavior ”at the edge” is truly an interfacial effect with the bulk of the two materials playing no role. This might be reminiscent of topological matter where unconventional states develop at the edge Ren et al. (2016).

Apart from these unexpected interfacial effects, ARPES studies have revealed that even the surface of (001) oriented STO, in the presence of oxygen vacancies, harbours a metallic state while the bulk remains insulating. It is to be noted that bulk STO is a band insulator in the absence of these vacancies. The metallic state or the 2D electron gas (2DEG) at the surface has a band structure which is very different from its 3D counterpart. It is not simply a 2D cut of the bulk band structure, as its shows essentially no dispersion in the direction perpendicular to the 2DEG Plumb et al. (2014).

In addition to the metallic state present at the surface, theoretical proposals advocating the presence of topological states in (111) oriented transition metal oxides have also been made for heterostructures Xiao et al. (2011); Liang et al. (2013); Rüegg and Fiete (2011); Rüegg et al. (2013) and subsequently surfaces, where a metallic state already exists Bareille et al. (2014); Rödel et al. (2014) . According to ARPES measurements, (111) oriented STO and KTaO3 (KTO) have a stongly confined 2DEG. The geometry of the conducting paths of the carriers in the 2DEG for this particular orientation is a honeycomb lattice. If its extent is limited to two unit cells, a theoretical analysis of the parity of eigenstates suggests a Z2 topological state provided that $E_{f}$ is close enough to the position of the topological band gap predicted. However, these two requirements, of limiting the gas to two unit cells and $E_{f}$ being close to the topological gap, are not met experimentally.

In this paper, we uncover a topological character of the conducting electronic states of the 2DEG present at the (001) oriented surface of STO. These features are specific to the surface and have no bearing on the bulk properties of the system. Accordingly, edge states are expected at the 1D boundary of this 2D layer.

While some of the basic ingredients known to be relevant to the physics of topological insulators are also seen here, certain specific features in this system need to be underlined, namely the occurence of several bands of different orbital character and the impact of the confining potential on the orbital ordering.

The scope of the present paper is to bring to light the existence of the 1D topological edge states, stemming from a non trivial Z2 number for the eigenstates in the Brillouin Zone (BZ) in the 2DEG observed in ARPES at the (001) surface of STO.

The first section of our paper summarizes the main experimental findings pertinent to the electronic structure of the 2DEG at the surface of STO which show a dominant $d$ character of the conduction bands, hint at the two dimensionality of the gas, show signatures of the confinement potential and the electric field associated to the Rashba spin-orbit coupling present at the surface. In the second section, we propose a realistic theoretical model to incorporate all the elements listed above, in an energy range of 300 meV around $E_{f}$, which causes us to neglect the higher Sr states and the $e_{g}$ states of Ti. After, we explicit each term which enters the model, confinement, bulk spin-orbit coupling and its combination with orbital mixing at the surface, giving rise to an effective Rashba coupling. We then ascertain the accuracy of the model by overlaying the theoretically determined energy and momentum dispersion curves onto the experimental ARPES data. Near $E_{f}$, the energy bands have a tilted, Dirac-like geometry in the vicinity of four symmetry-related points in the Brillouin Zone (BZ). We argue that these points could be viewed as effective $\Gamma$ points where again the orbitals are degenerate, as is observed in the cubic bulk structure without Spin-Orbit Coupling (SOC). In the next section, we present the eigenstates in the vicinity of the DP and we discuss the symmetries they posses and the winding of the phase of the wave-functions on constant energy contours. Branch cuts of the phases connecting pairs of time-reversed Dirac partners are also shown as well as the parity inversion occuring in the gapped region. An edge state analysis is provided in section four and we find linearly dispersing 1D edge states at the boundary of the 2D surface, by mapping our model to a simpler Bernevig-Hughes-Zhang (BHZ) model Bernevig et al. (2006). The last section includes the conclusions and perspectives where we discuss the potential impact of the magnetic states in the 2DEG. We also discuss possible experiments to provide a fingerprint of these 1D topological states. It is noteworthy to mention that since these states lie a few meV below the $E_{f}$, it is easier to access these states as compared to other transition metal oxides, notably (111) oriented KTO and STO.

II Experimental background

STO is known experimentally to be a trivial band insulator in the bulk, in the absence of defects. In the bulk for the cubic phase, the conduction $t_{2g}$ states $d_{xy},d_{yz},d_{xz}$ are degenerate at the $\Gamma$ point and are unoccupied, giving rise to its insulating behaviour. After the addition of bulk SOC, the $t_{2g}$ states are spin-orbit coupled and the $\Gamma_{7}^{+}$ spin split state is raised to about 25-30 meV above the degenerate $\Gamma_{8}^{+}$ states at the $\Gamma$ point. Below 100 K, there also exists a tetragonal distortion which might remove the degeneracy of the bands but it is largely dominated by SOC. However, away from $\Gamma$, the bands essentially recover their pure orbital character, irrespective of the perturbation. In 2011, interest was revived in STO, when a 2DEG was found to exist at its (001) surface Santander-Syro et al. (2011); Meevasana et al. (2011). Its origin is ascribed to oxygen vacancies confined to the boundary. To characterise it theoretically, the electronic structure, derived from spin-integrated ARPES experiments, was interpreted in terms of sub-bands. The sub-band structure was thought to be caused by the confinement potential, modelled by a triangular quantum well, which pins the 2DEG at the surface. This confinement potential also has another important effect; it overcomes the bulk SOC effect to lower the $d_{xy}$ states by about 200 meV below $E_{f}$, giving metallicity and raising the $d_{xz},d_{yz}$ states above the $d_{xy}$, thus, reversing the order of the bands. The two dimensionality of the gas was also revealed by the lack of dispersion in the $\hat{z}$ direction. The dominant orbitals in the conduction bands were shown to be the $t_{2g}$ $d$ states of Ti and those in the valence bands were shown to be O $2p$ states.

Subsequently, Spin-Polarised ARPES Santander-Syro et al. (2014) was performed and spin polarisation was found for the two sub-bands seen earlier in spin-integrated ARPES. The data was interpreted as spin split bands with oppositely winding spin chiralities at the Fermi surface. The spin textures observed were a consequence of the surface Rashba effect.

An issue raised with the detection of spin polarised bands was the presence of internal magnetism and spin-split bands Santander-Syro et al. (2014); McKeownWalker et al. (2016), which is one we will discuss further when we probe our system for possible topological states. In the next section, we will show the modelling of the bands.

III Modeling the 2DEG

In order to obtain an appropriate low-energy model for the electronic states of the 2DEG, we follow a two-step procedure. First, we construct a tight-binding Hamiltonian that includes the relevant terms capturing the physical properties of the material. In addition to an effective hopping between Ti atoms for each of the three $t_{2g}$ $d$ orbitals, we take into account the confinement close to the surface of STO and in a direction perpendicular to it, denoted by $z$. Confinement quantizes the dynamics along $z$. Furthermore, we include the SOC as well as orbital mixing contributions. Second, once the structure of the model is set by these physical considerations, we determine the values of the effective tight-binding parameters from a fitting to available ARPES spectra. In the insulating bulk, the conduction band manifold has a dominant $d$ orbital character of the charges which leads to an effective Ti $t_{2g}$ and O $2p$ Hamiltonian where hopping is between neighboring Ti and O sites. The directional anisotropy of the bondings implies that the hopping amplitude of $d_{yz}$ carriers is small in the $\hat{x}$ direction, denoted by t${}_{\text{h}}$, but large in the $\hat{y}\text{ and }\hat{z}$ directions denoted by t and is degenerate at the $\Gamma$ point, as can be seen in Fig. 1b) of Santander-Syro et al., 2011. In the presence of bulk SOC, the degeneracy of the bands at $\Gamma$ is lifted. To construct a low-energy Hamiltonian, one needs to take into consideration the orbitals which can play an important role at the surface near $E_{f}$. The reduction of the Hilbert space from a hybridization between O $2p$ states and Ti $t_{2g}$ states is a viable choice as the difference in energy between the $e_{g}$ and the $t_{2g}$ manifolds is of the order of 2.5 eV. The other states, for example, the core states of Ti, O and Sr lie far below the $E_{f}$ and other unoccupied bands of Ti, O and Sr lie much higher up in energy. Thus, a low energy Hamiltonian can be constructed from the reduction of the Hilbert space from a hybridization between O $2p$ states and Ti $t_{2g}$ states and then to an effective model with hopping from Ti to Ti sites. This bulk modelling can be used as a support when we move to the surface but it is important to note that at the surface, new terms in the Hamiltonian will completely change the allure of the bands. This modelling of the conduction bands is also in line with the dominant $d$ orbital charater of the conduction band as revealed in ARPES Santander-Syro et al. (2011); Meevasana et al. (2011) and confirmed in DFT Shen et al. (2012). Mapping a tight-binding Hamiltonian onto the experimental and simulated data is a standard procedure used by the ARPES and DFT communities to extract the relevant energies and effective masses. In order to quantify the interplay between the orbital and spin degrees of freedom, we describe the ARPES determined bands with the help of simple tight-binding surface models, inclusive of kinetic, bulk spin-orbit and surface orbital mixing terms.

The directions chosen for the surface are $\hat{x}$ and $\hat{y}$ while the perpendicular to the surface is in the $\hat{z}$ direction. The surface Hamiltonian $H_{\text{surf}}$ has the form

[TABLE]

$H_{\text{kin}}$ is composed firstly of tight-binding bands derived from the effective Ti-Ti hopping as in the bulk. Secondly, the motion of the carriers within the 2DEG is constrained along $\hat{z}$ by the confinement potential, caused by the oxygen vacancies. Confinement then yields two terms in the Hamiltonian – the first, a global energy offset -$V_{0}$, Santander-Syro et al. (2011) affects all three $t_{2g}$ bands equally and the second, $\epsilon_{1/2}$, represents the two relevant sub-bands due to the $k_{z}$-quantization.

At $\Gamma$, the sub-band value $\epsilon_{2}$ for $d_{xz}$ and $d_{yz}$ states is greater than the energy $\epsilon_{1}$ for the $d_{xy}$ state, reflecting the hierachy of hopping amplitudes along $\hat{z}$ in the bulk. The energies are thus,

[TABLE]

where the in-plane wave-vector components $k_{x}$ and $k_{y}$ are dimensionless (momenta times $a$, the lattice parameter). The $d_{xy}$ states are lowest in energy at $\Gamma$ and the orbital order is reversed as compared to that in the bulk. Consequently one observes several crossing points between one light ($d_{xy}$) and one heavy ($d_{xz}$ or $d_{yz}$) band in the BZ along $\Gamma$X and $\Gamma$Y. However, along $\Gamma$M and $\Gamma\bar{\text{M}}$, one gets special crossing points, where all three orbital energy dispersions are equal denoted by $\Gamma_{\text{dp}}$. The addition of confinement breaks the degeneracy of the orbital energies at $\Gamma$ that one observes in the bulk cubic phase, but in a certain sense this degeneracy is restored at $\Gamma_{\text{dp}}$ bringing a new symmetry at the surface in (and four new ”$\Gamma$” - like points). At $\Gamma_{\text{dp}}$, $\left\lvert k_{x}\right\rvert=\left\lvert k_{y}\right\rvert=k_{c}=2\arcsin{\left\{\sqrt{\left[\left(\epsilon_{2}-\epsilon_{1}\right)/\left(t-t_{h}\right)\right]}\right\}}$.

The next contribution in Hamiltonian (1) is simply the bulk SOC due to the Ti atoms which entangles the spin and orbital degrees of freedom, the expression of which is given by

[TABLE]

At the $\Gamma$ point, confinement has already lifted the degeneracy of the bands and the bulk SOC does not produce a significant effect. However it is an essential ingredient at the $\Gamma_{\text{dp}}$ points and results in making the band crossings partially avoided; one of the bands, $\Gamma_{7}^{+}$ becomes spin-orbit split off and which is raised by 25-30 meV above the other two degenerate $\Gamma_{8}^{+}$. The points at which the two degenerate bands still cross are the called Dirac Points (DP) situated at $\Gamma_{\text{dp}}$. They are four in total, one in each quadrant of the BZ and are related by C4 and time reversal symmetry (see Fig. 1). In the vicinity of the DP, the orbital dispersion resembles Dirac cones. Each of the Dirac cones is tilted and the tilt angle changes as one moves around $\Gamma_{\text{dp}}$. In the first quadrant, it is maximum along $\Gamma$M but goes to zero in the direction denoted by $\Delta$, perpendicular to $\Gamma{\text{M}}$ and passing through $\Gamma_{\text{dp}}$. The bulk SOC does not change the position of the crossing point and the Dirac cones are still found at $\Gamma_{\text{dp}}$. Because of time reversal symmetry, each spin orbital band is twofold degenerate in energy. We note that, while the points $\Gamma_{\text{dp}}$ are situated in directions of high symmetry, their crystal symmetry is lower than that of the special points $\Gamma$ or $M$ and $\bar{M}$. This feature allows for the existence of a tilt of the Dirac cones which would have been otherwise excluded Goerbig et al. (2008, 2009).

The last term in Hamiltonian (1) is the orbital mixing (OM) term, which originates from the electrostatic and the bond angle perturbations that occur near the surface. Given the value of the carrier density of the 2DEG, and considering that the excess electrons in the gas originate from the vacancies, we may consider that the vacancies are fairly homogenously distributed across the surface. The OM term gives rise to a Rashba term whose expression can be found in [Khalsa et al., 2013; Zhong et al., 2013]

[TABLE]

where $g$ is the strength of the OM and $\sigma=\uparrow,\downarrow$ . This produces new coupling terms among orbitals which were earlier forbidden in the bulk by symmetry. The experimentally relevant value of $g$ is at $E_{f}$. It is particularly large whenever $E_{f}$ lies close to the gaps of the band structure (i.e. near the avoided band crossings). A combination of tight-binding and DFT calculations estimate the size of the Rashba-like term to be 5-10 meV which is consistent also with the experimentally determined value. The OM term mixes the spin-orbital states derived from $H_{\text{kin}}$ and $H_{\text{SO}}$ and lifts their degeneracies. In particular gaps develop in the spectrum at $\Gamma_{\text{dp}}$ as shown in Fig 2.

In the basis $\Ket{d_{\beta};\sigma}$ ($\beta=(xy,yz,xz),\;\sigma=\uparrow,\downarrow$), with $\lambda^{{}^{\prime}}=\lambda/3$, $H_{\text{surf/2DEG}}$ reads

[TABLE]

It is to be noted that this Hamiltonian is for the same value of $k$ and its time-reversed block, degenerate in energy, also exists but will not be represented here for the sake of simplicity. Excellent fitting of the ARPES energy dispersion curves is obtained with $a=3.90$ Å, $t=0.388$ eV, $t_{h}=0.031$ eV, $\lambda=0.025$ eV, $g=0.02$ eV and $\epsilon_{1}-V_{0}=-0.205$ eV (1st $d_{xy}$ band) and $g=0.005$ eV and $\epsilon_{1}-V_{0}=-0.105$ eV (2nd $d_{xy}$ band), $\epsilon_{2}-V_{0}=-0.0544$ eV for $d_{yz}/d_{xz}$ bands. With these parameters and no further adjustment, we overlay the tight binding bands onto the momentum dispersion curves at $E_{f}$ and obtain a very good agreement. DFT calculations on slabs with vacancies confirm the values of the bulk SOC and the Rashba SOC. Within the resolution of ARPES, one sees these crossing points Plumb et al. (2014); Rödel et al. (2016). DFT calculations for slabs with vacancies also produce dispersion suggestive of the 3 band crossing Shen et al. (2012); Jeschke et al. (2015). A recent experimental report of the bands in LAO-STO Strocov et al. (2016) also shows two types of crossings but it is to be noted that the resolution of the ARPES measurement in this heterostructure is less than that achieved at the STO (001) surface. Gating of the samples is also a simple way to tune the energies of the surface states in order to exhibit $\Gamma_{\text{dp}}$ -type crossings. In the next section, we will detail the impact of all these terms on the surface spectrum and the surface eigenstates.

IV Nature of the eigenstates and analysis of the phase of the wavefunctions

With the above values of the parameters, the evolution of the spectrum under the addition of each term can be seen in Fig. 1 and Fig. 2. The $k$-independent bulk SOC terms and confinement produce twofold degenerate Dirac cones and another band, twofold degenerate itself, split off from the Dirac cones.

In order to probe the effect of orbital mixing, we examine the energies and the eigenstates when $H_{\text{om}}$ is included step by step in $H_{\text{surf/2DEG}}$, beginning from $g=0$.

a) If we set g=0, i.e at zero orbital mixing, the spectrum in the first quadrant of the BZ, near $\Gamma_{dp}$, is plotted in Fig. 1 showing three bands, each doubly degenerate. The three corresponding eigenvalues, denoted by $E_{b},E_{c},E_{d}$, are arranged in increasing order of the energy; the first two form the tilted Dirac cones positioned about 3 meV below Ef, while the last one lies 25-30 meV higher in energy above the DP and to first order, may be projected out when probing the low energy physics. We are, thus, left with twofold degenerate Dirac bands near $\Gamma_{\text{dp}}$. These cones can be further divided into an upper branch and a lower branch, Fig. 4. The Hamiltonian is then reduced to a block diagonal form with the upper block consisting of eigenvectors corresponding to energy $E_{b}$, i.e. $\Ket{B}$ and its conjugate $\Ket{B^{*}}$ and a lower block corresponding to the energy $E_{c}$ with eigenvectors $\Ket{C}$ and its conjugate $\Ket{C^{*}}$ . These kets are in fact, the eigen states of the Hamiltonian after addition of SOC and confinement in Eq. (1). After the addition of orbital mixing, we will couple these two blocks as we shall see later. $\Ket{B}$ and $\Ket{C}$ can be represented as a linear combination of all three $d$ orbitals and can be formally noted as

[TABLE]

where $a_{1}$, $b_{1}$, $c_{1}$, $a_{2}$, $b_{2}$ and $c_{2}$ are real functions of $k$ and have branch cuts. It is to be noted that the branch cuts do not occur at the same value of the tilt for the upper and lower block.

b) The orbital mixing g firstly contributes to a term which splits each branch of the doubly degenerate cones, upper as well as lower, into bonding and anti-bonding parts. The splitting is anisotropic and direction dependent.

The Hamiltonian in the basis of $\Ket{B}$,$\Ket{B^{*}}$,$\Ket{C}$ and $\Ket{C^{*}}$ reads

[TABLE]

where $R_{1}=4ga_{1}[i\sin(k_{x})b_{1}+\sin(k_{y})c_{1}]$ and $R_{2}=4ga_{2}[i\sin(k_{x})b_{2}+\sin(k_{y})c_{2}]$. The eigenvectors of the first block can be written in terms of the following orbitals.

[TABLE]

The eigenvectors are then given by

[TABLE]

corresponding to $E_{b}-\left\lvert R_{1}\right\rvert$, and

[TABLE]

for $E_{b}+\left\lvert R_{1}\right\rvert$. Similarly we have bonding and anti-bonding orbitals of the same character for the subspace spanned by $\ket{C},\ket{C^{*}}$,

[TABLE]

for $E_{c}-\left\lvert R_{2}\right\rvert$ and

[TABLE]

for $E_{c}+\left\lvert R_{2}\right\rvert$.

The tilted Dirac cones have a zero tilt in the $\Delta$ direction as can be seen in Fig. 4. The spectrum in the other directions with a non zero tilt will simply result in a rotation and shift about the DP. Focusing on the $\Delta$ direction for the sake of simplicity, we observe that the split branches of the cones cross at four points ABCD where we have a clear precursor of band inversion as the bonding state of the upper branch and the anti-bonding state of the lower branch intersect. This inversion is different in nature from the orbital order reversal caused by confinement.

The orbital mixing not only splits the cones but also introduces a quantized phase difference, $\psi$, in the expressions of the eigenvectors, which can be used to characterise the Berry phase of the bands. The phases $\psi_{1}\text{ and }\psi_{2}$, that enter the Eq. (14-17) of the eigenvectors are given by $\tan(\psi_{1})=c_{1}\sin(k_{y})/b_{1}\sin(k_{x})$ and $\tan(\psi_{2})=c_{2}\sin(k_{y})/b_{2}\sin(k_{x})$. The spin texture for a particular phase can then be modelled by a two component vector, as will be seen in Fig. 5. It is also seen that $\psi_{1}=-\psi_{2}$. The parity of the phase is opposite at points B and D as seen in Fig.4, making evident the band inversion that we get even at this first step.

c) The second step of the addition of g gives a term which gaps the spectrum at the crossing points of the anti-bonding and bonding branches, i.e. at all four points ABCD (Fig. 2). The phase presented above is now renormalised to $\psi^{{}^{\prime}}$ at the points ABCD, is no longer quantized and will show jumps wherever the band gap is reached. The Hamiltonian takes the following form in the basis of the eigenvectors explicited above-

[TABLE]

where $w1=-2g\beta_{ab}$, $v=-2g\beta_{ac}$, $a_{2}b_{1}+b_{2}a_{1}=\beta_{ab}$ and $a_{1}c_{2}+c_{1}a_{2}=\beta_{ac}$. With the matrix elements given in Hamiltonian (18) and that the unperturbed energies of the bonding and the anti-bonding orbitals, the perturbed energies can easily be calculated.

The phase winding after step b) can be seen on a constant energy contour obtained right after step a) for one value slightly larger (smaller) than that at A (at C) in Fig. 5 where the four DP, one in each quadrant, are shown for clarity. Two sets of contours are present- one centered around $\Gamma$ and four, oval shaped, lines surrounding the DP. $\psi$ is calculated around one DP in the first quadrant. For the contour above A, a branch cut extending from $\Gamma_{\text{dp}}$ to M is observed. Below point C a branch cut extending from $\Gamma_{\text{dp}}$ to $\Gamma$ is obtained. Even though the branch cut is not along the same direction for bands with energies above and below $\Gamma_{\text{dp}}$, the integral of the Berry phase along all the bands at a fixed value of $k$ still remains zero due to the contribution from the rest of the bands. The cones and the singularity in $\psi$ both evolve with C4 rotation. Two cones which are time reversed partners have a shift of this singularity by an additional $\pi$ between them. These branch cuts are at the origin of a rather unusual behavior of the Berry phase the quantization of which is normally imposed by its single-valuedness on a closed orbit around the Dirac points. Due to the branch cut and the associated $\pi$-jump, one therefore obtains integrated values of the Berry connection that correspond to non-integer “winding numbers” (here half-valued).

Below the DP, the energy contours have phases which wind in the reverse direction to that above $\Gamma_{\text{dp}}$, as can be expected. In between the points A and D, the phase also winds as expected along the energy cuts made in this region.

According to this two-tier sequence of events, the degeneracy for the cones is lifted by $H_{\text{om}}$ and at each step one only considers the mixing of two bands at a time. While this scheme becomes more accurate as the energy difference between the upper and the lower branches increases (i.e. far from $\Gamma_{\text{dp}}$), one should, a priori, consider all four degenerate bands instead of one pair at a time, in the vicinity of $\Gamma_{\text{dp}}$ (ABCD region). However, if we project the approximate eigenvectors obtained after the two step procedure onto the exact eigenvectors obtained by exact diagonalization of the Hamiltonian, we find an order one overlap which justifies the approximation. It is worthy to note that the parity of the anti-bonding orbitals is always odd while that of the bonding orbitals is always even. Thus, at the crossing points B and D where the band inversion takes place, one switches from a bonding to an anti-bonding orbital at D and vice versa at B, thus leading to an inversion in the parity of the eigenstates. With the form of the energies, which resemble gapped and inverted bands and the eigen-vectors which have an inversion of parity at the gapped points, we expect these topological properties to give rise to 1D edge states.

Next, we show that there is a convenient way to follow a variational procedure to model the Hamiltonian allowing us to describe the topological properties of the system in more detail.

V Edge states

We now use a variational approach to find a model Hamiltonian, closely capturing the properties of Eq. (1). The Hamiltonian, which is simply a mapping of the original problem unto a simpler and a better-understood model is a parabolic model, with breaking of the particle-hole symmetry, giving inverted bands with gap openings. Our choice of the parabolic simplification of our tight-binding bands is justified because we investigate the spectrum close to the DP, which is our new effecive ”$\Gamma$” point and in its vicinity, bands can be assumed to be parabolic. The model Hamiltonian reads

[TABLE]

where $\epsilon(k)=C-D_{x}k_{x}^{2}-D_{y}k_{y}^{2}\text{ and }M(k)=M-B_{x}k_{x}^{2}-B_{y}k_{y}^{2}$. The tilt of the Dirac cones, given by $(\epsilon_{d_{xy}}+\epsilon_{d_{xy}}+\epsilon_{d_{xy}})/3$, can be expanded in the vicinity of $\Gamma_{\text{dp}}$ to give $d_{0}(k)\mathbf{1}=(d_{x}k_{x}^{2}+d_{y}k_{y}^{2})\mathbf{1}$. It has been absorbed in the particle-hole-asymmetry term $\epsilon(k)$. $\alpha$ and $\gamma$ represent the anisotropy of the Rashba interaction and are themselves functions of $k$, which incorporate the branch cuts of the eigenvectors and the associated singularity of the Berry phase. The Hamiltonian $h(k)$ in Eq. (19) now looks like the two band BHZ model Bernevig et al. (2006) giving similar band inversion and gap opening. But there are also differences (i) the band gap is very small $\equiv 2$meV (ii) the Rashba coupling terms which are responsible for opening up the gap at the band crossings are anisotropic (iii) the spectrum does not have a rotational symmetry around the DP due to the tilt $d_{0}(k)$.

This model with band inversion at the 2D surface can then be probed for 1D edge states. Following the method detailed in Ref. [Zhou et al., 2008], we break translational symmetry in the $\hat{y}$ direction, while retaining it along $\hat{x}$. This leads to the spectrum (see appendix)

[TABLE]

where the $\lambda_{\pm}$ are the roots of the secular equation $H\Ket{\psi}=E\Ket{\psi}$. In the calculation of the edge-state spectrum, we have omitted the particle-hole asymmetry, i.e. we have used $D_{x}=D_{y}=0$. Indeed, the edge states are a consequence of the underlying topology of Hamiltonian (19) which is not affected by terms proportional to the unit matrix. The spectrum is linear in $k_{x}$ and closes the energy gap at $k_{x}=0$.

We note that the edge states that we have described here connect bands in the vicinity of a given $\Gamma_{\text{dp}}$ point. Time reversal symmetry implies that one also expects edge states connecting bands such that one of these is located near one $\Gamma_{\text{dp}}$ point and the other is near the time-reversed $\Gamma_{\text{dp}}$ partner in the BZ (e.g. points D and D’ in Fig. 6 of the appendix). A mathematical treatment of this type of edge state is beyond the scope of the present paper. The experimental detection of these edge states will be discussed below.

VI Discussion and Perspectives

Our calculation suggests that the Dirac gaps of the low energy spectrum of the 2DEG are located 3 meV below Ef; since the resolution that can be achieved in ARPES is less than this, $\Gamma_{\text{dp}}$ could actually be sitting at Ef, but even if the equality were not strictly met, gating the sample with a few volts in a side gate geometry or with a few tens of volts in a back gate geometry will suffice to bring $\Gamma_{\text{dp}}$ in coincidence with the Fermi energy. However, as we pointed out, in the case of the 2DEG at the (001) surface of STO, the states that we are considering are metallic. The topological edge states that we are predicting at the 1D boundary of the sample might not be easy to detect experimentally as they could be masked by the conducting sheet. One way around this difficulty is to couple the 2DEG to conventional s-wave superconductors. This stategy has been successfully implemented in the context of HgTe/CdTe heterostructures which are known to possess 1D topological edge states. Experiments carried out on HgTe/CdTe quantum wells measured the Josephson radiation from gapless Andreev Bound states Deacon et al. (2016). In out-of-equilibrium situations, missing Shapiro steps and emission at half of the Josephson frequency are predicted Fu and Kane (2009) and they were experimentally measured Wiedenmann et al. (2016). These are consequences of p-wave superconductivity being induced in the 2D topological insulator’s edge channels when it’s Josephson coupled to nearby conventional s-wave superconductors. A doublet of p-wave Andreev bound states arises and has an energy dispersion which is 4$\pi$ periodic in the Josephson phase. We suggest a similar Josephson junction setup for STO where the latter is sandwiched between two superconductors. The topological regime can be observed by gating the sample and tuning the Fermi energy to that of the edge states. A similar change in the frequency emitted should be observed. Optimizing the carrier concentration and the width of the STO junction to reduce the number of parallel channels would allow one to evidence the topological states.

The report of magnetism at the surface of STO raises the issue of the stability of the topological state to time reversal symmetry breaking but from an experimental standpoint, the claim of magnetism is still debatable Rödel et al. (2014)McKeownWalker et al. (2016)Lee et al. (2011)Taniuchi et al. (2016). If a purported ferromagnetic exchange splitting affected the $t_{2g}$ carriers, that would be detrimental to the edge states. The DFT study of [Altmeyer et al., 2016] reveals a dichotomy between localized magnetic moments of $e_{g}$ orbitals and spin textures for the $t_{2g}$ conduction states. Even if spin-polarized domains were present, one might nevertheless retain some of the topological character if the state of the 2DEG at the surface were spatially inhomogeneous. According to experimental reports in STO Taniuchi et al. (2016) and in LAO/STO Bert et al. (2011); Sulpizio et al. (2014); Salluzzo (2015), polarized domains coexist with patches of metallic regions. If the latter percolated across the entire surface of STO, the electronic bands of the 2DEG could still present a parity inversion leading to topological 1D edge states at the boundary.

In the present report, we have shown that a theoretical modeling of the ARPES spectra for (001) oriented STO leads us to uncover the existence of four special points in the 2D Brillouin zone of the 2DEG with energies close to $E_{f}$. In the vicinity of these points, the parity of the electronic wavefunctions becomes inverted. This leads to the appearance of topological 1D states at the boundary of the 2DEG and ways to detect those are suggested.

Acknowledgements.

Authors benefited from discussions with A.F. Santandar-Syro and E. Bocquillon. M. G. gratefully acknowledges support from the Institut Universitaire de France and from the French National Research Agency (ANR) (Project LACUNES No. 274 ANR-13-BS04-0006-01).

APPENDIX

Appendix A Details of edge states

The Hamiltonian 8 can be analysed either perturbatively or variationally to enable mapping onto a simpler model. The second approach consists of using a variational approach, mapping the four band model onto an effective two band model of a similar form as used in Ref. Bernevig et al. (2006). However, we take explicitly into account the particle-hole asymmetry, described by the parameters $D_{x}$ and $D_{y}$ into which we have also absorbed the tilt term $d_{0}(k)$. The Hamiltonian, as explicited in Eq. (3) of the main text (with C=0), is

[TABLE]

In this section, we are interested only in the topological properties of the systems which give rise to the particular edge-state structure. Since these topological properties are not affected by the particle-hole asymmetry, which yields a term proportional to the unit matrix, we set $D_{x}=D_{y}=0$ in the following discussion of the edge states. The model is solved in a finite strip geometry of width $L$ with periodic (open) boundary conditions in the x (y) direction. Here $k_{x}$ is still a good quantum number while $k_{y}$ is replaced by $k_{y}=-i\frac{\partial}{\partial y}$. The wave function is of a mixed real and $k$-space function of the form of $\phi(k_{x},y)=\phi(k_{x})e^{-\lambda y}$. With this form, we derive the secular equation to get four eigenvalues.

[TABLE]

The boundary conditions at the limits impose that the solutions of this equation become zero at $y=0$ , i.e, in the middle of the slab, and go to zero as $y\rightarrow\pm\infty$. The solutions are then of the form, for $y>0$

[TABLE]

and similarly for $y<0$. The condition of existence of such states implies that the $\mathbbm{R}(\lambda_{+})$ and $\mathbbm{R}(\lambda_{-})$ should be of the same sign (so as to avoid any solution which does not decay at $\pm\infty$). The symbol $\mathbbm{R}$ denotes the real part. Under this condition, the equation $H\ket{\phi_{B}}=E\ket{\phi_{B}}$ gives us

[TABLE]

which gives us the following spectrum

[TABLE]

and we see immediately that if we are to get an edge state, we should have the energy $E=0$ at $k_{x}=0$, i.e., the gap is closed at $k_{x}=0$. For the condition $k_{x}=0$, we have

[TABLE]

which has a solution $E=0$ for $MB_{y}>0$, i.e., $M$ and $B_{y}$ should have the same sign, which is the case if we have a band inversion. We also verify that in this case, we shall have $\mathbbm{R}(\lambda_{+})\mathbbm{R}(\lambda_{-})>0$ and thus $\mathbbm{R}(\lambda_{+})$ and $\mathbbm{R}(\lambda_{-})$ have the same sign which entails $0<M/B_{y}<1/4$. Additionally, there can be two kinds of spin-polarised edge states-those present around one Dirac point, i.e. between B and D in Fig. 6 and those spanning the BZ to their time reversed partners, D and D${}^{{}^{\prime}}$.