Downfolding the Su-Schrieffer-Heeger model

Charge-density waves are responsible for symmetry-breaking displacements of atoms and concomitant changes in the electronic structure. Linear response theories, in particular density-functional perturbation theory, provide a way to study the effect of displacements on both the total energy and the electronic structure based on a single ab initio calculation. In downfolding approaches, the electronic system is reduced to a smaller number of bands, allowing for the incorporation of additional correlation and environmental effects on these bands. However, the physical contents of this downfolded model and its potential limitations are not always obvious. Here, we study the potential-energy landscape and electronic structure of the Su-Schrieffer-Heeger (SSH) model, where all relevant quantities can be evaluated analytically. We compare the exact results at arbitrary displacement with diagrammatic perturbation theory both in the full model and in a downfolded effective single-band model, which gives an instructive insight into the properties of downfolding. An exact reconstruction of the potential-energy landscape is possible in a downfolded model, which requires a dynamical electron-biphonon interaction. The dispersion of the bands upon atomic displacement is also found correctly, where the downfolded model by construction only captures spectral weight in the target space. In the SSH model, the electron-phonon coupling mechanism involves exclusively hybridization between the low- and high-energy bands and this limits the computational efficiency gain of downfolded models.

The standard ab initio method for calculating the EPI is the density-functional perturbation theory (DFPT) [44]. The most important ingredients of this theory are the adiabatic Born-Oppenheimer approximation [45], density-functional theory (DFT) [46], and linearresponse theory. Briefly put, these state that it is possible to separate the dynamics of the electrons and ions, treat the electron in an effective one-body Schrödinger equation, and calculate the response of the electrons upon displacement of the nuclei within linear order, based only on the electronic density [47][48][49]. The resulting EPI simultaneously describes two sides of the same coin, namely how the electrons screen and renormalize the phonons and how the electronic structure will adjust to atomic displacements. For an overview of the historical development and the recent accomplishments of calculating the EPI from first principles, see Ref. [50].
Despite the unquestionable success of the current ab initio computational methods, another trend in the literature is to treat the important physical phenomena in correlated ma-terials with downfolding approaches. The central idea is to reduce the number of degrees of freedom compared to the full system by keeping only the relevant states in a low-energy theory. The other states are integrated out and determine the parameters of the downfolded system. The overarching purpose of this procedure is the application of more advanced and expensive computational techniques only to the low-energy space where correlations take place.
For phonon-related properties, the constrained density-functional perturbation theory (cDFPT) was introduced [51] and successfully applied to superconducting materials such as alkali-doped fullerides [52] and light elements [53]. Additionally it was applied to monolayer 1H-TaS 2 [29], where it was shown that the CDW in this material is induced by coupling between the longitudinal-acoustic phonons and the electrons from an isolated low-energy metallic band. With the help of cDFPT it is possible to extract unscreened or partially screened parameters such as the phonon frequency and the electron-phonon vertex from an ab initio calculation and use these as the basis for an effective low-energy model Hamiltonian. The usefulness of partially screened parameters lies in the fact that they get rid of the coupling between phonons and the high-energy electrons.
As discussed, the description of real physical phenomena that are tightly linked to the EPIs is frequently based on ab initio theories (DFPT, cDFPT) that involve substantial numerical and computational effort. The structure of the theory is not always transparent, and also obscured by details of the numerical implementation. To avoid these complications, a second branch in the literature is focused on model Hamiltonians. The most popular models of the EPI are the Fröhlich model [54] for polaron formation, the Holstein model [55] for optical phonons, and the Su-Schrieffer-Heeger (SSH) model [56] for CDWs.
For understanding the interplay of electronic structure and atomic displacements, the SSH model is the most instructive since it explicitly describes how the electronic band structure is renormalized by the displacements of the atoms. Previous investigations using this model have studied properties such as the effective mass [57,58] and the band structure [59,60], but also phonon-related properties [61]. In the model, a periodic displacement of the atoms can open a band gap and thereby lower the total energy of the system [56], leading to a CDW transition. In other words, electronic screening makes the CDW phonon go soft. This textbook example of a CDW transition [62] is appealing for the investigation of downfolding since it is possible to perform all calculations exactly once the Born-Oppenheimer approximation has been applied.
The origin of this extraordinary simplicity lies in the observation that the Born-Oppenheimer approximation makes the phonons classical and the remaining electronic degrees of freedom in the SSH model are noninteracting. Thus, given any fixed displacement, the resulting electronic Hamiltonian is easily diagonalized. In some sense, this is similar to the method employed in Hirsch-Fye Quantum Monte Carlo [63], where a Hubbard-Stratonovich transformation is used to generate a system of noninteracting electrons coupled to classical fields and the subsequent analysis only involves varying the classical field and evaluating the noninteracting electron system. Unlike in Hirsch-Fye Quantum Monte Carlo, here the classical field is directly observable and has a clear physical meaning.
We choose to study the SSH model here for its simplicity, acting as a minimal model for electron-phonon coupling. At the same time, this means that there are many relevant aspects of electron-phonon coupling and CDWs that are not captured by the SSH model. In particular, the SSH model neglects Coulomb interactions between the electrons, and these are responsible for important effects such as screening and entirely electronic CDWs without lattice displacement. Furthermore, in higher dimensions, the shape of the Fermi surface can play an important role, in the form of nesting and Van Hove singularities. Given the complexity of electron-phonon systems, studying simple models is a useful way to identify relevant effects and mechanisms.
In this work, we compare the direct calculation of properties of the SSH model in the Born-Oppenheimer approximation at finite displacement with a perturbative diagrammatic expansion around the undistorted state à la DFPT. In this model, the diagrammatic expansion can be evaluated analytically order by order and we show that it correctly captures how the electron-phonon coupling renormalizes the phonon frequency and the electronic structure. Then, in the spirit of downfolding, we move to an effective single-band model for the dimerization transition in the SSH model. The diagrammatic structure in this effective model differs substantially from the original model: an interaction between an electron and two phonons appears and this interaction turns out to be dynamical with a frequency set by the high-energy electrons that were integrated out. We show that this downfolded model faithfully reproduces the energy landscape and the CDW. Furthermore, we discuss a cDFPT-like approach to downfolding, which correctly describes the screening of the phonon frequencies. In the SSH model, the displacement-induced orbital reconstruction between target and rest space is the central aspect of the downfolding and there is no remaining electron-phonon coupling in the cDFPT low-energy model.

Model
In this work, we consider the SSH model [56] in the classical Born-Oppenheimer limit [64], i.e., we ignore the kinetic energy of the atoms. We consider spinless fermions in a one-dimensional lattice with Hamiltonian Here, u i is a (classical) variable describing the atomic displacements, with 0 ≤ i < N . We use the periodic boundary condition u N ≡ u 0 . The hopping t > 0 sets the electronic energy scale and the force constant k s > 0 that of the phonons. We consider dimerization, i.e., displacements of the form u i = (−1) i α/2, and double the unit cell to include entire dimers. This is illustrated in Fig. 1a. Using the notation a i = c 2i and b i = c 2i+1 , we obtain Performing a Fourier transform to momentum space, the Hamiltonian in matrix form reads with eigenvalues ε ± (k) = ±2t 1 + (α 2 − 1) sin 2 (k) = ±2t cos 2 (k) + α 2 sin 2 (k).
These give the dispersion shown in Fig. 1b. Note that the Brillouin zone is −π/2 ≤ k ≤ π/2, where k is made dimensionless by setting the atomic distance to unity. In the following, we assume that the electronic density n is smaller than 1 electron/dimer. Since the model is particle-hole symmetric, the case n > 1 follows by symmetry. The situation n = 1 (half-filling) is special and will be discussed in more detail below, see Sec. 7. At zero temperature, the electron density is proportional to the Fermi wave vector k f and independent of α: n = 2k f /π. The total electronic energy per dimer, in the thermodynamic limit N → ∞, is and the total energy per dimer is Note that in this model, displacements do not change the Fermi surface and the electronic energy E el depends on α only via Eq. (5), which will allow us to pull derivatives through the integral in Eq. (6). In Fig. 2a, we show how the total energy depends on α for fixed k s and n . The total energy is obviously symmetric in α, and the undistorted lattice at α = 0 is an extremum of the total energy. Without electrons, E bare = k s α 2 is a convex parabola with a minimum at α = 0. However, the coupling to the electrons can lead to a Peierls CDW phase transition where α = 0 turns into a local maximum and two global minima occur at α = ±α * . The finite α lowers the energy of the occupied states and thus the total electronic energy and this compensates for the gain in potential energy due to α.

Harmonic and anharmonic lattice potential
To analyze the phase transition, it is useful to perform a Taylor expansion of the lattice potential E(α) around α = 0.
Here, we have introduced the bare phonon frequency ω bare and the dressed phonon frequency ω. The difference ∆ω 2 , the electronic screening of the phonon, originates in the change in electronic structure in response to the lattice distortion. Screening lowers the phonon frequency, and the Peierls transition occurs when the dressed phonon frequency is equal to zero, i.e., ω = 0. In Fig. 2b, the Peierls transition is represented as the black line that separates the phases ω 2 < 0 (Peierls instability) and ω 2 > 0 (no instability). As we can see, a weak force constant k s and a density n close to half-filling is preferred for a Peierls instability. Beyond the Peierls transition, α = 0 is a local maximum of the potential and the higher-order terms, such as h (4) , are responsible for ensuring that E(α) has a minimum at some finite α. In Appendix A, we show that there can be at most two minima, symmetrically located around α = 0. Only even orders of α appear due to the symmetry of the system.

Electron-phonon coupling: Two-band model
In the previous section, we used our knowledge of the exact dependence of the electronic structureε on α to determine the potential-energy landscape. In ab initio calculations (e.g., DFPT), one will usually not have access to this. Instead, the only known quantities are the electronic structure of the undistorted structureε 0 and the electron-phonon coupling, the first derivative of the electronic structure with respect to the displacement. Access to the latter quantity is guaranteed by the 2n + 1 theorem [47][48][49]. Because of this, it is instructive to calculate the (approximate) potential-energy landscape of the SSH model-and in particular the screening of the phonon frequency-based just on these quantities in a perturbative expansion around α = 0. The Feynman rules can be read off from the Hamiltonian, Eq. (3), by writing it aŝ Here, f † is shorthand for the vector (a † , b † ). There is a single q = 0 phonon mode corresponding to dimerization, with frequency ω 2 bare = 2k s . This mode is entirely classical, since we are interested only in a Born-Oppenheimer potential-energy landscape. The electron-phonon coupling is a matrix in electronic space and is obtained asĝ = dε/dα evaluated at α = 0. In other words, it consists of the parts ofε that are proportional to α. Explicitly, Note that we are considering a single phonon mode at q = 0, so we do not need a q label on g. The lack of higher-order electron-phonon-coupling terms in Eq. (14) is a special property of the SSH model. To evaluate the Feynman diagrams, it is most convenient to express the electronic part of the Hamiltonian in the eigenbasis of the unperturbed electronic system. This basis transformation can be seen in Appendix B. The transformed electron-phonon coupling iŝ We observe that g couples the two bands and has no intraband component. In other words, to linear order in α around α = 0, distortions only change the orbital composition of the bands but not the dispersion of the bands. The vanishing diagonal elements of g can also be understood as a symmetry selection rule. The inversion symmetry of the system implies that ε(α) and ε(−α) have the same eigenvalues and this implies both Tr g = Tr dε dα = d dα Trε = 0, which holds in any basis, and n|ĝ |n = 0 for any α = 0 eigenvector |n , since these |n are eigenvectors of the inversion operator with eigenvalue ±1.

Leading diagram
We are interested in establishing the effective potential felt by the atoms, including electronic screening. Diagrammatically, this means that the phonon mode only appears as external lines, whereas internally the diagram consists of electronic propagators and electron-phonon vertices. All diagrams with n external lines need to be summed to obtain the α n coefficient in the potential E(α). 1 For all upcoming diagrams, we will use the electronic Green's function where1 is the identity matrix, the division denotes matrix inversion, andη k denotes the usual small imaginary constant that is positive (negative) for empty (occupied) states, respectively. For the phonon self-energy, i.e., with two external lines, there is only a single diagram, shown in Fig. 3a for one possible choice of the band indices, which corresponds to This allows us to simplify the result to This is consistent with Eq. (12). This shows that the harmonic energy landscape can be calculated entirely from the undistorted structure at α = 0, based on the electronic dispersion ε 0 and the electron-phonon couplingĝ.

Higher-order diagrams
It is also possible to calculate the energy landscape beyond the quadratic term. A special property of the SSH model is that the there are no higher-order electron-phonon vertices nor anharmonic bare phonon terms. Because of this, the entire perturbation theory is expressed in ε ± and g. For example, the diagram for the fourth order contribution α 4 is shown in Fig. 3b. This is the only connected diagram at this order. 2 Note that all external phonons have q = 0, so all electronic lines have the same momentum k and energy E. The band index of the electronic lines is alternating, since the electron-phonon coupling is entirely off-diagonal. The expression corresponding to this diagram is of the form which already includes a factor 2 accounting for the fact that there is a second way to assign the band indices. 3 The product of Green's functions can be reduced by repeated application of the relation AB = (B−A)/(A −1 −B −1 ) for A = B, which is helpful because G −1 ± (k, E) = E∓|ε 0 (k)|+iη k is very simple. Below, all G's have the same arguments k, E, which were dropped for notational convenience.
In the denominators we have already safely taken the limit η → 0. Now, the integral over E can be performed using dEG 2 ± (E) = 0 and dEG ± (E) = n(ε ± (k)). Here, n(ε ± (k)) is the occupation, which is unity for the − branch and |k| < k f and zero otherwise. This gives the same result as Eq. (13), Diagrams at higher order can be evaluated in the same way, by repeated simplification of products of Green's functions. An interesting aspect is that the entire potential-energy landscape E(α) can be calculated in this way (for n = 1, see Sec. 7) without ever determining how the band dispersion changes.

Change in electronic structure
The change in the electronic structure is given by the self-energy Σ(E, k) and can be obtained diagrammatically by considering the sum of all one-electron irreducible diagrams. Now, the electronic lines are amputated and the phononic ends of the vertices are connected to crosses representing α. This is similar to the way an external Zeeman magnetic field or scattering potential can be included in a diagrammatic theory. Note that due to the Born-Oppenheimer approximation, there is no true phonon propagator with two end points, which would represent the phonon dynamics. In the present model, it turns out that there is only a single, trivial diagram for the self-energy, with an equivalent diagram for Σ −+ . Together, they recover the exact electronic Green's functionĜ via the Dyson equation,

Single-band effective model
At n < 1, there is only one partially filled band and this motivates us to investigate the possibility of describing the CDW via a single-band effective model. Here, we construct a model consisting of the partially filled electronic band, the bare phonon, and the coupling between the two. Formally, such a model is obtained by integrating out the unoccupied band of the two-band model. The effective action of the single-band model contains (partially) renormalized, dynamically screened interactions between these electrons and the phonons. In fact, the interaction vertices in this effective model can and do have an entirely different structure compared to those of the original two-band model. Generally, the vertices in the effective theory are obtained by collecting all connected diagrams consisting of rest space (here: ε + ) internal lines with a particular number of external phonon and target space (here: ε − ) lines, and an infinite set of vertices can appear in this way. The only general constraints are the conservation of the fermion number and momentum conservation. Thus, the low-energy Hamiltonian can contain interactions of the form α m (c † c) n for arbitrary m and n. However, additional symmetries of the system can provide further constraints on the effective action.
Here, the single-band model is energetically completely symmetric in α ↔ −α and this implies that only even powers of α can appear in the effective action. In other words, only interaction vertices with an even number of phonon lines are allowed. 4 In fact, looking at the diagrammatic structure, it turns out that the single-band effective theory of the SSH model only contains one interaction vertex, shown in Fig 4a. This vertex has two phonon and two external electronic lines (one incoming, one outgoing) and takes the value Note that V depends explicitly on E; the screened interactions that enter the effective model are dynamical quantities. The effective model contains only a single fermion with dispersion ε − , so no further electronic band label is necessary. The downfolded SSH model has only a single effective interaction vertex. This happens because the electron-phonon coupling in the original SSH model only has a single external high-energy electron (blue line in Fig. 4a). On the other hand, if the original model had contained either electron-electron interactions in the high-energy band or electron-phonon coupling between different electronic states in the high-energy band, then the downfolding would be more involved, since more diagrammatic contributions would appear in the expression for the effective action.
For the energy E(α), the second-order contribution, shown in Fig. 4b, is which upon insertion of Eq. (27) is equal to the result we obtained in the two-band model. Similarly, the fourth-order contribution, also shown in Fig. 4b, is The denominator can be simplified using the same techniques as above and this gives the same final result as the earlier expression for h (4) .

Change in electronic structure in the single-band model
In the effective model, only the electronic target space is considered, corresponding to the lower band at zero displacement. The rest space has been integrated out. Σ is now a scalar quantity and it is once again given by a single diagram, In this case, Σ(k, E) is an explicit function of E and it is not possible to interpret it purely as a change in the dispersion. Since the true change in the electronic structure involves a change in the orbital composition of the bands and thus coupling between the bands and changes in the wave functions, it is not possible to capture this entirely in a single-band model. However, if we restrict ourselves to the vicinity of the lower band in terms of energy, we find which is equal to the exact second-order expansion of Eq. (5). At the same time, the self-energy of Eq. (30) has a pole at E = |ε 0 |, the energy of the upper band that has been integrated out. In the spectral function A(E, k) = − 1 π Im G(E, k), this shows up as interaction-induced spectral-weight transfer, as shown in Fig. 5. The original spectral weight of the noninteracting, i.e., undistorted downfolded model (grey peak) is distributed to the positions of the lower and upper band of the interacting, i.e., distorted model (orange peaks). Thus, even though it cannot represent the matrix structure of the electronic Green's function, the downfolded model has spectral weight at the right locations. Note that there is no imaginary part in the self-energy and thus no additional broadening of these peaks in the downfolded model; all broadening comes from the constant η = 0.05 used for plotting the spectrum.

Constrained density-functional perturbation theory
The downfolding procedure employed above is based on an explicit resummation of the diagrammatic series and is able to reproduce the screening from bare to dressed lattice potential exactly. This approach can be applied here, since we have full knowledge of the entire electronic structure and the electron-phonon coupling. In ab initio calculations, the downfolding is usually done somewhat differently. Indeed, cDFPT is a tool commonly used for downfolding electron-phonon systems onto an electronic target space and calculating corresponding partially screened phonon frequencies. In general, it evaluates a Feynman diagram similar to Fig. 3a, with the restriction that at least one of the two electronic propagators shall not be part of the target space.
In the SSH model, if the lower band is chosen as the target space, cDFPT includes the only relevant screening process, with one + and one − electron, in its calculation of the partially screened phonon frequency. In other words, Here Π m,n (k) is defined and used as in Eq. (19). In the SSH model, Π −− anyway does not contribute to the phonon renormalization, and as a result the cDFPT phonon frequency is identical to the fully screened phonon frequency. The cDFPT low-energy model then basically consists of the fully screened phonon, the lower electronic band, and no electron-phonon coupling, since g −− = 0. Because of this special property of the SSH model, there is no real distinction between the partially and fully screened phonon.

Breakdown of perturbation theory at half-filling
The series expansion of the potential E(α) around α = 0, performed either diagrammatically or by directly taking derivatives of ε(k, α), shows a regular pattern. Only even powers of α are allowed. For a given power α 2n , the diagrammatic contribution will be of the form (modulo prefactor) g 2n G n − G n + . The 2n electron-phonon vertices g contribute (2t) 2n sin 2n (k), whereas the Green's functions can be reduced to n(ε − (k))/(2ε 0 ) 2n−1 ∝ n(ε − (k))/ cos 2n−1 (k). The only role of the density is to determine the integration range, via k f . This becomes qualitatively important for n → 1, k f → π/2, since ε 0 (k f ) → 0. The denominator in the integral diverges and as a result the entire integral is no longer convergent. In other words, perturbation theory around α = 0 is not possible since E(α) is not an analytical function anymore. Physically, the dimerization at half-filling is a Peierls transition caused by the perfect nesting of the Fermi surface points ±π/2 with respect to the dimerization wave vector π (in the original Brillouin zone). Thus, at half-filling, dimerization will occur even at arbitrarily large force constant k s .

Conclusion and discussion
A key question in the investigation of coupled electron-phonon systems is the evolution of the total energy and electronic structure as a function of atomic displacement. In ab initio studies, it is desirable to gain (perturbative) access to this energy landscape starting from the undistorted structure and a small set of relevant electronic bands. In the SSH model, it is actually possible to perform this perturbative, diagrammatic expansion analytically and to trace the performance of effective models. This both provides a unique insight into "exact downfolding" and highlights the successes and possible failures of effective models.
The bare phonons in the SSH model are entirely harmonic by definition. Thus, all anharmonic effects in the potential energy have to be created by the (linear) coupling to the electrons and the resulting electronic screening. Due to the simple structure of the model, the screening can be calculated to arbitrary order in the displacement. It reduces the energetic cost of displacements and eventually leads to a CDW transition, i.e., the appearance of a new global minimum in the energy landscape at a finite displacement. In this model, all relevant quantities can be reduced to integrals over the occupied part of the Brillouin zone.
It is also possible to downfold onto a single-band model with only half the electronic degrees of freedom of the original system. The diagrammatic structure changes due to the downfolding; the electron-phonon coupling is now dynamical and quadratic in the displacement field. Still, the analytical evaluation of the diagrams determining the energy landscape is possible and agrees with the exact result. Regarding the electronic structure, the effective single-band model only has the ability to describe spectral-weight transfer and by construction does not have the ability to describe the changes in the orbital composition of the bands as the atoms move. In the cDFPT approach, as well as in the cRPA approach to Coulomb interactions, these changes in the electronic structure are usually not considered at all. This observation is potentially relevant for several two-dimensional transition-metal dichalcogenides. For example, monolayer 1H-TaS 2 has a single band crossing the Fermi level and this band consists of a combination of d 0,+2,−2 orbitals. It was already known that the electronic matrix structure is imprinted on the momentum structure of the electron-phonon coupling in ab initio downfolding [29] and that the resulting single-band electron-phonon model accurately describes the phonon frequencies (i.e., the energy landscape close to the undistorted structure). A similar situation, with a single composite band crossing the Fermi level, occurs in 1H-NbS 2 [66]. An open question is how these single-band effective models perform in the description of the true electronic structure of the distorted phase. If the distortions lead to hybridization between target and rest space, downfolded models can only capture the spectralweight transfer. On the other hand, downfolded approaches can fully describe processes that occur entirely in the target space. Thus, fluctuation diagnostics of the electron-phonon coupling [29] can provide an answer to this question.
The SSH model in the Born-Oppenheimer approximation-as studied here-is very much a simplification of the complex reality of electron-phonon-coupling and charge-density-wave physics. We assume that the lattice is one-dimensional, that the electronic hopping amplitudes and the bare restoring forces are linear in the displacement, that there is no electron-electron interaction, that there is a single relevant phonon mode (dimerization), and that the system is in the T = 0 ground state. Still, some general conclusions are possible from our work. It is possible to generate anharmonic phonon terms entirely electronically, from an initial Hamiltonian that has purely harmonic phonons. Diagrammatic expressions can be constructed for the electronic screening at and beyond the harmonic level; in the general case these will be infinite series of diagrams, but here there is only a single diagram at any order in the displacement. In the presence of multiple relevant phonons, see Appendix C, the Born-Oppenheimer energy landscape will include mode-mode coupling as well. Downfolding of the electronic space generates a new perturbation series, in which effective higher-order vertices appear naturally. Unlike in the original Hamiltonian, the vertices of the downfolded system are also dynamical (frequency-dependent). As a result, the self-energy is dynamical as well, leading to spectral-weight transfer in the downfolded model. We note that this happens even though the electrons are noninteracting. The magnitude of the self-energy in the low-energy band is approximately given by the electron-phonon coupling (between the target and the rest space) squared times the displacement squared divided by the energy separation between the low-energy and the high-energy band. This supports the natural strategy of including bands in the low-energy model that are close in energy and those that are strongly coupled to the target space via the relevant phonon modes.

A Number of minima of E(α)
The SSH model in the limit of large α is unlikely to be an accurate description of any real physics, but it is useful to establish some formal results. First of all, the triangle inequality provides us with bounds on the dispersion, Thus, in the limit of large α, ε(k) is roughly proportional to α sin(k). The total energy is then dominated by the purely lattice term proportional to k s α 2 . We conclude that the energy landscape E(α) is bounded from below, as it should be. Two types of energy landscape E(α) are discussed in the text, one with a single minimum at α = 0 and one with two minima at α = ±α * . In fact, we can proof that these are the only two possibilities, no further local minima are allowed.
First, we define the auxiliary function f (x) = − √ 1 + x 2 , so that We observe that the second derivative of f , f = −(1 + x 2 ) −3/2 , is monotonously increasing for x ≥ 0. This implies that d 2 ε − (k, α)/dα 2 is also monotonously increasing as a function of α for α ≥ 0 and the same holds for E(α), which is just a k-integral over ε − . Thus, there can be at most one α ≥ 0 where d 2 E(α)/dα 2 = 0. In E(α), local minima (d 2 E/dα 2 > 0) and local maxima (d 2 E/dα 2 < 0) alternate, so by the intermediate value theorem d 2 E/dα 2 must cross zero between every local optimum of E(α). This can happen only once for α ≥ 0, so there are at most two optima at α ≥ 0 and one of them is at α = 0 by symmetry. Since E(α) → +∞ for α → +∞, there is either a single global minimum at α = 0 or a local maximum at α = 0 and two global minima at ±α * .

B Basis transformation of the electron-phonon coupling
The electronic part is most conveniently expressed in the band basis ofε 0 , which isε evaluated at α = 0. The two eigenvalues ofε 0 are ε ±,0 = ±2t cos(k) with corresponding eigenvectors With the eigenvectors, we can form the transformation matrix which diagonalizesε 0 . This yields the electron-phonon coupling in the band basis,

C Beyond dimerization: 4-site unit cell
At half-filling, the dimerization is commensurate in the sense that 2k f = q dimerization . We have already shown that dimerization can also be energetically favorable away from half-filling, but so far we have not considered CDWs with other periodicities. In this appendix, we consider periodicity 4, which allows for the study of additional phonon modes. Because this doubling of the unit cell increases both the number of phonons and the number of electronic bands, it is more difficult to derive compact formulas and our treatment remains relatively brief, highlighting some similarities and differences to the 2-site unit cell. In this case, it is convenient to first consider the electronic dispersion as a function of the four hopping amplitudes t ij , as shown in Fig. 6. In the SSH model, these hopping parameters will be linear functions of the atomic displacements.
Here, −π/4 < k ≤ π/4 is the Brillouin zone corresponding to this unit cell. As for the dimerization transition, the total electronic energy is given by m dk ε m (k)n(ε m (k)). Now, in the SSH model, the hopping parameters depend linearly on the atomic displacements. We consider three phonon modes α 1 , α 2 , and α 3 defined by t 12 = t(1 + α 1 + α 2 ), α 1 is the dimerization mode studied in the main text, α 2 is sketched in Fig. 6, and α 3 is obtained from α 2 by translating the unit cell by one atom. They are eigenmodes at q = 0. Combining Eqs. (39) and (40), it is possible to calculate ε(k; α 1 , α 2 , α 3 ) and its derivatives with respect to α i . Using computer algebra, it is possible to evaluate these derivatives straightforwardly, although the expressions quickly become unwieldy. Below, we will briefly discuss the nonzero terms at the lowest orders. Finally, integrating these derivatives of the dispersion over the filled part of the Brillouin zone (for each band) then gives the terms in the Taylor expansion of E(α 1 , α 2 , α 3 ), as in Sec. 3 of the main text. The first derivative vanishes as expected, ∂ α 1 ε = ∂ α 2 ε = ∂ α 3 ε = 0. The second derivative is diagonal in the phonon index, ∂ α 1 ,α 2 ε = ∂ α 1 ,α 3 ε = ∂ α 2 ,α 3 ε = 0, so the only nonzero elements are ∂ α 1 ,α 1 ε and ∂ α 2 ,α 2 ε = ∂ α 3 ,α 3 ε. At the level of the third derivative, we find a finite term with mixed phonon labels, to be explicit: Here, the four components in the vector correspond to the bands from lowest to highest energy, and we have assumed k > 0. At fourth order, we find nonzero expressions only for the terms where the derivatives appear in pairs, e.g., ∂ α 1 ,α 1 ,α 2 ,α 2 ε. Symmetries and momentum conservation still ensure that many terms in the expansion vanish, but already at the third order we see that qualitatively new terms appear compared to energy landscape for the 2-site unit cell. In other words, the Feynman diagrams studied in the main text are all relevant in general, but diagrams that were "forbidden" in that simple system can play a role. It is difficult to make any statements about the sign and relative magnitude a priori; for a computational case study of nonlinear mode-mode coupling, see Ref. [67]. Similarly, the Hamiltonian can be written in terms of the bare dispersion and the electronphonon couplings, now as 4 × 4 matrices. In analogy to the main text, the terms in the expansion of E(α 1 , α 2 , α 3 ) can then be obtained diagrammatically.