# Spin-Peierls transition in the Heisenberg chain with finite-frequency phonons

\[

## Abstract

We study the spin-Peierls transition in a Heisenberg spin chain coupled to optical bond-phonons. Quantum Monte Carlo results for systems with up to $`N=256`$ spins show unambiguously that the transition occurs only when the spin-phonon coupling $`\alpha `$ exceeds a critical value $`\alpha _c`$. Using sum rules, we show that the phonon spectral function has divergent (for $`N\mathrm{}`$) weight extending to zero frequency for $`\alpha <\alpha _c`$. The equal-time phonon-phonon correlations decay with distance $`r`$ as $`1/r`$. This behavior is characteristic for all $`0<\alpha <\alpha _c`$ and the $`q=\pi `$ phonon does not soften (to zero frequency) at the transition.

\]

The $`S=1/2`$ Heisenberg spin chain is unstable towards dimerization (the spin-Peierls transition) when coupled to an elastic lattice . For phonons in the adiabatic limit, this transition has been predicted to occur for arbitrarily weak spin-lattice coupling . On the other hand, recent work for optical phonons in the anti-adiabatic (high-frequency) limit suggests a transition only above a critical coupling . Furthermore, the mechanism of the transition in this limit was suggested to be qualitatively different, with no softening of the $`q=\pi `$ phonon . A way to reconcile the results in the adiabatic and anti-adiabatic limits has been proposed within an improved mean-field (RPA) theory , with the result that a complete phonon softening occurs only for bare phonon frequencies $`\omega _0`$ less than a critical value. For higher $`\omega _0`$, a central peak appears in the phonon spectral function and the phonon branch remains gapped. Considering the manifestly uncontrolled nature of mean-field calculations in one dimension, non-perturbative results in the regime of phonon frequencies comparable to the magnetic exchange energy $`J`$ are required to test this novel scenario.

In this Letter, we address the issues of a critical spin-phonon coupling and the mechanism of the zero temperature dimerization transition in the strictly one-dimensional case, using quantum Monte Carlo (QMC) simulations to obtain numerically exact results for relatively large systems. The model we study is defined by the Hamiltonian

$$H=J\underset{i=1}{\overset{N}{}}(1+\alpha x_i)𝐒_i𝐒_{i+1}+\omega _0\underset{i=1}{\overset{N}{}}n_i,$$
(1)   
where

$$x_i=(a_i^++a_i)/\sqrt{2}$$
(2)   
is the phonon coordinate, and $`n_i=a_i^+a_i`$ is the phonon occupation number at bond $`i`$. We use a recently developed QMC method based on sampling the perturbation expansion in the interaction representation . For a finite lattice at finite inverse temperature $`\beta `$, the expansion converges for any decomposition of $`H=H_0+V`$ into diagonal ($`H_0`$) and perturbing ($`V`$) terms and can be used  as a basis for a “worldline” Monte Carlo algorithm in continuous imaginary time (i.e., without invoking the Trotter decomposition ). In a slight modification of the scheme introduced in Ref. , we here include only the bare phonons in the diagonal term; $`H_0=\omega _0_in_i`$. The updating of the spin degrees of freedom can then be carried out using a new and highly efficient “operator-loop” algorithm , which in particular allows for sampling of all winding number sectors and hence direct evaluation of the spin stiffness .

In this work we consider only an energy $`\omega _0=J/4`$ for the bare phonons and study the behavior for values of the spin-phonon coupling in the range $`0\alpha /J0.5`$. We have studied systems with $`N`$ up to $`256`$ at inverse temperatures $`\beta =J/T`$ sufficiently high to give ground state results. Typically, for the system sizes we have considered, $`\beta `$ as high as $`2N`$ is required to achieve convergence to the $`T=0`$ limit of all the quantities of interest. We have used at least $`\beta =4N`$ for all calculations presented here.

Note that in the model, Eq. (1), for $`\alpha >0`$ there is an energy gain associated with an average uniform phonon displacement $`x=(1/N)_ix_i>0`$, which leads to an increased average effective spin-spin coupling $`J_{\mathrm{eff}}=J+\alpha x>J`$. For $`\omega _0/J=0.25`$ at $`T=0`$, we find $`J_{\mathrm{eff}}/J=1.018,1.071,1.158,1.278`$, and $`1.430`$ for $`\alpha =0.1,0.2,0.3,0.4`$, and $`0.5`$, respectively. We will in some cases measure energies in units of $`J_{\mathrm{eff}}`$ instead of the bare exchange $`J`$.

The most direct signal of the dimerization that can be measured in our simulations is the approach of the staggered phonon-phonon correlation function $`(1)^rx_ix_{i+r}`$ to a non-zero value at long distances $`r`$. Our results for this quantity indicate a critical coupling $`0.1<\alpha _\mathrm{c}/J<0.35`$ for $`\omega _0/J=0.25`$ . In order to improve on the accuracy of this rough estimate, and to circumvent potential problems with detecting a very small dimerization (as would be the case for very weak coupling if the mean-field result $`\alpha _\mathrm{c}=0`$ would be correct), we have also considered several other quantities. It is particularly useful to study the effects of the dynamic phonons in the spin sector (in general, our simulation results for spin quantities have smaller statistical fluctuations than the phonon correlations). Thus we discuss here results for the spin stiffness and the staggered spin susceptibility.

If the dimerization transition occurs at some critical coupling $`\alpha _\mathrm{c}>0`$, it is expected to be of the Kosterliz-Thouless (KT) type . The spin stiffness $`\rho _\mathrm{s}`$, i.e., the ground state energy curvature with respect to a uniform twist $`\varphi `$ in the spin-spin interaction ,

$$\rho _\mathrm{s}=\frac{1}{N}\frac{^2E_0(\varphi )}{\varphi ^2},$$
(3)   
is then expected to exhibit a discontinuous jump from a finite value for $`\alpha \alpha _\mathrm{c}`$ to zero for $`\alpha >\alpha _\mathrm{c}`$ (reflecting the opening of a spin gap). For a finite system the jump will be smoothed. In Figure 1 we show results for the stiffness versus $`\alpha /J`$ for several system sizes. The behavior expected for a KT transition is seen clearly — $`\rho _s`$ rapidly approaches zero for $`\alpha /J0.4`$ but appears to converge to a finite value for $`\alpha /J0.2`$, indicating a critical coupling between these values, in agreement with our previous results for the dimerization. It is, however, not easy to extract an accurate value for $`\alpha _\mathrm{c}`$ using these results. The scaling behavior is complicated by logarithmic corrections, which we expect to be present for all $`\alpha \alpha _c`$ as in the case of the Heisenberg chain \[i.e., $`\alpha =0`$ in Eq. (1)\]. This is in contrast to the finite temperature KT transition in the two-dimensional XY model, where $`\rho _s`$ approaches its asymptotic value algebraically for $`T<T_\mathrm{c}`$ and logarithmically only exactly at $`T_\mathrm{c}`$ . An indication of the difficulties associated with the log corrections in the spin-phonon chain can be seen in our stiffness data for $`\alpha =0`$, for which the exact infinite-size value is known to be (in our units)  $`\rho _s=1/4`$; about $`8\%`$ lower than what we find for $`N=128`$.

Although the log corrections complicate the extraction of $`\alpha _\mathrm{c}`$ from the stiffness data, their presence in other quantities can in fact be useful in numerical calculations. The asymptotic behavior of the spin-spin correlation function of the Heisenberg chain is known from bosonization and conformal field theory ;

$$S_i^zS_{i+r}^z\frac{(1)^r}{r}\mathrm{ln}^{1/2}(r/r_0).$$
(4)   
We expect this form to apply for all $`\alpha <\alpha _c`$. The logarithmic correction should vanish at the critical point $`\alpha _c`$, as it is known to do, e.g., at the critical point of the frustrated $`J_1J_2`$ chain . We have calculated the staggered spin susceptibility

$$\chi _s(\pi )=\frac{1}{N}\underset{m,n}{}(1)^{nm}_0^\beta 𝑑\tau S_n^z(\tau )S_m^z(0),$$
(5)   
for which Eq. (4) and conformal invariance imply the finite-size scaling form 

$$\chi _s(\pi )N\mathrm{ln}^{1/2}(N/N_0).$$
(6)   
In Figure 2 we graph $`(\chi _s(\pi )/N)^2`$ vs $`\mathrm{ln}(N)`$ for $`\alpha /J`$ in the range $`0.10.3`$. For $`\alpha /J=0.1`$ and $`0.2`$ the linear behavior for the larger system sizes is consistent with the form (6) expected in the gapless phase, whereas in the $`\alpha /J=0.3`$ case there is a clear decrease with increasing $`N`$, corresponding to a finite asymptotic value for $`\chi _s(\pi )`$ and therefore the presence of a spin gap. For $`\alpha /J=0.23`$ the curve is flat within statistical errors for $`N64`$, implying that $`\chi _s(\pi )`$ diverges linearly with $`N`$ without log correction. Based on this result (and calculations for other $`\alpha /J`$ close to $`0.23`$) we conclude that $`\alpha _c/J=0.23\pm 0.01`$.

Having established a KT transition and the critical coupling, we now turn to the question of the behavior of the $`q=\pi `$ phonons at the transition. We consider the phonon spectral function

$$A(q,\omega )=\underset{m,n}{}\mathrm{e}^{\beta E_n}|m|x_q|n|^2\delta (\omega [E_mE_n]),$$
(7)   
where

$$x_q=\frac{1}{\sqrt{N}}\underset{j=1}{\overset{N}{}}\mathrm{exp}(iqj)x_j.$$
(8)   
This real-frequency dynamic quantity cannot be obtained directly in our simulations. In order to avoid the problems associated with numerically continuing imaginary time data to real frequency, we here study sum rules that relate $`A(q,\omega )`$ to quantities that can be directly calculated. Two useful integrals that can be easily obtained from Eq. (7) are

$`S_x(q)`$ $`={\displaystyle _0^{\mathrm{}}}𝑑\omega A(q,\omega )(1+\mathrm{e}^{\beta \omega }),`$ (10)   
$`\chi _x(q)`$ $`=2{\displaystyle _0^{\mathrm{}}}𝑑\omega A(q,\omega )\omega ^1(1\mathrm{e}^{\beta \omega }),`$ (11)   
where $`S_x(q)`$ and $`\chi _x(q)`$ are the static structure factor and susceptibility;

$`S_x(q)`$ $`=x_qx_q`$ (13)   
$`\chi _x(q)`$ $`={\displaystyle _0^\beta }𝑑\tau x_q(\tau )x_q(0).`$ (14)   
Using Eqs. (10) and (11) one can readily verify  that the ratio

$$R(q)=2S_x(q)/\chi _x(q)$$
(15)   
is an upper bound for the lowest phononic excitation of momentum $`q`$. For $`\alpha \alpha _c`$ we therefore expect $`R(\pi )0`$ as $`N\mathrm{}`$, reflecting the presence of two degenerate ground states with momenta $`0`$ and $`\pi `$ (linear combinations of the two possible real-space dimerized states). For $`\alpha =0`$, $`R(q)=\omega _0`$ for all $`q`$. A transition caused by a softening of the $`q=\pi `$ phonon would imply $`R(\pi )>0`$ for $`\alpha <\alpha _c`$ and $`R(\pi )0`$ as $`\alpha \alpha _c`$. This behavior is not seen in our results. Instead, $`R(\pi )`$ appears to approach zero as $`N`$ is increased even for $`\alpha `$ much smaller than the critical value, as shown in Figure 3. Hence the spectral weight extends to zero frequency also in the non-dimerized systems.

We also find that the total spectral weight, given by the structure factor according to Eq. (10), diverges with $`N`$. In Figure 4 we graph $`S_x(\pi )`$ versus the logarithm of the system size for values of $`\alpha `$ both below and above the critical coupling. We find a linear increase for $`\alpha <\alpha _c`$, indicating a logarithmic divergence and therefore an inverse distance decay of the real-space phonon-phonon correlation function. For $`\alpha =0.3>\alpha _c`$ the divergence is faster than logarithmic. The expected linear in $`N`$ behavior cannot be observed close to $`\alpha _c`$ for the system sizes we have studied, due to large short-distance contributions to $`S_x(\pi )`$. For $`\alpha 0.4`$ we do observe an almost linear divergence with $`N`$.

These results show that, in the thermodynamic limit, there is infinite $`q=\pi `$ phonon spectral weight also for $`\alpha <\alpha _c`$. This weight extends to zero frequency. The rate of decay of $`R(\pi )`$ with increasing $`N`$, seen in Fig. 3, shows that the low-frequency weight grows rapidly with $`N`$. The only plausible explanation for this is that the phonon spectral function has a central peak with infinite integral. We now elaborate on the reasons for this behavior.

In the absence of spin-phonon couplings, the low-lying excitations of the system are the two-spinon singlet and triplet states of the Heisenberg chain. Our results indicate that an arbitrarily weak coupling to the phonons induces an infinite phonon spectral weight into these states at the staggered momentum $`q=\pi `$. The asymptotic real-space staggered phonon correlation function has the same $`1/r`$ decay as the spin-spin correlation function (perhaps differing by multiplicative logarithmic corrections that cannot be detected in our results for the phonon correlations). This is not completely surprising, considering that the Heisenberg chain is also characterized by an inverse distance decay of the dimerization correlation function $`(S_iS_{i+1})(S_{i+r}S_{i+1+r})`$, as also recently noted by Gros and Werner . The corresponding susceptibility is therefore divergent and this leads to the spontaneous dimerization for arbitrarily weak spin-phonon couplings in the adiabatic case $`\omega _0=0`$. What we have shown here is that dynamic ($`\omega _0>0`$) phonons destroy the long-range order for weak spin-phonon couplings but nevertheless the spin and phonon excitations remain coupled in a manifestly non-perturbative fashion.

The most plausible scenario for the mechanism of the spin-Peierls transition is then the following: For weak coupling $`\alpha `$, the phonon spectral function at $`q=\pi `$ has a finite-weight peak close to the bare frequency $`\omega _0`$, as well as a central peak with divergent frequency integral \[$`\mathrm{ln}(N)`$ divergent as a function of system size\]. As $`\alpha `$ is increased the finite-frequency peak may shift slightly but remains at finite frequency. The central peak sharpens and at $`\alpha =\alpha _c`$ acquires a $`\delta `$ function component, corresponding to the development of static long-range order. For $`\alpha >\alpha _c`$ this $`\delta `$ function remains, and a gap develops to the remainder of what was the finite-width central peak for $`\alpha \alpha _c`$. This gap is the excitation energy of a lattice/magnetic soliton pair.

On general grounds, we find it unlikely that there would be any qualitative changes in the nature of the $`T=0`$ transition as $`\omega _0`$ is varied. A recent density matrix renormalization group study of a model closely related to the one we have studied also finds unambiguously that the critical coupling $`\alpha _c>0`$ for any $`\omega _0`$ . In addition, our results are consistent with the calculations in the anti-adiabatic limit  (although the infinite low-frequency $`q=\pi `$ phonon weight for $`\alpha <\alpha _c`$ was not noted there), even though our $`\omega _0=J/4`$ is closer to the adiabatic regime.

We have here discussed only the $`T=0`$ quantum phase transition in the strictly one-dimensional case. The finite $`T_c`$ in real materials such as CuGeO<sub>3</sub>  can be due to three-dimensional phonons, as well as interchain magnetic couplings. The non-softening nature of the quantum phase transition that we have found here clearly supports the suggestion  that the finite-$`T`$ transition may also be non-softening. In the improved RPA theory , softening occurs below a critical value of $`\omega _0`$. In future work, we plan to extend our simulations to two- and three-dimensional systems and hope to address this important issue.

We would like to thank Robert Bursill, Valeri Kotov, Ross McKenzie, Oleg Sushkov, and Johannes Voit for useful discussions. The numerical simulations were carried out at the NCSA facilities at the University of Illisois at Urbana-Champaign. This work is supported by the National Science foundation under Grant No. DMR-97-12765.