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

TL;DR

This paper investigates the phase diagram of the 1D Su-Schrieffer-Heeger model with spinless fermions and quantum phonons, revealing new phases and unconventional critical behavior through quantum Monte Carlo simulations.

Contribution

It provides a revised phase diagram including an extended charge-density-wave phase and links to deconfined quantum criticality, differing from previous studies.

Findings

Discovery of an extended charge-density-wave phase

Identification of a retardation-driven phase transition

Observation of unconventional power-law exponents at criticality

Abstract

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. In addition to Luttinger liquid and bond-order-wave (BOW) phases, we find an extended charge-density-wave (CDW) phase. Because of different broken symmetries, BOW and CDW phases are connected by a retardation-driven phase transition. Our results are consistent with the theory of the frustrated $XXZ$ chain, including unconventional power-law exponents at criticality, and an interpretation in terms of deconfined quantum criticality via proliferation of solitons.

Figures

Equations

$$\hat{H}$$

$$S_{\text{ret}}=-\frac{\lambda t}{2}{\iint }_0^{\beta }\text{d}\tau \text{d}{\tau }^{\prime }\sum _b{B}_b(\tau )P(\tau -{\tau }^{\prime })B_b({\tau }^{\prime }).$$

$$C_{\rho }(r)$$

$$C_{\text{b}}(r)$$

$$C_{\rho }(r)$$

$$C_{\text{b}}(r)$$

$$\displaystyle D_{\rho}=L\,\frac{\partial^{2}E(\phi)}{\partial\phi^{2}}\bigg{\rvert}_{\phi=0}\,,$$

$$\frac{D_{\rho }(L)}{D_{\rho }(\infty )}=1+\frac{g}{2\ln L+C}.$$

$${\chi }_{\alpha }(\pi )=\sum _r(-1{)}^r{\int }_0^{\beta }d\tau C_{\alpha }(r,\tau ).$$

$$\mathcal{H}=\frac{v}{2}\left[\frac{1}{\eta }({\partial }_x\theta {)}^2+\eta ({\partial }_x\varphi {)}^2+{\lambda }_{\varphi }\cos (\sqrt{8\pi }\varphi )\right].$$

$$\mathcal{L}=\overline{\mathbf{\psi }}\text{i}{\gamma }^{\mu }{\partial }_{\mu }\mathbf{\psi }+g_{\text{bow}}(\overline{\mathbf{\psi }}\mathbf{\psi }{)}^2+g_{\text{cdw}}(\overline{\mathbf{\psi }}\text{i}{\gamma }_5\mathbf{\psi }{)}^2,$$

$$Z=\sum _{n=0}^{\infty }\frac{1}{n!}\sum _{C_n}{w}_{{\nu }_1}\dots w_{{\nu }_n}\int \mathcal{D}(\bar{c},c)e^{-{\mathcal{S}}_0}{h}_{{\nu }_1}\dots h_{{\nu }_n}$$

$${\mathcal{S}}_1=-{\iint }_0^{\beta }d{\tau }_{1}d{\tau }_2{P}_{+}({\tau }_1-{\tau }_2)\sum _{a_1,a_2,b}{h}_{a_1{a}_2,b}({\tau }_1,{\tau }_2).$$

$$h_{11,b}({\tau }_1,{\tau }_2)=\frac{\lambda t}{2}{B}_b({\tau }_1)B_b({\tau }_2)$$

$$h_{10,b}({\tau }_1,{\tau }_2)$$

$$h_{01,b}({\tau }_1,{\tau }_2)$$

$$h_{00,b}({\tau }_1,{\tau }_2)=k{1}_b({\tau }_1)1_b({\tau }_2).$$

$$C_{\text{b}}(b_1,b_2;C_n)=\frac{1}{\beta t^2}\sum _p{I}_{b_1{b}_2}(p-1,p)K(p-1,p).$$

$$\displaystyle\rho_{\text{s}}=L\,\frac{\partial^{2}F(\phi)}{\partial\phi^{2}}\bigg{\rvert}_{\phi=0}\,.$$

$${\mathcal{S}}_{\varphi }={\mathcal{S}}_{\text{L}}{e}^{\mathrm{i}\varphi }+{\mathcal{S}}_{\text{R}}{e}^{-\mathrm{i}\varphi }+{\mathcal{S}}_{\text{LL}}{e}^{2\mathrm{i}\varphi }+{\mathcal{S}}_{\text{RR}}{e}^{-2\mathrm{i}\varphi }.$$

$$\displaystyle\rho_{\text{s}}=\frac{L}{\beta}\bigg{[}\left\langle\frac{\partial\mathcal{S}_{\phi}}{\partial\phi}\right\rangle^{2}+\bigg{\langle}{\frac{\partial^{2}\mathcal{S}_{\phi}}{\partial\phi^{2}}}\bigg{\rangle}-\bigg{\langle}{\bigg{(}\frac{\partial\mathcal{S}_{\phi}}{\partial\phi}\bigg{)}^{2}}\bigg{\rangle}\bigg{]}\bigg{\rvert}_{\phi=0}\,.$$

$$\displaystyle\bigg{\langle}\frac{\partial\mathcal{S}_{\phi}}{\partial\phi}\bigg{\rangle}\bigg{\rvert}_{\phi=0}$$

$$=-\mathrm{i}⟨\left(n_{\text{L}}-n_{\text{R}}\right)+2\left(n_{\text{LL}}-n_{\text{RR}}\right)⟩.$$

$$\displaystyle\bigg{\langle}\frac{\partial^{2}\mathcal{S}_{\phi}}{\partial\phi^{2}}\bigg{\rangle}\bigg{\rvert}_{\phi=0}$$

$$=⟨\left(n_{\text{L}}+n_{\text{R}}\right)+4\left(n_{\text{LL}}+n_{\text{RR}}\right)⟩$$

$$\displaystyle\bigg{\langle}{\bigg{(}\frac{\partial\mathcal{S}_{\phi}}{\partial\phi}\bigg{)}^{2}}\bigg{\rangle}\bigg{\rvert}_{\phi=0}=$$

$$=$$

$$+⟨\left(n_{\text{L}}+n_{\text{R}}\right)+4\left(n_{\text{LL}}+n_{\text{RR}}\right)⟩,$$

$${\rho }_s=\frac{L}{\beta }(⟨W^2⟩-⟨W{⟩}^2).$$

$$\frac{D_{\rho }(L)}{D_{\rho }(\infty )}=1$$

$$+\frac{a}{[2\ln L+C+\ln (C/2+\ln L){]}^2},$$

$${\chi }_{\nu }^2=\frac{1}{\nu }{\chi }^2=\frac{1}{\nu }\sum _{n=1}^N{\left(\frac{O_n-C_n}{{\sigma }_n}\right)}^2.$$

Full text

I Introduction

The rich phase diagrams of, e.g., dichalcogenides [1] or cuprates [2] motivate fundamental investigations of competing orders in strongly correlated quantum systems. Motivated by the significant complexity of real materials and justified by the concept of universality, most theoretical work is based on minimal and hence tractable models. A recent focus are Dirac systems [3; 4; 5; 6; 7; 8; 9]. The existence of two or more ordered phases also provides a route to study non-Landau deconfined quantum critical points (DQCPs), for which topological excitations of the order parameters play a central role [10]. For recent progress, see Refs. [11; 12; 13; 14; 9; 15; 16; 17; 18]. Different orders can arise either from different interactions, or from local and nonlocal components of the same Coulomb interaction [19]. The intricacy of such problems is reflected, e.g., in the debates surrounding the complex phase diagrams of extended Hubbard models on one-dimensional (1D) chains and 2D honeycomb lattices, see Refs. [20; 21] for reviews.

Retardation effects, which are negligible for Coulomb interactions, play a fundamental role in the context of electron-phonon coupling. As is known from the theory of superconductivity [22], phonon-mediated interactions have attractive and repulsive components. However, for a commensurate band filling, the different electron-phonon couplings are commonly associated with only a single type of order each. A charge-density-wave (CDW) state with a modulated electron site density follows from a Holstein coupling and is observed in molecular crystals [23]. A bond-order-wave (BOW) state with a dimerized kinetic energy emerges from a Su-Schrieffer-Heeger (SSH) coupling and is realized in conjugated polymers [23]. Both orders are illustrated in Fig. 1 for a 1D chain.

The SSH model [24] describes electrons coupled to quantum phonons. Originally introduced to study topological solitons in conjugated polymers [24], SSH models also have close relations with field theories of Dirac fermions [25]. The mean-field SSH model provides an important platform to explore interaction and nonequilibrium effects on 1D topological phases [26; 27; 28; 29; 30; 31; 32]. Novel experimental realizations include cold atoms [33] and resistor networks [34].

Despite the fundamental role of the SSH model, the retarded nature and momentum dependence of the phonon-mediated electron-electron interaction have so far prohibited a detailed understanding. Compared to simpler models of spins or fermions, theoretical and numerical approaches to electron-phonon problems face significant challenges [20]. Motivated by this situation, we explore the full phase diagram of the 1D spinless SSH model by quantum Monte Carlo (QMC) simulations and find important deviations from existing work and rich, previously undiscovered physics. We demonstrate the existence of not one [35] but two dimerized phases upon variation of the phonon frequency $\omega_{0}$: the familiar BOW Peierls phase with dimerized hopping and an extended CDW phase. This raises the central question of the nature of the BOW-CDW phase transition, which is entirely beyond the usual adiabatic and antiadiabatic approximations. Our findings, most notably unconventional exponents at criticality, are consistent with a field theory previously developed for frustrated spin chains [36]. They also establish 1D DQCP physics [16] in an electron-phonon model.

The paper is organized as follows. In Sec. II we define the SSH model and discuss its exact limits, in Sec. III we present our QMC results, in Sec. IV we discuss our results, and in Sec. V we conclude. The Appendices contain details on the QMC method and additional data.

II Model & exact limits

The Hamiltonian of the SSH model [24; 37],

[TABLE]

with $\hat{B}_{b}=\hat{c}^{\dagger}_{i(b)}\hat{c}^{\phantom{{\dagger}}}_{j(b)}+\hat{c}^{\dagger}_{j(b)}\hat{c}^{\phantom{{\dagger}}}_{i(b)}$ acting on bond $b$ between sites $i$ and $j=i+1$ and $\hat{H}_{\text{ph}}=\sum_{b}\left(\frac{1}{2M}\hat{P}^{2}_{b}+\frac{K}{2}\hat{Q}_{b}^{2}\right)$, describes spinless fermions coupled to optical bond phonons with momentum $\hat{P}_{b}$, displacement $\hat{Q}_{b}$, and frequency $\omega_{0}=\sqrt{K/M}$. For the present case of a half-filled band ($\langle\hat{n}_{i}\rangle=\langle\hat{c}^{\dagger}_{i}\hat{c}^{\phantom{{\dagger}}}_{i}\rangle=0.5$), numerics [38; 39] and field theory [35] suggest the same physics for optical phonons—which preclude a QMC sign problem—and the original acoustic phonons [24]. After integrating out the phonons, the partition function contains the retarded interaction

[TABLE]

The free phonon propagator $P(\tau)$ is local in space but its decay in imaginary time $\tau$ (here, $\beta=1/T$) is determined by $\omega_{0}$, with $P(\tau)\sim e^{-\omega_{0}\tau}$. The associated retardation effects are crucial for the phase diagram in Fig. 1. We use the coupling $\lambda=g^{2}/Kt$ and set $\hbar=k_{\text{B}}=1$.

For $\omega_{0}=0$, corresponding to classical phonons, mean-field theory is exact at $T=0$. Replacing $\hat{Q}_{b}$ with $\langle\hat{Q}_{b}\rangle=(-1)^{b}\Delta/g$ in Eq. (1) yields the fermionic hopping term $\hat{H}_{0}=-\sum_{b}[t+(-1)^{b}\Delta]\hat{B}_{b}$. The Peierls argument [40] implies that the bond dimerization $\Delta$ is nonzero for any $\lambda>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 $\lambda_{c}(\omega_{0})$ [35; 41; 38]. An exact solution (by the Bethe ansatz) is also possible in the opposite, antiadiabatic limit. For $\omega_{0}\to\infty$, the interaction (2) becomes instantaneous and Eq. (1) maps to the $t$-$V$ model $\hat{H}_{\infty}=-t\sum_{b}\hat{B}_{b}+V\sum_{i}\hat{n}_{i}\hat{n}_{i+1}$ with $V=\lambda$ [35] and a LL-CDW QPT at $V_{c}/t=\lambda_{c}=2$ [42].

BOW and CDW states, illustrated in Fig. 1, spontaneously break translation symmetry and are described by Ising order parameters that reflect the two possible BOW (CDW) dimerization patterns related by a shift by one lattice constant. The discrete Ising symmetry permits long-range order at $T=0$. CDW and BOW states can be distinguished by point group symmetries. CDW order breaks bond inversion symmetry but preserves invariance under site inversion, see Fig. 1. The opposite is true for the BOW phase. Given the different broken symmetries, a phase transition is expected between $\omega_{0}=0$ and $\omega_{0}=\infty$. In contrast, the influential Ref. [35] suggests a single dimerized phase mainly based on a continuum theory, details and extensions of which are discussed below. Similar conclusions were reached by (functional) RG calculations [43; 41] and for related spin-phonon models [44; 45; 46; 41]. A nonadiabatic mean-field approach yields BOW and CDW phases even for large $\omega_{0}$ [47]. These findings differ significantly from our results. Finally, the choice of spinless fermions is motivated by the absence of a LL phase [20] and only one symmetry-broken (BOW) phase in the spinful case; at small $t/\omega_{0}$, critical spin correlations are expected.

III Results

Our simulations were made possible by using a state-of-the-art QMC method based on the stochastic series expansion [48] and directed-loop updates [49]. The performance gain from extending the latter to retarded interactions [50] is essential to explore the phase diagram of the SSH model. The method has only statistical errors and relevant technical details are summarized in App. A. Results were obtained for periodic chains of $L$ sites, and for inverse temperatures $\beta t=2L$ representative of $T=0$.

III.1 Identification of the different phases

The three distinct regimes in Fig. 1 are revealed by the real-space correlation functions. Figure 2 shows the charge ($\rho$) and bond (b) correlators

[TABLE]

Only positive values are visible in the logarithmic representation. Using the conformal distance $x=L\sin\left(\pi r/L\right)$ minimizes the effects of the periodic boundaries. For $t/\omega_{0}=1/3$ and $\lambda=2$ [Fig. 2(a)], both exhibit a power-law decay of $q=2k_{\text{F}}$ correlations described by the LL expressions for 1D metals [51]

[TABLE]

with $K_{\rho}\approx 0.61$. At $\lambda=4$, corresponding to the BOW phase in Fig. 1, the saturation of $C_{\text{b}}(x)$ at large $x$ and the exponential decay of $C_{\rho}(x)$ are consistent with long-range bond order. In the nonadiabatic regime [$t/\omega_{0}=1/10$, Fig. 2(b)], we find behavior consistent with a LL at $\lambda=2$ and long-range CDW order at $\lambda=5$.

A general diagnostic to distinguish metallic and insulating 1D phases is the charge stiffness

[TABLE]

which measures the response of the system—via the ground-state energy $E$—to a magnetic flux $\phi$ [52]. In our simulations, it can be obtained from the estimator of the superfluid stiffness at low enough temperatures, see App. B. For $L\to\infty$, the charge stiffness is nonzero (zero) in metallic (insulating) phases [53; 54]. To illustrate the topology of the phase diagram, we show in Fig. 3 the stiffness $D_{\rho}(L)$ along a continuous path in parameter space: (a) we first increase $\lambda$ from 1 to 6 at fixed $t/\omega_{0}=1/10$, (b) then we tune $t/\omega_{0}$ from $1/10$ to $1/3$ at fixed $\lambda=6$, and (c) finally we decrease $\lambda$ again from 6 to 1 at fixed $t/\omega_{0}=1/3$. In accordance with Fig. 1, $D_{\rho}(L)$ remains nonzero for $\lambda<3$ (LL phase) and clearly vanishes for $L\to\infty$ for $\lambda\geq 4$ (CDW phase) at fixed $t/\omega_{0}=1/10$, see Fig. 3(a). The determination of critical values and the finite-size scaling at and near $\lambda_{c}$ ($\lambda_{c}\approx 3$ according to Fig. 1) will be discussed in detail below [see also Fig. 4(a)]. Similar behavior is observed in Fig. 3(c) for the LL-BOW transition at $t/\omega_{0}=1/3$. Remarkably, at fixed $\lambda=6$, the stiffness converges to a nonzero value only at intermediate $t/\omega_{0}$, see Fig. 3(b). It is shown in App. C.1 that the peak in $D_{\rho}(L)$ narrows with increasing $\lambda$ but its maximum value remains finite. Naively, one would assume that the insulating BOW and CDW phases are separated by an extended metallic region. However, numerical results and field-theory arguments discussed below provide evidence for a metallic line separating the ordered phases, as suggested by Fig. 1.

III.2 Finite-size scaling of the charge stiffness

To obtain the LL-BOW and LL-CDW critical values in Fig. 1, we analyzed the finite-size scaling of $D_{\rho}(L)$, see Fig. 4(a). For 1D metal-insulator transitions, a renormalization-group (RG) analysis of umklapp interactions predicts numerically challenging Berezinskii-Kosterlitz-Thouless (BKT) scaling with a critical value $K_{\rho}=1/2$ [51]. This scenario has been explicitly confirmed for $t/\omega_{0}=0$ [42]. Interestingly, a functional RG study of the SSH model [41] reported unconventional BKT physics with $K_{\rho}<1/2$. Large-scale simulations of classical frustrated 2D $XY$ models [55; 56]—relevant due to the usual quantum-classical mapping and a connection between retardation and frustration explained below—indicate a standard BKT transition [57; 55] (see, however, Ref. [58]), albeit with challenging crossover phenomena.

The RG calculations predict a characteristic logarithmic scaling exactly at a critical point due to marginally relevant operators. To first order [59],

[TABLE]

As demonstrated before for 2D classical $XY$ models [59; 60] and 1D quantum models [61; 62], critical values can be extracted by fitting the stiffness data to Eq. (6), treating the jump $D_{\rho}(\infty)$ [63] as well as $g$ and $C$ as parameters. Because the form (6) only holds at critical points, the latter can be identified as the minima in standard goodness-of-fit measures, such as the reduced chi-squared $\chi_{\nu}^{2}=\chi^{2}/\nu$ for $\nu$ degrees of freedom used here. Technical details of the stiffness fits are discussed in App. D. The universal stiffness jump at the critical point of the 2D $XY$ model [63] translates to $D_{\rho}(\infty)=t/2$ for the $t$-$V$ model (the SSH model with $t/\omega_{0}=0$) [42]. For $t/\omega_{0}>0$, we instead find nonuniversal stiffness jumps, with $D_{\rho}(L)<t/2$ even for small $L$ in, e.g., Fig. 5(a). Possible origins are discussed in App. E.

For $t/\omega_{0}=0$ (the $t$-$V$ model), such fits indeed yield a minimum of $\chi^{2}_{\nu}$ at the exact critical value $\lambda_{\text{c}}=2$ in Fig. 4(b). For $t/\omega_{0}=1/10$, we estimate $\lambda_{\text{c}}=3.0(3)$ from Fig. 4(c). Finite-size effects are larger for a smaller $\omega_{0}$ (longer interaction range in $\tau$), as visible in Fig. 4(d) for $t/\omega_{0}=1/3$ in terms of a larger minimal $\chi^{2}_{\nu}$; the critical value is determined as $\lambda_{\text{c}}=3.1(3)$. At fixed $\lambda=4$ [Fig. 4(d)], the stiffness fits exhibit two separate minima at $t/\omega_{0,\text{c1}}=0.15(3)$ and $t/\omega_{0,\text{c2}}=0.24(2)$, respectively, consistent with two critical points. For $\lambda=6$ and $8$, due to restrictions in system size, we estimated a single critical value from the peaks in $D_{\rho}(L)$ and $K_{\rho}(L)$ (see Fig. 5 and App. C) as $t/\omega_{0,\text{c}}=0.16(2)$ and $t/\omega_{0,\text{c}}=0.145(15)$, respectively. These critical values appear consistent with the finite-size behavior of other observables, as discussed in App. C.

III.3 CDW-BOW transition

Results for the retardation-driven CDW-BOW transition, a key feature of the SSH model, are presented in Fig. 5 for $\lambda=6$. The convergence of $D_{\rho}(L)$ in Fig. 5(a) and $K_{\rho}(L)=2\pi C_{\rho}(q_{1})/q_{1}$ (with $q_{1}=2\pi/L$) in Fig. 5(b) to a nonzero value around $t/\omega_{0,\text{c}}\approx 1/6$ indicates a metallic state. However, $K_{\rho}(L)$ falls below the typical lower bound $1/2$ for a LL [51]. The dotted line in Fig. 5(b) is $K_{\rho}=\overline{\eta}/2$, with $\overline{\eta}=0.79$ the exponent of the real-space charge and bond correlations at $t/\omega_{0}=1/6$ in Figs. 5(c) and 5(e). Combining the LL form of $C_{\alpha}(r)$, $\alpha\in\{\rho,\mathrm{b}\}$, given by Eq. (4) with conformal invariance implies $\chi_{\alpha}(\pi)/L\sim L^{1-\overline{\eta}}$ [64] for the staggered charge/bond susceptibility

[TABLE]

That this relation holds at $t/\omega_{0}=1/6$ [Figs. 5(d) and 5(f)] is evidence for a gapless state described by a conformal field theory. In particular, charge and bond susceptibilities have the same power-law exponent (within the available accuracy) and do not show any visible finite-size corrections even down to the smallest system sizes considered. Away from $t/\omega_{0}=1/6$, Figs. 5(c)–5(f) suggest either BOW or CDW order with one of the orders instantly becoming dominant and the other being suppressed. Our findings provide evidence for a direct quantum phase transition between CDW and BOW order with an intermediate metallic point instead of a narrow region. This interpretation is further supported by theoretical arguments given in Sec. IV. Data on the CDW-BOW transition for other $\lambda$ can be found in App C.

IV Discussion

The retardation effects that drive the BOW-CDW transition are generally difficult to capture analytically [65; 41; 39]. However, we find our numerical results to be fully consistent with a field theory previously derived for the antiferromagnetic $J_{1}$–$J_{2}$ $XXZ$ chain with Hamiltonian $\hat{H}_{J_{1}-J_{2}}=J_{1}\hat{H}_{XXZ}+J_{2}\sum_{i}{\bm{S}}_{i}{\bm{S}}_{i+2}$ and $\hat{H}_{XXZ}=\sum_{i}(\hat{S}^{x}_{i}\hat{S}^{x}_{i+1}+\hat{S}^{y}_{i}\hat{S}^{y}_{i+1}+{\Delta_{z}}\hat{S}^{z}_{i}\hat{S}^{z}_{i+1})$ [36]. Its bosonized continuum description takes the usual Tomonaga-Luttinger form [36; 16]

[TABLE]

Here, $v$ is the renormalized velocity, $\eta=1/\overline{\eta}=1/(2K_{\rho})$ determines the exponent of correlation functions, and the cosine term encodes umklapp scattering. The phase diagram of $\hat{H}_{J_{1}-J_{2}}$ has the same topology as Fig. 1 [36]; the Néel and dimer phases correspond to CDW and BOW phases, respectively.

A connection between the frustrated $XXZ$ model $\hat{H}_{J_{1}-J_{2}}$ and the SSH model (1) can be established via a mapping of phonon-mediated retarded interactions at small $t/\omega_{0}$ to frustrated spin interactions [66; 67; 46; 68], but is expected to contain additional terms not present in $\hat{H}_{J_{1}-J_{2}}$. Alternatively, a Hamiltonian of the form (8) can in principle be obtained for the SSH model via an RG treatment of the phonon-mediated interaction. In either case, in contrast to the $XXZ$ problem [36; 16], an explicit expression for $\lambda_{\phi}$ in terms of the SSH model parameters is not known. Here, we focus on comparing the predictions for Eq. (8) [36; 16] to our numerical data.

The LL-BOW and LL-CDW transitions can be attributed to the cosine umklapp term in Eq. (8) becoming relevant for $K_{\rho}<1/2$. In contrast to approximate functional RG results [41], we do not find evidence for a violation of the conventional $K_{\rho}=1/2$ at the critical point (see App. C.2 for details). BOW and CDW phases are associated with opposite signs of $\lambda_{\phi}$ and a pinning of the charge mode $\phi$ at different minima [51; 16]. The CDW phase is known in the Gross-Neveu literature as an Aoki phase [69]. The line of BOW-CDW transitions suggested by Fig. 1 for $\lambda\gtrsim 6$ mirrors the line of continuous dimer-Néel transitions of the frustrated $XXZ$ chain, along which $K_{\rho}$ varies continuously [36; 16]. Within the theory (8), and given $K_{\rho}<1/2$ [see Fig. 5(b)], a metallic state separating BOW and CDW phases requires $\lambda_{\phi}=0$. Since $\lambda_{\phi}$ depends on the independent SSH model parameters $t/\omega_{0}$ and $\lambda$, it will vanish at a single value $(t/\omega_{0})_{\text{c}}$ for a given $\lambda$. This scenario is fully supported by our observation that for a given $\lambda\gtrsim 6$, bond and charge correlations show exactly the same power-law exponents at a single $(t/\omega_{0})_{\text{c}}$ in Figs. 5(c)–5(f). For other $t/\omega_{0}$, umklapp scattering immediately gives rise to long-range order [Figs. 5(c)–5(f)], $K_{\rho}\to 0$ [Fig. 5(b)], and insulating behavior [Fig. 5(a)]. The clean power-law scaling observed on the BOW-CDW critical line in Figs. 5(c)–5(f) is consistent with the predicted absence of logarithmic corrections due to the vanishing of the umklapp term [70]. Finally, the absence of a CDW phase in functional RG results for the SSH model [41] may be the result of using a linearized spectrum. This approximation amounts to an infinite bandwidth $W$, whereas CDW order appears in Fig. 1 at $\omega_{0}\gtrsim W=4t$.

An intuitive physical picture of how BOW and CDW phases—breaking different symmetries—can be connected via a generically continuous phase transition is the scenario of a 1D DQCP [16]. It involves 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 $\bm{\varphi}=(\cos\varphi,\sin\varphi)$ [16], see the inset of Fig. 1, BOW (CDW) patterns correspond to $\varphi=0,\pi$ ($\varphi=\pm\pi/2$). For example, a defect in the BOW order connecting $\varphi=0,\pi$ contains a region with CDW order or $\varphi=\pi/2$. Simultaneous proliferation of BOW/CDW defects at $(t/\omega_{0})_{\text{c}}$ provides a mechanism for a continuous transition without fine-tuning.

Instead of the bosonized theory (8), the SSH model (1) can also be described in terms of a Gross-Neveu field theory of Dirac fermions [71]. Such a representation makes topological and symmetry aspects more transparent. While a Gross-Neveu theory was given in Ref. [35], our numerical results can only be captured by the more general form with two interactions [71]

[TABLE]

where ${\bm{\psi}}=(\psi_{A},\psi_{B})$ ($A$/$B$: sublattices, see Fig. 1). Similar to the umklapp term in Eq. (8), the interactions account for the lattice symmetries that distinguish BOW and CDW order [72]. In contrast, a mean-field approximation of Eq. (9) contains mass terms $g_{\text{bow}}m_{1}\overline{\bm{\psi}}{\bm{\psi}}$ and $g_{\text{cdw}}m_{2}\overline{\bm{\psi}}\,\text{i}\gamma_{5}{\bm{\psi}}$ [73]. It can be shown [25] that these mass terms anticommute and may be rotated into each other, establishing a chiral U(1) symmetry absent beyond the mean-field level. Combining the BOW and CDW masses into a vector $\bm{m}=(m_{1},m_{2})$, the mean-field spectrum $E(p)=\pm\sqrt{p^{2}+|\bm{m}|^{2}}$ [74]. Together with the chiral symmetry this form reveals the possibility of a continuous evolution between BOW and CDW order during which $\bm{m}$ changes its direction from $(m,0)$ to $(0,m)$ while $|\bm{m}|=m$ (the mean-field gap) remains nonzero. The absence of a CDW-BOW phase transition at the mean-field level reflects the fact that the CDW and BOW order parameters (characterized by different broken symmetries on the lattice, see Sec. II) become equivalent up to a chiral rotation in the continuum. The different symmetries are correctly captured by the bosonized theory (8) and also the Gross-Neveu theory (9) [72; 25]. In the absence of a chiral symmetry, the CDW-BOW transition involves a gap closing, $|\bm{m}|\to 0$ at the critical point. This agrees with our numerical results for the spinless SSH model (1), for which metallic behavior entails a vanishing single-particle gap.

Finally, it is interesting to contrast the BOW-CDW transition considered here with recent work on interaction-driven QPTs out of a topological band insulator [75; 31; 76] (see also Ref. [20]). In the latter, a static BOW mass term arises from a dimerized hopping (the mean-field SSH model). Soliton excitations are therefore only possible for the CDW order parameter and the critical behavior is significantly different. Instead of the deconfined scenario observed here, the BOW-CDW transition exhibits Ising criticality [31].

V Conclusions & outlook

We used an exact QMC method for retarded interactions to determine the phase diagram of the 1D SSH model with quantum phonons. In addition to the well-known metallic LL and the Peierls-ordered BOW phase we found an extended CDW phase at high phonon frequencies $\omega_{0}$ that has been absent in previous studies [35; 43; 41; 38]. We provided evidence that the CDW and BOW ordered phases are connected by a direct, continuous quantum phase transition with unconventional power-law exponents $K_{\rho}<1/2$ at the metallic critical point. Our findings are consistent with a bosonized field theory that was originally introduced for the frustrated $XXZ$ chain [36]. In analogy with frustrated spin systems [16], the CDW-BOW transition can be interpreted in terms of a 1D DQCP and the proliferation of solitons.

Our results demonstrate that competing orders can be generated from a single retarded interaction that originates, e.g., from an off-diagonal operator coupled to a phonon. Instead of having different competing interactions in an equal-time Hamiltonian, here the interaction range in imaginary time determines whether BOW or CDW order dominates. Competing phenomena therefore arise at a critical interaction range determined by the phonon frequency $\omega_{0}$. The study of retarded interactions might be a promising approach to generate complex phase diagrams also in higher dimensions. An interesting generalization of the work presented in this paper is the spinful 2D SSH model, which supports a phase transition between valence-bond and antiferromagnetic phases [77]. Moreover, the relation between the spinless SSH model and the still incompletely understood frustrated $XY$ model [57; 55; 58] should be explored.

Acknowledgements.

We acknowledge helpful discussions with F. Assaad, A. Furusaki, C. Mudry, and F. Pollmann. M.W. was supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES) under Award DE-FG02-08ER46542. F.P.T was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project number 414456783, M.H. via SFB 1170. We thank the Gauss Centre for Supercomputing (SuperMUC at the Leibniz Supercomputing Centre) for generous allocation of supercomputing resources.

Appendix A Quantum Monte Carlo method

We used the directed-loop QMC method for retarded interactions in the path-integral representation [50]. It is based on an interaction expansion of the partition function $Z=\int\mathcal{D}(\bar{c},c)\,e^{-\mathcal{S}_{0}-\mathcal{S}_{1}}$ around $\mathcal{S}_{0}=\int d\tau\sum_{i}\bar{c}_{i}(\tau)\,\partial_{\tau}\,c_{i}(\tau)$. A general interaction vertex $\mathcal{S}_{1}=-\sum_{\nu}w_{\nu}\,h_{\nu}$ can be written as a sum over vertex variables $\nu$, a weight $w_{\nu}$, and the Grassmann fields contained in $h_{\nu}$. The perturbation expansion becomes

[TABLE]

with sums over the expansion order $n$ and the ordered vertex list $C_{n}=\{\nu_{1},\dots,\nu_{n}\}$. For each time-ordered configuration of vertices, the expectation value over Grassmann fields can be represented by world lines. The trivial choice of $\mathcal{S}_{0}$ ensures that the imaginary-time evolution is entirely determined by the interaction vertices. Therefore, Eq. (10) is the path-integral equivalent of the stochastic series expansion (SSE) representation where $Z=\mbox{Tr}\,e^{-\beta H}$ is expanded in the total Hamiltonian [48; 78]. Accordingly, many algorithmic features, including the global directed-loop updates [49], directly transfer to the path-integral representation [78].

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\in\{1,2\}$ with local variables $\{a_{j},b_{j},\tau_{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

[TABLE]

It is important to note that the symmetrized phonon propagator $P_{\!+}(\tau)=\omega_{0}\cosh[\omega_{0}(\beta/2-\tau)]/[2\sinh(\omega_{0}\beta/2)]$ is included in the global weight $w_{\nu}$ of the vertex. Whereas the bond-bond interaction

[TABLE]

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.,

[TABLE]

This is possible because $\int_{0}^{\beta}d\tau_{2}\,P_{\!+}(\tau_{1}-\tau_{2})=1$. As the vertices (12) and (13) 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,

[TABLE]

With our choice of interaction vertices, we can formulate the diagonal and directed-loop updates similar to the SSE representation [49]. For the diagonal updates, we use the Metropolis algorithm to add and remove vertices $h_{00,b}(\tau_{1},\tau_{2})$ that do not change the world-line configurations but change the expansion order $n$. We propose time differences $\tau_{1}-\tau_{2}$ according to the phonon propagator using inverse-transform sampling. Because $P_{\!+}(\tau_{1}-\tau_{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. [50]. The constant $k$ in Eq. (14) 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}\leftrightarrow h_{10,b}/h_{01,b}\leftrightarrow 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. [49; 50].

The calculation of observables in the path-integral (interaction) representation is in many ways similar to the SSE representation. Sandvik et al. [78] 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_{\rho}(r)$ and the charge susceptibility $\chi_{\rho}(r)$, are simple to derive and given in Ref. [78]. 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 $\chi_{\text{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. [78] for the exact estimator. However, calculating the equal-time bond structure factor $C_{\text{b}}(r=b_{1}-b_{2})=\langle B_{b_{1}}B_{b_{2}}\rangle$ in the interaction representation is more involved. While a general derivation is outlined in Ref. [78], 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

[TABLE]

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 $\mathbb{1}_{b}$ as they were only introduced to simplify the Monte Carlo sampling. $I_{b_{1}b_{2}}(p-1,p)$ is zero unless bond operators $B_{b_{1}}(\tau_{p-1})$ and $B_{b_{2}}(\tau_{p})$ originating from the hopping terms $h_{10}/h_{01}$ appear at adjacent times; then $I_{b_{1}b_{2}}(p-1,p)=1$. An integral expression for $K(p-1,p)$ was derived in Ref. [78] and gives $K(p-1,p)=2/(\tau_{p+1}-\tau_{p-2})$ when 4 or more subvertices are present in a world-line configuration. The time difference $\tau_{p+1}-\tau_{p-2}\in[0,\beta]$ is defined by the two subvertices that surround the two bond operators under consideration. Note that $K(p-1,p)=2/\beta$ for 3 subvertices, $K(p-1,p)=1/\beta$ for 2 subvertices, and $K(p-1,p)=0$ for 0 or 1 subvertices. For further details, see Ref. [78].

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. [79; 80] that make use of the vertex distribution. In the following section, we use the framework outlined in Ref. [79] to show that the superfluid stiffness of an electron-phonon model can still be calculated from the winding number.

Appendix B Estimator for the superfluid stiffness

Consider a ring of length $L$ threaded by a magnetic flux $\phi$. At finite temperatures, the superfluid stiffness can be obtained from the free energy via [53]

[TABLE]

Because we study a 1D system [81] and our simulations at $\beta t=2L$ are essentially converged with respect to temperature, the measured values of $\rho_{\text{s}}$ are representative of the charge stiffness or Drude weight defined in Eq. (5).

Using $F=-\frac{1}{\beta}\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 $\hat{c}_{L+1}=e^{\mathrm{i}\phi}\hat{c}_{1}$. The boundary term of the action reads

[TABLE]

Here, $\mathcal{S}_{\text{L/R}}$ is the action of the hopping term (13) crossing the boundary to the left/right, whereas $\mathcal{S}_{\text{LL/RR}}$ corresponds to the bond-bond interaction (12) with both hopping operators going to the left/right. The superfluid stiffness can then be calculated as

[TABLE]

The first expectation value is given by

[TABLE]

For each Monte Carlo configuration, expectation values of terms $\mathcal{S}_{a}$ contained in the interaction vertex (11) can be obtained by counting the number of vertices $n_{a}$ [79]. For the Monte Carlo average we then obtain $\langle\mathcal{S}_{a}\rangle=-\langle n_{a}\rangle$. In the same way, the second term in Eq. (18) becomes

[TABLE]

and the third term is given by

[TABLE]

where we used $\langle\mathcal{S}_{a}\mathcal{S}_{b}\rangle=\langle n_{a}n_{b}\rangle-\delta_{ab}\langle n_{a}\rangle$. We get an additional shift for $a=b$ that cancels the contribution of (20). Our results are equivalent to calculating the winding number $W=n^{B}_{\text{L}}-n^{B}_{\text{R}}$ where $n^{B}_{\text{L/R}}$ counts the number of subvertices $B_{b}(\tau)$ crossing the boundary to the left/right. Here, $n_{\text{LL/RR}}$ contributes with a factor of 2 because each vertex contains two bond operators, whereas mixed contributions $n_{\text{LR}}$ drop out. Therefore, $\rho_{\text{s}}$ can be calculated in the same way for retarded interactions as for equal-time interactions [82], i.e.,

[TABLE]

Appendix C Additional data

C.1 CDW-BOW transition

Figure 6 shows $D_{\rho}(L)$ and $K_{\rho}(L)$ as a function of $t/\omega_{0}$ for $\lambda=4,6,8,12$. For all couplings, the data are consistent with a metallic region at intermediate $t/\omega_{0}$. Whereas the apparent narrowing of this region between $\lambda=4$ and $\lambda=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_{\rho}<1/2$. Values $K_{\rho}<1/2$ can be reconciled with metallic behavior by assuming $\lambda_{\phi}=0$ in Eq. (8) at the BOW-CDW critical point [36; 72]. Apart from $\lambda=6$, see Fig. 5(b), we also find evidence for $K_{\rho}<1/2$ at criticality for $\lambda=4$ [Fig. 6(e)], 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 the challenging finite-size scaling in the tricritical region of the phase diagram.

In contrast to Ginzburg-Landau theory, the BOW-CDW transition does not require fine-tuning of both $t/\omega_{0}$ and $\lambda$. For a fixed $\lambda$, $\lambda_{\phi}$ can be tuned to zero for a suitable value of $t/\omega_{0}$, giving rise to a line of critical points. Since $K_{\rho}<1/2$ at criticality, any nonzero $\lambda_{\phi}$ yields long-range BOW or CDW order. The theory hence excludes an extended metallic region (as opposed to a critical line) with $K_{\rho}<1/2$.

Previous work on the extended Hubbard model [83] suggests that a peak in $K_{\rho}(L)$ that narrows with increasing $L$ indicates a continuous transition, whereas the absence of a peak or a broadening with increasing $L$ signals a first-order transition. Figures 6(e)–(h) hence support continuous behavior, in accordance with theoretical expectations [36; 16].

C.2 LL-BOW and LL-CDW transition

The stiffness fits are based on the characteristic logarithmic scaling at the critical point. As is commonly done for non-BKT transitions in higher dimensions, the LL-BOW and LL-CDW critical phase boundaries can also be estimated from the real-space correlation functions in Eq. (3), evaluated at the maximum distance $L/2$. Results for $C_{\mathrm{b}}(L/2)$ at $t/\omega_{0}=1/3$ and $C_{\rho}(L/2)$ at $t/\omega_{0}=1/10$ are shown in Figs. 7(a) and 7(c), respectively. For both cases shown, the critical values can be identified visually as $\lambda_{\mathrm{c}}\approx 2.75$. They agree with our estimates from the stiffness fits within error bars, but suggest the possibility that the latter slightly overestimate the critical couplings. We notice that logarithmic corrections are expected to appear close to the transition. A similar analysis was done for the $t$-$V$ model in Ref. [64].

The critical points of the LL-BOW and LL-CDW phase transitions can further be deduced from the bond and charge susceptibilities defined in Eq. (7). In the LL phase, $\chi_{\mathrm{b}/\rho}(\pi)/L\sim L^{1-\overline{\eta}}$ with $\overline{\eta}=2K_{\rho}$. If we assume $K_{\rho}=1/2$ at the critical point—as predicted theoretically ($K_{\rho}<1/2$ is only expected at BOW-CDW critical points [36])—$\chi_{\mathrm{b}/\rho}(\pi)/L$ will converge to a constant at $\lambda=\lambda_{\mathrm{c}}$ but scale to zero (diverge) for $\lambda<\lambda_{\mathrm{c}}$ ($\lambda>\lambda_{\mathrm{c}}$). In contrast to the CDW-BOW transition considered in Sec. III.3, the presence of logarithmic corrections close to the critical point complicates a precise estimation of $\lambda_{\mathrm{c}}$. It was recently demonstrated for the $t$-$V$ model that $\chi_{\mathrm{b}/\rho}(\pi)/L$ at $\lambda=\lambda_{\mathrm{c}}$ increases slowly as a function of $L$ for the system sizes considered here [64]. Figures 7(b) and 7(d) are therefore in agreement with our estimates $\lambda_{\mathrm{c}}\approx 2.75$ from $C_{\mathrm{b}/\rho}(L/2)$. In principle, critical values can also be extracted from the correlation length, which shows a logarithmic scaling similar to the stiffness at $\lambda=\lambda_{c}$ [84].

The critical values from Fig. 7 are compatible with those from the stiffness fits. The assumption of $K_{\rho}=1/2$ at $\lambda_{\mathrm{c}}$ can give a more precise estimate of $\lambda_{\mathrm{c}}$ when using the susceptibility instead of the stiffness fits. However, the stiffness fits are more general and do not require $K_{\rho}=1/2$. The fact that all our estimates are consistent does not imply $K_{\rho}=1/2$ at $\lambda_{\mathrm{c}}$. On the other hand, we do not find evidence for $K_{\rho}<1/2$ at the LL-BOW and LL-CDW transitions, contrary to predictions from functional RG calculations for the LL-BOW transition in the SSH model [41].

Appendix D Stiffness fits

Standard BKT universality is predicted for the LL-BOW and LL-CDW transitions both in a general LL [51] and specifically for the frustrated $XXZ$ chain [36]. A detailed RG analysis [84] gives the finite-size scaling forms

[TABLE]

which provide the leading corrections to Eq. (6). However, in light of the observed nonuniversal jumps, functional RG predictions of $K_{\rho}<1/2$ at the LL-BOW transition [41], and $K_{\rho}<1/2$ at the BOW-CDW transition according to our data and theory [36], we determined the critical values in Fig. 1 using fits based on Eq. (6) with three parameters: $D_{\rho}(\infty)$, $g$, and $C$. In contrast, $g$ and $D_{\rho}(\infty)$ can be computed exactly for the classical 2D $XY$ model (see below), leaving only one free parameter. Specifically, for $\beta t\sim L_{y}=\infty$ (1D quantum chain at $T=0$), $g=1$ and $D_{\rho}=2/\pi$ ($D_{\rho}=t/2$) for the 2D $XY$ (1D $t$-$V$) model [60; 85]. 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 a fit involving $N$ data points, i.e., stiffness values $O_{n}$ for different system sizes $L_{n}\in\{L_{1},L_{2},\dots,L_{N}\}$ with corresponding statistical errors $\sigma_{n}$, the reduced $\chi^{2}$ is calculated from

[TABLE]

Here, the number of degrees of freedom $\nu$ is given by $N-M$, where $M$ is the number of fit parameters, and $C_{n}$ is the stiffness value predicted by the fit for system size $L_{n}$.

For the fits, we restricted the range of the jump to $0<D_{\rho}(\infty)<{2t}/{\pi}$, 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_{\text{min}}$ was imposed. Otherwise, the choice $g=0$ gives good fits throughout the LL phase and there would be no minimum of $\chi^{2}_{\nu}$ at the critical point. The exact value of $g_{\text{min}}$ does not significantly affect the results and was chosen as $0.25$. Finally, the allowed range of $C$ was $[0,\infty[$.

An important test case for the generalized, multi-parameter fit ansatz (6) 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 8(a)–(c) give a comparison of results based on Eqs. (6), (23), and (24). All three fit functions yield very similar and hence compatible minima of $\chi^{2}_{\nu}$ at the correct value $\lambda=2$. Figures 8(d)–(f) are based on fits that exploit the known values $g=1$ and $D_{\rho}(\infty)=t/2$. This additional information produces significantly sharper minima, in accordance with previous work on 2D $XY$ models [59]. At the same time, the first-order fit functions (6) and (23) do not fully capture the finite-size scaling on small system sizes, as manifested in $\chi^{2}_{\nu}\gg 1$ even at $\lambda=2$ in Figs. 8(d) and (e) for $L\geq 22$ and $L\geq 30$. Higher-order corrections are partially captured by varying $g$ and $D_{\rho}(\infty)$ [60], which explains the much better $\chi^{2}_{\nu}$ for the same range of $L$ in Figs. 8(a) and (b).

For the more challenging case of $t/\omega_{0}>0$, we focused on three-parameter fits based on Eqs. (6) and (23). Within the present accuracy, the results in Fig. 9 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 $L$ due to a reduced number of degrees of freedom. A similar picture arises for a fixed $\lambda=4$ in Figs. 9(e) and 9(f). For the present accuracy and range of system sizes, we cannot discriminate between the scaling forms (6), (23), and (24).

Appendix E 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}$ [60; 86]. For 1+1D quantum systems, $r=cL/\beta$, with $c$ a model-dependent constant. For example, $D_{\rho}(\infty)$ varies significantly as a function of $L_{x}/L_{y}$, covering the whole range from $2/\pi$ to [math] [86]. Similarly, $g$ varies between $1$ and $\infty$ as a function of $r$. In principle, a change of the range of the retarded interaction can mimic a change in $r$, leading to a dependence of $D_{\rho}(\infty)$ and $g$ on the phonon frequency.

There are several other known mechanisms for nonuniversal values of the stiffness jump. The bosonization relation $D_{\rho}=K_{\rho}u$ [51] implies that, even if $K_{\rho}=1/2$ at a QPT, $D_{\rho}(\infty)$ can change via the renormalized velocity $u$. For example, $u$ increases with $V$ in the $t$-$V$ model [42] but decreases with $\lambda$ in the Holstein model [80]. The stiffness can also be reduced by non-vortex excitations that are not captured by the standard BKT theory of the $XY$ model [87].