Renormalized phonon spectrum in the Su–Schrieffer–Heeger model

Motivated to understand phonon spectrum renormalization in the ground state of the half-filled Su–Schrieffer–Heeger model, we use the Born–Oppenheimer approximation together with the harmonic approximation to evaluate semi-analytically the all-to-all real-space ionic force constants generated through both linear and quadratic electron-phonon coupling. We then compute the renormalized phonon spectrum and the corresponding lattice zero-point energy (ZPE) as a function of the lattice dimerization. Crucially, the latter is included in the system’s total energy, and thus has a direct effect on the equilibrium dimerization. We find that inclusion of a small quadratic coupling leads to very significant changes in the predicted equilibrium dimerization, calling into question the use of the linear approximation for this model. We also argue that inclusion of the ZPE is key for systems with comparable lattice and electronic energies, and/or for finite size chains. Our method can be straightforwardly generalized to study similar problems in higher dimensions.


Introduction
The concept of Peierls instability in one-dimensional (1D) simple metals-the idea that 1D lattices are prone to dimerization in the presence of electron-phonon coupling-has been well established for over half a century. Notably, Su, Schrieffer and Heeger demonstrated it with the celebrated Su-Schrieffer-Heeger (SSH) Hamiltonian for polyacetylene [1], treating the lattice semiclassically.
In the SSH model, the effect of Peierls instability on the electronic spectrum is well-known. Yet the electron-phonon coupling responsible for the instability also affects the phonon spectrum, a fact that is less commonly discussed. The questions of how to properly account for lattice vibrations (phonons) and zero-point energy (ZPE) in systems with electron-phonon coupling, how strongly the phonon spectrum is modified, and how much the thus 'renormalized' ZPE influences the equilibrium lattice structure, remain open.
There are several available approaches to quantifying the effect of the electron-phonon interaction on the phonon spectrum. Early work by Ovchinnikov et al [2,3] focused on semiclassical techniques. Others, such as Nakahara and Maki [4] and Schulz [5], used Green's function methods and the random phase approximation to calculate the phonon spectrum renormalization, demonstrating phonon softening and phonon gap opening at the Brillouin zone edge (signatures of polaron formation as seen from effective mass renormalization have also been demonstrated using perturbative methods [6]). However, these studies did not explicitly sought to calculate the ground state dimerization self-consistently, by minimizing the total system energy including the ZPE; the equilibrium dimerization was obtained by minimizing only the electronic energy.
Specifically in the case of the SSH model, the questions highlighted above regarding the phonon spectrum and ZPE have a long history. Most early treatments neglected the quantum nature of the underlying crystal lattice: however, it soon became clear that the magnitude of the zero-point lattice fluctuations ∆u in the SSH model is comparable to the size of the dimerization u itself [7], raising the question whether the ground state is indeed dimerized. In the continuum SSH model these concerns were addressed with numerical renormalization group methods, which confirmed the presence of dimerization [4,7]-only mildly reduced by the lattice fluctuations-including for finite values of ionic mass M [8].
More recently, the dimerized ground state was also demonstrated in the lattice SSH model [9][10][11][12][13], with methods ranging from continuous-time quantum Monte Carlo to density matrix renormalization group. On the other hand, some recent studies combining ab initio density functional theory (DFT) calculations with the Fourier Grid Hamiltonian (FGH) method for the ionic Hamiltonian argued that the zero-point motion is sufficient to favor undimerized chains [14,15].
Separately, it was recognized that electron-electron interactions can play a significant role if added to the SSH model: for example, it was found that the phase diagram of the SSH model with Hubbard on-site repulsion can host a Mott-Hubbard insulator competing with the Peierls dimerized state [11,13]. It was also shown using the Gutzwiller approximation that electron-electron interactions generally tend to suppress the tendency towards dimerization by hardening the soft optical phonon and generally renormalizing the phonon spectrum to smaller energies [16], as well as (in the 2D case) potentially promoting nesting at other wavevectors [17]: at the same time, strong Hubbard interactions were seen to also encourage phase separation, just as in the Hubbard-Holstein model. However, few studies explicitly computed the phonon spectrum or investigated how it is affected by the lattice geometry.
Recently, a wide range of advanced variational methods have been developed to attack this problem, including the self-consistent ab initio lattice dynamics [18], time-dependent effective potential method [19], anharmonic lattice model from molecular dynamics [20], self-consistent phonon theory (SCPH) [21], and the stochastic self-consistent harmonic approximation [22] (the introduction of [22] offers a comprehensive recent review to the various methods). While these approaches offer incredible accuracy and are able to identify the equilibrium lattice configurations of a variety of challenging materials phases, the accuracy often comes at the expense of physical insight, and convergence can present difficulties.
In this paper, we study the phonon spectrum, ZPE and dimerization in the SSH model with linear and quadratic coupling. Combining the Born-Oppenheimer approximation and the harmonic approximation, we propose a direct, semi-analytical, self-consistent method of calculating the total system energy, accounting for both electronic and, crucially, lattice contributions. We show that the ZPE lattice contribution to the total energy is itself a function of dimerization. Minimizing the total energy allows us to self-consistently determine the equilibrium dimerization, and to evaluate when it is important to include quadratic coupling and/or the ZPE in the calculation. We acknowledge that a similar approach to calculating the renormalized phonon spectrum has appeared elsewhere [23], however to our knowledge we are the first to add quadratic coupling and also the first to include the resulting ZPE into total energy calculations and thus let it affect, self-consistently, the equilibrium lattice structure. Our method is demonstrated here for the SSH Hamiltonian, but it can be straightforwardly generalized to higher dimensions and multiple phonon and/or electron bands. The work presented here ignores electron-electron interactions. These can be introduced in the formalism described here at a mean-field level (we discuss possible effects of correlations on our results in more detail in section 5).
Our method can be straightforwardly integrated with existing DFT calculations, which presently do not account for ZPE during geometry optimization. Unlike more complex variational methods, it would only require adding ZPE-a quantity DFT methods already routinely evaluate with supercell calculations of the phonon spectrum-to the total energy during optimization; and this should not significantly affect convergence behavior of geometry optimization. The proposed method could be helpful in improving geometry optimization in systems with strong electron-phonon coupling-especially those with small effective bandwidth W and relatively large phonon frequencies Ω, such as in organic photovoltaics.
The work is organized as follows: in section 2, we describe our effective model Hamiltonian. In section 3 we describe the proposed method (with certain technical details kept in the appendices for clarity). Our results are discussed in section 4, and we conclude with some final remarks in section 5.

The model
We start from the Hamiltonian of a half-filled, undistorted, 1D chain with lattice constant a and one orbital per site, and include longitudinal lattice vibrations: Here c † nσ is a creation operator for an electron on site n with spin σ, t n,n+1 are nearest-neighbor hopping integrals, K is the stiffness of the bonds between neighboring ions, andû n =R n − na are the (operator) deviations of the ions from their undistorted equilibrium positions. This Hamiltonian, compared to an actual conjugated polymer chain it is often taken to model, has many approximations. First, it neglects electron-electron interactions. We note that in materials where these are important (such as polyacetylene [2,11,24]) one can incorporate them into the formalism we describe here at the Hartree-Fock level. Second, we include only nearest-neighbor hopping and only nearest-neighbor effective interactions between ions: however, both these approximations can be relaxed trivially. Third, as customary, we neglect the geometric complexity of real chain polymers, including the zig-zag structure of trans-polyacetylene, the presence of other ligands, out-of-plane bending and torsion, defects, etc [3]. All these geometric simplifications can be relaxed when a quantitative description of a specific material is sought, however the simple model of equation (1) suffices for our current purposes of illustrating the proposed method.
Hamiltonian (1) must be supplemented by a rule for the dependence of the hoppings t n,n+1 on the bond lengthsR n+1 −R n . It is customary to use the SSH linear approximation for the electron-phonon coupling: t n,n+1 = t − α(u n+1 − u n ), however we note that this is not consistent with the harmonic expansion of the Born-Oppenheimer potential energy surface, which is carried out to second order. For internal consistency, we expand the hopping amplitude also to second order, resulting in the quadratic coupling: We therefore consider the quadratic SSH model: In the following, we investigate a range of possible (t, α, β, K, M) choices, beginning with those relevant for polyacetylene. In all cases we demand that β||û n || ≪ α, α||û n || ≪ t and ||û n || ≪ a, so that the harmonic approximation remains valid.

The method
To simplify notation, the formulae below are for the linear coupling model (β = 0). Those corresponding to finite β are listed in appendix C.

Born-Oppenheimer decomposition
We treat our Hamiltonian within the well-established Born-Oppenheimer approximation, which is justified given the small ratio of the electron and ion masses. This allows us to sequentially solve the separate electronic and ionic problems, with the total electronic energy acting as a potential for the ions-that is, defining a Born-Oppenheimer potential energy surface. Specifically, the Hamiltonian of equation (2) where ⟨û n ⟩ ≡ ⟨R n ⟩ − na. For completeness, we sketch our Born-Oppenheimer procedure in appendix A.

The electronic problem
Keeping in line with previous work going back to Peierls [25], we anticipate a dimerization of the 1D chain. We single out this dominant static distortion from smaller lattice fluctuations by adopting the ansatẑ u n = (−1) n u +x n , ||x n || ≪ u.
We will check the self-consistency of this approximation a posteriori, by minimizing the total energy with respect to this variational parameter which sets the new equilibrium positions to R 0 n = na + (−1) n u.
Combining the static distortion term with the hopping term, we recover the standard electronic SSH model plus a perturbation (we denote ⟨x n ⟩ = x n ) The first term involves the usual conduction and valence band operators: where and for a chain with N sites (N/2 dimerized unit cells). Furthermore, the band energies are defined by This bandstructure is the usual SSH result, showing a gapped electronic spectrum whose size is controlled by the static dimerization u. The equilibrium contributions (x n = 0) to the ionic Hamiltonian, at half-filling, are: The next step is to find the dependence of the electronic energy on all x n = R n − R 0 n , so that we can then use E e ({R n }) to solve the ionic problem.

Harmonic approximation
To allow further analytical progress, we use the harmonic approximation to represent a region of the Born-Oppenheimer potential energy surface E e ({R n }). The zeroth order term is the (electronic) energy associated with the equilibrium lattice positions, E e ({R 0 n }), listed above. The first order term vanishes when we are expanding the electronic energy around the new equilibrium lattice positions. To second order in ionic displacements, we are left with: The second-order derivatives are commonly known as the dynamical matrix. We now proceed to calculate them.
Substituting equation (14) into equation (15), we find An immediate simplification can be made at this stage. The contributionV i−i = K/2 n ⟨û n+1 −û n ⟩ 2 does not depend explicitly on the electronic wavefunctions. This implies E e = ⟨Ĥ e ⟩ = ⟨Ĥ −V i−i ⟩+ K/2 n (u n+1 − u n ) 2 . Using the definition of x n = u n − (−1) n u, we arrive at To evaluate the derivatives, we use perturbation theory to calculate the expectation value F e ≡ ⟨Ĥ e −V i−i ⟩ to second order in x n . The only non-trivial dependence comes from theÛ el−ph electron-phonon coupling term in equation (6). To set the stage for the perturbative calculation, we rewriteÛ el−ph = n x nfn where the electronic operatorsf n are expressed in terms of the ν k,σ and χ k,σ conductance and valence band operators. Their form depends on whether n is even or odd. Specifically, using the shorthand sin(ka) ≡ s k , we find These expressions are cumbersome, but for a ground state calculation of a half-filled model we only need a single entry, as shown in appendix B.
The perturbative expansion is: e + . . . where the corrections are given by the usual quantum-mechanical expressions, namely Here |Ψ 0 ⟩ is the electronic Slater-determinant ground state of the half-filled SSH model, consisting of a full valence band and an empty conduction band. This leads to: Now we consider several cases depending on the even/odd character of n, m. The details are relegated to appendix B. The relevant results are the expressions for the dynamical matrix entries δK nm listed below (the factor of 2 is from the sum over spins).
Note that the sums over k, q run over the reduced Brillouin zone [−π/2a, π/2a] due to the dimerization ansatz.
While the expressions for the even-even and odd-odd sites turn out to be identical, there are subtle differences for the cross-terms that lift the phonon spectrum degeneracy at the Brillouin zone edge and split the optical phonon band from the acoustic band. In particular, there are two inequivalent terms δK 2n,2n+1 ≡ Y + 1 , δK 2n+1,2n+2 ≡ Y − 1 , which correspond to the fact that there are even and odd electronic operators that give rise to conduction and valence bands, It is remarkable that even though we started with only nearest-neighbor atomic force constants, through interactions with the extended electronic states we now have all-to-all force constants emerging. Because the equations (23)-(25) are closed-form expressions, we can easily calculate their values using computer integration, after going to the thermodynamic limit k → L 2π´d k (L = Na = (N/2) · 2a). The final step is to diagonalize the ionic Hamiltonian defined by these dynamical matrix elements.

The renormalized phonon spectrum
We adopt the following short-hand notation for the dynamical matrix elements: the onsite value is Z 0 ≡ 2K + δK n,n ; the two-sites away value is Z 2 ≡ δK n,n+2 ; the even and odd bond values are ; and the remaining longer-range values are Renaming the ionic operators asp n ,P n ,x n ,X n for the even and odd lattice sites of the nth unit cell, respectively (our unit cell is taken to start at the even site) we rewrite: Because the ionic motion is influenced by the dimerization and we find an optical phonon branch (in the folded Brillouin zone), even though we started with only a bare acoustic branch (in the full Brillouin zone). Longer-range terms do not affect the size of the unit cell, but they do further renormalize the phonon spectrum.
The resulting phonon spectrum (see appendix D for details on the diagonalization of equation (26)) is: This expression for the phonon spectrum demonstrates the appearance of an optical branch separated from the acoustic one when u ̸ = 0. Note that in the absence of electron-phonon coupling we have Z 0 = 2K, Z δ⩾2 = 0, Y ± 1 = −K, from which we readily recover the original acoustic spectrum folded in the smaller Brillouin zone: For a finite electron-phonon coupling, the number of Z δ to be included in the calculation, defined by the cut-off |δ| < δ max , depends on how quickly they decay. We select the cutoff δ max so as to ensure the convergence of the phonon spectrum and ground state energies. In practice, much higher δ max is needed for smaller u, because the energy difference denominator for small u is small, leading to difficulties in the integrand.
In the appropriate basis of phonon operators b jq , b † jq , we thus have: To complete the calculation and find the ground state, we need to find the u GS that minimizes the total energy of the system Furthermore, this value of u GS must lead to a real phonon spectrum ω jq , signifying that this is indeed a stable equilibrium point. Because the renormalization of the phonon spectrum ω jq depends on the value of u, the lattice ZPE contributes to determining the equilibrium value u GS , unlike in the traditional SSH approach where only the electronic contribution (second term) is considered. This is one key result of our work.

Results and discussion
We now discuss the impact of including the ZPE in total energy, as well as the effect of adding the quadratic coupling β ̸ = 0.

The renormalized phonon spectrum
When not otherwise stated, we adopt the typical polyacetylene parameters t = 2. First, we study the phonon spectrum of the model without any dimerization, u = 0. The traditional semiclassical calculation neglects any back-reaction of the electrons on the phonons and predicts an undisturbed, folded acoustic phonon band (equation (28)), shown in panel (a) of figure 1. By contrast, for the undimerized (u = 0) chain our method finds a massive Kohn anomaly in the optical phonon branch at q = 0 (the folded q = 2k F point), as expected for a system that is unstable to dimerization, see panel (b) of figure 1. The anomaly is so strong that it results in the phonon frequencies becoming imaginary near q = 0: for such values we take their magnitude and plot them as negative on the graph, in accordance with the usual convention for calculations near an unstable lattice configuration. Note that the solid line corresponds to β = 0 (linear coupling model), while the dashed line is for a finite β (quadratic coupling model). We follow this convention throughout.
An interesting observation is that a finite β renormalizes the phonon spectrum in figure 1(a) even without dimerization. The phonon bandwidth is increased slightly due to the constant term in the first-order perturbative expansion (see equation (C17)). In particular, for u = 0 we find δK nn ∼ 8β/π, δK n,n+1 ∼ −4β/π, boosting the bare spring constant K. This makes sense: even in the absence of dimerization, the quadratic coupling provides some additional 'glue' to the bonds.
The Kohn anomaly seen in panel (b) of figure 1 arises out of the many new force constants generated by the electron-phonon coupling, shown in figure 2. The case u = 0 is a difficult limit for our calculation due to the strong singularity in the denominators of equations (23)-(25) which makes the force constant corrections Z δ decay very slowly with δ, so that a δ max ∼ 160 is needed for convergence (the biggest obstacle to convergence is the q = 0 point of the acoustic spectrum which gathers in-phase contributions from all the spring constants). This slow convergence with neighbor number is not unreasonable for this case and has been noted by other researchers [26]. The force constants decrease as the distortion moves away from the singularity at u = 0, see panel (a) of figure 2; this improves the convergence drastically. Panel (b) shows that increasing the electron-phonon coupling strength increases the magnitude of the force constants, but does not affect their range significantly.
In general, the effect of finite β is to lower the force constants, with the most sizable effects on the first two force constants. This is because there are two kinds of β-generated corrections: those from the first order of perturbation theory (equation (C17)), which are relatively large but only nonzero for same-site and nearest neighbor spring constants, i.e. δ = 0, 1; and those ∼u 2 from second order of perturbation theory (equations (C29)-(C32)), which are much smaller because u ∼ 0.1 or smaller. The apparent increase of force constants with β for the blue curve with α = 1 eV Å −1 is only because we plot the absolute value Z: in reality, the force constants flip sign and become more negative, as the β corrections are subtracted off of the nearly-zero α ones.
Panel (c) of figure 1 shows the phonon spectrum corresponding to the equilibrium dimerization u GS (discussed below). As expected for a dimerized chain, a proper optical phonon branch emerges, separated from the acoustic branch by a gap of ∼10 meV at k F = π/a. The Kohn anomaly is lifted, its only remnant being the q = 0 softening of the optical phonon mode, by ∼60 meV if β = 0, and less for finite β. We note that now all phonon frequencies are real and positive, showing that this dimerized configuration represents a stable equilibrium.
The quadratic coupling β moderately increases the phonon gap: for β = 2 eV Å −1 , the gap changes from 10 meV to 12 meV-a 20% increase. The phonon gap is driven by the difference between the two inequivalent force constants, Y + − Y − , for the long-short bonds. This increase can be attributed in roughly equal parts to (a) the effect of β on force constants, and (b) the effect of β on increasing the equilibrium dimerization u gs . As discussed above, the former effect is because through equations (C29)-(C32), β contributes directly to Y + and Y − , and introduces further Y + δ /Y − δ terms which also contribute to the gap. The latter effect is more subtle: the quadratic SSH coupling induces a dimerization dependence into the effective hopping amplitude, t → t + 2βu 2 (we discuss this in more detail in section 4.2). This dependence widens the electronic bandwidth, favoring a larger equilibrium u gs . In turn, the larger u gs enhances the Y + /Y − asymmetry, so that the gap increases (roughly) linearly with u.
The quadratic coupling also has a significant effect on the optical phonon softening, which decreases by about 30 meV-a considerable 50% decrease. This is also primarily due to the increase in u gs discussed above: a larger dimerization supresses the softening because the latter arises from the singularity in the denominators of equations (23)-(25) (Kohn anomaly).

Comparison to other phonon spectrum calculations for polyacetylene
Before continuing the analysis, it is useful to compare our results to other existing predictions for the phonon spectrum, to gauge the performance of our method.    [27], and to a previous analytical calculation mentioned in their introduction [23]. It is unclear to us what parameters were used by these authors to generate the analytical results, hence we use the typical polyacetylene parameters listed above (with β = 0), and we scale the phonon spectra ω jq by the bare optical phonon frequency ω Q = 2K/M * . For the study by Miao et al, we estimate their ω Q = 1420 cm −1 from the available results.
We find reasonable agreement, despite the uncertainty about their parameter values. For the optical branch our method produces a slightly larger phonon softening and a smaller phonon gap. For the acoustic branch the agreement is perfect between the two analytical methods, while DFT predicts a strongly flattened  [28] (black dashed line). We use the same parameters t = 6.15 eV, α = 7.6 eV Å −1 , K = 81 eV Å −2 , M = 1.99 × 10 −26 kg but set u = 0.044 Å to reproduce the electronic band gap of 2.7 eV used in [28]. Convergence was achieved with δmax = 40. acoustic branch, similar to what is seen in other DFT studies of polyacetylene [29]. Miao et al hypothesize that this is because the SSH model treats only longitudinal phonons, whereas in DFT studies there is significant coupling between the long wavelength acoustic mode and transverse phonons. More work is needed to clarify this point.

Comparison to other spectrum calculations for carbyne
Carbyne-a pure carbon chain with alternating single and triple bonds-is also convenient for a comparison, given the similarly of its crystal structure to polyacetylene. However, carbyne has two degenerate electron orbitals, p y and p z , that can host delocalized π-bond electrons, compared to polyacetylene's single p z orbital. A simple way to address this difference, widely used in other semiempirical approaches, is to consider the orbitals as fully equivalent and account for them by adding an overall factor of 2 in equations (23)- (25).
In [28], carbyne was studied with a semi-empirical bond-bond polarization approach to construct a force field that gave good agreement with experimental measurements of Raman spectra. By adjusting model parameters to match the observed electronic band gap and Raman excitation frequencies, the authors adopted t = 6.15 eV, α = 7.6 eV Å −1 , K = 81 eV Å −2 , M = 1.99 × 10 −26 kg and found the dimerization u GS = 0.088 Å. Note that their approach did not explicitly minimize the ground state energy against dimerization or other crystal lattice parameters, and thus does not provide a theoretical origin for the ground state dimerization.
For these parameters, our method predicts an undistorted ground state, u GS = 0, because the spring constant is very stiff and the lattice deformation energy ∼Ku 2 grows too quickly. However, it is known from numerous previous experimental and ab initio studies that carbyne (and finite-length polyynes) are indeed dimerized. The fact that our method predicts an undimerized ground state may be due to the exclusion of electron-electron interactions from the model, whose importance was long ago recognized by Ovchinnikov [3].
To still be able to compare phonon spectra, we forego energy minimization and simply set our u = 0.044 Å, so as to reproduce the electronic band gap of 2.7 eV. Figure 4 shows a comparison of our corresponding phonon spectrum against results reported in figure 3 of [28]. We find strong agreement especially for the optical phonon branch, however the phonon gap is somewhat smaller in our approach.
Overall, we conclude that our quasi-analytical method works reasonably well in predicting the renormalized phonon spectrum.

The equilibrium dimerization
After checking the renormalized phonon spectrum generated by our method, we return to our calculation, and plot in figure 5 the total energy per site as a function of u (still for polyacetylene-like parameters). The curves labeled 'SSH' are obtained by the standard SSH calculation (disregarding the renormalization of the phonons and their ZPE), while our results are labeled 'BO+Har' .
For the linear coupling model (β = 0), inclusion of the renormalized ZPE contribution reduces the ground state dimerization value from 0.041 Å to 0.039 Å (vertical dot-dashed lines). This amounts to a 5% change-a small, but noticeable effect. Either with or without the ZPE, 2αu GS ≈ 0.33 eV ≪ t = 2.5 eV, so the linear approximation appears to be justified. Nevertheless, the dashed lines show that adding even a small The circles indicate the corresponding ground states, with vertical dot-dash lines as guides to the eye. Solid/dashed lines correspond to β = 0 and β = 2 eV Å −2 , respectively. The orange triangle is the total energy of the system in the absence of electron-phonon coupling, but accounting for ZPE. The convergence parameter is δmax = 160, the other parameters are as in figure 1. See text for more details. quadratic coupling has a very significant effect, changing the dimerization to u GS ≈ 0.055 Å. Note that even for this larger dimerization, 2βu GS ≈ 0.2 eV Å −1 ≪ α = 4.16 eV Å −1 , so a priori one would not expect the quadratic coupling to have such a significant effect.
This effect is primarily driven by the non-linear dependence on u introduced into the electronic bandwidth by the β term: t →t = t + 2βu 2 (see appendix C). This results in a smaller α/t ratio so naively one may expect a decrease in uGS. (Note that the observed decrease in the optical phonon softening near q = 0 at finite β, see figure 1(c), is consistent with the smaller α/t). Instead, the dependence of the bandwidth on u 2 drives the observed increase in u GS in the quadratic model. This shows that the effects of the quadratic coupling cannot be accounted for by a re-scaling of the parameters of the linear model. This significant dependence of u GS on β when βu GS ≪ α is one of our main results. It is very important because it is customary to choose α so as to obtain u GS equal to the value measured experimentally. The strong dependence of u GS on β shows that this choice of parameters is very dubious. Instead, one should include the quadratic coupling in the model and use another experimental data point (for instance, the softening of the optical mode) to find more appropriate values for the electron-phonon couplings α, β.

To dimerize or not to dimerize
The results in figure 5 also show that the typical ZPE is clearly larger than the Peierls dimerization barrier (the energy difference between the dimerized and undimerized structure). This has long been known, as is the fact that the magnitude of the zero-point fluctuations ∆u in the SSH model is comparable to the size of the dimerization u itself [7]. These observations called into question the validity of neglecting the quantum nature of the lattice in earlier work, and raised concerns whether the model's true ground state was dimerized at all.
These concerns were addressed for both continuum models, e.g. [4,7,8], and for lattice models, e.g. [9,[11][12][13], as mentioned in the Introduction. By using a variety of approaches, these works demonstrated that the half-filled, spinful, Hubbard-SSH model with quantum acoustic phonons admits two ground states: a Mott-Hubbard insulator and a Peierls (dimerized) insulator, but no metallic phase.
And yet in recent years several works have questioned this conclusion. Hudson and Allis [14] argue that zero-point motion is sufficiently strong to favor undimerized chains, invoking a two-step calculation: ab initio DFT calculations for the electronic energy, and a single-particle FGH method for the ionic Hamiltonian (with the ionic potential reconstructed from the DFT electronic energy). Their findings appear to be supported by other researchers [15]. Interestingly, recent experimental transport studies on short (N = 10 to 20) carbyne chains [30,31] found that un-strained chains are metallic. Meanwhile, polyacetylene's semiconducting behavior and Raman spectra [32][33][34] can be attributed to electron-electron interactions without recourse to a Peierls instability [35], as suggested by Ovchinnikov et al in the early 1980s [3,36].
We use our results to address the possible origin of (some of) this disagreement. The renormalized phonon dispersion in figure 1(c) shows that the modes most strongly affected by the electron-phonon coupling are those near q = 0. If only the q = 0 mode is considered, i.e. a Γ-point approximation is employed like in [14], the variation of the ZPE with dimerization u is overestimated. To exemplify this, in figure 6 we repeat figure 5 but assuming that the whole ZPE comes from the q = 0 modes scaled by the system size N/2 for a proper comparison (only β = 0 results are shown). Clearly, this additional approximation significantly Figure 6. Same as figure 5, but the Γ-point approximation is used to generate the ZPE contribution (i.e. only the energy of the q = 0 phonon modes is counted in the BO+Har approach, scaled by the system size). This additional approximation predicts a significantly reduced ground state dimerization-if it were not for the fact that the phonon spectrum becomes imaginary closer to u = 0 (energy not shown), the undimerized chain would be the predicted ground state. The two curves touch when the optical frequency vanishes at q = 0. changes the results. If it were not for the fact that the phonon spectrum becomes unstable near u = 0, the predicted ground state would now be undimerized. This emphasizes the importance of using the entire phonon spectrum when calculating the ZPE, and shows that employing the Γ approximation can overestimate the ZPE's contribution to the point of predicting an incorrect lattice structure.

Isotope effect
The results presented so far might suggest that calculating the renormalization of the phonon spectrum and including the ZPE in the total energy is an unwarranted complication, given that u GS is little affected. This is to be expected for polyacetylene, given the order of magnitude difference between the electronic and lattice energy scales. However, there may be materials where the energy scales are not so disparate-for instance, in a hypothetical hydrogen chain, or if the electron mass is strongly renormalized through interactions. In such cases, including the ZPE in the calculations should lead to noticeable changes to the predicted u GS .
To exemplify this point, we lower the mass parameter M as a simple way to boost the lattice ZPE energy scale. The results are shown in figure 7, where we compare the total energy curves for cases where the ionic mass M * = 0.1, 0.01M, M being the polyacetylene value; the other parameters are kept unchanged. Indeed, when M * decreases by two orders of magnitude (so that the characteristic phonon energy increases by one order of magnitude and becomes comparable to the electronic energy scale), the addition of the ZPE predicts a much less dimerized GS. We caution against a literal interpretation of the energies close to u = 0, where the phonon spectrum is unstable (and results are not shown). For contrast, we note that the standard SSH calculation remains as in figure 5, predicting u GS ≈ 0.4Å for any value of M * . As before, adding the quadratic coupling increases the equilibrium dimerization compared to that of the linear model, all else being equal (dashed lines).
Of course, the results in panel (b) are unphysical as at most M * ≈ 0.07M ≈ m p (for a hypothetical hydrogen 1D chain). However, the phonon energy scale can be increased not just by lowering M, but also by stiffening the spring constant K. The important point is that if K/M ∼ t, the ZPE contribution varies significantly as a function of the dimerization u, thus contributing nontrivially to setting u GS .
Comparing figures 5 and 7, we note that the finite β curve has lower energy in figure 5, but its energy seems to increase with decreasing mass parameter M * . This can be explained by changes to the magnitude of the zero-point lattice energy. Larger β leads to a larger ZPE contribution, mainly because the phonon bandwidth is mildly boosted by the first-order perturbative term (as we argued in section 4.1). This ZPE contribution varies relatively weakly with u and so generally just shifts the entire energy curve upward. As the mass is decreased in figure 7, the magnitude of the ZPE contribution increases: and it increases more for the case with finite β, leading to the change in relative positions of the energies.

Finite length rings
If we calculate the discrete sums in equations (23)-(25) instead of taking the thermodynamic limit N → ∞, we can study finite-size effects for closed rings with an even number of atoms. Here we fix δ max = N/2 where N is the total number of atoms. The physical point of comparison here is something like a benzene ring A significant reduction of the ionic mass boosts the characteristic phonon energy, and thus the importance of ZPE. As before, dashed lines represent results with β = 2 eV Å −2 . A significantly smaller lattice dimerization is favored once M * becomes sufficiently small: it would be zero were it not for the phonon spectrum becoming unstable-this happens near u = 0, and data for such u is not shown.
(N = 6). Cyclic polyenes are known to have a transition from all-equal (undimerized) to conjugated bonds at a critical size, typically pegged at N ∼ 30-50. In what follows, we always take the lowest energy state as the ground state, subject to the phonon spectrum being stable.
In figure 8 we plot the ground state dimerization u GS (panel (a)), and ground state energy E GS (panel (b)) as a function of chain length N using the polyacetylene model parameters and β = 0, calculated using the SSH and BO+Har approaches. For comparison, the corresponding infinite N limits are given by the dotted and dashed lines, respectively, and are denoted 'SSH inf ' and 'BO+Har inf ' in the legend. As discussed above, there is a small difference between the dimerization predicted in the infinite N SSH and BO+Har approaches, see figure 5. We split our dataset into chains whose length N is such that N/2 is even (blue triangles), while the rest are in the odd group (orange circles). We also show the SSH finite N results (red diamonds for both even and odd groups) in panel (a), to gauge the change in the ground state dimerization driven by the ZPE.
Convergence to the thermodynamic results is achieved for N ∼ 40. Interestingly, even group chains converge to the asymptotic limit from above (bigger dimerization, higher energy at small N), while odd group chains approach from below (smaller dimerization, lower energy at small N). This is a finite size effect, because in the even group, a pair of momentum points (one point for each of optical/acoustic, and valence/conduction bands) at the edge of the Brillouin zone is included, and it is here where the effect of dimerization on the electronic bands (and thus on the total energy) is strongest. For short rings with only a few allowed momenta, there is a big difference if the edge value is included or not.
This alternation is indeed seen in both the SSH and our approach. In organic chemistry, this is known as Hückel's '4n/4n+2' rule: the odd group has all electrons in 'bonding orbitals' , whereas the even group has a pair of electrons in the non-bonding orbital with k = ±π/2a. Figure 8(a) shows that for small N ≲ 20 values, the inclusion of the ZPE in the total energy can result in much more significant changes to the predicted u GS than is the case in the thermodynamic limit; this points out to the importance of the ZPE contribution in such cases.
We remark that open chains can also be studied with our method, but this leads to dimerizations that vary along the chain and thus to various complications whose study is left for future work.  (panel (a)) and ground state energy EGS (panel (b)) for finite-N chains with N/2 odd (orange circles) and even (blue triangles), as well as the SSH finite N results (red diamonds in panel (a)). The lines connecting the circles and triangles are guides to the eye. The SSH and BO+Har infinite N limits are given by the dotted and dashed lines, respectively, and are denoted 'SSH inf ' and 'BO+Har inf ' in the legend. The even/odd behavior is explained by the momentum points that are allowed for each N, and whether they fall on the BZ edge or not. Finite N results approach corresponding infinite N results for both the SSH and the BO+Har approach. The effect of ZPE on ground state dimerization is visible in panel (a), where the SSH predicts a higher uGS than the BO+Har for N > 14.

Conclusions
In this paper, we used the Born-Oppenheimer approximation together with a harmonic approximation of the Born-Oppenheimer potential energy surface, to propose a semi-analytical approximation that self-consistently accounts for the contribution of the zero-point lattice energy to the ground-state energy optimization. For consistency with the harmonic approximation, we also included quadratic terms in the electron-phonon coupling.
In contrast to variational techniques discussed in the introduction, in our approach the force constants generated through the electron-phonon coupling, are calculated directly instead of being taken to be variational parameters adjusted to minimize the energy; the static lattice dimerization u is our only variational parameter.
For a model with typical polyacetylene parameters, including the ZPE in the total energy led to a modest change of the average dimerization u GS , of around 5%. However, we argued that for systems whose lattice and electronic energies are more comparable, and/or for finite chains, including the ZPE in the total energy changes the predicted u GS much more significantly.
Adding even a very small quadratic coupling changes u GS for (otherwise) polyacetylene parameters by around 30%, and similarly big changes are observed for other parameters. This strongly calls into question the validity of using a linear approximation for the electron-phonon coupling in these models.
For polyacetylene-like model parameters, we find that the zero-point phonon energy contribution does not destroy the dimerized ground state of polyacetylene for an infinite chain. This agrees with several previous results and disagrees with others. Relevant to the latter, we point out the potential downfall of using the Γ-point approximation, as it overestimates the ZPE dependence on u and can thus lead to incorrect conclusions.
The key approximations we made-insofar as the Hamiltonian we use is meant to be a model for polyacetylene and other 1D organic chains-are to treat the ionic BO energy within the harmonic approximation, and to ignore electron-electron interactions. The former is fairly well-established: the utility of the harmonic approximation is that it allows calculating the force constants using perturbation theory, and thus an analytical expression for the phonon spectrum is available. In principle one can relax this approximation by including cubic and higher order anharmonic terms, whose coefficients can be calculated with the appropriate higher-order perturbation theory. However, additional approximations (e.g. mean-field) are needed to then ascertain the effect of these anharmonic terms on the phonon spectrum. Given its small lattice energy scale, we believe that for polyacetylene such a calculation is not warranted, but it might be relevant for a system with larger phonon frequencies Ω ∼ K/M.
We have completely ignored electron-electron interactions, which could be important for systems like polyacetylene. If added, such interactions would likely increase dimerization in the ground state [2,24] and further renormalize the phonon spectrum by reducing the phonon frequency at low and intermediate momenta [16]; on the other hand, such interactions are known to open the door to other ground states such as charge density wave order [11,12], phase separation [16], and even nesting instabilities at other wavevectors (in the 2D case) [17]. Weaker electron-electron interactions can be added into our formalism at the mean-field level, and it would be interesting to study the competition between the phonon spectrum renormalization effects of both electron-electron and electron-phonon interactions.
The remaining limitation of this study is that it only investigates T = 0 ground-state behavior. A finite-T generalization is in principle straightforward and would allow us to investigate the effects of the quadratic coupling on the transition temperature to an undimerized state [37]. We note that unusual behavior with T, attributed to quantum anharmonic effects, has been recently observed in carbyne chains by Romanin et al [38], wherein, counter-intuitively, dimerization actually increases with temperature: given that β deepens the potential energy well, it is possible that some of the behavior observed in [38] (obtained in ab initio), while unattainable with linear SSH coupling, could be manifested with quadratic coupling. We leave such calculations for future work.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).
The last term in equation (A5) is the term preventing the Born-Oppenheimer ansatz from being exact. An argument due to Slater [40] demonstrates this term to be smaller than the rest by a factor of m e /M: note that the last term is ℏ 2 /(2M)∂ 2 ψ/∂R 2 n ≈ ℏ 2 /(2M)∂ 2 ψ/∂r 2 n . The approximate equality holds because the wavefunction depends on the differences r n − R n of the coordinates and thus both derivatives remain of the same order. However, E e will have in it the kinetic energy of the electrons, which is order ℏ 2 /(2m e )∂ 2 ψ/∂r 2 n , with the difference being precisely a factor of m e /M. Therefore ignoring the last term in equation (A5) is an accurate approximation for m e ≪ M.

Appendix B. Perturbative calculation details
In this appendix we describe the details of getting from the general perturbative expressions in equation (22) to the results reported for δK nm in equations (23)- (25). First, we catalog the 'halves' of the perturbative expressions for the even-odd cases.
Even m: We arrived at this expression as follows: since on the right we have the ground state |Ψ 0 ⟩, the only term from the matrix product in equation (18) that can survive at half-filling is the χ † ν combination. As for the energy denominator, Given the reasoning above, we can immediately write down the odd version of the term Odd m: The corresponding left-hand sides of the overall perturbation expression are no different: now we keep the ν † χ combination, so that it does not annihilate the state on the left Even m: Connecting the left and right halves of the expressions forces q ′ = k, k ′ = q. Re-writing it for clarity The only thing that changes for odd m is the prefactor: Odd m: Now that we have all of the relevant expressions, we can begin to put them together. The first of equations (23)-(25) is derived as follows The other equations are obtained similarly.
All the definitions are as in appendix B, with the caveat that t → t + 2βu 2 .
As before, we are evaluating energy corrections where the terms are given by Introduce the notation Unlike in the conventional SSH model, the linear term from the perturbative expansion will have a non-zero contribution due toB el−ph in the Hamiltonian. That is the new calculation: the second-order term will be simply appended with the new coefficients α ± = α ± 2βu. There is no contribution fromB el−ph at second order because those will lead to terms cubic or quartic in atomic displacements, which will be set to zero in the harmonic approximation for the Born-Oppenheimer energy surface.
Start with the first-order term coming fromB el−ph . For |n − m| > 1, the first-order correction is zero, as is clear from the following calculation: Then: δK (1) 2n,2n = −β T 2n,2n+1 + T 2n−1,2n To finish the evaluation, we must compute the expectation values above. To do so, we should express the old c operators in terms of ν, χ (at half-filling, we will take the un-perturbed ground state to be the usual Fermi sea |Ψ 0 ⟩ = |k|<π/2a,σ ν † kσ |0⟩) Similarly, for the other expectation value, we find at half-filling Substituting into equation (C11), we find δK (1) 2n,2n = − 2β N/2 kqσ δ kq cos(ka)β k α q + cos(qa)β * q α k = −2 If n is odd, the matrix elements swap places, but the end result remains δK (1) 2n+1,2n+1 = −2 We may combine the two terms by writing It makes sense that the terms are the same for even and odd sites, as these are on-site corrections, and any site connects to a pair of short/long bonds. After more similar calculations, we find δK (1) nm = δ nm − 4β N/2 k ϵ k cos(ka) E k + (δ n+1,m + δ n−1,m ) 2β N/2 k ϵ k cos(ka) + (−1) n u∆ k sin(ka) E k . (C17) Evidently, there are two inequivalent spring constants, and thus an optical and an acoustic branch will emerge for a nonzero α. Now for the contribution fromÂ el−ph . As in the linear SSH model case, re-write it in terms of ν and χ operators,Â el−ph ≡ n x nĝn (C18) where we defined α ± (n) = α ± (−1) n 2βu, andĝ n = −α + (n)T n,n+1 + α − (n)T n,n−1 . Then α ± (2n) = α ± , andĝ 2n = (C20) The final step is to convert the operators to the ν, χ basis. Defining g k = 2i α sin(ka) − 4βu cos(ka) and h k = 2iα sin(ka) + 4βu cos(ka), we find As in the linear SSH model case, we only need a single entry from these matrices at half-filling. Moreover, defining δK (2) nm as the second-order corrections to the n, m-binding atomic spring, by analogy we also have δK (2) nm = It is here that we have to work our way through several cases depending on the even/odd character of n, m. First, catalog the halves of the perturbative expressions for the even-odd cases. Even m: Odd m: Now for corresponding left-hand sides. The main difference is that now we keep the ν † χ combination. We immediately have Even m: When connected with the appropriate right-hand side, the equality q ′ = k, k ′ = q is ensured. Hence The only thing that changes for odd m is the prefactor: Odd m: Now assemble the full expressions: Case 1: n even, m even.
Case 4: n odd, m even.
The subsequent calculation of the phonon spectrum is carried out exactly as for the linear case.