Correlation functions and excitation spectrum of the frustrated ferromagnetic spin-1
chain in an external magnetic field
T. Vekua,1 A. Honecker,2 H.-J. Mikeska,3 and F. Heidrich-Meisner4
Laboratoire de Physique Théorique et Modèles Statistiques,
Université Paris Sud, 91405 Orsay Cedex, France
Institut für Theoretische Physik, Universität Göttingen, 37077 Göttingen, Germany
Institut für Theoretische Physik, Universität Hannover, Appelstrasse 2, 30167 Hannover, Germany
Materials Science and Technology Division, Oak Ridge National Laboratory, Tennessee, 37831, USA and
Department of Physics and Astronomy, University of Tennessee, Knoxville, Tennessee 37996, USA
(Dated: April 5, 2007; revised: July 6, 2007)
Magnetic field effects on the one-dimensional frustrated ferromagnetic chain are studied by means
of effective field theory approaches in combination with numerical calculations utilizing Lanczos
diagonalization and the density matrix renormalization group method. The nature of the ground
state is shown to change from a spin-density-wave region to a nematic-like one upon approaching
the saturation magnetization. The excitation spectrum is analyzed and the behavior of the single
spin-flip excitation gap is studied in detail, including the emergent finite-size corrections.
I. INTRODUCTION
The interest in helical and chiral phases of frustrated
low-dimensional quantum magnets has been triggered by
recent experimental results. While many copper-oxide
based materials predominantly realize antiferromagnetic
exchange interactions, several candidate materials with
magnetic properties believed to be described by frus-
trated ferromagnetic chains have been identified,1,2,3,4,5,6
including Rb2Cu2Mo3O12 (Ref. 1), LiCuVO4 (Refs. 2,
3,4,5), and Li2ZrCuO4 (Ref. 6). The frustrated anti-
ferromagnetic chain is well-studied,7 but the magnetic
phase diagram of the model with ferromagnetic nearest-
neighbor interactions remains a subject of active theoret-
ical investigations.8,9,10,11
In this work we consider a parameter regime that is
in particular relevant for the low-energy properties of
LiCuVO4, corresponding to a ratio of J1 ≈ −0.3 J2 be-
tween the nearest neighbor interaction J1 and the frus-
trating next-nearest neighbor interaction J2 > 0. As
the interchain couplings for this material are an order
of magnitude smaller than the intrachain ones,3 we an-
alyze a purely one-dimensional (1D) model. Apart from
mean-field based predictions,8 the nature of the ground
state in a magnetic field h is not yet completely known.
Therefore, combining the bosonization technique with a
numerical analysis we determine ground-state properties
and discuss the model’s elementary excitations.
The Hamiltonian for our 1D model reads:
J1~Sx · ~Sx+1 + J2~Sx · ~Sx+2
Szx , (1)
where ~Sx represents a spin one-half operator at site x.
Bosonization has turned out to be the appropriate lan-
guage for describing the regime |J1| ≪ J2 of Eq. (1).
This result has been established by studying the magne-
tization process yielding a good agreement between field
theory and numerical data.9 The derivation of the ef-
fective field theory is summarized in Sec. II. Here, we
extend on such comparison of analytical and numerical
results and further confirm the predictions of field the-
ory by analyzing several correlation functions in Sec. III.
Then, in Sec. IV, we numerically compute the one- and
two-spin flip excitation gaps and compare them to field-
theory predictions. Finally, Sec. V contains a summary
and a discussion of our results.
II. EFFECTIVE FIELD THEORY
We start from an effective field theory describing the
long-wavelength fluctuations of Eq. (1). In the limit of
strong next-nearest neighbor interactions J2 ≫ |J1|, the
spin operators can be expressed as:
Szα(r) ∼ m+ c(m) sin
2kF r +
+ · · ·
S−α (r) ∼ (−1)re−iθα
π + · · · . (2)
kF = (
−m)π is the Fermi-wave vector and α = 1, 2 enu-
merates the two chains of the zig-zag ladder. In relation
with Eq. (1), note that ~S1(r) = ~S(x+1)/2 (~S2(r) = ~Sx/2)
for x odd (even). φα and θα are compactified quantum
fields describing the out-of-plane and in-plane angles of
fluctuating spins obeying Gaussian Hamiltonians:
H = v
(∂xφα)
2 +K(∂xθα)
, (3)
with [φα(x), θα(y)] = iΘ(y − x), where Θ(x) is the
Heaviside function. Sub-leading terms are suppressed in
Eq. (2). m is the magnetization of decoupled chains, re-
lated to the real magnetization M of the zig-zag system
M ≃ m
1− 2K(m)J1
πv(m)
. (4)
K(m) and v(m) are the Luttinger liquid (LL) parame-
ter and the spin-wave velocity of the decoupled chains,
http://arxiv.org/abs/0704.0764v2
respectively. The nonuniversal amplitude c(m) appear-
ing in the bosonization formulas (2) has been determined
from density matrix renormalization group (DMRG)
calculations.12 Note that in our notation M = 1/2 at
saturation.
Now we perturbatively add the interchain coupling
term to two decoupled chains, each of which is described
by an effective Hamiltonian of the form Eq. (3) and fields
φi and θi, i = 1, 2. For convenience, we transform to
the symmetric and antisymmetric combinations of the
bosonic fields φ± = (φ1±φ2)/
2 and θ± = (θ1±θ2)/
In this basis and apart from terms H±0 of the form (3),
the effective Hamiltonian describing low-energy proper-
ties of Eq. (1) contains a single relevant interaction term
with the bare coupling g1 ∝ J1 ≪ v:
Heff = H+0 +H
0 + g1
dx cos
, (5)
and the renormalized LL parametersK± are, in the weak
coupling limit:
K± = K
1∓ J1
. (6)
K+ is the Luttinger-liquid parameter of the soft mode
of the zig-zag ladder. The Hamiltonian (5) represents
the minimal effective low-energy field theory describing
the region J2 ≫ |J1| of the frustrated FM spin-1/2 chain
for M 6= 0.9,13 The relevant interaction term cos
opens a gap in the φ− sector. Since S
x+1 − Szx ∼ ∂xφ−,
relative fluctuations of the two chains are locked. This
implies that single-spin flips are gapped with a sine-
Gordon gap in the sector describing relative spin fluctua-
tions of the two-chain system.9 Gapless excitations come
from the ∆Sz = 2 channel, i.e. only those excitations
are soft where spins simultaneously flip on both chains.
DMRG results show that this picture applies to a large
part of the magnetic phase diagram.9
III. CORRELATION FUNCTIONS
We now turn to the ground state properties of Eq. (1)
as a function of magnetization, concentrating on several
correlation functions in order to identify the leading in-
stabilities. Note that our analysis is only valid if M 6= 0.
Apart from a term representing the magnetization M in-
duced by the external field, the longitudinal correlation
function shows an algebraic decay with distance r:
〈Szα(0)Szβ(r)〉 ≃ M2+
C1 cos(2kF r + (α− β)kF )
2π2rK+
8π2r2
The constants Ci, i = 1, 2, 3, appearing here and in
Eq. (9) will be determined through a comparison with
numerical results.
In contrast to Eq. (7), the transverse xy-correlation
functions decay exponentially reflecting the gapped na-
ture of the single spin-flip excitations. Here we do not
restrict ourselves to the equal-time expression only, be-
cause we will need non-equal time correlation functions
to extract the finite-size corrections to the gap later on.
We obtain:
〈S+α (0, 0)S−β (r, τ)〉 ≃
δα,β(−1)re−∆1(M)
τ2+r2/v2
(r2 + v2+τ
8K+ (r2 + v2−τ
where τ stands for the Euclidean time, ∆1(M) is the
∆Sz = 1 gap, and v± ∼ v ± J1/π in the weak cou-
pling limit. The Kronecker delta strictly applies to the
thermodynamic limit, while on the lattice an additional
contribution for α 6= β exists.
It is noteworthy that, different from Eq. (8), the in-
plane correlation functions involving bilinear spin combi-
nations decay algebraically. This stems from the gapless
nature of ∆Sz = 2 excitations. In fact, these are the
slowest decaying correlators close to the saturation mag-
netization:
〈S+1 (r)S
2 (r)S
1 (0)S
2 (0)〉 ≃
r1/K+
C3 cos(2kF r)
rK++1/K+
This result is reminiscent of a partially ordered state be-
cause the ordering tendencies in this correlation function
are more pronounced than those of the corresponding
single-spin correlation function Eq. (8). Therefore, we
call the correlator (9) ‘nematic’. Furthermore, we will
refer to a situation where Eq. (9) is the slowest decay-
ing one among all correlation functions as a ‘nematic-like
phase’.
By virtue of the exponential decay in (8), the correlator
(9) is proportional to:
〈(S+1 (r) + S
2 (r))
2 (S−1 (0) + S
2 (0))
2〉 . (10)
The term (Sα1 + S
2 appearing in the case of the S = 1
zig-zag ladder corresponds to the operator (Sα)2 in the
case of a S = 1 chain. One can think of an effective S = 1
spin formed from two neighboring S = 1
spins coupled
by the ferromagnetic interaction. A similar behavior of
correlation functions, namely the exponential decay of in-
plane spin components and the algebraic decay of their
bilinear combinations, is encountered also in the XY 2
phase of the anisotropic S = 1 chain14 and in the spin-1
chain with biquadratic interactions, see, e.g., Ref. 15.
The algebraic decay of the nematic correlator as op-
posed to the exponential decay of (8) suggests that there
are tendencies towards nematic ordering in this phase.
Depending on the value of K+ the dominant instabilities
are either spin-density-wave ones for K+ < 1 or nematic
ones for K+ > 1. From the result for K+ given in Eq. (6)
one can perturbatively evaluate the crossover value of J1:
|J1,cr| =
πv(m)
. (11)
For J1 < J1,cr the nematic correlator (9) is the slowest
decaying one, i.e. one is in the nematic-like phase. The
-0.015
-0.01
-0.005
0.005
0.001
0 4 8 12 16 20 24 28 32
fit, L=32
fit, L=48
fit, L=64
ED, L=32
ED, L=48
ED, L=64
FIG. 1: (Color online) Correlation functions at J1 = −J2 < 0,
and magnetization M = 3/8: (a) longitudinal component Sz
(b) transverse component S±
, (c) spin nematic S±
. x is
the distance in a single-chain notation. ED results for periodic
boundary conditions are shown by symbols, fits by lines. Note
the logarithmic scale of the vertical axis in panel (b).
behavior of the cross-over line can be read off from the
behavior of K(m): K(m) increases monotonically with
m, tends to K = 1 for m → 1/2, and satisfies K < 1
for m < 1/2 (see, e.g., Refs. 16,17). Therefore, we have
J1,cr = 0 for M = 1/2 with increasing ferromagnetic
|J1,cr| for decreasing M . This means that for J1 < 0
a regime opens at high M where nematic correlations
given by Eq. (9) dominate over spin-density-wave corre-
lations given by Eq. (7), in agreement with Chubukov’s
prediction.8
Now we check the correlation functions obtained within
bosonization against exact diagonalization (ED) results.
Numerical data obtained for J1 = −J2 < 0 and M = 3/8
on finite systems with periodic boundary conditions are
shown in Fig. 1. This parameter set allows for a clear
test of the above predictions, but represents the generic
behavior in the phase of two weakly coupled chains. To
take into account finite-size effects we use the observation
that for a conformally invariant theory, any power law
on a plane becomes a power law in the following variable
defined on a cylinder of circumference L:
x → L
. (12)
First we fit the nematic correlator given by Eq. (9),
which from bosonization is expected to be the leading
instability at high magnetizations. Using the part with
x ≥ 5 of the L = 64 data shown in Fig. 1c, we find
1/K+ = 0.904 ± 0.011, C2 = 0.143 ± 0.004, and C3 =
−0.326±0.013. Fig. 1c shows that all finite-size results for
the nematic correlator are nicely described by this fit with
the dependence on L taken into account by substituting
Eq. (12) for the power laws. Moreover, from K+ > 1, we
see that the system is indeed in the region dominated by
nematic correlations for M = 3/8 and J1 = −J2.
Now we turn to the longitudinal correlation function
which we fit to the bosonization result Eq. (7). Since
most numerical parameters have been determined by the
previous fit, only one free parameter is left which we de-
termine from the numerical results of Fig. 1a for L = 64
and x ≥ 14 as C1 = 0.060± 0.004. Predictions for other
system sizes are again obtained by substituting Eq. (12)
for the power laws. The agreement in Fig. 1a is not as
good as in Fig. 1c. However, it improves at larger dis-
tances x and system sizes L, indicating that corrections
omitted in Eq. (7) are still relevant on the length scales
considered here.
Finally, the xy-correlation function is shown in Fig. 1b
with a logarithmic scale of the vertical axis of this panel.
The exponential decay predicted by Eq. (8) is verified.
One further observes that correlations between the in-
plane spin-operators belonging to different chains (odd
x) are an order of magnitude smaller than on the same
chain (even x). This suppression of correlations between
different chains corresponds to the δ symbol in (8), which
strictly applies only in the thermodynamic limit and for
large distances.
We summarize the main result of this section: in-plane
spin correlators are exponentially suppressed for any fi-
nite value of the magnetization in the parameter region
|J1| < J2. The ground state crosses over from a spin-
density-wave dominated to a nematic-like phase with in-
creasing magnetic field, with the crossover line given by
Eq. (11).
IV. EXCITATIONS
We next address the excitation spectrum. Since the
gap to ∆Sz = 1 excitations should be directly accessible
to microscopic experimental probes such as inelastic neu-
tron scattering or nuclear magnetic resonance, we analyze
its behavior as a function of magnetization. Sufficiently
below the fully polarized state the gap can be calculated
analytically using results from sine-Gordon theory. In
addition, to leading order of the interchain coupling, one
can get qualitative expressions using dimensional argu-
ments for the perturbed conformally invariant model:
∆1(m) ∼
c2(m)|J1| sin(πm)
v(m)(1 − J1K(m)/πv(m))
, (13)
where ν = 2− 2K(m)
1+ J1K(m)/πv(m)
. m(h),K(h)
and v(h) can be determined numerically from the Bethe
ansatz integral equations.16,17,18,19
With this information and Eqs. (4) and (13) we de-
termine the qualitative behavior of the single-spin gap
∆1(M) as a function of M : it increases from zero at
zero magnetization, reaches a maximum at intermediate
magnetization values, then shows a minimum and, upon
approaching the saturation magnetization, it increases
again. As our formulas do not strictly apply at m = 0,
the notion of a vanishing gap at zero magnetization may
be a spurious result. Note that when the fully polarized
state is approached, the magnetization increases in an un-
physical fashion since in this limit bosonization becomes
inapplicable. At the point where the magnetization sat-
urates the exact value of the gap can be obtained from
the following mapping to hard-core bosons:8,11,20
Szi =
− a†iai , S
i = a
i . (14)
Comparing Eq. (14) with Eq. (2) one recognizes the lead-
ing terms in Haldane’s harmonic fluid transformation for
bosons.20 Using a ladder approximation which is exact in
the two-magnon subspace we arrive at:
4J22 − 2J1J2 − J21
2(J2 − J1)
J21 + 8J1J2 + 16J
J21 (3 J2 + J1)
J2 (J2 − J1)
. (15)
In Eq. (15) we have represented the gap as a difference
of two terms: the quantum and the classical instability
fields emphasizing its quantum origin.
In order to verify these field theory predictions, we
perform complementary numerical computations using
the DMRG method.21 Open boundary conditions are im-
posed and we typically keep up to 400 DMRG states.
From DMRG we obtain the ground-state energies E(Sz)
as a function of total Sz. For those values of Sz that
emerge as a ground state in an external magnetic field
we compute the single-spin excitation gap from
∆1(M) =
E(Sz + 1) + E(Sz − 1)− 2E(Sz)
. (16)
Fig. 2a shows numerical results for ∆1 at a selected value
of J1 = −0.3 J2 < 0 for the largest system sizes investi-
gated. We find that the finite-size behavior of the gap
∆1(M,L) for system sizes L ≥ 24 is well described by
a 1/L correction. This will be further corroborated by
field-theoretical arguments outlined below. Therefore, we
0 0.1 0.2 0.3 0.4 0.5
L=120
L=144
L=156
L=120
L=144
L=156
extrapolated
FIG. 2: (Color online) Density matrix renormalization group
results for the gaps at J1 = −0.3 J2 < 0 as a function of
magnetization M . Panel (a) shows the single-spin excitation
gap (16), panel (b) the finite-size gap (19) for two flipped
spins multiplied by the chain length L.
extrapolate it to the thermodynamic limit using a fit to
the form
∆1(M,L) = ∆1(M) +
+ · · · , (17)
allowing for an additional 1/L2 correction for those values
of M where at least 4 different system sizes are available.
This extrapolation is represented by the full circles in
Fig. 2a; errors are estimated not to exceed the size of the
symbols. Our extrapolation for ∆1 is consistent with a
vanishing gap at M = 0 in agreement with previous nu-
merical studies22 although bosonization predicts a non-
zero – possibly very small – gap.13,22,23 The behavior of
∆1(M) confirms the picture described above: the gap is
non-zero for M > 0, goes first through a maximum and
then a minimum and finally approaches ∆1/J2 ≈ 0.023
given by Eq. (15) for M → 1/2.
We further wish to point out that for chains with pe-
riodic boundary conditions, the coefficient a(M) of the
finite-size extrapolation Eq. (17) is determined by the
spin-wave velocity and the critical exponent of the soft
mode from the ∆Sz = 2 channel. Indeed, using Eq. (8)
where we can set r = 0, and use the conformal mapping
(12) to the cylinder, we see that the leading finite-size
correction to the gap is:
∆1(M,L) = ∆1(M) +
πv+(M)
4K+(M)
. (18)
0 π/2 π
wave vector
1.0 0.0
FIG. 3: (Color online) Numerical dispersion spectrum in the
subspaces of odd Sz computed for L = 24 and J1 = −J2 < 0.
The wave vector is given relative to the ground state wave
vector (0 for Sz = 0, 4, 8 and π for Sz = 2, 6).
Note that we have to replace sin with sinh in Eq. (12)
in order to extract a gap, since we are dealing with Eu-
clidean time. In addition we used the fact that in our ap-
proximation the effective Hamiltonian (5) is a direct sum
of symmetric and antisymmetric sectors. Moreover, it is
only the symmetric sector enjoying conformal invariance
and consequently we perform the replacement τ → sinh τ
only in the symmetric sector. The antisymmetric sector
has a spectral gap and its contribution to the finite-size
corrections of the single-spin flip excitation energy are
exponentially suppressed with system size.24 With this
method one cannot fix the amplitudes of the 1/L2 term
and beyond. Note furthermore that there may be ad-
ditional surface terms for open boundary conditions as
employed in the numerical DMRG computations. Nev-
ertheless there is a dominant 1/L correction in any case.
Next, we briefly look at the ∆Sz = 2 excitations. Their
finite-size gap is, in analogy to Eq. (16), computed with
DMRG from
∆2(M) =
E(Sz + 2) + E(Sz − 2)− 2E(Sz)
. (19)
Fig. 2b shows numerical results for L∆2(M,L) again at
the value J1 = −0.3 J2 < 0. One observes that the scaled
finite-size gaps collapse onto a single curve which shows
that ∆2(M,L) scales linearly to zero with 1/L, exactly
as expected for gapless excitations in 1D. Furthermore,
we observe that the scaled quantity L∆2(M,L) vanishes
as one approaches saturation M = 1/2 which indicates a
vanishing of the velocity of the corresponding excitations
at saturation.
We proceed by discussing the wave-vector dependence
of the ∆Sz = 1 excitation, while we remind the reader
that the low-energy excitations are in the ∆Sz = 2 sec-
tor. Fig. 3 shows representative ED results obtained for
rings with L = 24 and J1 = −J2 < 0. For ground
states with low Sz, the ∆Sz = 1 excitation spectrum
looks similar to the continuum of spinons. On the other
hand, close to saturation one has single-magnon excita-
tions with a minimum given by the classical value of
the wave vector kcl = arccos(|J1|/4J2).8,13,25 We read
off from Fig. 3 that upon lowering the magnetic field,
this minimum shifts from the classical incommensurate
value towards π/2, i.e. the value appropriate for two de-
coupled chains. This renormalization of the minimum of
the magnon excitations towards the value of decoupled
chains can be interpreted in terms of quantum fluctua-
tions, which are enhanced when the density of magnons
increases. A strong quantum renormalization of the pitch
angle from its classical value at zero magnetization was
previously observed by the coupled-cluster method and
DMRG calculations.26
V. SUMMARY
We have combined numerical techniques with analyti-
cal approaches and mapped out the ground state phase
diagram of the frustrated ferromagnetic spin chain in an
external magnetic field. We have established that with
increasing magnetic field, the ground state crosses over
from a spin-density-wave dominated to a nematic-like
phase. Single spin flip excitations are gapped, giving rise
to an exponential decay of in-plane spin correlation func-
tions in both regimes. We have studied the single- and
two-spin flip excitation energy numerically. Using tools
from conformal field theory we have further shown that
the amplitude of the leading 1/L correction term to the
single-spin flip gap is determined by the critical exponent
and the spin-wave velocity of the soft mode.
Finally, in order to apply our findings to the mate-
rial LiCuVO4, one should take into account interchain
interactions as well as anisotropies, which are expected
to be present in this system.3 At low fields, a helical state
has been observed experimentally.2,3 On the other hand,
for the purely one dimensional case, we have shown that
upon increasing the magnetic field there is a competition
between spin-density-wave and nematic-like tendencies.
Those are the two leading instabilities at high magnetiza-
tions and thus they are the natural candidates to become
long-range ordered in higher dimensions. The question
whether there are true phase transitions at high fields
in higher dimensions is beyond the scope of the current
work.
Acknowledgments
We thank A. Feiguin for providing us with his DMRG
code used for large scale calculations. Most of T.V.’s
work was done during his visits to the Institutes of
Theoretical Physics at the Universities of Hannover and
Göttingen, supported by the Deutsche Forschungsge-
meinschaft. The hospitality of the host institutions is
gratefully acknowledged. T.V. also acknowledges support
from the Georgian National Science Foundation under
grant N 06−81−4−100. LPTMS is a mixed research unit
8626 of CNRS and University Paris-Sud. A.H. is sup-
ported by the Deutsche Forschungsgemeinschaft (Project
No. HO 2325/4-1), and F.H.-M. is supported by NSF
grant No. DMR-0443144.
1 M. Hase, H. Kuroe, K. Ozawa, O. Suzuki, H. Kitazawa, G.
Kido, and T Sekine, Phys. Rev. B 70, 104426 (2004).
2 B. J. Gibson, R. K. Kremer, A. V. Prokofiev, W. Assmus
and G. J. McIntyre, Physica B 350, e253 (2004).
3 M. Enderle, C. Mukherjee, B. F̊ak, R. K. Kremer, J.-M.
Broto, H. Rosner, S.-L. Drechsler, J. Richter, J. Malek,
A. Prokofiev, W. Assmus, S. Pujol, J.-L. Raggazzoni, H.
Rakoto, M. Rheinstädter, and H. M. Rønnow, Europhys.
Lett. 70, 237 (2005).
4 M. G. Banks, F. Heidrich-Meisner, A. Honecker, H.
Rakoto, J.-M. Broto, and R. K. Kremer, J. Phys.: Cond.
Mat. 19, 145227 (2007).
5 N. Büttgen, H.-A. Krug von Nidda, L. E. Svistov, L. A.
Prozorova, A. Prokofiev, and W. Aßmus, Phys. Rev. B 76,
014440 (2007).
6 S.-L. Drechsler, O. Volkova, A. N. Vasiliev, N. Tristan, J.
Richter, M. Schmitt, H. Rosner, J. Málek, R. Klingeler, A.
A. Zvyagin, and B. Büchner, Phys. Rev. Lett. 98, 077202
(2007).
7 H.-J. Mikeska and A. K. Kolezhuk, Lect. Notes Phys. 645,
1 (2004).
8 A. V. Chubukov, Phys. Rev. B 44, 4693 (1991).
9 F. Heidrich-Meisner, A. Honecker, and T. Vekua, Phys.
Rev. B 74, 020403(R) (2006).
10 D. V. Dmitriev, V. Y. Krivnov, and J. Richter, Phys. Rev.
B 75, 0114424 (2007).
11 R. O. Kuzian and S.-L. Drechsler, Phys. Rev. B 75, 024401
(2007).
12 See, e.g., T. Hikihara and A. Furusaki, Phys. Rev. B 69,
064427 (2004).
13 D. C. Cabra, A. Honecker, and P. Pujol, Eur. Phys. J. B
13, 55 (2000).
14 H. J. Schulz, Phys. Rev. B 34, 6372 (1986).
15 A. Läuchli, G. Schmid, and S. Trebst, Phys. Rev. B 74,
144426 (2006).
16 K. Totsuka, Phys. Lett. A 228, 103 (1997).
17 D. C. Cabra, A. Honecker, and P. Pujol, Phys. Rev. B 58,
6241 (1998).
18 V. E. Korepin, N. M. Bogoliubov, and A. G. Izergin, Quan-
tum Inverse Scattering Method and Correlation Functions
(Cambridge University Press, Cambridge, England, 1993).
19 S. Qin, M. Fabrizio, L. Yu, M. Oshikawa, and I. Affleck,
Phys. Rev. B 56, 9766 (1997).
20 F. D. M. Haldane, Phys. Rev. Lett 47, 1840 (1981).
21 S. R. White, Phys. Rev. Lett. 69, 2863 (1992); Phys. Rev.
B 48, 10345 (1993).
22 C. Itoi and S. Qin, Phys. Rev. B 63, 224423 (2001).
23 A. A. Nersesyan, A. O. Gogolin, and F. H. L. Eßler, Phys.
Rev. Lett. 81, 910 (1998).
24 S. I. Matveenko and A. M. Tsvelik, private communication.
25 C. Gerhardt, K.-H. Mütter, and H. Kröger, Phys. Rev. B
57, 11504 (1998).
26 R. Bursill, G. A. Gehring, D. J. J. Farnell, J. B. Parkinson,
T. Xiang, and C. Zeng, J. Phys. Cond. Mat. 7, 8605 (1995).

  Magnetic field effects on the one-dimensional frustrated ferromagnetic chain
are studied by means of effective field theory approaches in combination with
numerical calculations utilizing Lanczos diagonalization and the density matrix
renormalization group method. The nature of the ground state is shown to change
from a spin-density-wave region to a nematic-like one upon approaching the
saturation magnetization. The excitation spectrum is analyzed and the behavior
of the single spin-flip excitation gap is studied in detail, including the
emergent finite-size corrections.

Correlation functions and excitation spectrum of the frustrated ferromagnetic spin-1
chain in an external magnetic field
T. Vekua,1 A. Honecker,2 H.-J. Mikeska,3 and F. Heidrich-Meisner4
Laboratoire de Physique Théorique et Modèles Statistiques,
Université Paris Sud, 91405 Orsay Cedex, France
Institut für Theoretische Physik, Universität Göttingen, 37077 Göttingen, Germany
Institut für Theoretische Physik, Universität Hannover, Appelstrasse 2, 30167 Hannover, Germany
Materials Science and Technology Division, Oak Ridge National Laboratory, Tennessee, 37831, USA and
Department of Physics and Astronomy, University of Tennessee, Knoxville, Tennessee 37996, USA
(Dated: April 5, 2007; revised: July 6, 2007)
Magnetic field effects on the one-dimensional frustrated ferromagnetic chain are studied by means
of effective field theory approaches in combination with numerical calculations utilizing Lanczos
diagonalization and the density matrix renormalization group method. The nature of the ground
state is shown to change from a spin-density-wave region to a nematic-like one upon approaching
the saturation magnetization. The excitation spectrum is analyzed and the behavior of the single
spin-flip excitation gap is studied in detail, including the emergent finite-size corrections.
I. INTRODUCTION
The interest in helical and chiral phases of frustrated
low-dimensional quantum magnets has been triggered by
recent experimental results. While many copper-oxide
based materials predominantly realize antiferromagnetic
exchange interactions, several candidate materials with
magnetic properties believed to be described by frus-
trated ferromagnetic chains have been identified,1,2,3,4,5,6
including Rb2Cu2Mo3O12 (Ref. 1), LiCuVO4 (Refs. 2,
3,4,5), and Li2ZrCuO4 (Ref. 6). The frustrated anti-
ferromagnetic chain is well-studied,7 but the magnetic
phase diagram of the model with ferromagnetic nearest-
neighbor interactions remains a subject of active theoret-
ical investigations.8,9,10,11
In this work we consider a parameter regime that is
in particular relevant for the low-energy properties of
LiCuVO4, corresponding to a ratio of J1 ≈ −0.3 J2 be-
tween the nearest neighbor interaction J1 and the frus-
trating next-nearest neighbor interaction J2 > 0. As
the interchain couplings for this material are an order
of magnitude smaller than the intrachain ones,3 we an-
alyze a purely one-dimensional (1D) model. Apart from
mean-field based predictions,8 the nature of the ground
state in a magnetic field h is not yet completely known.
Therefore, combining the bosonization technique with a
numerical analysis we determine ground-state properties
and discuss the model’s elementary excitations.
The Hamiltonian for our 1D model reads:
J1~Sx · ~Sx+1 + J2~Sx · ~Sx+2
Szx , (1)
where ~Sx represents a spin one-half operator at site x.
Bosonization has turned out to be the appropriate lan-
guage for describing the regime |J1| ≪ J2 of Eq. (1).
This result has been established by studying the magne-
tization process yielding a good agreement between field
theory and numerical data.9 The derivation of the ef-
fective field theory is summarized in Sec. II. Here, we
extend on such comparison of analytical and numerical
results and further confirm the predictions of field the-
ory by analyzing several correlation functions in Sec. III.
Then, in Sec. IV, we numerically compute the one- and
two-spin flip excitation gaps and compare them to field-
theory predictions. Finally, Sec. V contains a summary
and a discussion of our results.
II. EFFECTIVE FIELD THEORY
We start from an effective field theory describing the
long-wavelength fluctuations of Eq. (1). In the limit of
strong next-nearest neighbor interactions J2 ≫ |J1|, the
spin operators can be expressed as:
Szα(r) ∼ m+ c(m) sin
2kF r +
+ · · ·
S−α (r) ∼ (−1)re−iθα
π + · · · . (2)
kF = (
−m)π is the Fermi-wave vector and α = 1, 2 enu-
merates the two chains of the zig-zag ladder. In relation
with Eq. (1), note that ~S1(r) = ~S(x+1)/2 (~S2(r) = ~Sx/2)
for x odd (even). φα and θα are compactified quantum
fields describing the out-of-plane and in-plane angles of
fluctuating spins obeying Gaussian Hamiltonians:
H = v
(∂xφα)
2 +K(∂xθα)
, (3)
with [φα(x), θα(y)] = iΘ(y − x), where Θ(x) is the
Heaviside function. Sub-leading terms are suppressed in
Eq. (2). m is the magnetization of decoupled chains, re-
lated to the real magnetization M of the zig-zag system
M ≃ m
1− 2K(m)J1
πv(m)
. (4)
K(m) and v(m) are the Luttinger liquid (LL) parame-
ter and the spin-wave velocity of the decoupled chains,
http://arxiv.org/abs/0704.0764v2
respectively. The nonuniversal amplitude c(m) appear-
ing in the bosonization formulas (2) has been determined
from density matrix renormalization group (DMRG)
calculations.12 Note that in our notation M = 1/2 at
saturation.
Now we perturbatively add the interchain coupling
term to two decoupled chains, each of which is described
by an effective Hamiltonian of the form Eq. (3) and fields
φi and θi, i = 1, 2. For convenience, we transform to
the symmetric and antisymmetric combinations of the
bosonic fields φ± = (φ1±φ2)/
2 and θ± = (θ1±θ2)/
In this basis and apart from terms H±0 of the form (3),
the effective Hamiltonian describing low-energy proper-
ties of Eq. (1) contains a single relevant interaction term
with the bare coupling g1 ∝ J1 ≪ v:
Heff = H+0 +H
0 + g1
dx cos
, (5)
and the renormalized LL parametersK± are, in the weak
coupling limit:
K± = K
1∓ J1
. (6)
K+ is the Luttinger-liquid parameter of the soft mode
of the zig-zag ladder. The Hamiltonian (5) represents
the minimal effective low-energy field theory describing
the region J2 ≫ |J1| of the frustrated FM spin-1/2 chain
for M 6= 0.9,13 The relevant interaction term cos
opens a gap in the φ− sector. Since S
x+1 − Szx ∼ ∂xφ−,
relative fluctuations of the two chains are locked. This
implies that single-spin flips are gapped with a sine-
Gordon gap in the sector describing relative spin fluctua-
tions of the two-chain system.9 Gapless excitations come
from the ∆Sz = 2 channel, i.e. only those excitations
are soft where spins simultaneously flip on both chains.
DMRG results show that this picture applies to a large
part of the magnetic phase diagram.9
III. CORRELATION FUNCTIONS
We now turn to the ground state properties of Eq. (1)
as a function of magnetization, concentrating on several
correlation functions in order to identify the leading in-
stabilities. Note that our analysis is only valid if M 6= 0.
Apart from a term representing the magnetization M in-
duced by the external field, the longitudinal correlation
function shows an algebraic decay with distance r:
〈Szα(0)Szβ(r)〉 ≃ M2+
C1 cos(2kF r + (α− β)kF )
2π2rK+
8π2r2
The constants Ci, i = 1, 2, 3, appearing here and in
Eq. (9) will be determined through a comparison with
numerical results.
In contrast to Eq. (7), the transverse xy-correlation
functions decay exponentially reflecting the gapped na-
ture of the single spin-flip excitations. Here we do not
restrict ourselves to the equal-time expression only, be-
cause we will need non-equal time correlation functions
to extract the finite-size corrections to the gap later on.
We obtain:
〈S+α (0, 0)S−β (r, τ)〉 ≃
δα,β(−1)re−∆1(M)
τ2+r2/v2
(r2 + v2+τ
8K+ (r2 + v2−τ
where τ stands for the Euclidean time, ∆1(M) is the
∆Sz = 1 gap, and v± ∼ v ± J1/π in the weak cou-
pling limit. The Kronecker delta strictly applies to the
thermodynamic limit, while on the lattice an additional
contribution for α 6= β exists.
It is noteworthy that, different from Eq. (8), the in-
plane correlation functions involving bilinear spin combi-
nations decay algebraically. This stems from the gapless
nature of ∆Sz = 2 excitations. In fact, these are the
slowest decaying correlators close to the saturation mag-
netization:
〈S+1 (r)S
2 (r)S
1 (0)S
2 (0)〉 ≃
r1/K+
C3 cos(2kF r)
rK++1/K+
This result is reminiscent of a partially ordered state be-
cause the ordering tendencies in this correlation function
are more pronounced than those of the corresponding
single-spin correlation function Eq. (8). Therefore, we
call the correlator (9) ‘nematic’. Furthermore, we will
refer to a situation where Eq. (9) is the slowest decay-
ing one among all correlation functions as a ‘nematic-like
phase’.
By virtue of the exponential decay in (8), the correlator
(9) is proportional to:
〈(S+1 (r) + S
2 (r))
2 (S−1 (0) + S
2 (0))
2〉 . (10)
The term (Sα1 + S
2 appearing in the case of the S = 1
zig-zag ladder corresponds to the operator (Sα)2 in the
case of a S = 1 chain. One can think of an effective S = 1
spin formed from two neighboring S = 1
spins coupled
by the ferromagnetic interaction. A similar behavior of
correlation functions, namely the exponential decay of in-
plane spin components and the algebraic decay of their
bilinear combinations, is encountered also in the XY 2
phase of the anisotropic S = 1 chain14 and in the spin-1
chain with biquadratic interactions, see, e.g., Ref. 15.
The algebraic decay of the nematic correlator as op-
posed to the exponential decay of (8) suggests that there
are tendencies towards nematic ordering in this phase.
Depending on the value of K+ the dominant instabilities
are either spin-density-wave ones for K+ < 1 or nematic
ones for K+ > 1. From the result for K+ given in Eq. (6)
one can perturbatively evaluate the crossover value of J1:
|J1,cr| =
πv(m)
. (11)
For J1 < J1,cr the nematic correlator (9) is the slowest
decaying one, i.e. one is in the nematic-like phase. The
-0.015
-0.01
-0.005
0.005
0.001
0 4 8 12 16 20 24 28 32
fit, L=32
fit, L=48
fit, L=64
ED, L=32
ED, L=48
ED, L=64
FIG. 1: (Color online) Correlation functions at J1 = −J2 < 0,
and magnetization M = 3/8: (a) longitudinal component Sz
(b) transverse component S±
, (c) spin nematic S±
. x is
the distance in a single-chain notation. ED results for periodic
boundary conditions are shown by symbols, fits by lines. Note
the logarithmic scale of the vertical axis in panel (b).
behavior of the cross-over line can be read off from the
behavior of K(m): K(m) increases monotonically with
m, tends to K = 1 for m → 1/2, and satisfies K < 1
for m < 1/2 (see, e.g., Refs. 16,17). Therefore, we have
J1,cr = 0 for M = 1/2 with increasing ferromagnetic
|J1,cr| for decreasing M . This means that for J1 < 0
a regime opens at high M where nematic correlations
given by Eq. (9) dominate over spin-density-wave corre-
lations given by Eq. (7), in agreement with Chubukov’s
prediction.8
Now we check the correlation functions obtained within
bosonization against exact diagonalization (ED) results.
Numerical data obtained for J1 = −J2 < 0 and M = 3/8
on finite systems with periodic boundary conditions are
shown in Fig. 1. This parameter set allows for a clear
test of the above predictions, but represents the generic
behavior in the phase of two weakly coupled chains. To
take into account finite-size effects we use the observation
that for a conformally invariant theory, any power law
on a plane becomes a power law in the following variable
defined on a cylinder of circumference L:
x → L
. (12)
First we fit the nematic correlator given by Eq. (9),
which from bosonization is expected to be the leading
instability at high magnetizations. Using the part with
x ≥ 5 of the L = 64 data shown in Fig. 1c, we find
1/K+ = 0.904 ± 0.011, C2 = 0.143 ± 0.004, and C3 =
−0.326±0.013. Fig. 1c shows that all finite-size results for
the nematic correlator are nicely described by this fit with
the dependence on L taken into account by substituting
Eq. (12) for the power laws. Moreover, from K+ > 1, we
see that the system is indeed in the region dominated by
nematic correlations for M = 3/8 and J1 = −J2.
Now we turn to the longitudinal correlation function
which we fit to the bosonization result Eq. (7). Since
most numerical parameters have been determined by the
previous fit, only one free parameter is left which we de-
termine from the numerical results of Fig. 1a for L = 64
and x ≥ 14 as C1 = 0.060± 0.004. Predictions for other
system sizes are again obtained by substituting Eq. (12)
for the power laws. The agreement in Fig. 1a is not as
good as in Fig. 1c. However, it improves at larger dis-
tances x and system sizes L, indicating that corrections
omitted in Eq. (7) are still relevant on the length scales
considered here.
Finally, the xy-correlation function is shown in Fig. 1b
with a logarithmic scale of the vertical axis of this panel.
The exponential decay predicted by Eq. (8) is verified.
One further observes that correlations between the in-
plane spin-operators belonging to different chains (odd
x) are an order of magnitude smaller than on the same
chain (even x). This suppression of correlations between
different chains corresponds to the δ symbol in (8), which
strictly applies only in the thermodynamic limit and for
large distances.
We summarize the main result of this section: in-plane
spin correlators are exponentially suppressed for any fi-
nite value of the magnetization in the parameter region
|J1| < J2. The ground state crosses over from a spin-
density-wave dominated to a nematic-like phase with in-
creasing magnetic field, with the crossover line given by
Eq. (11).
IV. EXCITATIONS
We next address the excitation spectrum. Since the
gap to ∆Sz = 1 excitations should be directly accessible
to microscopic experimental probes such as inelastic neu-
tron scattering or nuclear magnetic resonance, we analyze
its behavior as a function of magnetization. Sufficiently
below the fully polarized state the gap can be calculated
analytically using results from sine-Gordon theory. In
addition, to leading order of the interchain coupling, one
can get qualitative expressions using dimensional argu-
ments for the perturbed conformally invariant model:
∆1(m) ∼
c2(m)|J1| sin(πm)
v(m)(1 − J1K(m)/πv(m))
, (13)
where ν = 2− 2K(m)
1+ J1K(m)/πv(m)
. m(h),K(h)
and v(h) can be determined numerically from the Bethe
ansatz integral equations.16,17,18,19
With this information and Eqs. (4) and (13) we de-
termine the qualitative behavior of the single-spin gap
∆1(M) as a function of M : it increases from zero at
zero magnetization, reaches a maximum at intermediate
magnetization values, then shows a minimum and, upon
approaching the saturation magnetization, it increases
again. As our formulas do not strictly apply at m = 0,
the notion of a vanishing gap at zero magnetization may
be a spurious result. Note that when the fully polarized
state is approached, the magnetization increases in an un-
physical fashion since in this limit bosonization becomes
inapplicable. At the point where the magnetization sat-
urates the exact value of the gap can be obtained from
the following mapping to hard-core bosons:8,11,20
Szi =
− a†iai , S
i = a
i . (14)
Comparing Eq. (14) with Eq. (2) one recognizes the lead-
ing terms in Haldane’s harmonic fluid transformation for
bosons.20 Using a ladder approximation which is exact in
the two-magnon subspace we arrive at:
4J22 − 2J1J2 − J21
2(J2 − J1)
J21 + 8J1J2 + 16J
J21 (3 J2 + J1)
J2 (J2 − J1)
. (15)
In Eq. (15) we have represented the gap as a difference
of two terms: the quantum and the classical instability
fields emphasizing its quantum origin.
In order to verify these field theory predictions, we
perform complementary numerical computations using
the DMRG method.21 Open boundary conditions are im-
posed and we typically keep up to 400 DMRG states.
From DMRG we obtain the ground-state energies E(Sz)
as a function of total Sz. For those values of Sz that
emerge as a ground state in an external magnetic field
we compute the single-spin excitation gap from
∆1(M) =
E(Sz + 1) + E(Sz − 1)− 2E(Sz)
. (16)
Fig. 2a shows numerical results for ∆1 at a selected value
of J1 = −0.3 J2 < 0 for the largest system sizes investi-
gated. We find that the finite-size behavior of the gap
∆1(M,L) for system sizes L ≥ 24 is well described by
a 1/L correction. This will be further corroborated by
field-theoretical arguments outlined below. Therefore, we
0 0.1 0.2 0.3 0.4 0.5
L=120
L=144
L=156
L=120
L=144
L=156
extrapolated
FIG. 2: (Color online) Density matrix renormalization group
results for the gaps at J1 = −0.3 J2 < 0 as a function of
magnetization M . Panel (a) shows the single-spin excitation
gap (16), panel (b) the finite-size gap (19) for two flipped
spins multiplied by the chain length L.
extrapolate it to the thermodynamic limit using a fit to
the form
∆1(M,L) = ∆1(M) +
+ · · · , (17)
allowing for an additional 1/L2 correction for those values
of M where at least 4 different system sizes are available.
This extrapolation is represented by the full circles in
Fig. 2a; errors are estimated not to exceed the size of the
symbols. Our extrapolation for ∆1 is consistent with a
vanishing gap at M = 0 in agreement with previous nu-
merical studies22 although bosonization predicts a non-
zero – possibly very small – gap.13,22,23 The behavior of
∆1(M) confirms the picture described above: the gap is
non-zero for M > 0, goes first through a maximum and
then a minimum and finally approaches ∆1/J2 ≈ 0.023
given by Eq. (15) for M → 1/2.
We further wish to point out that for chains with pe-
riodic boundary conditions, the coefficient a(M) of the
finite-size extrapolation Eq. (17) is determined by the
spin-wave velocity and the critical exponent of the soft
mode from the ∆Sz = 2 channel. Indeed, using Eq. (8)
where we can set r = 0, and use the conformal mapping
(12) to the cylinder, we see that the leading finite-size
correction to the gap is:
∆1(M,L) = ∆1(M) +
πv+(M)
4K+(M)
. (18)
0 π/2 π
wave vector
1.0 0.0
FIG. 3: (Color online) Numerical dispersion spectrum in the
subspaces of odd Sz computed for L = 24 and J1 = −J2 < 0.
The wave vector is given relative to the ground state wave
vector (0 for Sz = 0, 4, 8 and π for Sz = 2, 6).
Note that we have to replace sin with sinh in Eq. (12)
in order to extract a gap, since we are dealing with Eu-
clidean time. In addition we used the fact that in our ap-
proximation the effective Hamiltonian (5) is a direct sum
of symmetric and antisymmetric sectors. Moreover, it is
only the symmetric sector enjoying conformal invariance
and consequently we perform the replacement τ → sinh τ
only in the symmetric sector. The antisymmetric sector
has a spectral gap and its contribution to the finite-size
corrections of the single-spin flip excitation energy are
exponentially suppressed with system size.24 With this
method one cannot fix the amplitudes of the 1/L2 term
and beyond. Note furthermore that there may be ad-
ditional surface terms for open boundary conditions as
employed in the numerical DMRG computations. Nev-
ertheless there is a dominant 1/L correction in any case.
Next, we briefly look at the ∆Sz = 2 excitations. Their
finite-size gap is, in analogy to Eq. (16), computed with
DMRG from
∆2(M) =
E(Sz + 2) + E(Sz − 2)− 2E(Sz)
. (19)
Fig. 2b shows numerical results for L∆2(M,L) again at
the value J1 = −0.3 J2 < 0. One observes that the scaled
finite-size gaps collapse onto a single curve which shows
that ∆2(M,L) scales linearly to zero with 1/L, exactly
as expected for gapless excitations in 1D. Furthermore,
we observe that the scaled quantity L∆2(M,L) vanishes
as one approaches saturation M = 1/2 which indicates a
vanishing of the velocity of the corresponding excitations
at saturation.
We proceed by discussing the wave-vector dependence
of the ∆Sz = 1 excitation, while we remind the reader
that the low-energy excitations are in the ∆Sz = 2 sec-
tor. Fig. 3 shows representative ED results obtained for
rings with L = 24 and J1 = −J2 < 0. For ground
states with low Sz, the ∆Sz = 1 excitation spectrum
looks similar to the continuum of spinons. On the other
hand, close to saturation one has single-magnon excita-
tions with a minimum given by the classical value of
the wave vector kcl = arccos(|J1|/4J2).8,13,25 We read
off from Fig. 3 that upon lowering the magnetic field,
this minimum shifts from the classical incommensurate
value towards π/2, i.e. the value appropriate for two de-
coupled chains. This renormalization of the minimum of
the magnon excitations towards the value of decoupled
chains can be interpreted in terms of quantum fluctua-
tions, which are enhanced when the density of magnons
increases. A strong quantum renormalization of the pitch
angle from its classical value at zero magnetization was
previously observed by the coupled-cluster method and
DMRG calculations.26
V. SUMMARY
We have combined numerical techniques with analyti-
cal approaches and mapped out the ground state phase
diagram of the frustrated ferromagnetic spin chain in an
external magnetic field. We have established that with
increasing magnetic field, the ground state crosses over
from a spin-density-wave dominated to a nematic-like
phase. Single spin flip excitations are gapped, giving rise
to an exponential decay of in-plane spin correlation func-
tions in both regimes. We have studied the single- and
two-spin flip excitation energy numerically. Using tools
from conformal field theory we have further shown that
the amplitude of the leading 1/L correction term to the
single-spin flip gap is determined by the critical exponent
and the spin-wave velocity of the soft mode.
Finally, in order to apply our findings to the mate-
rial LiCuVO4, one should take into account interchain
interactions as well as anisotropies, which are expected
to be present in this system.3 At low fields, a helical state
has been observed experimentally.2,3 On the other hand,
for the purely one dimensional case, we have shown that
upon increasing the magnetic field there is a competition
between spin-density-wave and nematic-like tendencies.
Those are the two leading instabilities at high magnetiza-
tions and thus they are the natural candidates to become
long-range ordered in higher dimensions. The question
whether there are true phase transitions at high fields
in higher dimensions is beyond the scope of the current
work.
Acknowledgments
We thank A. Feiguin for providing us with his DMRG
code used for large scale calculations. Most of T.V.’s
work was done during his visits to the Institutes of
Theoretical Physics at the Universities of Hannover and
Göttingen, supported by the Deutsche Forschungsge-
meinschaft. The hospitality of the host institutions is
gratefully acknowledged. T.V. also acknowledges support
from the Georgian National Science Foundation under
grant N 06−81−4−100. LPTMS is a mixed research unit
8626 of CNRS and University Paris-Sud. A.H. is sup-
ported by the Deutsche Forschungsgemeinschaft (Project
No. HO 2325/4-1), and F.H.-M. is supported by NSF
grant No. DMR-0443144.
1 M. Hase, H. Kuroe, K. Ozawa, O. Suzuki, H. Kitazawa, G.
Kido, and T Sekine, Phys. Rev. B 70, 104426 (2004).
2 B. J. Gibson, R. K. Kremer, A. V. Prokofiev, W. Assmus
and G. J. McIntyre, Physica B 350, e253 (2004).
3 M. Enderle, C. Mukherjee, B. F̊ak, R. K. Kremer, J.-M.
Broto, H. Rosner, S.-L. Drechsler, J. Richter, J. Malek,
A. Prokofiev, W. Assmus, S. Pujol, J.-L. Raggazzoni, H.
Rakoto, M. Rheinstädter, and H. M. Rønnow, Europhys.
Lett. 70, 237 (2005).
4 M. G. Banks, F. Heidrich-Meisner, A. Honecker, H.
Rakoto, J.-M. Broto, and R. K. Kremer, J. Phys.: Cond.
Mat. 19, 145227 (2007).
5 N. Büttgen, H.-A. Krug von Nidda, L. E. Svistov, L. A.
Prozorova, A. Prokofiev, and W. Aßmus, Phys. Rev. B 76,
014440 (2007).
6 S.-L. Drechsler, O. Volkova, A. N. Vasiliev, N. Tristan, J.
Richter, M. Schmitt, H. Rosner, J. Málek, R. Klingeler, A.
A. Zvyagin, and B. Büchner, Phys. Rev. Lett. 98, 077202
(2007).
7 H.-J. Mikeska and A. K. Kolezhuk, Lect. Notes Phys. 645,
1 (2004).
8 A. V. Chubukov, Phys. Rev. B 44, 4693 (1991).
9 F. Heidrich-Meisner, A. Honecker, and T. Vekua, Phys.
Rev. B 74, 020403(R) (2006).
10 D. V. Dmitriev, V. Y. Krivnov, and J. Richter, Phys. Rev.
B 75, 0114424 (2007).
11 R. O. Kuzian and S.-L. Drechsler, Phys. Rev. B 75, 024401
(2007).
12 See, e.g., T. Hikihara and A. Furusaki, Phys. Rev. B 69,
064427 (2004).
13 D. C. Cabra, A. Honecker, and P. Pujol, Eur. Phys. J. B
13, 55 (2000).
14 H. J. Schulz, Phys. Rev. B 34, 6372 (1986).
15 A. Läuchli, G. Schmid, and S. Trebst, Phys. Rev. B 74,
144426 (2006).
16 K. Totsuka, Phys. Lett. A 228, 103 (1997).
17 D. C. Cabra, A. Honecker, and P. Pujol, Phys. Rev. B 58,
6241 (1998).
18 V. E. Korepin, N. M. Bogoliubov, and A. G. Izergin, Quan-
tum Inverse Scattering Method and Correlation Functions
(Cambridge University Press, Cambridge, England, 1993).
19 S. Qin, M. Fabrizio, L. Yu, M. Oshikawa, and I. Affleck,
Phys. Rev. B 56, 9766 (1997).
20 F. D. M. Haldane, Phys. Rev. Lett 47, 1840 (1981).
21 S. R. White, Phys. Rev. Lett. 69, 2863 (1992); Phys. Rev.
B 48, 10345 (1993).
22 C. Itoi and S. Qin, Phys. Rev. B 63, 224423 (2001).
23 A. A. Nersesyan, A. O. Gogolin, and F. H. L. Eßler, Phys.
Rev. Lett. 81, 910 (1998).
24 S. I. Matveenko and A. M. Tsvelik, private communication.
25 C. Gerhardt, K.-H. Mütter, and H. Kröger, Phys. Rev. B
57, 11504 (1998).
26 R. Bursill, G. A. Gehring, D. J. J. Farnell, J. B. Parkinson,
T. Xiang, and C. Zeng, J. Phys. Cond. Mat. 7, 8605 (1995).
