Role of long-range coupling on the properties of single polarons in models with dual electron-phonon couplings

We begin with a brief review of the extended Holstein model (EHM) $g(q)$ coupling both to establish a link to existing previous work,Feshke ; M1 and also because a similar model was used in Refs. ZX, ; Yao, for quantitative modeling of el-ph coupling in 1D CuO chains.

The extended Holstein coupling has the standard form Feshke ; M1 :

where $g$ is a dimensionless parameter quantifying the overall strength of the EH coupling, while the spatial dependence is described by:

The Heaviside function allows us to control the range of this model depending on the cutoff $M$. If $M$ is chosen so that the coupling is finite only for $j=0$, we regain the standard Holstein model (on-site coupling). Increasing $M$ allows us to study increasingly longer-range coupling, to see how this influences the results. In the following, we characterize the coupling range with (i) the label $1$ for the Holstein model ($M=0.5$); (ii) the label 3 when the extended coupling includes the three sites $j=0,\pm 1$ ($M=1.5$); (iii) the label 5 when the extended coupling includes the five sites $j=0,\pm 1,\pm 2$ ($M=2.5$); and (iv) the label 7 when the extended coupling includes the seven sites $j=0,\pm 1,\pm 2,\pm 3$ ($M=3.5$).

To make meaningful comparisons for these different ranges, we need to rescale the magnitude $g$ as a function of the cutoff $M$, so that the physical energy scale – the polaron formation energy in the single-site limit – remains the same. For EHM, the polaron formation energy (at $t=0$) for a given cutoff $M$ is

Comparing this energy to the ground-state energy $-2t$ of the bare electron defines the dimensionless electron-phonon coupling $\lambda=-{\epsilon_{p}}/{2t}$. For the Holstein model, then, $\lambda_{h}=g^{2}/(2t\Omega)$. In order to have the same overall $\lambda_{eh}=\lambda_{h}$ for different coupling ranges in the EHM, we rescale the value of $g$ such that

The $g(k,q)$ Peierls coupling is due to linear modulations of the nn hopping because of displacements of the atoms involved in the hopping:

Following previous work,Dominic we use the dimensionless effective coupling $\lambda_{p}=2{\alpha}^{2}/(\Omega t)$ to characterize the strength of the Peierls coupling.

For the Peierls coupling we do not consider a longer-range extension. That can be accomodated by VED, however it is not very physical given that such terms would arise from modulations of the longer-range hopping integrals. The hopping integrals are assumed to decrease exponentially with the distance between the sites, which is a much faster decrease than the power-law screening of the Coulomb interactions responsible for the $g(q)$ couplings. This is why we believe that the primary source of long-range coupling must be coming from the $g(q)$ model, and we only consider this option in our work.

The first dual model that we study combines the extended Holstein+Peierls couplings. Its Hamiltonian is:

and the strenghts of the couplings will be characterized by a $\lambda_{p}$ and a $\lambda_{eh}$, while the EHM range is set by $M$ as discussed above. We note that this extends the work of Ref. DominicPRB, where the combination of Peierls and Holstein models was investigated, providing another point of contact to previous work, and another way to validate our results. The comparison with those results will illustrate the effects of longer-range EHM coupling on the properties of the resulting polaron. On the other hand, contrasting these results with those of the EHM coupling will also allow us to understand the effects of adding Peierls coupling in the model.

However, although the EHM model was already suggested to be relevant for 1D chains while the Peierls coupling is generically expected to appear in any system with el-ph coupling, their combination is not very physical. This is because Holstein-like couplings are due to an internal distortion of the ’polar molecule’ assumed to be located at each lattice site, when an additional electron visits it. As such, this coupling is to the ’vibron’ describing the internal distortion of the molecule, whereas the Peierls coupling is to actual phonons describing the displacements of sites from their equilibrium position.

To mitigate for this issue, we propose another $g(q)$ model which does couple to the same actual phonons like the Peierls model, so that their combination is more physical. For completeness, we also note that for a chain with a strictly one-site basis, there is a single longitudinal acoustic phonon mode, there is no optical mode that could be modelled as an Einstein phonon, as we do here. Nevertheless, we continue to use the simplified Einstein phonon mode instead of a more realistic acoustic phonon mode. Physically, this is because both the Peierls coupling and the extended breathing-mode coupling described below vanish for $q=0$, i.e. there is no coupling to the gapless phonons at the center of the Brillouin zone. In contrast, the coupling is strong to the phonons at the Brillouin-zone boundary at $q=\pi$, which describe anti-phase motion of consecutive ions. In other words, the phonons to whom the coupling is strong are similar to optical phonons (gapped and describing anti-phase motion). The second reason to use an Einstein model is because this is what was used in all the previous work we contrast our results with. By using the same phonon model we insure that any differences in results are due to the changes in the coupling.

In this model, the electron-phonon coupling is due to the modulation of the screened Coulomb interaction $V(i-j)$ between the electron at site $i$ and the ion at site $j$, because of the motion of both sites. To linear order, this leads to the extended breathing-mode (EBM) coupling:

where $g_{j-i}=-g_{i-j}$ is an odd function. As a result, the total coupling to the phonons at the carrier site $j=i$ vanishes, and the EBM coupling can be written in a form similar to the EHM coupling:

however now $g(j)$ is an odd function (as opposed to an even one, for the EHM), so $g_{0}=0$. So while in the EHM the on-site coupling is the largest, for the EBM the on-site coupling vanishes.

To facilitate comparisons between results of the EBM and the EHM, in the following we take:

so that the long-range decay is similar. The coefficients $\chi_{j}$ are just signs. Specifically, we study three possible choices: (i) $\chi_{j}=j/|j|$, describing a case where there is overall attraction between the electron and the other ions in the lattice (we always take $g_{bm}>0$). We call this the ’EBM+’ case; (ii) $\chi_{j}=-j/|j|$, describing a case where there is overall repulsion between the electron and the other ions in the lattice. We call this the ’EBM-’ case; and (iii) $\chi_{j}=1$. This last choice is not leading to an odd $g_{j}$ function so it is not a proper breathing-mode coupling, instead it describes the EHM coupling without the on-site term. We will study it for comparison purposes, and we label it (somewhat improperly) ’EBM0’. Like for the EHM, we use the labels 1, 3, 5 or 7 to characterize the range of the extended coupling.

To characterize the strength of the EBM coupling, we proceed as follows. First, we note that we could use a formula similar to the one used for the $EBM$ coupling. However, that choice would mean that if $\lambda_{eh}=\lambda_{ebm}$, then $g_{bm}\gg g$ in order to compensate for the absence of the on-site coupling. This is not suitable for our purpose, which is to compare models with equal magnitudes of the longer range coupling, to see how important are the details of the chosen form.

This is why we will characterize $g_{bm}=g$ by saying (rather improperly) that $\lambda_{eh}=\lambda_{ebm}$. In other words, the magnitude $g_{bm}$ for a given value of $\lambda_{ebm}$ is taken to be equal to the magnitude $g$ of the EHM of equal range, when $\lambda_{eh}=\lambda_{ebm}$. Comparing EHM with EBM0 will allow us to gauge the importance of the on-site coupling (present for EHM and absent for EBM0) for equally strong longer-range coupling, while comparing EBM0 with EBM$\pm$ will reveal the importance of the even vs. odd coupling under inversion.

The second, more physical dual model that we investigate therefore combines the EBM+P couplings:

As mentioned just above, there are actually 3 distinct variants of the EBM depending on the choice of $\chi_{j}$. All three are investigated below.

First, we analyze the properties of polarons in the dual coupling model combining the extended Holstein + Peierls couplings. We set the energy scale to be $t=1$. We generated results for $\Omega=0.5,1$ and $3$, covering the crossover from weakly adiabatic to weakly anti-adiabatic regime. Figures 2 and 3 show representative results for a variety of coupling strengths and ranges, for the largest and smallest $\Omega$ values considered.

In all of these panels, the black dotted curves (label 1) are for the (on-site) Holstein+Peierls dual models, and are in great agreement with those generated in Ref. DominicPRB, using Bold Diagramatic Monte Carlo and the variational Momentum Average approximation. In particular, they show the expected change in the dispersion from one with a ground-state (GS) at $k_{gs}=0$ at weak Peierls coupling into one with a doubly degenerate GS with a $\pm k_{gs}\neq 0$ at strong Peierls coupling. The reason for this change is the dynamical generation of 2nd nearest-neighbor, phonon-mediated, hoppings in the presence of Peierls coupling. Dominic

As the range of the extended Holstein coupling is increased (labels 3, 5, 7 for couplings including nn, 2nd nn, and 3rd nn neighbors respectively) such that $\lambda_{eh}=0.5$ remains fixed, Fig. 2 shows a non-trivial quantitative change with the addition of the extended coupling to the nn sites, but then results are converged and no longer change if the range is further increased. By contrast, for the smaller $\Omega=0.5$ in Fig. 3, we see that dispersion continues to change as the range is increased. For weak Peierls coupling $\lambda_{P}\lesssim 0.5$, convergence is achieved at the 2nd nn range, but for the stronger Peierls coupling the dispersion is not yet converged at 3rd nn range. Furthermore, the quantitative changes are much more substantial than was the case for the results of Fig. 2.

At first sight, this strong sensitivity to the range of the extended coupling is rather unexpected, considering how relatively weak is the coupling to the 2nd and especially the 3rd nn sites, compared to the on-site one (see Fig. 1). However, it is known that the spatial size of the polaron cloud increases as $\Omega$ decreases, especially for moderate Holstein coupling such as the one used here, $\lambda_{h}=0.5$. One expects the structure of a more extended cloud to be quite sensitive to the details of the extended coupling, as the probabilities for various configurations could be reshuffled significantly for an extended vs. a local coupling. This appears to be consistent with the results reported so far.

However, the interplay between the effects of the extended range and those of a dual coupling are, in fact, even more intricate, as can be seen from Fig. 4. Here we consider how the GS properties of the polaron, specifically its energy $E(k_{gs})$, GS momentum $k_{gs}$, quasiparticle weight $Z_{gs}$ and inverse effective mass change with increasing Peierls coupling $\lambda_{P}$, at a fixed but large $\lambda_{h}=2$ and for a large $\Omega=3$. For these values the polaron is closer to the small polaron regime (note that $Z_{gs}<0.5$ at $\lambda_{P}=0$ for the Holstein model, label 1) and one would therefore anticipate a reduced sensitivity to range. The results show a significant change when coupling to the nn sites is turned on (label 3). This is reasonable, given that the Peierls coupling also involves phonons on the nn sites but with different signs, so ’interference’ effects can be expected (we further comment on this below). The less expected fact is that for larger $\lambda_{p}$ there is another sizable change when turning on the coupling to the 3rd nn sites (label 5), as most clearly seen in the value of the effective polaron mass. We find this hard to explain in any simplistic way.

Figure 5 provides a different way to analyze the sensitivity to the range of the coupling. Here we plot the critical value $\lambda^{*}_{P}$ at which $k_{gs}$ has the transition from 0 to a finite value, Dominic for several EHM strenghts $\lambda_{eh}=\lambda_{h}$ (curves of different colors in each panel). The various panels correspond to various ranges of the EHM coupling.

The black dotted curves are for Peierls only coupling ($\lambda_{h}=0$), thus they are the same in all panels. They agree with the known result from Ref. Dominic, and are provided as guides to the eye. The other three curves in each panel are for increasing values of the EHM coupling $\lambda_{h}$. For the dual Holstein+Peierls coupling (panel labelled 1), an increase of $\lambda_{h}$ results in a decrease of $\lambda^{*}_{P}$, in agreement with the results reported in Ref. DominicPRB, . Extending the Holstein coupling to nn sites (panel labelled 3) has a drastic effect: now $\lambda^{*}_{p}$ increases strongly with increasing $\lambda_{h}$ nearly everywhere in the parameter space. (The exception is for the smaller value $\lambda_{h}=0.5$ for a narrow range of phonon frequencies close to $\Omega\approx 0.5$, where $\lambda^{*}_{p}$ is slightly below the value for the pure Peierls model.) Further extending the range to 2nd and to 3rd nn sites (panels labelled 5 and 7, respectively) has less severe impact. For small $\Omega$ there is a visible decrease of $\lambda^{*}_{P}$ with increasing range for the strongest $\lambda_{h}=2$, but everywhere else the results for $\lambda^{*}_{P}$ appear to be converged – at least at this scale. This supports the view that the strongest sensitivity to the range of the coupling is observed as one moves towards the adiabatic limit. However, we emphasize again that even for parameters like those used in Fig. 4, which show little difference in the value of $\lambda^{*}_{P}$ between 2nd and 3rd nn ranges, the changes in $m^{*}$ are much more substantial.

Our conclusion is that insofar as a dual EHM+P coupling is concerned, the polaron properties’ sensitivity to the specific modeling of the extended coupling should certainly be a big concern in the adiabatic limit. As discussed, this is in line with expectations of a more extended polaron cloud in this limit. However, our results show that the sensitivity to the range of the coupling could be significant even rather far from the adiabatic limit, and that it varies for different properties. Furthermore, even in parts of the parameter space where the results converge fast with increasing range, we find that adding nn coupling makes a significant difference and the results are very different from those due to on-site coupling only.

The significant differences observed everywhere when adding coupling to nn sites motivates us to study how sensitive the results are to the particular $g(q)$ coupling used. To do this, in this section we replace the EHM with EBM coupling. We remind the reader that EBM has zero on-site strength, so the range here starts from coupling to nn sites and can be extended to coupling to 2nd and 3rd nn sites (labels 3, 5 and 7, respectively). Also, we study three non-equivalent possible overall signs for these couplings, as illustrated in Fig. 1.

In Figures 6 and 7 we show the shifted polaron dispersions energies $E_{k}-E_{0}$ vs. $k$ at $\Omega=3$ and $\Omega=0.5$, respectively. We remind the reader that $\lambda_{ebm}=0.5$ means that the coupling $g_{bm}$ is chosen to be equal to the EHM $g$ corresponding to the same $\lambda_{eh}=\lambda_{ebm}$, ie the magnitude of the couplings to nn (and further sites, if allowed) remains the same as for the corresponding EHM.

Before commenting on the sensitivity to the range of the coupling, it is interesting to note how peculiar are the shapes of some of these polaron dispersions, especially for the P+EBM+ model at weak and medium couplings in the adiabatic limit. The top left panel of Fig. 7, for instance, shows that the P+EBM- polaron dispersion has the usual expected shape at weak couplings, where the dispersion first roughly follows the bare dispersion and then flattens just below the polaron+one-phonon continuum that starts at $E_{0}+\Omega$. By contrast, the P+EBM+ polaron dispersion is located well below the continuum everywhere in the Brillouin zone, only reaching the continuum at $k=\pi$ for the weakest couplings considered. Given that the only difference is the overall sign of these EBM couplings, this major difference is clearly a consequence of the interference between the BM and the Peierls couplings. we discuss this in more details below.

Similar to the results obtained with EHM+P, we see little sensitivity of the polaron dispersion to the range for the larger $\Omega$, whereas as the adiabatic limit is approached this sensitivity appears again, especially at stronger $\lambda_{P}$. What is very striking, though, are the significant differences between the results for the three different EBM couplings (also from those for the corresponding EHM+P model). This is further demonstrated in Figs. 8-11, where we show the various GS properties of the polarons as a function of $\lambda_{P}$. All of these are quantitatively quite different for the three EBM flavors. They also show a different level of sensitivity to the range of the coupling, which depends not only on the ’flavor’ of EBM but also on the ground-state property considered. For instance, in Fig. 8, the dual P+EBM+ model shows significant differences in all panels for going from BM coupling to only nn sites, to EBM coupling to 2nd nn, and also to 3rd nn sites. The dual P+EBM- model only shows this sensitivity as the adiabatic limit is approached when $\lambda_{ebm}$ is increased, while the P+EBM0 coupling is the least sensitive to range. On the other hand, the top right panel of Fig. 10 reveals that for those parameters, the P+EBM- model has the larger sensitivity to range.

It is also interesting to note the unusual behavior of the P+EBM- model in the weak-coupling limit of $\lambda_{P}$, for the stronger $\lambda_{ebm}=2$ values. Here all ground-state quantities have rather unusual behavior, most strikingly $E_{gs}$ increases with $\lambda_{P}$ below the sharp transition. This is suggestive of a partial reciprocal cancellation of the effects of the Peierls and EBM+ couplings, consistent with previous work showing strong interference effects between the (short-range) BM and Peierls models.Roman In fact, we can follow the arguments in Refs. Dominic, ; Roman, to understand the behavior of these dual models in the strongly anti-adiabatic limit where $\Omega\gg\alpha,g_{bm}$. After projecting out the higher energy manifolds with one or more phonons, the resulting effective polaron dispersion is $E_{k}=\epsilon_{p}-2t^{*}\cos k-2t_{2}\cos 2k$. For all three flavors of EBM, we find the expected $\epsilon_{P}=-{4\alpha^{2}\over\Omega}-{g_{bm}^{2}\over\Omega}\sum_{|j|<M}f^{2}_{j}$ and the phonon-mediated 2nd nn effective hopping $t_{2}=-\alpha^{2}/\Omega$ which is the driver of the sharp transition in the polaron GS properties.Dominic ; DominicPRB However, we also find a renormalization of the nn hopping.Roman Specifically, $t^{*}=t$ for the P+EBM0 model, whereas $t^{*}=t\pm{2\alpha g_{bm}f_{1}\over\Omega}$ for the P+EBM$\pm$ couplings, demonstrating the interference between the Peierls and the nn BM couplings. For small $\alpha$ (weak $\lambda_{P}$) this renormalization of $t^{*}$ is the dominant change, explaining the already mentioned increase of $E_{gs}$ with $\lambda_{P}$ for the P+EBM- model. Of course, at larger $\lambda_{P}$, the cuadratic terms become more important, the dispersion becomes dominated by $t_{2}$ and the system transitions into the strong-coupling polaron regime. This argument also qualitatively explains the considerable differences in the values of $\lambda_{P}^{*}$ shown in Fig. 12. On the other hand, the sensitivity to the range of the EBM coupling, still clearly seen for the P+EBM+ model even at $\Omega=5$, is a higher-order effect whose explanation requires a more accurate approximation.