Strain-induced topological phase transition at (111) SrTiO$_3$-based heterostructures

The quasi-two-dimensional electronic gas at the (111) SrTiO$_3$-based heterostructure interfaces is described by a multi-band tight-binding model providing electronic bands in agreement at low energies with photoemission experiments. We analyze both the roles of the spin-orbit coupling and of the trigonal crystal field effects. We point out the presence of a regime with sizable strain where the band structure exhibits a Dirac cone whose features are consistent with \textit{ab-initio} approaches. The combined effect of spin-orbit coupling and trigonal strain gives rise to non-trivial spin and orbital angular momenta patterns in the Brillouin zone and to quantum spin Hall effect by opening a gap at the Dirac cone. The system can switch from a conducting to a topological insulating state \textit{via} modification of trigonal strain within a parameter range which is estimated to be experimentally achievable.

A general property of heterostructures, among which LAO/STO, is the possibility of controlling strain at the interfaces by realizing superlattices or using micromembranes [23][24][25]. The presence of strain allows the manifestation of phenomena absent in the bulk structure [26][27][28][29]. Moreover, these phenomena can be tuned by manipulating the strain itself. In particular, the possibility of strain-tuned topological phase transition has been theoretically predicted [30] and experimentally realized [31].
After the discovery of the presence of the q2DEGs on interfaces with different crystallographic cut directions [32], the attention was extended to the (111) LAO/STO interface, in which even more exotic phenomena have * mtrama@unisa.it † cataudella@na.infn.it ‡ carmine.perroni@unina.it been predicted [33]. Recently anisotropic magnetoresistance both in normal [34] and superconducting state [35,36] has been observed, suggesting the unconventional nature of superconductivity host by (111) LAO/STO interface. Further, recent studies [37] have discovered a striking similarity between this interface and graphene due to the presence of a Dirac point in the band structure. This similarity arises by virtue of the similar honeycomb lattice structure formed at the interface when STO is Ti-terminated. The analogy between graphene and the q2DEG in the (111) interface suggests the possibility of non-trivial topological effects [38]. In fact, since the Haldane seminal discovery [39], it has been known that q2DEG on honeycomb lattices in an insulating state can manifest Quantum Hall Effect (QHE) in the absence of magnetic field. The subsequent discovery of the Quantum Spin Hall Effect (QSHE) by Kane and Mele [40,41] showed that topologically protected edge states can appear in graphene even in presence of time reversal invariance, when SOC is included. This suggested that topological insulators (TIs) can be achieved in materials with strong SOC [42][43][44]: this has been later experimentally confirmed [45][46][47]. Further studies have highlighted the presence of TI in 3D bulk [48] and their nanostructures [49,50]. Finally, Xiao et al. [51] proposed the use of Transition Metal Oxides (TMOs) based heterostructures, in particular perovskites, in the (111) direction for realizing TIs. This possibility is offered by the presence of strong atomic SOC typical of the TMOs and the peculiar geometry of (111) perovskites. The authors used Tight Binding (TB) methodologies to show the presence of a QSHE in different materials. This methodology has been widely used to compute theoretically the electronic band structure for LAO/STO (111) [33,52,53].
In this work we analyze specifically the properties of the (111) STO interface electronic structure using an effective TB model, focusing on the LAO/STO one. Com-pared to the previous works which adopted TB methodologies we emphasize the role of the strain at interface and the atomic SOC in the band structure properties. Such an approach is useful to have a physical understanding of the mechanisms giving rise to the peculiar features of the system. We perform a comparison with the experimental data available on the (111) interface. In particular we fix the hopping parameters from the Angle-Resolved Photoemission Spectroscopy (ARPES) experiments on the bare (111) STO surface. We identify the region of the parameter space in which the Dirac point emerges within the electronic band structure; we find that this happens when a relevant term in the TB treatment is the strain. As already mentioned, the presence of a Dirac Point in the electronic band structure of the (111) LAO/STO interface was already predicted in [37], where Density Functional Theory (DFT) was adopted. The use of TB methodology allows an effective modeling of the band structure near the Dirac point including also the effects of the SOC with a small computational effort and, then, open up the way to study the transport properties and the topological invariants of this interface that can be very demanding for "first principle" approaches. We verify the concrete possibility of the regime with sizable strain studying the relation between the energy parameter in the TB model and the geometry of the physical structure. Furthermore we include an analysis of the SOC role in the properties of the energy spectra. We also investigate the presence of spin and orbital angular momentum patterns in the Brillouin Zone (BZ). In the last part of the work we perform a systematic scan of the parameter space to identify the region in which the system exhibits a QSHE. We verify the bulk-boundary correspondence by explicitly showing the presence of edge states on finite-thick zigzag ribbons. Our results suggest the possibility of a strain-induced topological transition for the system. This paper is organized as follows: in Section II the model Hamiltonian for the STO 2D bilayer is introduced, discussing every contribution we take into account. The methods used to estimate the hopping parameters is also presented. In Section III the properties of the energy spectra are explored, emphasizing the role of the trigonal strain; we also discuss the physical range in which the strain is expected to vary. Furthermore the presence of spin and angular momentum patterns in the BZ is shown. In Section IV the topological properties of the system are discussed together with the methods used to determine the topological invariant Z 2 . In Section V we report our conclusions.

II. MODEL AND METHODOLOGY
In this section we discuss the electronic band structure of the q2DEG using a minimal TB approach. This methodology allows us to encode the physical properties of the crystal in a small number of energy parameters in order to understand their role in the band structure formation. The projection of two different Ti atom layers along the (111) direction of the STO results in a honeycomb structure, similar to the graphene one, with a lattice constant ofã = 0.3905 2/3 nm as shown in Fig. (1). The conduction band is obtained by the superposition of the 3d orbitals of Ti atoms of the STO. Due to the approximate cubic geometry these orbitals split up into a triplet t 2g = {d yx , d zx , d xy } and a doublet e g = {d z 2 , d x 2 −y 2 } of orbitals, as described in Ref. [54]. The t 2g orbitals are responsible for the formation of the lower conduction band, so that in the following we will only consider them in the description of the energy structure. We consider the two dimensional structure originating from the projection on a single plane of only two layers of Ti, whose atoms are labeled using the symbols Ti 1 and Ti 2 . We denote by d iασ, k the annihilation operator of the electron of 2D dimensionless quasi-momenta k =ã K, where K is the quasi-momentum, occupying the orbital i = {xy, yz, zx} belonging to the layer α = {Ti 1 , Ti 2 } and of spin σ = ±1/2. With this notation we studied the following Hamiltonian We will discuss the different contributions separately. We start from the TB Hamiltonian which in the k-space takes the form where t αβ i ( k), considering also the next-to-nearest neighbour (NNN) hopping, contain three contributions: • the first and most important contribution will be the head-on overlap of two lobes of the same kind of orbitals for nearest neighbours. This happens between a Ti 1 and a Ti 2 , so that this term will describe an interlayer hopping. The amplitude for this hopping will be denoted by t 3 ; • the second contribution will come from the lateral overlap of two lobes of two orbitals of the same type between Ti 1 and Ti 2 . The amplitude of this interlayer term is denoted by t 2 ; • the third and smallest term will be a NNN intralayer hopping caused by the head-on overlap between two orbitals of two Ti on the same layer. The amplitude of this term will be denoted as t 1 . From best fit calculations this term will be neglected, as we can show in the following.
The matrix representation of the Hamiltonian (2) is given in Appendix A (Eq. (A1)). The values of the hopping parameters have been fixed by comparison with ARPES experiments [55] to the values t 3 = 0.5 eV; t 2 = 0.04 eV; t 1 = 0 eV, obtained as best fit parameters of a maximum likelihood analysis which is detailed in the Appendix A.
H T RI takes into account the possibility of having strain at the interface along the (111) direction. The physical origin of this strain is the possible contraction or dilatation of the crystalline planes along the (111) direction. The strain at the interface can be manipulating via deposition of other crystal, such as LAO. This coupling is taken into account by the following contribution, whose form is fixed by symmetry considerations [56] H T RI = ∆ 2 In the following we will refer to this term as the trigonal crystal field. More details about this term will be given in the next section. We initially point out that the sign of ∆ determines the kind of the applied strain along the (111) direction: if ∆ < 0 a contraction of the Ti atom planes occurs, while for ∆ > 0 a dilatation occurs. Another contribution taken into account is the atomic SOC, which is represented in the k-space by where ε is the Levi-Civita tensor, and σ k are the Pauli matrices. Finally we include a possible interaction with an external electric field orthogonal to the planes of the form where ξ T i1 = +1 and ξ T i2 = −1. This contribution will take into account the possibility of a small perturbation at the interface breaking the inversion symmetry between the two layers.

III. TRIGONAL CRYSTAL FIELD EFFECTS
The purpose of this section is to illustrate the role of the trigonal crystal field in determining the properties of system. We will show two extreme situations: a first one in which ∆ is a small perturbation compared to the TB and a second one in which ∆ is the dominant parameter. In order to exemplify the behaviour of the system for small values of ∆, we choose a benchmark ∆ = −0.005 eV, equal to the trigonal coupling for the unstrained interface [57]. For this choice the energy band structure and the density of states are shown in Fig. (2). At the K point of the Brillouin zone one can notice an energy dispersion which is linear in lower bands: this could erroneously suggest that a Dirac cone is present within the energy bands. In Appendix B we discuss a more complete characterization of the energy spectra near the K point, finding that the bands do not have the typical conical structure of the Dirac point in this region.
The band structure drastically changes when the parameter ∆ is assumed to be larger, exhibiting a Dirac cone in the low energy region as shown in Fig. (3). In fact when |∆/t 3 | 1 there is a splitting of the energy levels (neglecting the spin degeneracy) in a pair of bands of lower energy and four bands of higher energy. In the limit of |∆/t 3 | 1 these bands asymptotically become the eigenstates of the trigonal Hamiltonian called respectively a g and e π g [56], belonging to the corresponding symmetry group. The hopping contribution causes only a weak mixing between the two groups of states and gives rise to a Dirac cone at the K point in the Brillouin zone between the a g states. Our results are in agreement with the results obtained using DFT methodologies in Ref. [37] We can locally adopt a perturbative approach near the K point to obtain an effective Dirac Hamiltonian from which we can extract the Fermi velocity of the electron at the Dirac point, which is v F = ( 1 2 (t 2 + 2t 3 ))/ ≈ 2.5 × 10 5 m s −1 , and ∆ SO is the effective coupling induced by the SOC. It depends on the parameters as ∆ SO ∼ (λt 2 3 )/∆ 2 , since it originates from the perturbation to the third order in the TB and SOC contributions. A more detailed discussion is provided in Appendix C. The Dirac point is sensitive to the SOC, which opens a gap between the bands; this is described by the ∆ SO term in the Hamiltonian (6). The two dimensional carrier density n 2D needed to fill the bands up to the Dirac point is of the order of 10 14 cm −2 (for ∆ = −1 eV, n 2D = 7.6 × 10 14 cm −2 , corresponding to 0.34 electrons for each Ti ion). This is comparable with the typical electron densities injected at the LAO/STO interface, making it possible to perform an experimental investigation via field effect technique.
The presence of the trigonal distortion induces some peculiar behaviour of the angular and spin polarization between the two layers of Ti atoms. In order to show this we evaluated for a fixed filling µ the following quantity: where O α is the generic component of the angular momentum or spin momentum projected on the α-th layer, |E, k is a generic eigenstate of energy E and quasimomentum k and P α is a projection operator on the states localized on the α-th layer where we defined {|ψ iασ, k } as the vector of the original basis of the Hamiltonian. These quantities evaluated throughout the Brillouin zone are depicted in Fig. (4) for a benchmark choice of the parameters. This Figure  shows that a magnetic moment localized in k is formed between the two layers, in principle observable via spinpolarized ARPES experiments. A peculiar feature appears in Fig. (4): when the chemical potential is within the gap formed at the Dirac point, the out-of plane spin component is totally localized near the K point in the Brillouin zone, while the orbital angular momentum is almost zero. The latter can be understood by observing that in the limit −∆/t 3 1 the lower orbitals tend to the a g states which are characterized by L = 0. The localization of the spin, instead, can be understood by observing that far from the Dirac point the splitting among the first and second doublet is induced by TB alone, and therefore the mean value of the spin vanishes. On the other hand nearby the Dirac point Eq. (6) shows that the splitting induced by TB vanishes. The gap at the Dirac point is induced by the SOC, represented by the ∆ SO term in Eq. (6), leading to a non vanishing mean value of the spin on each of the two layers oriented along the (111) direction. Of course, in the absence of inversion symmetry, e.g. v = 0 in our case, the sum over the two layers of the projected mean values vanishes. This is reasonable and in agreement with literature [58,59]. Furthermore, since the system is time reversal invariant in the absence of a magnetic field, it does not show magnetic behaviour and a many body interaction is needed in order to recover ferromagnetic effects: for this reason, the total spin integrated over the whole Brillouin zone is always vanishing.

A. Connection between trigonal coupling and physical strain
We have shown that the possibility for the system to exhibit non trivial properties is critically connected with the strength of the strain parameter ∆. One can therefore wonder whether it is possible under realistic experimental conditions to achieve the necessary values for the trigonal coupling in order to have |∆/t 3 | 1. In Refs. [60,61] an analysis of the trigonal strain is performed in the case of a face-sharing octahedral structure. In such a case a trigonal distortion is a compression or a stretching of the face-shared planes of the oxygen octahedra around a metal ion. When the material is unstrained the octahedral oxygen cage split up the d-orbitals into t 2g and e σ g orbitals. When a strain is imposed the angle M-O-M between two metal ions and a oxygen in between changes, reducing the symmetry of the system from octahedral to trigonal and splitting up the states into three subgroups: e σ g , a g and e π g orbitals. The method of Refs. [60,61] allows a unified treatment of the energy splitting caused by the trigonal distortion and recovers the energy splitting of the octahedral crystal field when the M-O-M angle takes its undistorted value. The same calculation can be used for the case of a corner-sharing structure as in the STO case. The strain at the interface is translated into a rigid shift of the atoms along the (111) direction. This leads to a distortion of the octahedral oxygen cage structure around one Ti atom, analogous to the case of the face-sharing octahedra. This leads to the energy splitting between the a g , e π g and e σ g states depending on the distortion angle θ between the direction connecting the O ligand with the Ti ion line and the (111) direction as shown in Fig. (5). This energy splitting explicitly depends on the following parameters: • the energy splitting E 0 between the e σ g states and the triplet of states a g and e π g , which are degenerate when ∆ = 0. This splitting is typically expressed • the effective nominal charge z * which affects the electrons in the conduction band. This is evaluated by subtracting from the ion charge Z the charge of the electrons in the complete shells and of those electrons in the non-complete shells multiplied by a corrective factor (0.35), which takes into account the electron-electron screening.
Using these results for the STO structure, choosing E 0 = 2 eV as reported in Ref. [62], κ = 4.93, as determined for the d orbitals, and z * = 22 − 18 = 4, since Ti 4+ ions have no electrons in their external shells, we reproduce the energy dependence of the a g , e π g and e σ g orbitals as a function of the distortion angle θ, together with the behaviour of the strain parameter ∆ = (E ag − E e π g )2/3 in Fig. (6). For the latter we show the results also for κ = 12.5, corresponding to the muffin tin approximation [61]. The results show that for experimentally accessible distortion angles, ∆ reaches absolute values sufficiently large to cause the appearance of the Dirac cone in the band structure. In fact, compared to the equilibrium position θ 0 = arccos (1/ √ 3) ≈ 54.7 • for which there is no splitting between the a g and e π g states, for θ ≥ 58.7 • we find ∆ ≤ −0.6 eV. The muffin tin approximation leads to even smaller values of distortion angle necessary to reach the same value of ∆. For this reason for the rest of the paper, unless we specify otherwise, we consider as benchmark the value κ = 4.93 which leads to conservative estimates for ∆. We observe that our results do not take into account the correction due to the crystal field originated by the metal ions in the structure; the introduction of such a correction only increases the strength of the trigonal coupling in absolute value. Our results are therefore meant as a lower bound on the realizable coupling ∆, and they show that, already at this level, the trigonal dominated regime proposed in the previous sections seems to be experimentally achievable. These conclusions are further strengthened in Appendix D, in which we show that they are robust by accounting for the modification of the hopping parameters induced by the distortion angle. Let us finally discuss the change in the in-plane lattice constant induced by the biaxial strain along the (111) direction. As reported in Ref. [63], a biaxial (111) strain in the bulk STO could lead to transition to a ferroelectric phase if the strain is tensile. We first point out that, in the realistic case corresponding to ∆ = −0.65 eV and a distortion angle of 2 • evaluated using the muffin tin approximation, the compression of the lattice constant is about 5%. In order to maintain the volume of the lattice cell constant, the in-plane hexagonal lattice has to expand by about ∼ 2.5% of the lattice constant along the [110] and [112] directions. This would result in a dilation of the side of the hexagon by the same relative amount, preserving the C 3v symmetry of the 2D lattice. For this reason, we neglected this effect in our work. Furthermore, given that the considered strain for negative ∆ is compressive, we expect a paraelectric phase for the STO, according to Ref. [63], which does not modify the predicted results, up to a redefinition of the electric field coupling v.

IV. TOPOLOGICAL PROPERTIES
Due to the simultaneous presence of the SOC, the Dirac cone in the electronic band structure and the time reversal invariance, we expect that, at least for |∆/t 3 | 1, the system could exhibit QSHE. In this section we scan the parameter space of the system in order to determine which are the conditions for this to happen.  Since in the strained interface all the conditions for QSHE are met, we simulated a finite-thick (20-100 sites) zigzag ribbon, whose representation is shown in Fig. (7), obtaining numerical evidence of edge states. The latter are represented by the band inversion in the energy dispersion along the k X coordinate, the momentum in the or-thogonal direction to the finite length of the ribbon. The results for a 20 sites ribbon are represented in Fig. (8), where we highlight in color the edge states. Notice that the SOC is chosen to be 0.1 eV for ease of visualization (the typical values are smaller [52]). We also show the value of the Z 2 invariant computed for each of the Kramers doublet according to the method of the Wannier Charge Centers (WCCs). In this method we identify the fundamental direction of the reciprocal lattice, named λ 1 and λ 2 , and we compute the WCCs, namely the mean value of the position i ∂ ∂λ1 , for both of the states of a Kramers doublet as a function of the periodic variable λ 2 . For each doublet we subsequently count the number of discontinuities of the mean WCC along half a period of λ 2 . The sum modulo two of the number of discontinuities is equal to Z 2 invariant. The specific details of how to choose a continuous gauge along the BZ are discussed in more details in Refs. [64,65]. The edge states are always formed in between the non trivial bands. The localized nature of these states is confirmed by the rapid decreasing of the occupancy probability away from the edge for the red band. The choice of showing the results for a 20 sites ribbon is to highlight this rapid decreasing. For a 100 sites ribbon we find, as expected, the same qualitative behavior. A particularly interesting conclusion is the fact that edge states can form at the K point, in between the Dirac cone; this constitutes a further similarity between our heterostructure and graphene. The general topological characterization of the system has been obtained by computing the Z 2 topological invariant for each of the six Kramers doublets and subsequently determining the phase diagrams in the λ − ∆ plane. In the left panel of Fig. (9) we show the phase diagram for the lowest energy doublet. The diagrams for the other doublets can be found in the Appendix E. We can deduce from the figure that the first energy doublet is always non-trivial when λ and ∆ are non-zero. In order to achieve a QSHE the necessary requirements a chemical potential within the energy gap in the band structure and a value equal to 1 for the Z 2 invariant; the latter is the sum modulo 2 of the Z 2 of the filled doublet. A topologically non-trivial behaviour is therefore expected when the chemical potential lies between the first and second Kramers doublet. We emphasize the behaviour of the first doublet, since it is easier to access experimentally, but, referring to Fig. (8), a non trivial behaviour is expected also if the bands are filled up to the second doublet. Finally we studied the stability of the topological properties against an external electric field. In order to do this we evaluated the Z 2 invariant for a fixed value of λ = 0.1 eV, varying v and the trigonal parameter ∆ in a region in which the system exhibits an insulating behaviour. The results are depicted in middle panel of Fig. (9), showing how the topological structure for the lowest bands is preserved for electric potentials of v ∼ 0.004 eV for λ = 0.1 eV. This corresponds to an electric field between the two layers of Ti of the order of 10 7 V/m over a thick-  ness ofã ∼ 0.3 nm. Generally with the decreasing of the strain |∆| a larger field is needed to bring the system into the normal state. We studied the stability against weak perturbations also varying the spin-orbit coupling, since it could in principle be different from the value used in the previous analysis: we fixed the trigonal parameter to ∆ = −1 eV obtaining the phase diagram shown in the right panel of Fig. (9). In order to prove the possibility to induce a topological transition via modification of the strain we show in Fig. (10) the energy gap between the bands as a function of the strain parameter, together with the chemical potential for a fixed value of the surface carrier density n 2D . We have chosen this density value so that the chemical potential µ lies within the energy gap for ∆ = −0.65 eV and for benchmark values of λ = 0.1 eV and v = 0 eV. From Fig. (10) it is evident that decreasing ∆, and therefore increasing the strain, the chemical potential moves through the bands passing in the orange region. Here the two conditions of a non-trivial topological structure for the lower band and an insulating regime for the material are simultaneously met, so that the system can exhibit a QSHE.

V. CONCLUSIONS AND DISCUSSION
In this work we investigated the properties of the (111) STO-based heterostructure interface via TB methodologies. We took into account the contribution of the electron hopping in the bilayer of STO at the surface, the interfacial trigonal strain, the SOC and the electric field orthogonal to the bilayer. The values of the main hopping parameters were obtained by comparison with the ARPES experiments. The remaining parameters were systematically studied in order to identify the features of the electronic structure. We identified two different regimes: one in which the strain is a small perturbation compared to the hopping parameters and the other in which it is a relevant contribution. In the former situation our model is able to capture the main features of the band structure previously obtained in literature [53,55,66]. In the latter the system exhibits a Dirac point in the structure of the low energy bands, for experimentally accessible values of the carrier density. Compared with Ref. [37], which obtained the same behaviour with DFT methodologies, we accurately discussed the role of SOC and assessed the regime of the trigonal strain where a Dirac cone is expected. We found that the presence of the Dirac point requires values of the trigonal distortion such that |∆/t 3 | ≥ 1. Our assumption was supported by a systematical study of the dependence of the energetic parameter ∆ on the physical strain applied at the interface. We also showed the pres-ence of non trivial spin and orbital angular momentum patterns in the BZ, finding that in the presence of the trigonal distortion there can be a non vanishing magnetic dipole moment between the two layers of Ti at a fixed k in the Brillouin zone. In agreement with the expectation of time reversal invariance this magnetic moment vanishes when integrated over the whole Brillouin Zone; nevertheless these results are potentially accessible using spin-polarized ARPES experiments. Interestingly in the regime with sizable strain, when the chemical potential is within the gap at the Dirac point, the out-of plane components of the spin are completely localized near the K point in the Brillouin zone. These features suggest a non-trivial spin behaviour and open up the way to a systematic analysis of the spin properties of (111) STO heterostructures. We determined these properties in a minimal model using only two layers of Ti atoms. The presence of thin films over the STO, in particular LAO layers, are included as a reservoir of electrons, and its possible role in generating distortion is completely contained within the strain parameter ∆ and v. Moreover, we did not include in our study the Coulomb interaction. This allowed to obtain flexible results from the model with small computational effort: indeed we plan to use this model as a starting point for studying different phenomena. Furthermore, we expect the Coulomb interaction to be weak in the low filling region. Its effect should mainly be a renormalization of the energy spectra for the occupied bands [53], increasing locally the Fermi velocity. This conclusion should still be valid in the regime of interest |∆|/t 3 > 1 near the Dirac point, as shown in Ref. [67]. Here it is shown that the Coulomb interaction cannot open a gap in the band structure by itself, even though it widens the energy gap opened by SOC. Therefore, the main effects of Coulomb interactions that are not included in our work should be a renormalization of the Fermi velocity and the energy gap induced by SOC at the Dirac point. We point out that we did not investigate a poissible anti-ferrodistortive transition which can occur at low temperature in the bulk STO, and using (111) biaxial strain as predicted in Ref. [63]. Moreover we did not investigate the possible presence of charge ordered modes as predicted in [37] under suitable conditions. These questions can be investigated in future studies.
The low energy region in the regime with sizable strain is particularly interesting also due to its non-trivial topological properties. We computed the topological invariant of the system using the WCC method proposed by Soluyanov and Vanderbilt [64,65]. We found that the first Kramers doublet is characterized by a Z 2 = 1 in the presence of the inversion symmetry. This property is preserved in the presence of weak perturbations due to an external electric field. Finally, we showed that the system can exhibit a QSHE; this property has been investigated via finite-thick zigzag ribbon simulations, which show that edge states can be present under suitable conditions at the Dirac point. We found that a topological transition can occur since the system can be tuned from a conducting to an insulating state via strain manipulation.
In conclusion, as a result of a sizable strain, our TB scheme predicts a non-trivial topological phase for LAO/STO (111) including a gap at the Dirac point. The possibility to engineer the applied trigonal strain at the interface suggests that the transition to the non-trivial topological phase can be controlled in a parameter range which is experimentally achievable. In this appendix we discuss the methodology used to extract the values of the hopping parameters of the matrix t αβ i in Eq. (2). In literature STO bare surface ARPES data are available [55] which allow us to fix the hopping parameters t 1 , t 2 and t 3 , since the qualitative behaviour of the interface depends on the ratio between them and ∆. As a first approximation we can fix ∆ = 0, since its value is of the order of the meV for the bare STO [57], λ = 0 and v = 0. In this framework we can extrapolate analytically the eigenvalues of the system Hamiltonian, reported here in a matrix form for simplicity 1 : where the intralayer, diagonal, contributions are: 1 The double spin degeneracy is understood. and the interlayer contributions are: where k = k Xû 110 + k Yû 112 is the dimensionless in-plane k vector of the Brillouin Zone. We emphasize that the coupling among different t 2g orbitals are all equal by virtue of the C 3 symmetry of the interface. The possibility of exact diagonalization is evident due to the fact that the tight binding interaction does not mix the different d i orbitals with one another, but only the two layers. The matrix is therefore in a 2 × 2 block diagonal form.
The eigenvalues, which are doubly degenerate in the spin, are the following: where i runs over yz, zx, xy, γ i ( k) in which we have explicitly expressed the k dependence. In k = 0 (Γ point in the Brillouin Zone) the eigenvalues referring to the lower band (− sign) have a threefold degenerate minimum (neglecting the spin degeneracy) Near the bottom of the bands an approximate expression for the energies is .
These quantities represent the effective masses in the free particle approximation along the directions connecting respectively the Γ and K points and the Γ and M points in the Brillouin zone. The numerical values of the parameters have been fixed using the ARPES data [55] shown in Fig. (11): from the dashed curves we have extrapolated the values of the effective masses, together with a confidence range. In this way, a χ 2 analysis has been performed in order to search for the best fit combination of the three parameters. The best fit values are: show in Fig. (12) the χ 2 contours in the t 3 -t 2 plane. From the contours we see that the number of digits for the parameters is chosen coherently with their uncertainty on them. Analogous studies in literature [66] have found similar values for the parameters. In Fig. (11) (solid line) we show the effect of a weak trigonal strain on the pure TB curves. The strain changes the effective masses near the Γ point and induces a splitting between the first and the other two doublets of the order of 3/2∆. However, the splitting is not sufficient to induce an effective gap between the first two doublets and the third one.
with the sign referring to the lower or the upper bands, We have neglected the spin index for simplicity. Thus we can write the effective Hamiltonian as with λ = ± discriminating between the upper or the lower bands, and d λσ i,k is the annihilation operator referring to the states in Eqs. (B3). We conclude that the K point in this case is quite different from the Dirac point in graphene: in the latter, in fact, there are no states near the Dirac point, while in our case there is a large quantity of states corresponding to the intersection of the three plane bands. This is also seen from the density of states in Fig. (2) at the dashed line. In order to experimentally realize an electron filling up to the K point an electron surface density of n 2D = 1.59 × 10 15 cm −2 is needed, which is large compared with the typical charge carriers densities measured or predicted at the interface under examination. We will now discuss the stability of the degeneracy of the K point upon including the perturbations. First of all, using the trigonal crystal field coupling of ∆ = −0.005 eV the band structure changes as in Fig. (14). The degeneracy is partially split, without the formation of a gap. The presence of v does not influence the existence of the degeneracy: its effect is limited to pushing the upper and lower band further and further apart. The spin-orbit coupling, on its own, is not able to remove the degeneracy as shown in Fig. (15). On the other hand the simultaneous presence of the trigonal field and the spin orbit interaction removes completely the degeneracy which leads to the intersection of the bands. Of course, for sufficiently small values of λ and ∆, the effect is limited only to an extremely narrow region near K, and the behaviour returns to be quasi-linear as we go far from the K point. The case of the strained surface is particularly interesting due to the presence of the structure resembling the Dirac Point. Since ∆ = −1 eV it is possible to analyze the proximity of the point thanks to a perturbative approach in which the tight binding terms are treated as small. In fact we can consider the two lower bands as originating only from the a g Ti 1 , k and a g Ti 2 , k Bloch orbitals belonging to layer Ti 1 and Ti 2 : this would be an exact statement in the case of infinite ∆. Even though ∆ = −1 eV is not exactly in the range of applicability of perturbation theory (because t 3 = 0.5 eV), this would be enough to give us a semi-quantitative understanding of what happens at the K point. In the subspace of degeneracy constituted by the two states a g Ti 1 , k and a g Ti 2 , k the perturbation induced by the TB is Near the K point the latter can be expanded to the first order, due to the fact that the sum of the vanishes in that point: (C2) where we used the notation of δk X = (k X − 4π √ 27 ). We can write this Hamiltonian introducing the very compact notation of the Pauli spin matrices σ i : where δk is the k vector centered around the K point, v F = 1 2 (t 2 +2t 3 ) and we neglected the constants because we are interested in the qualitative behaviour. Using t 2 = 0.04 eV and t 3 = 0.5 eV we obtain v F ≈ 2.5 × 10 5 m s −1 which is smaller than the graphene one (v F ∼ 10 6 m s −1 ). The eigenstates of this Hamiltonian are where the phase e iθ k = δk X ∓ik Y | δk| , so that θ k coincides with the geometrical angle between δk and k X . The degeneracy near the K point is repeated six times throughout the hexagonal Brillouin zone: however, these six points are not all equivalent. In particular, near the K point the Hamiltonian is the same with a change in sign in σ x . If, for this second class of Dirac points, we reverse the ordering of the states in the matrix form of the Hamiltonian, we find that the form of the Hamiltonian is equal and opposite to the one near K. Thus, with the new ordering of states, we can write the effective Hamiltonian near K as The eigenstates of this Hamiltonian now become The degenerate points therefore group into three points equivalent to K and three points equivalent to K , forming two independent sublattices. This is again exactly what happens in graphene, where in the usual treatment a two-valued degree of freedom is introduced, generally named valley degree of freedom, corresponding to which class of Dirac point is under analysis. In the space of this new degree of freedom one can introduce new Pauli matrices, which we will denote by τ i , so that the full Hamiltonian can be written as where the state is now a four-spinor. We conclude that there is a resemblance between graphene's Dirac cone and (111) LAO/STO interface. An experimentally relevant difference is the fact that the Dirac point in the graphene stands exactly at the neutrality point of the system, between the conduction and valence band, while the Dirac point in the LAO/STO interface resides fully The electric field at the interface is fixed to v = 0; the spin-orbit coupling is taken as λ = 0.01 eV.
in the conduction band. Therefore, in order to investigate that point we need to fill sufficiently the |A g lower bands: the numerical predicted value (for our choice of parameters) is n 2D = 7.6 × 10 14 cm −2 . This value is comparable in order of magnitude with the typical charge carrier density predicted at these interfaces by the polar catastrophe model. We now discuss the presence of the effective emergent coupling due to the SOC at the Dirac point. The energy splitting induced by the SOC has to be treated carefully since, as we mention in the text, the mean value of L is equal zero on the a g states, leading to no induced splitting to the first order in the SOC. Moreover, both the SOC and the TB terms are to be treated as small perturbations to the larger trigonal strain term. We explicitly verified that a gap among the a g states at the Dirac point only opens to the third order in perturbation theory: in particular the gap opening is well-described by the perturbative term in the Hamiltonian with In these formulas for simplicity we set t 2 = 0. Therefore the complete Hamiltonian on the Dirac points is well described by: Appendix D: Strain influence on hopping parameters In this section we will discuss the effect of the strain on the hopping parameters. We will demonstrate that, even if the trigonal strain does change the hopping parameters, mainly due to the change in the bond lengths, the regime of |∆|/t 3 > 1 described in the text can always be reached. As we discussed in the main text, a trigonal distortion is a rigid shift of one layer with respect to the other along the (111) direction. First of all, this shift preserves the relative magnitude of the hopping parameters in different directions: since the trigonal strain preserves the C 3 symmetry, it does not give rise to new types of TB couplings violating C 3 , and therefore it only affects the numerical values of the parameters t 1 , t 2 , and t 3 . The numerical change is induced by two factors: the different angle in the projection of the orbitals, and the different bond lengths in the two centers Slater-Koster (SK) overlap integrals where ψ lm ( r) is the wave function of the electrons localized on r and V is the ionic potential which mediates the hopping between the two sites. We discuss the two dependencies in turn. On the one hand, the cosine of the angle δ between the Ti-Ti joining line and the z direction, which would be 1 in the equilibrium configuration (see Fig. (16)), changes in the distorted configuration. Since the (SK) integrals depend regularly on the cos(δ), this will affects the hopping parameters by powers of cos δ 1 − δ 2 /2, so that even for large angle of distortion, e.g. 10 • ∼ 0.17 rad, the correction of the hopping parameters is of the order of 1.5%. We will neglect the correction due to the distortion angle in comparison with the one due to the bond lengths. On the other hand, the distances between two Ti atoms are changed by the distortion of the crystalline planes. The unperturbed angle between one side of the cubic lattice and the (111) direction is θ 0 54.7 • . Under contraction of the crystalline planes, this angle increases to θ 0 + δ.The distance between two Ti atoms is decreased from a 0 to a 0 sin θ 0 / sin(θ 0 + δ). The SK integrals, and therefore the t 3 parameter, depend on the distance as d −(l+l +1) [26], and since we are in the subspace of the {t 2g } which have a quenched angular momentum, l = 1, l = 1 and thus t 3 ∝ d −3 . For a representative value δ = 10 • (which is large compared to the angles needed to reach the regime |∆| > t 3 ), the hopping parameter t 3 changes from 0.5 eV to 0.5 eV/(0.9) 3 0.7 eV. Therefore, due to the small values of the distortion angles needed to reach the regime of |∆|/t 3 > 1, the change in the hopping parameters is not expected to significantly influence our results. We emphasize in any case that our results are based on choices for the numerical parameters which are conservative: for example, for the constant κ we have taken a value κ = 4.93. The linearized muffintin orbitals approximation predicts larger values for κ, leading to |∆| which are even larger, therefore requiring even smaller distortion angles to enter the regime of |∆|/t 3 > 1. In Fig. (16) we show the dependence of t 3 and ∆ as a function of the distortion angle. We can see that, since t 3 grows slower than |∆|, it is possible to reach the |∆|/t 3 > 1 regime already for δ ∼ 3 • . In Fig. (17) we show the evolution of the band structure for decreasing ∆, taking into account the correction to the hopping parameters due to the distortion angle δ. The result shows that already for δ = 6 • the Dirac cone is present, corresponding to ∆ = −0.81 eV.