# Ground-state dispersion and density of states from path-integral Monte Carlo. Application to the lattice polaron

## I Introduction

Usually, quantum Monte Carlo (QMC) methods are used to study ground-state, thermodynamic, or static properties of quantum-mechanical systems. Some important dynamical characteristics can also be obtained through various forms of fluctuation-dissipation relations. Examples of these are the superfluid fraction of Bose-liquids , the Drude weight of conductors , the Meissner fraction of superconductors , and the effective mass of defects  and polarons . Beyond that dynamical calculations are less straightforward. For instance, calculation of the excitation spectrum normally requires measurement of the Green’s function at imaginary times and subsequent analytical continuation to real times.

However, there exists one special type of excitation spectrum that can be measured directly by QMC. This is the ground-state dispersion, i.e., the total energy of the system $`E_𝐏`$ as a function of the total momentum P. In a translationally-invariant system, P is a constant of motion, and the Hamiltonian does not mix subspaces with different P. Then, if a QMC is designed as to operate within a given P-subspace only, it may be able to access the ground state for the given P, thereby providing $`E_𝐏`$. Not for any physical system $`E_𝐏`$ is of interest. For a collection of identical particles, for instance, one has simply $`E_𝐏=𝐏^2/(2M)`$, $`M`$ being the total mass, which corresponds to free movement of the system as a whole. Positive examples include cases when the system can be divided into a tagged particle and an environment (usually bosons). The best known example of this kind is the polaron, i.e., an electron strongly interacting with phonons. In this case $`E_𝐏`$ is nothing but the polaron spectrum. The polaron spectrum will be the main subject of this paper.

There exist at least two different strategies of how to operate within a restricted P-subspace. The first one is to work in momentum space and to fix the total momentum of the system from outset. An example of this approach is the diagrammatic method of Prokof’ev and Svistunov . In this method QMC is used to sum the entire diagrammatic series for an imaginary-time Green’s function $`G(𝐏,\tau )`$. Since the total momentum is an external parameter of the series, it is possible to extract $`E_𝐏`$ from the $`\tau \mathrm{}`$ limit behavior of $`G`$, by fitting it to a single exponential $`e^{E_𝐏\tau }`$. This method is exact and universal but requires a separate simulation for each P-point.

The second strategy is to work in real space but to use Fourier-type projection operators to project on states with definite P. This amounts to the free boundary conditions in imaginary-time. Usually, the projections are used to calculate the second derivative of the energy with respect to momentum (effective mass) . In Ref.  the projection was applied for the first time to the whole polaron spectrum. In this scheme, the ground-state dispersion is measured directly, and all $`E_𝐏`$ are calculated simultaneously. Unfortunately, at non-zero P the weight of the polaron path is no longer positive definite, and one needs to deal with a sign-problem. It turns out, however, that the main idea can be reformulated in a way that does not require a division by the average sign, but only taking its logarithm. While the new formulation does not constitute a complete elimination of the sign-problem, it is more statistically stable and extends the parameter domain accessible in practical simulations. Below we derive the new formula, discuss its properties, and apply to the physically interesting example of the lattice polaron.

## II A formula for the ground-state dispersion

Up to our knowledge, the projection relations required for our purposes were first derived by Basile  . For completeness a derivation is given below. Let $`R`$ denote a many-body real-space configuration, and $`R+𝐫`$ a many-body configuration which is a result of the parallel transport of $`R`$ by a vector $`𝐫`$. (Note that the sum of $`R`$ and r is only symbolic. The dimensionality of $`R`$ is equal to the number of degrees of freedom, i.e., very large or infinite, while the dimensionality of r is the dimensionality of space, i.e. 1, 2, or 3.) States $`|R`$ form a complete orthogonal basis, $`𝐈=𝑑R|RR|`$, and $`R|R^{}=\delta (RR^{})`$. A different basis is formed by the states $`|n`$ which are characterized by the definite total momentum $`𝐏`$. One is interested in the projected partition function $`Z_𝐏`$ which includes only states with the given P:

$`Z_𝐏`$ $``$ $`{\displaystyle \underset{n}{}}n|e^{\beta H}|n\delta _{𝐏,𝐏_n}`$ (1)   
$`=`$ $`{\displaystyle 𝑑R𝑑R^{}R^{}|e^{\beta H}|RQ_𝐏},`$ (2)   
$$Q_𝐏=\underset{n}{}R|nn|R^{}\delta _{𝐏,𝐏_n}=R|𝐏𝐏|R^{},$$
(3)   
where $`\beta =(k_BT)^1`$ is the inverse temperature and $`H`$ is the full Hamiltonian. The meaning of Eq. (3) is that the two configurations, $`R`$ and $`R^{}`$, have to be projected on the given momentum $`𝐏`$. This is achieved as follows (below $`\mathrm{}=1`$ is set). Any arbitrary configuration $`R`$ generates a set of states $`|𝐏_R=V^{1/2}𝑑𝐫e^{i\mathrm{𝐏𝐫}}|R+𝐫`$, where $`V`$ is the volume. Inversely, $`|R+𝐫=V^{1/2}_𝐏e^{i\mathrm{𝐏𝐫}}|𝐏_R`$. Upon projection, only P-components of both configurations survive. As a result

$`Q_𝐏`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle 𝑑𝐫R+𝐫|𝐏𝐏|R^{}+𝐫}`$ (4)   
$`=`$ $`{\displaystyle \frac{1}{V^2}}{\displaystyle 𝑑𝐫\underset{𝐏^{}𝐏^{\prime \prime }}{}e^{i(𝐏^{}𝐏^{\prime \prime })𝐫}𝐏_{}^{}{}_{R}{}^{}|𝐏𝐏|𝐏_{}^{\prime \prime }{}_{R^{}}{}^{}}`$ (5)   
$`=`$ $`{\displaystyle \frac{1}{V}}𝐏_R|𝐏_R^{}`$ (6)   
$`=`$ $`{\displaystyle \frac{1}{V^2}}{\displaystyle 𝑑𝐫𝑑𝐫^{}R+𝐫|R^{}+𝐫^{}e^{i𝐏(𝐫𝐫^{})}}`$ (7)   
$`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle d(\mathrm{}𝐫)R+\mathrm{}𝐫|R^{}e^{i𝐏\mathrm{}𝐫}}`$ (8)   
$`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle d(\mathrm{}𝐫)e^{i𝐏\mathrm{}𝐫}\delta ((R+\mathrm{}𝐫)R^{})},`$ (9)   
where $`\mathrm{}𝐫=𝐫𝐫^{}`$. Substitution in Eq. (2) and integration over $`R^{}`$ yields

$`Z_𝐏`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle d(\mathrm{}𝐫)e^{i𝐏\mathrm{}𝐫}𝑑RR+\mathrm{}𝐫|e^{\beta H}|R}`$ (10)   
$`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle d(\mathrm{}𝐫)e^{i𝐏\mathrm{}𝐫}𝑑R\rho (R,R+\mathrm{}𝐫;\beta )},`$ (11)   
where $`\rho (R,R^{};\beta )`$ is the full many-body density matrix. Next, we assume that for each P the state with the lowest energy $`E_𝐏`$ is non-degenerate, and in the low-temperature limit the projected partition function is dominated by the contribution from this state, $`Z_𝐏\mathrm{exp}(\beta E_𝐏)`$. Now take the ratio of $`Z_𝐏`$ and $`Z_{𝐏=0}`$:

$`e^{\beta (E_𝐏E_0)}`$ $`=`$ $`\underset{\beta \mathrm{}}{lim}{\displaystyle \frac{Z_𝐏}{Z_{𝐏=0}}}`$ (12)   
$`=`$ $`\underset{\beta \mathrm{}}{lim}{\displaystyle \frac{d(\mathrm{}𝐫)e^{i𝐏\mathrm{}𝐫}𝑑R\rho (R,R+\mathrm{}𝐫;\beta )}{d(\mathrm{}𝐫)𝑑R\rho (R,R+\mathrm{}𝐫;\beta )}},`$ (13)   
where $`E_0`$ is the ground-state energy. The rhs is nothing but the average value of $`\mathrm{cos}𝐏\mathrm{}𝐫`$ taken over the distribution $`\rho `$. \[We have assumed that $`𝑑R\rho (R,R+\mathrm{}𝐫;\beta )`$ is an even function of $`\mathrm{}𝐫`$.\] A simple formula for $`E_𝐏`$ now follows

$$E_𝐏E_0=\underset{\beta \mathrm{}}{lim}\frac{1}{\beta }\mathrm{ln}\mathrm{cos}𝐏\mathrm{}𝐫,$$
(14)   
which is the main result of this section.

Eq. (14) shows that the ground-state dispersion can be obtained from the end-to-end distribution of many-body paths. It offers a direct way of evaluating the ground-state dispersion by QMC methods in cases when $`\rho (R,R^{};\beta )`$ is positive-definite. However, the QMC process must be organized in a special way, as apparent from Eq. (12). It must generate only such paths whose end configurations at imaginary time $`\tau =\beta `$ are exact images of the end configurations at $`\tau =0`$ except for a parallel transport by an arbitrary vector $`\mathrm{}𝐫`$. It is allowed to change $`\mathrm{}𝐫`$, make simultaneous changes of both end configurations, and to make arbitrary changes of paths at internal times $`0<\tau <\beta `$, but the end configurations must always be kept identical up to a shift. It is important that this restriction affects neither the ergodicity nor the applicability of the Metropolis algorithm.

Formula (14) involves only one measures quantity, $`\mathrm{cos}𝐏\mathrm{}𝐫`$, instead of two in the previous formulation. Moreover, it does not require division by the measured quantity, but only taking its logarithm. Additionally, Eq. (14) provides the difference between two large numbers, $`E_𝐏`$ and $`E_0`$, and a large cancellation of errors may occur. This makes Eq. (14) much more stable statistically than explicit energy estimators.

On the other hand, at small temperatures, the average cosine becomes exponentially small and it cannot be measured reliably. This reflects the fact that configurations with $`E_𝐏E_0k_BT`$ are very rare because of the Boltzmann factor. Thus, the present method is limited to excitation energies of the order of several $`k_BT`$.

## III Application to the lattice polaron

We now demonstrate the practical importance of Eq. (14) on the model problem of lattice polaron, which is often considered as a paradigmatic example of a particle strongly interacting with a boson field. We consider a hypercubic lattice with the nearest-neighbor hopping, dispersionless phonons, and the “density-displacement” electron-phonon interaction. The model Hamiltonian reads

$$H=t\underset{\mathrm{𝐧𝐧}^{}}{}c_𝐧^{}c_𝐧^{}\underset{\mathrm{𝐧𝐦}}{}f_𝐦(𝐧)c_𝐧^{}c_𝐧\xi _𝐦+\mathrm{}\omega \underset{𝐦}{}b_𝐦^{}b_𝐦.$$
(15)   
Here $`t`$ is the hopping amplitude (it will be used as the energy unit), $`\omega `$ is the phonon (oscillator) frequency, $`\xi _𝐦`$ is the internal coordinate of the mth oscillator, and $`f_𝐦(𝐧)`$ is the force between mth oscillator and the particle at site n ($`f`$ is a function of distance $`|𝐦𝐧|`$ only). The model is parametrized by the dimensionless frequency $`\overline{\omega }=\mathrm{}\omega /t`$ and by the dimensionless coupling constant $`\lambda =[_𝐦f_𝐦^2(0)]/(2M\omega ^2D)`$, where $`M`$ is the mass of the oscillator and $`D`$ is the half-bandwidth of the bare band. (For an isotropic band with nearest-neighbor hopping, $`D=zt`$, $`z`$ being the number of neighbors.)

For the polaron problem, a many-body configuration $`R`$ is specified by the position of the electron r and oscillator displacements $`\xi 𝐦`$. Making use of the Feynman’s idea of analytic integration over $`\xi 𝐦`$  the problem is reduced to a single-particle system with retarded self-interaction. The latter can be simulated exactly, using the continuous-time representation of polaron paths . The resulting algorithm  is very efficient and allow accurate determination of the ground-state energy and effective mass of the polaron for a wide class of models. It this paper, it will be shown that the method also produces accurate polaron spectra, when combined with Eq. (14). There are other reasons why the polaron is an ideal system to try formula (14). First, due to a constant phonon frequency, excited states are, at any P, separated from the restricted ground state by a finite energy gap $`\mathrm{}\omega `$. Therefore, instead of performing numerically the limit procedure to $`\beta =\mathrm{}`$, one can study the system at finite $`\beta `$, provided $`\mathrm{exp}(\beta \mathrm{}\omega )1`$ and the contribution from excited states is negligible. Second, by increasing the coupling constant $`\lambda `$ one can always decrease $`E_𝐏E_0`$, i.e., substantially increase $`\mathrm{cos}𝐏\mathrm{}𝐫`$, and stabilise the sumulations. Third, the polaron momentum P is not a parameter of simulations. This implies that statistics can be collected for all momenta simultaneously. In other words, the whole polaron spectrum is measured in a single QMC run. This will enable us to calculate for the first time exact polaron densities of states.

We begin with the simplest Holstein model which has local electron-phonon interaction, $`f_𝐦(𝐧)=\kappa \delta _{\mathrm{𝐦𝐧}}`$. In one dimension, the polaron spectrum has been extensively studied by exact diagonalisation , strong-coupling perturbation , and variational  techniques. Our QMC data for the one-dimensional Holstein model are shown in Fig. 1.

The most interesting feature of the spectrum is its non-cosine shape in the adiabatic regime $`\overline{\omega }1.0`$ (triangles and diamonds). At large momenta, the spectrum is more flat than at small momenta. The nature of this flattening was understood a long time ago . In the weak-coupling limit, the free-particle state hybridizes with the single-phonon state and creates a mixed ground state which is free-particle-like at small P and phonon-like at large P, hence the weak dispersion. With increasing $`\lambda `$, the free-particle state is replaced with a polaron state with an effective mass $`m^{}`$, which now hybridizes with the single-phonon state, still leading to a more flat dispersion at large momenta. The flattening effect weakens with growing $`\lambda `$ and $`\overline{\omega }`$ because both processes increase the energy separation of the two hybridizing states. In recent years the flattening of the polaron spectrum was observed in numerical studies .

Our Quantum Monte Carlo data fully confirm the previous analytical and numerical results, see Fig. 1. We found that the spectrum shape is more sensitive to the phonon frequency than to the coupling constant. At a small frequency $`\overline{\omega }=1.0`$, the increase of the coupling constant from $`\lambda =2.0`$ to $`\lambda =2.5`$ results in a 3.5-times increase of the effective mass, and in a 2.8-times drop of the bandwidth, yet the spectrum shape changes only slightly, cf. triangles and diamonds in Fig. 1. At the same time, a simultaneous 10-times increase of $`\lambda `$ and $`\overline{\omega }`$ results in a similar, 4.8 times, increase of effective mass but brings the spectrum shape very close to cosine, cf. triangles and squares in Fig. 1. Again, a doubling of $`\lambda `$ strongly affects $`E_0`$ and $`W`$ but not the spectrum shape, cf. circles and squares in Fig. 1.

It is instructive to compare the exact QMC results with the Lang-Firsov (LF) approximation  which is believed to be the correct description of the polaron in the antiadiabatic regime (high phonon frequency). The LF formula for the spectrum reads

$$E_PE_0=2te^{z\lambda /\overline{\omega }}(1\mathrm{cos}P),$$
(16)   
which also implies the relation between the renormalized mass and bandwidth

$$\frac{m^{}}{m_0}=\frac{W_0}{W}=e^{z\lambda /\overline{\omega }},$$
(17)   
where $`W_0=2zt`$ is the bare bandwidth and $`m_0=\mathrm{}^2/(2ta^2)`$ is the bare mass, $`a`$ being the lattice constant. For $`\overline{\omega }=10.0`$ and $`\lambda =10.0`$ QMC results are $`W=0.543(2)t`$ and $`m^{}=6.06(2)m_0`$ while LF yields $`W_{LF}=0.541t`$ and $`m_{LF}^{}=7.39m_0`$. For $`\overline{\omega }=10.0`$ and $`\lambda =20.0`$ one has $`W=0.0739(2)t`$ and $`m^{}=47.6(1)m_0`$ from QMC and $`W_{LF}=0.0733t`$ and $`m_{LF}^{}=54.6m_0`$ from LF. One can see that LF predicts very accurate values of the polaron bandwidth. This fact was established in the previous studies of the Holstein model . On the other hand, LF slightly overestimates the polaron mass. This is due to small deviations from the cosine shape, still present in the true spectrum at these model parameters. Still, the LF masses are reasonably close to the exact ones, and the agreement improves with the further increase of $`\overline{\omega }`$ and $`\lambda `$.

Consider now the two-dimensional Holstein model. The only exact polaron spectra published so far were calculated by Wellein, Fehske, and Loos with the exact diagonalization method . These authors found a flattening of the spectrum in the outer part of the Brillouin zone, even stronger than in the one-dimensional case. We checked that for the model parameters used in , formula (14) yields precisely the same values of $`E_𝐏E_0`$ as the exact diagonalization method. A definite advantage of the present method is that it allows simultaneous calculations at any desired number of P-points, while exact diagonalization studies are limited to a small number of P-points due to the finite size of the clusters. On the other hand, the QMC method is limited to the condition $`W\mathrm{}\omega `$, which prevents us from studying the weak-coupling regime and such an interesting phenomenon as the limit point of the polaron spectrum . (The latter is possible with the diagrammatic QMC .) Figure 2 shows our QMC data for a new set of parameters in the adiabatic regime, $`\overline{\omega }=1.0`$, $`\lambda =1.4`$, where 30 P-points have been used to represent the spectrum. One can see that the dispersion is indeed weak for $`|𝐏|>\pi /2`$, i.e. in the larger part of the Brillouin zone. Note, that at these parameters the polaron bandwidth is reduced by $`8t/0.12t=67`$ times, and the mass enhancement is 8.7, so we are already in the small polaron regime. Yet, the spectrum shape is profoundly non-cosine. With increasing $`\lambda `$, it will be approaching the cosine shape, but this is expected to happen only at such large $`\lambda `$ where the polaron is very heavy and is easily localized.

In two dimensions, the volume of the outer part of the Brillouin zone is larger than that of the inner one. Therefore the linear representation of the spectrum, like in Fig. 2, does not fully convey the changes in the band structure, caused by the flattening of the spectrum. The proper physical quantity which takes into account all the states of the Brillouin zone is the density of states (DOS). The QMC method, coupled with formula (14), provides the unique opportunity to calculate polaron DOS exactly, since it allows a simultaneous measurement of the whole spectrum. In this work, the two-dimensional Brillouin zone was divided in $`200^2`$ points at which the spectrum was measured. In the end, the total of $`\mathrm{40\hspace{0.17em}000}`$ polaron energies were distributed over 50 energy intervals between 0 and $`W`$. The resulting DOS for $`\overline{\omega }=1.0`$, $`\lambda =1.4`$ (the same parameters as in Fig. 2) is shown in Fig. 3. One can see that the effect of the spectrum flattening is indeed quite dramatic. The upper half-band is jammed into a narrow, $`0.015t`$-width, energy interval, thereby increasing DOS at the top of the band to $`50`$ times the DOS at the bottom of the band. The van Hove singularity is shifted from the middle to the top of the band. The lower half of the band contains only $`13\%`$ of all states. Overall, DOS looks qualitatively different from the free-particle one.

As in the one-dimensional case, the band structure approaches the free-particle-like as the phonon frequency increases. As an example, we considered a frequency equal to the bare bandwidth, $`\overline{\omega }=8.0`$, and $`\lambda =8.0`$. (This value of the coupling constant was chosen to have the polaron bandwidth, $`W`$, close to the previously considered adiabatic case.) The polaron spectrum and DOS are shown in Figs. 4 and 5, respectively. Deviations from the free-particle behavior are still visible, but they are small and not qualitative. DOS at the top of the band is just $`1.5`$ times larger than at the bottom, and the singularity appears close to the middle of the band. The “normalization” of the spectrum at large $`\overline{\omega }`$ is quite understandable. With increasing phonon frequency, the retardation effects become less important, the lattice deformation more readily follows the particle movement, and the whole complex behaves more like a free particle with a renormalized hopping integral. The Lang-Firsov formula (17) predicts $`W_{LF}=0.1465t`$ and $`m_{LF}^{}=54.6m_0`$ which is to be compared with the QMC results $`W=0.1510(3)t`$ and $`m^{}=38.4(1)m_0`$. Again, the LF approximation yields the correct bandwidth but overestimates the effective mass by some 40%.

The most spectacular transformation of DOS occurs in the three-dimensional Holstein model. In three dimensions, the volume of the outer part of the Brillouin zone is much larger than that of the inner part, and it should completely dominate the total DOS. We do not show the polaron spectrum which is not very informative. Density of states was calculated by measuring the polaron spectrum at $`60^3`$ points of the full Brillouin zone, and then distributing them among 50 energy intervals between 0 and $`W`$. DOS in the adiabatic regime, $`\overline{\omega }=1.0`$ and $`\lambda =1.2`$, is shown in Fig. 6. The states of the flat part of the spectrum form a massive peak at the top of the band. The width of the peak is about $`10\%`$ of the total bandwidth. DOS at the bottom of the band is negligible, the capacity of the lower half-band is less than $`1\%`$ of the total number of states. The two van Hove singularities are not visible, at least on the chosen level of energy resolution, they are absorbed into the peak. Overall, DOS is completely different from the free-particle one. Should the three-dimensional Holstein model with such parameters exist in nature, an extreme care would be necessary in interpreting experimental data. In any real material, the lowest states would likely be localized, and any response to an external perturbation would be dominated by the peak. Then, for instance, fitting to a free-particle-like form of DOS would lead to wrong estimates of the coupling constant and other errors.

As before, the band structure returns to the free-particle shape in the anti-adiabatic regime. We considered the case of the phonon frequency being equal to the bare bandwidth, $`\overline{\omega }=12.0`$, and coupling constant $`\lambda =10.0`$, when the polaron bandwidth is close to the just considered adiabatic case, see Fig. 7. Although still distorted, the DOS shape is close to the free-particle one, with the square-root behavior at the top and the bottom of the band, and with two van Hove singularities fully developed at the “right” places. The polaron bandwidth is $`W=0.0827(2)t`$, which is in good agreement with the Lang-Firsov value $`W_{LF}=0.0809t`$, while the polaron mass $`m^{}=112(1)m_0`$ is 24% lighter than the LF mass $`m_{LF}^{}=148m_0`$.

The polaron QMC algorithm of Ref.  is not limited to the Holstein model. In fact, it allows studies of arbitrary forms of the electron-phonon interactions (of the density-displacement type), and arbitrary forms of the particle kinetic energy. Combined with Eq. (14), it provides an efficient and exact way of calculating the band structure of a whole class of polaron models. As possibilities are numerous, we have chosen to illustrate the point on two particular examples.

The first example is the anisotropic two-dimensional Holstein model with $`\overline{\omega }=1.0`$, $`\lambda =1.4`$ (these parameters are the same as in Figs. 2 and 3), and the bare anisotropy ratio $`t_y/t_x=0.2`$. For a free particle with such an anisotropy, the saddle points at $`(\pm \pi ,0)`$ and $`(0,\pm \pi )`$ have different energies, which results in two singularities in DOS, positioned symmetrically with respect to the center and edges of the band. Polaron DOS, calculated by QMC, is shown in Fig. 8. The flattening effect creates a strong peak at the top of the band which absorbs the higher-energy singularity \[at $`(\pm \pi ,0)`$\]. At the same time, the second singularity is still clearly visible. Now it appears in the vicinity of the middle of the band.

The second example is the two-dimensional polaron model with long-range electron-phonon interaction \[combine with the Hamiltonian (15)\]:

$$f_𝐦(𝐧)=\frac{\kappa }{(|𝐦𝐧|^2+1)^{3/2}},$$
(18)   
where the distance $`|𝐦𝐧|`$ is measured in lattice constants. For this form of the force, $`\lambda =1.742\kappa ^2/(2M\omega ^2D)`$. The model describes a two-dimensional particle interacting with a parallel plane of ions vibrating perpendicular to the plane. It was proposed in Ref., where it was used to model the interaction of holes doped in copper-oxygen planes, with apical oxygens in the layered cuprates. It was found that this long-range (Fröhlich) polaron is much lighter than the short-range Holstein polaron. Here we present the density of states for $`\overline{\omega }=1.0`$ and $`\lambda =2.75`$, see Fig. 9. DOS shape is close to the free-particle one, with a single, well-developed singularity in the middle of the band. Note, that we are in the adiabatic regime, at the same frequency $`\overline{\omega }=1.0`$ where the two-dimensional Holstein polaron has a very distorted DOS, cf. Fig. 3. The comparison of Figs. 9 and 5 shows that a long-range electron-phonon interaction plays the same role as increasing phonon frequency, as far as the flattening effect is concerned.

## IV Discussion and conclusions

The main message of this paper is that a path-integral imaginary-time Quantum Monte Carlo is quite capable of direct measuring of real-time spectra. By relaxing the boundary conditions in imaginary time, one allows many-body paths to have arbitrary real-space shift $`\mathrm{}𝐫`$. Then the Fourier transform projects out configurations with a certain total momentum P. The result, expressed by the formula (14), is the ground-state energy as a function of the total momentum. The physical system, where such a ground-state dispersion is of interest, should be carefully chosen.

The role of temperature in this process is two-fold. On one hand, temperature should be made as low as possible to exclude the contribution from the excited states with the same P. On the other hand, only states with $`E_𝐏E_0k_BT`$ are excited in the system. This means that the corresponding configurations will be generated by the QMC process in amounts sufficient for a good statistics. For higher-energy states, configurations will be exponentially rare. This is the reason for the average cosine in Eq. (14) to become exponentially small in the low-temperature limit. In this case, the measurement process will be statistically unstable. Thus, the temperature should be of the order of the energy interval one is interested in. The two conditions on the temperature can be reconciled if the energy scale of the ground-state dispersion is much smaller than the energy gap between the ground state and the excited states within the same P-sector. This is realized in the polaron system where one has $`W\mathrm{}\omega `$ for a wide range of parameters. If this condition is not satisfied, Eq. (14) will be measuring the difference of projected free energies rather than that of ground-state energies.

Eq. (14) shows that in a many-body system there exists a general and simple relation between the ground-state dispersion and the end-to-end distribution of imaginary-time paths. It does not involve the energy estimator, although the latter is required for the separate calculation of $`E_0`$ and $`E_𝐏`$. This might be useful in cases when the evaluation of the energy estimator is computationally costly. There are other computational advantages. First, Eq. (14) involves only the logarithm of one measured quantity, the average cosine. Second, Eq. (14) calculates the difference of two energies both of which may be large. In the polaron problem, typical energies are of the order of a few $`t`$ but the bandwidth is $`W0.1t`$. Both $`E_0`$ and $`E_𝐏`$ can be calculated with typical accuracy $`0.30.5\%`$. This may result in a sizable error in their difference if the two energies are calculated separately and then subtracted. Eq. (14) produces much more stable energy differences because of the large cancellation of errors between $`E_𝐏`$ and $`E_0`$. In all the spectra presented in this paper, Figs. 1, 2, 4, the statistical errors are smaller that symbols representing the data. Finally, the whole dispersion (as well as $`E_0`$, and derivatives at any P-point) can be calculated during a single QMC run. This property of Eq. (14) also allows fast computation of the density of states.

To demonstrate the practical usefulness of Eq. (14), we have combined it with an exact continuous-time algorithm for the lattice polaron , and calculated first detailed polaron spectra and densities of states in two and three dimensions. Although our method of calculating the spectrum is limited to the condition $`W\mathrm{}\omega `$, i.e., to the intermediate and strong coupling, this is the most physically interesting regime. Together with the weak- and strong-coupling perturbation expansions, the exact diagonalization, density-matrix renormalization group , variational, and diagrammatic QMC techniques, the present method covers the whole parameter range of the polaron problems. One can say, that the problem of calculating the exact polaron spectrum has found its solution.

Apart from being a test for Eq. (14), the lattice polaron still has a considerable interest on its own. We have seen in this paper that the polaron spectrum in the adiabatic limit of the Holstein model is generically non-cosine, as was predicted theoretically and recently observed numerically. The flattening has been found to continue well into the strong-coupling regime, where the polaron mass is a few dozens and approaches a hundred. The same conclusion was reached previously in . For physical applications the most interesting regime is when polarons are not very heavy and can be mobile. We conclude that the non-cosine spectrum is typical for the Holstein model in the physically relevant parameter region, i.e., small phonon frequencies and intermediate couplings.

A surprise finding has been the great extent at which the flattening changes the band structure in high dimensions. Density of states is completely changed, van Hove singularities are shifted, low-lying states are almost irrelevant, relations between the effective mass and the bandwidth is broken. We have seen that the combination of dimensionality, phonon frequency, coupling strength, and anisotropy may produce densities of states of various shapes. In such a situation, one should be careful about the interpretation of any experimental data, like the estimation of $`\lambda `$ or polaron mass from the bandwidth or from the location of singularities. The use of a “naive” model band structure would lead to wrong conclusions.

We have also found that in all dimensions the polaron band structure becomes free-particle-like with increasing phonon frequency. In particular, the spectrum approaches the cosine shape and the effective mass approaches the inverse of the bandwidth. Moreover, numerical values of $`W`$ are well described by the Lang-Firsov approximation. Thus, our QMC data support the LF approximation as the right description of the polaron in the antiadiabatic regime. At the same time, we have found that LF could still overestimate the polaron mass by a few dozen percent in cases where the bandwidth is predicted correctly.

Finally, we considered a long-range electron-phonon interaction and found a free-particle-like band structure even in the adiabatic regime. In Ref.  it was found that in the adiabatic regime polaron masses are well described by the LF approximation. Two conclusions follow from these facts. First, all the unusual properties of the Holstein model caused by the flattening, may be specific to the local electron-phonon interaction and may not be generic polaron properties. Second, the long-range electron-phonon interaction on a lattice seems to have the same effect on the band structure as the increasing phonon frequency. Although it is clear intuitively that a long-range interaction leads to higher mobility of the lattice deformation, details of this mechanism are yet to be fully understood.

The author is grateful to A. S. Alexandrov, D. M. Ceperley, V. Elser, W. M. C. Foulkes, and V. V. Kabanov for useful discussions and communications. This work was supported by EPSRC under grant GR/L40113.