format: text

text:
# Peierls transition in the presence of finite-frequency phonons in the one-dimensional extended Peierls-Hubbard model at half-filling

## I Introduction

Almost half a century ago, Peierls demonstrated that a one-dimensional (1D) metal coupled to an elastic lattice could exhibit an instability towards a lattice distortion of wave vector $`q=2k_F`$.peierls  This leads to a gap in the electronic spectrum at the Fermi energy and, for the case of a half-filled band, the ground state is dimerized. Experimentally, the Peierls instability can be observed in a wide range of quasi-1D materials, e.g., conjugated polymers,pa  organic charge transfer salts,cts  MX salts,mx  and $`\text{CuGeO}_3`$.cugeo3  To explain quantitatively the properties of these materials, several different models extending beyond the original non-interacting Peierls model with a classical lattice have been proposed.review  These include the Su-Schrieffer-Heeger (SSH) model,ssh ; nakahara ; su ; fradkin-ssh ; zheng  the Holstein model,holstein ; fradkin-hol ; jeckelmann-hol  various Peierls-Hubbard phmodels ; takahashi ; fehske  and extended Peierls-Hubbard hirsch ; zimanyi ; takahashi ; jeckelmann-mi ; barford  models, as well as spin-Peierls models. pincus ; pytte ; cross ; sf ; uhrig ; gros ; wellein ; anders-sp ; robert-sp ; weisse ; raas

The Peierls instability is well understood in the static lattice limit (adiabatic phonons) in the absence of electron-electron (e-e) interactions. The ground state is a Peierls state for arbitrarily small electron-phonon (e-ph) coupling. However, the quantum lattice and e-e interaction effects are still not completely understood. For many quasi-1D materials the zero-point fluctuations of the phonon field are comparable to the amplitude of the Peierls distortion,su  and this has spurred a large number of studies of quantum fluctuations in SSH,nakahara ; fradkin-ssh ; zheng  and Holstein, fradkin-hol ; jeckelmann-hol  models. Recent numerical studies have shown that quantum fluctuations destroy the Peierls instability for small e-ph coupling and/or large phonon frequency in both the spinless sf  and spin-$`\frac{1}{2}`$ Holstein models jeckelmann-hol  at half-filling. For the SSH model, Fradkin and Hirschfradkin-ssh  carried out an extensive study of spin-$`\frac{1}{2}`$ ($`n=2`$) and spinless ($`n=1`$) fermions. In the anti-adiabatic limit (vanishing ionic mass), they mapped the system to an $`n`$-component Gross-Neveu model, which is known to have long-ranged dimerization for arbitrary coupling for $`n2`$ but not for $`n=1`$. For non-zero ionic mass, a renormalization group analysis showed that the low-energy behavior of the $`n=2`$ model is still governed by the zero mass limit of the theory. This implies that the spin-$`\frac{1}{2}`$ SSH model (but not the spinless model) has a dimerized ground state for arbitrarily weak e-ph coupling. Early numerical calculations are also consistent with this scenario.fradkin-ssh  It should be noted, however, that a dimerized state was also predicted for any non-zero e-ph coupling in the spin-$`\frac{1}{2}`$ Holstein model,fradkin-hol  for which more recent large-scale calculations have instead indicated a non-zero critical coupling.jeckelmann-hol

Independent-electron models, such as SSH and Holstein, are important from a theoretical standpoint but are not sufficient to account quantitatively for the experimentally observed properties of real materials. For that e-e interactions have to be included in addition to the e-ph couplings. pa  The interplay among the different interactions gives rise to a rich variety of broken-symmetry ground states as well as low-energy electron-lattice excitations like solitons, polarons, bipolarons, etc. pa  At half-filling, on-site (Hubbard) interactions open a charge gap,liebwu  and in the absence of e-ph couplings the system is then a Mott insulator with algebraically decaying ($`1/r`$ as a function of distance $`r`$) spin-spin correlations. Hence, the Peierls transition in this case is accompanied only by the opening of a spin gap, the charge gap already generated by the e-e interactions. Longer-range e-e interactions can destroy the Mott state, however, and hence must affect also the Mott-Peierls transition. Without phonons, even in the simplest half-filled extended Hubbard model with only on-site ($`U`$) and nearest-neighbor ($`V`$) interactions, solyom ; voit  some features of the phase diagram are still controversial. nakamura ; pinaki ; jeckelman-ehm ; jcomment  Adding e-ph interactions further increases the complexity of the problem, and the determination of the phase diagram remains a very challenging problem.

Zimanyi et al.zimanyi  have investigated models with both e-e and e-ph interactions using “g-ology” and RG techniques. They showed that the ground state has a spin gap if the combined backscattering amplitude $`g_1^T=g_1(\omega _0)+\stackrel{~}{g}_1(\omega _0)`$ is negative, where $`g_1(\omega _0)`$ is the contribution from e-e interactions and $`\stackrel{~}{g}_1(\omega _0)<0`$ comes from the e-ph interactions. Thus, in the extended Hubbard model, if the bare coupling $`g_1=U2V`$ is positive, then $`g_1^T0`$ for small values of the e-ph coupling and there is no spin gap. The transition to a Peierls state occurs only when the e-ph coupling exceeds a critical value and $`g_1^T`$ becomes negative. It should be noted, however, that the conventional scenario for the behavior close to the line $`U=2V`$ has recently been challenged nakamura ; pinaki ; jcomment  and a Peierls-like bond-ordered state most likely appears close to $`V=U/2`$ even in the absence of e-ph couplings. For the pure SSH model ($`U=V=0`$), $`g_1`$ is zero and $`g_1^T`$ is negative for any non-zero e-ph coupling. This implies a Peierls ground state for arbitrarily small e-ph coupling, in agreement with the earlier results of Fradkin and Hirsch.fradkin-ssh

In the limit of strong on-site e-e interactions, which inhibit doubly occupied sites, a half-filled system can be mapped to a spin-phonon model, which also can undergo a dimerization (spin-Peierls) transition. pincus ; cross  Extensive studies of spin-Peierls models in recent years uhrig ; gros ; wellein ; anders-sp ; robert-sp ; weisse ; raas  have largely been spurred by the discovery of a spin-Peierls transition at unusually high temperature (14 K) in $`\text{CuGeO}_3`$.cugeo3  Several different calculations, for different types of spin-phonon couplings, have shown that the transition occurs only above a finite spin-phonon coupling in the presence of finite-frequency phonons.uhrig ; anders-sp ; robert-sp ; weisse ; raas  This is in contrast to the adiabatic limit, where dimerization occurs for infinitesimal coupling.cross

Numerical studies of models with both e-ph and e-e interactions, in which the charge degrees of freedom are retained (Peierls-Hubbard and extended Peierls Hubbard models), have addressed the effect of interactions on the dimerization amplitude hirsch ; takahashi  and the excited states. barford  Detailed studies of the phase diagrams have in the past been limited by the small lattice sizes accessible when both e-e and e-ph interactions are included.fehske  The situation is rapidly improving, however, as modern quantum Monte Carlo sseloop ; anders-sp ; pinaki  and density matrix renormalization group robert-sp ; barford  methods can now access models with both e-e and e-ph interactions on chains with several hundred sites.

Here we consider a 1D extended Hubbard model with on-site ($`U`$) and nearest-neighbor ($`V`$) e-e interactions and couple it to dispersionless optical bond phonons via modulation of the electron kinetic energy. We study the transition from a Mott insulating state, with dominant spin-spin correlations, to a Peierls (dimerized) spin-gapped state. Since the parameter space of this model is rather large, with a bare phonon frequency ($`\omega _0`$) and an e-ph coupling ($`\alpha `$) in addition to the e-e couplings, we have limited our study to a physically reasonable ratio $`U/V=4`$ of the e-e parameters.

In Sec. II we define the model and the various physical quantities that we have calculated using a quantum Monte Carlo method (stochastic series expansion sseloop ). In Sec. III we discuss several signals that we have used to detect the Mott-Peierls transition. In Sec. IV we present the phase diagram in the ($`\omega _0,\alpha `$)-plane for a fixed value of the on-site interaction $`U`$ that has previously been used in models of polyacetylene. We also discuss the effects of varying the e-e interaction strength at fixed $`\omega _0`$. We map the model to a spin-Peierls model for large $`U,V`$ and compare the Mott-Peierls boundary with known spin-Peierls results. In Sec. V we summarize our results and discuss some future prospects.

## II Model and observables

The Hamiltonian is given by

$`H`$ $`=`$ $`t{\displaystyle \underset{i,\sigma }{}}(1+\alpha [a_i^{}+a_i])(c_{i+1,\sigma }^{}c_{i,\sigma }+c_{i,\sigma }^{}c_{i+1,\sigma })`$ (1)
$`+\mu {\displaystyle \underset{i}{}}n_i+U{\displaystyle \underset{i}{}}(n_{i,}{\displaystyle \frac{1}{2}})(n_{i,}{\displaystyle \frac{1}{2}})`$
$`+V{\displaystyle \underset{i}{}}(n_i1)(n_{i+1}1)+\omega _0{\displaystyle \underset{i}{}}a_i^{}a_i,`$
where $`a_i^{}(a_i)`$ creates(annihilates) a phonon on the bond between sites $`i`$ and $`i+1`$, $`c_{i,\sigma }^{}(c_{i,\sigma })`$ is the spin-$`\sigma `$ electron creation(annihilation) operator, and $`n_i=n_{i,}+n_{i,}`$. For the half-filled band that we study here, the chemical potential $`\mu =0`$. We set the single-electron hopping $`t`$ to unity. For the e-e interactions, we first take the values $`U=2.5`$, $`V=U/4`$, which have previously been used as values corresponding approximately to what is expected in polyacetylene.jeckelmann-mi  We will also consider other values of $`U`$, keeping the ratio fixed at $`U/V=4`$. We study the system as a function of the bare phonon frequency $`\omega _0`$ and the e-ph interaction $`\alpha `$.

The dispersionless optical phonons we use are different from the bare SSH phonons, which have vanishing energy for momentum $`q0`$. However, these acoustic phonons decouple from the electronic low-energy states involved in the Peierls instability and therefore only the optical phonons close to $`q=\pi `$ need to be kept.fradkin-ssh ; zimanyi  Hence, in this regard we expect the optical phonons in (1) to be equivalent to fully quantum mechanical SSH phonons. In the non-interacting limit ($`U,V0`$), the ground state should therefore be a dimerized Peierls insulator for any non-zero $`\alpha `$.fradkin-ssh

To obtain numerically exact ground state results we have used the stochastic series expansion (SSE) quantum Monte Carlo method sseloop  for periodic chains with up to $`N=256`$ sites. The SSE method is a finite-temperature technique based on importance sampling of the diagonal elements of the Taylor expansion of $`e^{\beta H}`$, where $`\beta `$ is the inverse temperature; $`\beta =t/T`$. Ground state expectation values can be obtained using sufficiently large $`\beta `$, and there are then no approximations beyond the statistical errors. Typically, $`\beta =2N`$ or $`4N`$ was sufficient for the quantities presented here to have converged to their ground state values. Using the recently developed “operator loop” update, sseloop  the electronic degrees of freedom are treated in the same manner as in the absence of phonons.pinaki  The phonons are also treated in the occupation number basis directly with the SSE representation pinakithesis  (i.e., slightly different from the interaction picture used for the phonons in in Refs. anders-sp,  and irsse, ). At the (low) energy scales that we are interested in here, the number of phonons per bond is small (typically $`<10`$) and there are no problems in using a truncated basis (the truncation can be arbitrarily large in the SSE).

The Mott and Peierls phases can be characterized using the static spin (S) and bond (B) structure factors and susceptibilities, as will be further discussed in Sec. III. The structure factors are defined by

$`S_S(q)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{k,l}{}}e^{iq(kl)}S_k^zS_l^z,`$ (2)
$`S_B(q)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{k,l}{}}e^{iq(kl)}K_kK_l,`$ (3)
where $`K_j=_{\sigma =,}(c_{j+1,\sigma }^{}c_{j,\sigma }+\mathrm{h}.\mathrm{c}.)`$ is the kinetic energy operator on the $`i`$th bond. The corresponding static susceptibilities are given by

$`\chi _S(q)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{k,l}{}}e^{iq(kl)}{\displaystyle _0^\beta }𝑑\tau S_k^z(\tau )S_l^z(0),`$ (4)
$`\chi _B(q)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{k,l}{}}e^{iq(kl)}{\displaystyle _0^\beta }𝑑\tau K_k(\tau )K_l(0).`$ (5)
Direct evidence for the presence or absence of the spin and charge gaps can also be obtained from spin and charge stiffnesses $`\rho _c`$ and $`\rho _s`$, which are defined as the second derivative of the internal energy per site with respect to phase factors multiplying the kinetic energy; $`\rho _{c,s}=^2E(\varphi _{c,s})/\varphi _{c,s}^2`$.kohn  The SSE estimators for all these observables can be found in Ref. pinaki, .

## III Detecting the Mott-Peierls transition

For our choice of $`V=U/4`$, the ground state in the limit of zero e-ph interaction ($`\alpha =0`$) is a Mott insulator with no spin gap (but finite charge gap) and is characterized by a $`1/r`$ decay of the staggered spin-spin correlations.voit ; pinaki  The transition to a dimerized Peierls state is marked by the development of a staggered kinetic energy modulation and will hence be signalled by divergent peaks at $`q=\pi `$ in the bond-order structure factor, Eq. (3), and susceptibility, Eq. (5). The dimerization is accompanied by the opening of a spin gap, the charge gap remaining finite. Hence the $`q=\pi `$ peak in the spin structure factor and susceptibility, Eqs. (2) and (4), which diverge in the Mott phase, become non-divergent in the Peierls state.

In the adiabatic limit, the system is dimerized for any $`\alpha >0`$. We here present several results showing that the Peierls transition occurs at a critical coupling $`\alpha _c>0`$ when $`\omega _0=1`$ and $`U=2.5`$. The phase diagrams discussed in Sec. IV are based on the same signals for the transition at other e-e couplings and $`\omega _0`$.

Since the charge gap is finite in both the Mott and Peierls states, the charge stiffness $`\rho _c`$ should vanish in the thermodynamic limit for all $`\alpha `$. The upper panel of Fig. 1 shows $`\rho _c`$ as a function of $`\alpha `$ for several system sizes $`N`$. As $`N`$ grows, $`\rho _c`$ indeed rapidly converges to zero for all $`\alpha `$, in agreement to our expectations. The Mott state has no spin gap (finite spin stiffness) whereas the Peierls state has a finite spin gap (zero spin stiffness). If the Peierls transition occurs at a critical coupling $`\alpha _c>0`$, it should be of the Kosterliz-Thouless type, voit  where the spin stiffness changes discontinuously from a finite value for $`\alpha \alpha _c`$ to zero for $`\alpha >\alpha _c`$. The spin stiffness graphed in the lower panel of Fig. 1 shows a jump developing with increasing $`N`$, indicating a critical coupling $`\alpha _c0.3`$.

The spin stiffness data do not easily yield a more accurate estimate of the critical coupling. As discussed in the context of the spin-Peierls model, anders-sp  logarithmic corrections lead to large finite-size effects for $`\alpha \alpha _c`$. A more accurate estimate can be obtained from the scaling behavior of the finite-size staggered bond and spin susceptibilities.anders-sp  It is known from bosonization studies that in the Mott phase the equal-time staggered spin and bond correlations both decay with distance $`r`$ as $`1/r`$, up to multiplicative logarithmic corrections.voit  At the Mott-Peierls phase boundary, the log-corrections can be expected to disappear,eggert1  and this can be used as a criterion for the phase transition. In the dimerized Peierls phase the bond correlation function approaches a constant at long distances, whereas the spin correlations decay exponentially. It is convenient to study the associated static susceptibilities defined in Eqs. (4) and (5), which in a critical state scale with one power of $`N`$ higher than the structure factors Eqs. (2) and (3). In the Peierls phase $`\chi _S(\pi )/N`$ should converge to $`0`$ and $`\chi _B(\pi )/N`$ should diverge, whereas in the Mott phase $`\chi _S(\pi )/N`$ should diverge logarithmically and $`\chi _B(\pi )/N`$ should approach zero logarithmically (the log-corrections for spin and bond correlations are differentvoit ). Fig. 2 shows both quantities versus $`\text{ln}(N)`$ for several values of $`\alpha `$. The expected behavior is indeed observed, and within statistical errors both $`\chi _S(\pi )/N`$ and $`\chi _B(\pi )/N`$ are independent of $`N`$ for the largest chains when $`\alpha 0.26`$.

Additional confirmation of the critical coupling is obtained by studying the behavior of $`S_S(q)/q`$ as $`q0`$. It has been shown voit ; torsten  that if the ground state is gapless in the spin sector, then $`S_S(q)/q1/\pi `$ as $`q0`$, whereas if there is a spin gap, $`S_S(q)/q0`$. Even a very small spin gap can be detected this way, since it is in practice sufficient to see that $`\pi S_S(q)/q`$ decays below $`1`$ for small $`q`$ to conclude that a spin gap must be present. In Fig. 3 we present results for different values of $`\alpha `$. The curves for $`\alpha 0.26`$ are above 1 for all $`q`$, and the decay towards 1 is very slow. The asymptotic approach to 1 can be expected to be logarithmic.eggert2  On the other hand, the $`\alpha 0.30`$ curves drop below 1. From these results we estimate $`\alpha _c`$=0.28$`\pm `$0.01, which is compatible with the $`q=\pi `$ quantities in Fig. 2. In general, we have found that $`S_S(q)/q`$, which indirectly signals the opening of a spin gap, is the easiest and most reliable way to detect a bond-ordered state (see also Ref. pinaki, ).

## IV phase diagrams

The critical coupling depends on the parameters of the Hamiltonian, in particular, the bare phonon frequency $`\omega _0`$. Using the above criteria for distinguishing the Mott and Peierls phases, we have calculated $`\alpha _c`$ as a function of $`\omega _0`$, keeping $`U=2.5`$, $`V=0.625`$. As shown in Fig. 4, $`\alpha _c`$ decreases linearly to zero for small $`\omega _0`$. This phase diagram is hence consistent with the known $`\alpha _c=0`$ for adiabatic phonons.

For polyacetylene Fradkin and Hirsch fradkin-ssh  used rescaled phonon parameters, which in our units correspond to $`\omega _0=0.067`$, $`\alpha =0.052`$. This point is indicated in Fig. 4, and, in accordance with the strong dimerization of polyacetylene, is well within the Peierls phase.

As argued above, the e-ph coupling in the present model is similar to that in the SSH model for the purpose of studying the dimerization transition. This implies that in the limit of $`(U,V)0`$, we should be able to reproduce previous SSH results. In particular, according to Fradkin and Hirsch,fradkin-ssh  $`\alpha _c`$ should be zero even for finite frequency phonons. To verify this, we have studied $`\alpha _c`$ as a function of $`(U,V)`$, keeping a fixed ratio of $`U/V=4`$. With this ratio the ground state for all $`U`$ is a Mott insulator with zero spin gap in the limit of vanishing e-ph interaction.pinaki  We have studied only a single phonon frequency, $`\omega _0=1`$. The resulting phase diagram is presented in Fig. 5. The critical e-ph coupling decreases monotonically with decreasing $`(U,V)`$, but the smallness of the spin gap as $`(U,V)0`$ makes it hard to obtain reliable results below $`U=0.4`$. We can therefore not make a definite statement about this limit. Nevertheless, our results are consistent with a power-law behavior $`\alpha _cU^\gamma `$ with $`\gamma 0.3`$, but a logarithmic form for $`U0`$ can also not be excluded.

For large $`U,V`$, and with the ratio $`V=U/4`$, the extended Hubbard model can be mapped onto the spin-$`\frac{1}{2}`$ Heisenberg chain with exchange coupling $`J=4t^2/(UV)`$ (this mapping becomes invalid for $`VU/2`$, where phase transitions to bond-ordered and charge-ordered phases occur voit ; nakamura ; pinaki ). Carrying out this transformation for the full electron-phonon model (1), the phonon-modulated exchange is

$$J(x_i)=\frac{4t^2(1+\alpha x_i)^2}{UV},$$
(6)
where $`x_i=a_i^{}+a_i`$. Under the assumption that the non-linear term $`(\alpha x_i^2)`$ can be neglected, which is not a priori clear when $`\alpha `$ is not small, we obtain exactly the spin-Peierls model considered in, e.g., Refs. anders-sp,  and raas, :

$$H_{\mathrm{SP}}=\underset{i}{}(J+gx_i)𝐒_i𝐒_{i+1}+\omega _0\underset{i}{}a_i^{}a_i,$$
(7)
with

$$J=\frac{4t^2}{UV},g=\frac{8\alpha t^2}{UV}.$$
(8)
For the model (7), an analytic expression for the critical spin-phonon coupling $`g`$ has been obtained for the whole range of bare phonon frequencies $`\omega _0/J`$.raas  The form, Eq. (12) of Ref. raas, , is expected to be exact in the anti-adiabatic limit, $`\omega _0/J\mathrm{}`$, which here corresponds to $`U\mathrm{}`$. It is in good agreement with numerical (SSE) results anders-sp  for the spin-Peierls transition even for frequencies as low as $`\omega _0/J=0.25`$. In Fig. 5 we compare our SSE results for the extended Peierls-Hubbard model with the spin-Peierls form for $`U`$ up to $`10`$. The transition curve crosses the Mott-Peierls transition curve at $`U8`$, and is not in good agreement away from this point. The spin-Peierls critical $`\alpha `$ increases linearly with $`U`$ as $`U\mathrm{}`$, whereas the Mott-Peierls boundary has a slower increase with $`U`$. The poor agreement for small $`U`$ is not surprising, as the mapping to the spin chain, with coupling $`J`$ given by (8), can only be expected to be good for large $`U`$, and the form used for the spin-Peierls transition curve is not expected to be quantitatively accurate for very small $`\omega /J`$ (corresponding here to small $`U`$).raas  Considering, however, that the analytical form is accurate for the spin-Peierls model with $`\omega _0/J=0.25`$, raas ; anders-sp  and that the nonlinear spin-phonon coupling should be negligible for small $`\alpha `$, the disagreement for the small-$`U`$ region in Fig. 5 must be due to the poor correspondence between the full electron model and the spin chain with the lowest-order $`J`$ in Eq. (8). The poor agreement for $`U>8`$ indicates that the nonlinear spin-Phonon term $`(\alpha x_i)^2`$ does becomes important as $`U\mathrm{}`$.

One effect of the nonlinear coupling term is to renormalize $`J`$: Writing $`x_i^2=x_i^2+\mathrm{\Delta }(x_i^2)`$, the renormalized J is given by

$$J_{\mathrm{eff}}=\frac{4t^2(1+\alpha ^2x_i^2)}{UV},$$
(9)
i.e., $`J_{\mathrm{eff}}>J`$. Evaluation the spin-Peierls transition curve using $`J_{\mathrm{eff}}`$ instead of $`J`$ clearly would reduce $`\alpha _c`$ for given $`U`$ and bring the result closer to the actual Mott-Peierls curve in Fig. 5. However, we have not evaluated $`x_i^2`$ for a quantitative test of this effect. In any case, if $`x_i^2`$ is large there is also no reason to expect that the remaining nonlinear coupling $`\mathrm{\Delta }(x_i^2)`$ can be neglected if $`\alpha `$ is not small. This issue clearly deserves further study.

## V summary

In summary, we have studied several aspects of the phase diagram of an extended 1D Hubbard model coupled to optical bond phonons. We have demonstrated that the stochastic series expansion technique sseloop ; pinaki  can be used for large electron-phonon chains (here up to 256 sites) to compute several different quantities that signal the opening of the spin gap at the Mott-Peierls transition. The spin gap boundary is also in good agreement with direct probes of the bond order (i.e., kinetic-energy correlation functions and susceptibilities at $`q=\pi `$).

Our phase diagrams are in agreement with what is generally expected, but to our knowledge they have not been computed quantitatively before. For large e-e couplings, we have pointed out the relevance of an effective nonlinear spin-phonon coupling in the mapping of the Hubbard model to a spin chain. Because of this, standard spin-phonon models, where the nonlinear term is not included, cannot be expected to reproduce fully the phase diagrams of electron-phonon models.

The methods that we have used here should also be applicable to systems away from half-filling. We plan such calculations aimed at studying the stability of the soliton lattice, which is formed in doped systems in the adiabatic limit,ssh ; salkola ; jeckelmann-mi  in the presence of finite-frequency phonons. A previous quantum Monte Carlo study has indicated that the soliton lattice is stable.takahashi  However, open boundary conditions were used for relatively small chains, and therefore the issue of whether this is a stable phase on an infinite lattice remains to be clarified.

###### Acknowledgements.

This work was supported by the NSF under Grants No. DMR-99-86948 (PS) and DMR-97-12765 (DKC), and by the Academy of Finland, project 26175 (AWS). The numerical calculations were carried out at the NCSA in Urbana, Illinois.


text:
# Baryons from the lattice: past, present and future

## 1 Introduction

With new exciting experiments ongoing and planned at several facilities, distributed over three continents, the past few years have witnessed the emergence of a new arena of research on the borderline between particle and nuclear physics: hadron physics. These developments have been complemented by several break-throughs in our theoretical understanding of the long distance, non-perturbative regime of hadronic matter, directly from QCD. Here I will focus on the most recent lattice results.

Progress in the field has been slow but steady. Several “revolutions” in the lattice technology, for instance non-perturbative renormalisation and improvement programmes of action and operators, were often hardly identifiable for the outside observer, a feature that we share with many other fields of research. In the case of strong QCD the failure of immediate satisfaction is due to the highly non-trivial phenomenology emerging from a Lagrangian as simple as that of QCD and yet so colourful and complex that even now it keeps experimentalists and theoreticians very busy.

A major but all to often neglected application of lattice results is to establish validity range and applicability of phenomenological models or effective field theories directly from QCD and to calculate low energy parameters that appear in these models or approximations. In quite a few situations, lattice predictions have matured to a degree that allows direct confrontation with experiment. Triggered by new experimental programmes and improved interaction with other theoreticians, an increasing number of lattice practitioners became interested in what has now become the domain of hadronic physics. This more focused approach resulted in an accelerated development of the field over the past two years with more to come in the future. I will go through the present state of lattice studies, discuss the two major sources of error, and review results on spectrum and structure of baryons, in particular $`N\gamma \mathrm{\Delta }`$ transition form factors, charge distributions, three-body forces and the hyperon spin content.

## 2 The state of the lattice

Simulation results are obtained at finite lattice spacings and have to be extrapolated to the continuum limit. Moreover, the space-time box size should be big enough to comfortably accommodate the hadronic scales of interest. The dominant systematic uncertainties of present-day lattice results however tend to be related to unrealistically heavy $`u`$ and $`d`$ quark masses, which typically have to be chirally extrapolated over a long distance from $`m_\pi ^2>0.25m_\rho ^2`$ to the physical limit $`m_\pi ^20.03m_\rho ^2`$, often combined with the quenched approximation, i.e. neglecting the polarisation effects of sea quarks on the QCD vacuum.

The latter approximation, which easily reduces the computational effort by a factor of $`10^3`$, is justified in the limit of the number of colours $`N_c\mathrm{}`$, however, it is not a priori clear whether $`1/31`$. Phenomenology suggests that this assumption is not as outrageous as it might sound at first. This seems substantiated by simulations of pure $`SU(N_c)`$ gauge theories for $`N_c=2,\mathrm{},6`$ as well as by a comparison of quenched lattice predictions of the light hadron spectrum with experiment: ratios of light hadron masses on the lattice indeed have at last been found to be inconsistent with the observed spectrum. The differences are typically smaller than 10 %.

However, there are cases in which “quenching” clearly matters. Quantum mechanical perturbation theory tells us that an $`𝒪(ϵ^2)`$ correction to an energy eigenvalue corresponds to an $`𝒪(ϵ)`$ correction on wave functions which hence can easily be changed by as much as 30 %. This will affect quantities like decay constants and parton distributions. The situation with for instance bottomonium fine structure splittings or electroweak decay rates is even worse, as these are roughly proportional to the square of the wave function at the origin. Strong decays do not exist in this approximation as quarks and anti-quarks cannot annihilate, and therefore one might hesitate to expect very broad resonances to be modelled in a reasonable way. The same also holds true for any valence quark based model. We shall discuss some other shortcomings of quenching below.

Fortunately, with several new supercomputer installations on the horizon by early 2004 we are at the brink of the fourth generation of simulations including sea quarks. These certainly will bring us down in the quark mass to about the point where the $`\rho `$ meson becomes unstable and will at the same time allow for controlled continuum limit extrapolations.

Recently a way has been found to implement the analogue of exact chiral symmetry without “doublers” on the lattice, known as overlap, Neuberger or Ginsparg-Wilson (GW) fermions. One particular (truncated) realisation of the Neuberger operator are domain wall fermions. GW fermions are the only known consistent way of regularising a quantum field theory with massless fermions (like e.g. the standard model) and hence obviously a major theoretical break-through in particle physics as a whole. The famous Nielsen-Ninomiya no-go theorem is circumvented by relaxing the definition of chiral symmetry. I illustrate the yearly number of citations of these two important articles published in 1981 and 1982, respectively, in Figure 1 and leave judgement to the reader which paper would have better impacted on Research Assessment Exercises. In practical terms GW fermions mean that realistically light quark masses will become accessible to future numerical simulations, first quenched and subsequently un-quenched.

Going to light quark masses is expensive as the linear spatial lattice extent $`L`$ has to be large relative to the inverse pion mass, to avoid strong finite size effects from pions exchanged “around” the boundaries. For instance if we intend to simulate a physical pion with say $`L>4/m_\pi `$ this implies $`L>5.7`$ fm. The main limitation however is that the effort of inverting the lattice Dirac operator scales with a large power of the inverse quark mass. This will improve with GW quarks that do not suffer from un-physical zero or near-zero modes: while at $`m_\pi >0.5m_\rho `$ GW fermions are easily two orders of magnitude more expensive to simulate than “standard” formulations, there will be a reward at small quark masses. In the quenched approximation this improvement factor has already turned out to be infinite, rendering a pion lighter than 200 MeV from impossible to possible. Chiral fermions also greatly simplify operator mixing and renormalisation.

## 3 Chiral extrapolations and “quenching”

Based on the assumptions that QCD bound states are mesons and baryons, that there is a mass gap and spontaneous chiral symmetry breaking at zero quark mass, an effective low energy chiral Lagrangian can be derived in the spirit of the Born-Oppenheimer approximation. This Lagrangian will, to leading order, describe interactions between the (fast moving) massless Nambu-Goldstone pions and other hadrons. In nature quarks and thus pions are not massless and the leading mass corrections are formally of order $`m_\pi /\mathrm{\Lambda }_{\chi SB}`$ where $`\mathrm{\Lambda }_{\chi SB}4\pi F_\pi >1`$ GeV. Of course in an effective field theory the number of terms explodes at higher orders and predictive power is eventually lost, unless an early truncation is possible.

Third generation lattice simulations with sea quarks were limited to unrealistically heavy pions, heavier than about 400 MeV. Only recently masses as low as 180 MeV have become possible in the quenched approximation, albeit not yet at sufficiently large volumes. To allow for a first principles connection with phenomenology it is mandatory to establish an overlap region between lattice simulations and (leading order) chiral perturbation theory. In general the size of this window will depend on the observable in question. GW fermions help a lot since only in this case we have a variant of chiral symmetry at finite lattice spacings while strictly speaking in other formulations (without doublers) a sensitive comparison can only be performed after extrapolation of lattice results to the continuum limit. While even $`400`$ MeV $`\mathrm{\Lambda }_{\chi SB}`$ (modulo the ambiguity of “$``$vs.$`<`$”) such a pion is still doomed to “see” some of the internal structure of the proton, with an inverse charge radius of about 250 MeV: it is doubtful that the quark and gluon nature of QCD can completely be ignored with such a “hard” pion probe. With broadening baryonic wave functions such issues are likely to be more critical in QCD with light sea quarks than they already are in the quenched approximation.

Naïvely, hadron masses are a polynomial in the quark mass, $`m_qm_\pi ^2`$. However, pion loops give rise to non analytic (NA) functional dependence on the quark mass. For instance the nucleon mass is given by,

$$m_N(m_\pi )=m_N(0)+c_2m_\pi ^2+c_3m_\pi ^3+\mathrm{}.$$
(1)
The coefficients of the NA terms are related to phenomenological low energy constants. For instance $`c_3=3g_A^2/(32\pi F_\pi ^2)`$. In the quenched approximation, which neglects disconnected fermion lines, the leading NA (LNA) term is not proportional to $`m_\pi ^3`$ anymore but to $`m_\pi ^2\mathrm{ln}m_\pi `$, known as a chiral log. Only very recently convincing signs of such logarithmic behaviour have been found, albeit not yet with the expected coefficients.

Recently, corrections to standard LNA terms due to massive pions (resulting in a different pion cloud) have been modelled, in particular by the Adelaide group. Unlike the chiral Lagrangian itself, such models cannot be derived in any limit directly from QCD, and the predictive power has to be digested with some caution. Figure 2 reveals that we have not yet made contact with chiral perturbation theory (dotted curve, $`c_3`$ fixed). Nonetheless the data are well fitted by Eq. (1), with $`c_3`$ as a free fit parameter (dashed curve), and the extrapolated value has dangerously tiny statistical errors. The Adelaide cloudy bag model (CBM) suggests a different functional form that screens the LNA term at pion masses close to a bag parameter, $`\mathrm{\Lambda }`$. With a fixed phenomenological value of $`\mathrm{\Lambda }`$ the solid curve is obtained while a three parameter fit yields the dashed-dotted curve. By treating the difference as a systematic uncertainty such models allow us to arrive at more realistic error estimates. Alternatives to the CBM include NJL based as well as Bethe-Salpeter based ideas.

The quenched approximation is sick in several ways. I already mentioned that many unstable particles become stable. It is not always entirely clear when this “bug” is a virtue and in which cases not. A more serious flaw is that unitarity is violated since quarks are not included in a consistent way. This defect, which becomes visible in particular at small quark masses and in the scalar sector, also reflects onto and can be understood in terms of the chiral expansion. For instance chiral logs originate from this source. Another related problem is the absence of the axial anomaly since quarks do not feed back onto the glue. This leads to the spontaneous breaking of a greater flavour symmetry group to start with and results in $`n_F^2`$ instead of $`n_F^21`$ pions: the $`\eta ^{}`$ becomes heavy only at small sea quark masses, with possible implications onto the phenomenology of excited baryon states: the mass splitting between the nucleon ($`N=N^{1/2^+}`$) and the $`N^{}=N^{1/2^{}}`$ for instance is a consequence of chiral symmetry breaking and its size is related to transitions $`N^{}N+\pi N^{}`$. The number of pions should naturally be expected to seriously affect such observables. (Note that also in the quenched approximation there is a pion cloud since valence quarks can travel forwards and backwards and forwards again in time. Pion exchange is possible too, as valence quarks are swapped.)

To this end it has recently been suggested to model quenched LNA terms and to subtract these from the corresponding lattice simulation results, replacing them by the expected un-quenched counterparts. Needless to say that many assumptions are required, in particular relations between low energy parameters appearing in quenched and un-quenched chiral Lagrangians and in the CBM. However, in the absence of lattice results with light sea quarks, playing around with such models provides us with valuable rough ideas of possible sizes of quenching effects. On the other hand, once quality un-quenched data become available this will positively feed back onto such models, providing insights impossible to gain from experiment alone where quark masses and the number of active flavours are fixed.

Determinations of nucleon parton distributions seem to suffer from extremely large uncertainties in the chiral extrapolation; in particular this is so for the leading moments. In a recent study of the LHPC/SESAM Collaborations of quark distributions including sea quarks almost perfect agreement was found with quenched results, obtained both by them and by QCDSF, down to $`m_\pi 600`$ MeV. However, a naïve polynomial extrapolation of $`x_{ud}`$ was found to overestimate the phenomenological value by as much as a factor 5/3. Such a factor can indeed be accounted for by modelling pion mass effects on the LNA contribution in a sensible way, indicating that the quark mass is the dominant (and a huge) source of uncertainty.

## 4 Spectrum

In the past two years there has been a vibrant flow of lattice publications on the spectrum of excited nucleon states. All but one have been based on the quenched approximation so far. Only the LHPC/UKQCD/QCDSF Collaborations, using the Wilson-Clover action, attempted a continuum limit extrapolation. Other strategies were implementations of improved actions like the $`D_{234}`$ action by Lee et al. or the FLIC action by Melnitchouk et al.. The Riken-BNL group employed domain wall fermions to improve the chiral properties.

One feature that all these articles share is a level inversion between the positive parity Roper $`N^{}(1440)`$ resonance and the negative parity $`N^{}(1535)`$, relative to experiment. It was not clear whether this was a quenching effect, hence shared with many models, most of which are essentially quenched too. But maybe the resonance observed in nature did something funny, to which the quark-only operator used on the lattice to create the Roper did not couple well. A pessimist might also argue that such resonances are quite unstable which could very well affect their position in ways that are not easily foreseeable, prior to doing the real thing.

However, new results obtained with GW fermions by the Kentucky-Washingon group at lighter pion masses, $`m_\pi >180`$ MeV, than ever realised before indicate an exciting behaviour: as the pion becomes lighter than about 550 MeV the slope of the Roper data bends downwards, heading straight in the direction of the experimental value (cf. Figure 3). While at the lightest quark mass $`Lm_\pi 2.7`$ and finite size effects are likely, everything with $`m_\pi ^2>0.1\text{GeV}^2`$ is pretty safe. There is overlap with other lattice studies in the region $`m_\pi ^2>0.25\text{GeV}^2`$ and one might wonder why the onset of this behaviour has not been detected before. Inspection shows, however, that within the statistical errors an effect of the expected size would not have appeared significant in these cases.

Certainly, once the independence of the Kentucky results from Bayesian fitting priors and lattice volume has been firmly established, these are very exciting news that tell model builders how transitions between heavy and light quark regions can look like: a prime example for the cross fertilisation between lattice QCD, models and experiment in this thriving area of research. The Kentucky group has investigated many additional baryonic channels and it is astonishing how good quenched QCD seems to explain the positions of these strongly decaying resonances, once the quark masses become sufficiently light, another unsolved miracle of nature.

## 5 Structure

Exciting results on the transition matrix element for electromagnetic $`N`$ to $`\mathrm{\Delta }`$ transitions have been reported. The relevant operator is usually decomposed into three Sachs type form factors, accompanying magnetic dipole transitions, $`G_{M1}(q^2)`$ as well as electric and Coulomb/scalar quadrupole transitions, $`G_{E2}(q^2)`$ and $`G_{C2}(q^2)`$. The quadrupole transitions ratios $`R_{EM}=G_{E2}/G_{M1}`$ and $`R_{SM}=G_{C2}/G_{M1}`$ are related to the question of nucleon deformation. Recent experimental results, in particular from the CLAS and OOPS Collaborations indicate the ratios $`R_{EM}0.02`$ and $`R_{SM}0.07`$ around $`q^20.53`$, which is the lowest momentum transfer that could be realised in a recent lattice study with two flavours of sea quarks, not much lighter than the strange. These first quantitative lattice results seem to get the sign right and point towards values $`R_{EM}=0.03(1)`$ and $`R_{SM}=0.03(2)`$, consistent with experiment. Interestingly, the quenched reference result $`R_{EM}=0.009(8)`$ is compatible with zero, indicating that a more realistic pion cloud might be essential for the effect. Obviously, lattice simulations are still miles away from computing, say, octupole form factors for transitions between orbitally excited nucleons.

In another lattice investigation the electric charge and matter distributions of $`\pi `$, $`\rho `$, $`N`$ and $`\mathrm{\Delta }`$ have been studied at rest. One can easily calculate the charge radius from such a distribution and the Fourier transform is related to form factors at zero momentum transfer. While the (cut-off and scheme dependent) matter distributions of all investigated hadrons were found to be remarkably similar there are differences in the charge distributions. In the quenched approximation the pion is only approximately half as wide as the $`\mathrm{\Delta }`$, with $`N`$ and $`\rho `$ almost of the same size as the $`\mathrm{\Delta }`$. Even with a pion mass of almost 600 MeV un-quenching significantly broadened the latter three distributions, and in particular the $`\mathrm{\Delta }`$, while the $`\pi `$ remains completely unaffected. The quenched proton charge radius comes out to be $`r_p=0.59(4)`$ fm, after polynomial extrapolation to the physical limit. This has to be compared with the value $`r_p0.81`$ fm from dipole fits to experimentally determined electromagnetic form factors. The above direct lattice result is in agreement with $`r_p`$ 0.6 fm obtained indirectly from a dipole fit to electric and magnetic proton form factors calculated on the lattice by the QCDSF Collaboration. However, from recent JLAB Hall A results, showing an energy dependence of the ratio $`G_e/G_m`$, we know that the dipole ansatz is not the whole truth, which stresses the importance of determining quantities like charge radii in more than one way. The same QCDSF work suggests $`r_p0.7`$ fm for two flavours of sea quarks, closer to the phenomenological value.

For the first time deformations of the hadron wave functions have been investigated. Obviously, a spectroscopic quadrupole moment becomes only possible for an angular momentum of at least one, such that deviations from spherical symmetry are ruled out for the nucleon. However, the $`\rho `$ was reported to be somewhat prolate with the $`\mathrm{\Delta }`$ being slightly oblate and the deformations appeared to increase when including sea quarks. In neither of the cases was the deformation of the $`\mathrm{\Delta }`$ statistically significant though.

It is not clear whether QCD forces are stringy with a three body potential that depends on the length of the shortest path connecting the three quarks to one point (Mercedes or $`Y`$ law) or if they merely can be understood in terms of sums of two body interactions ($`\mathrm{\Delta }`$ law). Leaving aside the question whether an instantaneous Born-Oppenheimer interaction picture is appropriate for a situation with light quarks, the nucleon wave function seems to be better described by a Schrödinger equation with a potential that depends on the “$`\mathrm{\Delta }`$-distance” than by the $`Y`$ configuration. This finding is also supported by recent lattice calculations of $`SU(3)`$ and $`SU(4)`$ gauge theory as well as of the $`Z_3`$ Potts model. The results indicate that for interquark distances of physical relevance, i.e. $`r<1`$ fm, the $`\mathrm{\Delta }`$-picture yields a very accurate description while, long after the “string” of QCD with sea quarks is broken, the $`Y`$ picture is eventually approached.

Recently, spin and transversity of the quark contributions to the $`\mathrm{\Lambda }`$ hyperon have been calculated by QCDSF. The fraction $`\mathrm{\Delta }q_\mathrm{\Lambda }`$ of the spin carried by quarks of flavour $`q`$ is given by the forward matrix element of the axial vector current, $`\mathrm{\Lambda }(p^{},s)|\overline{q}\gamma _\mu \gamma _5q|\mathrm{\Lambda }(p,s)=2\left(p^{}p\right)_\mu \mathrm{\Delta }q_\mathrm{\Lambda },`$ i.e. by the lowest moment of the helicity distribution. Assuming $`SU(3)`$ flavour symmetry one can express the $`\mathrm{\Delta }q_\mathrm{\Lambda }`$s in terms of proton spin fractions $`\mathrm{\Delta }q_p`$. Converting the lattice results into the $`\overline{MS}`$ scheme at $`\mu =2`$ GeV yields an axial charge $`g_A=\mathrm{\Delta }u_p\mathrm{\Delta }d_p`$ that falls short of the experimental value $`g_A1.26`$ by 10-20 %, due to quenching and/or due to uncertainties in the chiral extrapolation, such that we cannot expect the result $`\mathrm{\Delta }u_\mathrm{\Lambda }=\mathrm{\Delta }d_\mathrm{\Lambda }=0.02(4)`$ to be consistent with the real world either.

The lattice data on proton and $`\mathrm{\Lambda }`$ spin fractions turn out to be in excellent agreement with $`SU(3)`$ flavour symmetry, an observation that might be independent of sea quarks and certainly consistent with the success in understanding semileptonic decay rates of the $`\mathrm{\Lambda }`$. Combining experimental data obtained for the proton helicity with the systematic error on $`SU(3)_F`$ symmetry suggested by this lattice study, yields the predictions $`\mathrm{\Delta }u_\mathrm{\Lambda }=\mathrm{\Delta }d_\mathrm{\Lambda }=0.17(3)`$ and $`\mathrm{\Delta }s_\mathrm{\Lambda }=0.63(3)`$, again in the $`\overline{MS}`$ scheme at $`\mu =2`$ GeV.

## 6 Outlook

Experimental efforts are paralleled by lattice studies of baryon spectrum and structure. The main limitation that has to be overcome in the near future is the use of unrealistically heavy quarks in quenched simulations, as well as in simulations with sea quarks. Recent theoretical developments and advances in computer technology mean that optimism is justified.

## Acknowledgments

This work is supported by PPARC grants PPA/A/S/2000/00271, PPA/G/O/1998/00559 and EU grant HPRN-CT-2000-00145.


text:
# Susceptibility of the 2D 𝑆=1/2 Heisenberg antiferromagnet with an impurity

## Abstract

We use a quantum Monte Carlo method (stochastic series expansion) to study the effects of a magnetic or nonmagnetic impurity on the magnetic susceptibility of the two-dimensional Heisenberg antiferromagnet. At low temperatures, we find a log-divergent contribution to the transverse susceptibility. We also introduce an effective few-spin model that can quantitatively capture the differences between magnetic and nonmagnetic impurities at high and intermediate temperatures.

Static impurities can be introduced into the CuO2 planes of the high-$`T_c`$ cuprates as a means of probing their electronic correlations and excitations; some unwanted impurities and defects are also always present in “pure” systems. Impurities coupled to a two-dimensional (2D) host system hence constitute an important class of quantum many-body problems. When Cu ions in hole-doped CuO2 planes are substituted with Zn, experiments indicate that local magnetic moments are induced at the Cu sites neighboring the nonmagnetic Zn impurities mahajan . This behavior can be reproduced in spin-gapped Heisenberg antiferromagnets ng , e.g., in systems of weakly coupled spin-$`\frac{1}{2}`$ ladders sandvik ; martins ; sachdev . Although spin models cannot address questions specifically related to a metallic or superconducting host, static impurities in various gapped (paramagnetic) and gapless (antiferromagnetic) Heisenberg systems are important limiting models for understanding the physics of quantum spin defects. They also have direct experimental realizations in the cuprates doping ; greven .

Several studies have addressed static vacancies ng ; sandvik ; sachdev ; bulut ; nagaosa  and added spins igarashi ; kotov  in 2D $`S=1/2`$ Heisenberg models. The localized moments forming in spin-gapped systems are now well understood. Recently, a universal behavior was predicted sachdev  for the $`T>0`$ impurity effects in systems that have long-range order at $`T=0`$ (i.e., in the “renormalized classical” regime chn ), as well as in systems close to a quantum-critical point. Some predictions for quantum-critical systems have been confirmed numerically sachdev2 , whereas other aspects of the theory remain unsettled troyer ; sushkov . The predictions in the renormalized-classical regime have not yet been tested.

In this Letter, we consider the 2D square-lattice spin-$`\frac{1}{2}`$ Heisenberg antiferromagnet with (i) a vacancy (nonmagnetic impurity) or (ii) an off-plane added spin (magnetic impurity) coupled to a single host spin. We perform quantum Monte Carlo calculations, using the stochastic series expansion (SSE) technique sse , and determine the effects of the two different impurities on the uniform magnetic susceptibility. Our results confirm quantitatively a classical-like longitudinal Curie contribution resulting from alignment of the impurity moment with the local Néel order sachdev . However, we also find evidence of a logarithmic divergence as $`T0`$ in the transverse component of the impurity susceptibility, instead of the predicted $`T`$ independence at low $`T`$. We point out qualitative differences between the vacancy and the added spin at high and intermediate $`T`$ and explain these in terms of an effective few-spin model.

We begin by defining three spin-$`\frac{1}{2}`$ Hamiltonians describing the models investigated:

$`H_0`$ $`=`$ $`J{\displaystyle \underset{i,j}{}}𝐒_i𝐒_j,`$ (1a)
$`H_{}`$ $`=`$ $`J{\displaystyle \underset{\genfrac{}{}{0.0pt}{}{i,j}{i,j0}}{}}𝐒_i𝐒_j,`$ (1b)
$`H_+`$ $`=`$ $`J{\displaystyle \underset{i,j}{}}𝐒_i𝐒_j+J_{}𝐒_a𝐒_0,`$ (1c)
where $`i,j`$ denotes a pair of nearest-neighbor sites on a periodic $`L\times L`$ lattice. $`H_0`$ is the standard Heisenberg model. In $`H_{}`$ the spin $`𝐒_0`$ has been removed, creating a vacancy. We also study systems with two vacancies at maximum separation on different sublattices. In $`H_+`$ an added spin-$`\frac{1}{2}`$, $`𝐒_a`$, is coupled to one of the spins, $`𝐒_0`$, in the plane. We here only consider $`J_{}=J`$.

We have calculated the total susceptibilities

$$\chi _k^z=\frac{1}{T}\left(\underset{i=1}{\overset{N}{}}S_i^z\right)^2,$$
(2)
with $`k=0,,+`$ corresponding to $`H_0`$, $`H_{}`$, and $`H_+`$. We have used square lattices with $`L=4,8,16,32`$, and $`64`$, at temperatures down to $`T=J/32`$. The number of spins $`N=L^2`$ for $`k=0`$ and $`L^2\pm 1`$ for $`k=\pm `$. In order to determine the effects of the impurities, we follow Ref. sachdev,  and define the impurity susceptibilities

$$\chi _{\text{imp}}^{z,\pm }=\chi _\pm ^z\chi _0^z.$$
(3)
In the case of two vacancies, the impurity susceptibility is further normalized by a factor $`\frac{1}{2}`$, so that for $`L\mathrm{}`$ it should be the same as for a single vacancy.

The temperature dependence of $`\chi _{\text{imp}}^z`$ was discussed by Sachdev et al. sachdev  on the basis of quantum field theory. They concluded that at $`T=0`$, the component $`\chi _{\text{imp}}^{}`$ parallel to the direction of the Néel order vanishes, while the perpendicular component $`\chi _{\text{imp}}^{}=C_3/\rho _s`$, where $`C_3`$ is a constant and $`\rho _s`$ is the spin stiffness of the host. For $`T>0`$ the correlation length $`\xi `$ grows exponentially as $`T0`$ chn . For a general impurity spin, a moment of magnitude $`S`$ was then predicted to result from the alignment of $`𝐒`$ with a large Néel-ordered domain, i.e., it does not have the quantum mechanical magnitude $`\sqrt{S(S+1)}`$ which might have been naively expected. The full low-$`T`$ impurity susceptibility is hence $`\chi _{\text{imp}}^z=\frac{1}{3}\frac{1}{T}S^2+\frac{2}{3}\chi _{\text{imp}}^{}`$ sachdev .

In order to have a simple system which reproduces this expected low-$`T`$ behavior, and also accounts for non-universal behavior at higher $`T`$ for different types of impurities, we introduce a simple effective few-spin model. The idea is to model a large Néel-ordered region by a classical vector $`𝐍`$. Three Hamiltonians corresponding to the original models (1) are defined:

$`H_0^{\mathrm{eff}}`$ $`=`$ $`\alpha 𝐒_0𝐒_e+r𝐍𝐒_eh_zM_0^z,`$ (4a)
$`H_{}^{\mathrm{eff}}`$ $`=`$ $`r𝐍𝐒_eh_zM_{}^z,`$ (4b)
$`H_+^{\mathrm{eff}}`$ $`=`$ $`\alpha 𝐒_0𝐒_e+r𝐍𝐒_e+J_{}𝐒_a𝐒_0h_zM_+^z.`$ (4c)
The effective Hamiltonian for the vacancy model, $`H_{}^{\mathrm{eff}}`$, contains a single spin $`𝐒_e`$ representing a remnant spin-$`\frac{1}{2}`$ due to the sublattice asymmetry caused by the vacancy. It is coupled to the unit vector $`𝐍`$, representing the orientation of the local Néel order of the host antiferromagnet. The magnitude of this order is absorbed into the effective coupling $`r`$. In $`H_0^{\mathrm{eff}}`$ we “reinsert” the spin $`𝐒_0`$ that was removed in the vacancy model, and couple it to $`𝐒_e`$. Further, including the Heisenberg interaction between the host spin $`𝐒_0`$ and the added spin $`𝐒_a`$, we arrive at the effective Hamiltonian for the added spin model $`H_+^{\mathrm{eff}}`$.

We determine the impurity susceptibilities (3) of the effective models in the same way as for the original models. In Eqs. (4) we have explicitly indicated how an applied external field, which defines the $`z`$ direction, couples to the systems. The magnetization operators are $`M_{}^z=S_e^z`$, $`M_0^z=S_0^z+S_e^z`$, and $`M_+^z=S_a^z+S_0^z+S_e^z`$. The susceptibilities for $`k=0,,+`$, corresponding to $`H_0^{\mathrm{eff}}`$, $`H_{}^{\mathrm{eff}}`$, and $`H_+^{\mathrm{eff}}`$, can be written in the form

$$\chi _k^z=\frac{}{h_z}M_k^z=\frac{1}{3}\chi _k^{}+\frac{2}{3}\chi _k^{},$$
(5)
where $``$ and $``$ refer to the directions parallel and perpendicular to the vector $`𝐍`$. After straight-forward diagonalization of the effective Hamiltonians in the $``$ basis, the two components can be evaluated using

$`\chi _k^{}`$ $`=`$ $`{\displaystyle \frac{1}{T}}(M_k^{})^2,`$ (6a)
$`\chi _k^{}`$ $`=`$ $`{\displaystyle _0^{1/T}}𝑑\tau M_k^{}(\tau )M_k^{}(0).`$ (6b)
In (5) an average over all orientations of $`𝐍`$ relative to the fixed $`z`$ axis has been taken. The coupling of the impurity spin to $`𝐍`$, therefore, and in accordance with Ref. sachdev, , leads to a classical-like low-$`T`$ divergence; $`\chi _{\text{imp}}^{z,\pm }\frac{1}{12}\frac{1}{T}+\frac{2}{3}\chi _\pm ^{}`$, where the second term is constant at low $`T`$ (instead of having the $`r=0`$ form $`\frac{1}{6}\frac{1}{T}`$).

We do not attempt to derive values for the couplings $`\alpha `$ and $`r`$ \[however, $`J_{}=J`$ as in (1c)\]. Moreover, at least $`r`$ should in principle have some $`T`$ dependence. All results presented will be for constant $`\alpha /J=2.1`$ and $`r/J=1.75`$, which give a reasonable over-all agreement with the SSE calculations. We will show that although the effective model is highly simplified, it captures some of the differences between the vacancy and the added spin.

Next, we present the results of the SSE calculations for the models (1) and compare with the corresponding effective models. The impurity susceptibilities of interest, Eq. (3), are defined as differences between two extensive quantities. Although improved estimators evertz  were utilized in the SSE, the relative statistical noise grows very rapidly with increasing $`L`$. We are therefore currently limited to $`L64`$ and $`TJ/32`$.

The single-vacancy impurity susceptibility multiplied by $`4T`$ is shown in Fig. 1. At high $`T`$ the data for different $`L`$ are indistinguishable, while at low $`T`$ finite-size effects are clearly seen for $`L16`$. In the range of $`T`$ considered, all finite-size effects are eliminated within statistical errors for $`L=64`$. The observed behavior at high $`T`$ is easily understood as the total susceptibility is then just the sum of Curie contributions of each independent spin; $`\chi _{\text{imp}}^{z,}(T\mathrm{})=(L^21)/4TL^2/4T=\frac{1}{4}\frac{1}{T}`$. At low $`T`$ the $`S=\frac{1}{2}`$ ground state of $`H_{}`$ and $`S=0`$ of $`H_0`$ lead to the observed Curie behavior $`\chi _{\text{imp}}^{z,}=+\frac{1}{4}\frac{1}{T}`$ for small $`L`$. As $`L`$ grows this finite-size effect vanishes, and we observe the predicted $`\frac{1}{12}\frac{1}{T}`$ contribution sachdev  arising from the longitudinal part, in analogy with the effective model. Hence the nearly constant $`4T\chi _{\text{imp}}^{z,}\frac{1}{3}`$ for $`L=32`$ and $`64`$. The effective model reproduces the behavior reasonably well, as shown in the inset of Fig. 1.

Looking more carefully at the SSE data, the predicted low-$`T`$ constant behavior in the quantity $`\chi _{\text{imp}}^{z,}\frac{1}{12}\frac{1}{T}`$, which should reduce to $`\frac{2}{3}\chi _{\text{imp}}^,`$ as $`T0`$, is not observed. In Fig. 2 (solid circles) we present results for $`\chi _{\text{imp}}^{z,}\frac{1}{12}\frac{1}{T}`$ for $`L=64`$ (data for smaller $`L`$ indicate that there are no significant finite-size effects for $`L=64`$). It shows an apparent logarithmically divergent behavior, roughly from the onset $`T`$ of renormalized-classical behavior in the pure 2D Heisenberg susceptibility, $`T/J0.3`$ kim . This is also approximately where the corresponding effective-model result converges to a constant.

Results for two vacancies are shown in Fig. 3. The high-$`T`$ behavior has the same explanation as for the single vacancy. At low $`T`$, the moments due to the two vacancies, which reside on different sublattices, are pinned by the Néel order antiparallel to each other, resulting in a vanishing $`\chi _{\text{imp}}^,`$. Hence, in this case $`\chi _{\text{imp}}^{z,}`$ does not diverge for small $`L`$ and we do not, therefore, multiply our results by $`T`$. The inset in Fig. 3 shows a comparison between one and two vacancies for $`L=16`$. The $`T`$ at which the two curves deviate corresponds to a correlation length of the same order as the separation between the two impurities, $`\xi L/2`$, i.e., above this $`T`$ the two impurities couple to different Néel domains and behave as independent vacancies. Since $`\xi `$ diverges exponentially, the point of deviation moves very slowly towards $`T=0`$ as $`L`$ grows. At low $`T`$ the resulting $`\chi _{\text{imp}}^{z,}`$ for two vacancies shows little $`T`$ dependence for large $`L`$. However, no sign of convergence of the plateau value is seen as the system size grows. If $`\chi _{\text{imp}}^,`$ is finite as $`T0`$, we would expect such a convergence at low $`T`$ (with a peak at intermediate $`T`$ for very large $`L`$, where the cut-off of the $`\frac{1}{12}\frac{1}{T}`$ divergence occurs at low $`T`$). It should be noted that in a finite system there will always be some interactions also between the $``$ components of the two impurity spins at low $`T`$. We can therefore not expect the behavior for two vacancies to be given in a straight-forward way by the single-vacancy results in Fig. 2. The roughly $`\mathrm{ln}(L)`$ divergence of the plateau height in Fig. 3 is, however, fully in line with a log-divergent $`\chi _{\text{imp}}^,`$ for a single vacancy.

Results for $`\chi _{\text{imp}}^{z,+}`$ for systems with an off-plane added spin are shown in Fig. 4. The high-$`T`$ behavior, as well as the low-$`T`$ behavior for small systems ($`L=4,8`$), is understood by the same arguments as for the vacancy. Note that in this case the high-$`T`$ impurity susceptibility is $`+\frac{1}{4}\frac{1}{T}`$, instead of $`\frac{1}{4}\frac{1}{T}`$ for the vacancy. Again, we believe that all finite-size effects are eliminated for $`L=64`$ down to $`T=J/32`$. Contrary to the behavior for the vacancy in Fig. 1, the $`T`$ at which $`4T\chi _{\text{imp}}^{z,+}`$ assumes an almost constant value $`\frac{1}{3}`$ has not yet been reached at $`T=J/32`$. Moreover, we note the substantial differences between $`\chi _{\text{imp}}^{z,+}`$ and $`\chi _{\text{imp}}^{z,}`$ at intermediate $`T`$. In particular, the shoulder seen in Fig. 4 for $`T0.21.2`$ has no counterpart in Fig. 1. This feature is clearly due to the internal structure of the added spin impurity and is reproduced very well by the effective model. In Fig. 2 (open circles) we present the results for $`\chi _{\text{imp}}^{z,+}\frac{1}{12}\frac{1}{T}`$ for $`L=64`$. As in the vacancy case, it appears to be log-divergent at low $`T`$. Note also that the effective model describes the added spin impurity very well down to quite low temperatures. It captures the behavior better than for the vacancy, but a better agreement for the vacancy can also be achieved by using slightly different $`\alpha `$ and $`r`$.

Summarizing our results, we have confirmed the predicted sachdev  contribution $`\frac{1}{12}\frac{1}{T}`$ to the impurity susceptibility. However, we also find a logarithmic divergence of $`\chi _{\text{imp}}^{z,\pm }\frac{1}{12}\frac{1}{T}`$, instead of the constant $`\frac{2}{3}\chi _{\text{imp}}^{}`$ predicted for this quantity at low $`T`$ (and reproduced with our effective model). We cannot, of course, completely exclude an approach to a constant at still lower temperatures, but such a slow convergence had also not been anticipated. We note that the separation of $`\chi _{\text{imp}}^z`$ into transverse and longitudinal components, Eq. (5), is strictly correct for the 2D Heisenberg model only at $`T=0`$. However, the very sudden cutoff of the divergence seen in the two-vacancy data in Fig. 3 supports the notion of a component aligning very strongly to the local Néel order (which becomes the global order at the $`L`$-dependent crossover $`T`$) and justifies the relation $`\frac{2}{3}\chi _{\text{imp}}^{}=\chi _{\text{imp}}^z\frac{1}{12}\frac{1}{T}`$ also at relatively high temperatures. Hence, our results are most naturally interpreted as a log divergent $`\chi _{\text{imp}}^{}`$. This is in fact in line with a Green’s function calculation by Nagaosa et al. nagaosa . They found that for a system with a vacancy, the frequency dependent impurity susceptibility at $`T=0`$, $`\chi _{\text{imp}}^{}(T=0,\omega )`$, was log divergent when $`\omega 0`$. In view of the renormalized-classical chn  picture, this is consistent with exact results harris  for the classical 2D Heisenberg model. An anomalous perpendicular susceptibility was also recently noted for the $`S=1/2`$ model at finite impurity concentration chernyshev .

The results presented here call for a reexamination of the field theory sachdev  of quantum impurities in the renormalized-classical regime. Very recent efforts to explain the log divergence, motivated by our numerical findings, have indicated that the impurity moment acquires a previously unnoticed correction of order $`T\mathrm{ln}(1/T)`$ to its leading order value $`S`$ in the renormalized-classical regime sachdev3 ; sushkov2 , and that this can be interpreted as a log-divergent contribution to $`\chi _{\mathrm{imp}}^{}`$ as proposed here (results in the quantum-critical regime sachdev  remain unchanged).

Apart from the log divergence, the effective few-spin model that we have introduced here gives a good description of the impurity susceptibility at high and intermediate $`T`$. In particular, it captures very well the nonmonotonic $`T`$ dependence that we have found in the case of an added spin.

We would like to thank Oleg Sushkov, and, in particular, Subir Sachdev for very useful correspondence and comments on the manuscript. This work is supported by the Academy of Finland, project No. 26175.