Origin of the crossover from polarons to Fermi liquids in transition metal oxides

Transition metal oxides host a wealth of exotic phenomena ranging from charge, orbital and magnetic order to nontrivial topological phases and superconductivity. In order to translate these unique materials properties into device functionalities these materials must be doped; however, the nature of carriers and their conduction mechanism at the atomic scale remain unclear. Recent angle-resolved photoelectron spectroscopy investigations provided insight into these questions, revealing that the carriers of prototypical metal oxides undergo a transition from a polaronic liquid to a Fermi liquid regime with increasing doping. Here, by performing ab initio many-body calculations of angle-resolved photoemission spectra of titanium dioxide, we show that this transition originates from non-adiabatic polar electron–phonon coupling, and occurs when the frequency of plasma oscillations exceeds that of longitudinal-optical phonons. This finding suggests that a universal mechanism may underlie polaron formation in transition metal oxides, and provides a pathway for engineering emergent properties in quantum matter.

lucidating the nature of charge carriers in doped transition metal oxides (TMOs) is key to understanding the mechanism of electrical conduction in these multifunctional materials. In conducting oxides the infrared-active vibrations can couple strongly to electrons, leading to the formation of polarons 1 . Polarons are electrons dressed by a phonon cloud 2 , and represent a paradigmatic example of emergent state in condensed matter. Depending on their mass and size, polarons exhibit widely different conduction mechanisms, from band-like transport to thermally activated hopping transport 3,4 . Despite being central to the science and technology of oxides, little is known about the properties of polaronic states.
The interest in electron-phonon coupling and polaronic quasiparticles in TMOs has been reinvigorated by recent angleresolved photoelectron spectroscopy (ARPES) experiments [5][6][7][8][9] . The signature of polaronic behaviour in ARPES spectra is the appearance of satellites below the conduction band, at integer multiples of the optical phonon energy. This is reported in Fig. 1a,b for the paradigmatic case of doped anatase TiO 2 (ref. 5). These pioneering measurements showed that by increasing the carrier concentration, polaronic satellites gradually evolve into the photoemission kinks observed in metals and superconductors 10 (see Fig. 1c). It was proposed that this crossover reflects the evolution of charge carriers from polarons to a Fermi liquid 5,8 . In order to clarify the origin of this transition without making any a priori assumption about the underlying mechanism, first principles calculations are urgently called for. However, the investigation of polaronic features in ARPES spectra from first principles and their evolution with doping is exceptionally challenging and has never been reported before.
In the following we focus on the prototypical example of anatase TiO 2 . On top of its well-known applications in solar energy harvesting 11,12 and superhydrophilic technology 13,14 , this material is also being investigated in the quest for transparent conducting oxides based on non-toxic and Earth-abundant elements 15,16 . Despite its pivotal role in a broad range of technologies, the nature of the charge carriers in anatase is still controversial 17 . Here we address these issues by calculating ARPES spectra and polaron wavefunctions entirely from first principles. We develop a theoretical and computational framework that allows us to investigate polarons and Fermi liquid quasiparticles on the same footing, and without resorting to any empirical parameters. Using this approach, we show how the interplay between the dynamical screening of the electron plasma and the Fröhlich electron-phonon coupling is responsible for the transition between polaronic and Fermi liquid states. We propose that the mechanism identified in this work may be universal, and also applies to other oxides such as SrTiO 3 and ZnO.

Results
Angle-resolved photoemission spectra. Our calculated ARPES spectra are shown in Fig. 1d-f, for the same doping levels as in the measurements of ref. 5, reproduced in Fig. 1a-c. These maps show the bottom of the conduction band of n-doped anatase TiO 2 , for three doping levels in the range 10 18 -10 20 cm À 3 . All the spectra exhibit a bright parabolic band, whose size increases with doping. This reflects the rise of the Fermi energy inside the conduction band as the electron density increases. Besides this bright feature, panels a-b (experiments) and d-e (calculations)  Figure 1 | Ab initio ARPES spectra of n-doped anatase TiO 2 . (a-c) ARPES spectra of anatase TiO 2 measured by Moser et al. 5 . The measurements were taken at 20 K on samples with 5 Â 10 18 cm À 3 (a), 3 Â 10 19 cm À 3 (b) and 3.5 Â 10 20 cm À 3 (c). The zero of the energy is set to the Fermi level. The electron momentum k x is along the GS line of the anatase Brillouin zone (see j). Reproduced with permission from ref. show each a pair of satellites, a bright one at a binding energy around 0.1 eV, and a dim one near 0.25 eV. These features are identified as polaronic effects 5 . Moving on to higher doping in panels c and f, the satellites disappear and are replaced by band structure kinks near 0.1 eV. Overall, our calculated ARPES spectra are in remarkable agreement with the experiments of ref. 5. In order to achieve this unprecedented level of precision without any adjustable parameters, we developed an innovative computational framework.
In our calculations the photoelectron intensity maps are obtained using the single-particle spectral function A(k, o), where :k and :o are the crystal momentum and binding energy, respectively, and : is the reduced Planck constant. The spectral function is calculated using the state-of-the-art cumulant expansion, which has been employed recently to successfully describe plasmon satellites [18][19][20][21] . The cumulant expansion can naturally be applied to investigate the spectral properties of a polaronic system since the theory stems from the exact solution of the Fröhlich electron-boson coupling Hamiltonian 22,23 . In this formalism the spectral function is expressed as ref. 24: where n is the electron band index, e nk is the electron eigenvalue in absence of many-body interactions, and t is the time variable. The quantity C nk is the so-called cumulant function, and can be calculated by using the standard electron-phonon self-energy S nk as a seed (see Methods and Supplementary Note 1): Here g mnn (k, q) is the electron-phonon vertex, and describes the probability amplitude for an electron in the initial state nk j i to be scattered into the final state mk þ q j i by a phonon with momentum :q and energy :o qn in the branch n (ref. 25). The terms n qn and f mk þ q denote the Bose-Einstein and Fermi-Dirac occupations, respectively; O BZ is the Brillouin zone volume, and Z a positive infinitesimal. The numerical evaluation of equations (1) and (2) is very challenging owing to the singular nature of the Fröhlich interaction at long wavelengths 26 . To overcome this challenge we use the Wannier function technique of refs 27,28, as implemented in the EPW code 29 .
The strength of the electron-phonon interaction is most commonly expressed in terms of a single parameter l, which describes the enhancement of the electron mass from the band effective mass m b to the polaron mass m* via m* ¼ m b (1 þ l), as well as the renormalization of the Fermi velocity 30 . To determine this parameter we first extract the band structures underlying Fig. 1d-f, by tracking the maxima in the energy distribution curves. This procedure yields a set of parabolic bands and their satellites in Fig. 1g,h, and a distorted parabola in Fig. 1i (blue curves). For comparison we also report the electronic bands in absence of electron-phonon interactions (red curves). For carrier concentrations of 5 Â 10 18 and 3 Â 10 19 cm À 3 , polaron satellites are clearly visible in Fig. 1g,h, whereas at 3.5 Â 10 20 cm À 3 we see a band structure kink but no satellites. From these band structures we obtain the mass enhancement parameter as the ratio between the Fermi velocity of the bare band, v 0 F , and that of the dressed band, v F , as indicated in Fig. 1i: l k ¼v 0 kF =v kF À 1 (ref. 30), where we explicitly included the dependence on the wavevector at the Fermi surface. As we move from the lowest to the highest doping level we obtain l k ¼ 0.73, 0.70 and 0.20, respectively. Our calculated value at intermediate doping is in excellent agreement with the mass enhancement determined in experiments, l ¼ 0.7 (ref. 5).
In Fig. 1k we show the spectrum calculated at intermediate doping for electron momenta along GX and GZ, as well as the Fermi surface cut. Here we see that the Fermi surface pocket is elongated along the GZ direction. This elongation reflects the anisotropic character of the band effective masses, which we calculate to be m Origin of satellites and kinks. Having established that our calculations can accurately reproduce experimental spectra without adjustable parameters, we now proceed to identify the mechanisms that drive the formation of polaron satellites by selectively turning off individual components of the calculations. The energy separations between the quasiparticle bands and the first satellite in Fig. 1d,e are 106 and 124 meV, respectively. In the case of Fig. 1f the kink appears at a binding energy of 100 meV. These energy scales are compatible with a Fröhlich-type coupling to the longitudinal-optical (LO) E u phonon at 109 meV (see phonon dispersions in Supplementary Fig. 1). This vibrational mode corresponds to the stretching of the Ti-O bonds in the ab plane, as shown in Fig. 2e. Another candidate bosonic mode is the c-axis A 2u phonon at 88 meV (ref. 31); however, the energy of this mode appears too small to account for the satellites and kinks in Fig. 1. In order to quantify the importance of the E u phonon, in Fig. 2a we compare two calculations: the complete spectrum at intermediate doping and a calculation where the coupling to modes with energy above 100 meV is artificially suppressed. We see that, upon removing high-energy phonons, the intensity of the first satellite decreases and the second satellite disappears. The effective mass is also visibly lower; in fact, the analysis of the mass enhancement parameter yields l ¼ 0.3, to be compared to the total coupling l ¼ 0.7. It follows that the E u phonon contributes B60% of the total coupling, hence it represents the primary mechanism behind the satellites.
In Fig. 2b we test the importance of many-body correlations beyond the one-shot Migdal approximation. The spectrum at intermediate doping, obtained from the cumulant expansion, is shown together with the result of a calculation within the one-shot Migdal approximation (see Methods). The Migdal approximation is obtained from the complete electron-phonon self-energy by neglecting the three-point vertex 2 . In this approximation the electron-phonon self-energy contains only non-crossing diagrams. The additional approximation adopted in first principles calculations is to replace the fully renormalized electron and phonon propagators by those evaluated within density functional theory (DFT) 25 . We name this choice one-shot Migdal approximation to emphasize the lack of self-consistency in the electron propagator. From Fig. 2b we notice that the one-shot Migdal approximation fails twice: firstly, the separation between the quasiparticle band and the first satellite is too large as compared to experiment (151 meV instead of B100 meV); secondly, the dim satellite around 0.25 eV is completely missing. This test highlights the crucial role of high-order electronphonon correlations in the description of polarons in conducting oxides. It would be interesting to see whether the deficiencies of the one-shot Migdal approximation could be avoided by performing a fully self-consistent Migdal calculation 32 , and to establish the importance of crossing diagrams which are included in the cumulant expansion; however this test is currently out of reach.
We now move to discuss the origin of the crossover between satellites and kinks in the spectra as a function of doping. The Fröhlich theory of polarons 26 considers a single electron added to a polar insulator. This description does not take into account that, as more electrons are added to the system, the polar electronphonon coupling is weakened by the electronic screening of the charge carriers. Following ref. 2, we treat this issue by screening the electron-phonon vertex by the frequency-dependent Lindhard dielectric function with the calculated effective mass and dielectric permittivity of anatase TiO 2 (see also Methods). For completeness we show in Supplementary Fig. 2 the dielectric screening as a function of carrier density and how the Fröhlich electron-phonon vertex is influenced by doping. In order to illustrate the importance of carrier screening, in Fig. 2c we compare two scenarios: the spectrum calculated at high doping by accounting for electronic screening and the same system, but this time ignoring the screening of the electron-phonon coupling by doped carriers. The result is striking: in the absence of screening one obtains a sharp polaron satellite, in stark disagreement with experimental evidence. This comparison indicates that a correct description of the electronic screening is absolutely crucial to capture the photoemission kink at high doping. On the contrary, when we repeat the comparison of Fig. 2c for the cases of low and intermediate doping, we find a completely different picture: at these doping levels the electronic screening does not play any significant role.
These conflicting observations can be rationalized by inspecting the timescales of lattice vibrations and electronic screening. The E u phonon vibrates with a period T ph ¼ 38 fs. The characteristic response time of the electronic screening is set by the plasma frequency of doped carriers, T el ¼ 2p/o p . In the case of n-doped TiO 2 the electrons occupy a singly degenerate parabolic band minimum; therefore o p ¼ ne 2 where n is the electron density, E 0 the permittivity of vacuum and E 1 the high-frequency dielectric constant of TiO 2 (ref. 2). For the electron densities considered in Fig. 1a-c we calculate T el ¼ 122, 50 and 15 fs, respectively. From these values we deduce that, at the lowest doping, the carriers are too slow to screen the long-range electric field generated by the oscillation of the E u phonon (122 fs versus 38 fs). In this case the screening is ineffective and the Fröhlich interaction dominates the spectrum. On the contrary, at the highest doping the electrons oscillate faster than the LO phonon (15 fs versus 38 fs). In this case the electronic screening is almost complete and the Fröhlich coupling is largely suppressed. In this regime the strength of the kink depends critically on the carrier concentration; with increasing doping the coupling to the LO phonons is gradually suppressed, and the ARPES spectrum is dominated by the weaker coupling of carriers to non-polar phonons 8 . These considerations are summarized in Fig. 2d, where we compare the energy of the E u phonon with the plasma energy :o p of the carriers, and we monitor the evolution of the coupling strength with doping. We can identify two regions: a polaronic regime, corresponding to the situation o p oo ph , and a Fermi liquid regime, corresponding to o p 4o ph . In the polaronic region the electronic screening is ineffective, we see satellites in the spectra, and the electronphonon mass enhancement is not sensitive to the doping level. In the Fermi liquid region the Fröhlich coupling is strongly suppressed, polaron satellites are gradually replaced by photoemission kinks, and the coupling strength decreases. To further validate this trend, we performed additional calculations for a doping concentration of 1 Â 10 20 cm À 3 in the transition region. We found that one satellite is still present in the spectrum (Supplementary Fig. 3). A careful investigation of the mass renormalization parameter yields l ¼ 0.34, thus confirming that the electron-phonon coupling is weakened by the electronic screening. The present analysis reveals that the origin of the crossover from polarons to a Fermi liquid in the ARPES spectra of doped TiO 2 is to be found in a form of electron-phonon coupling, which we refer to as a non-adiabatic Fröhlich interaction.
Given the qualitative change in the band structures at the crossover from the polaronic to the Fermi liquid regime, it is natural to ask how this evolution is reflected in the wavefunctions. In order to explore this aspect we calculate the polaron wavefunctions by generalizing the perturbation theory of ref. 2 to ab initio calculations. Figure 2f shows the square moduli of the polaron wavefunctions at the bottom of the conduction band near the maximum of the polaron wavefunction and a few unit cells away: in blue the intermediate doping level, which is inside the polaronic region of Fig. 2d; in yellow the high doping level, in the Fermi liquid region. The corresponding envelope functions are also shown. Here we see that the wavefunction of an electron in the Fermi liquid region is essentially akin to a periodic Bloch function. On the contrary, the wavefunction of an electron in the polaronic region shows spatial localization. To quantify the size of the polaron in this latter case we define a polaron radius r p using the half width at half maximum. From Fig. 2f we obtain r p ¼ 5.7 nm. This result indicates that, despite the qualitative changes in the ARPES spectra, we are in the presence of large polarons throughout the entire doping range, and that polarons in TiO 2 are considerably more delocalized than previously thought 5 . Since thermally activated hopping transport corresponds to r p of the order of the lattice constant 33 , our finding supports the notion that electrical conduction in anatase TiO 2 takes place via standard band-like transport. To avoid ambiguities we emphasize that the present result refers to intrinsic mobile polarons, not to localized electronic defect states such as those associated with O vacancies and which are not mobile 34 .

Discussion
The non-adiabatic Fröhlich mechanism identified here is simple enough that it is likely to play a role in many other conducting oxides. In order to test this hypothesis we estimate the critical density for the crossover in doped SrTiO 3 (ref. 8) and doped ZnO (ref. 9). In SrTiO 3 the plasma energy for an electron concentration of 10 20 cm À 3 is :o p B64 meV (ref. 35); therefore, to match the 100 meV LO phonon of SrTiO 3 one would need a doping level of 2.6 Â 10 20 cm À 3 . In the experiments of ref. 8 the carriers are confined in a thin surface layer corresponding to approximately three unit cells; therefore, we estimate the critical carrier density for the two-dimensional electron liquid in the range 3 Â 10 13 cm À 2 . This value is remarkably close to the critical density determined experimentally, 4 Â 10 13 cm À 2 (ref. 8). More accurate calculations will need to take into account the two-dimensional screening at the surface of SrTiO 3 and its quantum paraelectric nature. A similar analysis can be carried out for ZnO (ref. 9). In this case the energy of the highest LO phonon is 72 meV, while the plasma energy for a surface electron concentration of 7.5 Â 10 13 cm À 2 is 320 meV (ref. 36). The transition would then appear around 3.8 Â 10 12 cm À 2 , slightly below the density reported in ref. 9. In agreement with our estimate, the spectra of ref. 9 exhibit a behaviour which is intermediate between kinks and satellites.
In sum, our findings indicate that the electron-phonon coupling in TiO 2 is more complex than previously thought, and that the non-adiabatic Fröhlich coupling could be the unifying mechanism behind the transition from polarons to Fermi liquids in conducting oxides; explicit ab initio calculations will be required to confirm this point. Looking further ahead, our work suggests that a fine control of the interplay between lattice vibrations and plasma oscillations may offer a pathway for investigating emergent states in quantum materials, and provide opportunities in the development of quantum technologies based on TMOs.

Methods
Ground-state calculations. Ab initio calculations were carried out for anatase TiO 2 (space group I4 1 /amd), using the experimental lattice parameter a ¼ 3.784 Å (ref. 37). We used DFT within the generalized gradient approximation of Perdew, Burke and Ernzerhof 38 . The core-valence interaction was described by means of norm-conserving pseudopotentials, with the semicore Ti-3s and Ti-3p states explicitly taken into account. Electron wavefunctions were expanded in a planewaves basis set with kinetic energy cutoff of 200 Ry, and the Brillouin zone was sampled using a 6 Â 6 Â 6 Monkhorst-Pack mesh. Lattice-dynamical properties were calculated using density functional perturbation theory (DFPT). All DFT and density functional perturbation theory calculations were performed using the Quantum ESPRESSO 39 package.
Electron-phonon coupling. Calculations of electron-phonon couplings were performed using the EPW code 29 , the cumulant expansion was performed separately (see Supplementary Note 1). The Fröhlich electron-phonon matrix element was calculated using the method of ref. 28. An accurate description of the Fröhlich vertex as described in ref. 28 is essential: ignoring this effect leads to a severe underestimation of the mass enhancement 40 . To evaluate equation (2) we computed electronic and vibrational states as well as the scattering matrix elements on a 4 Â 4 Â 4 Brillouin-zone grid. These quantities were interpolated with ab initio accuracy onto a fine grid with 2 Á 10 6 random q-points using EPW. The positive infinitesimal Z was set to 10 meV. Temperature effects were accounted for by including the Fermi-Dirac occupation in the spectral function, corresponding to the experimental temperature of 20 K.
Doping. Doping was included using the rigid-band approximation since the system is degenerate. In fact, by considering the Mott criterion for the metal-insulator transition 41 Þis the effective Bohr radius, we obtain a value n c ¼ 1.3 Â 10 18 cm À 3 for the critical density that is below the doping levels investigated. The screening of the electron-phonon interaction arising from the doped carriers was taken into account by computing the dielectric function E(q, o) in the random-phase approximation 2,42 , for a homogeneous electron gas with the calculated effective mass m b and dielectric permittivity E 1 of anatase TiO 2 (refs 2,42) (see Supplementary Note 3 for a rationale). In this expression the plasma frequency directly reflects the carrier density; therefore, the influence of doping on the dielectric screening in anatase TiO 2 (ref. 43) is included in the calculations. The resulting non-adiabatic matrix element is g NA mnn ðk; qÞ¼g mnn ðk; qÞ=E q; o qn þ i=t nk À Á , where t nk is the electron lifetime near the band edge. Here we use :/t nk ¼ 55 meV from ref. 5. We checked that this choice does not affect our conclusions ( Supplementary Fig. 4). We note that, owing to the strong anisotropy of the Fermi surface in n-doped anatase TiO 2 , it is possible that the nominal doping levels estimated in ref. 5 using the Fermi momentum may underestimate the actual carrier densities in the measured samples.
Angle-resolved photoemission spectra. The cumulant function used in the calculation of the ARPES spectra is obtained from the self-energy of equation (2) as follows 18 : C nk ðtÞ ¼ðp' Þ À 1 R 1 À 1 do y o À e nk =' ð Þ e ðio þ ZÞt Im S nk e nk =' À o ð Þ = o À iZ ð Þ 2 , where y is the Heaviside function. More details on the cumulant expansion are provided in Supplementary Note 1. The calculations within the one-shot Migdal approximation shown in Fig. 2b were performed by using directly equation (2) inside the spectral function 25 : Aðk; oÞ¼2p À 1 P n jIm ' o À e nk À S nk ð Þ À 1 j. The spectra were computed without including photon matrix elements effects 10 . This approximation is justified since the variation of the orbital character over the B0.1 eV binding energy range shown in Fig. 1 is negligible 10 . In order to extract the bands shown in Fig. 1g-i we proceeded as follows. For the quasiparticle bands, which are the brightest features in Fig. 1d-f, we calculate the poles of the electron Green's function as E nk ¼ e nk þ Re S nk (E nk =' ), where e nk and E nk are the bare and the renormalized dispersions, respectively. For the satellites in Fig. 1d,e the previous procedure is not applicable since these features are not poles of the Green's function. In this case we directly identify the maxima in the energy dispersion curves 10 , and obtain the pairs of parabolic bands which are visible in Fig. 1g,h at binding energies below 0.1 eV.
Polaron wavefunction. In order to calculate the wavefunction of the polaron we generalized the perturbation theory approach of ref. 2 to the case of multiple electron bands and phonon branches, as well as ab initio electronphonon matrix elements. We foundc nk ðr; fs kp gÞ¼Cc nk ðrÞÂ½1 þ P n R dq XBZ g NA nnn ðk; qÞ 2 e iqÁr = e nk À e nk þ q À ' o qn À Á 2 i j0 p i, where c nk andc nk are the electron wavefunction without and with electron-phonon interactions, respectively. s kp È É indicate the nuclear coordinates, 0 p is the phonon ground state and C is a normalization constant. The details of the derivation are given in Supplementary Note 4. The wavefunctions in Fig. 2f were obtained by setting nk to the conduction band edge; the envelope functions were calculated by retaining only the term inside the square brackets in the above expression. A random grid with 2 Â 10 6 points was used for the Brillouin-zone integration, and a small imaginary part Z ¼ 10 meV was added to the denominator.
Data availability. The data that support the findings of this study are available from the corresponding author on request, as is the modified version of the EPW package to include screening and the custom code for the cumulant expansion.