text:
# Holstein polaron in two and three dimensions by quantum Monte Carlo

## I Introduction

In recent years, a lot of experimental evidence has been given for the importance of electron-phonon interaction in strongly correlated systems such as the cupratesBar-Yam et al. (1992) or the manganites.Edwards (2002) Although considerable theoretical progress has been made in understanding and describing many details of the physics of these compounds, the quantum nature of the phonons has often been neglected in actual calculations.Millis et al. (1996); Dagotto et al. (1999) However, a quantum-mechanical treatment of the lattice degrees of freedom has been shown to give significantly different results in certain parameter regimes.Edwards (2002)

Due to the complexity of the models for the abovementioned classes of materials, numerical methods have been used extensively. One very powerful approach is the quantum Monte Carlo (QMC) method, which allows for simulations on relatively large lattices and gives quasiexact results (i.e., exact apart from statistical errors which can, in principle, be made arbitrarily small) also at finite temperature. The latter point represents a major advantage over, e.g., the density matrix renormalization group (DMRG) or variational methods—which are restricted to the calculation of ground-state properties—since fascinating phenomena such as high-temperature superconductivity and colossal magnetoresistance can be investigated. Exact diagonalization (ED) can be performed at finite temperatures (see, e.g., Ref. Aichhorn et al., 2003), but for the case of electron-phonon systems such calculations would be restricted to very small systems and unrealistic parameters. Consequently, all existing work based on ED is for zero temperature only. Additionally, for coupled electron-phonon systems, the infinite-dimensional Hilbert space associated with the boson degrees of freedom represents a substantial difficulty for ED and DMRG, in contrast to QMC. Nevertheless, QMC methods are often limited by (1) the minus-sign problem, which restricts simulations to high temperatures and/or small systems, (2) the fact that the calculation of dynamical properties, such as the one-electron Green function, requires analytic continuation to the real-time axis which is an ill-posed problem, and (3) by strong autocorrelations and large statistical errors.

In a recent paper,Hohenadler et al. (2004) from here on referred to as I, we have proposed a new QMC approach to the Holstein model, which is based on the canonical Lang-Firsov transformationLang and Firsov (1962) and a principal component representation of the phonon degrees of freedom. The resulting algorithm has been shown to completely avoid the problem of autocorrelations, strongly reduce statistical errors and to allow for very efficient simulations.Hohenadler et al. (2004) The work here and in I is restricted to the case of a single electron which is often called the “Holstein polaron problem”. Obviously, in connection with materials such as the cuprates or manganites, we are interested in studying the many-electron system. However, the one-electron case represents an important first step towards an understanding of more general situations. Moreover, due to the coupling to the phonons, even the model with one electron constitutes a complex many-body problem which has received a lot of attention in the past. An important feature of our approach is the fact that it can be generalized to the more realistic many-electron case encountered in transition metal oxides.

In this paper, we extend the QMC method proposed in I to square (cubic) lattices in two (three) dimensions. Although the one-dimensional case considered in I is of great theoretical interest, real strongly correlated materials usually require two or three dimensional models. Consequently, before turning to the more challenging many-electron case, it is crucial to examine the influence of dimensionality on our algorithm. To this end, we investigate the dependence of the minus-sign problem, which has been found to appear in simulations of the Lang-Firsov transformed model, on the dimensionality D of the lattice, as well as on the parameters of the model. We then use our method to study in detail the well-known small-polaron cross over, which occurs with increasing electron-phonon coupling strength, and its dependence on phonon frequency, temperature, system size and dimensionality. We remove the Trotter error, which results from the discretization of imaginary time, by extrapolating our results to the limit of continuous time, thereby significantly improving the accuracy of the data.

The paper is organized as follows. In Sec. II we review the physics of the Holstein model with one electron, and in Sec. III we briefly discuss the QMC method introduced in I. Results are presented in Sec. IV and Sec. V contains our conclusions.

## II The Holstein model

The Holstein HamiltonianHolstein (1959) with dimensionless phonon variables readsHohenadler et al. (2004)

$`H`$ $`=`$ $`t{\displaystyle \underset{ij}{}}c_i^{}c_j^{}+{\displaystyle \frac{\omega }{2}}{\displaystyle \underset{i}{}}\left(\widehat{p}_i^2+\widehat{x}_i^2\right)\alpha {\displaystyle \underset{i}{}}\widehat{n}_i\widehat{x}_i`$ (1)
$`=`$ $`K+P+I,`$
where we have introduced the abbreviations $`K`$, $`P`$ and $`I`$ for the kinetic, potential and interaction terms, respectively. Spin indices have been suppressed in the notation, owing to the spin symmetry of the one-electron problem. In Eq. (1), $`c_i^{}`$ ($`c_i`$) creates (annihilates) a spinless electron at lattice site $`i`$, $`\widehat{x}_i`$ and $`\widehat{p}_i`$ denote the displacement and momentum operator of a harmonic oscillator at site $`i`$, and $`\widehat{n}_i=c_i^{}c_i^{}`$. The coupling term $`I`$ describes the local interaction of the single electron considered here with dispersionless Einstein phonons. In the first term, the symbol $`ij`$ denotes a summation over all hopping processes $`ij`$ and $`ji`$ between neighboring lattice sites $`i,j`$. The parameters of the model are the hopping integral $`t`$, the phonon energy $`\omega `$ ($`\mathrm{}=1`$), and the electron-phonon coupling constant $`\alpha `$. As in I, we define the dimensionless coupling constant $`\lambda =\alpha ^2/(\omega W)`$, where $`W=4t\text{D}`$ is the bare bandwidth in D dimensions. We shall also use the dimensionless phonon frequency $`\overline{\omega }=\omega /t`$, also called the “adiabatic ratio”, and express all energies in units of $`t`$. Consequently, the model can be described by two independent parameters, $`\overline{\omega }`$ and $`\lambda `$. While the bulk of results in this paper has been calculated assuming periodic boundary conditions in real space, we shall also discuss the effect of the boundary conditions on the sign problem in Sec. IV.1.

Since the literature on the Holstein polaron problem is vast, we restrict our discussion to work in more than one dimension. A brief review of the one-dimensional model can be found in I. Moreover, we will focus on recent progress in the field, and on numerical methods. A comprehensive overview of earlier analytical work can be found, e.g., in the books of Alexandrov and MottAlexandrow and Mott (1995) and Mahan.Mahan (1990) The Holstein polaron in $`\text{D}>1`$ has been studied using a large variety of numerical techniques. In contrast to many perturbative approaches,Alexandrow and Mott (1995); Mahan (1990) the latter can also accurately describe the physically most interesting regime of small but finite phonon frequency ($`0<\overline{\omega }<1`$), and intermediate electron-phonon coupling ($`\lambda 1`$). Much information about the Holstein polaron has been obtained using QMC. De Raedt and LagendijkDe Raedt and Lagendijk (1982, 1983, 1984) applied Feynman’s path-integral technique to the lattice problem of Eq. (1). In this approach, the phonon degrees of freedom are integrated out analytically, and one is left with a single-particle system with a retarded self-interaction, which can be simulated using the Monte Carlo method. The only approximation consists of discretizing the imaginary time using the Suzuki-Trotter decomposition.De Raedt and Lagendijk (1983) As the phonons have been completely eliminated from the problem, very efficient simulations on large lattices even in higher dimensions are possible. Moreover, the method is not restricted to the Holstein Hamiltonian (1). It can be extended to include long-range electron-phonon coupling as well as dispersive phonons.De Raedt and Lagendijk (1984) KornilovitchKornilovitch (1997) later employed the same method to study the whole range of values of the adiabatic ratio $`\overline{\omega }`$, and extrapolated the results to remove the error due to the Trotter approximation. He also pointed out the connection of the algorithm to world-line methods.Kornilovitch (1997) Extending this work, Kornilovitch developed a QMC approach which is formulated in continuous imaginary time, and which is capable of directly measuring the polaron band dispersion $`E(𝒌)`$ and the density of states in one to three dimensions, by sampling world lines with open boundary conditions in imaginary time.Kornilovitch (1998, 1999) Although the method gives very accurate results and can also be applied to more general models with, e.g., long-range interaction, it is limited to the regime of intermediate and strong electron-phonon coupling as well as $`\overline{\omega }1`$ by a minus-sign problem.Kornilovitch (1998, 1999) Consequently, it is not as universally applicable as the algorithm of de Raedt and Lagendijk.De Raedt and Lagendijk (1982, 1983, 1984); Kornilovitch (1997) As the above discussion reveals, the QMC methods of Refs. De Raedt and Lagendijk, 1982, 1983, 1984; Kornilovitch, 1997, 1998, 1999 are very well suited to study the one-electron problem. Moreover, de Raedt and Lagendijk also extended their approach to the bipolaron problem of two electrons with opposite spins.De Raedt and Lagendijk (1986) However, as pointed out by Kornilovitch,Kornilovitch (1997) these are all world-line methods. Consequently, despite the possibility of integrating out the phonon degrees of freedom even in the many-electron case, they face a serious sign problem in more than one dimension, for two or more fermions of the same spin, and can therefore not be used to study many-particle systems. In the context of superconductivity, the Holstein and the Holstein-Hubbard model with many electrons have been investigatedBlankenbecler et al. (1981); Scalapino and Sugar (1981); Levine and Su (1990, 1991); Niyaz et al. (1993) using the grand-canonical determinant method of Blankenbecler et al. Blankenbecler et al. (1981) (see also I). However, as discussed in I, due to strong autocorrelations, these simulations were restricted to rather large phonon frequencies $`\overline{\omega }1`$, while, e.g., the cuprates and the manganites fall into the adiabatic regime $`\overline{\omega }1`$. This is exactly the point where our new approach enters. As we shall see below, it completely avoids the problem of autocorrelations.

Apart from QMC, several other methods have been successfully applied to the Holstein polaron. This includes ED in combination with a truncation of the phonon Hilbert space,Fehske et al. (1997); Wellein et al. (1996); Wellein and Fehske (1997) finite-cluster strong coupling perturbation theory (SCPT),Stephan (1996) cluster perturbation theory (CPT),Hohenadler et al. (2003) DMRG,Jeckelmann and White (1998) and a variational diagonalization method.Ku et al. (2002) As pointed out before, standard ED is restricted to rather small numbers of lattice sites and/or phonon states, while DMRGJeckelmann and White (1998) and the variational method of Ref. Ku et al., 2002 are applicable over a wider range of values of phonon frequency and electron-phonon coupling, and are much less influenced by finite-size effects. The same is true for SCPT and CPT, which exactly diagonalize small clusters—for which enough phonon states can be included in the calculation—and extrapolate the results to the thermodynamic limit by treating the rest of the system as a perturbation.Stephan (1996); Hohenadler et al. (2003) Nevertheless, as mentioned before, all this work was limited to zero temperature. Finally, the Holstein polaron has been investigated recently using weak- and strong-coupling perturbation theory,Romero et al. (1999a, b) and a variational wave function which is a superposition of Bloch states for the weak- and the strong-coupling regime.Cataudella et al. (1999)

From all this work, many properties of the Holstein polaron are well understood. Similar to one dimension,Hohenadler et al. (2004) there is a cross over from a quasiparticle with slightly increased effective mass to a heavy small polaron as the electron-phonon coupling strength increases. As pointed out by Fehske et al. Fehske et al. (1997) and Capone et al.,Capone et al. (1997) the two conditions for the existence of a small polaron are $`\lambda >1`$ and $`\lambda \text{D}/\overline{\omega }>1`$. We will first discuss the weak-coupling state. Concerning its nature in higher dimensions, there exist two different views. From calculations based on the adiabatic approximation, i.e., taking the limit $`\overline{\omega }0`$, one expects a qualitatively different behavior in one dimension compared to $`\text{D}>1`$, also for $`\overline{\omega }>0`$. Wellein et al. Wellein et al. (1996) distinguish between the adiabatic ($`\overline{\omega }<1`$) and the nonadiabatic ($`\overline{\omega }>1`$) regime. In the adiabatic case, and for $`\text{D}>1`$, the electron is expected to remain quasifree with an almost unchanged effective mass and kinetic energy for $`\lambda <1`$, corresponding to a quasiparticle with infinite radius. This contrasts strongly with the one-dimensional case, in which the electron is always self-trapped by the surrounding lattice distortion and forms a polaron with finite radius for any $`\lambda >0`$. On the other hand, in the nonadiabatic regime $`\overline{\omega }>1`$, the behavior is very similar in all dimensions, and a very gradual decrease of the kinetic energy (or increase in effective mass) is observed as the coupling strength increases. Due to the large energy of the phonon excitations, only the zero-phonon state contributes significantly. Moreover, the phonons are fast and react almost instantly to the motion of the electron. Consequently, a lattice distortion only persists in the immediate vicinity of the electron, and this rather small quasiparticle is sometimes called a “nonadiabatic Lang-Firsov polaron”.Wellein et al. (1996) Romero et al. Romero et al. (1999b) take on a slightly different viewpoint. They argue that the large polaron state is essentially the same in any dimension, and that the only effect of increasing the dimension of the system is given by the observed sharpening of the transition to a small polaron in $`\text{D}>1`$ compared to 1D. For $`\lambda `$ larger than a critical value $`\lambda _\text{C}`$, determined by the aforementioned conditions $`\lambda >1`$ and $`\lambda \text{D}/\overline{\omega }>1`$, where the cross over occurs, both views agree on the existence of a a small-sized so-called “Holstein polaron”.Wellein et al. (1996) The latter is a heavy quasiparticle with a strongly reduced mobility. The cross over at $`\lambda _\text{C}`$ is very sharp, especially for $`\overline{\omega }1`$, but it does not represent a real phase transition.Lowen (1988) Even though the electron is trapped in the potential well originating from the response of the lattice to its motion, the ground state is still Bloch-like. For simplicity, in the remainder of this paper, we shall always refer to the weak-coupling state as a large polaron, either with finite or infinite radius, depending on which of the abovementioned viewpoints one holds. In our opinion, a definite decision about which of the two alternatives is correct cannot be made using QMC, which is restricted to finite clusters and, more importantly, finite temperatures. As a consequence, there will always be a contribution from excited states, making it difficult to reveal the true nature of the ground state in the weak-coupling regime.

## III  Quantum Monte Carlo method

Since the extension of the QMC algorithm to higher dimensions is straight forward, here we shall merely give an overview of the method which has been discussed in detail in I.

### III.1 QMC algorithm

The cornerstone of the new approach is the canonical Lang-Firsov transformation,Lang and Firsov (1962) which separates the polaron effects, due to the electron-phonon interaction, from the zero-point fluctuations of the harmonic oscillators in Eq. (1). The transformed model with one electron takes the formHohenadler et al. (2004)

$$\stackrel{~}{H}=\stackrel{~}{K}+PE_\text{P},\stackrel{~}{K}=t\underset{ij}{}c_i^{}c_j^{}e^{\mathrm{i}\gamma (\widehat{p}_i\widehat{p}_j)}$$
(2)
with $`P`$ as defined in Eq. (1), and the polaron binding energy $`E_\text{P}=\lambda W/2`$. The parameter $`\gamma `$, which corresponds to the displacement of the harmonic oscillator in the presence of an electron,Hohenadler et al. (2004) is given by $`\gamma ^2=2E_\text{P}/\omega `$. The method employs a Trotter decomposition of the imaginary time axis into $`L`$ intervals of size $`\mathrm{\Delta }\tau =\beta /L`$, where $`\beta =(k_\text{B}T)^1`$ is the inverse temperature. The partition function, obtained by integrating out the phonon coordinates,Hohenadler et al. (2004) is given by

$$𝒵_L=\text{const.}𝒟pw_\text{b}w_\text{f},$$
(3)
where $`𝒟p`$ denotes the $`L\times N^\text{D}`$ dimensional integral over all phonon momenta $`p`$, and $`N`$ is the linear size of the lattice in D dimensions. The bosonic weight $`w_\text{b}`$ is defined as $`e^{\mathrm{\Delta }\tau S_\text{b}}`$ with the bosonic action

$$S_\text{b}=\underset{i=1}{\overset{N}{}}𝒑_i^\text{T}A𝒑_i,$$
(4)
$`𝒑_i=(p_{i,1},\mathrm{}p_{i,L})`$, and a tridiagonal $`L\times L`$ matrix $`A`$ with nonzero elements

$$A_{jj}=\frac{\omega }{2}+\frac{1}{\omega \mathrm{\Delta }\tau ^2},A_{j,j\pm 1}=\frac{1}{\omega \mathrm{\Delta }\tau ^2},$$
(5)
and periodic boundary conditions $`L+11`$. As discussed in I, the electronic weight $`w_\text{f}`$ is given by

$$w_\text{f}=\text{Tr}_\text{f}\mathrm{\Omega },\mathrm{\Omega }=\underset{\tau =1}{\overset{L}{}}e^{\mathrm{\Delta }\tau \stackrel{~}{K}_\tau }.$$
(6)
Here $`\stackrel{~}{K}_\tau `$ is $`\stackrel{~}{K}`$ with the phonon operators $`\widehat{p}_i`$ replaced by the values $`p_{i,\tau }`$ on the $`\tau `$th Trotter slice. The exponential of the hopping term can be written as

$`e^{\mathrm{\Delta }\tau \stackrel{~}{K}_\tau }`$ $`=`$ $`D_\tau \kappa D_\tau ^{}`$ (7)
$`\kappa _{jj^{}}`$ $`=`$ $`\left(e^{\mathrm{\Delta }\tau th^{\text{tb}}}\right)_{jj^{}},(D_\tau )_{jj^{}}=\delta _{jj^{}}e^{\mathrm{i}\gamma p_{j,\tau }},`$
where $`h^{\text{tb}}`$ is the $`N^\text{D}\times N^\text{D}`$ tight-binding hopping matrix for the lattice under consideration. In fact, this is the only nontrivial change compared to the one-dimensional case. As pointed out in I, for a single electron, the fermionic trace can easily be calculated from the matrix representation of $`\mathrm{\Omega }`$ as $`w_\text{f}=_i\mathrm{\Omega }_{ii}`$. Due to the complex-valued hopping term, $`w_\text{f}`$ is not strictly positive, which gives rise to the minus-sign problem discussed below. We would like to mention that, in contrast to some determinant QMC methods,von der Linden (1992) the $`L`$-fold matrix product involved in the calculation of the matrix $`\mathrm{\Omega }`$ is well conditioned also for large systems at low temperatures, so that a time-consuming numerical stabilization is not necessary.

In I, we have introduced the so-called principal component representation for the phonon degrees of freedom. In terms of the latter, the bosonic weight takes the simple Gaussian form

$$w_\text{b}=e^{\mathrm{\Delta }\tau _i𝝃_i^\text{T}𝝃_i}$$
(8)
with the principal components $`𝝃_i=A^{1/2}𝒑_i`$. It is strictly positive. Finally, the efficiency of the QMC algorithm can be greatly improved by the use of a reweighting of the probability distribution. In our case, this amounts to transferring all the influence of the electronic degrees of freedom—which are treated exactly—to the observables calculated via

$$O=\frac{Ow_\text{f}_\text{b}}{w_\text{f}_\text{b}}$$
(9)
with the expectation values with respect to the bosonic weight $`w_\text{b}`$ defined as

$$O_\text{b}=\frac{𝒟pw_\text{b}O(p)}{𝒟pw_\text{b}}.$$
(10)
As outlined in I, the expectation values defined in Eq. (10) can be determined during the QMC simulation. The applicability of the reweighting method has been discussed in detail in I, and here we only mention that we have carefully checked that it also works reliably in higher dimensions. This could be expected, since the physics of the Holstein polaron is rather similar in all dimensions (see Sec. II). Details about the calculation of observables within our formalism can be found in I. Due to the analytic integration over the phonon coordinates $`x`$ used here, interesting observables such as the correlation functions $`_i\widehat{n}_i\widehat{x}_{i+\delta }`$ are difficult to measure accurately. Other quantities such as the quasiparticle weight, and the closely related effective mass,Wellein et al. (1996) can be determined from the one-electron Green function at long imaginary times,Brunner et al. (2001) but results but would not be as accurate as existing work on the one-electron case considered here (see, e.g., Refs. Jeckelmann and White, 1998; Romero et al., 1999b, and Ku et al., 2002). Therefore, we have restricted ourselves to the kinetic energy of the electron, which contains a lot of information about the small-polaron cross over. For the more demanding many-electron case, to which our method can be extended (see Sec. V), other methods produce far less reliable data.

We are now in a position to discuss the actual QMC procedure which has been explained more thoroughly in I. From the above discussion and I, it is obvious that we only simulate the phonons, while the electronic degrees of freedom have been integrated out analytically. This is equivalent to the method proposed by Blankenbecler et al. Blankenbecler et al. (1981) for the many-electron case. However, here the Monte Carlo sampling is independent of the fermionic weight $`w_\text{f}`$ given by Eq. (6), with the latter entering the calculations only via the reweighting procedure \[see Eq. (9)\]. Thereby, we avoid the problem of nonpositive weights during the QMC updates, since we do not use $`w_\text{f}`$ as a probability for accepting or rejecting new phonon configurations. Nevertheless, the sign problem still manifests itself in terms of statistical errors, as can be seen from Eq. (9). Similar to other occurrences of the minus-sign problem, e.g., for the case of the Hubbard model away from half filling,von der Linden (1992) simulations would become very difficult if $`w_\text{f}_\text{b}`$ should tend to zero.

Owing to the Gaussian form of the bosonic action, when written in terms of the principal components $`\xi `$, the latter can be sampled exactly by drawing random numbers from a normal distribution. In contrast to usual Markov chain Monte Carlo simulations, every new configuration is accepted, and no autocorrelations between successive values of the $`\xi `$ exist. This is in strong contrast to conventional QMC methods for the Holstein model, also with many electrons (see discussion in I). Except for situations in which the phonons are integrated out analytically, simulations become extremely difficult at low temperatures and for small phonon frequencies due to strongly increasing autocorrelations. We regard the complete absence of such correlations as a major advantage of our method. After the $`\xi `$ have been updated for each site of the $`L\times N^\text{D}`$ space-time lattice, a transformation back to the momenta $`p`$ is performed using the matrix $`A^{1/2}`$. Then, for each observable of interest, $`Ow_\text{f}`$ as well as the fermionic weight $`w_\text{f}`$ are calculated for the current phonon configuration, after which the next update can begin. At the end of the program, results for observables are obtained according to Eq. (9) (see also I).

### III.2 Sign problem

Before we come to a discussion of the sign problem in the approach presented here, we would like to give a quick review of its occurrence in other QMC methods. The situation is best illustrated for the case of the world-line algorithm (see, e.g., Ref. Batrouni and Scalettar, 1998). The use of the latter to simulate systems of interacting fermions is restricted to one dimension by the Fermi statistics of the electrons. This is a consequence of the negative matrix elements $`w`$, which appear when two fermion world lines wind around each other one or more times as they traverse the space-time lattice. With other methods, despite the occurrence of the sign problem, simulations can still be carried out in many situations. Since the QMC sampling requires nonnegative weights, one uses $`|w|`$ instead of $`w`$. As a consequence, the sign of the fermionic weight, $`\text{sgn}(w)`$, is treated as part of the observables. The latter then have to be calculated via

$$O=\frac{O\text{sgn}(w)_{|w|}}{\text{sgn}(w)_{|w|}},$$
(11)
and it is obvious that simulations will become extremely difficult if the denominator in Eq. (11) tends to zero. In fact, if a sign problem is present, the number of measurements during a QMC run has to be increased by a factor $`\text{sign}^2`$ to obtain results of the same accuracy.von der Linden (1992) As discussed in Ref. Chandrasekharan and Wiese, 1999, using $`|w|`$ instead of $`w`$ corresponds to simulating an effective model of hard-core bosons. The average sign can be written as

$$\text{sign}=e^{\beta V(f_wf_{|w|})},$$
(12)
where $`f_w`$ and $`f_{|w|}`$ denote the free energy per site of the fermionic and bosonic model, respectively, and $`V`$ is the volume of the system.Chandrasekharan and Wiese (1999) From Eq. (12), it is obvious that $`\text{sign}`$ decreases exponentially with increasing $`\beta `$ and $`V`$.

The auxiliary-field methodvon der Linden (1992) for the Hubbard model faces similar problems. Here the weight of a configuration is given by the product of determinants for $``$ and $``$ electrons, respectively. The product is strictly positive only for half filling, whereas simulations for other particle densities become very demanding at low temperatures and/or for large systems. Since at large $`U`$ the determinant is similar to the weight of world lines of Hubbard-Stratonovitch variables,Hirsch (1986) this similarity of the dependence of $`\text{sign}`$ on the parameters of the system is not surprising. Finally, there is no sign problem in determinantal grand-canonical simulations of the Holstein model at any filling, since the coupling of electrons and phonons is the same for both spin directions.Loh Jr. et al. (1990)

The sign problem in the current approach clearly has a fundamentally different origin, since there is only a single electron in the system, so that no winding of world lines around each other can take place. In fact, the negative fermionic weights are a result of the complex phase factor in the transformed hopping term \[Eq. (2)\]. We will see below that, in contrast to other methods, the sign problem in our approach diminishes, i.e., $`\text{sign}1`$, with increasing system size, and a detailed investigation of its dependence on the parameters of the model and the choice of boundary conditions will be given in Sec. IV.

### III.3 Numerical details and performance

As discussed in I, the error due to the Trotter decomposition is proportional to $`(\mathrm{\Delta }\tau )^2`$. We perform simulations at different values of $`\mathrm{\Delta }\tau `$, typically $`\mathrm{\Delta }\tau =0.1`$, $`0.075`$ and $`0.05`$, and exploit the linear dependence of the results on $`(\mathrm{\Delta }\tau )^2`$ to extrapolate to $`\mathrm{\Delta }\tau =0`$. This is a common procedure in the context of discrete-time QMC methods,von der Linden (1992) and allows one to remove the Trotter error if the values of $`\mathrm{\Delta }\tau `$ are sufficiently small. Moreover, a given accuracy can often be reached using larger values of $`\mathrm{\Delta }\tau `$, compared to calculations which do not use the abovementioned extrapolation, thereby saving computer time.

We conclude this section by a comparison of our approach with other QMC methods for the Holstein polaron. As mentioned in Sec. I, the methods of de Raedt and LagendijkDe Raedt and Lagendijk (1982, 1983, 1984) and KornilovitchKornilovitch (1997, 1998, 1999) are based on an analytic integration over the phonon degrees of freedom. This separation of electronic and bosonic degrees of freedom greatly reduces the statistical noise due to phonon fluctuations, which are induced by the electron-phonon interaction. The fluctuations increase noticeably with decreasing phonon frequency, decreasing temperature and increasing electron-phonon coupling. Although the Lang-Firsov transformation used here performs a very similar task, namely, to separate polaron effects from the zero-point and thermal fluctuations of the free oscillators (see I), the integral over the bosonic degrees of freedom is calculated with Monte Carlo, thereby leaving us with a residual influence of the phonons. In fact, the numerical effort for calculations with our approach is proportional to the $`LN^{3\text{D}}`$, similar to the grand-canonical determinant method for the Holstein model.Blankenbecler et al. (1981) This could be improved to a computer time $`N^{2\text{D}}`$ by linearizing the exponential of the hopping term to a tridiagonal matrix, which has not been done in the work presented here. Although our algorithm is not as efficient as the methods of Refs. De Raedt and Lagendijk, 1982, 1983, 1984; Kornilovitch, 1997, 1998, 1999, in which the numerical effort is independent of $`N`$, proportional to $`L^2`$ and depends only linearly on D, we will see in the following section that it is well suited for accurate simulations on large clusters in one and two dimensions, and on reasonably large clusters in 3D. Moreover, as we are interested in developing a QMC method that can also be applied to the many-electron system in the adiabatic regime—in order to study, e.g., quantum phonon effects in the manganites (see I)—such a decrease in performance when compared to the world-line methodsDe Raedt and Lagendijk (1982, 1983, 1984); Kornilovitch (1997, 1998, 1999) is acceptable, since the latter cannot be applied to the many-electron case in more than one dimension. Additionally, the advantage of the analytical integration over the phonons is biggest in the one-electron case, in which they represent the majority of the degrees of freedom in the path-integral representation of the partition function.Kornilovitch (1997) This is no longer true at finite band filling, where the contributions of fermions and bosons are comparable.

## IV  Results

This section is divided into two parts. First, in Sec. IV.1, we investigate in detail the dependence of the aforementioned sign problem on electron-phonon coupling, dimensionality, boundary conditions, system size, phonon frequency, and temperature. In Sec. IV.2, we present results for the electronic kinetic energy to study the small-polaron cross over discussed in Sec. II.

### IV.1 Sign problem

Here we are interested in the average sign of the fermionic weight $`w_\text{f}`$, which is given by

$$\text{sign}=\frac{w_\text{f}_\text{b}}{|w_\text{f}|_\text{b}},$$
(13)
where the expectation value $`\mathrm{}_\text{b}`$ has been defined in Eq. (10). In Fig. 1, we show its dependence on electron-phonon coupling and system size. We would like to point out that while the extrapolation to $`\mathrm{\Delta }\tau =0`$, discussed in Sec. III, has been used for the results shown in the following subsection, the calculations for the average sign have been performed for a single, fixed value $`\mathrm{\Delta }\tau =0.05`$, for which the Trotter error is smaller than statistical errors. Some other results for $`\text{sign}`$ in one dimension have already been reported in I.

From the general discussion in Sec. III, it is clear that the sign problem encountered in the present approach is of a different nature than in, for example, QMC simulations of the Hubbard model. As reported in I, for the Holstein polaron problem under consideration, it is most pronounced for small systems, low temperatures and small phonon frequencies $`\overline{\omega }1`$. Therefore, the bulk of results presented below will be for such a set of “worst case” parameters, including $`N=4`$, $`\beta t=10`$, and $`\overline{\omega }=0.1`$. Figure 1 shows that, in one dimension, the average sign of $`w_\text{f}`$ in the critical region of intermediate electron-phonon coupling increases quickly as the system size increases from $`N=4`$ to $`N=16`$, which is in strong contrast to Eq. (12). This increase of the minimum as a function of $`N`$ is also shown in the inset of Fig. 1. As we have only calculated $`\text{sign}`$ for a finite number of $`\lambda `$-values, an approximation for the minimum has been determined using a spline interpolation.

The minimum of $`\text{sign}`$ occurs near $`\lambda =1`$, where the transition from a large to a small polaron takes place (see Sec. II). The increase of statistical errors in this regime as a consequence of the sign problem is similar to the situation encountered with simulations of the untransformed model (see I). However, the use of the transformed model still gives significantly more accurate results for the same number of measurements, in particular for low temperatures and small phonon frequencies. Finally, one may be tempted to explain the unusual system-size dependence of the sign problem by ascribing its origin to the periodic boundary conditions in real space. If the latter were indeed the source of the sign problem, the boundary effects would decrease with increasing system size, in accordance with the results of Fig. 1. This possibility has been investigated, and we shall see below that the sign problem persists also for the case of open boundary conditions.

We now come to the Holstein model in two dimensions. In Fig. 2(a), we show results for $`\text{sign}`$ for different lattice sizes, again starting with a very small linear dimension $`N`$. All other parameters are the same as before, in particular $`\beta t=10`$ and $`\overline{\omega }=0.1`$. Obviously, for the smallest cluster size shown, the minimum of the average sign has diminished to a value of approximately 0.1, so that large numbers of measurements are necessary. However, similar to one dimension, $`\text{sign}`$ increases with increasing system size, and for the largest system size shown ($`N=12`$), we find a rather uncritical minimum value of about 0.5.

The results for open boundary conditions, shown in Fig. 2(b), reveal that for small clusters, the average sign increases compared to the case of periodic boundary conditions. However, with increasing system size, $`\text{sign}`$ quickly converges to the same values, independent of the boundary conditions. This is just what one would expect, since with increasing $`N`$, the effect of the choice of boundary conditions on the properties of the system diminishes. Moreover, we can conclude that the negative weights do not simply result from hopping processes of the electron across the periodic boundaries, since in that case we would expect $`\text{sign}1`$ for open boundary conditions, in contrast to Fig. 2(b). Similar behavior has been found in one and three dimensions.

The influence of the phonon frequency $`\omega `$ on the average sign is shown in Fig. 3(a). Clearly, the sign problem is most noticeable for small values of $`\overline{\omega }`$, while it diminishes quickly as we increase the adiabatic ratio (see inset). This is very similar to the small-polaron cross over. As discussed in Sec. II, the latter sharpens significantly with decreasing $`\overline{\omega }`$, while an abrupt transition is completely absent in the nonadiabatic regime $`\overline{\omega }>1`$. Moreover, the coupling $`\lambda _\text{C}`$ at which the minimum of $`\text{sign}`$ occurs increases with $`\overline{\omega }`$, in agreement with the aforementioned small-polaron condition $`\lambda \text{D}/\overline{\omega }>1`$.

In Fig. 3(b), we report the average sign as a function of $`\lambda `$, and for different inverse temperatures $`\beta `$. Again we have taken $`N=4`$ and $`\overline{\omega }=0.1`$, the parameters for which the sign problem is most noticeable. While for $`\beta t=10`$, the minimum of $`\text{sign}`$ lies below 0.1, the situation is much better already for $`\beta t=5`$, as shown in the inset. At even higher temperature $`\beta t=1`$, the fermionic weight is always positive so that we have $`\text{sign}=1`$ for all $`\lambda `$. The dependence of the sign problem on temperature is therefore similar to other QMC methods \[see Eq. (12)\], although we do not find a simple exponential relation.

Finally, we also present results for $`\text{sign}`$ in three dimensions, for lattices of different linear size $`N`$, and again for the parameters $`\overline{\omega }=0.1`$ and $`\beta t=10`$ as a function of $`\lambda `$. Figure 4 reveals that the minimum of $`\text{sign}`$ in 3D has an even more pronounced form than in two dimensions. The sign problem diminishes slightly as we increase the system size from $`N=4`$ to $`N=6`$. However, accurate simulations in this regime are still quite demanding. As we will see in Sec. IV.2, a temperature of $`\beta t=5`$ is sufficient to study the low-temperature properties of the Holstein model. For this higher temperature, and all other parameters unchanged, $`\text{sign}`$ is close to 1 even for $`N=4`$ (see inset in Fig. 4), so that accurate results can be obtained even for small phonon frequencies in three dimensions (see Sec. IV.2).

In summary, the investigation of the sign problem has shown that our method works well for a large range of values of phonon frequency, electron-phonon coupling and temperature, as long as the system size is large enough. This contrasts, for example, with the world-line algorithm of Kornilovitch,Kornilovitch (1998, 1999) which is restricted to intermediate and strong coupling, as well as low temperatures and rather large phonon frequencies (see also Sec. II), also by a minus-sign problem. Although we have investigated the influence of all important parameters on the sign problem, a physical interpretation of its origin has not emerged. Nevertheless, it is clear that the negative fermionic weights are a result of the phase factors in the Lang-Firsov transformed hopping term in Eq. (2). This is similar to the sign problem which occurs, e.g., in simulations of electron-phonon models in an external magnetic field,Loh Jr. et al. (1990) but here the phonon fields $`p_{i,l}`$ vary with time ($`l`$) and position ($`i`$), and are coupled in imaginary time by the bosonic action \[see Eqs. (4) and (5)\]. Moreover, the dependence of the sign problem on $`\overline{\omega }`$, $`\lambda `$, $`\beta `$ and $`N`$ bears striking resemblance to the influence of these parameters on the properties of the Holstein polaron. In particular, its reduction with increasing system size may be a consequence of the dilution of the system (the particle density $`n0`$ as $`N\mathrm{}`$, in the one-electron case). Finally, it is interesting to note that the large statistical fluctuations—resulting from the sign problem in the case of the transformed model—occur at exactly the same points in parameter space as in the untransformed model. This suggests a strong correlation between the minimum in $`\text{sign}`$ and the transition to a small polaron, which both occur at or near $`\lambda _\text{C}`$, similar to simulations of other models, in which the sign problem occurs exactly where the most interesting physics is going on, i.e., in the vicinity of phase transitions.von der Linden (1992)

### IV.2 Holstein polaron

The most interesting observable, which is easily accessible with our method (see discussion in Sec. III.1), and allows us to investigate the small-polaron cross over, is the one-electron kinetic energy $`E_\text{k}=K`$, given by the expectation value of the first term in Hamiltonian (1). In order to compare results for different dimensions, we define the normalized quantity

$$\overline{E}_\text{k}=E_\text{k}/(2t\text{D})$$
(14)
with $`\overline{E}_\text{k}=1`$ for $`T=0`$ and $`\lambda =0`$. Due to the large amount of work that has been devoted to the Holstein polaron in the past (see Sec. II), a lot is known about the transition from a large to a small polaron with increasing electron-phonon coupling. In this work, we therefore concentrate on those aspects which have not been studied in a systematic way so far. To this end, we exploit the advantages of our QMC method which allows us to investigate, in particular, finite-size and finite-temperature effects. The latter have only been touched upon briefly by previous authors employing QMC,De Raedt and Lagendijk (1982, 1983, 1984); Kornilovitch (1997, 1998, 1999) who focused on very largeDe Raedt and Lagendijk (1982, 1983, 1984); Kornilovitch (1997) or infinite systems,Kornilovitch (1998, 1999) in the ground stateKornilovitch (1998, 1999) or at two different temperatures.De Raedt and Lagendijk (1982, 1983) However, since a large amount of work has been done using ED or other methods based on diagonalization of small clusters, it is essential to study the convergence of the results with increasing system size. Additionally, we shall also present a comparison of the kinetic energy in one, two and three dimensions. As pointed out before, the results shown here have been obtained by extrapolating to $`\mathrm{\Delta }\tau =0`$, thereby eliminating the error due to the Trotter discretization.

#### IV.2.1 Two dimensions

To study the small-polaron cross over, previous authors focused on the effective mass of the electronKornilovitch (1998); Jeckelmann and White (1998); Romero et al. (1999b); Ku et al. (2002) and the quasiparticle weight.Wellein and Fehske (1997); Cataudella et al. (1999); Ku et al. (2002) Unfortunately, as pointed out in Sec. III.1, these observables cannot be calculated directly with our QMC algorithm. Nevertheless, we can examine many interesting aspects of the cross over by calculating the kinetic energy of the electron, and comparison will be made to existing work.De Raedt and Lagendijk (1982, 1983, 1984); Kornilovitch (1997); Jeckelmann and White (1998)

We begin with the dependence of the cross over on the phonon frequency. To this end, we present in Fig. 5 results for the kinetic energy calculated for $`N=12`$, $`\beta t=10`$ and different values of $`\overline{\omega }`$. The large range of the adiabatic ratio in Fig. 5, $`0.1\overline{\omega }8.0`$, shows the ability of our method to give accurate results for almost arbitrary values of the phonon frequency, especially in the adiabatic regime $`\overline{\omega }1`$. This is in contrast to the method of Kornilovitch,Kornilovitch (1998, 1999) which is restricted to $`\overline{\omega }1`$ and $`\lambda 1`$ by a severe minus-sign problem. In our case, the only limitations regarding the accessible values of $`\overline{\omega }`$ are the moderate sign problem, discussed in Sec. IV.1, which for small systems and low temperatures gives rise to a noticeable increase of statistical errors as $`\overline{\omega }0`$, and the increasing Trotter error as $`\overline{\omega }\mathrm{}`$, which requires the use of more and more time slices (see discussion in I).

Figure 5 also shows the well known fact that the transition from large to small polaron near $`\lambda =1`$ sharpens considerably with decreasing phonon frequency. While there is an abrupt decrease in $`\overline{E}_\text{k}`$ for $`\overline{\omega }=0.1`$, the cross over is very smooth for $`\overline{\omega }1`$. Additionally, we see from Fig. 5 that the cross over position $`\lambda _\text{C}`$ increases with $`\overline{\omega }`$ in the nonadiabatic regime. The physics of the transition to a small polaron has been discussed, e.g., by Capone et al. Capone et al. (1997) In the adiabatic regime ($`\overline{\omega }1`$), the cross over is entirely determined by the balance of kinetic and electron-phonon coupling- or displacement energy, where the latter is given by the polaron binding energy $`E_\text{P}`$. As soon as the gain in displacement energy outweighs the loss in kinetic energy, the electron localizes in a potential well and forms a polaron. The parameter $`\lambda `$ is defined as the ratio of these two contributions and may be written as $`\lambda =E_\text{P}/(W/2)`$ ($`W/2`$ is the kinetic energy of a free electron at $`T=0`$). Therefore, in the adiabatic regime, the cross over occurs at $`\lambda _\text{C}=1`$. With increasing $`\overline{\omega }`$, the lattice energy becomes more and more important, since more energy is required to excite phonons. As a consequence, the distortions of the lattice around the position of the electron—giving rise to the large effective mass and low mobility in the small polaron regime—are much smaller, and the local oscillators will predominantly be in their ground state. Thus, even for $`\lambda >1`$, where a trapped state is energetically favored in the adiabatic regime, the electron remains mobile. The decrease of $`\overline{E}_\text{k}`$ with increasing $`\lambda `$ is a result of the exponentially decreasing overlap of a displaced and an undisplaced harmonic oscillator in its ground state, which reduces the hopping matrix element between neighboring lattice sites. In the nonadiabatic regime, a small polaron is formed if $`\gamma ^2=E_\text{P}/\omega >1`$ (see Ref. Capone et al., 1997), where $`\gamma `$ is the induced lattice distortion at the site of the electron which follows from the Lang-Firsov transformation.Hohenadler et al. (2004) This is identical to the condition $`\lambda \text{D}/\overline{\omega }>1`$ given above. The larger lattice energy also gives rise to the more gradual decrease of $`\overline{E}_\text{k}`$ for intermediate and large values of $`\overline{\omega }`$. In particular, the kinetic energy is much larger for $`\lambda >1`$ and $`\overline{\omega }>1`$ than for $`\overline{\omega }<1`$.

Finally, like in one dimension,Hohenadler et al. (2004) $`\overline{E}_\text{k}`$ remains small but finite even at very strong coupling, which is a consequence of the fast and undirected motion of the electron inside the polaron. As pointed out by Kornilovitch,Kornilovitch (1997) the kinetic energy in the small polaron regime is therefore not related to the effective mass of the electron, the latter exhibiting an exponential decrease as a function of $`\lambda `$ (Refs. Kornilovitch, 1998Jeckelmann and White, 1998Ku et al., 2002, and Romero et al., 1999b).

To address the issue of finite-size effects, we have calculated $`\overline{E}_\text{k}`$ for $`\beta t=10`$, $`\overline{\omega }=0.1`$ and different linear lattice sizes $`N=4`$ – 12 \[see Fig. 6(a)\]. The choice of $`\overline{\omega }=0.1`$ is reasonable since the large polaron, which exists for $`\lambda <\lambda _\text{C}`$, is most extended for small phonon frequencies, as discussed in Sec. II, so that finite-size effects can be expected to be largest. To illustrate this point, we also present results for a larger phonon frequency $`\overline{\omega }=2.0`$ \[Fig. 6(b)\]. In the latter case, the local oscillators can respond very quickly to the motion of the electron, and the extension of the phonon cloud or lattice distortion surrounding the electron is much smaller. Since for $`\overline{\omega }>1`$ the transition to a small polaron happens at larger $`\lambda `$, the polaron will be larger for intermediate values $`1\lambda 2`$ compared to the adiabatic regime, in which $`\lambda _\text{C}=1`$. In total, we therefore expect to have smaller finite-size effects for $`\overline{\omega }=2.0`$ than for $`\overline{\omega }=0.1`$ as long as we are in the weak-coupling regime $`\lambda 1`$, while the opposite should be true for $`1\lambda 2`$. For $`\lambda >2`$, a highly immobile polaron state exists in both cases and results should therefore be virtually independent of $`N`$. All this is well confirmed by the results shown in Fig. 6.

Figure 6 also reveals that the results begin to saturate for $`N8`$, as pointed out previously by Kornilovitch.Kornilovitch (1997) However, in contrast to the one-dimensional case,Hohenadler et al. (2004) where we performed simulations for $`N`$ as large as 32, we find a nonnegligible dependence on $`N`$ up to the largest system size ($`N=12`$). This is better illustrated in Fig. 7, in which we show $`\overline{E}_\text{k}`$ as a function of $`1/N`$, and for several values of $`\lambda `$. To allow for a better representation, some curves have been shifted, as indicated in the legend. From Fig. 7, we see that $`\overline{E}_\text{k}`$ changes very little for $`N>4`$ for the case of strong coupling $`\lambda =2`$ and, in fact, remains constant within the error bars, similar to the 1D results in I. For smaller values of $`\lambda `$, no such saturation is found on the scale of Fig. 7, and the behavior of the kinetic energy for large $`N`$ is almost linear when plotted as a function of $`1/N`$. We used a linear fit of the data for $`N=8`$, 10, and 12 to obtain an approximation to the thermodynamic limit. The results are shown in Fig. 6. Obviously, $`\overline{E}_\text{k}`$ for $`N=\mathrm{}`$ has decreased noticeably for small values of $`\lambda `$ (including $`\lambda =0`$), while it remains almost unchanged in the small polaron regime. The decrease of $`\overline{E}_\text{k}`$ for $`\lambda <\lambda _\text{C}`$ can easily be understood if we consider the fact that our method works at a finite temperature $`1/\beta `$. As a consequence, for very small $`N`$, the energy gap between the ground state with $`𝒌=0`$ and the first excited state with $`𝒌0`$ is larger than the thermal energy $`(\beta t)^1`$. With increasing system size, thermal population of excited states becomes possible. For $`N`$ large enough ($`N20`$ for $`\beta t=10`$ and $`\lambda =0`$), results converge to those for $`N=\mathrm{}`$ and, in fact, the extrapolated data for $`\lambda =0`$, shown in Fig. 7, agree well with the results for a free electron on an infinite lattice. We ascribe the smallness of this finite-temperature effect in the strong-coupling regime above $`\lambda _\text{C}`$ to the very small width of the polaron band (see, e.g., Ref. Hohenadler et al., 2003). Consequently, the low-energy coherent states with different $`𝒌`$ have very similar energies.

Figure 6(a) also shows that for small phonon frequencies, the influence of the lattice size is largest near $`\lambda _\text{C}`$. The cross over is sharper for small systems than for the case of large $`N`$. A similar behavior has been found by Marsiglio for the one-dimensional model.Marsiglio (1995)

Finally, as pointed out by Kornilovitch, Kornilovitch (1997) despite the good convergence of the results for quantities such as the kinetic energy even for small $`N4`$, other observables, such as the effective mass, rely on the knowledge of the energy of states with infinitely small momenta. Therefore, they will show a more pronounced dependence on $`N`$ when calculated on small clusters accessible with, e.g., ED.

We also investigated the effect of temperature on our results, again for $`\overline{\omega }=0.1`$ and $`\overline{\omega }=2.0`$. The results, shown in Fig. 8, indicate that $`\overline{E}_\text{k}`$ is more affected by the finite temperature of the simulation in the adiabatic case $`\overline{\omega }=0.1`$ \[Fig. 8(a)\]. This is a consequence of the fact that calculations at finite temperatures only give ground-state-like results when $`\beta \omega 1`$. Clearly, this condition is much more difficult to meet for $`\overline{\omega }=0.1`$, and requires larger values $`\beta t>10`$.

The changes of $`\overline{E}_\text{k}`$ with temperature result from an interplay of several effects. For $`\lambda =0`$, the kinetic energy approaches its full noninteracting value of $`2t\text{D}`$ (i.e., $`\overline{E}_\text{k}=1`$) as $`T0`$. At finite $`T`$, however, states with nonzero total quasimomentum $`𝒌`$ will contribute and thereby lead to a decrease of $`\overline{E}_\text{k}`$. As discussed above, this effect of temperature on $`\overline{E}_\text{k}`$ is expected to decrease with increasing $`\lambda `$, and to be extremely small in the strong-coupling regime. In the adiabatic case \[Fig. 8(a)\], the transition at $`\lambda _\text{C}=1`$ is smeared out at high temperatures. For both $`\overline{\omega }=0.1`$ and $`\overline{\omega }=2.0`$, a qualitative change in behavior occurs near $`\lambda _\text{C}`$. $`\overline{E}_\text{k}`$ decreases with increasing temperature for $`\lambda <\lambda _\text{C}`$, whereas the opposite is true for $`\lambda >\lambda _\text{C}`$. The behavior above $`\lambda _\text{C}`$ can be understood by considering the electronic hopping amplitude, given by the overlap of the wave functions of a displaced and an undisplaced oscillator at neighboring sites. While the latter is exponentially reduced with increasing electron-phonon coupling at $`T=0`$ (see above), it increases with temperature since the oscillators can occupy excited states, corresponding to wave functions which are more spread out than the ground state. Moreover, the thermal energy allows the electron to overcome the potential barrier more easily. This thermally activated hopping—which shows up more clearly for small phonon frequencies since, in this case, phonon excitations require less energy—has been studied already a long time ago by Holstein.Holstein (1959)

It is interesting to notice that the dependence of $`\overline{E}_\text{k}`$ on temperature, as shown in Fig. 9, is almost linear at low temperatures and for $`\overline{\omega }=2.0`$ \[see Fig. 9(b)\], while in the adiabatic case \[Fig. 9(a)\] this is only true for strong electron-phonon coupling. However, as pointed out before, for $`\overline{\omega }=0.1`$, an inverse temperature $`\beta t=10`$ is not sufficient to obtain well-converged results. Therefore, a linear dependence, as for $`\overline{\omega }=2.0`$, may still be found at lower temperatures. Such calculations, with the accuracy of the results presented here, would be more time-consuming. The data for $`N=12`$, $`\beta =10`$, $`\mathrm{\Delta }\tau =0.05`$, and $`\lambda =1.0`$, i.e., in the cross over regime where statistical errors are largest (about 0.5%), required about 28 days of CPU time on an Intel Xeon 2600 MHz computer. However, as has been mentioned previously, the numerical effort may be reduced by a factor of order $`N^\text{D}`$ by linearizing the exponential of the hopping matrix, denoted as $`\kappa `$ (see Sec. III.3.) Similar to the finite-size scaling performed above, we also linearly extrapolated the data for $`\beta t=10`$ and 15 to the zero-temperature limit $`\beta t=\mathrm{}`$. The results for $`\overline{\omega }=2.0`$, for which such a linear scaling is reasonable, are shown in Fig. 8(b). While the general trend agrees well with our expectations based on the finite-temperature data shown, the scaling procedure clearly overestimates the temperature effects, thereby leading to spurious values $`\overline{E}_\text{k}>1`$ at $`\lambda =0`$. However, this can easily be understood keeping in mind that results saturate at low-enough values of $`\beta t`$, so that the linear extrapolation used here becomes insufficient.

Finally, we would like to mention that the results of de Raedt and LagendijkDe Raedt and Lagendijk (1982, 1983, 1984) were also given for $`\beta t=5`$ (and $`\overline{\omega }=1`$), but the small number of Trotter slices ($`L=32`$ in their case) gives rise to relatively large systematic errors.Kornilovitch (1997) This is not the case for the results of Ref. Kornilovitch, 1997, in which the same extrapolation to $`\mathrm{\Delta }\tau =0`$ was employed as here. Finally, the method of KornilovitchKornilovitch (1997, 1998) is free of such errors, but only permits one to calculate ground-state properties for a restricted range of $`\overline{\omega }`$ and $`\lambda `$.

#### IV.2.2 Three dimensions

In contrast to the two-dimensional case discussed above, less work has been done in three dimensions.De Raedt and Lagendijk (1982, 1983, 1984); Kornilovitch (1998, 1999); Romero et al. (1999a, b); Ku et al. (2002) In fact, we are only aware of one calculation of the kinetic energy, which is by de Raedt and Lagendijk (DRL).De Raedt and Lagendijk (1983) To compare with their work, we chose the same values for the phonon frequency $`\overline{\omega }=1.0`$ and temperature $`\beta t=5`$. As pointed out in Sec. III, the numerical effort for calculations with our method, which is proportional to $`N^{3\text{D}}`$ for the algorithm in the form used here but could be reduced to $`N^{2\text{D}}`$, restricts us to smaller systems than those considered by DRL.De Raedt and Lagendijk (1982, 1983, 1984) For simplicity, we have therefore limited ourselves to a maximum of $`N=6`$, for which results can easily be obtained within a reasonable amount of computer time, while the data presented in Refs. De Raedt and Lagendijk, 1982, 1983, 1984 is for $`N=32`$. To be more specific, our calculations for one value of $`\lambda `$, for $`N=6`$ and $`\mathrm{\Delta }\tau =0.05`$, took about 10 h on an Intel Xeon 2600 MHz computer. Due to the relatively small system size in our work, it is important to study to what extent the results are converged with respect to $`N`$. To this end, in Fig. 10, we present $`\overline{E}_\text{k}`$ as a function of $`\lambda `$ for $`N=4,5`$, and 6. Surprisingly, the results are already satisfactorily converged, keeping in mind the rather small values of $`N`$. There is a maximal change of less than 20% in the transition region at $`\lambda =1`$, while $`\overline{E}_\text{k}`$ remains almost constant for small and large $`\lambda `$, as the linear size increases from $`N=4`$ to $`N=6`$. Thus, increasing $`N`$ further will not change the results qualitatively, although the finite temperature of our simulations will manifest itself in a way similar to the two-dimensional case. Our findings agree well with the results of DRL.De Raedt and Lagendijk (1983) The main difference is that for weak coupling, our results are closer to the zero-temperature values (e.g., $`\overline{E}_\text{k}=1`$ at $`\lambda =0`$). The reason for this discrepancy—despite the fact that we have used the same temperature—is the smaller lattice size in our calculations. In contrast to the two-dimensional case considered above, we have not performed a scaling to $`N=\mathrm{}`$ in 3D, since the clusters under consideration are too small to reveal a systematic power law dependence as a function of $`1/N`$.

Finally, we wish to investigate the effect of dimensionality on the small-polaron cross over. Therefore, we compare $`\overline{E}_\text{k}`$ in one, two and three dimensions using $`\beta t=5`$, $`N=6`$, and $`\overline{\omega }=1.0`$. The dependence on D, which is shown in Fig. 11, is in perfect agreement with previous work. The transition from a mobile large polaron to a small polaron, moving in a very narrow but still coherent band, sharpens considerably with increasing dimension of the system, and while $`\overline{E}_\text{k}`$ only displays a gradual decrease in 1D—without any signs of an abrupt change at $`\lambda =1`$—we find a sharp and well-defined transition in three dimensions.

We conclude this section by comparing the accuracy of our results with the QMC methods of DRL,De Raedt and Lagendijk (1982, 1983, 1984); Kornilovitch (1997) and Kornilovitch.Kornilovitch (1998, 1999) As discussed in Sec. III, their main advantage over our method is the fact that they allow one to obtain data which are essentially free of finite-size effects, in any dimension $`\text{D}=1`$ – 3, and with modest computational effort. However, we have seen above that even in three dimensions, where the limitation of our algorithm is most noticeable, results are reasonably converged. While the approach of Refs. Kornilovitch, 1998 and Kornilovitch, 1999 is limited to $`T=0`$, DRL’s method as well as the current approach can, in principle, be used to study any temperature. Apart from the sign problem discussed in Sec. III.2, the only limitation which occurs is the fact that one has to increase the number of Trotter slices as $`\beta \mathrm{}`$, so as to keep the Trotter error smaller than the statistical errors. This situation can be greatly improved by extrapolating results to $`\mathrm{\Delta }\tau =0`$, as has been done in this work, while the continuous-time algorithm of KornilovitchKornilovitch (1998, 1999) is free of any such discretization errors. The accuracy of the results presented here depends on $`\overline{\omega }`$, $`\beta t`$, $`N`$ and $`\lambda `$. Away from $`\lambda 1`$, error bars are use usually smaller than the line width, corresponding to relative errors of less than 0.5%. This is comparable to the accuracy of the results given by KornilovitchKornilovitch (1998) and significantly more accurate than the original results by DRL.De Raedt and Lagendijk (1982, 1983)

## V Conclusions

Extending the work of Ref. Hohenadler et al., 2004 to the case of the two- and three dimensional Holstein model with one electron, we have shown that our QMC approach allows accurate calculations with modest computational effort for a large range of parameters. In particular, the minus-sign problem, due to the Lang-Firsov transformation, has been found to diminish quickly with increasing system size, and not to have a significant effect on simulations. We have presented a detailed study of the small-polaron cross over and its dependence on the parameters of the system. In particular, we have focused on finite-size and finite-temperature effects, which have not been investigated systematically before.

As discussed above, our approach is not as fast as other QMC methods for the Holstein polaron,De Raedt and Lagendijk (1982, 1983, 1984); Kornilovitch (1997, 1998, 1999) the main limitation being the restriction to smaller but still reasonably large lattices. This difference in performance is acceptable keeping in mind that it can be extended to the many-electron case (see below), in contrast to the world-line methods of Refs. De Raedt and Lagendijk, 1982, 1983, 1984; Kornilovitch, 1997, 1998, 1999, which work best for the case of a single electron or two electrons of opposite spin, and face a severe sign problem in more than one dimension, similar to other world-line methods.von der Linden (1992) Finally, despite the sign problem, we can perform accurate simulations for virtually all physically interesting values of $`\overline{\omega }`$ and $`\lambda `$, in contrast to the method of Refs. Kornilovitch, 1998 and Kornilovitch, 1999.

Motivated by the promising results of this work and Ref. Hohenadler et al., 2004, our next objective will be a generalization to the many-electron case. To this end, it is important to notice the striking similarity of the QMC method presented here with the determinant method first introduced by Blankenbecler et al.,Blankenbecler et al. (1981) which has been successfully used to study superconductivity and charge-density-wave formation in the Holstein model.Blankenbecler et al. (1981); Scalapino and Sugar (1981); Levine and Su (1990, 1991); Niyaz et al. (1993) In fact, the main difference between our one-electron algorithm and the corresponding grand-canonical method is the calculation of the fermionic weight, which in case of the latter is given by the determinant of $`\widehat{1}+\mathrm{\Omega }`$. Here $`\widehat{1}`$ denotes the unity operator, and $`\mathrm{\Omega }`$ is the corresponding generalization of Eq. (6) to many electrons. Consequently, the numerical effort is expected to be similar to the simulations with one electron. However, an important open question concerns the dependence of the sign problem on the number of electrons. This work is currently in progress.

###### Acknowledgements.

This work has been supported by the Austrian Science Fund (FWF), project No. P15834. One of us (M.H.) is grateful to the Austrian Academy of Sciences for financial support. The majority of the results presented here have been obtained using the Gigabit-Cluster at Graz University of Technology.


format: text

text:
# A class of (ℓ-dependent) potentials with the same number of (ℓ-wave) bound states

## Abstract

We introduce and investigate the class of central potentials

$$V_{\text{CIC}}(g^2,\mu ^2,\mathrm{},R;r)=\frac{g^2}{R^2}\left(\frac{r}{R}\right)^4\mathrm{}\left\{\left[1+\left(\frac{1}{2\mathrm{}+1}\right)\left(\frac{r}{R}\right)^{2\mathrm{}+1}\right]^21+\mu ^2\right\}^2,$$

which possess, in the context of nonrelativistic quantum mechanics, a number of $`\mathrm{}`$-wave bound states given by the ($`\mathrm{}`$-independent !) formula

$$N_{\mathrm{}}^{\text{(CIC)}}(g^2,\mu ^2)=\left\{\left\{\frac{1}{\pi }\sqrt{g^2+\mu ^21}\left(\sqrt{\mu ^21}\right)^1\mathrm{arctan}\left(\sqrt{\mu ^21}\right)\right\}\right\}.$$

Here $`g`$ and $`\mu `$ are two arbitrary real parameters, $`\mathrm{}`$ is the angular momentum quantum number, and the double braces denote of course the integer part. An extension of this class features potentials that possess the same number of $`\mathrm{}`$-wave bound states and behave as $`\left(a/r\right)^2`$ both at the origin ($`r0^+`$) and at infinity ($`r\mathrm{}`$), where $`a`$ is an additional free parameter.

journal: Physics Letters Athanks: FNRS Postdoctoral Researcher

Recently we revisited  the classical problem of establishing upper and lower limits for the number $`N_{\mathrm{}}`$ of $`\mathrm{}`$-wave bound states possessed by a central potential $`V(r)`$ (vanishing as $`r\mathrm{}`$) in the context of nonrelativistic quantum mechanics. As it is well known (see, for instance, ), this number $`N_{\mathrm{}}`$ coincides with the number of zeros, in the interval $`0<r<\mathrm{}`$, of the zero-energy radial wave function $`u_{\mathrm{}}(r)`$, namely of the solution $`u_{\mathrm{}}(r)`$, characterized by the boundary condition $`u_{\mathrm{}}(0)=0`$, of the zero-energy (stationary, radial) Schrödinger equation

$$u_{\mathrm{}}^{\prime \prime }(r)+\left[V(r)+\frac{\mathrm{}\left(\mathrm{}+1\right)}{r^2}\right]u_{\mathrm{}}(r)=0.$$
(1)
\[Here and throughout we use of course the standard quantum-mechanical units such that $`\mathrm{}=2m=1,`$ where $`m`$ is the mass of the particle under consideration, and $`\mathrm{}`$ (a nonnegative integer, $`\mathrm{}=0,1,2,\mathrm{}`$) is the angular momentum quantum number\]. In this context we discovered a class of ($`\mathrm{}`$-dependent) potentials, which (depend on two dimensionless parameters and) have the remarkable property to possess a number $`N_{\mathrm{}}`$ of $`\mathrm{}`$ -wave bound states that (is given by a neat formula and) does not depend on the angular momentum quantum number $`\mathrm{}`$. These potentials, to which we assigned the name $`V_{\text{CIC}}(r)`$ because we found them while visiting the Centro Internacional de Ciencias (CIC) in Cuernavaca (see Acknowledgements), read as follows:

$`V_{\text{CIC}}(g^2,\mu ^2,\mathrm{},R;r)`$ $``$ $`{\displaystyle \frac{g^2}{R^2}}f_{\text{CIC}}(\mu ^2,\mathrm{};{\displaystyle \frac{r}{R}}),`$ (2a)
$`f_{\text{CIC}}(\mu ^2,\mathrm{};x)`$ $`=`$ $`x^4\mathrm{}\left[\left(1+{\displaystyle \frac{x^{2\mathrm{}+1}}{2\mathrm{}+1}}\right)^21+\mu ^2\right]^2.`$ (2b)
Here and throughout $`R,g^2,\mu ^2`$ are three arbitrary positive parameters, the first of which has the dimensions of a length, and the other two are dimensionless and must satisfy the inequality $`g^2+\mu ^2>1`$ (which is necessary for the existence of bound states). The corresponding number $`N_{\mathrm{}}^{\text{(CIC)}}`$ of $`\mathrm{}`$-wave bound states is indeed given by the following ($`\mathrm{}`$-independent !) neat formula (proven below):

$$N_{\mathrm{}}^{\text{(CIC)}}(g^2,\mu ^2)=\left\{\left\{\stackrel{~}{N}_{\mathrm{}}^{\text{(CIC)}}(g^2,\mu ^2)\right\}\right\},$$
(3)
where the double braces denote of course the integer part and

$`\stackrel{~}{N}_{\mathrm{}}^{\text{(CIC)}}(g^2,\mu ^2)={\displaystyle \frac{1}{2\pi }}\sqrt{g^2+\mu ^21}\left(\sqrt{1\mu ^2}\right)^1\mathrm{log}\left[{\displaystyle \frac{1+\sqrt{1\mu ^2}}{1\sqrt{1\mu ^2}}}\right]`$ (4a)
$`\text{if}\mu ^21,`$
$$\stackrel{~}{N}_{\mathrm{}}^{\text{(CIC)}}(g^2,\mu ^2)=\frac{g}{\pi }\text{if}\mu ^2=1,$$
(4b)
$`\stackrel{~}{N}_{\mathrm{}}^{\text{(CIC)}}(g^2,\mu ^2)={\displaystyle \frac{1}{\pi }}\sqrt{g^2+\mu ^21}\left(\sqrt{\mu ^21}\right)^1\mathrm{arctan}\left(\sqrt{\mu ^21}\right)`$ (4c)
$`\text{if}\mu ^21.`$
This remarkable property evokes a mathematical and pedagogical interest, and it is moreover of some applicative relevance because these potentials can be used as comparison potentials: clearly any potential $`V(r)`$ that is more “attractive” than an $`\mathrm{}`$-wave CIC potential, see (2), $`V(r)V_{\text{CIC}}(g^2,\mu ^2,\mathrm{},R;r)`$ for $`0r<\mathrm{},`$ shall possess a number $`N_{\mathrm{}}`$ of $`\mathrm{}`$ -wave bound states at least as large as $`N_{\mathrm{}}^{\text{(CIC)}}(g^2,\mu ^2)`$, $`N_{\mathrm{}}N_{\mathrm{}}^{\text{(CIC)}}(g^2,\mu ^2)`$, and conversely any potential $`V(r)`$ that is less “attractive” than an $`\mathrm{}`$-wave CIC potential, see (2), $`V(r)V_{\text{CIC}}(g^2,\mu ^2,\mathrm{},R;r)`$ for $`0r<\mathrm{},`$ shall possess a number $`N_{\mathrm{}}`$ of $`\mathrm{}`$-wave bound states not larger than $`N_{\mathrm{}}^{\text{(CIC)}}(g^2,\mu ^2)`$, $`N_{\mathrm{}}N_{\mathrm{}}^{\text{(CIC)}}(g^2,\mu ^2)`$ . This motivated us to write the present paper, in order to advertise this finding and to elaborate on it by providing some information on the shape $`f_{\text{CIC}}(\mu ^2,\mathrm{};x)`$ of these potentials, see (2b).

First of all let us note the following qualitative features. Clearly this function is negative, $`f_{\text{CIC}}(\mu ^2,\mathrm{};x)<0,`$ for all positive values of $`x`$, $`0<x<\mathrm{}`$, and it vanishes proportionally to the inverse $`\left[4\left(\mathrm{}+1\right)\right]`$-power of $`x`$ at large $`x`$,

$$\underset{x\mathrm{}}{lim}\left[x^{4(\mathrm{}+1)}f_{\text{CIC}}(\mu ^2,\mathrm{};x)\right]=\left(2\mathrm{}+1\right)^4.$$
(5)
For $`\mathrm{}=0`$, this function has its minimum at $`x=0`$,

$$\underset{0x<\mathrm{}}{\mathrm{min}}\left[f_{\text{CIC}}(\mu ^2,0;x)\right]=f_{\text{CIC}}(\mu ^2,0;0)=\mu ^4,$$
(6)
and it increases monotonically from this minimum value to its vanishing value at $`x=\mathrm{}`$ (namely, $`f_{\text{CIC}}^{}(\mu ^2,0;x)>0`$ for $`0<x<\mathrm{}`$; more specifically, if $`0<\mu ^2<6`$ the second derivative of $`f_{\text{CIC}}(\mu ^2,0;x)`$ is everywhere negative, $`f_{\text{CIC}}^{\prime \prime }(\mu ^2,0;x)<0`$ for $`0<x<\mathrm{}`$, while if $`\mu ^2>6`$, the second derivative changes sign, namely $`f_{\text{CIC}}^{\prime \prime }(\mu ^2,0;x)>0`$ for $`0<x<\stackrel{~}{x}(\mu ^2),f_{\text{CIC}}^{\prime \prime }(\mu ^2,0;x)<0`$ for $`\stackrel{~}{x}(\mu ^2)<x<\mathrm{}`$, with $`\stackrel{~}{x}(\mu ^2)=\sqrt{\frac{\mu ^21}{5}}1`$).

For positive $`\mathrm{}`$, $`\mathrm{}=1,2,\mathrm{}`$, this function vanishes at the origin proportionally to $`x^4\mathrm{}`$ and it has a single minimum at $`x=x_{\text{min}}(\mu ^2,\mathrm{})`$,

$$x_{\text{min}}(\mu ^2,\mathrm{})=\left[\frac{\left(2\mathrm{}+1\right)}{2\left(\mathrm{}+1\right)}\xi (\mu ^2)\right]^{\frac{1}{2\mathrm{}+1}},$$
(7)
where it attains the value

$`\underset{0x<\mathrm{}}{\mathrm{min}}\left[f_{\text{CIC}}(\mu ^2,\mathrm{};x)\right]`$ $`=`$ $`f_{\text{CIC}}[\mu ^2,\mathrm{};x_{\text{min}}(\mu ^2,\mathrm{})]`$ (8)
$`=`$ $`\left(\mathrm{}+1\right)^2\left[{\displaystyle \frac{2(\mathrm{}+1)}{2\mathrm{}+1}}\right]^{\frac{2}{2\mathrm{}+1}}{\displaystyle \frac{\left[\xi (\mu ^2)\right]^{\frac{4\mathrm{}}{2\mathrm{}+1}}}{\left[\xi (\mu ^2)+2\mu ^2(\mathrm{}+1)\right]^2}}.`$
Here we use the short-hand notation $`\xi \left(\mu ^2\right)=\sqrt{1+4\mathrm{}\left(\mathrm{}+1\right)\mu ^2}1`$. Note that $`x_{\text{min}}(\mu ^2,\mathrm{})=1`$, and $`x_{\text{min}}(\mu ^2,\mathrm{})`$ grows as $`\left[\frac{\mathrm{}(2\mathrm{}+1)^2\mu ^2}{\left(\mathrm{}+1\right)}\right]^{\frac{1}{2(2\mathrm{}+1)}}`$ at large $`\mu `$; while the minimum of $`f_{\text{CIC}}(\mu ^2,\mathrm{};x)`$, see (8), decreases as $`\mathrm{}^2\left(1+\mu \right)^2`$ at large $`\mathrm{}`$, and increases as $`L(\mathrm{})\mu ^{\frac{4(\mathrm{}+1)}{2\mathrm{}+1}}`$ at large $`\mu `$ with $`L(\mathrm{})=\left(\mathrm{}\right)^{\frac{2\mathrm{}}{2\mathrm{}+1}}\left(\mathrm{}+1\right)^{\frac{2(\mathrm{}+1)}{2\mathrm{}+1}}`$ $`\left(2\mathrm{}+1\right)^{\frac{2}{2\mathrm{}+1}}`$.

The formulas (7) and (8), as well of course as the expression of the number of bound states, see (4b), simplify considerably for $`\mu ^2=1`$:

$$x_{\text{min}}(1,\mathrm{})=\left[\frac{\mathrm{}\left(2\mathrm{}+1\right)}{\left(\mathrm{}+1\right)}\right]^{\frac{1}{2\mathrm{}+1}},$$
(9)
$$\underset{0x<\mathrm{}}{\mathrm{min}}\left[f_{\text{CIC}}(1,\mathrm{};x)\right]=f_{\text{CIC}}[1,\mathrm{};x_{\text{min}}(1,\mathrm{})]=\mathrm{}^{\frac{4\mathrm{}}{2\mathrm{}+1}}\left(\frac{\mathrm{}+1}{2\mathrm{}+1}\right)^{\frac{4(\mathrm{}+1)}{2\mathrm{}+1}}.$$
(10)
Another remarkable property of the function $`f_{\text{CIC}}(\mu ^2,\mathrm{};x)`$ is displayed by the formula (that can be easily verified by direct integration)

$$_0^{\mathrm{}}𝑑x\sqrt{f_{\text{CIC}}(\mu ^2,\mathrm{};x)}=\frac{\pi g}{\sqrt{g^2+\mu ^21}}\stackrel{~}{N}_{\mathrm{}}^{\text{(CIC)}}(g^2,\mu ^2),$$
(11)
see (4).

Graphs of the function $`f_{\text{CIC}}(\mu ^2,\mathrm{};x)`$ (appropriately renormalized) are presented in Figs. 1 and 2 for $`\mathrm{}=0`$ and $`\mathrm{}=2`$ and for various values of $`\mu `$.

Let us now prove the claim made above, namely the validity of (3) with (4). This is done (albeit without providing here any explanation of how this result was discovered ) by relating as follows the solution $`u_{\mathrm{}}(r)`$ of the zero-energy radial Schrödinger equation (1) (or rather its logarithmic derivative) to a “phase function” $`\eta _{\mathrm{}}(r)`$ that provides a convenient tool to count the zeros of $`u_{\mathrm{}}(r)`$:

$$\sqrt{V(r)}\mathrm{cot}\left[\eta _{\mathrm{}}(r)\right]\frac{V^{}(r)}{4V(r)}=\frac{u_{\mathrm{}}^{}(r)}{u_{\mathrm{}}(r)}.$$
(12)
It is indeed plain that the boundary condition $`u_{\mathrm{}}(0)=0`$ entails $`\eta _{\mathrm{}}(0)=0,`$ while the definition $`\eta _{\mathrm{}}(\mathrm{})=\stackrel{~}{N}_{\mathrm{}}\pi `$ entails that the number $`N_{\mathrm{}}`$ of zeros of $`u_{\mathrm{}}(r)`$ in $`0<r<\mathrm{}`$ is given by the integer part of $`\stackrel{~}{N}_{\mathrm{}}`$. It is moreover easy to verify that the zero-energy radial Schrödinger equation (1) yields for $`\eta _{\mathrm{}}(r)`$ the simple first-order ODE

$$\eta _{\mathrm{}}^{}(r)\left\{1+\frac{\mu ^21}{g^2}\left[\mathrm{sin}\eta _{\mathrm{}}(r)\right]^2\right\}^1=\sqrt{V(r)},$$
(13)
provided the potential $`V(r)`$ satisfies the nonlinear ODE

$$\frac{\mathrm{}\left(\mathrm{}+1\right)}{r^2}+\frac{5}{16}\left(\frac{V^{}(r)}{V(r)}\right)^2\frac{V^{\prime \prime }(r)}{4V(r)}+\frac{\mu ^21}{g^2}V(r)=0.$$
(14)
And it is as well easy to verify that the potential $`V_{\text{CIC}}(g^2,\mu ^2,\mathrm{},R;r)`$, see (2), does indeed satisfy this ODE, (14).

It is then plain, by taking advantage of the boundary values of $`\eta _{\mathrm{}}(r)`$ at $`r=0`$ and at $`r=\mathrm{}`$ (see above) and of the relation (11) with (2), that the integration from $`r=0`$ to $`r=\mathrm{}`$ of the ODE (13) yields the expression (4) for $`\stackrel{~}{N}_{\mathrm{}}^{\text{(CIC)}}`$. Q.E.D.

Finally let us point out a consequence of the obvious identity of the zero-energy Schrödinger equation which obtains from (1) by replacing in it $`V(r)`$ with $`W(r),`$ with the zero-energy Schrödinger equation which obtains from (1) by replacing in it $`\mathrm{}`$ with $`\lambda `$, provided

$$W(r)=\left(\frac{a}{r}\right)^2+V(r),$$
(15)
$$\lambda \lambda (a^2,\mathrm{})=\frac{1}{2}+\sqrt{\left(\mathrm{}+\frac{1}{2}\right)^2+a^2}.$$
(16)
This of course entails that the number $`N_{\mathrm{}}`$ of $`\mathrm{}`$-wave bound states possessed by the potential

$`W_{\text{CIC}}(g^2,\mu ^2,a^2,\mathrm{},R;r)`$ $``$ $`{\displaystyle \frac{1}{R^2}}F_{\text{CIC}}(g^2,\mu ^2,a^2,\mathrm{};{\displaystyle \frac{r}{R}})`$ (17a)
$`F_{\text{CIC}}(g^2,\mu ^2,a^2,\mathrm{};x)`$ $`=`$ $`\left({\displaystyle \frac{a}{x}}\right)^2+g^2f_{\text{CIC}}(\mu ^2,\lambda ;x),`$ (17b)
where $`a`$ is an arbitrary real constant and $`f_{\text{CIC}}(\mu ^2,\lambda ;x)`$ is of course defined by (2b) with (16), is still given by (3) with (4) (this formula provides now the number of bound-states with angular momentum quantum number $`\lambda `$ rather than $`\mathrm{}`$; but this number does not depend on this quantum number, see (4); which of course also entails that $`\lambda `$ is not required here to be an integer). This remark entails that the potential $`W_{\text{CIC}}(g^2,\mu ^2,a^2,\mathrm{},R;r)`$, which depends now on the three dimensionless (arbitrary, positive) constants $`g^2,\mu ^2,a^2`$, provides an additional, more flexible tool to assess, by comparison techniques, the number $`N_{\mathrm{}}`$ of $`\mathrm{}`$-wave bound states possessed by a given central potential $`V(r)`$. Graphs of the function $`F_{\text{CIC}}(g^2,\mu ^2,a^2,\mathrm{};x)`$ are presented in Figs. 3, 4 and 5 for $`g=10,`$ $`\mu =1`$, $`\mathrm{}=0`$, $`\mathrm{}=2`$ and $`\mathrm{}=10`$ and for various values of $`a`$. In Fig. 6, graphs of the function $`F_{\text{CIC}}(g^2,\mu ^2,a^2,\mathrm{};x)`$ are presented for $`a=1`$, $`\mathrm{}=5`$ and for various values of $`g`$ and $`\mu `$ such as the number of $`\mathrm{}`$-wave bound states is always equal to three, see (4).

Let us end this paper by emphasizing that, because of the availability of two, or even three, free parameters (for every given value of the angular momentum quantum number $`\mathrm{}`$), the family of CIC potentials considered herein is quite flexible (as also shown by the figures); these makes these potentials – which have a known number of bound states – quite convenient to be used as comparison potentials, whenever one wishes to assess quickly and easily the number of ($`\mathrm{}`$-wave) bound states associated with a given potential.

Acknowledgments

The results reported in this paper have been obtained during the Scientific Gathering on Integrable Systems held from November 3rd to December 13th, 2002, at the Centro Internacional de Ciencias (CIC) in Cuernavaca. It is a pleasure to thank professor Thomas Seligman, director of CIC, and his collaborators and colleagues for the pleasant hospitality and the fruitful working environment provided by CIC.


text:
# 𝐵_𝑐 photoproduction at HERA collider11footnote 1The short talk given at Weizmann Institute Workshop on HEAVY QUARK PHYSICS AT HERA II (Rehovot, Israel, October 19-22, 2003).

## 1 Calculation technique

The main features of the $`ep`$-production of $`B_c`$in the frame work of pQCD had been discussed in our previous work. The dominant subprocess $`\gamma gB_c+X`$ is described by 24 Feynman diagrams. We suppose that the amplitude of the production process can be expressed through the amplitude $`T`$ of hard production of four heavy quarks $`b\overline{b}c\overline{c}`$ followed by the soft fusion of $`\overline{b}`$ and $`c`$ quarks into $`B_c`$ meson due to the meson wave function $`\mathrm{\Psi }`$:

$$A_{\gamma gB_c+b+\overline{c}}T_{\gamma gb\overline{b}c\overline{c}}^{\text{hard}}(p_i,\stackrel{}{q})\stackrel{~}{\mathrm{\Psi }}_q^{}(\stackrel{}{q})d^3q,$$

where $`p_i`$ are four momenta of $`B_c`$, $`b`$ and $`\overline{c}`$; $`\stackrel{}{q}`$ is the three momentum of $`\overline{b}`$ quark in the $`B_c`$ meson rest frame.

One can neglect the $`\stackrel{}{q}`$ dependence in $`T_{\gamma gb\overline{b}c\overline{c}}^{\text{hard}}`$ and get for $`S`$ wave states production the following equation:

$$A_{\gamma gB_c+b+\overline{c}}\mathrm{\Psi }(0)T_{\gamma gb\overline{b}c\overline{c}}^{\text{hard}}(p_i,0),$$

where $`\mathrm{\Psi }(0)`$ is the value of $`B_c`$ wave function value at origin (coordinate space).

In such approach $`\overline{b}`$ and $`c`$ quarks in $`B_c`$ have the same velocities. Thus:

$$p_{\overline{b}}=\frac{m_b}{M_{B_c}}P_{B_c},$$

$$p_c=\frac{m_c}{M_{B_c}}P_{B_c}.$$

If one need to get the yields of $`B_c`$ and $`B_c^{}`$ separately, the product of spinors $`v_{\overline{b}}\overline{u}_c`$, corresponding to the $`\overline{b}`$ and $`c`$ quarks in the $`T_{\gamma gb\overline{b}c\overline{c}}^{\text{hard}}`$ amplitude, should be substituted by the projection operator

$$𝒫(\mathrm{\Gamma })=\sqrt{M}\left(\frac{\frac{m_b}{M}\widehat{P}_{B_c}m_b}{2m_b}\right)\mathrm{\Gamma }\left(\frac{\frac{m_c}{M}\widehat{P}_{B_c}+m_c}{2m_c}\right),$$

where $`\mathrm{\Gamma }=\gamma ^5`$ for $`s=0`$, or $`\mathrm{\Gamma }=\widehat{\epsilon }^{}(P_{B_c},s_z)`$ for $`s=1`$, where $`\epsilon (P_{B_c},s_z)`$ is the polarization vector for the spin-triplet state.

Nevertheless, one can easily show that $`\mathrm{\Gamma }`$ can be expressed through the $`v_{\overline{b}}`$ and $`\overline{u}_c`$:

$$𝒫(\gamma ^5)=\sqrt{\frac{2M}{2m_b2m_c}}\frac{1}{\sqrt{2}}\{v_b(p_b,+)\overline{u}_c(p_c,+)v_b(p_b,)\overline{u}_c(p_c,)\},$$

$$𝒫(\widehat{\epsilon }^{}(P_{B_c},1))=\sqrt{\frac{2M_{B_c}}{2m_b2m_c}}v_b(p_b,+)\overline{u}_c(p_c,),$$

$$𝒫(\widehat{\epsilon }^{}(P_{B_c},0))=\sqrt{\frac{2M_{B_c}}{2m_b2m_c}}\frac{1}{\sqrt{2}}\{v_b(p_b,+)\overline{u}_c^{}(p_c,+)+v_b(p_b,)\overline{u}_c(p_c,)\},$$

$$𝒫(\widehat{\epsilon }^{}(P_{B_c},+1))=\sqrt{\frac{2M_{B_c}}{2m_b2m_c}}v_b(p_b,)\overline{u}_c(p_c,+).$$

The more detailed description of the calculation technique one can find in the papers , where the gluonic and photonic $`B_c`$ production cross sections had been calculated (see also ).

To get the cross section of $`B_c`$ production in $`ep`$ collision one should covert the cross section of the subprocess $`\gamma gB_c+X`$ with the distribution function of gluon in the initial proton (we use CTEQ6 parametrization of PDF ) and Weizsäcker-Williams photon function.

## 2 Results

The parameter values have been chosen as follows:

$$\alpha _s=0.2,$$

$$m_c=1.5\text{GeV},$$

$$m_b=4.8\text{GeV}$$

$$|\mathrm{\Psi }(0)|^2=0.116\text{GeV}^3.$$

The total production cross section of $`B_c^+`$, $`B_c^+`$, $`B_c^{}`$ and $`B_c^{}`$ mesons is about 5 pb. Higher excitations yield about 30% to the total cross section. Thus for integrated HERA luminosity 1 $`\text{fb}^1`$ about $`710^3`$ $`B_c`$ mesons will be produced. It is rather small amount of $`B_c`$ to find these particles at HERA, because of small branching ratios of $`B_c`$ decay mode, which are interesting for the extraction of $`B_c`$ signal. Nevertheless, one can hope to find these particles at HERA collider.

## 3 Remarks

The most popular model of $`B`$ meson production is the fragmentation one. In accordance to this mechanism at large $`p_T>p_T^\text{f}`$ quark $`b`$ quark is produced in the hard process followed by the soft process of heavy quark fragmentation. For example: $`\gamma gb\overline{b}`$ followed by $`bB`$. The meson production cross section in the frame work of fragmentational model have a very simple form:

$$\frac{d\sigma (\gamma gB+X)}{dz}=\sigma (\overline{b}b)D_{bB}(z),$$

where $`z=2E_{B_c}/\sqrt{s_{\gamma g}}`$.

But can we estimate the numerical value "large" $`p_T^\text{f}`$?

For $`B_c`$ production we can do it in the frame of model discussed in Chapter 1 and 2 of this paper:

$$p_T^\text{f}30\text{GeV}5M_{B_c}$$

At $`p_T<5M_{B_c}`$ nonframentational terms (recombination terms) give the dominant contribution into $`B_c`$ production cross section.

## Acknowledgements

We would like to thank L. Gladilin for the fruitful discussion. The work was supported in part by grants of RF Education Ministry E00-3.3-62, CRDF M0-011-0, Russian Scientific School 1303.2003.