The Effect of Next-Nearest Neighbour Hopping in the One, Two, and Three Dimensional Holstein Model

Allowing a single electron to hop to next-nearest neighbours (NNN) in addition to the closest atomic sites in the Holstein model, a modified Trugman method is applied to exactly calculate the effect on the polaronic effective mass in one, two, and three dimensions, building on the previous study of the one-dimensional NNN Holstein model. We also present perturbative calculations and a heuristic scaling factor for the coupling strength and ion frequency to nearly map the NNN Holstein model back onto the original Holstein model. When account is taken of the modified electronic bandwidth near the electron energy, we find that including NNN hopping effectively increases the polaron effective mass.

In the realm of BCS theory, it is well known that electron-phonon interactions in solid materials are integral to the emergence of superconductivity, as they are responsible for the effective attraction that leads to the formation of Cooper pairs 1 . On the other hand, the importance of electron-phonon interactions in high temperature superconductivity is not yet clear 2 . Since polarons are simply quasiparticles consisting of electrons dressed with the net effect of these electron-phonon interactions, it is important to understand this basic building block to fully understand conventional superconductivity, and possible extensions to nonconventional superconductors 3 . To this end, the problem of a single electron in the conduction band of a crystal lattice has been extensively studied 4 . Specifically, a numerically exact algorithm for solving the Holstein model with tight-binding electron bands in the thermodynamic limit was formulated in ref. 5, and now that problem is effectively solved. Several extensions were subsequently reported, including ones to better manage disparate electron (t) and phonon (ω E ) energy scales (in particular, ω E ≪ t) 6,7 , higher dimensionality [8][9][10] , extended interaction range 11,12 and inclusion of next-nearest neighbour (NNN) single-particle hopping amplitude 13 . In this last study it was found that including NNN hopping in the one dimensional Holstein model altered significantly the electron's effective mass in strong coupling.
The purpose of this paper is to follow up on this study. Thus far studies of polaron properties within the Holstein model have revealed that the effective mass becomes very large with rather modest electron-phonon coupling strength. This is incompatible with experiment, specifically with the evidence that some conventional superconductors have a large electron-phonon coupling strength, and yet show almost no sign of single-electron polaronic behaviour in the normal state 14 . Chakraborty et al. 13 found that including NNN hopping in the one-dimensional Holstein model could decrease the polaron effective mass significantly, particularly at strong coupling. This is potentially very important since this is a means for lowering the polaron effective mass to a realistic level, such that an Eliashberg treatment 15-19 makes sense.
To more fully understand the effects due to NNN electron hopping we will first present our perturbation theory calculations for the energy and effective mass of the NNN Holstein model in one, two, and three dimensions. We use square and simple cubic lattices for two and three dimensions, respectively. These results agree for sufficiently low coupling strength with our exact numerical calculations using our previously refined algorithm for the Holstein model 12 extended to include NNN interactions. We note that quantitative agreement with perturbation theory extends over a surprisingly limited range of electron-phonon interaction strength, even in three dimensions, which is the most applicable to bulk superconductors. A low phonon frequency approximation to the perturbation theory results suggests a scaling of the phonon frequency with the low energy effective bandwidth, which explains the results obtained as a function of NNN electron hopping. We also note an additional scaling Scientific RepoRts | 6:32591 | DOI: 10.1038/srep32591 factor that accounts for the results with non-zero NNN hopping with respect to those with nearest-neighbour hopping only, over a more extended coupling strength range.
Since including NNN electron hopping also modifies the 'effective' electronic bandwidth (to be defined more precisely below), we should account for this in using the appropriate phonon frequency. That is, since altering the adiabatic ratio ω E /t, even in the case with NN hopping only, is known to lead to changes in the polaronic effective mass for the same coupling strength, then we should be careful to use an appropriately scaled phonon frequency.
After a brief introduction we use perturbation theory to determine the polaron effective mass in weak coupling. Since these expressions are analytical, they are well-suited to examine the various scaling factors. We then present exact solutions, in one, two, and three dimensions, to examine the effect on polaron mass over the entire coupling range. We also note a heuristic scaling, found numerically, that very accurately maps the parameters with NNN hopping back to those without, before closing with a summary.

Model and Methods
Holstein model. The Holstein model 20 is perhaps the simplest model for describing electron-phonon interactions; it treats (optical) phonons as local ion vibrations, and assumes that each atomic site oscillates with the same characteristic frequency ω E . With NNN hopping included, the Hamiltonian that describes such a system is: 3 is a dimensionless measure of the electron-phonon coupling strength, with M being the atomic mass and α being the coupling strength as defined in real space.
In order to diagonalize this Hamiltonian, we transform into k-space, according to the equation: where R j points to lattice site j, and  k is a wave vector summed over the First Brillouin Zone (FBZ). The relation for ˆ † c j may be obtained simply by taking the Hermitian conjugate of the above expression, and the bosonic Fourier transforms are defined almost identically. In the FBZ there are N distinct k values within π π − a a ( / , / ) in each direction. The transformed Hamiltonian then becomes and holds for all dimensions with dispersion relations   k ( ) given in Table 1.
In this paper, we will examine various properties of the ground state, which for > − ′ t t 1 4 is at zero total crystal momentum p. For all dimensions this results in a low energy dispersion E p quadratic in p = |p|, so that the ground state effective mass of the electron is given by: Weak-coupling perturbation theory. Beginning with the perturbative approach (in the weak coupling regime), we consider the electron-phonon interaction to be the perturbation, so that the unperturbed energy is The unperturbed ground state for arbitrary p is therefore Here, 0 is simply the electron-phonon vacuum state. Unperturbed excited states include all states with a single electron and any number of phonons such that the total crystal momentum still adds to p. It is easy to check that given these definitions, , independent of the choice of total crystal momentum. Under these conditions, the energy correction to second order (in α or g) is: In this case, V pert is the electron-phonon interaction term of the Hamiltonian in Eq. 3. Note that only one phonon is considered in the unperturbed excited states because V pert only creates (annihilates) one phonon. Upon evaluating this sum, we may apply Eq. 4 to find the effective mass. In order to do so, we convert the sum in Eq. 6 to an integral over k-values, since the thermodynamic limit (N → ∞ ) implies a continuum of k values between − π/a and π/a.
In the one-dimensional case, the corrected energy according to second order perturbation theory (in g) is and β ≡ t t / 2 . The above result, substituted into Eq. (4), gives an expression for the effective ground-state electron mass m* at p = 0: , where W is the electronic bandwidth. This definition for λ applies for 3D as well, though in 2D, λ ≡ More generally, evaluating for the second-order energy correction in two or three dimensions proves tedious, and has a cumbersome, unenlightening answer (as in eq. (8)). For these cases, we have also integrated the result numerically to check our analytical results. In the figures that follow, these are referred to as "numerically integrated perturbation theory. " For an approximate analytical result, we observe that the integrand in the energy correction expression decays more or less to 0 by some k 0 < π for most small values of the parameters β and ω E . Cutting the integral off at this k 0 and making the approximation that  k 1 for k ≤ k 0 , we achieve the analytic approximations shown in Table 2. Unfortunately, these approximations prove to be rather crude in 2D and 3D, which limits their usefulness. Regardless, we present them alongside our numerically integrated results for completeness in Figs 1, 2 and 3. The (1 4 ) 3/2 2 8 Table 2. Approximate ground state effective mass for small ω ω ≡ t / E E and small β ≡ t t / 2 , from weakcoupling perturbation theory.

Figure 1. Numerically integrated perturbation theory, approximate perturbation theory, and exact numerical solution in 1D.
Note that the approximate perturbation theory is on top of the numerically integrated perturbation theory so here the approximation in Table 2 works very well. We used parameter values of ω E /t = 0.1 and t 2 /t = 0.025.

Figure 2. Numerically integrated perturbation theory, approximate perturbation theory, and exact numerical solution in 2D.
Here the approximate perturbation calculation fails quite badly even in the very small perturbative regime from g 2 = 0 to g 2 = 1. We used parameter values of ω E /t = 0.2 and t 2 /t = 0.025. Scientific RepoRts | 6:32591 | DOI: 10.1038/srep32591 approximations do work better for smaller ω E ; however we have tested then in a reasonably physically representative regime of small ω E 12 and even here the agreement is poor. Secondly, it is important to note that even without approximations, the range of coupling strengths over which the numerical perturbation theory is valid is very small. This feature is similar to our results for the standard Holstein model 12 so we do not recommend using the perturbation calculations for physical predictions but only as a check on more powerful numerical calculations such as the Trugman Method.

Modified Trugman method for exact numerical solutions. The single polaron problem is solved
here with the variational exact diagonalization method described in Bonča et al. 5 and revised by the authors as described in refs 6,12 to account for a rapidly growing Hilbert space from the additional terms in the Hamiltonian. For the data included in these plots 20-100 preliminary diagonalizations were performed with the most strongly contributing basis states selected at each iteration to seed the next iteration. All results were converged for the effective mass. While our algorithm for generating basis states differs significantly from that of Chakraborty et al. 13 , our 1D results are in agreement with theirs and our algorithm permits efficient calculations in higher dimensions and for realistic values of ω E /t.
Note that the introduction of non-zero t 2 leaves the bandwidth invariant; however it does alter the curvature of the dispersion relation near k = 0. The result is a different bare electron mass (see Eq. (9)) and, correspondingly, a different effective bandwidth. The different bare electron mass and different normalization of the effective mass has been noted before 13 , but the different effective bandwidth was not considered. If we take some small  k a 1/ 0 , and calculate the bandwidth in the region [− k 0 , k 0 ] up to second order accuracy in k 0 , then we find the ratio of next-nearest neighbour bandwidth to nearest neighbour bandwidth to be: and: This suggests the phonon energy scale ω E also should be rescaled by the same factors to keep the ratio of phonon energy to effective bandwidth constant. The interaction strength parameter g is dimensionless and thus remains unchanged, though the electron-phonon interaction term in the Hamiltonian is rescaled since it is proportional not simply to g, but to ω g E . It can be seen from the result in Table 2 that rescaling ω E by the effective bandwidth change would transform the NNN approximate effective mass onto that for NN hopping only. In other words, if we use a renormalized phonon frequency with the same value with respect to the effective electronic bandwidth, then the addition of next nearest neighbour hopping has no effect on the effective mass (according to Table 2). However these are only approximate perturbation calculations and the exact results show that the NNN effective mass is substantially different from the NN effective mass even when the proper scalings have been taken into account.
On the other hand, this rescaling of ω E (which was not done by Chakraborty et al. 13 ), definitely reduces the effect of the NNN hopping on the polaron effective mass. In Fig. 4 we compare the data with and without the scaling of ω E to the original Holstein model in 1D with = ± . t t / 01 2 . The 1D case has no crossover coupling strength in the standard Holstein model, where the polaron effective mass suddenly begins to increase exponentially with coupling strength, and while t 2 > 0 did decrease the effective mass somewhat, no crossover point was found for the NNN Holstein model in 1D either. Our results are in agreement with Chakraborty et al. 13 for the same value (i.e. unscaled) of ω E . However, including the correction described in the preceding paragraph decreases the effective mass change and even the direction of the change. The 1D case is not very realistic for bulk materials so we continue with 2D and 3D calculations.
In two and three dimensions we find that for small NNN hopping parameters the effective mass deviates from the standard Holstein model only slightly until the crossover coupling strength in 2D and 3D is reached. At this point the effective masses increases sharply for both models. However, the introduction of non-zero t 2 changes the crossover point slightly. We have included the unscaled results with the ω E /t scaled results for the sake of completeness in Figs 5 and 6 for 2D and 3D, respectively.
While our approximate perturbation theory suggests that the scalings of ω E /t given in Table 2 would map the NNN effective masses onto the NN effective mass, Figs 4, 5 and 6 make it clear that this scaling does not work very well except in the 2D case. The approximate perturbation theory only agrees exactly with the numerical  Note that the scaling of ω E maps the NNN calculations back to the original. This is in agreement with our approximate perturbation theory results, though perhaps serendipitously since the approximate perturbation theory was not very close to the exact numerical perturbation theory. perturbation theory in 1D and there is no reason to trust even the numerical perturbation theory for the effective mass outside of very small coupling strength as we have remarked on other occasions 12 . In 2D the approximate ω E scaling suggested in Table 2 mapped the effective masses back close to that obtained with no NNN hopping, but in 1D and 3D it over-corrected the change in the effective mass. Also, not surprisingly, in all three dimensions a simple scaling does not work well in the very strong coupling regime. So in agreement with Chakraborty et al. 13 , NNN hopping does introduce changes in the properties of the polaron, and leads to a decreased effective mass for some additional (positive) t 2 hopping if no phonon frequency scaling is introduced. But NNN hopping leads to an increased effective mass if the phonon frequency is also increased to account for the increase in the 'local' bandwidth. Most importantly for our understanding of the conventional framework for superconductivity, the inclusion of NNN hopping changes the critical coupling strength in 3D at which the effective mass increases sharply towards infinity (see the large coupling regime of Fig. 6).
Heuristic scaling. In the course of our investigations we further found a heuristic scaling of the coupling strength that, combined with the bandwidth-inspired scaling, works very well; however, the underlying physical motivation is still rather unclear. We introduce a scaling factor in the dimensionless interaction parameter ⋅ = g B g original scaled such that the new term in the Hamiltonian is:

Summary and Conclusions
We have performed weak coupling perturbation calculations for the Holstein model with next nearest neighbour (NNN) hopping and compared them with exact calculations in one, two, and three dimensions. We confirmed previous results obtained in one dimension by Chakraborty et al. 13 . However, we point out that a more appropriate comparison, at least in weak coupling, requires a change of the phonon frequency, ω E , as given in Table 2, for the case with NNN hopping. With this change accounted for we find that the effect of NNN hopping on the effective mass is opposite to the change without a phonon frequency change. Including t 2 with the same sign as t reduces the polaron effective mass when the same phonon frequency is used, whereas including a change in the phonon frequency according to the changes in effective bandwidth as in Eqs (10) and (11) increases the effective mass in 1D and 3D. We feel that including the change in phonon frequency is the more physically correct procedure. In two dimensions NNN hopping has very little effect on the effective mass in weak coupling. In 1D and 3D we have also found a heuristic scaling factor for g, the dimensionless electron-phonon coupling strength, that maps the results for the polaron effective mass of the NNN Holstein model back onto the standard model without NNN hopping. While the physical reason for this effect is not known, this heuristic scaling allows us to crudely estimate how increasing t 2 impacts the coupling strength at which the (sharp) crossover occurs for polaronic behaviour. For small values of t 2 /t the crossover remains in the regime of moderate electron-phonon coupling. Therefore it remains difficult to reconcile the fairly strong coupling attributed to some real metals/ superconductors with the diverging effective mass predicted for a single polaron at the same coupling strength.