Single polaron properties for double-well electron-phonon coupling

# Single polaron properties for double-well electron-phonon coupling

Clemens P.J. Adolphs Department of Physics and Astronomy, University of British Columbia, Vancouver, British Columbia, Canada, V6T 1Z1    Mona Berciu Department of Physics and Astronomy, University of British Columbia, Vancouver, British Columbia, Canada, V6T 1Z1 Quantum Matter Institute, University of British Columbia, Vancouver, British Columbia, Canada, V6T 1Z4
July 29, 2019
###### Abstract

We show that in crystals where light ions are symmetrically intercalated between heavy ions, the electron-phonon coupling for carriers located at the light sites cannot be described by a Holstein model. We introduce the double-well electron-phonon coupling model to describe the most interesting parameter regime in such systems, and study it in the single carrier limit using the momentum average approximation. For sufficiently strong coupling, a small polaron with a robust phonon cloud appears at low energies. While some of its properties are similar to those of a Holstein polaron, we highlight some crucial differences. These prove that the physics of the double-well electron-phonon coupling model cannot be reproduced with a linear Holstein model.

###### pacs:
71.38.-k, 71.38.Ht, 63.20.kd, 63.20.Ry

## I Introduction

When charge carriers couple to phonons, magnons, or other bosonic excitations, the resulting dressed quasiparticles – the polarons – often behave drastically different from the free carriers. This is why understanding the consequences of carrier-boson coupling is important for many materials such as organic semiconductors,organic_1 (); Ciuchi and Fratini (2011) cuprates,Lanzara et al. (2001); Shen et al. (2004); Reznik et al. (2006); Lee et al. (2006); Gadermaier et al. (2010); Gunnarsson and Rosch (2008) manganites,Mannella et al. (2005) two-gap superconductors like MgB,Yildirim et al. (2001); Liu et al. (2001); Mazin and Antropov (2003); Choi et al. (2002) and many more. To describe them, many models of varying complexity have been devised and studied. The simplest is the Holstein model for electron-phonon coupling,Holstein (1959) where carriers couple to a branch of dispersionless optical phonons through a momentum-independent coupling . Physically, it describes a modulation of the on-site potential of the carrier due to the deformation of the “molecule” hosting it. Longer-range coupling that modulates the carrier’s on-site potential leads to couplings that depend on the boson’s momentum, such as the FröhlichFrÃ¶hlich (1954) or the breathing-mode models.bm () If the bosons modulate the hopping of the carrier, the coupling depends on the momenta of both carrier and boson, as is the case in the Su-Schrieffer-Heger (SSH) modelSSH (); extra () or for a hole coupled to magnons in an antiferromagnet, as described by a model.tJmod ()

All these electron-phonon coupling models assume that the coupling is linear in the lattice displacements. This is a natural assumption because if the displacements are small, the linear term is the most important contribution. However, the coefficient of the linear term may vanish due to symmetries of the crystal. In such cases, the most important contribution is the quadratic term.

Here we introduce, motivate and study in detail a Hamiltonian describing such quadratic electron-phonon (e-ph) coupling relevant for many common crystal structures, consisting of intercalated sublattices of heavy and light atoms. We focus on the single carrier limit and the parameter regime where the carrier dynamically changes the effective lattice potential from a single-well to a double-well; hence, we call this the double-well e-ph coupling. We use the momentum-average approximationBerciu (2006); Goodvin et al. (2006) to compute the properties of the resulting polaron with high accuracy. We find that although the polaron shares some similarities with the Holstein polaron, it also differs in important aspects. Indeed, we show that the physics of the double-well e-ph coupling model cannot be described by a renormalized linear Holstein model.

To the best of our knowledge, this is the first systematic, non-perturbative study of such a quadratic model. PreviouslyAdolphs and Berciu (2013) we studied the effect of quadratic (and higher) corrections added to a linear term. Weak, purely quadratic coupling was studied using perturbation theory in Refs. Entin-Wohlman et al., 1983, Entin-Wohlman et al., 1985. Other works considered complicated non-linear lattice potentials and couplings but treated the oscillators classically,Zolotaryuk et al. (1998); Raghavan and Kenkre (1994); Kenkre (1998) or discussed anharmonic lattice potentials but for purely linear coupling.Voulgarakis and Tsironis (2000); Maniadis et al. (2003) Away from the single-carrier limit, the Holstein-Hubbard model in infinite dimensions was shown to have parameter regions where the effective lattice potential has a double-well shape;Freericks et al. (1993); Meyer et al. (2002); Mitsumoto and Ōno (2005) this was then used to explain ferroelectricity in some rare-earth oxides.[See][andrelatedreferences]double_well_2 However, the effect of a double-well e-ph coupling on the properties of a single polaron were not explored in a fully quantum-mechanical model on a low-dimensional lattice.

This work is organized as follows: in Section II we introduce the Hamiltonian, motivate its use for relevant systems, and discuss all approximations made in deriving it. In Section III we review the theoretical means by which we study our Hamiltonian. In Section IV we present our results, and in Section V we give our concluding discussion and an outlook for future work.

## Ii The model

The crystal structures of interest are illustrated in Fig. 1(a) for 1D, and Fig. 1(b) for 2D cases. The 3D crystal would have a perovskite structure but we do not discuss it explicitly because, as we show below, dimensionality plays no role in determining the polaron properties.

The undoped compound is an insulator made of light atoms, shown as filled circles, intercalated between heavy ones, shown as empty circles. To zeroth order, the vibrations of the heavy atoms can be ignored while those of the light atoms are described by independent harmonic oscillators , where annihilates a phonon at the light atom. (We set the mass of the light ions , and also ). In reality there is weak coupling between these oscillators giving rise to a dispersive optical phonon branch. However, the dispersion can be ignored if its bandwidth is small compared to all other energy scales. We do so in the following.

Consider now the addition of a carrier. If it occupies orbitals centered on the heavy atoms, its coupling to the oscillations of the light atoms is described by breathing-mode coupling models.bm () Here we are instead interested in the case where the carrier is located on the light atoms. Such is the situation for a CuO plane as shown in Fig. 1(b), since the parent compound is a charge-transfer insulatorZSA () so that upon doping, the holes reside on the light O sites (of course, there are additional complications due to the magnetic order of the Cu spins; we ignore these degrees of freedom in the following). The carrier moves through nearest-neighbor hopping between light atoms: , where is the carrier annihilation operator at light atom .

Given the symmetric equilibrium location of the light ion hosting the carrier between two heavy ions, it is clear that the e-ph coupling cannot be linear in the displacement of that light ion: the sign of the displacement cannot matter. Thus, e-ph coupling in such a material is not described by a Holstein model. This assertion is supported by detailed modelling. For simplicity, we assume that the interactions with the neighboring heavy atoms are dominant (longer-range interactions can be easily included but lead to no qualitative changes). There are, then, two distinct contributions to the e-ph coupling:

#### Electrostatic coupling:

The carrier changes the total charge of the light ion it resides on. If the distance between adjacent light and heavy ions is , and if is their additional Coulomb interaction due to the carrier, then the potential increases by . This is an even function and thus has no linear (or any odd) terms in . The coefficient of the quadratic term can be either positive or negative, depending on the charge of the carrier (electron or hole).

#### Hybridization:

Even though charge transport is assumed to take place in a light atom band, there is always some hybridization allowing the carrier to hop onto an adjacent heavy ion. If is the corresponding energy increase, assumed to be large, then the carrier can lower its on-site energy by through virtual hopping to a nearby heavy ion and back. The hopping depends on the distance between ions; for small displacements where is some material-specific constant. Because the light ion is centered between two heavy ions, such contributions add to The potential is again even in . In this case, the coefficient of the quadratic term is always negative.

Given that , it follows that the largest (quadratic) contribution to the e-ph coupling for such a crystal has the general form:

 H(2)e-ph=g2∑ic†ici(bi+b†i)2

where all prefactors have been absorbed into the energy scale , and the sum is over all light ions. From the analysis above we know that may have either sign.

Physically, shows that the presence of a carrier modifies the curvature of its ion’s lattice potential, and thus changes the phonon frequency at that site from to . If then , making phonon creation more costly. As we show in Appendix C, this leads to a rather uninteresting large polaron with very weakly renormalized properties. This is why in the following we focus on the case with .

For sufficiently negative , vanishes or becomes imaginary, i.e. the lattice is unstable. This is unphysical; in reality the bare ionic potential contains higher order terms that stabilize the lattice. This means that for we must include anharmonic (quartic) terms in the phonon Hamiltonian and, for consistency, also in the e-ph coupling, so that

 Hph =Ω∑ib†ibi+Θ∑i(b†i+bi)4 H(4)el-ph =∑n∈{2,4}gn∑ic†ici(b†i+bi)n,

where is the scale of the anharmonic corrections. In physical situations and , or the Taylor expansions would not be sensible starting points.

The anharmonic terms in make the total Hamiltonian unwieldy, because the phonon vacuum is no longer the undoped ground-state, and the new undoped ground state has no simple analytical expression. In order to be able to proceed with an analytical approximation, we argue that these terms can be absorbed into the e-ph coupling; this is a key approximation of the model. The reasoning is as follows: At those lattice sites that do not have a carrier, the quartic terms have little effect if . This statement is verified by exact diagonalization of . Results are shown in Fig. 2 where we plot the overlap (per site) between the undoped ground-states with and without anharmonic corrections, as well as the average number of phonons at a site of the undoped lattice. Even for unphysically large values , the overlap remains close to 1 while , showing that the undoped ground-state has not changed significantly in the presence of anharmonic corrections. From now we ignore these corrections at sites without an additional carrier.

However, for sites that have a carrier present, we cannot ignore the anharmonic term: As discussed, it is crucial for stabilizing the lattice. Since this term is similar to the quartic term in the e-ph coupling, they can both be grouped together, resulting in the approximate Hamiltonian for our crystal:

 H=^T+Ω∑ib†ibi+g2∑ic†ici(b†i+bi)2+(g4+Θ)∑ic†ici(b†i+bi)4 (1)

with an effective quartic e-ph coupling term , which from now on we will simply call . This is the Hamiltonian that we investigate in this work.

Before proceeding, let us review what we are neglecting when we discard the anharmonic corrections at the unoccupied sites. Besides ignoring the change in the undoped ground state from to (which is a reasonable approximation if , as discussed), we also assume that only the e-ph coupling can change the number of phonons in the system, whereas in the full model the phonon number is also changed by anharmonic corrections. This latter approximation is valid if the timescale for anharmonic phonon processes is much longer than the characteristic polaron timescale , where is the effective polaron mass.

Let us briefly summarize the basic properties of the lattice potential, which equals for sites without an extra carrier, and for sites with one carrier. If , the first term describes a harmonic well with frequency and describes a single well centered at . If , however, becomes purely imaginary. In this case, becomes a double-well potential with a local maximum at . The two wells are centered at For we obtain , which locally describes a harmonic well of frequency . Interestingly, this is independent of , whose only role is to control the location and depth of the two wells (they are further apart and deeper for smaller ).

## Iii Formalism

We want to find the single particle Green’s function , where is the resolvent of Hamiltonian (1). From this, we can obtain all the polaron’s ground state properties as well as its dispersion.(Goodvin et al., 2006)

Grouping terms in the Hamiltonian according to how they affect the phonon number, we rewrite with and do not change the number of phonons, while and change it by and , respectively. The constant in is absorbed into in the following derivations, but plots of the spectral weight will show actual energies.

One important property of this Hamiltonian is that it preserves the phonon number parity on each site: because its terms only change the number of phonons by multiples of two, any eigenstate is a sum of basis states having only even (or only odd) number of phonons. The Hilbert space can thus be divided into an even and an odd (phonon number) sector, which can be diagonalized separately. We emphasize that this symmetry is different from the parity symmetry under a global lattice inversion . The latter has been studied extensively for the linear Holstein model,linear_symmetry () where it was shown that polaron states with total momentum have well defined (spatial) parity. The phonon number parity, on the other hand, corresponds to a unitary transformation , i.e., a local inversion of the phonon coordinates. The number parity symmetry also correlates with the local spatial parity of the ions, since the spatial parity operator for site can be written as .

### iii.1 The even sector

We compute the Green’s function via the same continued matrix fractions methodMÃ¶ller et al. (2012) previously used by us to compute the Green’s function of a generalized Holstein model with linear and higher-order termsAdolphs and Berciu (2013) within the framework of the momentum average (MA) approximation. This approximation was shown to be highly accurate for models with Holstein coupling.Berciu (2006); Goodvin et al. (2006) The reasons for this (such as obeying exact sum rules) can be verified to hold for this model, too. To be specific, here we implement the MA flavor which allows us to also locate the continuum lying above the polaron band.Berciu and Goodvin (2007)

We begin our derivation by dividing the Hamiltonian into with . Using Dyson’s identity , where , we obtain

 G(k,ω)=G0(k,ω)[1+∑ieikRi√N(g2+6g4)F1(k,ω;i,i)+g4F2(k,ω;i,i)] (2)

with being the generalized propagator for a system with phonons in total, of them on site with the other two on site . The difference between MA and the original MA, which we also call MA, is that for we also use its exact equation of motion (EOM),

 F1(k,ω;i,j)=G(k,ω;j)G0(j−i,ω−2Ω)(2g2+12g4)+F1(k,ω;j,j)G0(j−i,ω−2Ω)(4g2+36g4)+8g4F2(k,ω;j,j)G0(j−i,ω−2Ω)+∑lG0(l−i,ω−2Ω)[F2(k,ω;l,j)(g2+6g4)+F3(k,ω;l,j)g4]. (3)

which is obtained by applying Dyson’s identity again, and introducing and . The equations of motion for the propagators with are approximated by replacing the free propagator , where is the momentum averaged free propagator. At low energies this is a good approximation because decays exponentially with the distance if in dimensions. This is also justified by the variational meaning of the MA approximations, discussed at length elsewhere.bar (); Berciu and Goodvin (2007) (Basically, MA assumes that all phonons in the cloud are at the same site but also allows for a pair of phonons to be created at a site away from the cloud).

The resulting EOMs are different depending on whether or . If we define and for , we obtain

 F=n(k,ω;i)=¯g0(ω−2nΩ)[F=n−2(2n)¯4g4+F=n−1((g2+6g4)(2n)¯2+4g4(2n)¯3)+(4ng2+12ng4+24n2g4)F=n+(g2+6g4+8ng4)F=n+1+g4F=n+2]. (4)
 F≠n(k,ω;i,j)=¯g0(ω−2nΩ)[g4(2n−2)¯4F≠n−2+((g2+6g4)(2n−2)¯2+(2n−2)¯3⋅4g4)F≠n−1+[2(2n−2)g2+12(n−1)g4+6(2n−2)2g4]F≠n+[g2+6g4+4(2n−2)g4]F≠n+1+g4F≠n+2] (5)

where we use the notation . We also omitted the arguments from the appearing on the right hand sides, as they remain unchanged.

These EOMs connect generalized Green’s functions with and . We reduce this to a first order recurrence relationAdolphs and Berciu (2013) by introducing vectors and analogously for . Below, we write without the index or for results that apply to both and . By inserting the EOMs into the definition of , we obtain a matrix EOM for the ,

 γnWn=αnWn−1+βnWn+1. (6)

The coefficients of these matrices are read off from the EOM for the . They are listed in appendix A.1.

Using the fact that we can showAdolphs and Berciu (2013) that with . By introducing a suitably large cut-off where we set , we can compute all with as continued matrix fractions. Knowledge of allows us to express and in terms of and . Following a series of steps presented in appendix A.2, we obtain a closed equation for in terms of , which we then finally use to compute . The end result of these manipulations is the self energy

 Σ(ω)=(g2+6g4+A=1|12g4)~g0(ω)a=01−~g0(ω)(a=1−a≠)+g4A=1|11.

with ) and the other coefficients defined in appendix A.2. The independence of the self-energy on momentum is the consequence of the local form of the coupling and of the non-dispersive phonons, similar to the MA results for the Holstein model.Berciu and Goodvin (2007) Momentum-dependence would be acquired in a higher flavor of MA, but is likely to be weak. Finally, the Green’s function is:

 G(k,ω)=1ω−ϵk−Σ(ω)+iη. (7)

One can now use the matrices to generate the generalized propagators , which allow one to reconstruct the entire polaron wavefunction (within this variational space).Can () For the quantities of interest here, however, the single-particle Green’s function suffices.

### iii.2 The odd sector

Here we calculate the Green’s function for a state that already has a phonon in the system. Since the phonon number can only change by or , this single phonon can never be moved to another site, so it is natural to compute the Green’s function in real space. The most general such real space Green’s function is:

Applying the Dyson identity leads to the EOM

We then split the sum over all lattice sites into a term where the electron is on the same site as the extra phonon, and a sum over all the other sites. The subsequent steps are very similar to those for the even-sector Green’s function. We summarize them in Appendix B, where we also discuss how various propagators that enforce translational symmetry – i.e. propagators defined in momentum space – can be obtained from these real-space Green’s functions.

The end result for the real-space Green’s functions is where . The coefficients and are listed in appendix B.

## Iv Results

### iv.1 Atomic limit: t=0

We begin our analysis with the atomic limit since it is a good starting point for understanding the properties of the small polaron, which is the more interesting regime. However, we note an important distinction between the Holstein model and our double-well model. In the former, the atomic limit is the infinite-coupling limit. In the latter, sets an additional energy scale. Thus, the atomic limit is not the same as the strong coupling limit; the latter also requires that be small.

Before doing any computations, we can describe some general features of the spectrum. As already discussed, the phonon component of the wavefunctions has either even or odd phonon number parity. Since this is due to the spatial symmetry in the local ionic displacement, in any eigenstate the ion is equally likely to be found in either well. As usual, the ground state has even symmetry since it has no nodes in its wavefunction. Subsequent eigenstates always have one more node than the preceding eigenstate, so states with even and odd parity alternate. The exception is the limit of infinite well separation, , where the and eigenstates become degenerate. The system can then spontaneously break parity to have the ion definitely located in the left or in the right well, like in a ferroelectric. For a finite this is not possible in the single carrier limit, but it can be achieved at finite carrier concentration through spontaneous symmetry breaking.

As discussed, our results are obtained with MA. In the atomic limit MA is exactBerciu (2006) because for the free propagator is diagonal in real-space so the terms ignored by MA vanish. Thus, MA results must be identical here to those obtained by other exact means. To check our implementation of MA, we used exact diagonalization (ED) with up to a few thousand phonons; this suffices for an accurate computation of the first few eigenstates. ED and MA results agree, as required.

Figure 3 shows the ground-state quasiparticle weight (the overlap between the polaron ground-state and the non-interacting carrier ground-state), and the ground-state average number of phonons in the cloud, , as a function of , for various values of . has an interesting behavior. At it is slightly below because of the quartic terms. As is increased, first rises towards a value close to and then sharply drops. This turnaround is caused by the terms that involve both and , i.e. from and from . Starting from and making it increasingly more negative will at first decrease these coefficients, thereby renormalizing the ground state less. For even more negative , however, decreases sharply as the absolute value of these coefficients increases; this is paralleled by a strong increase in . Based on this argument, the peak in should occur for , which is indeed the case. The strong-coupling limit of a small polaron (corresponding to small , large values) is therefore reached either by increasing or by lowering .

While this allows us to conclude that in the atomic limit the crossover into the small polaron regime occurs at , it also illustrates the difficulty in defining an effective coupling for this model. For the Holstein model, the dimensionless effective coupling is the ratio between the ground-state energies in the atomic limit and in the free electron limit; the crossover to the small polaron regime occurs at . For the double-well model the introduction of an effective coupling is not as straightforward, because the atomic limit has vastly different properties depending on the ratio , so comparing the energy in this limit to that of a free electron is not sufficient. (Moreover, there is no analytic expression for the ground state energy of the double well potential). For these reasons, we continue to use the bare coupling parameters and to characterize our model.

For strong coupling, we can accurately estimate the ground state energy by using the barrier depth and effective harmonic frequency of the double-well potential, . Fig. 4 shows the relative error of this estimate, which indeed decreases as parameters move deeper into the small polaron regime. Since here the tunnelling between the two wells also becomes increasingly smaller, one may think that we can describe this regime accurately by assuming that the carrier becomes localized in one of the wells (thus breaking parity), i.e. that we can approximate the full lattice potential as being a single harmonic well centered at either or . Of course, the latter situation can be modelled with a linear Holstein model.

It turns out that this is not the case. In the standard Holstein model, the charge carrier cannot change the curvature of the lattice potential and thus cannot account for the difference between and . To account for the change in the curvature of the well, one would have to consider at least a Holstein model with both linear and quadratic e-ph coupling terms. Although it is possible to find effective parameters , and so that the resulting lattice potential in the presence of the carrier has the same location and curvature as one of the wells of the double-well potential, the corresponding quasi-particle weight severely underestimates . This is because the single well approximation severely overestimates the lattice potential at , thereby reducing the overlap between the ground state of the shifted well and that of the original well. We conclude that the double-well coupling cannot be accurately described by a (renormalized) Holstein coupling even in this simplest limit.

### iv.2 Finite Hopping

We focus on results from the even sector because it describes states accessible by injecting the carrier in the undoped ground-state. The odd sector is accessed only if the carrier is injected into an excited state with an odd number of phonons present in the undoped system; we briefly discuss this case at the end of the section.

We begin by plotting the ground-state values of and , for 1D and 2D lattices, in Figs.  5 and 6 respectively. Since the MA self-energy is local, the effective polaron mass , where is the free carrier mass; we therefore do not plot separately. Apart from , the parameters are like in Fig. 3. Note that the kinks in the curves for are not physical; they arise from numerical difficulties in resolving the precise location of the ground state peak when .

Qualitatively, the polaron properties show the same dependence on as in the atomic limit, but the shape and location of the turnarounds is slightly modified: As one would expect, the presence of finite hopping counteracts the formation of a robust polaron cloud and increases the quasi-particle weight for any given and when compared to the atomic limit.

The results in one and two dimensions are strikingly similar. The 2D is slightly larger than the 1D , and in 2D is slightly lower than in 1D. This is expected because in higher dimensions, the polaron formation energy is competing against a larger carrier kinetic energy. These results suggest that dimensionality is not playing a key role for the double-well model, similar to the situation for the Holstein model. This is why we did not consider 3D systems explicitly.

We now move on to discuss the evolution of the spectral weight with increasing , at a fixed value of . This is shown in Fig. 7 for 1D, and in Fig. 8 for 2D. Because the evolution is again qualitatively similar in the two cases, we analyze in more detail the 1D results. Here, at small quadratic coupling , we observe the appearance of a polaron band below a continuum of states. This continuum begins at , and consists of excited states comprising the polaron plus two phonons far away from it. (In our MA approximation, the continuum actually begins at , not at , for reasons detailed in Ref. Berciu and Goodvin, 2007).

Note that due to the parity-preserving nature of the Hamiltonian there is no analog of the polaron+one-phonon continuum starting at , which is observed in all linear coupling models. Trying to mimic the results of the double-well coupling with a linear model will, therefore, lead to a wrong assignment for the value of .

At small , the polaron band flattens out just below the polaron+two-phonon continuum. With increasing , its bandwidth decreases as the polaron becomes heavier, and additional bound states appear below the continuum. This is similar to the evolution of the spectrum of a Holstein polaron when moving towards stronger effective coupling.Goodvin et al. (2006) However, as already discussed, this does not mean that the two Hamiltonians can be mapped onto one another.

For completeness, let us also discuss some of the features of the odd sector. In particular, we focus on the local Green’s function , which can be written as

 Giii(ω)=¯g0(~ω)+¯g0(~ω)2(a=−a≠)1−¯g0(~ω)(a=−a≠)

with . One can verify that equals the MA self-energy for the even sector, up to a shift by of its frequency. The equation for then shows that the odd sector spectral function comprises two parts: (i) the first term is just the momentum-averaged spectral function of the even-sector, shifted in energy by due to the presence of the extra phonon. One can think of these as states where the even-sector polaron does not interact with the extra phonon. This contribution therefore has weight starting from ; (ii) the second part describes interactions between the polaron and the extra phonon. An interesting question is whether these can lead to a bound state, i.e. to a new polaron with odd numbers of phonons in its cloud.

This question is answered in Fig. 9 where we plot for different values of and , , , in one dimension. The vertical bars indicate the position of , where indeed a continuum begins, as expected from the previous discussion. At sufficiently strong coupling we find a discrete bound state below that continuum, showing that the polaron can bind the extra phonon. In fact, it is more proper to say that the extra phonon (which is localized somewhere on the lattice) binds the polaron to itself and therefore localizes it. One can think of this as an example of “self-trapping”, except here there is an external trapping agent in the form of the extra phonon.

One might wonder whether this localized bound state in the odd sector could ever be at an energy below the polaron ground-state energy of the even sector, i.e. become the true ground-state. This is not the case; as explained above, in the atomic limit the ionic states alternate between even and odd symmetry. Introducing a finite hopping allows the polaron to further lower its energy by delocalizing, but this is only possible in the even sector. Thus, we always expect the even-sector polaron to have an energy below that of this localized state.

As stated before, the two subspaces with even and odd phonon number are never mixed, at least at zero temperature. At finite temperature, the extra charge is inserted not into the phonon vacuum but into a mixed state containing a number of thermally excited phonons. We therefore expect the resulting spectral function to show features of both the even and odd sectors. To be more precise, some spectral weight should be shifted from the even-sector spectral weight to the odd-sector spectral weight as increases and there is a higher probability to find one or more thermal phonons in the undoped state. We plan to study the temperature depend properties of this double-well coupling elsewhere.

## V Summary and Discussions

Here we introduced and motivated a model for purely quadratic e-ph coupling, relevant for certain types of intercalated lattices, wherein the carrier dynamically changes the on-site lattice potential from a single well into a double well potential. All the approximations made in deriving this model were analyzed. In particular, we argued that ignoring the anharmonic lattice terms at the sites not hosting the carrier should be a good approximation. However, a more in-depth numerical analysis might be needed to further validate this assumption.

We used the momentum average approximation to obtain the model’s ground state properties and its spectral function in the single polaron limit, in one and two dimensions. We found that for sufficiently strong quadratic coupling a small polaron forms. Although the polaron behaves somewhat similarly to the polaron of the linear Holstein model, the double-well model cannot be mapped onto an effective linear model: apart from the difference in the location of the continuum in the even sector, the double-well model also has an odd sector that should be visible at finite , and which is entirely absent in the Holstein model. This is due to the double-well potential model’s invariance to local inversions of the ionic coordinate; this symmetry is not found in the Holstein model. The polaron in this odd sector is also qualitatively different from the Holstein polaron, in that it is localized near the additional phonon present in the system when the carrier is injected. Of course, if the assumption of an Einstein mode is relaxed, then the phonon acquires a finite speed and this polaron would become delocalized, as expected for a system invariant to translations. However, this would still be qualitatively different than a regular polaronic solution because this polaron’s dispersion would be primarily controlled by the phonon bandwidth, not the carrier hopping.

Our results suggest that researchers interpreting their measurements from, e.g., angular-resolved photoemission spectroscopy, must carefully consider the nature of their system’s e-ph coupling: if they assume linear coupling where the lattice symmetry calls for a quadratic one, the parameters extracted from fitting to such models will have wrong values.

While we have laid here the basis for a thorough investigation of the properties of the double-well e-ph coupling model, much work remains. We believe that adjusting already existing numerical schemes such as diagrammatic Monte Carlo to this model is straightforward and look forward to a comparison of numerically exact results with our MA results. In addition, there are certain ranges of parameters for which MA is not well-suited, such as the adiabatic limit at weak coupling, or systems with finite carrier densities. We anticipate that these regimes will be explored with a range of numerical and analytical tools, especially the finite carrier regime which should be relevant for modelling ferroelectric materials.

We plan to extend our study of the double-well e-ph coupling beyond the single-polaron limit. We deem especially interesting the parameter range where the lattice potential remains a single well if only one carrier is present, but changes into a double well when a second charge is added. In this case, we anticipate the appearance of a strongly bound bipolaron while the single polarons are relatively light. Such states are not possible in the Holstein model.

Finally, extending our MA treatment to finite temperature should yield interesting insights into the interplay between the two symmetry sectors revealed by the spectral weight.

###### Acknowledgements.
We thank NSERC and QMI for financial support.

## Appendix A Details for the even-sector

### a.1 Coupling matrices

The matrices appearing in Eq. (6) are:

 γ=n|11 =1−¯g0(ω−4nΩ)(8ng2+24ng4+96n2g4) γ=n|12 =−¯g0(ω−4nΩ)(g2+6g4+16ng4) γ=n|21 =−¯g0(ω−(4n+2)Ω)((g2+6g4)(4n+2)¯2+ =4g4(4n+2)¯3) γ=n|22 =1−¯g0(ω−(4n+2)Ω)((8n+4)g2+ =(24n+12)g4+24(2n+1)2g4)
 α=n|11 =¯g0(ω−4nΩ)(g4(4n)¯4) α=n|12 =¯g0(ω−4nΩ)((g2+6g4)(4n)¯2+4g4(4n)¯3) α=n|21 =0 α=n|22 =¯g0(ω−(4n+2)Ω(g4(4n+2)¯4)
 β=n|11 =¯g0(ω−4nΩ)g4 β=n|12 =0 β=n|21 =¯g0(ω−(4n+2)Ω)(g2+6g4+(16n+8)g4) β=n|22 =¯g0(ω−(4n+2)Ω)g4

The matrices for sector are the same if we substitute everywhere except in the argument of .

### a.2 Manipulation of the EOMs

We can rewrite the EOM of by inserting the matrices and and collecting terms. This results in

 F1(ij)=G0(j−i,ω−2Ω)[a=0G(j)+a=1F=1(j)]+∑l≠jG0(l−i,ω−2Ω)a≠F≠1(lj). (8)

where we omit the arguments and for shorter notation. We give expressions for the various coefficients below. For now, we rewrite the EOM as

 F1(ij)=G0(j−i,ω−2Ω)×[a=0G(j)+(a=1−a≠1)F=1(j)]+∑lG0(l−i,ω−2Ω)a≠F1(lj). (9)

Defining , we can write this as a matrix product:

 ∑l[δil−a≠G0(ω−2Ω)il]F1(lj)=G0(ω−2Ω)ij[a=0G(j)+(a=1−a≠1)F1(jj)].

We multiply this from the left with and obtain

 ∑l[G−10(ω−2Ω)rl−a≠δrl]F1(lj)=δrj[a≠0G(j)+(a=1−a≠)F1(jj)].

Next, we use the fact that , so subtracting from this just shifts its frequency to obtain . As a result:

 F1(ij)=G0(ω−2Ω−a≠)ij×[a=0G(j)+(a=1−a≠1)F1(jj)].

Since in the EOM for we only require , we solve for that diagonal element and obtain

 F1(jj)=¯g0(ω−2Ω−a≠)a=0G(j)1−¯g0(ω−2Ω−a≠)(a=1−a≠1).

The coefficients are obtained by just inserting the appropriate matrices into the EOM and collecting terms:

 a=0 =2g2+12g4+(g2+14g4)A=1|11+g4A=1|21 a=1 =4g2+36g4+(g2+14g4)A=2|12+g4A=2|22 a≠ =(g2+6g4)A≠1|12+g4A≠1|22

Finally, are used in Eq. (2) to obtain .

## Appendix B Details for the odd-sector

### b.1 Equations of Motion

Starting from the EOM for , we let act on the states in those sums, to find for the diagonal state:

 H1c†lb†l|0⟩=(2g2+12g4)c†lb†l|0⟩+(g2+10g4)c†lb†,3l|0⟩+g4c†lb†,5l|0⟩

while for the off-diagonal ones:

 H1c†i′b†l|0⟩=(2g2+6g4)c†i′b†l|0⟩+(g2+6g4)c†i′b†,2lb†l|0⟩+g4c†i′b†,4lb†l|0⟩.

We now define the generalized Green functions as:

so we always have the extra phonon at site . The equation of motion for then becomes: Again, we start by separating the cases and . The resulting equations of motion for are like those of the even-sector with , while those for are like those of the even-sector with .

In the spirit of MA, only the EOM for , which already has one phonon present, is kept exact, while in the EOMs for all the with we approximate . We introduce matrices . Again we obtain an equation like Eq. (6), where now:

 γ=11 =1−¯g0(ω−(4n−1)Ω)((4n−1)(2g2+6g4 =+6g4(4n−1)) γ=12 =−¯g0(ω−(4n−1)Ω)(g2+6g4+4g4(4n−1)) γ=21 =−¯g0(ω−(4n+1)Ω)((4n+1)¯2(g2+6g4) =+(4n+1)¯3⋅4g4) γ=22 =1−¯g0(ω−(4n+1)Ω)(4n+1) =×(2g2+6g4+6g4(4n+1))
 α=11 =¯g0(ω−(4n−1)Ω)(4n−1)¯4g4 α=12 =¯g0(ω−(4n−1)Ω) =×((4n−1)¯2(g2+6g4)+(4n−1)¯3⋅4g4) α=21 =0 α=22 =¯g0(ω−(4n+1)Ω)(4n+1)¯4g4
 β=11 =¯g0(ω−(4n−1)Ω)g4 β=12 =0 β=21 =¯g0(ω−(4n+1)Ω)(g2+6g4+4g