Giant oscillating thermopower at oxide interfaces

Understanding the nature of charge carriers at the LaAlO3/SrTiO3 interface is one of the major open issues in the full comprehension of the charge confinement phenomenon in oxide heterostructures. Here, we investigate thermopower to study the electronic structure in LaAlO3/SrTiO3 at low temperature as a function of gate field. In particular, under large negative gate voltage, corresponding to the strongly depleted charge density regime, thermopower displays high negative values of the order of 104–105μVK−1, oscillating at regular intervals as a function of the gate voltage. The huge thermopower magnitude can be attributed to the phonon-drag contribution, while the oscillations map the progressive depletion and the Fermi level descent across a dense array of localized states lying at the bottom of the Ti 3d conduction band. This study provides direct evidence of a localized Anderson tail in the two-dimensional electron liquid at the LaAlO3/SrTiO3 interface.

The simulated resistance shows oscillations in correspondence of the Fermi level crossing the localized levels, as does the Seebeck curve. However, these oscillations are not observed experimentally. We first point out that, from the experimental point of view, we carried out the resistance measurements with the highest possible level of accuracy and reliability, trying different configurations. In particular, we measured several samples in the maximum possible range of V g , compatible with input impedance of nanovoltmeters, R sheet measurements were carried out using both d.c. (either applying alternating polarity or fixed polarity) and a.c. (lock in) techniques, we measured both voltage drops with a current bias and flowing current with a voltage bias, we took care of avoiding pinch off condition, we changed measuring parameters in such a way that the voltage drop in the transport measurement was similar to or even smaller than the voltage drop in the thermopower measurement, we measured data points closely spaced in the V g axis. In all our attempts, oscillatory behavior was never observed in the resistance curves of any of the investigated samples. We rule out the slight intrinsic spatial inhomogeneity in these systems as a cause for the absence of oscillations in R sheet curves, as it should have a smearing effect on S and R sheet curves alike. On the other hand, similar discrepancy between R sheet (V g ) and S(V g ) curves is not noticed here for the first time. In other systems, e.g. GaAs/AlGaAs [1,2] and MoS 2 [3] , oscillations and divergence of S at negative V g were observed, without corresponding oscillations in R. The key difference between Seebeck and resistance lies in the respective nature of these two measurements, namely the former is an open-circuit measurement with no net flow of electric charge, where the effect of the applied thermal gradient is counterbalanced by a voltage drop, while the latter is a closed-circuit measurement with injected current. We suggest that these differences, combined with strong non ohmic behavior of LAO/STO in the depletion regime, widely observed in literature [4] , play a decisive role in washing out features associated to localized level filling in the resistance measurements. Similar non linear effects may be also responsible for the absence of oscillations in several resistivity measurements reported in literature [ 1 , 3 ] , where diverging R and S in the depletion regime are likely associated to non ohmic transport mechanisms. On the contrary, literature works reporting oscillations of S versus V g , associated to very low and nondiverging resistance versus V g curves, evidence oscillations both in Seebeck and resistance curves [5,6] .
Oppositely, the measurement of S is unobtrusive and permits an ultra-refined scan of the unperturbed, closely-spaced polaronic densities lying below the band edge. Hence, the oscillations of S should be thus ascribed to a paramount advantage of thermoelectric spectroscopy: the capability to probe extremely small changes in the charge density, without introducing significant perturbations. For this aspect, thermopower is far superior not only to resistivity measurements, but also to most of nowadays available spectroscopic techniques.  [7] to metal-insulator transition at a specific charge threshold, magnetic correlations [8] , and phase separation [9] . Furthermore, there is a number of spectroscopic measurements (see e.g. HAXPES spectra in [10] ) where a peak of Ti 3+ character is clearly detected, attributed to localized electrons occupying a small fraction of Ti sites.
All these evidences consistently suggest the presence of localization in some form in the highly depleted regime (typically below few 10 12 cm -2 ). On the other hand, what we define a 'direct' evidence has to do with a detailed scanning of the electronic structure of these states. This aspect has been, to our knowledge, lacking so far, and in fact, we claim that the Seebeck oscillations measured in this work are the first direct evidence of this sort. There are two major reasons at the basis of this unprecedented achievement: a) the specific characteristics of the thermoelectric measurement, sensitive to charge density fluctuations as small as 10 10 cm -2 which translates, in terms of energies, to a resolution of the order of meV and below. This is a resolution that few spectroscopies can achieve, in practice.
b) the induction of a strongly-depleted regime: as long as the metallic regime is retained, no technique could be capable to distinguish the presence of tiny densities of polaronic levels. The complete depletion of the conduction band is 'conditio-sine-qua-non' for the emergence of localized behavior in transport, thermoelectric, or even electric and optical measurements.
Consider for example ARPES measurements, which is the prototypical approach to explore the electronic structure of a system: experiments carried out so far on LAO/STO have employed soft-X ray ARPES which at its best is performed with a resolution of 50 meV [11,12,13] , i.e. a resolution way too low to distinguish levels separated a few meV one from each other and from the band bottom, as those described in our work. Furthermore, the well known "waterfall effects" creating blurred trails below the band bottom would make any attempt to distinguish features from these localized states totally hopeless. Finally, we point out that UV-ARPES with less than 10 meV resolution could be done, but the electron escape depth would be as short as few Angstrom, thus not enough to overcome the LAO barrier. More generally for what concern photoemission measurements, a long list of works can be found for LAO/STO with different techniques (HAXPES, RIXS, XPS, XAS, etc.) but after a careful exploration of this massive literature, we can safely conclude that none of those present the suited characteristics in terms of energy and charge density resolution to reveal the electronic structure at the same level of detail as that inferred from our Seebeck measurements.
A possible suitable technique to detect such localized states, could be scanning tunneling spectroscopy. With this technique, operating at low temperature (5K) it could be possible in principle to resolve electronic levels that are even more closely spaced than the ones observed in our experiment.
However the design of such an experiment should face several difficulties, mostly related to the buried interface and to the contributions to the DOS coming also from surface and defect states. To the best of our knowledge, only one paper in literature shows scanning tunnel spectroscopy on LAO/STO but, in that case, data were averaged on 75 meV [14] , so that features on the meV scale were not visible.

Supplementary Note 4: Phonon-drag model
Introduction -Starting from the general expression for phonon-drag [15,16,17,18,19,20] , we develop a specific formulation adapted to our means. We use a formalism based on 3D electronic structure coupled with 3D phonons, instead of the more common 2D electronic structure coupled with 3D phonons. Use of 3D formulas for the electrons is motivated by the need to treat on equal footing the bands and localized states, without a-priori assumption on the dimensionality. Of course, the electronic structure is fully anisotropic, and the 2D limit is easily recovered by imposing large effective masses in the interface-orthogonal direction. The multi-band effective mass modeling includes both delocalized conduction states and localized polaronic states (the latter having renormalized mass [21] and finite bandwidth) as shown in Figure 3 of the main text.
Since we are only interested in the low-T limit, only the interaction of electrons with acoustic phonons will be considered as a source of phonon-drag. The interaction is treated according to the deformation potential approach plus piezoelectric scattering to account for the polar character of the system. For simplicity, only one acoustic phonon branch and intra-band electron-phonon scattering is considered. These approximations may be seriously detrimental for detailed quantitative predictions, but are sufficient for the purpose of reproducing the order of magnitude and major features of the observed behavior.
From the coupled Boltzmann equations for electrons and phonons, the phonon-drag in the j direction is expressed as: where the factor 2 accounts for spin degeneracy, e is the electron charge,  j the electron conductivity, and V the volume; q and  q are phonon wavevector and frequency, respectively; the electron-phonon relaxation time, whereas includes all the other relevant phonon scattering processes (phonon-phonon scattering, boundary scattering, impurity scattering, etc.). j V is a velocity factor: and ) ( k n v j phonon and electron velocities, respectively. From Equation (1) it is apparent that phonon-drag essentially depends on two dominant factors: the inverse electric conductivity and the ratio of electron-phonon to total phonon scattering.
Clearly, phonon-drag only matters if EPR is non-discardable with respect to the other phonon scattering processes. Following Cantrell notations [16] , the inverse electron-phonon relaxation time (i.e. electron-phonon scattering frequency) is: where k n f and N q are Fermi-Dirac and Phonon equilibrium distributions, respectively, A(q) the electron-acoustic phonon coupling amplitude, and the two delta functions impose energy and momentum conservation for the absorption process (emission is accounted by the velocity factor).
Equation (1) implicitly assumes the relaxation time approximation for phonons: the decrease of phonon population Ph N q due to phonon scattering is: Velocity factor -In our model: The minus sign derives from the fact that for positive band curvature (thus electrons) the band velocity increases with k.
Energy-conserving delta function -To handle non isotropic effective masses, it is useful to introduce the following change of variables: where  is the angle formed by vectors K and Q. We then operate another change of variable from  to X defined as: The integral over X is only non vanishing for X 0 included in the integral limits , that is: the inequality can be rewritten: which has a simple interpretation: after absorption the carrier energy must be higher (lower) than the band energy corresponding to antiparallel (parallel) K and Q orientation. For conduction (valence) electrons, only the right (left) inequality matters.
Sum over crystal momentum -The sum over k is transformed first into a 3D-integral, then changed to an integral over K in polar coordinates. We integrate over the equatorial angle of K space, and take the azimuth angle as the  between K and Q, and then change  to X. Solving for the energy delta function we finally obtain: N is the number of bands. The minus sing comes from the velocity factor. Finally, K can be changed with the energy (Equation (10) where 0 nk  and W n are band bottom and bandwidth, respectively, and we have assumed, as customary, that the electronic relaxation time only depends on the wavevector through the energy.
Integration over phonon wavevector -The modulus of the rescaled variable Q in Equation (16) includes the band masses according to Equation (8); this would require to solve a 3D integration over the phonon wavevector. To simplify the calculation, we assume the following approximation: Here the phonon wavevector is integrated up to the Debye frequency  Debye temperature (500 K for STO [22] ), and: Electron-phonon coupling -In the simplest treatment, the electron-acoustic phonon scattering for an isotropic phonon distribution can be taken as the sum of two terms [23,24] : where D is the deformation potential,  the mass density, 0  and  vacuum permittivity and dielectric constant, em K the 3D-averaged electromechanical coupling, and q 0 the Debye screening length. The first term ("deformation potential scattering") describes the coupling of electrons with longwavelength longitudinal acoustic waves treated as an homogeneous strain. The second ("piezoelectric scattering") is the additional contribution due to the coupling with the electric field produced by the strain. For a non-polar system ( em K =0) or in the limit of large doping concentration (i.e. strong Debye screening) the second term vanishes and only the deformation potential contributes to the acoustic scattering. In case of small screening (q 0 =0) and highly ionic compounds, on the other hand, the piezoelectric contribution (1/q) may become dominant at small q. For some parameters appearing in Equation (22) the appropriate values to be used for our LAO/STO model are not easily determined, since these measurements may crucially depend on structural and chemical composition of samples, as well as doping type and concentration, thus they should be considered as merely indicative. We take em K = 0.37, i.e. the average of planar and transversal components measured for a solid solution of Tibased ceramics [25] . The dielectric constant under large negative gate is also difficult to be quantified reliably. It has been demonstrated [26] that below T=50 K a large negative gate field induces a polar transition in STO which sweeps away the quantum paraelectric state and largely reduces the screening capability. Tentatively, we use 300  k , that is the room-T value for STO bulk. Furthermore, we can assume a very long Debye length so that carrier screening is discardable ( ). This is indeed a basic hypothesis for having strong electron-phonon scattering [ 27 ] , i.e. in the low-density chargelocalized limit the carrier mobility is so small that the piezoelectric interaction becomes unscreened.
Finally for LAO/STO we use D=8.2 eV (used for GaAs in Ref. [16] ) and (that is the 3Daverage sound speed at low temperature [28] ).Equation (22) which is finally integrated numerically.

Supplementary note 5: Phonon relaxation time
For the electron-acoustic phonon scattering we start from: For what concerns other phonon scattering processes, a simple parametrization is proposed in the seminal article of Callaway [29] : where A describes the scattering by point impurities, B the phonon-phonon scattering including both normal (crystal momentum-conserving) and umpklapp scattering, and L v / s the boundary scattering (L is a characteristic sample length; in our model L=1 mm). In the following numerical examples we will use values for A and B given in Ref. [29] . While probably not appropriate for quantitative predictions, Equation (29) is sufficient for the purpose of qualitative description of general trends.

Supplementary Note 6: Diffusive Seebeck model
The formulation of diffusive Seebeck for a multiband effective-mass model is based on the Boltzmann Transport Equation in relaxation time approximation, and band energies modeled by nonisotropic effective masses. The model is described in detail in Ref. [30] where it was first applied to LAO/STO. Here we can just recall the final expression. The DC conductivity (in 3D) in direction j is given by: The Seebeck associated to the n th band is Finally, the diffusive Seebeck is obtained as:  (21)) is real and does not cancel at low temperature, thus no such limitation exists on the attainable amplitude.

Supplementary Note 7: Electronic relaxation time
Since we are mainly interested in the very low-T regime, we can consider a minimal model for the electronic relaxation time including acoustic-phonon scattering (AP) and impurity scattering (IS). AP is treated within elastic deformation potential approach: where the electron energy is relative to the CBB. For IS we adopt the well known Brooks-Herring formula: where n I is the impurity concentration, Z the impurity charge. For localized states, we assume the following hopping frequency expression:  is a characteristic hopping time, and the two exponentials represent the probability of hopping by thermal excitation (for energies lower than the mobility edge CBB  ) and by tunneling across an energy barrier E T , respectively. For the LAO/STO well simulation we use h0  = 10 -11 sec, and E T = 3 meV (that is in our model the energy separation between two localized states). the inverse Bose occupancy, and the integral over the electron-hole occupancy (see Equations (28), (29) and (30)). For what concern the 1 ep   dependence on the integral, in the T=0 limit the phonon can only be absorbed by electrons in the energy range ] , [

Supplementary Note 8: Test cases
. Thus, the number of available electronic states, and so the scattering rate, grows linearly with q   and saturates at ) ( In Supplementary Figure 2(B) we report the total phonon scattering frequency has a strong T 3 dependence, it becomes preponderant over electron-phonon at increasing T, causing the phonon-drag to quickly fade away with increasing T. Indeed, only at a very low temperature the profile of total scattering is visibly affected by electron-phonon. The phonon-drag T-behavior can be understood considering the relative scattering ratio (see Equation (1) and (3) Bulk SrTiO 3 n-doped (band plus localized states) -It is well known that the transport properties of STO are characterized by polaronic behavior and electron localization. Similar behavior has been also reported for LAO/STO [31] . A simple way to model these states is through poorly-dispersed massrenormalized bands laying below the CBB, and assuming that conduction across these states can only occur by hopping. In the example shown in Supplementary Figure 3  For what concern the fundamental understanding of the localized states, in literature they are typically related to two possible sources, which however do not exclude each other and can occur simultaneously: a) Anderson localization (disorder) which is known to affect the LAO/STO interface; (the presence of an Anderson tail of localized states just below the conduction edge of oxide heterostructures is well documented in literature [32] ); b) Mott localization, invoked to explain other correlation phenomena observed in this system (see e.g. capacitance enhancement described in [ 7 ] ).
Indeed, at very low charge density, the charge can localize in Ti 3+ sites, eventually helped by structural deformations (polarons) even in absence of structural disorder.
To mimic the experiment, the simulation assumes that at zero voltage the lowest conduction band is occupied by a mobile carrier density (n 2D 1.210 13 cm -2 ) typical for LAO/STO. Then, it is assumed that the leading effect of a negative V g is the progressive charge depletion, thus mapping S(E F ) is substantially equivalent to mapping S(V g ). Seebeck and resistivity are then calculated as a function of the progressively decreasing E F . Additional effects of negative V g , such as a possible change in the STO dielectric permittivity at the interface [26] , are inessential since our simplified model only includes by construction electronic states confined in a single interface layer. We underline that specific quantitative features of the model (energies, effective masses, bandwidths) are set to reproduce as closely as possible the characteristics of the S measurements in terms of number of oscillations, oscillation amplitude, S vs. V g absolute value. Of course these characteristics change, within a certain extent, from sample to sample, and so must do the parameters of the model. In particular the values in Supplementary Table I  The resulting phonon-drag, diffusive Seebeck, and resistivity are displayed in Figure 3 of the main text. Here we add (Supplementary Figure 4) the analysis of their state-by-state contributions, which is useful to highlight some key aspects of their behavior.
Consider first the diffusive Seebeck of each individual state (Supplementary Figure 4B)): each of them displays a negative-S region followed, after crossing the S=0 axis, by a positive bell-like feature.
Indeed, according to the Mott formula [33]   Notice also that for a model made of multiple overlapping bands shifted from each other in energy at the bottom, E F would not cross any DOS peaks, and in turn S d could display non-zero centered oscillations corresponding to the E F crossing the onset of each band bottom. However, as extensively discussed in the main article, in the band regime the amplitude of diffusive Seebeck at low T is totally discardable on the scale of the measured Seebeck. In our intensive attempts to simulate the experiment on the basis of a multi-band modeling, S d could never be found to exceed a few tents of V/K at low T, no matter how large the DOS slope (i.e. the effective masses) could be. On the basis of this analysis, we can confidently exclude that diffusive Seebeck could ever cause the huge, furiously oscillating thermopower observed in the measurements. behavior (Equation (22)). We see that starting from low doping (red solid curve), the scattering first increases with E F due to the progressive increase of W and, in turn, of available electronic transitions. Then, a regime change occurs (highlighted by the red-to-blue color change) and the scattering start to decrease with E F in consequence of the progressive decrease of the effective masses. For each curve, the phonon frequency region where the scattering is effective roughly follows the DOS profile of the electronic state crossed by E F at that doping.