Competing Orders and Unconventional Criticality in the Su-Schrieffer-Heeger Model

The phase diagram of the one-dimensional Su-Schrieffer-Heeger model of spinless fermions coupled to quantum phonons is determined by quantum Monte Carlo simulations. It differs significantly from previous work and features Luttinger liquid, bond-order-wave (BOW), and charge-density-wave (CDW) phases. The retarded phonon-mediated interaction gives rise to physics known from the frustrated XXZ chain. Most notably, a continuous BOW-CDW quantum phase transition characterized by unconventional power-law exponents can be interpreted in terms of deconfined quantum criticality involving proliferation of solitons (Mudry et al. [1]).

acting on bond b between sites i and j = i + 1 andĤ ph = b Equation (1) describes itinerant spinless fermions coupled to optical bond phonons with momentumP b , dis-placementQ b , and frequency ω 0 = K/M . For the present case of a half-filled band ( n i = c † i c i = 0.5), numerics [32,33] and field theory [34] suggest the same physics for optical phonons-which avoid a QMC sign problem-and the original acoustic phonons [18]. We use the dimensionless coupling λ = g 2 /Kt and set , k B = 1.
After integrating out the phonons, the partition function for Eq. (1) contains a retarded interaction The free phonon propagator P (τ ) is local in space but its decay in imaginary time τ (here, β = 1/T ) is determined by ω 0 , P (τ ) ∼ e −ω0τ . The associated retardation effects are crucial for the phase diagram in Fig. 1.
(1) yields the fermionic hopping termĤ 0 = − b [t + (−1) b ∆]B b studied in the context of topological insulators. The Peierls argument  [35] implies that the dimerization ∆ is nonzero for any λ > 0 and opens a gap at the Fermi level. Quantum lattice fluctuations can destroy long-range order at sufficiently weak coupling and thereby allow for a Luttinger liquid (LL) to BOW quantum phase transition (QPT) at a finite λ c (ω 0 ) [32,34,36]. An exact solution (by the Bethe ansatz) is also possible in the opposite, antiadiabatic limit. For ω 0 → ∞, the interaction (2) becomes instantaneous and Eq. (1) maps to the t-V model H ∞ = −t bB b + V in ini+1 with V = λ [34] and an LL-CDW QPT at V c /t = λ c = 2 [37]. The BOW and CDW states spontaneously break translation symmetry. Long-range order at T = 0 is described by Ising order parameters that reflect the two possible BOW (CDW) dimerization patterns related by a shift by one lattice site. CDW and BOW states can be distinguished by point group symmetries. As illustrated in Fig. 1, CDW order breaks bond inversion symmetry but preserves invariance under site inversion. The opposite is true for BOW order. Because different symmetries of Hamiltonian (1) are spontaneously broken, a Ginzburg-Landau theory would correctly suggest the absence of an adiabatic connection between the ground states at ω 0 = 0 and ω 0 = ∞, an aspect ignored in previous work. Apart from its origin in retardation effects, a generically continuous BOW-CDW QPT is captured by the theory of Refs. [38,39] for a frustrated XXZ chain, see below.
Method.-To connect the two exact limits, we used a state-of-the-art QMC method based on the stochastic series expansion [40] and directed-loop updates [41]. The performance gain from extending the latter to retarded interactions [42] is essential to explore the phase diagram of the SSH model. The method has only statistical errors and relevant technical details are summarized in the Supplementary Material (SM) [43]. All results were obtained for periodic chains of L sites at inverse temperatures βt = 2L.
In the nonadiabatic regime [t/ω 0 = 1/10, Fig. 2(b)], we find behavior consistent with an LL at λ = 2 and longrange CDW order at λ = 5. Results for the BOW-CDW transition at λ = 6 are shown in Fig. 5(c). The charge stiffness D ρ converges to a nonzero value (vanishes) for L → ∞ in a metallic (insulating) phase [45,46]. In 1D and at sufficiently low temperatures [47], it can be measured using the winding number estimator for the superfluid stiffness [43,48]. In conjunction with Fig. 2 Fig. 3(a) hence reveals a metal-insulator LL-CDW transition with increasing λ at t/ω 0 = 1/10. If we instead start in the CDW phase at λ = 4, an increase of t/ω 0 appears to drive two consecutive QPTs separated by a metallic region [ Fig. 3(b)]. The correlators in Fig. 2 identify the critical points as CDW-LL and LL-BOW QPTs, respectively. Comparing Fig. 3 Fig. 5(a) reveals that a potential intermediate phase narrows with increasing coupling, but a peak in D ρ can be observed even at λ = 12 [43].
A renormalization-group (RG) analysis of umklapp interactions would suggest Berezinskii-Kosterlitz-Thouless (BKT) LL-BOW and LL-CDW transitions with a critical value K ρ = 1/2 and numerically challenging logarithmic scaling due to a marginally relevant operator [44]. This has been explicitly confirmed for t/ω 0 = 0 [37]. A functional RG study of the SSH model [36] reported unconventional BKT critical behavior with a critical K ρ < 1/2. We will argue below that key aspects of the SSH model are captured by a theory for the frustrated XXZ chain, which predicts conventional BKT LL-BOW and LL-CDW transitions but K ρ < 1/2 along the BOW- CDW critical line [38]. Large-scale simulations of classical frustrated 2D XY models [49,50] indicate a standard BKT transition [49,51] (see, however, Ref. [52]), albeit with challenging crossover phenomena.
To obtain the phase boundaries in Fig. 1, we analyzed the finite-size scaling of D ρ (L). The universal stiffness jump at the critical point of the 2D XY model [53] translates to D ρ (∞) = t/2 for the t-V model (the SSH model with t/ω 0 = 0) [37]. For t/ω 0 > 0, we instead find nonuniversal stiffness jumps, with D ρ (L) < t/2 even for small L in, e.g., Fig. 3(b). Possible origins are discussed in the SM [43] and include velocity renormalization, geometry effects, and frustration. We used the first-order scaling ansatz [54] with fit parameters D ρ (∞), g, and C [43]. As for the 2D XY model, fits to Eq. (4) reveal the critical point as a minimum in goodness-of-fit measures (here: the reduced chi-squared χ 2 ν = χ 2 /ν for ν degrees of freedom) [54,55]; for applications to quantum models see Refs. [56,57].
Discussion.-BOW order at t/ω 0 = ∞ and CDW order at t/ω 0 = 0 were established exactly in Ref. [34] by the arguments outlined after Eq. (2). However, the phase diagram was predicted to have a single, monotonic BKT phase boundary separating a LL from a dimerized phase. Similar conclusions were reached by (functional) RG calculations [36,59] and for related spin-phonon models [36,[60][61][62]. A nonadiabatic mean-field approach yields BOW and CDW phases even for large ω 0 [63]. These findings differ significantly from ours.
A phenomenological description of the different phases in Fig. 1 is provided by the Goldstone-Wilczek theory of Dirac fermions ψ = (ψ A , ψ B ) (A/B: sublattice components, see Fig. 1) coupled to classical phonon fields (we drop a phonon term L ϕ ) [64], The masses can be combined into m = (m bow , m cdw ) so that the spectrum E(p) = ± p 2 + |m| 2 [65]. For |m| = 0, L has a U(1) chiral symmetry generated by γ 5 [19]. The soliton charge is e/2 for m cdw = 0 [66] but can be arbitrary for m cdw > 0 [64,67]. Spontaneous mass generation is described by the Gross-Neveu model [68] The umklapp interactions reduce the chiral symmetry to a discrete Ising symmetry and are directly linked to the commensurately filled lattice [39]. In the context of Gross-Neveu theories, the CDW phase observed here for the SSH model is known as an Aoki phase [69]. The BOW phase constitutes an interaction-generated topological insulator [13], adiabatically connected to topological band insulators in the so-called BDI class [70].
A bosonized theory that captures spontaneous BOW and CDW order from a single interaction-as appropriate for Eq. (1)-as well as other key aspects of our findings is that of the frustrated J 1 -J 2 XXZ chain [38,39], The cosine term is irrelevant for K ρ > 1/2 (LL phase) and on the BOW-CDW transition line (where λ φ = 0 [1]). It is a relevant perturbation in the BOW and CDW phases, which are associated with opposite signs of λ φ and a pinning of the charge mode φ at different minima [1,44]. A relation between the SSH electron-phonon model and the frustrated J 1 -J 2 XXZ model emerges from the mapping of phonon-mediated retarded interactions to frustrated spin interactions [62,[71][72][73].
Within the theory (7), QPTs between the LL and symmetry-broken phases are conventional BKT transitions [38]. The unusual critical behavior at the BOW-CDW transition in Fig. 5 mirrors that along the line of continuous dimer-Néel transitions of the frustrated XXZ chain, with a continuously varying exponent η = 2K ρ < 1 [38]. Apart from λ = 6, see Fig. 5, we also find evidence for K ρ < 1/2 at criticality for λ = 4 [43], consistent with a location on the BOW-CDW transition line. Therefore, the two separate critical points (with significant uncertainty) in Fig. 1, inferred from Fig. 4(e), may be an artifact of challenging finite-size scaling in the tricritical region of the phase diagram.
A physical picture of how BOW and CDW phasescharacterized by different broken symmetries-can be connected via a generically continuous phase transition is provided by the scenario of 1D deconfined quantum criticality of Ref. [1]. It is based on solitons in the CDW (BOW) order parameter that can be added in pairs and interpolate between the two degenerate CDW (BOW) configurations. Parameterizing the phase of the order parameter by ϕ = (cos ϕ, sin ϕ) [1], see inset of Fig. 1, BOW (CDW) patterns correspond to ϕ = 0, π (ϕ = ±π/2). For example, a defect in the BOW order connecting ϕ = 0, π contains a region with CDW order or ϕ = π/2. Simultaneous proliferation of BOW/CDW defects at ω 0,c provides a mechanism for a continuous transition without fine-tuning. Whereas the mean-field theory (5) suggests a BOW-CDW 'transition' via rotation of m while keeping the gap |m| open, the theory of Refs. [1,38] yields |m| → 0 at the critical point together with an emergent U(1) symmetry. This agrees with our numerical results for the spinless SSH model (1), for which metallic behavior entails a vanishing single-particle gap. Interactiondriven QPTs out of a topological band insulator were recently addressed in Refs. [20,26,74] (see also Ref. [14]), where BOW solitons are excluded and the BOW-CDW transition exhibits Ising criticality [26].
The BOW-CDW transition is driven by retardation, which is difficult to capture analytically [33,36,75]. An effective Hamiltonian of the form (7) has to be obtained from an RG treatment of the phonon-mediated interaction, accounting for the different signs of λ φ in the BOW and CDW phases and λ φ = 0 at the transition. The onset of CDW order at ω 0 W suggests that the free bandwidth W = 4t is an important energy scale. Hence, the standard bosonization approach based on a linearization of the electronic dispersion around the Fermi points (equivalent to W → ∞) may be insufficient. A CDW phase is absent in functional RG results [36].
Outlook.-Generalizations of our work include spinful fermions and higher dimensions. For the 1D spinful SSH model, a transition from BOW to a critical spin-densitywave state is expected, whereas a 2D SSH model supports valence-bond and antiferromagnetic phases [76].
We acknowledge helpful discussions with F. Assaad, A.

Quantum Monte Carlo Method
We used the directed-loop QMC method for retarded interactions in the path-integral representation [42]. It is based on an interaction expansion of the partition function Z = D(c, c) e −S0−S1 around S 0 = dτ ic i (τ ) ∂ τ c i (τ ). A general interaction vertex S 1 = − ν w ν h ν can be written as a sum over vertex variables ν, a weight w ν , and the Grassmann fields contained in h ν . The perturbation expansion becomes with sums over the expansion order n and the ordered vertex list C n = {ν 1 , . . . , ν n }. For each time-ordered configuration of vertices, the expectation value over Grassmann fields can be represented by world lines. The trivial choice of S 0 ensures that the imaginary-time evolution is entirely determined by the interaction vertices. Therefore, Eq. (S1) is the path-integral equivalent of the stochastic series expansion (SSE) representation where Z = Tr e −βH is expanded in the total Hamiltonian [40,77]. Accordingly, many algorithmic features, including the global directed-loop updates [41], directly transfer to the path-integral representation [77]. The retarded interaction of the SSH model includes two bond operators acting at different imaginary times. Therefore, a compatible interaction vertex must contain two subvertices j ∈ {1, 2} with local variables {a j , b j , τ j } labeling the operator type, bond, and time of each operator. For the SSH model with a coupling to optical bond phonons we have b 1 = b 2 = b. The interaction vertex of the SSH model becomes It is important to note that the symmetrized phonon propagator P + (τ ) = ω 0 cosh[ω 0 (β/2 − τ )]/[2 sinh(ω 0 β/2)] is included in the global weight w ν of the vertex. Whereas the bond-bond interaction is already nonlocal in time, the single hopping terms of the kinetic energy are promoted to retarded interactions by including unit operators with a second time variable, i.e., This is possible because β 0 dτ 2 P + (τ 1 − τ 2 ) = 1. As the vertices (S3) and (S4) both contain off-diagonal hopping operators, we have to include a purely diagonal term in the interaction vertex. The simplest choice is a constant shift of the action, With our choice of interaction vertices, we can formulate the diagonal and directed-loop updates similar to the SSE representation [41]. For the diagonal updates, we use the Metropolis algorithm to add and remove vertices h 00,b (τ 1 , τ 2 ) that do not change the world-line configurations but change the expansion order n. We propose time differences τ 1 −τ 2 according to the phonon propagator using inverse-transform sampling. Because P + (τ 1 − τ 2 ) appears as a global weight in front of each vertex, it drops out of the directed-loop equations. The latter can be solved for each vertex similarly to the original approach, see the Supplemental Material of Ref. [42]. The constant k in Eq. (S5) has to be chosen such that every weight in the loop assignments is positive. During the propagation of the directed loop, unit operators can be transformed into bond operators and vice versa, leading to local updates h 00,b ↔ h 10,b /h 01,b ↔ h 11,b . Note that the vertices are constructed in such a way that each subvertex can be changed individually while the other subvertex remains unchanged. For details on the updating schemes, we refer to Refs. [41,42].
The calculation of observables in the path-integral (interaction) representation is in many ways similar to the SSE representation. Sandvik et al. [77] systematically compared estimators for electronic correlation functions derived in the two representations. Estimators that only include diagonal operators, such as the charge structure factor C ρ (r) and the charge susceptibility χ ρ (r), are simple to derive and given in Ref. [77]. Estimators including off-diagonal operators can often be recovered from the vertex distribution if there is a vertex that only includes this operator. Measuring the static or dynamic correlations functions of two bond operators at arbitrary bonds b 1 and b 2 is only possible when considering the hopping vertices h 10,b /h 01,b . It turns out that the bond susceptibility χ b (r) has a very simple estimator where only the total number of hopping vertices at bonds b 1 or b 2 has to be computed, see Ref. [77] for the exact estimator. However, calculating the equal-time bond structure factor C b (r = b 1 − b 2 ) = B b1 B b2 in the interaction representation is more involved. While a general derivation is outlined in Ref. [77], we only state the final estimator for the SSH model. For a Monte Carlo configuration C n , the bond structure factor can be estimated from In principle, the sum over p runs over the time-ordered list of all subvertices contained in a world-line configuration. However, we can exclude the unit operators 1 b as they were only introduced to simplify the Monte Carlo sampling. I b1b2 (p − 1, p) is zero unless bond operators B b1 (τ p−1 ) and B b2 (τ p ) originating from the hopping terms h 10 /h 01 appear at adjacent times; then I b1b2 (p − 1, p) = 1. An integral expression for K(p − 1, p) was derived in Ref. [77] and gives K(p−1, p) = 2/(τ p+1 −τ p−2 ) when 4 or more subvertices are present in a world-line configuration. The time difference τ p+1 − τ p−2 ∈ [0, β] is defined by the two subvertices that surround the two bond operators under consideration. Note that K(p − 1, p) = 2/β for 3 subvertices, K(p − 1, p) = 1/β for 2 subvertices, and K(p − 1, p) = 0 for 0 or 1 subvertices. For further details, see Ref. [77]. The Monte Carlo configurations do not give direct access to observables containing phonon fields because the latter have been integrated out to obtain a retarded fermionic interaction. However, bosonic observables can be recovered from electronic correlation functions using generating functionals. In particular, we derived efficient estimators for the total energy, specific heat, fidelity susceptibility, and phonon propagator in Refs. [78,79] that make use of the vertex distribution. In the following section, we use the framework outlined in Ref. [78] to show that the superfluid stiffness of an electron-phonon model can still be calculated from the winding number.

QMC estimator for the superfluid stiffness
Consider a ring of length L threaded by a magnetic flux φ. At finite temperatures, the superfluid stiffness can be obtained from the free energy via [45] Because we consider a 1D system [47] and our simulations at β = 2L are essentially converged with respect to temperature, the measured values of ρ s are representative of the charge stiffness or Drude weight defined as [80] where E is the ground-state energy. Using F = − 1 β ln Z, the stiffness is directly related to the action of the SSH model. The magnetic flux can be incorporated by imposing twisted boundary conditionsĉ L+1 = e iφĉ 1 . The boundary term of the action reads Here, S L/R is the action of the hopping term (S4) crossing the boundary to the left/right, whereas S LL/RR corresponds to the bond-bond interaction (S3) with both hopping operators going to the left/right. The superfluid stiffness can then be calculated as The first expectation value is given by For each Monte Carlo configuration, expectation values of terms S a contained in the interaction vertex (S2) can be obtained by counting the number of vertices n a [78]. For the Monte Carlo average we then obtain S a = − n a . In the same way, the second term in Eq. (S10) becomes = − (S L + S R ) + 4 (S LL + S RR ) = (n L + n R ) + 4 (n LL + n RR ) (S12) and the third term is given by where we used S a S b = n a n b − δ ab n a . We get an additional shift for a = b that cancels the contribution of (S12). Our results are equivalent to calculating the winding number W = n B L − n B R where n B L/R counts the number of subvertices B b (τ ) crossing the boundary to the left/right. Here, n LL/RR contributes with a factor of 2 because each vertex contains two bond operators, whereas mixed contributions n LR drop out. Therefore, ρ s can be calculated in the same way for retarded interactions as for equal-time interactions, i.e., ρ s = L( W 2 − W 2 )/β [48]. Figure S1 shows D ρ (L) and K ρ (L) as a function of t/ω 0 for λ = 4, 6, 8, 12. For all couplings, the data are consistent with a metallic region at intermediate t/ω 0 . Whereas the apparent narrowing of this region between λ = 4 and λ = 6 matches the phase boundaries in Fig. 1, the theory discussed in the main text suggests that the BOW-CDW transition involves a gap closing and hence metallic behavior only at a single point. At this transition, the LL parameter K ρ < 1/2, in accordance with Fig. S1(e). Values K ρ < 1/2 can be reconciled with metallic behavior by assuming λ φ = 0 in Eq. (7) at the BOW-CDW critical point [38,39].

Additional data
In contrast to Ginzburg-Landau theory, the BOW-CDW transition does not require fine-tuning of both t/ω 0 and λ. For a fixed λ, λ φ can be tuned to zero for a suitable value of t/ω 0 , giving rise to a line of critical points. Since K ρ < 1/2 at criticality, any nonzero λ φ yields long-range BOW or CDW order. The theory hence excludes an extended metallic region (as opposed to a critical line) with K ρ < 1/2, even though this is difficult to verify numerically.
Previous work on the extended Hubbard model [81] suggests that a peak in K ρ (L) that narrows with increasing L indicates a continuous transition, whereas the absence of a peak or broadening with increasing L signals a first-order transition. Figures S1(e)-(h) hence support continuous behavior, in accordance with theoretical expectations [1,38].

Stiffness fits
Standard BKT universality is predicted for the LL-BOW and LL-CDW transitions both in a general LL [44] and specifically for the frustrated XXZ chain [38]. A detailed RG analysis [82] gives the finite-size scaling forms which provide the leading corrections to Eq. (4). However, in the light of the observed nonuniversal jumps, functional RG predictions of K ρ < 1/2 at the LL-BOW transition [36], and K ρ < 1/2 at the BOW-CDW transition according to our data and theory [38], we determined the critical values in Fig. 1 using fits based on Eq. (4) with three parameters: D ρ (∞), g, and C. In contrast, g and D ρ (∞) can be computed exactly for the classical 2D XY model (see below), leaving only one free parameter. Specifically, for βt ∼ L y = ∞ (1D quantum chain at T = 0), g = 1 and D ρ = 2/π (D ρ = t/2) for the 2D XY (1D t-V ) model [55,83]. As expected and demonstrated below, multi-parameter fits provide less accurate, but nonetheless fully consistent, critical values (shallower minima, stronger dependence on the range of L) than single-parameter fits. This is particularly relevant for the analysis of quantum systems such as the SSH model, where the range and number of system sizes are limited.
For the fits, we restricted the range of the jump to 0 < D ρ (∞) < 2t/π, using the known value of the noninteracting case. To discriminate between the logarithmic scaling at the critical point and the very weak finite-size dependence at weak coupling [see Fig. 4(a)], a nonzero lower bound g min was imposed. Otherwise, the choice g = 0 gives good fits throughout the LL phase and there would be no minimum of χ 2 ν at the critical point. The exact value of g min does not significantly affect the results and was chosen as 0.25. Finally, the allowed range of C was [0, ∞[. An important test case for the generalized, multi-parameter fit ansatz (4) was the LL-CDW transition of the t-V model, for which the critical value is known. We used the same range of system sizes as for the SSH model. Figures S2(a)-(c) give a comparison of results based on Eqs. (4), (S14), and (S15). All three fit functions yield very similar and hence compatible minima of χ 2 ν at the correct value λ = 2. Figures S2(d)-(f) are based on fits that exploit the known values g = 1 and D ρ (∞) = t/2. This additional information produces significantly sharper minima, in accordance with previous work on 2D XY models [54]. At the same time, the first-order fit functions (4) and (S14) do not fully capture the finite-size scaling on small system sizes, as manifested in χ 2 ν 1 even at λ = 2 in Figs. S2(d) and (e) for L ≥ 22 and L ≥ 30. Higher-order corrections are partially captured by varying g and D ρ (∞) [55], which explains the much better χ 2 ν for the same range of L in Figs. S2(a) and (b). For the more challenging case of t/ω 0 > 0, we focused on three-parameter fits based on Eqs. (4) and Eq. (S14). Within the present accuracy, the results are compatible with each other but slightly less systematic than for the t-V model. In particular, the fits become less robust upon increasing the smallest value of L due to a reduced number of degrees of freedom. A similar picture arises for a fixed λ = 4 in Figs. S3(e) and (f). For the present accuracy and range of system sizes, we cannot discriminate between the scaling forms (4), (S14), and (S15).

Nonuniversal stiffness jumps
For the XY model, the stiffness jump and the constant g can be computed from a series for a given aspect ratio r = L x /L y [55,84]. For 1+1D quantum systems, r = cL/β, with c a model-dependent constant. For example, D ρ (∞) varies significantly as a function of L x /L y , covering the whole range from 2/π to 0 [84]. Similarly, g varies between 1 and ∞ as a function of r. In principle, a change of the range of the retarded interaction can mimic a change in r, leading to dependence of D ρ (∞) and g on the phonon frequency.
There are several other known mechanisms for nonuniversal values of the stiffness jump. The bosonization relation D ρ = K ρ u [44] implies that, even if K ρ = 1/2 at a QPT, D ρ (∞) can change via the renormalized velocity u. For example, u increases with V in the t-V model [37] but decreases with λ in the Holstein model [79]. The stiffness can also be reduced by non-vortex excitations that are not captured by the standard BKT theory of the XY model [85].