A semi exact solution for a metallic phase in a Holstein-Hubbard chain at half filling with Gaussian anharmonic phonons

The nature of phase transition from an antiferromagnetic SDW polaronic Mott insulator to the paramagnetic bipolaronic CDW Peierls insulator is studied for the half-filled Holstein-Hubbard model in one dimension in the presence of Gaussian phonon anharmonicity. A number of unitary transformations performed in succession on the Hamiltonian followed by a general many-phonon averaging leads to an effective electronic Hamiltonian which is then treated exactly by using the Bethe-Ansatz technique of Lieb and Wu to determine the energy of the ground state of the system. Next using the Mott–Hubbard metallicity condition, local spin-moment calculation, and the concept of quantum entanglement entropy and double occupancy, it is shown that in a plane spanned by the electron–phonon coupling coefficient and onsite Coulomb correlation energy, there exists a window in which the SDW and CDW phases are separated by an intermediate phase that is metallic.

www.nature.com/scientificreports/ Feshke et al. 14 have also implemented the DMRG method and established the occurrence of the metallic regime between the two insulating phases. They have also proved that intermediate metallic regime widens as the phonon frequency increases. Several other studies using the renormalization group (RG) technique 16 , Monte-Carlo simulations 17 , exact numerical diagonalization and cluster perturbation theory 18 etc. have also shown the evidence of the intervening metallic phase between the SDW and CDW states. Tezuku et al. 19 have studied the HH model with strong e-e and strong e-p interactions with DMRG. They have considered the region between the anti-adiabatic and adiabatic limits and observed that when the e-e and e-p interactions are of comparable strength, the pairing and CDW correlations are degenerate. They have furthermore shown that when the phonon energy scale is much higher than the e-p interaction scale and the electron-hole symmetry is broken, the on-site superconducting phase overlaps with CDW state and then the transition from SDW to CDW does not require any intermediate phase. In a modified study, Tezuka et al. 20 have calculated the correlation functions with the real space dynamics. Though for the pure (un-doped) HH model, they have found a metallic region in between the SDW-CDW phases, for the doped HH model with broken electron-hole symmetry, they have found the pairing correlation to be more dominant.
Tam et al. 21 have studied the 1D HH model at half-filling using renormalization group (RG) method. They have considered the e-p interaction and e-e interaction on the same footing and considered all the possible retardation effects of the phonon dynamics. Their study using mean-field RG suggests a direct transition from the CDW state to the SDW state.
The above studies suggest that more detailed analytical investigations are required to predict unequivocally the existence of the intermediate metallic phase at the CDW-SDW cross-over region. To that end, Chatterjee and collaborators [22][23][24][25][26][27][28] have carried out a few improved variational calculations and have shown that with every improvement of the variational wave function, the width of the intermediate metallic phase increases. This certainly lends a fair amount of credence to the original prediction of TC.
Chatterjee and Takada (CT) 28 have also examined the problem in the presence of lattice anharmonicity. They have considered cubic and quartic phonon anharmonicities and have shown that the metallic phase becomes broader in the presence of lattice anharmonicity and thus the conjecture on the presence of the intervening metallic regime in the HH system is strengthened in the presence of anharmonic phonons. This work is of much importance because lattice anharmonicity has been found to play a crucial role in high T c superconductors. In fact, it has been observed that apex oxygen has a substantial anharmonic motion in the cuprates and also the phonon anharmonicity makes a significant impact on the electronic structure of these systems [29][30][31][32][33] . Konior 34 has explained the importance of Gaussian phonon anharmonicity in the context of high T c superconductors and have shown that in the presence of this anharmonicity, the hopping parameter reduces at a slow rate causing an enhancement in the polaron mobility and the polaron bandwidth, which is a favourable condition for the phonon mechanism to stake a claim for inducing pairing. Lavanya, Sankar and Chatterjee (LSC) 25 have recently re-examined the work of CT with Gaussian anharmonic potential by applying in succession a number of unitary transformations followed by an averaging with a general many-phonon state and the Bethe ansatz technique. This gives a wider metallic phase.
The principal aim of the present paper is to further modify the variational wave function of the phonon sub-system used by LSC for the anharmonic HH system to obtain a better solution for the GS energy and the SDW-CDW phase diagram. This calculation can be considered as semi-exact as we have included rigorously all possible phonon processes including coherence and correlations while treating the phonon subsystem and solved the electronic part exactly with the help of the Bethe ansatz technique and LW's solution. The GS energy, local spin moment, von Neumann entropy, the double occupancy parameter and the phase diagram at the SDW-CDW transition region have been obtained.

The Model
A one dimensional HH system with Gaussian phonon anharmonicity may be described by the Hamiltonian with where H e describes the Hubbard Hamiltonian, H p is the phonon Hamiltonian and H ep is the Holstein e-p interaction. In Eq. (2), the parameter t is the nearest-neighbour hopping integral, the operator c † iσ (c iσ ) creates (annihilates) a spin-σ electron at the i-th site, n iσ = c † iσ c iσ being the corresponding electron occupation number and U gives the onsite e-e interaction energy. In Eq. (3), b † i (b i ) represents an operator that creates (annihilates) an optical phonon at site i with dispersionless frequency ω 0 , ap and γ measure respectively the strength and range www.nature.com/scientificreports/ of the phonon anharmonicity. In Eq. (4), g is the on-site e-p coupling strength which can be written as: g = √ αω 0 , where dimensionless α is referred to as the e-p coupling constant.

Formulation
GS energy. In order to solve the Hamiltonian (1) we choose to seek a variational solution. First of all, we apply the modified Lang-Firsov transformation (LFT) 35 with the generator where η is the variational parameter that carries the information of the polaronic structure. For strong e-p interaction, η → 1 , and Eq. (5) generates the usual LFT 36 and gives a reasonable approximation for the antiadiabatic region in which the ion motion is much faster than the electron motion. This should work well for the narrow-band materials which are essentially strongly correlated systems. Because of the above transformation, the Hamiltonian H transforms to H 1 = e R 1 He −R 1 . To deal with the adiabatic regime where the electron motion is much faster than the ion motion, we perform the Takada-Chatterjee (TC) transformation 10 with the generator: where we assume, h i = h , as all sites are equivalent. After the second transformation, the transformed Hamiltonian becomes: H 2 = e R 2 H 1 e −R 2 . The above two transformations together can be generated by: With η = 1, the transformation (7) represents the conventional LFT which gives exact results in the antiadiabatic limit, while for η = 0 , it takes care of the adiabatic limit. Thus both the anti-adiabatic and the adiabatic regions can be studied by considering: 0 < η < 1. Thus η can be called an adiabaticity parameter. It is important to note that Eq. (7) assumes the phonons associated with the electron to be in a coherence state. This is essentially a semi-classical approximation in which it is assumed that the phonons in the polaron cloud are independent of each other satisfying a Poissonian distribution. In other words, the phonons emitted or absorbed by the electrons are completely uncorrelated and in that sense, the present transformation is equivalent to the Hartree approximation.
The presence of Gaussian anharmonicity in the system introduces anharmonicity to infinite order and results in a finite lifetime of the phonons through phonon-phonon interactions. Furthermore, an electron undergoes a recoil motion while emitting a phonon. While undergoing a recoil motion, if the electron emits another phonon, then these two successively emitted virtual phonons will be correlated. Both the anharmonicity and the correlation effects of the phonons can be taken into account by considering squeezed phonons. Squeezing of the phonon vacuum state can be accomplished by the celebrated Bogoliubov transformation with the generator 10 : where the squeeze parameter α s is to be obtained variationally. The transformed Hamiltonian is now given by: The variational parameter α s has been considered by all investigators as independent of the electron concentration until 2019, when Malik, Mukhopadhyay and Chatterjee (MMC) 26 considered the squeezing of the phonon state to be partly dependent on electron density. According to MMC, the correlation between phonons emitted by the electrons may depend on the number of electrons available at a particular lattice site. Thus, to be more realistic, we next apply, a squeezing transformation with the generator: where α d is the variational parameter. Correspondingly, the transformed Hamiltonian can be written as: It may be pointed out that in Eq. (8), phonon correlation and anharmonicity have been included at a mean-field level while Eq. (9) incorporates the fluctuations. The principal idea of performing a series of canonical transformation on H is to disentangle the electrons and phonons and if the disentanglement is accomplished, then the Hilbert space containing the electronic and phonon variables will become separable. In a simpler language, the purpose of the transformations is to obtain a transformed Hamiltonian in which the electron and phonon variables will be decoupled. However, the Hamiltonian H is still not exactly soluble but it is quite reasonable to assume that after performing the four canonical transformations with generators (5), (6), (8) and (9), the electrons and phonons are weakly entangled in the Hamiltonian H and therefore the eigenstate of H can be approximately written as the product of an electronic state |ψ el � (hitherto unknown) and a phonon state | ph � which has to be judiciously chosen. Thus the total wave function corresponding to the full Hamiltonian H can be written as: www.nature.com/scientificreports/ As we have already mentioned, both adiabatic and anti-adiabatic effects have been incorporated in Eq. (10). In order to project out the effective electronic Hamiltonian, we have to eliminate the phonons from the problem. To accomplish that we need to take the expectation value of the Hamiltonian H in a suitable phonon state | ph � . Usually a zero-phonon state is chosen for the averaging phonon state | ph � . To make the calculation most accurate, we choose the averaging phonon state as: where ϕ n (x) is the n− th excited-state oscillator eigen function in 1D and the coefficients r n 's are to be obtained variationally. The idea is to start the calculation with M = 0 and then keep on increasing the value of M till the energy converges. It may be noted that the canonical transformations performed on Hamiltonian H by the generators given by Eqs. (5), (6), (8) and (9), followed by the averaging of the transformed Hamiltonian H with respect to the phonon state (11) is same as taking the expectation value of the Hamiltonian H with respect to the state which is the variational state for the phonon sub-system. We choose units in which = ω 0 = 1, ω 0 being the dispersionless phonon frequency. The effective electronic Hamiltonian can be defined as: where H k y and H l y are the Hermite polynomials of degree ' k ' and ' l ' , respectively, and The GS energy of the system described by the Hamiltonian H eff can be obtained exactly at half-filling with the help of the Bethe ansatz method first implemented by LW 11 . The LW solution has however been obtained for U eff > 0. We modify the solution to include the results for U eff ≤ 0 10, 38 . The GS energy per electron (ε) is finally obtained as www.nature.com/scientificreports/ where ε is finally minimized numerically with respect to the variational parameters to obtain the GS energy.

Average lattice displacement ( x i ), Entanglement Entropy (EE) and Local Spin Angular
Momentum ( L 0 ). In this sub-section we attempt to calculate a few interesting quantities such as the Average lattice displacement, Entanglement entropy and Local spin angular momentum L 0 for the system under consideration. These quantities give in general the physical properties of the system. The Entanglement entropy and Local spin angular momentum, in particular, provide information about the different phases that the system may possess. The numerical results for these quantities with respect to different physical parameters of the system will be presented in "Numerical results and discussion". The average lattice displacement is calculated with respect to the variational state |ψ ph � . The result is obtained as In order to study the role of quantum correlation in phase transition, we calculate the von Neumann Entanglement entropy for the 1D HH Hamiltonian. Considering a set of four available states |0�, |↑�,|↓� and | ↑↓� , the single-site entanglement entropy is calculated as: where ρ r is called the reduced density operator and is given by where Using the Hellmann-Feynman theorem, we get Thus all the occupation numbers can be calculated and the corresponding von Neumann entanglement entropy (E ϑ ) is evaluated.
The mean-square spin angular momentum per site ( L 0 ) can be defined as: where S i is the electron spin at site i , and S 2 n i↓ , n 2 i↑ = n i↑ , n 2 i↓ = n i↓ , we obtain, S 2 i = 3 1 − 2�n i↑ n i↓ � /4 and consequently where use has been made of Eq. (26). L 0 gives a measure of the spin magnetic moment and will be loosely referred to as the spin moment. For a completely un-correlated electron gas, we can write: �n i↑ n i↓ � = �n i↑ ��n i↓ � , and the average spin moment per site ( L 0 ) becomes equal to 0.375.

Numerical results and discussion
Ground State (GS) Energy. The single-site GS energy is determined by varying ε in the space of the variational parameters. The minimum value of ε gives the GS energy. The results are shown in Fig. 1. The harmonic TC results ( = 0&γ = 0) 10 are also plotted to show the effect of anharmonicity. It is clear that the lattice anharmonicity has a significant effect on the GS energy. We also show the results of the anharmonic case reported recently by Lavanya et al. (LSC) 25 . The LSC results have been obtained by using three canonical transformations. The present work uses an additional transformation namely the transformation given by Eq. (9). For = 0.1 and γ = 0.05 , we find that the new transformation (9) has only a marginal effect on the GS energy at small U . www.nature.com/scientificreports/ However, at large U, some observable effect is evident. More importantly, as we will show later, it has a more discernible effect on the phase diagram. Here we would like to mention an important point. We know that in a variational calculation, an error in the trial wave function of order δ will give an error of the order δ 2 in the energy. So, even a reasonable improvement in the wave function may not provide a sizable change in the energy. More importantly, here we are not so much interested in the energy but in the phase diagram. Our interest lies in examining whether with the improvement in the variational wave function, the width of the metallic phase broadens or shortens. If an improved variational calculation which obviously lowers the energy (by whatever amount) narrows the metallic phase, then the claim that an intermediate metallic phase exists at the crossover region of the CDW-SDW phases may be discounted. Our aim has been to show that with every improved variational calculation (which might improve the ground sate energy only marginally) the intermediate metallic Broadening of the metallic phase. The variations of the effective hopping integral (t eff ) and the effective onsite e-e interaction energy (U eff ) are respectively studied in Fig. 3a,b with respect to the onsite Coulomb  www.nature.com/scientificreports/ energy ( U) for the different strengths of the e-p coupling strength (α) . As expected, for α = 0 , t eff becomes equal to the bare Hubbard hopping parameter t and U eff becomes equal to the Hubbard U. As α increases, t eff decreases and with the increase in U , it gradually increases and saturates to the Hubbard value. At small values of U , the effective attractive onsite e-e interaction induced by the e-p interaction overcomes the repulsive Coulomb interaction U and therefore U eff becomes negative, i.e., attractive. The lattice is then unstable against the Peierls transition in which bound states of singlet bipolarons form on every alternate site leading to an insulating phase called the CDW state. On the contrary, when U is larger compared to α , the repulsive onsite electronic interaction wins, and the polarons cannot hop from one site to the other and consequently the GS of the system is given by the AFM mott-insulator state which is also known as SDW. Figure 3a shows that in the weak e-p interaction regime, the variation of t eff and U eff with U is continuous. But for higher values of α, discontinuous jumps occur in the behaviour of t eff and U eff . The discontinuous jump corresponds to a direct CDW-SDW transition. In order to examine the nature of the transition between the two phases observed in Fig. 3a,b, the quantity, ( dt eff /dU) is plotted in Fig. 4 with respect to U for α = 0.5, 0.8 and 1.0 . The double-peak structure in (dt eff /dU) is clearly evident. One can also observe that the peaks grow in height and shift towards larger values of U with increasing α. Furthermore, the new transformation used in the present calculation broadens the width between the two peaks in Fig. 4. Corresponding to the peak values of the (dt eff /dU) vs U plot, the phase diagram is drawn in the (α − U) plane. This is shown in Fig. 5. The intermediate phase satisfies the condition: 4t eff ≥ U eff , which is the signature of a metallic or a conducting phase. The metallic phase is flanked by the SDW phase on the left and the CDW phase on the right. The figure shows that the intermediate metallic phase appearing at the CDW-SDW cross-over region is now wider compared to that predicted by LSC 25 . It is important to emphasize that it is not important by how much the present modified variational wave function broadens the width of the intermediate  www.nature.com/scientificreports/ region, that the improved variational calculation widens the metallic phase is itself a result of great significance. The reason is simple. If a modified variational wave function predicts a narrower intermediate phase, it will have a disastrous effect on the prediction of the existence of a metallic region (MR) because one may then argue that the metallic phase may as well collapse if a more improved variational calculation is performed. That with every improved variational calculation, the metallic phase widens is indeed an encouraging result. In Fig. 6, we plot (dt eff /dU) with U for larger values of α. We find that the double peak structure almost disappears as α increases and one can observe from Fig. 6 that for α > 2 , only a single peak structure appears.
This indicates the absence of any intermediate phase for large α. Figure 7a,b illustrate respectively the behaviour of t eff and U eff as a function of α . As α → 0, t eff → t . Thus at low α, the system GS is in the SDW phase. As α increases, t eff gradually decreases and finally falls off to zero. Figure 7b tells us that corresponding U eff becomes maximally negative. This indicates the formation of massive singlet bipolarons giving rise to the CDW phase. Here also we see that for large U, SDW-CDW transition is again direct.
The phase boundaries of the (α − U) phase diagram can also be verified by plotting ( dt eff /dα) with respect to the e-p interaction coefficient α . This is shown in Fig. 7c. The peak-like structure is again obtained in the ( dt eff /dα) vs α plot for different U values and the peak width for the present result is broader than that of the LSC result. This again proves the intermediate metallic phase is widened.  www.nature.com/scientificreports/ Double occupancy (d) and EE (E υ ). In Fig. 8a,b, we plot respectively the double occupancy (d) and quantum EE (E ν ) as a function of α. The entanglement entropy (EE) gives a measure of the accessible states the system can have. Obviously then, the maximum in EE would correspond to a conducting state. It is observed that for certain combinations of α and U , E ν has maxima and for other values E ν becomes very small. Small values of EE correspond to insulating states. When the e-p interaction becomes strong compared to the e-e interaction, the electrons form pairs and the double occupancy parameter d reaches the maximum value of 0.5 driving the system to the CDW state. For d < 0.5 , the formation of a polaronic SDW state takes place. Similar behaviour is observed in Fig. 9a,b when we plot the double occupancy (d) and quantum EE (E ν ) with respect to U for different values of e-p interaction ( α).
In order to unravel the effect of e-p and e-e interaction simultaneously, we plot 3D graphs in Fig. 10. Figure 10a shows that between the SDW and CDW phases, there lies a region where the value of d neither corresponds to the SDW region with d = 0 nor to the CDW region with d = 0.5. Therefore, the effect of e-p and e-e interactions has been found simultaneously on d and E ν . This intermediate cross-over region corresponds to the metallic phase. In Fig. 10b, the peak of EE ( E ν ) lies over the metallic region (MR) in the (α − U) plane. Therefore, the peak denotes the metallic phase.

Mott-Hubbard (MH) criterion in 3D plots.
We have already emphasized that for a metallic state the bandwidth follows the criterion: 2zt eff ≥ U eff . In Fig. 11, we present a 3D representation of |U eff | and 4t eff with respect to U and α . The figure displays a region of (α, U) where the condition: 4t eff ≥ U eff is satisfied. This is the metallic phase. There are two other regions in the (α − U) plane where this condition is not satisfied and those are insulating phases. Among them the phase where U eff > 0 corresponds to the SDW phase and the one where U eff < 0 corresponds to the CDW state. In order to look into the metallic region by the Mott-Hubbard criterion more directly, 4t eff is plotted in Fig. 12, with respect to U eff and the region that satisfies the Mott-Hubbard condition is indicated by the dotted line. It is observed that the region satisfying the Mott-Hubbard criterion is

Conclusion
The nature of SDW-CDW transition has been studied in a 1D half-filled Holstein-Hubbard model with Gaussian phonon anharmonicity by improving the variational calculation of Lavanya et al. 10 . Using a number of unitary transformations performed in succession followed by a generalized many-phonon averaging an effective electronic Hamiltonian is obtained. The phonon-subsystem has been treated in a semi-exact way. The effective electronic Hamiltonian has been solved exactly using the Bethe ansatz technique to obtain the ground state