Spectroscopic Properties of Polarons in Strongly Correlated Systems by Exact Diagrammatic Monte Carlo Method

TL;DR

This paper introduces advanced Diagrammatic Monte Carlo methods for accurately studying the spectral and dynamical properties of polarons and related strongly correlated electron-phonon systems, providing unbiased results.

Contribution

It develops novel, approximation-free Monte Carlo techniques for calculating Green functions and spectral functions in complex electron-phonon coupled systems.

Findings

Exact ground state properties of various polarons obtained.

Spectral functions and optical conductivities computed without bias.

Results applicable to multiple strongly correlated models.

Abstract

We present recent advances in understanding of the ground and excited states of the electron-phonon coupled systems obtained by novel methods of Diagrammatic Monte Carlo and Stochastic Optimization, which enable the approximation-free calculation of Matsubara Green function in imaginary times and perform unbiased analytic continuation to real frequencies. We present exact numeric results on the ground state properties, Lehmann spectral function and optical conductivity of different strongly correlated systems: Frohlich polaron, Rashba-Pekar exciton-polaron, pseudo Jahn-Teller polaron, exciton, and interacting with phonons hole in the t-J model.

Figures

Equations

$${\hat{H}}_0^{\text{mbox}par}=\sum _{\mathbf{k}}{\epsilon }_a(\mathbf{k})a_{\mathbf{k}}^{†}{a}_{\mathbf{k}}+\sum _{\mathbf{k}}{\epsilon }_h(\mathbf{k})h_{\mathbf{k}}{h}_{\mathbf{k}}^{†}$$

$${\hat{H}}_{\text{mbox}a-h}=-N^{-1}\sum _{\mathbf{pk}{\mathbf{k}}^{\prime }}\mathcal{U}(\mathbf{p},\mathbf{k},{\mathbf{k}}^{\prime })a_{\mathbf{p}+\mathbf{k}}^{†}{h}_{\mathbf{p}-\mathbf{k}}^{†}{h}_{\mathbf{p}-{\mathbf{k}}^{\prime }}{a}_{\mathbf{p}+{\mathbf{k}}^{\prime }}.$$

$${\hat{H}}_{\text{mbox}par-bos}=i\sum _{\kappa =1}^Q\sum _{\mathbf{k},\mathbf{q}}(b_{\mathbf{q},\kappa }^{†}-b_{-\mathbf{q},\kappa })$$

$$\left[{\gamma }_{aa,\kappa }(\mathbf{k},\mathbf{q})a_{\mathbf{k}-\mathbf{q}}^{†}{a}_{\mathbf{k}}+{\gamma }_{hh,\kappa }(\mathbf{k},\mathbf{q})h_{\mathbf{k}-\mathbf{q}}^{†}{h}_{\mathbf{k}}+{\gamma }_{ah,\kappa }(\mathbf{k},\mathbf{q})h_{\mathbf{k}-\mathbf{q}}^{†}{a}_{\mathbf{k}}\right]+h.c.$$

$${\hat{H}}_{\text{mbox}bos}=\sum _{\kappa =1}^Q\sum _{\mathbf{q}}{\omega }_{\mathbf{q},\kappa }{b}_{\mathbf{q},\kappa }^{†}{b}_{\mathbf{q},\kappa }.$$

$${\hat{H}}_0^{\text{mbox}PJT}=\sum _{\mathbf{k}}\sum _{i=1}^T{ϵ}_i(\mathbf{k})a_{i,\mathbf{k}}^{†}{a}_{i,\mathbf{k}},$$

$${\hat{H}}_{\text{mbox}par-bos}=i\sum _{\kappa }^Q\sum _{\mathbf{k},\mathbf{q}}\sum _{i,j=1}^T{\gamma }_{ij,\kappa }(\mathbf{k},\mathbf{q})(b_{\mathbf{q},\kappa }^{†}-b_{-\mathbf{q},\kappa })a_{i,\mathbf{k}-\mathbf{q}}^{†}{a}_{j,\mathbf{k}}+h.c.$$

$$L_{\mathbf{k}}(\omega )=\sum _{\nu }\delta (\omega -E_{\nu }(\mathbf{k}))|⟨\nu |a_{\mathbf{k}}^{†}|\text{mbox}vac⟩{|}^2$$

$${\hat{H}}_{\text{mbox}0}=\sum _{\mathbf{k}}ϵ(\mathbf{k})a_{\mathbf{k}}^{†}{a}_{\mathbf{k}}+\sum _{\mathbf{q}}{\omega }_{\mathbf{q}}{b}_{\mathbf{q}}^{†}{b}_{\mathbf{q}},$$

$${\hat{H}}_{\text{mbox}int}=\sum _{\mathbf{k},\mathbf{q}}V(\mathbf{k},\mathbf{q})(b_{\mathbf{q}}^{†}-b_{-\mathbf{q}})a_{\mathbf{k}-\mathbf{q}}^{†}{a}_{\mathbf{k}}+h.c..$$

$$G_{\mathbf{k}}^{{\mathbf{pp}}^{\prime }}(\tau )=⟨\text{mbox}vac\mid a_{\mathbf{k}+{\mathbf{p}}^{\prime }}(\tau )h_{\mathbf{k}-{\mathbf{p}}^{\prime }}(\tau )h_{\mathbf{k}-\mathbf{p}}^{†}{a}_{\mathbf{k}+\mathbf{p}}^{†}\mid \text{mbox}vac⟩.$$

$$G_{\mathbf{k},ij}(\tau )=⟨\text{mbox}vac\mid a_{i,\mathbf{k}}(\tau )a_{j,\mathbf{k}}^{†}\mid \text{mbox}vac⟩,i,j=1,2.$$

$$G_{\mathbf{k}}(\tau )=⟨\text{mbox}vac\mid a_{\mathbf{k}}(\tau )a_{\mathbf{k}}^{†}\mid \text{mbox}vac⟩.$$

$$G_{\mathbf{k}}(\tau )={\int }_0^{\infty }d\omega L_{\mathbf{k}}(\omega )e^{-\omega \tau },$$

$$G_{\mathbf{k}}(\tau \gg \max \left[{\omega }_{\mathbf{q},\kappa }^{-1}\right])\to Z^{(\mathbf{k})}\exp [-E_{\text{mbox}g.s.}(\mathbf{k})\tau ],$$

$$\begin{array}[]{l}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;G_{{\bf k}}(n,\tau;\,{\bf q}_{1},\ldots,{\bf q}_{n})\,=\\ \langle\mbox{vac}|\,b_{{\bf q}_{n}}({\tau})\cdots b_{{\bf q}_{1}}({\tau})\,a_{\bf p}(\tau)a^{{\dagger}}_{\bf p}\,b^{{\dagger}}_{{\bf q}_{1}}\cdots b^{{\dagger}}_{{\bf q}_{n}}|\mbox{vac}\rangle\;,\;\;{\bf p}\,=\,{\bf k}-\sum_{j=1}^{n}{\bf q}_{j}\;.\end{array}$$

$${\Psi }_{\text{mbox}g.s.}(\mathbf{k})=\sum _{i=1}^T\sum _{n=0}^{\infty }\sum _{{\mathbf{q}}_1...{\mathbf{q}}_n}{\theta }_i(\mathbf{k};{\mathbf{q}}_1,...,{\mathbf{q}}_n)c_{i,\mathbf{k}-{\mathbf{q}}_1...-{\mathbf{q}}_n}^{†}{b}_{{\mathbf{q}}_1}^{†}...b_{{\mathbf{q}}_n}^{†}\mid \text{mbox}vac⟩$$

$$Z^{(\mathbf{k})}(n)\equiv \sum _{i=1}^T\sum _{{\mathbf{q}}_1...{\mathbf{q}}_n}\mid {\theta }_i(\mathbf{k};{\mathbf{q}}_1,...,{\mathbf{q}}_n){\mid }^2$$

$$⟨N⟩\equiv ⟨{\Psi }_{\text{mbox}g.s.}(\mathbf{k})\mid \sum _{\mathbf{q}}{b}_{\mathbf{q}}^{†}{b}_{\mathbf{q}}\mid {\Psi }_{\text{mbox}g.s.}(\mathbf{k})⟩=\sum _{n=1}^{\infty }n{Z}^{(\mathbf{k})}(n)$$

$${\Psi }_{\text{mbox}g.s.}(\mathbf{k})=\sum _{\text{mbox}\mathbf{p}}{\xi }_{\text{mbox}\mathbf{kp}}(g.s.)a_{\text{mbox}\mathbf{k}+\mathbf{p}}^{†}{h}_{\text{mbox}\mathbf{k}-\mathbf{p}}^{†}\mid \text{mbox}vac⟩.$$

$$G_{\mathbf{k}}^{{\mathbf{p}=\mathbf{p}}^{\prime }}(\tau \to \infty )=\mid {\xi }_{\text{mbox}\mathbf{kp}}(g.s.){\mid }^2{e}^{-E_{\text{mbox}g.s.}(\mathbf{k})\tau }.$$

$$L_{\mathbf{k}}(\omega )={\hat{\mathcal{F}}}_{\omega }^{-1}\left[G_{\mathbf{k}}(\tau )\right].$$

$$\mathcal{I}(\omega )={\hat{\mathcal{F}}}_{\omega }^{-1}\left[\sum _{{\mathbf{pp}}^{\prime }}{G}_{\mathbf{k}=0}^{{\mathbf{pp}}^{\prime }}(\tau )\right].$$

$${\sigma }_{\beta \delta }(\omega )=\pi {\hat{\mathcal{F}}}_{\omega }^{-1}\left[⟨J_{\beta }(\tau )J_{\delta }⟩\right]/\omega .$$

$$G_{\mathbf{k}}(\tau )={⟨\text{mbox}vac\left|T_{\tau }\left[a_{\mathbf{k}}(\tau )a_{\mathbf{k}}^{†}(0)\exp \left\{-{\int }_0^{\infty }{\hat{H}}_{\text{mbox}int}({\tau }^{\prime })d{\tau }^{\prime }\right\}\right]\right|\text{mbox}vac⟩}_{\text{mbox}con};\tau >0.$$

$$G_{\mathbf{k}}(\tau )=\sum _{m=0,2,4\dots }^{\infty }\sum _{{\xi }_m}\int d{x}_1^{\prime }\cdots d{x}_m^{\prime }{\mathcal{D}}_m^{({\xi }_m)}(\tau ;\{x_1^{\prime },\dots ,x_m^{\prime }\}).$$

$${\tilde{G}}_{\mathbf{k}}^{(0)}(\tau )=\frac{1}{\pi }{\int }_{-\infty }^{\infty }d\omega e^{-\omega \tau }\frac{Im{\Sigma }_{\text{mbox}ret}(\mathbf{k},\omega )}{{\left[\omega -ϵ(\mathbf{k})-Re{\Sigma }_{\text{mbox}ret}(\mathbf{k},\omega )\right]}^2+{\left[Im{\Sigma }_{\text{mbox}ret}(\mathbf{k},\omega \right]}^2}$$

$${\mathcal{D}}_2(\tau ;\{{\tau }_2^{\prime },{\tau }_1^{\prime },\mathbf{q}\})=$$

$${\left|V(\mathbf{k},\mathbf{q})\right|}^2{D}_{\mathbf{q}}^{(0)}({\tau }_2^{\prime }-{\tau }_1^{\prime })G_{\mathbf{k}}^{(0)}({\tau }_1^{\prime })G_{\mathbf{k}-\mathbf{q}}^{(0)}({\tau }_2^{\prime }-{\tau }_1^{\prime })G_{\mathbf{k}}^{(0)}(\tau -{\tau }_2^{\prime }).$$

$$M=\frac{{\mathcal{D}}_m^{({\xi }_m)}(\tau ;\{x_1^{\prime },\dots ,x_l^{(new)},\dots ,x_m^{\prime }\})/W(x_l^{(new)})}{{\mathcal{D}}_m^{({\xi }_m)}(\tau ;\{x_1^{\prime },\dots ,x_l^{(old)},\dots ,x_m^{\prime }\})/W(x_l^{(old)})}.$$

$$M=\frac{p_{\mathcal{A}}}{p_{\mathcal{B}}}\frac{{\mathcal{D}}_{m+2}^{({\xi }_{m+2})}(\tau ;\{x_1^{\prime },\dots ,x_m^{\prime };{\mathbf{q}}^{\prime },{\tau }^{\prime },{\tau }^{\prime \prime }\})}{{\mathcal{D}}_m^{({\xi }_m)}(\tau ;\{x_1^{\prime },\dots ,x_m^{\prime }\})W({\mathbf{q}}^{\prime },{\tau }^{\prime },{\tau }^{\prime \prime })}.$$

Full text

11institutetext: CREST, Japan Science and Technology Agency (JST), AIST, 1-1-1, Higashi, Tsukuba 305-8562, Japan. 22institutetext: Russian Research Centre “Kurchatov Institute”, 123182 Moscow, Russia. 33institutetext: Department of Applied Physics, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113, Japan.

Spectroscopic Properties of Polarons in Strongly Correlated

Systems by Exact Diagrammatic Monte Carlo Method

A. S. Mishchenko 1122   

N. Nagaosa 33

1 Introduction

Theoretical study of polarons in the strongly correlated system is like an attempt to view contents of a Pandora box embedded into another, even more sinister and obscure, container of riddles, enigmas and mysteries. This desperate situation occurs because solution is not known even for the simplest polaron problem, i.e. when a perfectly stable quasiparticle (QP) with momentum as a single quantum number interacts with a well defined bath of bosonic elementary excitations. To the contrary, the definition of the strongly correlated system implies that QPs might be highly unstable and the very notion of QPs, both in electronic and bosonic subsystems, is under question. Thus, one faces the problem of an interplay between ill defined objects and it is crucial to solve the problem without approximations. Further difficulty, pertinent to realistic systems, is an interplay of the momentum and other quantum numbers characterizing internal states of a QP.

The problem of polaron originally emerged as that of an electron coupled to phonons (see Appel ; Pekar51 ). In the initial formulation a structureless QP is characterized by the only quantum number, momentum, which changes due to interaction of the QP with phonons Landau33 ; Frohlich . Later, depending on what can be called “particle” and “environment”, and how they interact with each other, the polaron concept was related to extreme diversity of physical phenomena. There are many other objects which, having nothing to do with phonons, are isomorphic to simple polaron DevEnciclop , as, e.g. an exciton-polaron in the intraband scattering approximation Ansel ; UKKTH86 ; Toyo ; Toyobook . Another example is the problem of a hole in the antiferromagnet which is closely related to polaron since hole movement is accompanied by the spin flips which, in the spin wave approximation, are equivalent to creation and annihilation of magnons Kane ; Izymov97 .

The concept of polaron was further generalized to include internal degrees of freedom which, interacting with environment, change their quantum numbers. Example of a complex QP is the Jahn-Teller polaron, where electron-phonon interaction (EPI) changes quantum numbers of degenerate electronic states Kanamori ; Abragam ; Kugel . This generalization is important due to it’s relevance to the colossal magnetoresistance phenomena in the manganese oxides Millis95 ; AlexBrat99 . Another example is the pseudo Jahn-Teller polaron, where EPI is inelastic and leads to transitions between close in energy electronic levels of a QP RashPJT ; ToyPJT ; BersuBook . Further generalization is a system of several QPs which interact both with each other and environment. For example, effective interaction of two electrons through exchange by phonons can overcome the Coulomb repulsion and form a bound state, bipolaron Vinetskii ; AnderBP ; HiramToy ; AlexRunn ; AlexRun2 . On the other hand, coupling of attracting hole and electron to the lattice vibrations Haken56 ; Bassani75 ; Pollman77 can create a lot of qualitatively different objects: localized exciton, weakly bound pair of localized hole and localized electron, etc. Sumi77 ; UKKTH86 . Scattering by impurities introduces additional complexity to the polaron problem because interference of impurity potential with lattice distortion, which accompanies the polaron movement, can contribute either constructively or destructively to the localization of a QP on impurity Primesi1 ; Primesi2 ; UKKTH86 .

In addition, a bare QP and bosonic bath can not be considered as well defined in the correlated systems. Angle Resolved Photoemission Spectra (ARPES), revealing the Lehmann Function (LF) of quasiparticle, demonstrate broad peaks in many correlated systems: cooper oxide high-temperature superconductors ZX95 ; Shen_03 ; DichotomyZhou , colossal magnetoresistive manganites Dessau_98 ; Manella_04 ; Manella_05 , quasi-one-dimensional Peierls conductors Perf_01 ; Perf_02 , and Verwey magnetites Schrupp_05 . Besides, phonons are also broadened in many correlated systems, e.g. in high-temperature semiconductors Pintsch and mixed-valent materials Mook79 ; Mook82 . One of possible reasons for these broadenings is the interaction of the QPs with the lattice degrees of freedom. However, in many realistic cases other subsystems, not explicitly included into the polaron Hamiltonian, are responsible for the decay of QP and phonons, e.g., another electronic bands, phonon anharmonicity, interaction with nuclear spins, etc. Then, if this auxiliary broadening is known in some approximation, one can formulate an ambitious goal to study spectral response when “bare” quasiparticle with known damping interacts with “broadened” bosonic excitations.

No one of traditional numerical methods, to say nothing of analytical ones, can give approximation free results for measurable quantities of polaron, such as optical conductivity or angle resolved photoemission spectra, for in macroscopic system of arbitrary dimension. Besides, we are not aware of any numerical method which can incorporate in an approximation free way the information on the damping of QP and bosonic bath. Below we describe basics of recently developed Diagrammatic Monte Carlo (DMC) method for numerically exact computation of Green functions and correlation functions in imaginary time for few polarons in a macroscopic system Zhetf ; PS ; MPSS ; Exciton ; Aus02 ; Optics ; tJpho ; UFN_05 ; JPSJ_exc . Analytic continuation of imaginary time functions to real frequencies is performed by a novel approximation free approach of stochastic optimization (SO) MPSS ; UFN_05 ; JPSJ_exc , circumventing difficulties of popular Maximal Enthropy method. Finally we focus on results of application of the DMC-SO machinery to various problems MN01 ; Hole01 ; RP02 ; Isotope ; tJphoSpain ; FrCondon

The basic models, related to the polaronic objects in correlated systems, which can be solved by DMC-SO methods, are stated in the next Sect. It is followed in Sect. 1.2 by description of stumbling blocks encountered by analytic methods. Sect. 2 concerns the basics of DMC-SO methods. However, those who are not interested in the details of the methods can briefly look through the definitions in the introduction for Sect. 2 and turn to Sect. 3 where LF and optical conductivity of Fröhlich polaron are discussed (see also D ). Results of studies of the self-trapping phenomenon are presented in Sect. 4 and application of DMC-SO methods to the exciton problem can be found in Sect. 5. The chapter is completed by Sect. 6 devoted to studies of ARPES of high temperature superconductors.

1.1 Formulation of a General Model with Interacting Polarons

In general terms, the simplest problem of a complex polaronic object, where center-of-mass motion does not separate from the rest of degrees of freedom, is introduced as system of two QPs

[TABLE]

($a_{{{\bf k}}}$ and $h_{{{\bf k}}}$ are annihilation operators, and $\varepsilon_{a}({\bf k})$ and $\varepsilon_{h}({\bf k})$ are dispersions of QPs), which interact with each other

[TABLE]

($N$ is the number of lattice sites) through the instantaneous Coulomb potential and the scattering by bosons

[TABLE]

($\gamma_{[aa,ah,hh],\kappa}$ are interaction constants) where quanta of Q different branches of bosonic excitations are created or annihilated, which are described by

[TABLE]

In general, each QP can be a composite one with internal degree of freedom represented by $T$ different states

[TABLE]

which quantum numbers can be also changed due to nondiagonal part of particle-boson interaction

[TABLE]

Complicated model (1)-(6) is still too far from the cases encountered in strongly correlated systems. Due to coupling of QPs (1) and (5) and bosonic fields (4) to additional degrees of freedom, these excitations are not well defined from the onset. Namely, the dispersion relation of the QP spectrum $\epsilon({\bf k})$ in realistic system is ill-defined. One can speak of a Lehmann Function (LF) AGD ; Mahan ; JG of a QP

[TABLE]

,which is normalized to unity $\int_{0}^{+\infty}d\omega L_{\bf k}(\omega)=1$ and can be interpreted as a probability that a QP has momentum ${\bf k}$ and energy $\omega$. (Here $\{|\nu\rangle\}$ is a complete set of eigenstates of Hamiltonian $\hat{H}$ in a sector of given momentum ${\bf k}$: $H\,|\nu({\bf k})\rangle=E_{\nu}({\bf k})\,|\nu({\bf k})\rangle$.) Only for noninteracting system the LF reduces to delta function $L_{\bf k}^{\mbox{\scriptsize NONINT}}(\omega)\,=\,\delta(\omega-\epsilon({\bf k}))$ and, thus, sets up dispersion relation $\omega=\epsilon({\bf k})$.

Specific cases of model (1)-(6) describe enormous variety of physical problems. Hamiltonians (1) and (2), in case of attractive potential ${\cal U}({\bf p},{\bf k},{\bf k}^{\prime})>0$, describe an exciton with static screening Knox ; Egri . Besides, expressions (1)-(4) describe bipolaron for repulsive interaction Vinetskii ; AnderBP ; HiramToy ; AlexRunn ; AlexRun2 ${\cal U}({\bf p},{\bf k},{\bf k}^{\prime})<0$ and exciton-polaron otherwise Haken56 ; Bassani75 ; Pollman77 . The simplest model for exciton-phonon interaction, when only two ($T=2$) lowest states of relative electron-hole motion are relevant (e.g. in one-dimensional charge-transfer exciton Haarer75 ; HaPhi75 ; ElsWeis85 ), is defined by Hamiltonians (4)-(6)). The same relations (4)-(6) describe the problems of Jahn-Teller [all $\epsilon_{i}$ in Hamiltonian (5) are the same] and pseudo Jahn-Teller polaron. The problem of a hole in an antiferromagnet in spin-wave approximation is expressed in terms of Hamiltonians (4)-(6) with $Q=1$ and $T=1$. When hole also interacts with phonons, one has to take into account one more bosonic branch and set $Q=2$ in (4) and (6). Finally, the simplest nontrivial problem of a polaron, i.e. of a structureless QP interacting with one phonon branch, is described by noninteracting Hamiltonians of QP $\hat{H}_{\mbox{\scriptsize par}}$ and phonons $\hat{H}_{\mbox{\scriptsize ph}}$

[TABLE]

and interaction term

[TABLE]

The simplest polaron problem, in turn, can be subdivided into continuous and lattice polaron models.

1.2 Limitations of Analytic Methods in Problem of Polarons

Analytic solution for the problem of exciton in a rigid lattice is available only for small radius Frenkel regime Frenkel and large radius Wannier regime Wannier . However, even limits of validity for these approximations are not known. Random phase approximation approaches Knox ; Egri , are capable of obtaining some qualitative conclusions for intermediate radius regime though its’ quantitative results are not reliable due to uncontrolled errors. The situation is similar with the problem of structureless polaron, where analytic solutions are known only in the weak and strong coupling regimes. Besides, reliable results for these regimes are available only for ground state properties.

Although several novel methods, capable of obtaining properties of excited states, were developed recently, variational coherent-states expansion DeFilCat05 and free propagator momentum average summation Berciu as a few examples, all of them, to provide reliable data in a specific regime, need either comparison with exact sum rules DevreLeRo76 ; KornilEL02 or with exact numerical results.

Application of variational methods to study of excitations is a tricky issue since, strictly speaking, they are valid only for the ground state. As an example for the importance of sum rules in variational treatment, we refer to the problem of the optical conductivity of the Fröchlich polaron. Possibility of existence of Relaxed Excited State (RES), which is a metastable state where lattice deformation has adjusted to the electronic excitation rendering stability and narrow linewidth of the spectroscopic response, was briefly mentioned by S. I. Pekar in early 50’s. Then, conception of RES was rigorously formulated by J. T. Devreese with coworkers and has been a subject of extensive investigations for years DevEnciclop ; Dev66 ; KED69 ; Dev72 ; DeRev72 ; Dev_twopho_73 ; Optics ; FrCondon . Calculations of impedance Dev72 in the framework of technique Fey62 supported the existence of a narrow stable peak in the optical conductivity. However, even the authors of Dev72 were skeptical about the fact that the width of RES in the strong coupling regime appeared to be more narrow than the phonon frequency, i.e. inverse time which is, according to the Heisenberg uncertainty principle, is required for the lattice readaptation. In consequent paper Dev_twopho_73 they realized the importance of many-phonon processes and studied two-phonon contribution to optical conductivity. Importance of many phonon processes was confirmed when variational results Dev72 were compared with exact DMC simulations Optics . Variational result well reproduced the position of the peak in exact data though failed in description of the peak width in the strong coupling regime Optics . Finally, when approach Dev72 was modified and several sum rules were accurately introduced into variational model FrCondon , both position and width of the peak were quantitatively reproduced. Studies FrCondon (see Sect. 3.1), do not address rather philosophical question whether RES exists or not, though inevitably prove that, in contrast to the foregoing beliefs, there in no stable excited state of the Fröhlich polaron in the strong coupling regime. Note that sometimes excited states can not be handled by analytic methods even for weak couplings: perturbation theory expression for LF of the Fröhlich polaron model diverges at the phonon energy $\omega_{\mbox{\scriptsize ph}}$ [See (34) in Sect. 3.1.] and more elaborate treatment is necessary.

Difficulties of semianalytic methods enhance in the intermediate coupling regime where results are sometimes wrong even for ground state properties. For example, the variatioanl approach LLP53 , which has been considered as an intermediate coupling theory, appeared to be valid only in the weak coupling limit MPSS . Special interest to the methods, giving reliable information on excited states, is triggered by the self-trapping phenomenon which occurs just in the intermediate coupling regime. This phenomenon is a dramatic transformation of QP properties when system parameters are slightly changed Landau33 ; UKKTH86 ; Toyobook ; R82 . In the intermediate coupling regime “trapped” QP state with strong lattice deformation around it and “free” state with weakly perturbed lattice may hybridize and resonate because of close energies at some critical value of electron-lattice interaction $\gamma_{c}$. It is clear that, to study self-trapping, one has to apply a method giving reliable information on excited states in the intermediate coupling regime.

2 Diagrammatic Monte Carlo and Stochastic Optimization Methods

In this section we introduce definitions of exciton-polaron properties which can be evaluated by DMC and SO methods. An idea of DMC approach for numerically exact calculation of Green functions (GFs) in imaginary times is presented in Sect. 2.1, and a short description of SO method, which is capable of making unbiased analytic continuation from imaginary times to real frequencies, is given in Sect. 2.2. Using combination of DMC and SO, one can often circumvent difficulties of analytic and traditional numerical methods. Therefore, a brief comparative analysis of advantages and drawbacks of DMC-SO machinery is given in Sect. 2.3.

To obtain information on QPs it is necessary to calculate Matsubara GF in imaginary time representation and make analytic continuation to the real frequencies Mahan . For the two-particle problem (1)-(4) the relevant quantity is the two-particle GF Exciton ; Aus02

[TABLE]

(Here $h_{{\bf k-p}}(\tau)=e^{\hat{H}\tau}h_{{\bf k-p}}e^{-\hat{H}\tau}$, $\tau>0$.) In the case of exciton-polaron, vacuum state $\mid\mbox{vac}\rangle$ is the state with filled valence and empty conduction bands. For the bipolaron problem it is a system without particles. In the simpler case of a QP with two-level internal structure described by (4)-(6) the relevant quantity is the one-particle matrix GF MN01 ; Aus02

[TABLE]

For a structureless polaron the matrix (11) reduces to one-particle scalar GF

[TABLE]

Information on the response to an external weak perturbation (e.g. optical absorption) is contained in the current-current correlation function $\langle J_{\beta}(\tau)J_{\delta}\rangle$ ($\beta$/$\delta$ are Cartesian indexes).

Lehmann spectral representation of $G_{{\bf k}}(\tau)$ Mahan ; JG at zero temperature

[TABLE]

with the Lehmann function (LF) $L_{\bf k}(\omega)$ given in (7), reveals information on the ground and excited states. Here $\{|\nu\rangle\}$ is a complete set of eigenstates of Hamiltonian $\hat{H}$ in a sector of given momentum ${\bf k}$: $H\,|\nu({\bf k})\rangle=E_{\nu}({\bf k})\,|\nu({\bf k})\rangle$. The LF $L_{\bf k}(\omega)$ has poles (sharp peaks) on the energies of stable (metastable) states of particle. For example, if there is a stable state at energy $E({\bf k})$, the LF reads $L_{\bf k}(\omega)=Z^{({\bf k})}\,\delta(\omega-E({\bf k}))\,+\,\ldots$, and the state with the lowest energy $E_{\mbox{\scriptsize g.s.}}({\bf k})$ in a sector of a given momentum ${\bf k}$ is highlighted by asymptotic behavior of GF

[TABLE]

where $Z^{({\bf k})}$-factor is the weight of the state. Analyzing the asymptotic behavior of similar $n$-phonon GFs MPSS ; MN01

[TABLE]

one obtains detailed information about lowest state. For example, important characteristics of the lowest state wave function

[TABLE]

are partial $n$-phonon contribution

[TABLE]

which is normalized to unity $\sum_{n=0}^{\infty}Z^{({\bf k})}(n)\equiv 1$, and the average number of phonons

[TABLE]

in polaronic cloud. Another example is the wave function of relative electron-hole motion of exciton in the lowest state in the sector of given momentum

[TABLE]

The amplitudes $\xi_{\mbox{\scriptsize{\bf k p}}}(g.s.)$ of this wave function can be obtained Exciton from asymptotic behavior of the following GF (10)

[TABLE]

Information on the excited states is obtained by the analytic continuation of imaginary time GF to real frequencies which requires to solve the Fredholm equation $G_{{\bf k}}(\tau)=\hat{\cal F}\left[L_{\bf k}(\omega)\right]$ (13)

[TABLE]

The equation (13) is a rather general relation between imaginary time GF/correlator and spectral properties of the system. For example, the absorption coefficient of light by excitons ${\cal I}(\omega)$ is obtained as solution of the same equation Exciton

[TABLE]

Besides, the real part of the optical conductivity $\sigma_{\beta\delta}(\omega)$ is expressed Optics in terms of current-current correlation function $\langle J_{\beta}(\tau)J_{\delta}\rangle$ by relation

[TABLE]

2.1 Diagrammatic Monte Carlo Method

DMC Method is an algorithm which calculates GF (10)-(12) without any systematic errors. This algorithm is described below for the simplest case of structureless polaron MPSS , and generalizations to more complex cases can be found in consequent references111Generalization of described below technique to the case of exciton (1-2) is given in Exciton and its modification for pseudo-Jahn-Teller polaron (4-6) is developed in MN01 ; Aus02 . Method for evaluation of current-current correlation function can be found in Optics and a case of a polaron interacting with two kinds of bosonic fields is considered in tJpho .. DMC is based on the Feynman expansion of the Matsubara GF in imaginary time in the interaction representation

[TABLE]

Here $T_{\tau}$ is the imaginary time ordering operator, $|\mbox{vac}\rangle$ is a vacuum state without particle and phonons, $\hat{H}_{\mbox{\scriptsize int}}$ is the interaction Hamiltonian in (9). Symbol of exponent denotes Taylor expansion which results in multiple integration over internal variables $\{\tau_{1}^{\prime},\tau_{2}^{\prime},\ldots\}$. Operators are in the interaction representation $\hat{A}(\tau)=\exp[\tau(\hat{H}_{\mbox{\scriptsize par}}+\hat{H}_{\mbox{\scriptsize ph}})]\hat{A}\exp[-\tau(\hat{H}_{\mbox{\scriptsize par}}+\hat{H}_{\mbox{\scriptsize ph}})]$. Index “con” means that expansion contains only connected terms where no one integral over internal time variables $\{\tau_{1}^{\prime},\tau_{2}^{\prime},\ldots\}$ can be factorized.

Vick theorem expresses matrix element of time-ordered operators as a sum of terms, each is a factor of matrix elements of pairs of operators, and expansion (24) becomes an infinite series of integrals with an ever increasing number of integration variables

[TABLE]

Here index $\xi_{m}$ stands for different Feynman diagrams (FDs) of the same order $m$. Term with $m=0$ is the GF of the noninteracting QP $G_{\bf k}^{(0)}(\tau)$. Function ${\cal D}_{m}^{(\xi_{m})}(\tau;\{x_{1}^{\prime},\ldots,x_{m}^{\prime}\})$ of any order $m$ can be expressed as a factor of GFs of noninteracting quasiparticle, GFs of phonons, and interaction vortexes $V({{\bf k},\bf q})$. For the simplest case of Hamiltonian system expressions for GFs of QP $G_{{\bf k}}^{(0)}(\tau_{2}-\tau_{1})=\exp\left[-\epsilon({\bf k})(\tau_{2}-\tau_{1})\right]$ ($\tau_{2}>\tau_{1}$) and phonons $D_{\bf q}^{(0)}(\tau_{2}-\tau_{1})=\exp\left[-\omega_{\bf q}(\tau_{2}-\tau_{1})\right]$ ($\tau_{2}>\tau_{1}$) are well known.

An important feature of the DMC method, which is distinct from the row of other exact numerical approaches, is the explicit possibility to include renormalized GFs into exact expansion without any change of the algorithm. For example, if a damping of QP, caused by some interactions not included in the Hamiltonian, is known, i.e. retarded self-energy of QP $\Sigma_{\mbox{\scriptsize ret}}({\bf k},\omega)$ is available, renormalized GF

[TABLE]

can be introduced instead of bare GF $G_{{\bf k}}^{(0)}(\tau)$. Explicit rules for evaluation of ${\cal D}_{m}^{(\xi_{m})}$ do not depend on the order and topology of FD. GFs of noninteracting QPs $G_{{\bf k}}^{(0)}(\tau_{2}-\tau_{1})$ (or $\widetilde{G}_{{\bf k}}^{(0)}(\tau_{2}-\tau_{1})$) with corresponding times and momenta are ascribed to horizontal lines and noninteracting GFs of phonon $D_{\bf q}^{(0)}(\tau_{2}-\tau_{1})$ (multiplied by the factor of corresponding vortexes $V({\bf k}^{\prime},{\bf q})V^{*}({\bf k}^{\prime\prime},{\bf q})$) are attributed to phonon propagator arch (see Fig. 1a). Then, ${\cal D}_{m}^{(\xi_{m})}$ is the factor of all GSs. For example, expression for the weight of the second order term (Fig. 1b) is the following

[TABLE]

The DMC process is a numerical procedure which, basing on the Metropolis principle Metro ; LanBin , samples different FDs in the parameter space $(\tau,m,\xi_{m},\{x_{m}^{\prime}\})$ and collects statistics of external variable $\tau$ in a way that the result of this statistics converges to exact GF $G_{\bf k}(\tau)$. Although sampling of the internal parameters of one term in (25) and switch between different orders is performed within the the framework of one and the same numerical process, it is instructive to start with the procedure of evaluation of a specific term ${\cal D}_{m}^{(\xi_{m})}(\tau;\{x_{1}^{\prime},\ldots,x_{m}^{\prime}\})$.

Starting from a set $\{\tau;\{x_{1}^{\prime},\ldots,x_{m}^{\prime}\}\}$, an update $x_{l}^{(old)}\rightarrow x_{l}^{(new)}$ of an arbitrary chosen parameter is suggested. This update is accepted or rejected according to Metropolis principle. After many steps, altering all variables, statistics of external variable converges to exact dependence of the term on $\tau$. Suggestion for new value of parameter $x_{l}^{(new)}=\hat{S}^{-1}(R)$ is generated by random number $R\in[0,1]$ with a normalized distribution function $W(x_{l})$ in a range $x_{l}^{(min)}<x_{l}<x_{l}^{(max)}$. There are only two restrictions for this otherwise arbitrary function. First, new parameters $x_{l}^{(new)}$ must not violate FD topology, i.e., for example, internal time $\tau_{1}^{\prime}$ in Fig. 1c must be in the range $[x^{(min)}=0,x^{(max)}=\tau_{3}^{\prime}]$. Second, the distribution must be nonzero for the whole, allowed by FD topology, domain. This ergodicity property is crucial since it is necessary to sample the whole domain for convergence to exact answer. At each step, update $x_{l}^{(old)}\rightarrow x_{l}^{(new)}$ is accepted with probability $P_{acc}=M$ (if $M<1$) and always otherwise. The ratio $M$ has the following form

[TABLE]

For uniform distribution $W=\mbox{const}=\left(x_{l}^{(max)}-x_{l}^{(max)}\right)^{-1}$, the probability of any combination of parameters is proportional to the weight function ${\cal D}$. However, for better convergence the distribution $W(x_{l}^{new})$ must be as close as possible to the actual distribution given by function ${\cal D}_{m}^{(\xi_{m})}(\{\ldots,x_{l}^{(new)},\ldots,\})$.

For sampling over FDs of all orders and topologies it is enough to introduce two complimentary updates. Update ${\cal A}$ transforms FD ${\cal D}_{m}^{(\xi_{m})}(\tau;\{x_{1}^{\prime},\ldots,x_{m}^{\prime}\})$ into higher order FD ${\cal D}_{m+2}^{(\xi_{m+2})}(\tau;\{x_{1}^{\prime},\ldots,x_{m}^{\prime};\;{\bf q}^{\prime},\tau^{\prime}_{3},\tau^{\prime}_{4}\})$ with extra phonon arch, connecting some time points $\tau^{\prime}_{3}$ and $\tau^{\prime}_{4}$ by phonon propagator with momentum ${\bf q}^{\prime}$ (Fig. 1c). Note that the ratio of weights ${\cal D}_{m+2}^{(\xi_{m+2})}/{\cal D}_{m}^{(\xi_{m})}$ is not dimensionless. Dimensionless Metropolis ratio

[TABLE]

contains normalized probability function $W({\bf q}^{\prime},\tau^{\prime},\tau^{\prime\prime})$, which is used for generating of new parameters222The factor $p_{\cal A}/p_{\cal B}$ depends on the probability to address add/remove processes.. Complementary update ${\cal B}$, removing the phonon propagator, uses ratio $M^{-1}$ MPSS .

Note that all updates are local, i.e. do not depend on the structure of the whole FD. Neither rules nor CPU time, needed for update, depends on the FD order. DMC method does not imply any explicit truncation of FDs order due to finite size of computer memory. Ever for strong coupling, where typical number of phonon propagators $N_{ph}$, contributing to result, is large, influence of finite size of memory is not essential. Really, according to Central Limit Theorem, number of phonon propagators obeys Gauss distribution centered at $\bar{N}_{ph}$ with half width of the order of $\sqrt{\bar{N}_{ph}}$ Sandvik . Hence, if a memory for at least $2\bar{N}_{ph}$ propagators is reserved, diagram order hardly surpasses this limit.

2.2 Stochastic Optimization Method

The problem of inverting of integral equation (13) is an ill posed problem. Due to incomplete noisy information about GF $G_{\bf k}(\tau)$, which is known with statistic errors on a finite number of imaginary times in a finite range $[0,\tau_{\mbox{\scriptsize max}}]$, there is infinite number of approximate solutions which reproduce GF within some range of deviations and the problem is to chose “the best one”. Another problem, which is a stumbling block for decades, is the saw tooth noise instability. It occurs when solution is obtained by a naive method, e.g. by using least-squares approach for minimizing deviation measure

[TABLE]

Here $\tilde{G}_{\bf k}(\tau)$ is obtained from approximate LF $\tilde{L}_{\bf k}(\omega)$ by applying of integral operator $\tilde{G}_{\bf k}(\tau)={\cal F}\left[\tilde{L}_{\bf k}(\omega)\right]$ in (13). Saw tooth instability corrupts LF in the ranges where actual LF is smooth. Fast fluctuations of the solution $\tilde{L}_{\bf k}(\omega)$ often have much larger amplitude than the value of actual LF $L_{\bf k}(\omega)$. Standard tools for saw tooth noise suppression are based on the early 60-es idea of Fillips-Tikhonov regularization method regularization ; regcm ; Phillips ; Tikhonov . A nonlinear functional, which suppresses large derivatives of approximate solution $\tilde{L}_{\bf k}(\omega)$, is added to the linear deviation measure (30). Most popular variant of regularization methods is the Maximal Entropy Method JG .

However, typical LF of a QP in a boson field consists of $\delta$-functional peaks and smooth incoherent continuum with a sharp border MPSS ; RP02 . Hence, suppression of high derivatives, as a general strategy of the regularization method, fails. Moreover, any specific implementation of the regularization method uses predefined mesh in the $\omega$ space, which could be absolutely unacceptable for the case of sharp peaks. If the actual location of a sharp peak is between predefined discrete points, the rest of spectral density can be distorted beyond recognition. Finally, regularization Maximal Entropy approach requires assumption of Gauss distribution of statistic errors in $G_{\bf k}(\tau)$, which might be invalid in some cases JG .

Recently, a Stochastic Optimization (SO) method, which circumvents abovementioned difficulties, was developed MPSS . The idea of the SO method is to generate a large enough number $M$ of statistically independent nonregularized solutions $\{\tilde{L}_{\bf k}^{(s)}(\omega)\},s=1,...,M$, which deviation measures $D^{(s)}$ are smaller than some upper limit $D_{u}$, depending of the statistic noise of the GF $G_{\bf k}(\tau)$. Then, using linearity of expressions (13), (30), the final solution is found as the average of particular solutions $\{\tilde{L}_{\bf k}^{(s)}(\omega)\}$

[TABLE]

Particular solution $\tilde{L}_{\bf k}^{(s)}(\omega)$ is parameterized in terms of sum

[TABLE]

of rectangles $\{P_{t}\}=\{h_{t},w_{t},c_{t}\}$ with height $h_{t}>0$, width $w_{t}>0$, and center $c_{t}$. Configuration

[TABLE]

which satisfies normalization condition $\sum_{t=1}^{K}h_{t}w_{t}=1$, defines function $\tilde{G}_{{\bf k}}(\tau)$. The procedure of generating particular solution starts from stochastic choice of initial configuration ${\cal C}_{s}^{\mbox{\scriptsize init}}$. Then, deviation measure is optimized by a randomly chosen consequence of updates until deviation is less than $D_{u}$. In addition to updates, which do not change number of terms in the sum (32), there are updates which increase or decrease number $K$. Hence, since the number of elements $K$ is not fixed, any spectral function can be reproduced with desired accuracy.

Although each particular solution $\tilde{L}_{\bf k}^{(s)}(\omega)$ suffers from saw tooth noise at the area of smooth LF, statistical independence of each solution leads to self averaging of this noise in the sum (32). Note that suppression of noise happens without suppression of high derivatives and, hence, sharp peaks and edges are not smeared out in contrast to regularization approaches. Therefore, saw tooth noise instability is defeated without corruption of sharp peaks and edges. Moreover, continuous parameterization (32) does not need predefined mesh in $\omega$-space. Besides, since the Hilbert space of solution is sampled directly, any assumption about distribution of statistical errors is not necessary.

SO method was successfully applied to restore LF of Fröhlich polaron MPSS , Rashba-Pekar exciton-polaron RP02 , hole-polaron in $t$-$J$-model Hole01 ; tJpho , and many-particle spin system Aplesnin03 . Calculation of the optical conductivity of polaron by SO method can be found in Optics . SO method appeared to be helpful in cases when GF’s asymptotic limit, giving information about ground state, can not be reached. For example, sign fluctuations of the terms in expansion (25) for a hole in the $t$-$J$-model lead to poor statistics at large times Hole01 , though, SO method is capable of recovering energy and $Z$-factor even from GF known only at small imaginary times Hole01 .

2.3 Advantages and Drawbacks of DMC-SO Machinery

Among numerical methods, capable of obtaining quantitative results in the problem of exciton (1) and (2), one can list time-dependent density functional theory Onida02 , Hanke-Sham technique of correcting particle-hole excitation energy HanSham65 ; Benedict98 , and approaches directly solving Bethe-Salpeter equation Alb98 ; RohLo98 ; Marini03 . The latter ones provide rather accurate information on the two-particle GF. However, usage of finite mesh in direct/reciprocal space, which is avoided in DMC method, leads to its’ failure in Wannier regime RohLo98 .

In contrast to DMC method, none of the traditional numeric methods can give reliable results for measurable properties of excited states of polaron at arbitrary range of electron-phonon interaction for the macroscopic system in the thermodynamic limit. Exact diagonalization method Ste96 ; WelFe97 ; FeLoWe97 ; FeLoWe97a can study excited states though only on rather small finite size systems and results of this method are not even justified in the variational sense in the thermodynamic limit BoTruBa99 . There is a batch of rather effective variational “exact translation” methods BoTruBa99 ; ChungTrug ; ShaTrug ; Barisic02 ; Barisic04 where basis is chosen in the momentum space and, hence, the variational principle is applied in the thermodynamic limit. Although these methods can reveal few discrete excited states, its fail for long-range interaction and for dispersive, especially acoustic phonons due to catastrophic growth of variational basis. A non perturbative theory, which is able to give information about spectral properties in the thermodynamic limit at least for one electron, is Dynamical Mean Field Theory DMF1 ; DMF2 ; DMF3 ; DMF4 . However it gives an exact solution only in the case of infinite dimension which does not correspond to a realistic system and can be considered only as a guide for extrapolation to finite dimensions Chuki97 .

Recently developed cluster perturbation theory, where exact diagonalization of a cluster is further improved by taking into account inter-cluster interaction Sene ; Sene1 ; Hohen ; Hohen1 ; Hohen2 , is applicable for study of the excited states, but limited to one-dimensional lattices or two-dimensional systems with short-range interaction. Traditional density-matrix renormalization group method Jeckel ; Jeckel1 ; Jeckel2 ; Jeckel3 ; Jeckel4 is very effective though mostly limited to one-dimensional systems and ladders. Finally, recently developed path integral quantum Monte Carlo algorithm Kornilo ; Kornilo1 ; Kornilo2 ; HagueKonilo06 is valid for any dimension and properly takes into account quasi long-range interactions AlexKornilo99 . Path integral method is capable of obtaining the density of states Kornilo ; Kornilo1 and isotope exponents Kornilo2 ; KorniloAlex04 . However calculations of measurable characteristics of excited states, such as ARPES or optical conductivity, by this method were never reported.

In conclusion, none of methods, except DMC-SO combination, can obtain at the moment approximation-free results for measurable physical quantities for a few QPs interacting with a macroscopic bosonic bath in the thermodynamic limit. Indeed, there are limitations of the DMC and SO methods. DMC method does not work in many-fermion systems due to sign problem and SO method fails at high temperatures, comparable to energies of dominant spectral peaks, because even very small statistical noise of GFs turns Fredholm equation (13) into essentially “ill defined” problem regularization .

3 Spectral Properties of the Fröhlich Polaron

Before development of DMC-SO methods, the information on the excited states of polaron models, especially the Fröhlich one, was very limited. Knowledge of LF was based on results of infinite-dimensions approximation CiPaFe95 , exact diagonalization AleKaRay94 ; WelFe97 ; FeLoWe97 ; FeLoWe97 , or strong coupling expansion AlexRan . No one of the above techniques was capable of obtaining the LF of polaron without approximations, especially for long-range interaction where difficulties of traditional numerical methods dramatically increase. In a similar way, optical conductivity (OC) of Fröhlich model was known only in strong coupling expansion approximation LandPekar , within the framework of the perturbation theory Gurevich , or was based on the variational Feynman path integral technique Dev72 . In this sect. we consider exact DMC-SO results on LF MPSS and OC Optics ; FrCondon of Fröhlich polaron model.

3.1 Lehmann Function of the Fröhlich Polaron

The perturbation theory expression for the high-energy part ($\omega>0$) of the LF for arbitrary interaction potential $V(\mid{\bf q}\mid)$ reads MPSS (frequency of the optical phonon $\omega_{ph}$ is set to unity)

[TABLE]

Low-energy part of the LF for the short-range interaction $V({\mid\bf q}\mid)\,=\,i\,\left(2\sqrt{2}\alpha\pi\right)^{1/2}\,(q^{2}+\kappa^{2})^{-1/2}$ , reducing to the Fröhlich one when $\kappa\to 0$, is

[TABLE]

Comparison of low-energy parts of the LF of the Fröhlich model, obtained by DMC-SO and taken from (35), shows perfect agreement for $\alpha=0.05$: the accuracy for the polaron energy and Z-factor is about $10^{-4}$. On the other hand, high-energy part of numeric result (Fig. 2) significantly deviates from that of the analytic expression (35). This is not surprising since for Fröhlich polaron the perturbation theory expression is diverging as $\omega\to\omega_{ph}$ and, therefore the perturbation theory breaks down. When perturbation theory is obviously valid, e.g. for the case of finite $\kappa=1$, there is a perfect agreement between analytic expression and DMC-SO results (Fig. 2b). Note that the high-energy part of $L_{{\bf k}=0}(\omega)$ is successfully restored by SO method despite the fact that the total weight of the feature for $\alpha=0.05$ is less than $10^{-2}$.

The main deviation of the actual LF from the perturbation theory result is the extra broad peak in the actual LF at $\omega\sim 3.5$. To study this feature $L_{{\bf k}=0}(\omega)$ was calculated for $\alpha=0.5$, $\alpha=1$, and $\alpha=2$ (Fig. 2c-e). The peak is seen for higher values of the interaction constant and its weight grows with $\alpha$. Near the threshold, $\omega=1$, LF demonstrates the square-root dependence $\sim\sqrt{\omega-1}$ (Fig. 2f).

To trace the evolution of the peak at higher values of $\alpha$ the LF was calculated MPSS for $\alpha=4$, $\alpha=6$, and $\alpha=8$ (Fig. 3). At $\alpha=4$ the peak at $\omega\sim 4$ already dominates. Moreover, a distinct high-energy shoulder appears at $\alpha=4$, which transforms into a broad peak at $\omega\sim 8.5$ in the LF for $\alpha=6$. The LF for $\alpha=8$ demonstrates further redistribution of the spectral weight between different maxima without significant shift of the peak positions.

3.2 Optical Conductivity of the Fröhlich Polaron: Validity

of the Franck-Condon Principle in the Optical Spectroscopy

The FC principle Franck ; Condon and its validity have been widely discussed in studies of optical transitions in atoms, molecules Bertran ; Urbain , and solids Lax ; Toyobook . Generally, the FC principle means that if only one of two coupled subsystems, e.g. an electronic subsystem, is affected by an external perturbation, the second subsystem, e.g., the lattice, is not fast enough to follow the reconstruction of the electronic configuration. It is clear that the justification for the FC principle is the short characteristic time of the measurement process $\tau_{mp}\ll\tau_{ic}$, where $\tau_{mp}$ is related to the energy gap between the initial and final states, $\Delta E$, through the uncertainty principle: $\tau_{mp}\simeq\hbar/(\Delta E)$ and $\tau_{ic}$ is the time necessary to adjust the lattice when the electronic component is perturbed. Then, the spectroscopic response considerably depends on the value of the ratio $\tau_{mp}/\tau_{ic}$ For example, in mixed valence systems, where the ionic valence fluctuates between configurations $f^{5}$ and $f^{6}$ with characteristic time $\tau_{ic}\approx 10^{-13}$s, spectra of fast and slow experiments are dramatically different Khomsky ; Falicov . Photoemission experiments with short characteristic times $\tau_{mp}\approx 10^{-16}$s (FC regime), reveal two lines, corresponding to $f^{5}$ and $f^{6}$ states. On the other hand slow Mössbauer isomer shift measurements with $\tau_{mp}\approx 10^{-9}$s show a single broad peak with mean frequency lying between signals from pure $f^{5}$ and $f^{6}$ shells. Finally, according to paradigm of measurement process time, magnetic neutron scattering with $\tau_{mp}\approx\tau_{ic}$ revealed both coherent lines with all subsystems dynamically adjusted and broad incoherent remnants of strongly damped excitation of $f^{5}$ and $f^{6}$ shells Aleks95 ; KiMi95 . Actually, the meaning of the times $\tau_{ic}$ and $\tau_{mp}$ varies with the system and with the measurement process.

To study the interplay between measurement process time $\tau_{mp}$ and adjustment time $\tau_{ic}$, the OC of the Fröhlich polaron was studied in FrCondon from the weak to the strong coupling regime by three methods. DMC method gives numerically exact answer which is compared with memory function formalism (MFF), which is able to take dynamical lattice relaxation into account, and strong coupling expansion (SCE) which assumes FC approach. It was found that near critical coupling $\alpha_{c}\approx 8.5$ a dramatic change of the OC spectrum occurs: dominating peak of OC splits into two satellites. In this critical regime the upper (lower) one quickly decreases (increases) it’s spectral weight as the value of coupling constant increases. Besides, while OC follows prediction of MFF at $\alpha<\alpha_{c}$, its dependence switches to that predicted by SCE for larger couplings. It was concluded that, for the OC measurement of polaron, the adjustment time $\tau_{ic}\approx\hbar/{\cal D}$ is set by typical nonadiabatic energy ${\cal D}$. Nonadiabaticity destroys FC classification at $\alpha<\alpha_{c}$ while FC principle rapidly regains its validity at large couplings due to fast growth of energy separation between initial and final states of optical transitions.

Comparison of exact DMC-SO data for OC with existing results of approximate methods showed Optics that the Feynman path integral technique Dev72 of Devreese, De Sitter, and Goovaerts, where OC is calculated starting from the Feynman variational model Fey55 , is the only successfully describing evolution of the energy of the main peak in OC with coupling constant $\alpha$ (see D ). However, starting from the intermediate coupling regime this approach fails to reproduce the peak width. Subsequently, the path integral approach was rewritten in terms of MFF PeetersDev83 . Then, in FrCondon the extended MFF formalism, which introduces dissipation processes fixed by exact sum rules, was developed Cataudella .

As shown in Fig. 4a, in the weak coupling regime, the MFF, with or without dissipation, is in very good agreement with DMC data, showing significant improvement with respect to weak coupling perturbation approach Gurevich which provides a good description of OC spectra only for very small values of $\alpha$. For $1\leq\alpha\leq 8$, where standard MFF fails to reproduce peak width (Fig. 4b-d) and even the peak position (Fig. 4c), the damping, introduced to extended MFF scheme, becomes crucial. Results of extended MFF are accurate for the peak energy and quite satisfactory for the peak width (Fig. 4b-e). Note that the broadening of the peak in DMC data is not a consequence of poor quality of analytic continuation procedure since DMC-SO methods is capable of revealing such fine features as 2- and 3-phonon thresholds of emission (Fig. 4b).

However, a dramatic change of OC occurs around critical coupling strength $\alpha_{c}\approx 8.5$. The dominating peak of OC splits into two ones, the energy of lover one corresponding to the predictions of SCR expansion and that of upper one obeying extended MFF value (Fig. 5a). The shoulder, corresponding to dynamical extended MFF contribution, rapidly decreases it’s intensity with increase of $\alpha$ and at large $\alpha$ (Fig. 5b-c) the OC is in good agreement with strong coupling expansion, assuming FC scheme. Finally, comparing energies of the peaks, obtained by DMC, extended MFF and FC strong coupling expansion (Fig. 5d), we conclude that at critical coupling $\alpha_{c}\approx 8.5$ the spectral properties rapidly switch from dynamic, when lattice relaxes at transition, to FC regime, where nuclei are frozen in initial configuration.

In order get an idea of the FC breakdown authors of FrCondon consider the following arguments. The approximate adiabatic states are not exact eigenstates of the system. These states are mixed by nondiagonal matrix elements of the nonadiabatic operator ${\cal D}$ and exact eigenstates are linear combinations of the adiabatic wavefunctions. Being interested in the properties of transition from ground ($g$) to an excited ($ex$) state, whose energy correspond to that of the OC peak, it is necessary to consider mixing of only these states and express exact wavefunctions as a linear combinations Brovman ; KiMi of ground and excited adiabatic states. The coefficients of superposition are determined from standard techniques Brovman ; KiMi where nondiagonal matrix elements of the nonadiabatic operator Brovman are expressed in terms of matrix elements of the kinetic energy operator $M$, the gap between excited and ground state $\Delta E=E_{ex}-E_{g}$ and the number $n_{\beta}$ of phonons in adiabatic state:

[TABLE]

The extent to which lattice can follow transition between electronic states, depends on the degree of mixing between initial and final exact eigenstates through the nonadiabatic interaction. If initial and ground states are strongly mixed, the adiabatic classification has no sense and, therefore, the FC processes have no place and lattice is adjusted to the change of electronic states during the transition. In the opposite limit adiabatic approximation is valid and FC processes dominate. The estimation for the weight of FC component $I_{FC}$ FrCondon is equal to unity in the case of zero mixing and zero in the case of maximal mixing. The weight of adiabatically connected (AC) transition $I_{AC}=1-I_{FC}$ is defined accordingly. Non-diagonal matrix element $M$ is proportional to the root square of $\alpha$ with a coefficient $\eta$ of the order of unity. In the strong coupling regime, assuming that $\Delta E\approx\gamma\alpha^{2}$ ($\gamma\approx 0.1$ from MC data), and $n_{\beta}\approx\Delta E$ ($n_{\beta}\gg 1$), one gets

[TABLE]

where $\tau_{mp}=1/\Delta E$ and $\tau_{ic}=1/D$. For $\eta$ of the order of unity one obtains qualitative description of a rather fast transition from AC- to FC-dominated transition, when $I_{FC}$ and $I_{AC}$ exchange half of their weights in the range of $\alpha$ from 7 to 9. The physical reason for such quick change is the faster growth of energy separation $\Delta E\sim\alpha^{2}$ compared to that of the matrix element $M\sim\alpha^{1/2}$. Finally, for large couplings, initial and final states become adiabatically disconnected. The rapid AC-FC switch has nothing to do with the self-trapping phenomenon where crossing and hybridization of the ground and an excited states occurs. This phenomenon is a property of transition between different states and related to the choice whether lattice can or can not follow adiabatically the change of electronic state at the transition.

4 Self-Trapping

In this section we consider the self-trapping (ST) phenomenon which, due to essential importance of many-particle interaction of QP with bosonic bath of macroscopic system, was never addressed by exact method before. We start with a basic definition of the ST phenomenon and introduce the adopted criterion for it’s existence. Then, generic features of ST are demonstrated on a simple model of Rashba-Pekar exciton-polaron in Sect. 4.1. It is shown in Sect. 4.2 that the criterion is not a dogma since even in one dimensional system, where ST is forbidden by criterion of existence, one can observe all main features of ST due to peculiar nature of electronic states.

In general terms UKKTH86 ; R82 , ST is a dramatic transformation of a QP properties when system parameters are slightly changed. The physical reason of ST is a quantum resonance, which happens at some critical interaction constant $\gamma_{c}$, between “trapped” (T) state of QP with strong lattice deformation around it and “free” (F) state. Naturally, ST transition is not abrupt because of nonadiabatic interaction between T and F states and all properties of the QP are analytic in $\gamma$ Gerlach . At small $\gamma<\gamma_{c}$, ground state is an F state which is weakly coupled to phonons while excited states are T states and have a large lattice deformation. At critical couplings $\gamma\approx\gamma_{c}$ a crossover and hybridization of these states occurs. Then, for $\gamma>\gamma_{c}$ the roles of the states exchange. The lowest state is a T state while the upper one is an F state.

First, and up to now the only quantitative criterion for ST existence was given in terms of the ground state properties in the adiabatic approximation. This criterion considers stability of the delocalized state in undistorted lattice $\Delta=0$ with respect to the energy gain due to lattice distortion $\Delta^{\prime}\neq 0$. ST phenomenon occurs when completely delocalized state with $\Delta=0$ is separated from distorted state with $\Delta^{\prime}\neq 0$ by a barrier of adiabatic potential. One of these states is stable while another one is meta-stable. The criterion of barrier existence is defined in terms of the stability index

[TABLE]

where $d$ is the system dimensionality. Index $l$ determines the range of the force $\lim_{q\to 0}\psi(q)\sim q^{-l}$, where $\psi({\bf R})$ is the kernel of interaction $U({\bf R}_{n})=\psi({\bf R}_{n}-{\bf R}_{n^{\prime}})\nu({\bf R}_{n^{\prime}})$ connecting potential $U({\bf R}_{n})$ with generalized lattice distortion $\nu({\bf R}_{n^{\prime}})$ UKKTH86 . The barrier exists for $s>0$ and does not exist for $s<0$. The discontinuous change of the polaron state, i.e. ST, occurs in the former case while does not happen in the latter case. When $s=0$, this scaling argument alone can not conclude the presence or absence of the ST and more detailed discussion for each model is needed.

4.1 Typical Example of the Self-Trapping: Rasba-Pekar Exciton-Polaron

Classical example of a system with ST phenomenon is the three dimensional continuous Rasba-Pekar exciton-polaron in the approximation of intraband scattering, i.e. when polar electron-phonon interaction (EPI) with dispersionless optical phonons $\omega_{\mbox{\scriptsize ph}}=1$ does not change the wave function of internal electron-hole motion. System is defined as a structureless QP with dispersion $\epsilon({\bf k})=k^{2}/2$ and short range coupling to phonons RP02 ; UKKTH86 . General criterion of the existence of ST is satisfied for three dimensional system with short range interaction RP02 ; UKKTH86 ; UFN_05 and, thus, one expects to observe typical features of the phenomenon.

It is shown RP02 that in the vicinity of the critical coupling $\gamma_{c}\approx 18$ the average number of phonons $\langle N\rangle$ in (18) and effective mass $m^{*}$ quickly increase in the ground state by several orders of magnitude (Fig. 6b-c). Besides, a quantum resonance between polaronic phonon clouds of F and T state is demonstrated. Distribution of partial $n$-phonon contributions $Z^{({\bf k}=0)}(n)$ in (17) has one maximum at $n=0$ in the weak coupling regime, which corresponds to weak deformation, and one maximum at $n\gg 1$ in the strong coupling regime, which is the consequence of a strong lattice distortion. However, due to F-T resonance there are two distinct peaks at $n=0$ and $n\gg 1$ for $\gamma\approx\gamma_{c}$ (Fig. 6d).

Near the critical coupling $\gamma_{c}$ the LF of polaron has several stable states (Fig. 7 a-b) below the threshold of incoherent continuum $E_{\mbox{\scriptsize gs}}+\omega_{\mbox{\scriptsize ph}}$. Any state above the threshold is unstable because emission of a phonon with transition to the ground state at ${\bf k}=0$ with energy $E_{\mbox{\scriptsize gs}}$ is allowed. On the other hand, decay is forbidden by conservation laws for states below the threshold. Dependence of the energies of ground and excited resonances on the interaction constant resembles a picture of crossing of several states interacting with each other (Fig. 7c).

According to the general picture of the ST phenomenon, lowest F state in the weak coupling regime at ${\bf k}=0$ has small effective mass $m^{*}\approx m$ of the order of the bare QP mass $m$. To the contrary, the effective mass of excited state $m^{*}\gg m$ is large. Hence, below the critical coupling the energy of the F state, which is lowest at ${\bf k}=0$, has to reach a flat band of T state at some momentum. Then, F and T state have to hybridize and exchange in energy. DMC data visualize this picture (Fig. 7 d-e). After F state crosses the flat band of excited T state, the average number of phonons increases and dispersion becomes flat.

It is natural to assume that above the critical coupling the situation is opposite: ground state is the T state with large effective mass while excited F state has small, nearly bare, effective mass. Indeed, this assumption was confirmed in the framework of another model which is considered in Sect. 6.1. Moreover, it was shown that in the strong coupling regime excited resonance inherits not only bare effective mass around ${\bf k}=0$ but the whole dispersion law of the bare QP tJpho .

4.2 Degeneracy Driven Self-Trapping

According to the criterion (38), ST phenomenon in one-dimensional system does not occur. Although this statement is probably valid for the case of single band in relevant energy range, it is not the case for the generic multi-band cases. This fact has been unnoticed for many years, which prevented the proper explanation of puzzling physics of quasi-one-dimensional compound Anthracene-PMDA, although it’s optical properties HaPhi75 ; BrPh80 ; Ha74 ; Ha77 ; ElsWeis85 ; KuGo94 directly suggested resonance of T and F states. The reason is that in Anthracene-PMDA, in contrast to conditions at which criterion (38) is obtained, there are two, nearly degenerate exciton bands. Then, one can consider quasi-degenerate self-trapping mechanism when ST phenomenon is driven by nondiagonal interaction of phonons with quasidegenerate exciton levels MN01 . Such mechanism was already suggested for explanation of properties of mixed valence systems KiMi though it’s relevance was never proved by an exact approach.

The minimal model to demonstrate the mechanism of quasi-degenerate self-trapping involves one optical phonon branch with frequency $\omega_{ph}=0.1$ and two exciton branches with energies $\epsilon_{1,2}(q)=\Delta_{1,2}+2[1-\cos(q)]$, where $\Delta_{1}=0$ and $\Delta_{2}=1$. Presence of short range diagonal $\gamma_{22}$ and nondiagonal $\gamma_{12}$ interactions (with corresponding dimensionless constants $\lambda_{22}=\gamma_{22}^{2}/(2\omega)$ and $\lambda_{12}=\gamma_{12}^{2}/(2\omega)$) leads to classical self-trapping behavior even in one-dimensional system MN01 (see Fig. 8).

5 Exciton

Despite numerous efforts over the years, there has been no rigorous technique to solve for exciton properties even for the simplest model (1)-(2) which treats electron-electron interactions as a static renormalized Coulomb potential with averaged dynamical screening. The only solvable cases are the Frenkel small-radius limit Frenkel and the Wannier large-radius limit Wannier which describe molecular crystals and wide gap insulators with large dielectric constant, respectively. Meanwhile, even the accurate data for the limits of validity of the Wannier and Frenkel approximations have not been available. As discussed in Sects. 1.2 and 2.3, semianalytic approaches has little to add to problem when quantitative results are needed whereas traditional numerical methods fail to reproduce them even in the Wannier regime. To the contrary, DMC results do not contain any approximation.

To study conditions of validity of limiting regimes by DMC method, electron-hole spectrum of three dimensional system was chosen in the form of symmetric valence and conduction bands with width $E_{c}$ and direct gap $E_{g}$ at zero momentum Exciton . For large ratio $W=E_{c}/E_{g}$, when $W>30$, exciton binding energy is in good agreement with Wannier approximation results (Fig. 9a) and probability density of relative electron-hole motion corresponds (Fig. 9c) to hydrogen-like result. The striking result is the requirement of rather large valence and conduction bandwidths ($W>20$) for applicability of Wannier approximation. For smaller values of $W$ the binding energy and wave function of relative motion (Fig. 9d) deviate from large radius results. In the similar way, conditions of validity of Frenkel approach are rather restricted too. Moreover, even strong localization of wave function does not guarantee good agreement between exact and Frenkel approximation result for binding energy. At $1<W<10$ the wave function is already strongly localized though binding energy considerably differs from Frenkel approximation result. For example, at $W=0.4$ relative motion is well localized (Fig. 9e) whereas the binding energy of Frenkel approximation is two times larger than exact result (Inset in Fig. 9a).

A study of conditions necessary for formation of charge transfer exciton in three dimensional systems is crucial to finalize protracted discussion of numerous models concerning properties of mixed valence semiconductors CuKiMi . A decade ago unusual properties of SmS and SmB6 were explained by invoking the excitonic instability mechanism assuming charge-transfer nature of the optically forbidden exciton KiMiZhe88 ; KiMi90 . Although this model explained quantitatively the phonon spectra Aleks89 ; KiMi91 , optical properties TraWa84 ; Le95 , and magnetic neutron scattering data KiMi95 , it’s basic assumption has been criticized as being groundless Kas94_1 ; Kas94_2 . To study excitonic wavefunction, dispersions of valence and conduction bands were chosen as it is typical for mixed valence materials: almost flat valence band is separated from broad conduction band, having maximum in the centre and minimum at the border of Brillouin zone Exciton . Results presented in Fig. 9b support assumption of KiMiZhe88 ; KiMi90 since wave function of relative motion has almost zero on-site component and maximal charge density at near neighbors.

6 Polarons in Undoped High Temperature Superconductors

It is now well established that the physics of high temperature superconductors is that of hole doping a Mott insulator ManoRev94 ; Dag94 ; Lee06 . Even a single hole in a Mott insulator, i.e. a hole in an antiferromagnet in case of infinite Hubbard repulsion $U$, is substantially influenced by many-body effects Kane because it’s jump to a neighboring site disturbs antiferromagnetic arrangement of spins. Hence a thorough understanding of the dynamics of doped holes in Mott insulators has attracted a great deal of recent interest. The two major interactions relevant to the electrons in solids are electron-electron interactions (EEI) and electron-phonon interactions (EPI). The importance of the former at low doping is no doubt essential since the Mott insulator is driven by strong Hubbard repulsion, while the latter was considered to be largely irrelevant to superconductivity based on the observations of a small isotope effect on the optimal $T_{c}$ isotope and an absence of a phonon contribution to the resistivity (for review see resistivity ).

On the other hand, there are now accumulating evidences that the EPI plays an important role in the physics of cuprates such as (i) an isotope effect on superfluid density $\rho_{s}$ and Tc away from optimal doping Keller , (ii) neutron and Raman scattering Pint99 ; Raman ; Raman1 experiments showing strong phonon softening with both temperature and hole doping, indicating that EPI is strong Khal97 ; Gunnars2 . Furthermore, the recent studies of cuprates by the angle resolved photoemission spectroscopy (ARPES), which spectra are proportional to the LF (7) Shen_03 , resulted in the discovery of the dispersion ”kinks” at around 40-70meV measured from the Fermi energy, in the correct range of the relevant oxygen related phonons Lanzara01 ; Cuk2004 ; Dev2004 . These particular phonons - oxygen buckling and half-breathing modes are known to soften with doping Egami ; Pint99 and with temperature Cuk2004 ; Dev2004 ; Egami ; Pint99 ; Raman ; Raman1 indicating strong coupling. The quick change of the velocity can be predicted by any interaction of a quasiparticle with a bosonic mode, either with a phonon Cuk2004 ; Dev2004 or with a collective magnetic resonance mode ChuNor ; EschNor ; EschNor1 . However, the recently discovered “universality” of the kink energy for LSCO over the entire doping range Zhou03 casts doubts on the validity of the latter scenario as the energy scale of the magnetic excitation changes strongly with doping.

Besides, measured in undoped high $T_{c}$ materials ARPES revealed apparent contradiction between momentum dependence of the energy and linewidth of the QP peak. On the one hand the experimental energy dispersion of the broad peak in many underdoped compounds ZX95 ; KyleDop obeys the theoretical predictions Xiang_96 ; Kyu_Fer , whereas the experimental peak width is comparable with the bandwidth and orders of magnitude larger than that obtained from theory of Mott insulator Hole01 . Early attempts to interpret this anomalously short lifetime of a hole by an interaction with additional nonmagnetic bosonic excitations, e.g. phonons Pothuizen97 , faced generic question: is it possible that interaction with media leaves the energy dispersion absolutely unrenormalized, while, induces a decay which inverse life-time is comparable or even larger than the QP energy dispersion? A possibility of an extrinsic origin of this width can be ruled out since the doping induces further disorder, while a sharper peak is observed in the overdoped region.

In order to understand whether phonons can be responsible for peculiar shape of the ARPES in the undoped cuprates, the LF of an interacting with phonons hole in Mott insulator was studied by DMC-SO tJpho . The case of the LF of a single hole corresponds to the ARPES in an undoped compound. For a system with large Hubbard repulsion $U$, when $U$ is much larger than the typical bandwidth $W$ of noninteracting QP, the problem reduces to the $t$-$J$ model ChaSpaOles77 ; GJRice87 ; ManoRev94 ; Izymov97

[TABLE]

Here $c_{j\sigma}$ is projected (to avoid double occupancy) fermion annihilation operator, $n_{i}$ ($<2$) is the occupation number, ${\bf S}_{i}$ is spin 1/2 operator, $J$ is an exchange integral, and $\langle ij\rangle$ denotes nearest-neighbor sites in two dimensional square lattice. Different theoretical approaches revealed ManoRev94 ; Assad00 ; Hole01 basic properties of the LF. The LF has a sharp peak in the low energy part of the spectrum which disperses with a bandwidth $W_{J/t}\sim 2J$ and, therefore, the large QP width in experiment can not be explained. More complicated $tt^{\prime}t^{\prime\prime}$-$J$ model takes into account hoppings to the second $t^{\prime}$ and third $t^{\prime\prime}$ nearest neighbors and, hence, dispersion of the hole changes Chernyshov96_1 ; Chernyshov96_2 ; Tohyama_00 ; Xiang_96 ; Kyu_Fer ; Shen_03 . However, for parameters, which are necessary for description of dispersion in realistic high $T_{c}$ superconductors ZX95 ; Xiang_96 , peak in the low energy part remains sharp and well defined for all momenta Zaanen_95 .

After expressing spin operators in terms of Holstein-Primakoff spin wave operators and diagonalizing the spin part of Hamiltonian (39) by Fourier and Bogoliubov transformations SVR88 ; Kane ; Liu_91 ; Liu_92 , $tt^{\prime}t^{\prime\prime}$-$J$ Hamiltonian is reduced to the boson-holon model, where hole (annihilation operator is $h_{\bf k}$) with dispersion $\varepsilon({\bf k})=4t^{\prime}\cos(k_{x})\cos(k_{y})+2t^{\prime\prime}[\cos(2k_{x})+\cos(2k_{y})]$ propagates in the magnon (annihilation operator is $\alpha_{\bf k}$) bath

[TABLE]

with magnon dispersion $\omega_{\bf k}=2J\sqrt{1-\gamma_{\bf k}^{2}}$, where $\gamma_{\bf k}=(\cos k_{x}+\cos k_{y})/2$. The hole is scattered by magnons as described by

[TABLE]

with the scattering vertex $M_{\bf k,q}$. Parameters $t$, $t^{\prime}$ and $t^{\prime\prime}$ are hopping amplitudes to the first, second and third near neighbors, respectively. If hopping integrals $t^{\prime}$ and $t^{\prime\prime}$ are set to zero and bare hole has no dispersion, the problem (40-41) corresponds to $t$-$J$ model.

Short range interaction of a hole with dispersionless optical phonons $\hat{H}^{\mbox{\scriptsize e-ph}}=\Omega_{0}\sum_{\bf k}b_{\bf k}^{\dagger}b_{\bf k}$ of the frequency $\Omega_{0}$ is introduced by Holstein Hamiltonian

[TABLE]

where $\sigma$ is the momentum and isotope independent coupling constant, $M$ is the mass of the vibrating lattice ions, and $\Omega_{0}$ is the frequency of dispersionless phonon. The coefficient in front of square brackets is the standard Holstein interaction constant $\gamma=\sigma/\sqrt{(}2M\Omega_{0})$. In the following we characterize strength of EPI in terms of dimensionless coupling constant $\lambda=\gamma^{2}/4t\Omega_{0}$. Note, if interaction with magnetic subsystem (41) is neglected and hole dispersion $\varepsilon({\bf k})$ is chosen in the form $\varepsilon({\bf k})=2t[\cos(k_{x})+\cos(k_{y})]$, the problem (40), (42) corresponds to standard Holstein model where hole with near neighbor hopping amplitude $t$ interacts with dispersionless phonons.

We consider the evolution of ARPES of a single hole in $t$-$J$-Holstein model (40)-(42) from the weak to the strong coupling regime and dispersion of the LF in the strong coupling regime in Sect. 6.1. It occurs that properties of the LF in the strong coupling regime of the EPI explain the puzzle of broad lineshape in ARPES in underdoped high $T_{c}$ superconductors. Therefore, in order to suggest a crucial test for the mechanism of phonon-induced broadening, we present calculations of the effect of the isotope substitution on the ARPES in Sect. 6.2.

6.1 Spectral Function of a Hole Interacting with Phonons

in the $t$-$J$ Model: Self-Trapping and Momentum Dependence

Previously, the LF of $t$-$J$-Holstein model was studied by exact diagonalization method on small clusters Bauml_98 and in the non-crossing approximation (NCA)333NCA is equivalent to self-consistent Born approximation (SCBA) for both phonons and magnons Ramsak_92 ; Kyung_96 . However, the small system size in exact diagonalization method implies a discrete spectrum and, therefore, the problem of lineshape could not be addressed. The latter method omits the FDs with mutual crossing of phonon propagators and, hence, is an invalid approximation for phonons in strong and intermediate couplings of EPI. This statement was demonstrated by DMC, which can sum all FDs for Holstein model both exactly and in the NCA tJpho . Exact results and those of NCA are in good agreement for small values $\lambda\leq 0.4$ and drastically different for $\lambda>1$. For example, for $\Omega_{0}/t=0.1$ exact result shows a sharp crossover to strong coupling regime for $\lambda>\lambda_{H}^{c}\approx 1.2$ whereas NCA result does not undergo such crossover even at $\lambda=100$. On the other hand, NCA is valid for interaction of a hole with magnons since spin S=1/2 can not flip more than once and number of magnons in the polaronic cloud can not be large. Note that the $t$-$J$-Holstein model is reduced to problem of polaron which interacts with several bosonic fields (3)-(4).

DMC expansion in tJpho takes into account mutual crossing of phonon propagators and, in the framework of partial NCA, neglects mutual crossing of magnon propagators, to avoid sign problem. NCA for magnons is justified for $J/t\leq 0.4$ by good agreement of results of NCA and exact diagonalization on small clusters SVR88 ; Kane ; Marsiglio91 ; MartHor_91 ; Liu_92 . Recently results of exact diagonalization were compared in the limit of small EPI for $t$-$J$-Holstein model, boson-holon model (40-42) without NCA, and boson-holon model with NCA RoGu_SCBA . Although agreement is not so good as for pure $t$-$J$ model, it was concluded that NCA for magnons is still good enough to suggest that one can use NCA for a qualitative description of the $t$-$J$-Holstein model.

Figures 10a-e show low energy part of LF in the ground state at ${\bf k}=(\pi/2,\pi/2)$ in the weak, intermediate, and strong coupling regimes of interaction with phonons. Dependence on the coupling constant of energies of resonances (Fig. 10f), $Z^{{\bf k}=(\pi/2,\pi/2)}$-factor of lowest peak (Fig. 10g), and average number of phonons in the polaronic cloud $\langle N\rangle$ (Fig. 10h) demonstrates a picture which is typical for ST (see R82 ; RP02 and Sect. 4). Two states cross and hybridize in the vicinity of critical coupling constant $\lambda_{\mbox{\scriptsize t-J}}^{c}\approx 0.38$, $Z^{{\bf k}=(\pi/2,\pi/2)}$-factor of lowest resonance sharply drops and average number of phonons in polaronic cloud quickly rises. According to the general understanding of the ST phenomenon, above the critical couplings $\lambda>\lambda_{\mbox{\scriptsize t-J}}^{c}$ one expects that the lowest state is dispersionless while the upper one has small effective mass. This assumption is supported by the momentum dependence of the LF in the strong coupling regime (Fig. 11a-e). Dispersion of upper broad shake-off Franck-Condon peak nearly perfectly obeys relation

[TABLE]

which describes dispersion of the pure $t$-$J$ model in the broad range of exchange constant $0.1<J/t<0.9$ Marsiglio91 (Fig. 11f). Note that this property of the shake-off peak is general for the whole strong coupling regime (Fig. 11f).

Momentum dependence of the shake-off peak, reproducing that of the free particle, is the direct consequence of the adiabatic regime. Actually, phonon frequency $\Omega_{0}$ is much smaller than the coherent bandwidth $2J$ of the $t$-$J$ model, giving the adiabatic ratio $\Omega/2J=1/6\ll 1$. Besides, as experience with the OC of the Fröhlich polaron (Sect. 3.2) shows, there is one more important parameter in the strong coupling limit. Namely, the ratio between measurement process time $\tau_{mp}=\hbar/\Delta E$ where $\Delta E$ is the energy separation of shake-off hump from the ground state pole, and that of characteristic lattice time $\tau\approx 1/\Omega_{0}$ is much less than unity. Hence, fast photoemission probe sees the ions frozen in one of possible configurations Misko97 . The LF in the FC limit is a sum of transitions between a lower $E_{\mbox{\scriptsize low}}(Q)$ and an upper $E_{\mbox{\scriptsize up}}(Q)$ sheets of adiabatic potential, weighted by the adiabatic wave function of the lower sheet $\mid\psi_{\mbox{\scriptsize low}}(Q)\mid^{2}$ RoGu2005 . If EPI is absent both in initial $E_{\mbox{\scriptsize low}}(Q)=Q^{2}/2$ and final $E_{\mbox{\scriptsize up}}(Q)={\cal D}+Q^{2}/2$ states, the LF is peaked at the energy ${\cal D}$. Then, if there is EPI $\Delta E_{\mbox{\scriptsize up}}(Q)=-\lambda Q$ only in the final state, i.e. when hole is removed from the Mott insulator, the upper sheet of adiabatic potential $E_{\mbox{\scriptsize up}}(Q)={\cal D}-\lambda^{2}/2+(Q-\lambda)^{2}/2$ has the same energy ${\cal D}$ at $Q=0$. Since the probability function $\mid\psi_{\mbox{\scriptsize low}}(Q)\mid^{2}$ has maximum at $Q=0$, the peak of the LF broadens but it’s energy does not shift RoGu2005 (Fig. 11g).

Behavior of the LF is the same as observed in the ARPES of undoped cuprates. The LF consists of a broad peak and a high energy incoherent continuum (see Fig. 11a). Besides, dispersion of the broad peak “c” in Figs. 11 reproduces that of sharp peak in pure $t$-$J$ model (Fig. 11b-f). The lowest dispersionless peak, corresponding to small radius polaron, has very small weight and, hence, can not be seen in experiment. On the other hand, according to experiment, momentum dependence of spectral weight $Z^{({\bf k})^{\prime}}$ of broad resonance exactly reproduces dispersion of $Z^{({\bf k})}$-factor of pure $t$-$J$ model. The reason for such perfect mapping is that in adiabatic case $\Omega_{0}/2J\ll 1$ all weight of the sharp resonance in $t$-$J$ model without EPI is transformed at strong EPI into the broad peak. This picture implies that the chemical potential in the heavily underdoped cuprates is not connected with the broad resonance but pinned to the real quasiparticle pole with small $Z$-factor. This conclusion was recently confirmed experimentally KyleDop .

Comparing the critical EPI for a hole in the $t$-$J$-Holstein model (40-42) $\lambda_{\mbox{\scriptsize t-J}}^{c}\approx 0.38$ and that for Holstein model $\lambda_{H}^{c}\approx 1.2$ with the same value of hopping $t$, we conclude that spin-hole interaction accelerates transition into the strong coupling regime. The reason for enhancement of the role of EPI is found in RoGu_SCBA . Comparison of the EPI driven renormalization of the effective mass in $t$-$J$-Holstein and Holstein model shows that large effective mass in the $t$-$J$ model is responsible for this effect. The enhancement of the role of EPI by EEI takes place at least for a single hole at the bottom of the $t$-$J$ band. Had the comparison been made with half-filled model, the result would have been smaller enhancement or no enhancement at all Sangio_Gunn_06 . On the other hand, coupling constant of half-breathing phonon is increased by correlations RoGuPRB_04 . Finally, we conclude that effect of enhancement of the effective EPI by EEI is not unambiguous and depends on details of interaction and filling. However, this effect is present for small filling in the $t$-$J$-Holstein model.

6.2 Isotope Effect on ARPES in Underdoped High-Temperature

Superconductors

The magnetic resonance mode and the phonon modes are the two major candidates to explain the “kink” structure of the electron energy dispersion around 40-70 meV below the Fermi energy, and the isotope effect (IE) on ARPES should be the smoking-gun experiment to distinguish between these two. Gweon et al. Gweon_04 performed the ARPES experiment on O18-replaced Bi2212 at optimal doping and found an appreciable IE, which however can not be explained within the conventional weak-coupling Migdal-Eliashberg theory. Namely the change of the spectral function due to O18-replacement has been observed at higher energy region beyond the phonon energy ($\sim 60$meV). This is in sharp contrast to the weak coupling theory prediction, i.e., the IE should occur only near the phonon energy. Hence the IE in optimal Bi2212 remains still a puzzle. On the other hand, the ARPES in undoped materials, as described in Sect. 6.1, has recently been understood in terms of the small polaron formation tJpho ; RoGuZX ; RoGu2005 . Therefore, it is essential to compare experiment in undoped systems with presented in this Sect. DMC-SO data, where theory can offer quantitative results.

In addition to high-$T_{c}$ problem, strong EPI mechanism of ARPES spectra broadening was considered as one of alternative scenarios for diatomic molecules Zawatzky_89 , colossal magnetoresistive manganites Dessau_98 , quasi-one-dimensional Peierls conductors Perf_01 ; Perf_02 , and Verwey magnetites Schrupp_05 . Therefore, exact analysis of the IE on ARPES at strong EPI is of general interest for conclusive experiments in a broad variety of compound classes.

Dimensionless coupling constant $\lambda=\gamma^{2}/4t\Omega$ in (42) is an invariant quantity for the simplest case of IE. Indeed, assuming natural relation $\Omega\sim 1/\sqrt{M}$ between phonon frequency and mass, we find that $\lambda$ does not depend on the isotope factor $\kappa_{\mbox{\scriptsize iso}}=\Omega/\Omega_{0}=\sqrt{M_{0}/M}$, which is defined as the ratio of phonon frequency in isotope substituted ($\Omega$) and normal ($\Omega_{0}$) systems. We chose adopted parameters of the $tt^{\prime}t^{\prime\prime}$-$J$ model which reproduce the experimental dispersion of ARPES Xiang_96 : $J/t=0.4$, $t^{\prime}/t=-0.34$, and $t^{\prime\prime}/t=0.23$ . The frequency of the relevant phonon Shen_03 is set to $\Omega_{0}/t=0.2$ and the isotope factor $\kappa_{\mbox{\scriptsize iso}}=\sqrt{16/18}$ corresponds to substitution of O18 isotope for O16.

To sweep aside any doubts of possible instabilities of analytic continuation, we calculate the LF for normal compound ($\kappa_{\mbox{\scriptsize nor}}=1$), isotope substituted ($\kappa_{\mbox{\scriptsize iso}}=\sqrt{16/18}$) and “anti-isotope” substituted ($\kappa_{\mbox{\scriptsize ant}}=\sqrt{18/16}$) compounds. Monotonic dependence of LF on $\kappa$ ensures stability of analytic continuation and gives possibility to evaluate the error-bars of a quantity ${\cal A}$ using quantities ${\cal A}_{\mbox{\scriptsize iso}}-{\cal A}_{\mbox{\scriptsize nor}}$, ${\cal A}_{\mbox{\scriptsize nor}}-{\cal A}_{\mbox{\scriptsize ant}}$, and $({\cal A}_{\mbox{\scriptsize iso}}-{\cal A}_{\mbox{\scriptsize ant}})/2$.

Since LF is sensitive to strengths of EPI only for low frequencies Isotope , we concentrate on the low energy part of the spectrum. Figure 12 shows IE on the hole LF for different couplings in nodal and antinodal points, respectively. The general trend is a shift of all spectral features to larger energies with increase of the isotope mass ($\kappa<1$). One can also note that the shift of broad FCP is much larger than that of narrow real-QP peak. Moreover, for large couplings $\lambda$ the shift of QP energy approaches zero and only decrease of QP spectral weight $Z$ is observed for larger isotope mass. On the other hand, the shift of FCP is not suppressed for larger couplings. Except for the LF in nodal point at $\lambda=0.62$ (Fig. 12a, b), where LF still has significant weight of QP $\delta$-functional peak, there is one more notable feature of the IE. With increase of the isotope mass the height of FCP increases. Taking into account the conservation law for LF $\int_{-\infty}^{+\infty}L_{\bf k}(\omega)=1$ and insensitivity of high energy part of LF to EPI strength Isotope , the narrowing of the FCP for larger isotope mass can be concluded. To understand the trends of the IE in the strong coupling regime we analyze the exactly solvable independent oscillators model (IOM) Mahan . The LF in IOM is the Poisson distribution

[TABLE]

where $\xi_{0}=\gamma^{2}_{0}/\Omega_{0}^{2}=4t\lambda/\Omega_{0}$ is dimensionless coupling constant for normal system and ${\cal G}_{\kappa,l}(\omega)=\delta[\omega+4t\lambda-\Omega_{0}\kappa l]$ is the $\delta$-function. The properties of the Poisson distribution quantitatively explain many features of the IE on LF444Cautions should be made about approximate form of EPI (42). Strictly speaking, actual momentum dependence of the interaction constant $\sigma$ RoGu_04 ; Ishihara_04 can slightly change the obtained differences between nodal and antinodal points though the general trends have to be left intact because ST is caused solely by the short range part of EPI R82 ..

The energy $\omega_{\mbox{\scriptsize QP}}=-4t\lambda$ of the zero-phonon line $l=0$ in (44) depends only on isotope independent quantities which explains very weak isotope dependence of QP peak energy in insets of Fig. 12. Besides, change of the zero-phonon line weight $Z^{(0)}$ obeys relation $Z^{(0)}_{\mbox{\scriptsize iso}}/Z^{(0)}_{\mbox{\scriptsize nor}}=\exp\left[-\xi_{0}(1-\kappa)/\kappa\right]$ in IOM. These IOM estimates agree with DMC data within 15% in the nodal point and within 25% in the antinodal one. IE on FCP in the strong coupling regime follows from the properties of zero $M_{0}=\int_{-\infty}^{+\infty}L(\omega)d\!\omega=1$, first $M_{1}=\int_{-\infty}^{+\infty}\omega L(\omega)d\!\omega=0$, and second $M_{2}=\int_{-\infty}^{+\infty}\omega^{2}L(\omega)d\!\omega=\kappa\xi_{0}\Omega_{0}^{2}$ moments of shifted Poisson distribution (44). Moments $M_{0}$ and $M_{2}$ establish relation ${\cal D}=h^{\mbox{\scriptsize FCP}}_{\mbox{\scriptsize iso}}/h^{\mbox{\scriptsize FCP}}_{\mbox{\scriptsize nor}}=1/\sqrt{\kappa}\approx 1.03$ between heights of FCP in normal and substituted compounds. DMC data in the antinodal point perfectly agree with the above estimate for all couplings. This is consistent with the idea that the anti-nodal region remains in the strong coupling regime even though the nodal region is in the crossover region. In the nodal point DMC data well agree with IOM estimate for $\lambda=0.75$ (${\cal D}\approx 1.025$) whereas at $\lambda=0.69$ and $\lambda=0.62$ influence of the ST point leads to anomalous values of ${\cal D}$: ${\cal D}\approx 1.07$ and ${\cal D}\approx 0.98$, respectively. Shift of the low energy edge at half maximum $\Delta_{1/2}$ must be proportional to change of the root square of second moment $\Delta_{\sqrt{M_{2}}}=\sqrt{\xi_{0}}\Omega_{0}[1-\sqrt{\kappa}]$. As we found in numeric simulations of (44) with Gaussian functions555 Results are almost independent on the parameter $\eta$ of the Gaussian distribution ${\cal G}_{\kappa,l}(\omega)=1/(\eta\sqrt{2\pi})\exp(-[\omega+4t\lambda-\Omega_{0}\kappa l]/(2\eta^{2}))$ in the range $[0.12,0.2]$. ${\cal G}_{\kappa,l}(\omega)$, relation $\Delta_{1/2}\approx\Delta_{\sqrt{M_{2}}}/2$ is accurate to 10% for $0.62<\lambda<0.75$. Also, simulations show that the shift of the edge at one third of maximum $\Delta_{1/3}$ obeys relation $\Delta_{1/3}\approx\Delta_{\sqrt{M_{2}}}$. DMC data with IOM estimates are in good agreement for strong EPI $\lambda=0.75$ (Fig. 13). However, shift of the FCP top $\Delta_{p}$ and $\Delta_{1/2}$ are considerably enhanced in the self-trapping (ST) transition region. The physical reason for enhancement of IE in this region is a general property regardless of the QP dispersion, range of EPI, etc. The influence of nonadiabatic matrix element, mixing excited and ground states, on the energies of resonances essentially depends on the phonon frequency. While in the adiabatic approximation ST transition is sudden and nonanalytic in $\lambda$ R82 , nonadiabatic matrix elements turn it to smooth crossover Gerlach . Thus, as illustrated in Fig. 13a, the smaller the frequency the sharper the kink in the dependence of excited state energy on the interaction constant

In the undoped case the present results can be directly compared with the experiments. It is found that the IE on the ARPES lineshape of a single hole is anomalously enhanced in the intermediate coupling regime while can be described by the simple independent oscillators model in the strong coupling regime. The shift of FCP top and change of the FCP height are relevant quantities to pursue experimentally in the intermediate coupling regime since IE on these characteristics is enhanced near the self trapping point. In contrast, shift of the leading edge of the broad peak is the relevant quantity in the strong coupling regime since this value increases with coupling as $\sqrt{\lambda}$. These conclusions, depending on the fact whether self trapping phenomenon is encountered in specific case, can be applied fully or partially to another compounds with strong EPI Dessau_98 ; Perf_01 ; Perf_02 ; Schrupp_05 .

6.3 Conclusions and Perspectives

In this article, we have focused mainly on the polaron problem in strongly correlated systems. This offers an approach from the limit of low carrier concentration doped into the (Mott) insulator, which is complementary to the conventional Eliashberg-Migdal approach for the EPI in metals. In the latter case, we have the Fermi energy $\varepsilon_{F}$ as a relevant energy scale, which is usually much larger than the phonon frequency $\Omega_{0}$. In this case, the adiabatic Migdal approximation is valid and the vertex corrections, which correspond to the multi-phonon cloud and are essential to the self-trapping phenomenon, are suppressed by the ratio $\Omega_{0}/\varepsilon_{F}$. Therefore an important issue is the crossover from the strong coupling polaronic picture to the weak coupling Eliashberg-Migdal picture. This occurs as one increases the carrier doping into the insulator. As is observed by ARPES experiments in high temperature superconductors, the polaronic states continue to survive even at finite doping KyleDop . This suggests a novel polaronic metallic state in underdoped cuprates, which is common also in CMR manganites Manella_05 and is most probably universal in transition metal oxides. In the optimal and overdoped region, the Eliashberg-Migdal picture becomes appropriate Cuk2004 ; Dev2004 , but still a nontrivial feature of the EPI is its strong momentum dependence leading to the dichotomy between the nodal and anti-nodal regions. It is an interesting observation that the highest superconducting transition temperature is attained at the crossover region between the two pictures above, which suggests that both the itinerancy and strong coupling to the phonons are essential to the quantum coherence. It should be noted that this crossover occurs in a nontrivial way also in the momentum space, i.e., the nodal and anti-nodal regions behave quite differently as discussed in Sect. 6.2. However, the relevance of the EPI to the high $T_{c}$ superconductivity is still left for future investigations.

We hope that this article convinces the readers the vital role of ARPES experiments and numerically exact solutions to the EPI problem, the combination of them offers a powerful tool for the momentum-energy resolved analysis of these rather complicated strongly correlated electronic systems. This will pave a new path to the deeper understanding of the many-body electronic systems.

We thank Y. Toyozawa, Z. X. Shen, T. Cuk, T. Devereaux, J. Zaanen, S. Ishihara, A. Sakamoto, N. V. Prokof ev, B. V. Svistunov, E. A. Burovski, J. T. Devreese, G. de Filippis, V. Cataudella, P. E. Kornilovitch, O. Gunnarsson, N. M. Plakida, and K. A. Kikoin, for collaborations and discussions.