Orbital and spin order in oxide two-dimensional electron gases
aa r X i v : . [ c ond - m a t . s t r- e l ] A ug Orbital and spin order in oxide two-dimensional electron gases
John R. Tolsma, ∗ Marco Polini, and Allan H. MacDonald Department of Physics, The University of Texas at Austin, Austin Texas 78712, USA Istituto Italiano di Tecnologia, Graphene Labs, Via Morego 30, I-16163 Genova, Italy (Dated: October 10, 2018)We describe a variational theory of multi-band two-dimensional electron gases that captures theinterplay between electrostatic confining potentials, orbital-dependent interlayer electronic hoppingand electron-electron interactions, and apply it to the d-band two-dimensional electron gases thatform near perovskite oxide surfaces and heterojunctions. These multi-band two-dimensional electrongases are prone to the formation of Coulomb-interaction-driven orbitally-ordered nematic ground-states. We find that as the electron density is lowered and interaction effects strengthen, spontaneousorbital order occurs first, followed by spin order. We compare our results with known properties ofsingle-component two-dimensional electron gas systems and comment on closely related physics insemiconductor quantum wells and van der Waals heterostructures.
PACS numbers: 73.20.-r,71.15.-m,71.10.-w, 73.21.-b,75.10.-b
I. INTRODUCTION
Orbital and spin order are relatively ubiquitous instrongly interacting systems, and have long been stud-ied in ABO bulk crystals with the perovskite crys-tal structure. In perovskites B is typically a transi-tion metal and A is either an alkaline or rare earthmetal. Among these materials, those with partially filledtransition metal t g d-orbitals are particularly interest-ing because orbital and spin order appear almost simul-taneously near the magnetic transition temperature. In crystals for which the number of t g electrons pertransition-metal is close to one, the probability of twoelectrons occupying the same site is substantial, andshort-range Hubbard-type interaction models capture themost important parts of the interaction physics. Gen-eralizations of the Hubbard model apply for larger inte-ger values of the number of electrons per metal. As aconsequence of relatively recent advances it is nowpossible to add electrons to quantum wells formed by d perovskites, most commonly SrTiO . The resulting sys-tems are two-dimensional metals with far fewer than oneconduction electron per transition-metal site. Because ofits long range, the typical Coulomb interaction energyof an individual electron in these systems drops to zeroonly as two-dimensional density n / , in contrast to the ∝ n behavior of the Hubbard model. Therefore, in thesesystems the full long-range Coulomb interaction must beretained in order to achieve a realistic description ofelectron-electron interaction effects. In this Article weaddress orbital and spin order, and their interplay in lowcarrier-density perovskite quantum well systems.The possibility of spin order in electron gases with asingle parabolic band was raised very early in the his-tory of the quantum theory of solids. Although manyexperiments have hinted at ferromagnetic instabilities inthese single-band electron gas systems, unambiguousevidence for magnetism has so far been lacking. Theo-retically, order is expected only at extremely low carrierdensities, at least in the absence of disorder.
Our study captures the peculiarities of perovskite multi-bandtwo-dimensional electron gas systems which open up aricher range of experimental possibilities. d b) y -k y kE k xy k c)a) dielectric η TiO FIG. 1: (Color online) A schematic summary of the t g . a) At left is a top-gated SrTiO het-erojunction which hosts t g electrons distributed across two-dimensional TiO layers highlighted in orange. At right isa sketch of the squared amplitude of the subband wavefunc-tions, | η | , for the t g xz and yz subband wavefunctions are shown in red andthe xy subband wavefunction in blue, illustrating the differ-ence between their interfacial confinements. The energy off-set, equal to 2( t − t ′ ) in the single-layer limit, between theanisotropic xz and yz (red) band edges and the isotropic xy (blue) band edge is illustrated in b) and the anisotropy of theelliptical xz and yz Fermi surfaces (red) is contrasted withthe xy band’s circular Fermi surface (blue) in c). t and t ’ arerespectively the large and small hopping amplitudes for t g orbitals discussed in the main text. In multiband two-dimensional electron gas heterojunc-tions and quantum wells, the ground-state energy de-pends sensitively on the density in each spin-resolvedband. We propose a modified variational theory whichcaptures the interplay between orbital-dependent inter-layer electronic hopping, the electrostatic confining po-tential, and electron-electron interactions. Our theory at-tempts to capture physics that is missing in most density-functional-theory calculations which typically em-ploy extended version of the local-density approximation(LDA) for the exchange-correlation (XC) energy func-tional, E LDAXC (cid:2) n ( r ) (cid:3) = Z d r n ( r ) ξ XC (cid:0) n ( r ) (cid:1) , (1)where ξ XC (cid:0) n ( r ) (cid:1) is the XC energy-per-electron of thesingle-band isotropic electron gas model with homoge-neously distributed electron density n ( r ). Despite itsgreat success, it is generally known that the LDA and itsgeneralizations, ( e.g. the generalized gradient approxi-mation ) do not in all cases adequately describe theXC energy. In covalent semiconductors, this shortcom-ing is highlighted by a significant underestimation of theband-gap energy. This limitation can be traced back toEq. (1), where we see that the LDA XC energy functionalis sensitive to the local electron density but not to theorbital symmetries of the occupied electronic states .Motivated by an interest in orbital and spin order in t g two-dimensional electron gases, we have developeda broadly applicable variational theory for the ground-state energy which is sensitive to the electron density ineach band, and therefore is sensitive to electron orbitalsymmetry in systems in which each band has a dominantorbital character.Our article is organized as follows. In Section II webriefly summarize the t g and extend it to account explicitlyfor the concerted action of orbital-dependent electronichopping along the confinement direction and electron-electron interactions. In Section III we describe a varia-tional approximation for the ground-state that accountsself-consistently for interfacial confinement forces andtheir influence on XC energies. As a first applicationof this method, we calculate the phase diagram for the t g t g instabilities. We comment on the ex-perimental signatures of these states and suggest possibleextensions of the calculations presented here. Finally, inSection V we summarize our results and comment on theapplicability of our variational approximation to semi-conductor quantum wells, van der Waals heterostruc-tures, and other multi-component 2D electron systems. II. t g ELECTRON GAS MODEL
We begin by summarizing the simple model for thelow-energy conduction bands in SrTiO quantum wellsand heterojunctions that we study. We first separate thenon-interacting part of the Hamiltonian into two sets ofterms: those which describe electronic motion parallel ( k )to the confinement direction, and those which describeelectronic motion perpendicular ( ⊥ ) to the confinementdirection, which we take to be the ˆ z direction. The pur-pose of this separation will become clear in Section III.The band structure of SrTiO electron gases hasbeen thoroughly discussed in the literature. Here webriefly describe some details required for the present anal-ysis. The low-energy conduction bands of SrTiO quan-tum wells and heterojunctions are formed from the xy , yz , and xz d-orbitals ( the t g orbitals) of Ti. SrTiO hasa perovskite crystal structure, but the sublattice formedby the Ti atoms is simple cubic and described by the lat-tice vector R = a ( n x ˆ x + n y ˆ y + n z ˆ z ) where a = 3 . . Orbital symmetry dic-tates that each t g electron hops mainly between stateswith the same orbital character, with a large hopping am-plitude t in the two directions in the plane of the orbital,and a smaller hopping amplitude t ′ in the third direc-tion. The weak hopping directions for xy , xz , and yz orbitals are ˆ z , ˆ y , and ˆ x , respectively. Nearest neighborhopping on the simple-cubic Ti sublattice is described bythe non-interacting Hamiltonian H = X ασ X RR ′ c † R ασ h α ( R − R ′ ) c R ′ ασ , (2)where h α ( R − R ′ ) = − t δ ∆ R , ± a ˆ x − t δ ∆ R , ± a ˆ y − t ′ δ ∆ R , ± a ˆ z for α = xy − t δ ∆ R , ± a ˆ x − t ′ δ ∆ R , ± a ˆ y − t δ ∆ R , ± a ˆ z for α = xz − t ′ δ ∆ R , ± a ˆ x − t δ ∆ R , ± a ˆ y − t δ ∆ R , ± a ˆ z for α = yz , (3)where the electron spin is labelled by σ and we have defined ∆ R = R − R ′ .2DEGs have been realized in δ -doped SrTiO , inionic-liquid gated SrTiO , at heterojunctions betweenSrTiO and polar perovskites, and in electrostat-ically gated SrTiO . For definiteness, we consider aSrTiO heterojunction extending from z = 0 to z = a ( N l +1) with an electrical gate placed at z <
0. Here, N l is the number of TiO layers (see Figure 1) with appre-ciable occupation. Electron motion in the plane perpen-dicular to the ˆ z -direction is unbounded and Bloch bandsare formed from t g -orbital hopping on the square Ti 2Dsublattices. Electron creation operators in real and 2Dmomentum space are related by c † R ασ = 1 √ N X k e − i k · R ⊥ c † k R z ασ , (4)where R ⊥ ranges over the 2D square lattice projection ofthe Ti simple-cubic lattice, and k is a crystal momentumin the corresponding 2D Brillouin-zone.With this notation, the non-interacting Hamiltoniancan then be written as H = T k + T ⊥ , where the inter -layer hopping term is T k = X ασ X k R z R ′ z c † k R z ασ h k α ( R z − R ′ z ) c k R ′ z ασ , (5)with h k α ( R z − R ′ z ) = − t ′ δ R z − R ′ z , ± a ˆ z for α = xy − t δ R z − R ′ z , ± a ˆ z for α = xz − t δ R z − R ′ z , ± a ˆ z for α = yz, (6)and the intra -layer hopping term is T ⊥ = X ασ X k R z c † k R z ασ h ⊥ α ( k ) c k R z ασ (7)with h ⊥ α ( k ) = ~ k m L − t for α = xy ~ k x m L + ~ k y m H − t + t ′ ) for α = xz ~ k x m H + ~ k y m L − t + t ′ ) for α = yz . (8)In using Eq. (8) we are assuming that the number ofelectrons per Ti atom is much less than one, a condi-tion that is violated only in very narrow quantum wellsystems, so that only states near the Γ point of the2D Brillouin zone are occupied. The energy dispersionof the conduction bands can be described using effectivemasses, which we obtain by expanding about Γ to lead-ing order in k . Eq. (8) shows that the xz and yz bandshave one heavy mass, m H , and one light mass, m L , whilethe xy band is isotropic and has a light mass m L in boththe x and y directions. In terms of the original hoppingamplitudes, m H = ~ / (2 t ′ a ) and m L = ~ / (2 ta ). In heterojunction systems, several xy , yz , and xz typesubbands are expected to be occupied even at moder-ate electron densities. However, since &
75% of theelectron density is contained in the lowest xy , yz and xz subbands , we address the case in which only one sub-band of each orbital type is retained. This truncationis sufficiently realistic to account for the most interest-ing properties of heterojunction 2DEGs. The assumptioncan easily be relaxed when there is an interest in quan-titatively describing the properties of particular 2DEGsystems that have more occupied subbands. With thisassumption Eq. (7) can be written as T ⊥ = X k ασ c † k ασ h ⊥ α ( k ) c k ασ , (9)where α is an orbital label. Eq. (8) gives the bare banddispersions for the three bands retained in our model,which are plotted in Figure 1. According to recent tight-binding fits to Shubnikov-de Haas measurements ofbulk n-type SrTiO , t = 236 meV, and t ′ = 35 meV sothat m H m = Ry t ′ (cid:16) a B a (cid:17) ∼ m L m = Ry t (cid:16) a B a (cid:17) ∼ , (11)where m is the bare electron mass, a B is the Bohr radius,and Ry is the Rydberg energy unit.The isolated-layer bare-electron Hamiltonian in Eq. (7)has several important qualitative features that areschematically depicted in Fig. 1. First, the xz and yz bands have elliptical Fermi surfaces, while the xy bandis isotropic. Second, the xz and yz band edge is higherin energy than the xy band edge by 2( t − t ′ ) ∼
402 meVwhen carriers are confined to a single atomic layer. Thisoffset will decrease as the electrons spread out in the con-finement direction, but the quantum confinement energywill always be lower for xy oribitals. Because the xy bandelectrons have a much weaker hopping amplitude in theˆ z -direction ( t ′ ) than the xz and yz electrons ( t ), theyare more strongly confined to interfaces in heterojunc-tion systems. We note that it is possible to crudelyaccount for the difference in interfacial confinement byintroducing a parameter d , which is defined as the dif-ference between xz , yz subbands and the xy subbandaverage ˆ z -direction positions (see Fig. 1), however in thecurrent article we do not make this simplifying approxi-mation.Next we add electron-electron interactions to ourHamiltonian. When the electron density per Ti atomis much smaller than one, the Fermi surface occupies asmall fraction of the Brillouin zone and the probability oftwo electrons simultaneously occupying the same Ti siteis very small. In this limit, including only the Hubbardpart of the full electron-electron interaction misses themost important Coulomb interactions. Because of itslong range, the typical Coulomb interaction energy of anindividual electron vanishes only as two-dimensional den-sity n / , in contrast to the ∝ n scaling of the Hubbardmodel. We therefore include in our model calculationsthe full Coulomb interaction V = 12 X i = j e /κ q ( r ⊥ ,i − r ⊥ ,j ) + ( z i − z j ) . (12)Here, κ is an effective dielectric constant. For the typ-ical electron-electron interaction transition energies, therelevant dielectric constant does not include the difficult-to-model soft-phonon contribution , but depends on thedielectric environment on both sides of the relevant het-erojunction or surface. In order to take advantage oftranslational symmetry in the x - y plane, it will be conve-nient to Fourier transform the Coulomb interaction withrespect to the in-plane electron position operators r ⊥ ,i : V = 12 A X i = j X q v q ( z i − z j ) e i q · ( r ⊥ ,i − r ⊥ ,j ) (13)where v q ( x ) = 2 πe exp( − q | x | ) / ( κq ). We then sepa-rate the q = 0 term in Eq. (13), which can be regu-larized by combining it with a remote neutralizingbackground. This procedure leads, up to an irrelevantconstant, to the following Hamiltonian contribution: V k = − A X i = j πe κ | z i − z j | , (14)which we refer to below as the Hartree contribution. Be-cause this term in the Hamiltonian depends only on asmall number of macroscopic observables, namely thenumber of electrons in each layer, its contribution to theground-state energy is given exactly by the correspondingclassical mean-field energy, as discussed in Section III B.The remaining portion of the electron-electron interac-tion, V ⊥ = 12 A X i = j X q = v q ( z i − z j ) e i q · ( r ⊥ ,i − r ⊥ ,j ) , (15)gives an energy contribution that is dependent on elec-tronic exchange and correlation effects. III. VARIATIONAL THEORY FOR THEGROUND-STATE ENERGY
In this Section we develop a variational approximationfor the ground-state energy that aims to treat interfa-cial confinement on an equal footing with electrostaticand XC interaction energies. Our strategy is to firstassign confinement wavefunctions to each t g subband,to then use a two-dimensional random-phase (GW) ap-proximation to account for exchange and correlation with a fixed confinement-wavefunction constraint, and finallyto optimize the confinement wavefunctions by minimiz-ing the total energy including electrostatic energies andconfinement-direction hopping. The orbital-dependentsubband wavefunctions can be expanded in the form | ψ ασ i = N l X R z =1 η ασR z | R z α σ i , (16)where | R z α σ i is a single-particle confinement directionlattice state. The single-particle Hilbert space is con-structed as a direct product of the | R z α σ i space andtwo-dimensional momentum space, | ψ k ασ i = | k i ⊗ | ψ ασ i .In Sect. III A we derive an expression for the contribu-tion to the ground-state energy from the single-particleand interaction terms that are sensitive to planar spatialcorrelations: H ⊥ = T ⊥ + V ⊥ (see below.) Because thenumber of electrons in each orbital is a good quantumnumber, the energy can be expressed as a function ofthe subband wavefunction spinors, η ασ , and the electrondensity associated with each spin-resolved orbital, n ασ .In Sect. III B we derive an equation for the wavefunctionsattached to each spin-resolved band that is based on atotal-energy minimization principle. Finally, in Sect III Cwe describe the self-consistent solution of this equation.Using this procedure we can determine the distributionof density amongst the spin-resolved bands, allowing forthe possibility of spontaneous orbital and spin order. A. Two-Dimensional Electron Gas Exchange andCorrelation
In this Section we will apply the random phase approx-imation (RPA) to a 2DEG Hamiltonian constructed byfixing confinement-direction subband wavefunctions: H ⊥ = X ασ c † k ασ h ⊥ α ( k ) c k ασ + 12 A X q =0 X kk ′ αα ′ σσ ′ V ασ,α ′ σ ′ ( q ) c † k + q ασ c † k ′ − q α ′ σ ′ c k ′ α ′ σ ′ c k ασ . (17) H ⊥ accounts fully for lateral hopping and for the com-ponent of Coulomb interactions sensitive to the electroncoordinates perpendicular to the confinement direction.It is simply the second-quantization version of Eq. (15)combined with Eq. (9). The dispersion of the bands re-tained in our model is given in Eq. (8) and plotted inFig. 1. The orbital- and spin-dependent effective interac-tion in Eq. (17) is constructed by attaching the appropri-ate subband wavefunction to each orbital and averagingover the confinement direction V ασ,α ′ σ ′ ( q ) = X R z ,R z | η ασR z | πe e − q | R z − R z | κq | η α ′ σ ′ R z | . (18)In the numerical calculations for the t g i.e. η ασ → η α , allowingthe spin indices ( σ, σ ′ ) in Eq. (18) can be dropped. Theground-state energy of the 2DEG Hamiltonian in Eq. (17)has kinetic energy E ⊥ K , exchange energy E ⊥ X , and corre-lation energy E ⊥ C contributions that depend on a set ofdensities that specify the occupation of each spin-resolvedband, { n ασ } . The exchange and correlation energiesalso depend on the set of transverse subband wavefunc-tions through the expansion coefficients, { η ασ } definedin Eq. (16). Here and throughout the remainder of thisArticle we use the convention { x ασ } ≡ { x xy ↑ , x xy ↓ , x yz ↑ , x yz ↓ , x xz ↑ , x xz ↓ } (19)to express sets of variables that depend on band α andspin σ quantum numbers.The kinetic energy contribution to the ground-stateenergy of Eq. (17) is simply the sum of contributionsfrom each occupied spin-resolved band, i.e. E ⊥ K = X ασ E ασ K . (20)For each spin-resolved isotropic xy band E ασ K = π ~ An ασ m L , (21)and for each spin-resolved elliptical xz and yz band E ασ K = π ~ An ασ m DOS , (22) where A is the 2D sample area and we have defined thedensity-of-states mass m DOS = √ m L m H . The offset be-tween the xz and yz band-edge energy and the xy band-edge energy is accounted for in the analysis of H k pre-sented in Sect. III B.The exchange energy is given by first-order perturba-tion theory. Since the Coulomb interaction vertices arediagonal in band index, the total exchange energy is givenby the sum of independent contributions from each oc-cupied band E ⊥ X = X ασ E ασ X . (23)The exchange energy of electrons in band α and with spin σ is given by E ασ X = − A X k n (F) ασ ( k ) X q =0 V ασ,ασ ( q ) n (F) ασ ( k + q ) (24)where n (F) ασ ( k ) is the usual Fermi-Dirac distribution func-tion. At zero temperature the integral over k can be per-formed analytically and the remaining two-dimensionalintegral over q can be evaluated numerically: E ασ X = − A π Z ∞ q dq Z π dφ X R z ,R z | η ασR z | πe e − qF ( ζ,φ ) | R z − R z | κqF ( ζ, φ ) | η ασR z | Θ (cid:16) √ πn ασ − q (cid:17) × " n ασ − q π r πn ασ − (cid:16) q (cid:17) − π arcsin (cid:18) q √ πn ασ (cid:19) , (25)where Θ ( x ) is the Heaviside step-function and we havedefined F ( ζ, φ ) = q ζ − / cos ( φ ) + ζ / sin ( φ ) . (26)The parameter ζ = m H /m L in Eq. (25) describes thedegree of band-anisotropy present in the xz and yz bands.Eq. (25) applies to the isotropic xy bands after setting ζ to one, in which case the integral over φ becomes trivial.In the numerical calculations presented in Sect. IV wespecifically consider thin-films of SrTiO . Throughoutthis Article we define the “thin-film limit” quantitativelyvia the condition N l a p n/g ≪
1, where p n/g is an approximation for the 2D Fermi wavevector of each ofthe g = 6 spin-resolved bands in our t g N l a is the quantum well thickness ( i.e. the max-imum distance separating any two electrons in the z -direction). Qualitatively this limit describes the casewhen the typical inter-electron distance which is perpen-dicular to the confinement direction, | r ⊥ ,i − r ⊥ ,j | , is muchgreater than the inter-electron distance parallel to theconfinement direction, | z i − z j | , and the latter can be ne-glected. In this limit the Coulomb interaction loses itsorbital dependence, and we therefore replace V ασ,ασ ( q )with v q ≡ πe / ( κq ). Subsequently, Eq. (25) for theexchange energy of a single spin-resolved band can beevaluated analytically E ασ X = − e n / ασ A π / κ ζ / K(1 − ζ ) , (27)where K ( x ) is the complete elliptic integral of the firstkind. Eq. (27) applies to both the elliptical xz and yz bands and, by setting ζ →
1, also to the circular xy band.We note that the exchange energy is always negative ,that the exchange interaction acts only between electronsin the same spin and orbitally resolved band, and thatthe dependence of exchange energy on density in a sin-gle band is superlinear. It follows that exchange favorsstates in which electrons occupy fewer bands. This is thephysical mechanism for the type of ferromagnetism envi-sioned by Bloch. In 2DEGs where different bands havedifferent orbital content (e.g. xy, yz or xz ), the exchangeinteraction can lead to orbitally-ordered states, as wellas spin-ordered states. We discuss these types of insta-bilities in detail in Section IV. It follows from Eq. (27)that the exchange energy is weakly reduced by bandanisotropy. This property works against exchange-drivenspontaneous symmetry breaking in SrTiO of combining coupling-constant-integration with the fluctuation-dissipation theorem.This method begins by considering the coupling-constantof our Hamiltonian ( i . e . the electron’s charge) as a tun-able parameter: we replace − e with − e √ λ and allow λ to vary between zero and one. The ground-state wave-function Ψ ( λ ) and total ground-state energy E ( λ ) bothevolve as λ is varied. At any particular value of λ ,the contribution to the total ground-state energy fromthe coupling-constant-rescaled Hamiltonian in Eq. (17)is found by taking the expectation value E ⊥ ( λ ) = h Ψ ( λ ) |H ⊥ ( λ ) | Ψ ( λ ) i = hT ⊥ i λ + hV ⊥ ( λ ) i λ . (28)After taking the derivative with respect to λ and usingthe Hellman-Feynman theorem this becomes dE ⊥ ( λ ) dλ = 1 λ hV ⊥ ( λ ) i λ . (29)The XC energy is then obtained by integrating over thecoupling-constant E ⊥ XC ≡ E ⊥ (1) − E ⊥ (0) = Z dλλ hV ⊥ ( λ ) i λ . (30)The instantaneous structure factor describes corre-lations between density-fluctuations at wavevector − q inband ( α, σ ), and density-fluctuations at wavevector + q in band ( α ′ , σ ′ ), and is defined by S ( λ ) ασ,α ′ σ ′ ( q ) = h n ασ − q n α ′ σ ′ q i λ /N , (31) where n ασ − q = P k c † k + q ασ c k ασ and N is the total numberof electrons. Eq. (30) can then be written E ⊥ XC = 12 A Z dλλ X q =0 X αα ′ σσ ′ V ( λ ) ασ,α ′ σ ′ ( q ) × h N S ( λ ) ασ,α ′ σ ′ ( q ) − δ αα ′ δ σσ ′ N ασ i , (32)where N ασ is the number of electrons in the spin-resolvedband labelled by ( α, σ ), and V ( λ ) ασ,α ′ σ ′ is given by Eq. (18)after replacing − e with − e √ λ . Next we make use of thefluctuation-dissipation theorem relating the instanta-neous structure factor to the density-response function S ( λ ) ασ,α ′ σ ′ ( q ) = − ~ nπ Z ∞ χ ( λ ) ασ,α ′ σ ′ ( q , iω ) dω . (33)The density-response function, χ ( λ ) ασ,α ′ σ ′ ( q , iω ), is the pro-portionality factor between the change in electron densityat − q in the band labelled by ( α, σ ), when an arbitrar-ily weak perturbation couples to the density at q in theband labelled by ( α ′ , σ ′ ). We have analytically contin-ued the integration path to the complex frequency axisin Eq. (33) to avoid plasmon poles. We can now writethe XC energy as E ⊥ XC = 12 A Z dλλ X q =0 X αα ′ σσ ′ V ( λ ) ασ,α ′ σ ′ ( q ) × (cid:20) − ~ Aπ Z ∞ χ ( λ ) ασ,α ′ σ ′ ( q , iω ) dω − δ αα ′ δ σσ ′ N ασ (cid:21) . (34)Since the exchange energy is equivalent to first-orderperturbation theory in homogeneous electron gases, itis clear that we can obtain an alternative expres-sion for the exchange energy by replacing the fulldensity-response function, χ ( λ ) ασ,α ′ σ ′ ( q , iω ), with the non-interacting density-response function, χ (0) ασ,α ′ σ ′ ( q , iω ).Note that the non-interacting density-response functionis diagonal in the combined spin-band index ( α, σ ), i.e. χ (0) ασ,α ′ σ ′ ( q , iω ) = δ αα ′ δ σσ ′ χ (0) ασ ( q , iω ). We can subtractthe coupling-constant-integration expression for the ex-change energy from Eq. (34) to obtain the correlationenergy E ⊥ C = − ~ π Z ∞ dω Z dλλ X q =0 X αα ′ σσ ′ V ( λ ) ασ,α ′ σ ′ ( q ) × h χ ( λ ) ασ,α ′ σ ′ ( q , iω ) − χ (0) ασ,α ′ σ ′ ( q , iω ) i . (35)It is, in fact, necessary to remove the exchange-energycontribution from the coupling-constant formulation ofthe XC energy in order to obtain a convergent frequencyintegral.The density-response functions appearing in Eq. (35)form a matrix, χ ( λ ) ( q , ω ), with rows labelled by ( α, σ )and columns labelled by ( α ′ , σ ′ ). In the RPA, the density-response matrix is related to the diagonal matrix ofnon-interacting density-response functions χ (0) ( q , ω ) ≡ diag( χ (0) α ↑ ( q , ω ) , χ (0) α ↓ ( q , ω ) , χ (0) α ↑ ( q , ω ) , χ (0) α ↓ ( q , ω ) , . . . )via (cid:2) χ ( λ ) ( q , iω ) (cid:3) − = (cid:2) χ (0) ( q , iω ) (cid:3) − − V ( λ ) ( q ) , (36)where the elements of the matrix V ( λ ) ( q ) are given byEq. (18) after replacing − e with − e √ λ . Analytic ex-pressions for Re χ (0) xyσ ( q , iω ) and Im χ (0) xyσ ( q , iω ) for theisotropic 2D band case can be found, for example, inRef. 28. The subband wavefunctions, η ασ , do notchange these functions from the formulas for a per-fectly two-dimensional system. The analytic expressionsfor the elliptical-band response functions χ (0) xzσ ( q , iω )[ χ (0) yzσ ( q , iω )] can be easily identified by rescaling mo-menta, k x → k x p m L /m DOS and k y → k y p m H /m DOS [ k x → k x p m H /m DOS and k y → k y p m L /m DOS ] inEq. (8). Since this rescaling maps the elliptical bandonto a circular band, the elliptical band density-responsefunctions χ (0) xzσ ( q , iω ) and χ (0) yzσ ( q , iω ) can be written interms of χ (0) xyσ ( q , iω ): χ (0) xzσ ( q , iω ) = χ (0) xyσ ( q ′ , iω ; m DOS ) (cid:12)(cid:12)(cid:12) q ′ → √ q x ζ / + q y ζ − / (37)and χ (0) yzσ ( q , iω ) = χ (0) xyσ ( q ′ , iω ; m DOS ) (cid:12)(cid:12)(cid:12) q ′ → √ q x ζ − / + q y ζ / (38)where we have again defined ζ = m H /m L . Eq. (35) andEq. (36) combine to give a general expression for the RPAcorrelation energy which can be applied to 2DEGs withan arbitrary number of subbands. When the number ofbands is sufficiently small, Eq. (36) can be inverted an-alytically and the integral over the coupling-constant λ can be performed by hand. Furthermore, if each band isisotropic then the non-interacting density-response func-tions depend only on | q | and the integral over wavevectororientation angle is trivial.In applying the coupling-constant integration algo-rithm to the t g model for SrTiO , we assume that allelectrons in the elliptical xz and yz bands can be de-scribed by a single transverse wavefunction, η E , whileall electrons in the circular xy band can be described bya second transverse wavefunction, η C . This approxima-tion is motivated by Eq. (6), which says that electrons inthe elliptical bands hop between layers in the ˆ z -directionwith amplitude t , while electrons in the circular band hopwith amplitude t ′ . There are then only three distinct in-teractions contained in Eq. (18), i.e. V C , C = V xyσ,xyσ ′ , V E , E = V xzσ,xzσ ′ = V yzσ,yzσ ′ = V xzσ,yzσ ′ = V yzσ,xzσ ′ ,and V C , E = V E , C = V xyσ,xzσ ′ = V xyσ,yzσ ′ = V xzσ,xyσ ′ = V yzσ,xyσ ′ where we have omitted the q dependence of theinteractions for brevity. Since there are only two inde- pendent components, Eq. (36) becomes (cid:2) χ ( λ ) ( q , iω ) (cid:3) − = ( χ (0)C ) −
00 ( χ (0)E ) − ! − V ( λ )C , C V ( λ )C , E V ( λ )E , C V ( λ )E , E ! . (39)Here, χ (0)C ( q , iω ) is the density response function of thecircular bands, χ (0)C ( q , iω ) = X σ χ (0) xyσ ( q , iω ) , (40)and χ (0)E ( q , iω ) is the density response function of theelliptical bands, χ (0)E ( q , iω ) = X σ X α = xz,yz χ (0) ασ ( q , iω ) . (41)Finally, as in the calculation of the exchange energypresented above, we consider the thin-film limit in whichthe Coulomb-interaction entering the RPA correlationenergy loses its orbital dependence, and we can replace V ασ,α ′ σ ′ ( q ) with v q . Eq. (39) simplifies to a single-component expression1 χ ( λ ) ( q , ω ) = 1 χ (0) ( q , ω ) − v ( λ ) q (42)where we define χ (0) ( q , ω ) = P ασ χ (0) ασ ( q , ω ). In thissame limit Eq. (35) reduces to E ⊥ C = − ~ π Z ∞ dω Z dλλ X q =0 v q × h χ ( λ ) ( q , iω ) − χ (0) ( q , iω ) i , (43)which after carrying out the integration over the couplingconstant gives E ⊥ C = − ~ π Z ∞ dω X q =0 × h v q χ (0) ( q , iω ) + ln (cid:16) − v q χ (0) ( q , iω ) (cid:17)i . (44)As explained in more detail in Appendix A, we findthat the (always negative) correlation energy per elec-tron is decreased in magnitude by increasing the bandanisotropy. Furthermore, the correlation energy islargest in magnitude when the electron density is dis-tributed equally amongst all of the spin-resolved bands.This property implies that the correlation energy worksagainst the type of spontaneous symmetry breaking fa-vored by exchange. We discuss this competition in Sec-tion IV, but first must account for the contribution tothe ground-state energy from the q = 0 electron-electroninteraction term ( i.e. the Hartree potential), the elec-trostatic external potential, and from hopping betweenlayers.In this section we have obtained expressions for thecontribution to the total ground-state energy from termsin the t g E ⊥ = E ⊥ K + E ⊥ X + E ⊥ C . (45)In the next section we will give expressions for the to-tal ground-state energy of the t g N l a p n/g ≪
1, nor our restriction to only two distincttransverse wavefunctions, ( η C , η E ), are critical to the ap-plicability of the variational theory. B. Total Energy
Next we calculate the contribution to the total ground-state energy from terms in the t g H k = T k + V k + V k ext . This is obtained by taking the expectation value of H k in the many-body state defined by h Ψ { n ασ , η ασ } |H ⊥ | Ψ { n ασ , η ασ } i = E ⊥ . (46)Although E ⊥ was not calculated in the previous sectionby appealing to an explicit many-body wavefunction, theappropriate N-electron wavefunction can always be writ-ten as a linear combination of Slater determinants ofsingle-particle states | ψ k ασ i : | Ψ { n ασ , η ασ } i = X i c i n X P ( − P ˆP (cid:2) N Y j =1 | ψ k j α j σ j i (cid:3)o , (47)where ˆP is the permutation operator, P is the number ofpermutations in each term, and we are summing over allpermutations of the single-particle state labels amongstthe N electrons. Here, c i is a complex coefficient andthe index i runs over all Slater determinants with a fixednumber of electrons in each spin-resolved band and fixedtransverse wavefuctions; we choose the sets { n ασ } and { η ασ } as good quantum numbers of our many-body vari-ational wavefunction | Ψ { n ασ , η ασ } i .The first contribution to E k is from interlayer hoppingdescribed by Eq. (6). Taking the expectation value of T k we find that E k K = A X ασ n ασ (cid:0) η ασ † · t α · η ασ (cid:1) , (48)where t xz = t yz = − t δ R z ,R z ± a + 2( t − t ′ ) δ R z ,R ′ z and t xy = − t ′ δ R z ,R z ± a are N l × N l matrices in layer index.The diagonal contribution 2( t − t ′ ) to t xz and t yz ac-counts for the energy offset between the band edge of the elliptical xz and yz bands and the circular xy band whenconfined to a single TiO layer (see Fig. 1). This expres-sion is independent of the correlations among transversedegrees of freedom, approximated above using the RPA.Similarly, the Hartree energy can be obtained by tak-ing the expectation value of the portion of the Coulombinteraction sensitive to electron positions along the con-finement direction, V k : E k H = A X ασ n ασ (cid:0) η ασ † · V H · η ασ (cid:1) = − A πe κ X R z ,R ′ z n R z | R z − R ′ z | n R ′ z , (49)where n R z is the total electron density in layer R z . Thematrix V H entering Eq. (49) is diagonal in TiO layerindex, with elements given by (cid:2) V H (cid:3) R z ,R ′ z = − δ R z ,R ′ z (cid:18) πe κ (cid:19) × X R ′′ z | R z − R ′′ z | X ασ n ασ | η ασR ′′ z | . (50)For definiteness, we will assume in this article the com-mon case in which there is only one electrostatic gate(which we place at z = 0) with total charge + e n =+ e P ασ n ασ (see Figure 1). The electronic energy fromthe external gate’s confinement potential is then E k ext = A X ασ n ασ (cid:0) η ασ † · V ext · η ασ (cid:1) , (51)where V ext is also a diagonal matrix in TiO layer index,with matrix elements given by (cid:2) V ext (cid:3) R z ,R ′ z = (cid:18) πe nR z κ (cid:19) δ R z ,R ′ z . (52)Combining terms, we write the contribution from theparallel part of the Hamiltonian to the total ground-stateenergy as E k = E k K + E k H + E k ext . (53)Notice that since each term in H ⊥ depends only ona small number of macroscopic observables ( e.g. thedensity in each layer and/or the density in each spin-resolved band), which are completely determined fromthe good quantum numbers of our variational wavefunc-tion Eq. (47), each terms contribution to the ground-state energy is given exactly by the corresponding classi-cal mean-field energy.When E k is combined with the contributions to theground-state energy from the perpendicular part of theHamiltonian, we obtain the total ground-state energy asa functional of the density in each spin-resolved band andthe transverse subband wavefunction spinors: E (cid:2) { n ασ , η ασ } (cid:3) = E k + E ⊥ . (54) C. Minimizing the total energy
To solve for the ground-state values of { n ασ } and { η ασ } , we minimize the energy functional in Eq. (54)subject to the constraint equations X ασ n ασ = n X R z | η ασR z | = 1 . (55)The second equation is applied for each α and σ . Forthese constraints we introduce Lagrange multipliers, µ and { ν ασ } , respectively. Minimization of the total energyfunctional with respect to each particular n ασ yields anequation µ ⊥ ασ + η ασ † · (cid:0) t α + V H + V ext (cid:1) · η ασ = µ , (56) where we have defined µ ⊥ ασ = A − dE ⊥ /dn ασ . In obtain-ing the results presented in Sect. IV, we have evaluated µ ⊥ ασ numerically by combining the above expressions for E ⊥ with finite difference formulas. Minimization of thetotal energy functional with respect to each componentof the spinor η ασ † yields an eigenvalue equation (cid:0) X + Y ασ + t ασ + V H + V ext (cid:1) · η ασ = (cid:18) ν ασ A n ασ (cid:19) η ασ (57)for each value of α and σ . Here, the diagonal matrices X and Y ασ are found from differentiating the correla-tion energy and exchange energy, respectively, and theirmatrix elements are given by (cid:2) X (cid:3) R z ,R ′ z = δ R z ,R ′ z (cid:18) − ~ π (cid:19) Z ∞ dω Z dλλ X q =0 X αα ′ σσ ′ × M ασ,α ′ σ ′ ( R z , q ) n χ ( λ ) ασ,α ′ σ ′ − χ (0) ασ,α ′ σ ′ o + V ( λ ) ασ,α ′ σ ′ ( q ) X ββ ′ ττ ′ χ ( λ ) ασ,βτ M βτ,β ′ τ ′ ( R z , q ) χ ( λ ) β ′ τ ′ ,α ′ σ ′ (58)and (cid:2) Y ασ (cid:3) R z ,R ′ z = δ R z ,R ′ z (cid:18) − A π (cid:19) Z ∞ q dq Z π dφ M ασ,ασ ( R z , q F [ ζ, φ ]) Θ (cid:16) √ πn ασ − q (cid:17) × " n ασ − q π r πn ασ − (cid:16) q (cid:17) − π arcsin (cid:18) q √ πn ασ (cid:19) , (59)where we have defined M ασ,α ′ σ ′ ( R z , q ) = X R ′′ z πe e − q | R z − R ′′ z | κq (cid:16) | η ασR ′′ z | + | η α ′ σ ′ R ′′ z | (cid:17) . (60)In Eq. (58) we have ommitted the ( q , iω ) dependence ofthe density response functions for brevity. Eqs. (55)-(57)provide enough independent linear relationships to solvesolve self-consistently for the various densities, { n ασ } ,and transverse wavefunction spinors { η ασ } , introducedin our model.As already mentioned above, the numerical results wepresent in Sect. IV have been obtained for thin films ofSrTiO , which we identify by the condition N l a p n/g ≪
1. While this approximation is not essential to our vari-ational approach, by making it we simplify the numeri-cal problem to be solved in this first application of our method in which we study orbital and spin ordering insta-bilities. Specifically, in the thin-film limit the Coulombinteraction, Eq. (18), loses its orbital dependence and theexchange and correlation energies are no longer directlydependent on the transverse wavefunction spinors (seeEq. (27) and Eq. (44)). As a result the matrices X and Y ασ vanish. After also making the approximation thatonly two distinct transverse wavefunction spinors exist,0 η C and η E , the constraints in Eq. (55) can be written as X β =E , C n β = n X R z | η C R z | = 1 X R z | η E R z | = 1 (61)and Eq. (56) and Eq. (57) become µ ⊥ β + η β · (cid:0) t β + V H + V ext (cid:1) · η β = µ (62)and (cid:0) t β + V H + V ext (cid:1) · η β = (cid:18) ν β A n β (cid:19) η β , (63)where β = E or C. Eq. 62 guarantees that all occupiedsubbands have the same chemical potential, while Eq. 63chooses the subband wavefunction as the eigenvalue ofthe mean-field Hamiltonian. Appropriate expressions formatrices t β and V H , as well as for µ ⊥ β follow in an obviousmanner from the more general definitions in Sect. III B.We note in passing that although the exchange and cor-relation energies in this approximation are no longer di-rectly dependent on the transverse wavefunction spinors,they are still sensitive to confinement effects. For exam-ple, the spinors, η C and η E are determined by solvingEq. (63) and are therefore sensitive to the band-offsetbetween the xz and yz bands and the xy bands. Theband-offset obviously influences the equilibrium valuesof n β , which in turn affects the values of the exchangeand correlation energies defined in Eq. (27) and Eq. (44).We now outline one method for obtaining the self-consistent solution. The first step in the determinationof the densities ( n E , n C ), and wavefunctions (cid:0) η E , η C (cid:1) ,is to make an initial choice for them which satisfies theequations of constraint. In step two, use these to cal-culate V H , µ ⊥ E and µ ⊥ C . Eq. (62) with β = C can becombined with Eq. (62) with β = E to eliminate µ , androot-finding numerical algorithms can be applied to solvefor new values of n E and n C . These values are then usedto re-evaluate the Hartree potential V H . Eq. (63) is aneigenvalue problem with a real, symmetric, and tridiag-onal matrix. It follows that the eigenvalues are real andthe eigenvectors can be chosen to have only real com-ponents, i.e. (cid:0) η β (cid:1) † = (cid:0) η β (cid:1) T . We solve this eigenvalueproblem numerically and select the eigenvectors, η E and η C , which correspond to the smallest eigenvalues. Theseeigenvalues corresponds to the band-edge energies of thesubbands in our model. We keep the wavefunctions cor-responding to the smallest energy eigenvalues becausewe are only interested in keeping the lowest subband ofeach orbital type xy , xz and yz . We next return to theaforementioned “step two”, and repeat until convergenceis obtained. Adaptable numerical algorithms exist in theliterature to assist in solving the eigenvalue problem, finding spline interpolations for the correlation energy,finding roots, and evaluating multi-dimensional integrals.After a self-consistent solution has been obtained, theeigenvalues of Eq. (63) give an approximation for theband-edge energies of the bands retained in the t g χ ( λ ) ( q , iω ). Whilebeing exact at high electron density, the RPA neglectsvertex corrections which become more important at verylow density. Various methods exist for including vertex-correction physics in the matrix of density-response func-tions (see Ref. 28 and references therein), and we proposethat in future calculations the relative success of these ap-proximations can be gauged by comparing the predictedband-edge energy offsets against those measured directlyin angle-resolved photoemission experiments. IV. SPIN AND ORBITAL ORDER
In this Section we present the results of our variationalcalculation for the ground-state energy as a function of n E , n C , η E , and η C . Specifically, we consider three typesof fully broken symmetry states characterized by the dis-tribution of electrons in the four spin-resolved ellipticalbands: spin-ordered ( e.g. n xz ↑ = n yz ↑ = n E / n xz ↓ = n yz ↓ = 0), orbital-ordered ( n xz ↑ = n xz ↓ = n E / n yz ↑ = n yz ↓ = 0), and orbital- and spin-ordered( n xz ↑ = n E and n xz ↓ = n yz ↑ = n yz ↓ = 0). In allof these states we keep the density in the isotropic xy band equally distributed amongst the two spin-species, n xy ↑ = n xy ↓ = n C /
2. The distribution of electrons inthese spontaneously broken symmetry states contrastswith the distribution in the normal phase, where all fourspin-resolved elliptical bands are either equally occupied(above some critical density which depends on the thick-ness of SrTiO ), or all empty. Each broken symmetrystate has is own set of n E , n C , η E , and η C , which wevary to find the ground-state energy of that phase. Byselecting the state with the lowest energy at each valueof total 2DEG density, n , and number of TiO layers, N l ,we have have constructed the phase diagrams shown inFig. 2.Examining Fig. 2a), we find that spontaneous orbital-order occurs first, and is joined by an accompanying spin-order as the density is decreased. Several striking fea-tures appear here. First, the lower critical density forsymmetry breaking decreases as the thickness of SrTiO increases. Rather than being an interaction effect, thisis caused by a reduction in the intrinsic band-edge off-1 n (cm -2 )n (cm -2 )12345678910N l l a)b) FIG. 2: (Color online) Phase diagram for the t g plotted for varying total density, n , and number ofTiO layers, N l . In all phases shown there is a finite electrondensity in the isotropic xy band which is equally distributedbetween spins. The colors distinguish phases with differentdistributions of electrons amongst the four spin-resolved el-liptical bands ( xz, ↑ , xz, ↓ , yz, ↑ , and yz, ↓ ). Purple identifiesthe normal low-density state in which no electrons occupy theelliptical bands. Green identifies a orbital- and spin-orderedstate in which only one of the spin-resolved elliptical bandscontains electrons (e.g. xz, ↑ ). Light blue identifies an orbital-ordered state in which only one elliptical band is occupied(e.g. xz, ↑ and xz, ↓ ). Yellow identifies a spin-ordered statein which both elliptical bands have electrons but only of onespin species (e.g. xz, ↑ and yz, ↑ ). Red identifies the nor-mal high-density state in which the density in all four spin-resolved elliptical bands is equal. In panel a) correlationsare included within the RPA. In panel b) correlations are ne-glected and only the exchange energy is included. Althoughthe orbital-ordered state is degenerate with the spin-orderedstate in b), correlations prefer an orbital-ordered state andthe spin-ordered state is absent in a). set 2( t − t ′ ), see e.g. Fig. 1. Specifically, the elliptical xz and yz bands have a large hopping amplitude par-allel to confinement, t , and the isotropic xy band has a small hopping amplitude, t ′ . Therefore the xy band ismore tightly localized to the interface than the xz and yz bands (a result that is quantitatively contained in thedifferences between η C and η E , respectively), and thedifference in band-edge energies decreases relative to thesingle-layer value. This mechanism becomes more pro-nounced as the sample thickness is increased, leading toa smaller effective band-edge offset at each given totaldensity. Accurately capturing the change in band-edgeoffset as the total density is varied, and perhaps moreimportantly as density is allowed to shift between the el-liptic bands to the circular bands near an instability, iskey to finding the correct phase diagram. In the next Sec-tion we elaborate on the role of correlations in selectingthe most energetically favorable broken symmetry state. A. Plasmon-pole approximation: Correlationsprefer orbital order
In the isotropic single-band 2DEG, the competitionbetween kinetic and exchange energies qualitatively cap-tures the transition from the normal to the ferromagneticstate. The role of correlations is to lower the critical den-sity for symmetry breaking . Surprisingly, in our modelfor thin-film SrTiO heterojunctions, we find that thecorrelation energy plays a more qualititave role. Morespecifically, at the exchange-only level the spin-orderedstate is energetically degenerate with the orbital-orderedstate, as shown by the presence of the diagonally-hatchedblue and yellow region in Figure 2b). The effectiveexchange-interaction only acts between electrons in thesame spin-resolved band and therefore the ground-stateenergy is insensitive to whether the xz, ↑ and xz, ↓ bandsare occupied, as in a orbital-ordered state, or whether the xz, ↑ and yz, ↑ bands are occupied, as in a spin-orderedstate; the total exchange energy is given in Eq. (23), andis simply the sum of the exchange energy from each sepa-rate spin-resolved band. Once correlations are included,however, our calculations reveal that the system prefersorbital-order to spin-order. Indeed, the orbitally orderedphase (blue) remains in Fig. 2a), while the spin orderedphase (yellow) is absent.This result can be understood by applying the plasmon-pole approximation (PPA) to the RPA correla-tion energy. For simplicity we consider a model withoutthe xy band, and only consider the four remaining bands( i.e. α = xz, yz and σ = ↑ , ↓ ) of our model. The plasmonfrequency of the spin polarized (SP) and orbitally polar-ized (OP) states is given by the zeroes of their respectivedielectric functions ǫ ( q , Ω i ) = 1 − v q ˜ χ i ( q , Ω i ) = 0 , (64)where i = OP or SP, and ˜ χ i ( q , ω ) is the proper densityresponse function of the i ’th broken symmetry state. Inthe RPA the proper density response function is given bysumming the non-interacting density response function2 Spin orderOrbital order0.0 0.2 0.4 0.6 0.8 1.0q/k
FDOS -0.48-0.46-0.44-0.42-0.40-0.38-0.36 E c o rr ( q ) / N ( e V A ) . . . FIG. 3: (Color online) The contribution to the correlation en-ergy per electron from transitions with wavevector q is plottedfor the orbital-ordered (solid brown) and spin-ordered (dashedgreen) broken symmetry states. (Compare to Eqs. (74) andEq. (75) respectively.) The correlation energy is calculatedwithin the plasmon-pole approximation and in a simplifiedmodel in which the circular xy band is absent. We have takenthe total density to be 5 × cm − . of each occupied spin-resolved band. For the orbital-ordered state˜ χ OP ( q , ω ) = χ (0) xz, ↑ ( q , ω ) + χ (0) xz, ↓ ( q , ω ) , (65)and for the spin-polarized state˜ χ SP ( q , ω ) = χ (0) xz, ↑ ( q , ω ) + χ (0) yz, ↑ ( q , ω ) . (66)To make analytic progress we will consider plasmons inthe long-wavelength limit (i.e. q → ω > v F q ). Thenon-interacting density response function of an isotropicband (e.g the xy band we are currently neglecting) in thislimit is lim q → χ (0) xyσ ( q, ω ) = n xyσ q m L ω . (67)Expressions for the long-wavelength non-interacting den-sity response functions of the elliptical xz and yz bandscan be found from Eq. (37) and Eq. (38). After insert-ing the relevant proper density response function intoEq. (64), we find the plasmon frequency of the spin-ordered and orbital-ordered state isΩ SP = r v q n m q (68)and Ω OP = r v q nm DOS q ζ q x + ζ − q y , (69)respectively. Here, n is the total density of the xz and yz bands only since we are not considering the xy band. Thedensity-of-states mass is m DOS = √ m L m H , the reduced mass is ¯ m = m L m H / ( m L + m H ). Interestingly, the SPplasmon becomes isotropic in the long-wavelength limit.In the PPA, the imaginary part of the density responsefunction is expanded about the plasmon-pole of the sys-tem . Using the long-wavelength form of the densityresponse function we obtainIm χ PPA i ( q, ω ) = − π Ω i v q δ ( ω − Ω i ( q )) . (70)The fluctuation-dissipation theorem allows us to thenfind the instantaneous structure factor S PPA i ( q ) = − ~ πn Z ∞ Im χ PPA i ( q , ω ) dω (71)and then the interaction energy per electron in the i ’thbroken symmetry state is given by E i C N = 18 π Z d q v q (cid:2) S PPA i ( q ) − (cid:3) . (72)Although the long-wavelength PPA static structure fac-tor gives a formally divergent interaction energy, the con-tribution to the interaction energy at each wavevector | q | has meaning in the long-wavelength limit | q | ≪ k F45 . Wedefine the correlation energy per electron from transitionsat wavevector | q | , E i C ( q ), by the relationship E i C N = Z ∞ E i C ( q ) N dq. (73)We find that the contribution to the interaction energyper electron at each wavevector | q | is E OPC ( q ) N = ~ q π n r v q nm DOS (cid:20) ζ / Ξ (cid:18) ζ − ζ (cid:19)(cid:21) − qv q π (74)in the orbital-ordered state, and in the spin-ordered stateis given by E SPC ( q ) N = ~ q πn r v q n m − qv q π , (75)where Ξ( x ) is the complete elliptic integral of the sec-ond kind. The last term on the right hand side in bothEq. (74) and (75) is proportional to the classical self-interaction energy of the plasmon density fluctuation atwavelength 2 π/q . The first term on the right hand sideis the zero-point-energy of each plasmon mode .In Fig. 3 we plot Eq. (75) and Eq. (74) against | q | inunits of an effective density-of-states Fermi wavevectordefined by k FDOS = p πn/g , where g = 2 is the num-ber of occupied spin-resolved bands in either the spin-or orbital-ordered states. The orbitally polarized phaseis of lower energy as a result of the smaller zero-pointplasmon energy at wavevector | q | .3 V. SUMMARY AND DISCUSSION
Orbital order is an interesting possibility in t g electronsystems, even in the itinerant limit of doped d materi-als with a small value for the number of electrons permetal site. We refer to the low-density t g systems as t g electron gases. The orbitally ordered t g electron gas isan electron nematic, a translationally invariant metallicphase with a spontaneously generated spatial anisotropy.In this article we have developed and applied a theoryof a two-dimensional t g electron gases confined at a sur-face or interface by an external electric field, which ac-counts for transverse correlations induced by Coulombicelectron-electron interactions using a random phase ap-proximation, and optimizes the subband wavefunctionsof the three t g two-dimensional bands by minimizing to-tal energy. In the present calculation we have assumedthat the t g t g electron gases, for example the difference betweenthe exchange and correlation energies of electrons in el-liptical vs. circular bands.In the ideal isotropic single-band electron gas, the RPAapproximation that we employ for transverse correlationsis quantitatively unsuccessful in predicting the Blochferromagnetic instability, drastically overestimatingthe density at which it occurs. However the RPA canbe viewed as the leading term in an expansion in N − f ,where N f is the number of fermion bands, and we canexpect that it becomes more accurate when the numberof Fermion species increases. This observation suggeststhat the RPA is a reasonable starting point for treatingCoulomb correlations in the six band t g t g electrons areconfined to an interface and also confined in a quantumwell with a finite number of metal layers, N l . Below acritical density that decreases with increasing N l only xy bands are occupied because they have smaller confine-ment energies. We find that orbital order occurs whenthe yz, xz elliptical band density is below a critical value,and that it occurs over a broader range of total density forwider quantum wells. At very low elliptical band densityspin order occurs in additional to orbital order.Anisotropic magnetoresistance measurements ofLaAlO /SrTiO heterojunctions indicate uniaxialanisotropy which spontaneously appears above a criticaldensity at which the xz and yz bands start to be occu-pied, and disappears when the xz and yz bands’ densitybecomes large. The magnetic properties of these systemsis not perfectly understood to date, but some theoreticalanalysis based on Hubbard-like model interactions andKondo models has recently appeared. We point outthat both the orbitally ordered as well as the orbitallyand spin ordered phases we have found would introducea uniaxial anisotropy in the system, and may explain these measurements.Several other systems have shown indications of theBloch-like instability which drives the transition to or-bital and spin order in the t g electron gas. Experi-mentally, some evidence for a ferromagnetic instabilityhas been observed in transport studies of very clean low-density Si 2DEGs . Spin-ordering in these itinerantelectron 2DEGs is supported by a large enhancement inthe spin-susceptibility as the density is lowered . Mean-while, in the multivalley semiconductor AlAs, where thevalley degree-of-freedom is closely analogous to our t g . In addition, exchange-correlation en-ergies are understood to play an important quantitativerole in valley occupancy symmetry breaking in Si in-version layers . Lastly, numerical studies of single-band 2DEG systems which include both disorder andthe Coulomb interaction have found a magnetic orderedground-state . Despite these studies, Bloch-like bro-ken symmetry states have not been unambiguously ob-served to date. Our calculations suggest that Bloch-likeinstability to an orbital-ordered state may be present inthin-film SrTiO heterojunctions.The variational approach described in this Articlecan be applied to wide semiconductor quantum wells inwhich there is also a subtle interplay between correlationsand subband wavefunctions, and to Van der Waals het-erostructures . Although multiple orbitals are not oftenpresent near the Fermi energy, these layered materialshave a weak electronic hopping amplitude in the confine-ment direction, and as a result many subbands are oc-cupied even in relatively thin multilayers. Additionally,many van der Waals materials have been observed to bestrongly interacting, often exhibiting diverse phenomenalike charge-density wave instabilities , superconductiv-ity , and negative compressibility. Acknowledgments
JRT thanks Jeremy Levy, Lu Li, Jos´e Lorenzana andFengcheng Wu for interesting discussions. Work inAustin was supported by the DOE Division of Mate-rials Sciences and Engineering under grant DE-FG02-ER45118 and by the Welch foundation under grantTBF1473. MP was supported by Fondazione IstitutoItaliano di Tecnologia and the European Union’s Horizon2020 research and innovation programme under GrantAgreement No. 696656 “GrapheneCore1”.
Appendix A: The RPA correlation energydependence on number of bands and bandanisotropy
In this Appendix we explore the dependence of theRPA correlation energy on the number of bands with4distinct orbital character, and on the degree of bandanisotropy. In Figure 4 we examine how the correla-tion energy changes by the addition of a second iden-tical isotropic band of effective mass m eff . We plot thecorrelation energy in units of effective Rydberg Ry . = m eff ry ./κ where ry . = 13 . r s in this figure is defined in terms ofthe total 2DEG density, n , which is distributed equallyamongst either one or two bands, r s = m eff a B κ √ πn . (A1)Since the dielectric constant κ and the effective mass m eff enter into the definition of both Ry . and r s , the curves inFigure 4 are independent of the exact numerical value ofboth κ and m eff . We find that the RPA correlation en- ergy per electron is larger for the system with two bands.In Figure 5 we examine how the correlation energy ofa single 2D band changes as the degree of anisotropyvaries. We take the band to have an elliptical Fermisurface and a heavy mass, m H , in the x -direction anda light mass, m L , in the y -direction (exactly like the yz band of our model for SrTiO in the main text). Weplot the correlation energy in units of effective RydbergRy . = √ m H m L ry ./κ . The parameter r s in this figure isdefined in terms of the total 2DEG density n , r s = √ m H m L a B κ √ πn . (A2)We find that the correlation energy per electron is re-duced as the degree of band anisotropy increases. ∗ Electronic address: [email protected] J. B. Goodenough, Phys. Rev. , 564 (1955). E. O. Wollan and W. C. Koehler, Phys. Rev. , 545(1955). J.-G. Cheng, Y. Sui, J.-S. Zhou, J. B. Goodenough and W.Su, Phys. Rev. Lett. , 087205 (2008). J.-Q. Yan, J.-S. Zhou and J. B. Goodenough, Phys. Rev.Lett. , 235901 (2004). J.-S. Zhou, J. B. Goodenough, J.-Q. Yan, J.-G. Cheng, K.Matsubayashi, Y. Uwatoko and Y. Ren Phys. Rev. B. ,224422 (2009). T. Mizokawa and A. Fujimori, Phys. Rev. B. , 5368(1996). M. Mochizuki and M. Imada, J. Phys. Soc. Jpn. , 1833(2004). J. Mannhart and D. G. Schlom, Science , 1607 (2010). A. P. Kajdos, D. G. Ouellette, T. A. Cain, and S. Stemmer,Appl. Phys. Lett. , 082120 (2013). Y. Kozuka, M. Kim, H. Ohta, Y. Hikita, C. Bell, and H.Y. Hwang, Appl. Phys. Lett. , 222115 (2010). S. Stemmer and S. J. Allen, Annu. Rev. Mater. Res. ,151 (2014). J. A. Sulpizio, S. Ilani, P. Irvin, and J. Levy, Annu. Rev.Mater. Res. , 117 (2014). J. R. Tolsma, A. Principi, R. Asgari, M. Polini, and A. H.MacDonald, Phys. Rev. B. , 045120 (2016). F. Bloch, Z. Physik , 545 (1929). D. P. Young, D. Hall, M. E. Torelli, Z. Fisk, J. L. Sarrao,J. D. Thompson, H.-R. Ott, S. B. Oseroff, R. G. Goodrich,and R. Zysler, Nature , 412 (1999). S. A. Vitakalov, H. Zheng, K. M. Mertes, M. P. Sarachik,and T. M. Klapwijk, Phys. Rev. Lett. , 086401 (2001). A. A. Shashkin, S. V. Kravchenko, V. T. Dolgopolov, andT. M. Klapwijk, Phys. Rev. Lett. , 086801 (2001). N. Teneh, A. Yu. Kuntsevich, V. M. Pudalov, and M.Reznikov, Phys. Rev. Lett. , 226403 (2012). K. Takashina, Y. Niida, V. T. Renard, B. A. Piot, D. S.D. Tregurtha, A. Fujiwara, and Y. Hirayama, Phys. Rev.B , 201301(R) (2013). V. T. Renard, B. A. Piot, X. Waintal, G. Fleury, D.Cooper, Y. Niida, D. Tregurtha, A. Fujiwara, Y. Hirayama, FIG. 4: (Color online) A plot of the the RPA correlationenergy per electron (in units of effective Rydbergs) as theparameter r s is varied from one to forty. Eq. (A1) defines r s in terms of the total electron density. The blue curve isfor a single band isotropic 2DEG with equal populations ofboth spin species. The red curve is for a 2DEG with twoidentical isotropic bands with equal spin populations in bothbands. Note that the bands are of different orbital content,and therefore interband scattering is negligeable. The cor-relation energy per electron increases in magnitude with thenumber of occupied bands.and K. Takashina, Nature Commun. , 7230 (2015). Y. Zhang and S. Das Sarma, Phys. Rev. B. , 115317(2005). C. Attaccalite, S. Moroni, P. Gori-Giorgi and G.B.Bachelet, Phys. Rev. Lett. , 256601 (2002). B. Tanatar and D. Ceperley, Phys. Rev. B. , 5005 (1989). P. Hohenberg and W. Kohn, Phys. Rev. , B864 (1964). W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965). D. C. Langreth and J. P. Perdew, Phys. Rev. B. , 5469(1980). FIG. 5: (Color online) A plot of the RPA correlation energyper electron (in units of effective Rydberg) for single band2DEGs of varying band anisotropy defined by the ratio of themass in the x -direction, m H , and the mass in the y -direction, m L . The correlation energy per electron is reduced in magni-tude for increasing band anisotropy. J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996). G. F. Giuliani and G. Vignale,
Quantum Theory of theElectron Liquid (Cambridge University Press, Cambridge,2005). R. O. Jones and O. Gunnarsson, Rev. Mod. Phys. , 689(1989). Z. S. Popovic, S. Satpathy, and R. M. Martin, Phys. Rev.Lett. , 256801 (2005). G. Khalsa and A.H. MacDonald, Phys. Rev. B. , 125121(2012). S.E. Park and A.J. Millis, Phys. Rev. B. , 205145 (2013). K. Ueno, S. Nakamura, H. Shimotani, A. Ohtomo, N.Kimura, T. Nojima, H. Aoki, Y. Iwasa, and M. Kawasaki,Nature Mater. , 855 (2008). A. Joshua, S. Pecker, J. Ruhman, E. Altman, and S. Ilani,Nature Commun. , 1129 (2012). Santosh Raghavan, S. James Allen, and Susanne Stemmer,Appl. Phys. Lett. , 212103 (2013). Y. Lee, C. Clement, J. Hellerstedt, J. Kinney, L. Kinnis-chtzke, X. Leng, S. D. Snyder, and A. M. Goldman, Phys.Rev. Lett. , 136809 (2011). P. Moetakef, C. A. Jackson, J. Hwang, L. Balents, S.J. Allen, and S. Stemmer, Phys. Rev. B. , 201102(R)(2012). R. Chen, S. Lee, and L. Balents, Phys. Rev. B. ,161119(R) (2013). S. J. Allen, B. Jalan, S. Lee, D. G. Ouellette, G. Khalsa,J. Jaroszynski, S. Stemmer, and A. H. MacDonald, Phys.Rev. B. , 045114 (2013). Although the long-wavelength low-temperature static di-electric constant of bulk SrTiO is in the tens of thou-sands [A. S. Barker Jr. and M. Tinkham, Phys. Rev. ,1527 (1962)] because of the presence of a soft LO phononmode near the Γ point, the effective dielectric constantwhich screens electron-electron interactions is expected to be substantially smaller. The electronic transitions whichcontribute to the ground-state energy in the t g model areon the scale of a few hundred meV, much larger than thesoft-phonon mode energy which is closer to a few meV [Y.Yamada and G. Shirane, J. Phys. Soc. Jpn. , 396 (1969)].Screening by this mode is therefore weak. Furthermore,the LO phonon mode is soft only in close vicinity to the Γpoint. Even static electric fields are ineffectively screenedby this mode unless they are constant over distances whichgreatly exceed a lattice constant. Based on these consid-erations, we think that the effective dielectric constant tobe included in electron-electron interaction calculations iscloser to κ ∼
15, similar to small-gap covalent semiconduc-tors. At liquid helium temperatures the dielectric constantof bulk SrTiO is κ ∼
10 in the range 2 −
40 meV and κ ∼ , 2124 (1972)]. I. S. Gradshteyn and I. M. Ryzhik,
Table of Integrals, Se-ries, and Products (Elsevier Academic Press, 2007). L. Zheng and A. H. MacDonald, Phys. Rev. B. , 5522(1994). E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel,J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling,A. McKenney, and D. Sorensen,
LAPACK Users’ Guide (Society for Industrial and Applied Mathematics, Philadel-phia, 1999). W. Press, S. A. Teukolsky, W. T. Vetterling, and B. P.Flannery,
Numerical Recipes in C (Cambridge UniversityPress, Cambridge, 2002). D. Pines and P. Nozi´eres,
The Theory of Quantum Liquids (W.A. Benjamin, Inc., New York, 1966). E. Fradkin, S. A. Kivelson, M. J. Lawler, J. P. Eisensteinand A. P. Mackenzie, Annu. Rev. Condens. Matter Phys. ,153 (2010). M. Ben Shalom, C. W. Tai, Y. Lereah, M. Sachs, E. Levy,D. Rakhmilevitch, A. Palevski and Y. Dagan, Phys. Rev.B. , 140403(R) (2009). A. Annadi, Z. Huang, K. Gopinadhan, X. Renshaw Wang,A. Srivastava, Z. Q. Liu, H. Harsan Ma, T. P. Sarkar,T. Ventatesan and Ariando, Phys. Rev. B. , 201102(R)(2012). A. Joshua, J. Ruhman, S. Pecker, E. Altman and S. Ilani,Proc. Nat. Acad. Sci. , 9633 (2013). S. Banerjee, Onur Erten, and M. Randeria, Nature Phys. ,626 (2013). J. Ruhman, A. Joshua, S. Ilani, and E. Altman, Phys. Rev.B. , 125123 (2014). V. M. Pudalov, M. E. Gershenson, H. Kojima, N. Butch, E.M. Dizhur, G. Brunthaler, A. Prinz, and G. Bauer, Phys.Rev. Lett. , 196404 (2002). O. Gunawan, Y. P. Shkolnikov, K. Vakili, T. Gokmen, E.P. De Pootere, and M. S. Shayegan, Phys. Rev. Lett. ,186404 (2006). W. L. Bloss, S. C. Ying, and J. J. Quinn, Phys. Rev. B. ,1839 (1980). W. L. Bloss, L. J. Sham, and V. Vinter, Phys. Rev.Lett. , 1529 (1979). G. Benenti, G. Caldara, and D. L. Shepelyansky, Phys.Rev. Lett. , 5333 (2001). A. K. Geim, and I. V. Grigorieva, Nature , 419 (2013). X. Xi, L. Zhao, Z. Wang, H. Berger, L. Forro, J. Shan, andK. F. Mak Nature Nanotech. , 765 (2015). J. M. Lu, O. Zheliuk, I. Leemakers, N. F. Q. Yuan, U.Zeitler, K. T. Law, and J. T. Ye, Science , 1353 (2015). S. Larentis, J. R. Tolsma, B. Fallahazad, D. C. Dillen, K.Kim, A. H. MacDonald, and E. Tutuc, Nano. Lett.4