Polaron mobility in the"beyond quasiparticles"regime

In a number of physical situations, from polarons to Dirac liquids and to non-Fermi liquids, one encounters the"beyond quasiparticles"regime, in which the inelastic scattering rate exceeds the thermal energy of quasiparticles. Transport in this regime cannot be described by the kinetic equation. We employ the Diagrammatic Monte Carlo method to study the mobility of a Fr\"{o}hlich polaron in this regime and discover a number of non-perturbative effects: a strong violation of the Mott-Ioffe-Regel criterion at intermediate and strong couplings, a mobility minimum at $T \sim \Omega$ in the strong-coupling limit ($\Omega$ is the optical mode frequency), a substantial delay in the onset of an exponential dependence of the mobility for $T<\Omega$ at intermediate coupling, and complete smearing of the Drude peak at strong coupling. These effects should be taken into account when interpreting mobility data in materials with strong electron-phonon coupling.

Mobility of non-degenerate charge carriers in ionic semiconductors with strong electron-phonon coupling is a long-standing problem [1][2][3]. Textbook treatment starts with the notion of a transport scattering time τ controlling momentum relaxation in the equation of motioṅ p = E − p/τ , where E is the external electric field (we set e = 1,h = 1, and k B = 1, for brevity). The underlying assumption is that all effects originating from coupling the particle to its environment (in this work we consider a single particle coupled to a translationallyinvariant bath) are reduced to an effective friction force. The equation of motion leads to a familiar result for the frequency-dependent Drude mobility, with M being the bare particle mass. Deducing the scattering time from the DC mobility by using τ = µ(0)M ≡ µM faces a problem because coupling to the environment necessarily leads to mass renormalization, M → M * . One might argue that substituting M * instead of M into Eq. (1) would solve the issue. However, this procedure implicitly assumes that M * can be measured or calculated separately in a setup where particles propagate coherently and relaxation processes can be neglected at the relevant energy scales. In other words, one has to require that Eτ (E, T ) 1 (in general, τ may depend on the particle energy, E), or if E is replaced with the typical thermal energy (d/2)T , with d being the spatial dimension. Here, τ (T ) should be understood as the inelastic scattering time. For elastic scattering, either by impurities or phonons above the Bloch-Grüneisen temperature, condition (2) is not relevant [4]. Likewise, it can be strongly violated in lattice models in the hopping mobility regime [5,6] when the local rather than momentum representation is more adequate. Condition (2) can be reformulated as λ T , where is the mean free path and λ T is the particle de Broglie wavelength. In this form, it coincides with the "thermal" version of the Mott-Ioffe-Regel (MIR) criterion for the validity of the kinetic equation approach, which has attracted a lot of attention recently in the context of bad and strange metals [7][8][9].
In this Letter, we analyze the violation of the MIR criterion for a tractable yet non-trivial system, namely the Fröhlich polaron model [10]. Apart from being a canonical polaron problem [11,12], it plays the main role in understanding charge transport in ionic semiconductors. It is described by the Hamiltonian H = H e +H ph +H e−ph where with V (q) = i 2 3/2 πΩ 3/2 α/M 1/2 q 2 . Coupling between electrons and optical phonons with energy Ω can no longer be assumed weak when the dimensionless constant α > 1. Perturbation theory predicts that the polaron Zfactor in the ground state is given by Z = 1 − α/2 to leading order. As revealed by stochastically exact Diagrammatic Monte Carlo (diagMC) simulations [13,14], the actual Z-factor is about 0.3 already at α = 2 and drops down to ≈ 0.01 at α = 6. As far as the Fröhlich polaron mobility at finite temperature is concerned, rigorous analytical results can only be obtained within some version of perturbation theory, based either on weak coupling or the Migdal theorem. In the anti-adiabatic regime (T Ω), one can only rely on weak coupling (α 1). This limit was considered by Kadanoff and Langreth [2], who calculated all diagrams for µ up to order O(1/α) + O(1) with the result The O(1) correction accounts for effective mass renormalization, M * = M/(1 − α/6), and is legitimate. There is also a number of variational results; a widely used one is the Low and Pines formula [1], which has the same form as (5) with M * = M (1 + α/6) 3 /f (α), where f (α) ≈ 5/4 for 3 < α < 6. In the adiabatic limit (T Ω), one can neglect vertex corrections when computing the scattering rate and find the mobility from the kinetic equation, whose solution is radically simplified by the fact that scattering is quasi-elastic. The result is a mobility that slowly decreases with T [15], The 1/ √ T scaling of µ is a combination of the T scaling of the phonon occupation numbers and the √ T scaling of the thermal velocity. In contrast to Eq. (5), α is not is not required to be small for Eq. (6) to be valid.
For α 1 both limits predict that at T ∼ Ω one enters the "beyond quasiparticles" regime, in which τ (T ) -defined as M µ(T )-is shorter than the Planckian bound. As far as we know, the mobility of Fröhlich polarons has never been computed with controlled accuracy in this parameter regime.
In this Letter, we address this unsolved problem by extending the diagMC method [13,14] to finite temperatures in combination with analytic continuation methods [14,16,17] and computing the mobility from the Matsubara current-current correlation function. We find that for α > 1 the MIR condition, formulated as µ 1/M T , is indeed violated. Contrary to expectations based on perturbation theory, the exponential increase of the mobility with decreasing temperature, Eq. (5), is observed only at temperatures significantly below Ω. Even more surprising is a mobility minimum developing at T ∼ Ω for large values of α. We interpret both effects as a competition between the decreasing number of thermal excitations and a strong enhancement of the polaron mass, both of which taking place as temperature is lowered. Method. As far as the diagMC method is concerned, our finite temperature calculation (in which we used the units M = Ω = 1) is similar to that employed for calculating the optical conductivity in Ref. [18] for T = β −1 → 0. The difference is that now phonon propagators in imaginary time are temperature dependent, D(τ ) = [e −Ωτ + e Ωτ −Ωβ ]/[1 − e −Ωβ ], and we do not assume that e −Ωβ is small. A typical diagram contributing to the current-current correlation function A typical Feynman diagram contributing to the current-current correlation function. Solid (dashed) lines stand for electron (phonon) propagators G (D), dots are representing interaction vertices |V |, and crosses mark velocity measurements.
By using Matsubara frequencies, ω m = 2πmT (with m integer), we obtain an efficient unbiased Monte Carlo estimator for its Fourier transform, The value of k m for a diagram of order N is readily found from the electron momenta {k i } in imaginary time intervals i = 1, . . . , 2N ranging from τ i to τ i+1 (cf. Fig 1), There are several exact and asymptotic relations relations that the correlation function has to satisfy. Any of them can be used as an independent check on the accuracy of the data. For the equal-time correlator we find that T m C m = k 2 /dM 2 . The sum rule for the mobility translates into C m=0 = 1. The asymptotic behavior in the limit of large |m| comes from diagrams containing very short time intervals τ i+1 − τ i = ∆τ → 0 with large particle momenta k i ∼ 1/ √ ∆τ → ∞. Dressing such intervals by vertex corrections is not necessary because it would result in additional factors of √ ∆τ . Therefore, one can use perturbation theory to compute the highfrequency limit with the final result To determine the mobility, one needs to perform the analytic continuation numerically, which amounts to inverting the Kramers-Kronig transformation  (5) and (6), respectively. Below the dashed red line, µ = 1/M T , the MIR criterion is violated.
The asymptotic behavior implied by C m immediately yields the asymptotic behavior of µ(ω) which can be used to subtract the leading tail contribution when performing the analytic continuation. The nature of Eq. (10) is such that even tiny error bars on C m translate into large uncertainties on integrals of µ(ω) over physically relevant intervals. All errors on the µ(ω) function itself are conditional and depend on constraints for allowed behavior [16]. In this work, Monte Carlo data for the lowest Matsubara frequencies were accurate at the level of five to six significant digits. Our analytic continuation method for solving Eq. (10) is similar to the one used in Ref. [5]. The main difference is that we extract µ(ω) from the data parameterized by the Matsubara frequency rather than imaginary time. We also employ more conservative protocols for computing both the mobility and its errorbars from multiple solutions of the stochastic optimization with the consistent constraints method (see Supplemental material) [14,16].
Results. To begin with, we verified that our results for the mobility are consistent with the predictions (5) and (6) based on perturbation theory for small α = 0.25. In this case, the MIR criterion is not violated even at T ∼ Ω, i.e., the mobility remains large compared to 1/M T . When the coupling strength is increased to α = 2.5, the situation changes. As Fig. 2 shows, the MIR criterion is violated over a broad temperature interval 0.2 < T /Ω < 10. The data for T > Ω approach the limiting form in Eq. (6), (color online). Mobility of a Fröhlich polaron as a function of temperature at strong coupling (α = 6). All notations are identical to those in Fig. 2. Eq. (5) is plotted using the exact result M * = 7.3 instead of its perturbative approximation; otherwise the curve would go unphysical. The Low-Pines formula [1] is shown by a dotted line [19].
as expected, because Eq. (6) is valid for any α. On the other hand, one should not expect the low-temperature formula (5) to be valid beyond weak coupling. Indeed, it is clear from Fig. 2 that the slow temperature dependence extends down to at least T ∼ Ω/2. At lower temperatures, our data are consistent with the exponential increase of the mobility, but uncertainties amplified by the analytic continuation procedure become too large for a meaningful analysis.
A delay (at T < Ω) in the onset of the exponential dependence may be anticipated already from perturbation theory. Indeed, matching Eqs. (5) and (6), we find that the crossover temperature with λ ≡ M * /M is (logarithmically) suppressed compared to Ω for λ 1. One might argue that the MIR criterion cannot be specified down to a well-defined numerical factor and, thus, having µM T = 1/4 for α = 2.5 is merely a borderline situation. However, the strong coupling case (α = 6) (cf. Fig. 3) shows that there is no obvious limit on how small the value of µM T can be. Even more surprising is the finding that the mobility develops a minimum at T /Ω < ∼ 1 and that the high-temperature limit is recovered only after a maximum at T /Ω ∼ α, see Fig. 3. The µM T value at the minimum is about 1/30.
The mobility minimum for lattice polarons is a wellestablished phenomenon (see, e.g., Refs. [5,6,[20][21][22]), known as "activated hopping" between lattice sites. However, we are not aware of any work on the mobility minimum for Fröhlich polarons because the notion of hopping is ill-defined in the continuum. To understand the origin of the non-monotonic behavior, we note that for α = 6 the effective mass renormalization at low temperatures and small momenta is very strong, M * /M ≈ 7.3 [14], while the physics at high temperatures is still captured by the Migdal theorem. The mobility minimum then emerges from the competition between the decreasing number of thermal excitations and the strong particle mass renormalization: The former dominates at T /Ω 1, while the latter is more important at T /Ω > ∼ 1. If renormalization of the effective mass, which is temperature and energy dependent, does play such an important role, then the dependence of µ on ω should be radically different from the prediction of the Drude model in Eq. (1). To check this assertion we present the real part of µ(ω, T ) in Fig. 4. At the two lowest temperatures (T = 0.125 Ω and T = 0.25 Ω), µ(ω, T ) exhibits a narrow Drude peak at low frequencies and a broad maximum at intermediate frequencies.
In the context of small lattice polarons, the broad maximum in the mobility is commonly attributed to ionization of bound states [23]. In a broader context, however, this maximum is referred to as "Holstein band" and arises because emission of phonons is activated for ω above the characteristic phonon frequency [24,25]. The Fröhlich model is just one example; another one is the Holstein model the result for which is shown in the inset in Fig. 4.
For T /Ω > ∼ 0.5 the Drude peak in the mobility disappears. Moreover, for T /Ω > ∼ 2.0 the broad maximum becomes barely distinguishable. Any attempt to extract the scattering time from the frequency dependence of the mobility obviously fails in this regime. Our interpretation instead is to use the physically appealing relation τ −1 ∼ E ∼ T in order to conclude that the effective mass defined through M * ∼ τ /µ undergoes strong renormalization as the temperature is lowered down to T < Ω (which is further confirmed by calculating the effective mass at T = 0).
Even if the mobility retains a Drude peak but the MIR criterion is violated, the scattering time can not be extracted from the Drude formula. This is supported by an example of a heavy particle with a large transport scattering cross-section embedded into a Fermi sea, e.g., an electron bubble in 3 He. Its mobility is described by Eq. (1), cast in the form µ(ω) = 1/M (Γ − iω), with a width Γ ≈ const down to temperatures exponentially lower than Γ [26]. This result is also well known in the context of Ohmic dissipative models [27,28]. Interpreting Γ as the scattering rate at temperatures T Γ implies that the particle mass is not renormalized. However, if the Planckian bound on the scattering rate 1/τ (T ) is to be obeyed, one is forced to conclude that the effective mass has to diverge as M * ∝ M Γ/T in order to maintain a constant DC mobility µ = 1/M Γ ∼ τ (T )/M * (T ). That the last interpretation is physically correct is revealed by considering the superfluid state at T T C Γ, when the undamped particle motion is controlled by the strongly renormalized effective mass M * ∝ M Γ/T C [26].
In conclusion, we have studied the problem of electron mobility in ionic semiconductors in the temperature regime where the MIR criterion for the applicability of the kinetic-equation approach is violated. We found that the mobility can be orders of magnitude smaller than the MIR value of 1/M T . This result is consistent with recently observed MIR violation for non-degenerate charge carriers in doped SrTiO 3 [29]. At strong coupling, the mobility has a minimum at a temperature comparable to the optical mode frequency and increases with a further increase of the temperature despite the fact that the number of thermally excited phonons grows linearly with T . We ascribe this behavior to the "undressing" of polarons at higher temperatures. After going through a maximum at T Ω, the mobility follows the µ ∼ T −1/2 law predicted by the kinetic equation. cussions with K. Behnia and T. Timusk. This work was initiated and partially performed at Aspen Center for Physics supported by the National Science Foundation grant PHY-1607611. We acknowledge support by the Simons Collaboration on the Many Electron Problem (N.P.), the National Science Foundation under grant DMR-1720465 (N.V.P.) and DMR-1720816 (A.K. and D.L.M.), the ImPACT Program of the Council for Sci-ence, Technology and Innovation (Cabinet Office, Government of Japan) and JST CREST Grant Number JP-MJCR1874, Japan (A.S.M.), and by FP7/ERC Consolidator grant QSIMCORR No. 771891, the Nanosystems Initiative Munich and the Munich Center for Quantum Science and Technology (L.P.). Open Data for this project can be found at https://gitlab.lrz.de/ QSIMCORR/mobilityfroehlich where τ depends both on ω and T . After passing a Drude peak of width ∼ 1/τ (T, 0), the mobility enters the "collisionless" regime in which For any kind of electron-phonon interaction, 1/τ (ω, T ) increases with ω faster then ω 2 , and thus µ (ω, T ) increases with ω. This increase is curbed off at ω larger than the characteristic phonon frequency, when 1/τ (ω, T ) becomes independent of ω and µ (ω, T ) drops off with ω. This is the origin of the maximum in µ (ω, T ). We illustrate this behavior on the example of Holstein model, 1 which describes a short-range electron-phonon interaction with optical phonons. To lowest order in the electron-phonon coupling, the imaginary part of self-energy is given by where g is the electron-phonon coupling constant, is the retarded Green's function of free electron gas with chemical potential µ, is the Green's function of optical phonons with frequency Ω, n F (ε) = e ε/T + 1 −1 is the Fermi function with µ = 0, and N 0 (Ω ) is the Bose function. For non-degenerate electrons, it is convenient to introduce a new energy variable E ≡ ε + µ, in terms of which Here, Σ 0 = πg 2 ν(Ω)/2, f 0 (E) is the Fermi function in the Maxwell-Boltzmann limit given by Eq. (7), ν(E) is the electron density of states, and θ(x) is the Heaviside step-function. For a short-range electron-phonon interaction, the vertex corrections to the current-current correlation function are absent, and the Kubo formula for the real part of the mobility to leading order in g is reduced to where n is the number density of electrons, M is the electron mass, and ω is chosen to be positive without loss of generality. At weak coupling, |Σ (E)| E. The spectral functions in Eq. (15) are sharply peaked only for E > 0, which effectively limits the range of integration over E to interval (0, ∞). Solving the integral over p with condition |Σ (E)| E taken into account and using Eq. (7) for f 0 (E), we obtain The remaining integral over E is solved numerically. The resultant frequency dependence µ is shown in the inset of Fig. 4 of the main text for T /Ω = 0.125 and Σ 0 /Ω = 0.5. For non-degenerate electrons (T µ ≈ E F , where E F is the Fermi energy), one only has to replace E by E F in the square-root factors and by ε in the Fermi functions in Eq. (14). Assuming also the adiabatic limit (Ω E F ), we then obtain for the imaginary part of the self-energỹ whereΣ 0 = πg 2 ν(E F )/2. In the Kubo formula, one replaces the integration measure d 3 p/(2π) 3 by ν(E F )dξ p dO/4π, where O is the solid angle and ξ p = p 2 /2m − E F , and takes the factor of p/M ≈ v F ≡ 2E F /M outside the integral. After integration over ξ p , we then find The resultant µ(ω) is shown in Fig. 1 for T /Ω = 0.125 andΣ 0 /Ω = 0.1.

C. Analytic continuation
To extract the mobility and its uncertainty from solutions of Eq. (10) in the main text we employed the stochastic optimization with consistent constraints (SOCC) method 2,3 . First, we find N solutions, {µ i (ω)}, i = 1, . . . , N , that would reproduce Matsubara frequency data within their error bars (our accuracy was better than five significant digits for the lowest frequencies). The set of mobilities {µ i } was obtained from the limiting values of frequency-dependent solutions, µ i = µ i (ω → 0), for which we considered several protocols: (i) The value of µ i (ω) at the lowest frequency, (ii) linear, and (iii) quadratic fits+extrapolation of µ i (ω) data to ω = 0 (from the range of small frequencies where µ i (ω) is structureless). However, differences in final results between these protocols were found to be much smaller than the statistical error bars described below.
Despite having highly accurate data in the Matsubara frequency domain the distribution of possible mobility values turns out to be rather broad, see Fig. 2, which is a typical scenario for such problems. Given a typical probability distribution as shown in Fig. 2  its error bars: (i) An unrestricted average over all mobility values, µ =μ = (1/N ) N i=1 µ i ; (ii) An average over a restricted set of N * solutions with µ i = 0 (some solutions have zero spectral density at ω = 0 and they are accounted for by the SOOC method), µ =μ (r) = (1/N * ) N i=1 µ i Θ(µ i + 0 + ); (iii) The most probable value; i.e., the location of the maximum of the probability distribution function shown in Fig. 2. We denote it as µ = µ (max) .
Since the distribution of mobilities is asymmetric (limited from below and unlimited from above) we estimate both the upper and lower error bars. The upper(lower) error bars are then defined though asymmetric dispersions computed from statistics of N > (N < ) solutions with µ i > µ (µ i < µ): We find that in all cases considered the error bars σ (u) and σ (d) are much larger than differences between the alternative definitions of µ. For example, for α = 6 and T /Ω = 8, see Fig. 2, the distribution is such thatμ = 0.037,μ (r) = 0.040, µ (max) = 0.036, and σ (d) = 0.016, σ (u) = 0.018. All data presented in the main text are based on averages over restricted sets of N * solutions; i.e. µ =μ (r) and