Field-controlled spin current in frustrated spin chains

We study states with spontaneous spin current, emerging in frustrated antiferromagnetic spin-$S$ chains subject to a strong external magnetic field. As a numerical tool, we use a non-Abelian symmetry realization of the density matrix renormalization group. The field dependence of the order parameter and the critical exponents are presented for zigzag chains with S=1/2, 1, 3/2, and 2.


Introduction
In usual magnets, spin degrees of freedom order ferro-or antiferromagnetically, while those two basic orders may be viewed as special cases of a general spin density wave order parameter. It had been known for a long time [ 1,2] that other "exotic" spin orderings with a vanishing expectation value of the spin density are theoretically possible, but only recently such orders have been observed in numerical studies of realistic models [ 3,4,5,6,7]. One prominent example of unconventional order parameter is the spin current (the so-called vector chirality [ 1], or p-type nematic [ 2] order parameter). In a spin chain, the vector chirality operator κ(n) ≡ κ n,n+1 can be defined as a vector product of two adjacent spins: where n labels the lattice sites along the chain. Obviously, the vector chiral order is nonzero in a state which has a classical helical magnetic order, and the projection of κ onto the helix axis distinguishes between the left and right spirals. For low-dimensional magnets, such unconventional order parameters may become dominating, since quantum fluctuations have a tendency to destroy the usual magnetic order. A state with vector chiral order may be viewed as the result of phase fluctuations destroying the helical order, but still keeping the preferred sense of rotation contained in a classical spiral [ 8]. The interpretation of κ as the spin current can be understood from the equation of motion, which for a Hamiltonian of the form H = nm J m S n · S n+m readsh(∂ S n /∂t) = m J m { κ n−m,n − κ n,n+m }. One can show that in the ground state the net spin current should be zero, which for the above Hamiltonian translates into the condition m mJ m κ n,n+m = 0.
Spin current states have been relatively well studied for anisotropic frustrated chains [ 9,3,10,4,11]. recently, it has been shown [ 12,6,13,7] that a strong external magnetic field can be used as a control parameter to drive the spin current in isotropic frustrated spin chains. The interest to this topic is further boosted by the fact that there are currently several quasi-onedimensional magnetic materials that are intensely studied as possible candidates for manifestation of unconventional orders, particularly LiCuVO 4 [ 14], Li 2 ZrCuO 4 [ 15], and Cu 2 Cl 4 -H 8 C 4 SO 2 (Sul-Cu 2 Cl 4 ) [ 16]. We have previously studied [ 6] the spin current correlations in frustrated chains with the spin S = 1/2 and 1, and have shown that the behavior of S = 1/2 and S = 1 chains was very different. It was not clear whether this can be attributed to some general difference in the behavior of chains with integer and half-integer S. The goal of the present paper is to answer this question. For that purpose, we numerically study the spin current in isotropic spin-S zigzag chains with the spin up to S = 2. It is shown that chains with S ≥ 1 behave similarly and exhibit a state with a finite spontaneous spin current in the entire region of finite magnetization. In Sect. 2 we give a brief overview of the theory of spin current correlations, Sect. 3 presents the numerical results, and Sect. 4 contains a short summary.

Frustrated chain in magnetic field: spin current correlations
We consider a frustrated antiferromagnetic spin chain described by the following Hamiltonian: where S n are spin-S operators at the n-th site, J 1 > 0 and J 2 > 0 are the nearest and next-nearest neighbor exchange constants, respectively, and H is the external magnetic field, assumed to be applied along the z axis. At zero field, if the frustration parameter α ≡ J 2 /J 1 is large enough, the ground state is generally expected to have a finite spectral gap ∆, both for integer and halfinteger S. When the applied field exceeds the critical value H c = ∆, the system acquires finite magnetization. Further increase of the field beyond the saturation field value H s brings the system into a fully polarized state.
Let us first assume that α ≫ 1, so that one deals with the limit of two weakly coupled spin-S subchains. The behavior of spin chains in magnetic field have been extensively studied for S = 1/2 [ 17,18,19] and S = 1 [ 20,21,22]. We will assume that for H c < H < H s the low-energy physics of a single chain is well described in terms of the effective Tomonaga-Luttinger liquid (TLL) theory with the following Hamiltonian: Here K is the so-called TLL parameter, v ∝ J 2 is the Fermi velocity, φ is the compact bosonic field (ϕ ≡ ϕ + √ π), and θ is its dual satisfying the commutation relations [ϕ(x), θ(y)] = iΘ(y − x), where Θ(x) is the Heaviside function.
We assume further, by the analogy with the S = 1/2 and S = 1 case [ 17,20,22], that in the continuum limit the most relevant low-energy part of spin operators can be represented as Here a = 1, 2 labels the two subchains, the space coordinate x is defined at the middle of every bond along the original zigzag chain, x a = x+ a− 3/2 define the sites of the subchains, M S is the ground state magnetization per spin, k F has a meaning of the Fermi momentum and is connected to M (e.g., k F = π(1+M )/2 for S = 1/2), A i are nonuniversal amplitudes. All parameters (K, v, M , A i ) should be understood as functions of the field H which generally have to be extracted by comparing the numerical results with those following from the TLL description [ 17,18,19,20,21,23]. Introducing symmetric and antisymmetric combinations of the bosonic fields ϕ ± = (ϕ 1 ± ϕ 2 )/ √ 2, θ ± = (θ 1 ±θ 2 )/ √ 2, one obtains two coupled TLLs, and the nature of the leading interaction term strongly depends on the value of the TLL parameter K. The longitudinal (S z S z ) part of the zigzag exchange leads to a splitting of the TLL parameter values for the symmetric and antisymmetric sectors, which for α ≫ 1 are found as 2K/(πvα)) 1/2 . At the same time, the S z S z part of the interaction gives rise to the interaction term proportional to cos √ 8πϕ − − πM whose scaling dimension is 2K − and which favors n-type spin-nematic (or XY2 in the nomenclature of Ref. [ 24]) correlations [ 25]. The transversal part of the zigzag exchange yields the so-called "twist term" sin √ 2πθ − (∂ x θ + ) that has a nonzero conformal spin and the formal scaling dimension 1 + 1/(2K − ). The twist term favors states with spontaneous spin current [ 9]. If K < 1, as is the case for S = 1/2, the two interaction terms compete, which leads to the Ising-type transition between the nematic and spin-current state [ 12]. For S = 1 it is known that K > 1 [ 20,21,23], thus the nematic-favoring term is irrelevant and the spin-current favoring term dominates. For higher spins, to our knowledge, there is no numerical data on the behavior of the TLL parameter K. However, since the analysis based on a mapping to the nonlinear σ-model [ 20] can be extended to higher S, one may expect that for all spin-S chains with S ≥ 1, the TLL parameter K is greater than 1, at least in the region of small (but finite) magnetization M .
Inside the phase with a spontaneous spin current, the twist term is the most relevant perturbation, and the mean-field treatment along the lines of Ref. [ 9] leads to the conclusion [ 12] that both ∂ x θ + and sin √ 2πθ − should be nonzero, and the antisymmetric sector should be massive. Expressing the spin current κ(n) as defined in (1) through the bosonic fields, one obtains [ 6] where we have omitted the contributions from massive fields and operators with higher scaling dimensions. At large distances x ≫ ξ, where ξ is the largest correlation length determined by the smallest gap in the neglected massive fields, the asymptotics of the leading spin-current correlations takes the form where the incommensurate wave vector Q is given by Q = π/2 + π/2 ∂ x θ + and C i are some numbers. The asymptotics (6), derived by means of bosonization, is expected to be valid in the limit α ≫ 1 and for H sufficiently far from any of the critical fields H c , H s (close to H c or H s the bosonization approach becomes inapplicable since the effective bandwidth goes to zero and the interactions cease to be small compared to the bandwidth). In a close vicinity of the saturation field H s , large-S analysis [ 12] maps the system to a two-component dilute Bose gas with repulsive interaction, and the external magnetic field drives the system into a condensed state, playing the role of the chemical potential. The interspecies repulsion turns out to be strong enough to satisfy the phase separation condition, so only one of the species condenses and the other condensate is depleted. At the end, one deals with the one-component pseudo-condensate whose physics is again described by a (one-component) TLL, and the asymptotic form of the spin current correlators for H close to H s is given by [ 12]: where , K ′ > 1 is another TLL parameter (different from K considered above), and the wave vector Q ′ is just the pitch of the classical helix, Q ′ = ±(π − arccos(1/4α)).

Results of numerical study
We have studied numerically the frustrated chain model (2) for systems with the spin S = 1/2, 1, 3/2, and 2 and the frustration parameter α = J 2 /J 1 was fixed at α = 1 (so the chains may be considered as stripes of the triangular lattice). We have used the density matrix renormalization group (DMRG) method [ 26,27] in its matrix product state formulation [ 28], making full use of the non-Abelian SU (2) symmetry [ 29,28]. Although the external field breaks SU (2)  fact that the Zeeman energy term commutes with the rest of the Hamiltonian makes it possible to include the effect of the magnetic field by calculating the ground state in a sector with a finite total spin S tot . The use of SU (2) symmetry allows one to reduce substantially the number of states m which is necessary to describe the system: essentially, one treats the multiplet of states of the same total spin as a single representative state. A slight disadvantage is that the non-Abelian method allows to compute only reduced matrix elements (in the sense of the Wigner-Eckart theorem). Since the spin current is a vector, the available correlator that is easily available is the rotationally invariant scalar product κ(n) · κ(n ′ ) , which is a mixture of the longitudinal and transversal spin current correlations. This complicates somewhat the analysis of the data: from the theoretical analysis in Sect. 2 it follows that for typical values of K ∼ 1 the longitudinal correlations κ z (x)κ z (0) decay to their asymptotic value much faster than the transversal ones. Thus by analyzing the scalar correlator κ(n) · κ(n ′ ) one can only estimate the critical exponent η = η ⊥ characterizing the decay of the transversal spin current correlations. We have studied zigzag chains of several lengths L ranging from 64 to 256. It turned out that a relatively large number m of representative states was necessary to reach good convergence; depending on S, L and S tot , we have used m ranging from 400 to 1500. The scalar spin current correlator κ(n) · κ(n ′ ) has been calculated for a large number of ground states in sectors with different total spin S tot .
The correlator has been averaged over the starting and final positions n, n ′ , and contributions with n or n ′ being closer as as a fixed "cutoff" (taken here to be 20 sites) to the chain ends were discarded. The DMRG data for the correlator has been fitted to the power-law form suggested by (6), (7); the introduction of a finite phase shift δ helps to suit the open boundary conditions. From those fits we have extracted the behavior of the equilibrium spin current κ 0 and the exponent η as functions of the chain magnetization M = S tot /L, shown respectively in Fig. 1 and Fig. 2. The wave vector q extracted from the fit depends only weakly on the magnetization, for all S studied, so we do not present the corresponding dependences. The quality of fits is generally at its best in the intermediate M region, and is deteriorating at M → 0 or M → 1 for the reasons discussed in Ref. [ 6]. For small M it was only possible to extract the asymptotic value κ 2 0 of the spin current correlator, but not the critical exponent η.

Summary and discussion
We have studied spin-S isotropic antiferromagnetic frustrated chains in strong magnetic fields, described by the model (2) with the frustration parameter α = 1 and the spin S ranging from 1/2 to 2, by means of the matrix product density matrix renormalization group technique. Existence of a phase with field-induced spontaneous spin current (vector chiral, or p-type nematic order) is established for all S studied. For S ≥ 1, the behavior of the order parameter and its correlations as functions of the magnetization qualitatively agrees with the theory proposed in Ref. [ 12]. For S = 1/2, the existence of the vector chiral phase and the transition between the chiral and XY2 phases are captured well by the large-α bosonization approach of Ref. [ 12].
However, the behavior of S = 1/2 chain close to the saturation M = 1 does not fit into the common large-S scheme of [ 12]: there is a finite non-chiral region in the immediate vicinity of M = 1 which, as shown in [ 13], belongs to the two-component Tomonaga-Luttinger liquid (TLL) phase. This peculiarity might be attributed to the fact that the large-S approach in [ 12] is based on identifying the two emerging particle species as bosons, while the underlying particles for S = 1/2 are the Jordan-Wigner fermions. From the Bethe-ansatz results for the Hubbard model [ 30] it is known that pure density-density interaction does not destroy the two-component TLL in the low particle density limit. Thus, for S = 1/2 the nature of the transition from the two-component TLL into the chiral phase close to the saturation field still remains to be explained.