---
abstract: 'We show that spin-$S$ chains with SU(2)-symmetric, ferromagnetic nearest-neighbor and frustrating antiferromagnetic next-nearest-neighbor exchange interactions exhibit metamagnetic behavior under the influence of an external magnetic field for small $S$, in the form of a first-order transition to the fully polarized state. The corresponding magnetization jump increases gradually starting from an $S$-dependent critical value of exchange couplings and takes a maximum in the vicinity of a ferromagnetic Lifshitz point. The metamagnetism results from resonances in the dilute magnon gas caused by an interplay between quantum fluctuations and frustration.'
author:
- 'M. Arlego'
- 'F. Heidrich-Meisner'
- 'A. Honecker'
- 'G. Rossini'
- 'T. Vekua'
date: 'May 27, 2011; revised November 22, 2011'
title: Resonances in a dilute gas of magnons and metamagnetism of isotropic frustrated ferromagnetic spin chains
---

Introduction {#sec:intro}
============

Quantum spin systems at low temperatures share many of the macroscopic quantum behaviors with systems of bosons such as Bose-Einstein condensation, superfluidity (see Ref.  and references therein), or macroscopic quantum tunneling, [@chudnovsky88] and may indeed be viewed as quantum simulators of interacting bosons.[@giamarchi08] Spin systems with frustration, in particular, realize exotic phases of strongly correlated bosons, such as spin liquids [@balents10] or supersolids.[@liu73] Experimentally, such systems can be studied both in low-dimensional quantum magnets (see, e.g., Ref. ) and in ultra-cold atomic gases.[@struck11; @jo11] In the latter case, interest in the many-body physics has been excited by the extraordinary control over interactions between bosons via Feshbach resonances. [@bloch08] These can, in particular, be used to tune interactions from repulsive to attractive. In the attractive regime, the collapse of an ultra-cold gas of bosons was observed in experiments.[@gerton00] In this work, we discuss a mechanism by which the same can be achieved in spin systems.

In the large-$S$ limit, spins map onto bosons with a finite but large Hilbert space, justifying a description in terms of soft-core bosons. Yet, in many cases, similarities between spins and soft-core bosons may even exist for $S=1/2$. In three dimensions (3D), systems of spins and bosons resemble each other the more the smaller the density of bosons is, whereas 1D spin-$S$ antiferromagnets close to saturation behave as spinless fermions, or hard-core bosons. In spin systems, the external magnetic field $h$ tunes the density of magnons. The limit of a dilute gas of magnons is then realized as the fully polarized state (the vacuum of magnons) is approached from below.

In this work, we argue that resonances can play a crucial role in frustrated quantum spin systems largely determining their low-energy behavior in an external magnetic field. We consider the frustrated ferromagnetic (FM) spin-$S$ chain and show that upon changing system parameters such as the coupling constants or $h$, one can tune the effective interaction between magnons from repulsive to attractive by exploiting the existence of resonances. As a main result, we demonstrate that, close to resonances, where the scattering length is much larger than the lattice constant, and in the case of attractive effective interactions, the intrinsic hard-coreness of spins does not play a significant role in the limit of a dilute gas of magnons. A behavior resembling the collapse of attractively interacting bosons [@nozieres82; @Mueller] therefore exists in such spin systems close to their fully polarized state. The thermodynamic instability of collapsed states causes jumps in the magnetization curve just below saturation. This has to be contrasted with the scattering length being of the order of a few lattice sites. In that case, which is realized for spin 1/2, the formation of mutually repulsive, multi-magnon bound states can be observed (see, e.g., Refs. ).

Our work demonstrates that the mapping of a [*purely 1D*]{} spin system close to saturation to an effective theory of a dilute Bose gas properly accounts for the physics of the model. This is accomplished by connecting the scattering vertices of the microscopic lattice model with the coupling constants of the effective theory. Furthermore, despite the purely 1D nature of our problem, a $1/S$ expansion is a valuable tool and yields the correct physics.

Concretely, we study the following system: $$H_S = \sum_{i=1}^L \left \lbrack J\vec{S}_i\cdot \vec{S}_{i+1}+ J'   \vec{S}_i\cdot \vec{S}_{i+2}-h  S^z_i \right\rbrack\,.
\label{eq:ham}$$ $\vec{S}_i=(S_i^x, S^y_i, S^z_i)$ is a spin-$S$ operator acting on site $i$ and $L$ is the number of sites. $J < 0$ is the FM, nearest-neighbor exchange interaction and $J'=1$ is the antiferromagnetic (AFM), next-nearest-neighbor exchange interaction setting the energy scale. In the absence of an external field $h$, $H_S$ has a FM ground state for $J\le-4$ for all $S$ (see Ref. ). $J=-4$ is a ferromagnetic Lifshitz point where the quadratic term in the dispersion of the magnons vanishes.

Systems with competing FM and AFM interactions are of timely interest, [@shannon06] in particular, the spin-$1/2$ version of Eq. ,[@chubukov91a; @kolezhuk05; @hm06a; @vekua07; @kecke07; @hikihara08; @sudan09] motivated by the experimental realizations in, e.g., LiCuVO$_4$ (Refs. ) and Li$_2$ZrCuO$_4$ (Ref. ). The 1D case of $J<0$ and $S > 1/2$ is largely unexplored (for $J>0$ and $S>1/2$, see Refs.  and ).

We are mainly interested in the region $-4<J<0$ and magnetization $M=S^z/{(S L)}$ ($S^z=\sum_i \langle S^z_i \rangle $) close to saturation $M=1$, for general spin $S$. We will proceed in three steps: First, in Sec. \[sec:two\_magnon\], we discuss the solution of the two-magnon problem and introduce the scattering length. Second, in Sec. \[sec:dilute\], we map the low-energy limit of Eq. , close to saturation, to a dilute 1D gas of two species of bosons interacting via an effective short-range interaction.[@Batyev84; @nikuni95] We then calculate the interaction vertices in this effective theory using a $1/S$ expansion. Finally, in Sec. \[sec:dmrg\], we compare the analytical results with exact numerical ones using the density matrix renormalization group (DMRG) method [@white92b; @schollwoeck05] and exact diagonalization (ED). We put a particular focus on the case of $S=1$. A summary of our results is presented in Section \[sec:summary\], while technical details of the mapping to a dilute gas and of the $1/S$ expansion are given in Appendix \[app:dilute\]. A comparison of DMRG results for open boundary conditions vs.results for periodic boundary conditions is shown in Appendix \[app:dmrg\].

Two-magnon problem {#sec:two_magnon}
==================

Solution of the two-magnon problem
----------------------------------

We now solve the interacting two-magnon problem, starting with the thermodynamic limit (for $S=1/2$, see Ref. ). On a chain of finite length $L$ with periodic boundary conditions the total momentum $K$ is a good quantum number due to the translational invariance of the Hamiltonian $H_S$ in Eq. . Thus, it is convenient to use a basis separating momentum subspaces $$|K,r\rangle=\sum_{l=1}^{L}e^{iK(l+r/2)}S_{l}^{-}S_{l+r}^{-}|F\rangle \, ,
\label{eq:mixed-basis}$$ where $|F\rangle$ is the fully polarized state, $K=2q\pi/L$ ($q=0,1,\cdots,L-1$) and $r$ is the relative distance of two magnons. The allowed values of $r$ depend on $S$ and the parity of $L$ and $q$. For instance, in the case of $S>1/2$ and $L$ even, $r=0, 1,..,L/2-1,(L/2)$ for $q$ odd (even).

We expand a general two-magnon state with momentum $K$ into the (unnormalized) basis of Eq. (\[eq:mixed-basis\]) as $$|\Psi_{2M}\rangle=\sum_{r}C_{r}|K,r\rangle$$ and determine $C_{r}$ analytically by solving the two-magnon Schrödinger equation $$H_S|\Psi_{2M}\rangle=E_{2M}|\Psi_{2M}\rangle.$$ This leads to the recurrence relations $$\begin{aligned}
\label{eq:recurrence}
\Omega_0C_0&=&\frac{S}{\sqrt{S(2S-1)}}(\zeta_{1}C_1 +\zeta_{2}C_2  )\nonumber\\
(\Omega_0-J)C_1&=&\frac{(2S-1)^{3/2}}{{S^{3/2}}}\zeta_{1}C_0+\zeta_{1}C_2 \nonumber\\
&&+\zeta_{2}(C_1+C_3) \nonumber\\
(\Omega_0-1)C_2&=& \frac{(2S-1)^{3/2}}{{S^{3/2}}}\zeta_{2}C_0+\zeta_{2}C_4\nonumber\\
&&+ \zeta_{1} (C_1+C_3)\nonumber\\
\Omega_{0}C_{r} & = & \zeta_{1}\left(C_{r+1}+C_{r-1}\right)\nonumber\\
&&+  \zeta_{2}\left(C_{r+2}+C_{r-2}\right),\quad\mathrm{for}\,\,r\geq3,\end{aligned}$$ where $\zeta_{1}=2SJ\cos{(K/2)}$, $\zeta_{2}=2S\cos{(K)}$. When $|\Psi_{2M}\rangle$ is a bound state, $\Omega_0=E_b-4S(1+J^2/8)$ where $E_b$ is the (negative) binding energy (defined as the bound-state energy minus the energy of the minimum of the two-magnon scattering states).

The (unnormalized) two-magnon bound states for a given $K$ are constructed with the ansatz $$C_{r}=e^{-\kappa_{-}r}+ve^{-\kappa_{+}r}\quad(r\geq1),\label{eq:Ansatz}$$ which, inserted in Eq. (\[eq:recurrence\]), leads to a characteristic quartic equation for $r\geq 3$ $$\Omega_{0}z^{2}-\zeta_{1}(z^{3}+z)-\zeta_{2}(z^{4}+1)=0 \, ,$$ $z$ being any of $e^{-\kappa_{\pm}}$ with $\text{Re}[\kappa_{\pm}]>0$. The remaining unknown quantities $C_{0}$, $v$ and $E_b$ are determined from the remaining relations listed in Eq. (\[eq:recurrence\]).

Scattering length in the lattice problem
----------------------------------------

For $S>1/2$, bound states with energies below the minimum of the two-magnon scattering continuum exist only for $K\simeq \pm 2k_{cl}$ and only in a finite window of couplings $$\label{range_J1c}
    -4 < J < J_{cr}(S),$$ with $S$-dependent critical values $J_{cr}(S)$, as illustrated in Table \[tab:Jcr\]. The critical value $J_{cr}(S)$, which is $J_{cr}\approx -2.11$ for $S=1$, quickly approaches $J_{cr}(S) \simeq -4$ with increasing $S$ (see Tab. \[tab:Jcr\]). In fact, the $1/S$ analysis to be presented in Sec. \[sec:dilute\] suggests the existence of an $S_{cr}$ beyond which this window disappears completely.

       $S$             1            3/2            2            5/2
  -------------- ------------- ------------- ------------- -------------
   $-J_{cr}(S)$   2.11 (2.95)   3.31 (3.42)   3.68 (3.66)   3.84 (3.80)

  : Critical exchange couplings $J_{cr}(S)$ for the existence of metamagnetism in Eq.  derived from solving the two-magnon problem (values in parenthesis: Results from the $1/S$ expansion of Sec. \[sec:dilute\]).[]{data-label="tab:Jcr"}

We define the *scattering length* of bound states, in the thermodynamic limit $L \rightarrow \infty$, from their spatial extent (in analogy to the continuum problem of particles interacting via a short-range, attractive potential): $$\label{boundstatescatteringlength}
a_S =\frac{1}{\mathrm{min} \{\text{Re}[\kappa_{\pm}]\}}.$$ The binding energy takes its lowest value ([*i.e.*]{}, the largest absolute value) for $K=K^*\simeq \pm 2k_{cl}=\pm 2\arccos{(-J/4 )}$ and this quantity, with extremely high accuracy, is related to the scattering length by $$\label{bindingenergy}
E_b(K^*)\simeq-\frac{1}{ma_S^2}\, ,$$ where $m$ is the one-magnon *mass*, $$\label{magnonmass}
m=\frac{2}{S(4-J)(4+J)}.$$ The relation Eq. (\[bindingenergy\]) between the binding energy and the scattering length that holds for our microscopic lattice model is typical for a 1D Bose gas in the continuum interacting via an attractive contact potential, the Lieb-Liniger model.[@LL]

![(Color online) 1D scattering length $a_S/a$ ($a$: lattice spacing) for $S=1/2,1,3/2$ (solid, dashed, dot-dashed line). []{data-label="fig:scatt"}](figure1.eps){width="0.9\columnwidth"}

The scattering length $a_S$ can also be determined from the scattering problem of two magnons for general $S$ in the thermodynamic limit. We have solved this problem, with the momenta of the two magnons (which participate in scattering) taken in the vicinity of the same dispersion minimum, $k_1=k_{cl}+k$ and $k_2=k_{cl}-k$. From the asymptotic form of the two magnon scattering state wavefunction we extract the scattering phase shift $\delta_S (k)$ for any $S$, $$\lim_{r\to \infty}C_r \sim \cos{(rk+\delta_S(k))}.$$ To extract the scattering length from the scattering phase shift we use the same relation as in a 1D continuum model of particles interacting via a short-range potential $$\label{scatteringstatescatteringlength}
 a_S=\lim_{k\to 0} \frac{\cot{(\delta_S(k))}}{k}.$$ This allows us to calculate the scattering length in the repulsive regime $a_S < 0$ as well (when two-magnon bound states are not formed below the minimum of the scattering continuum). In the attractive regime $a_S > 0$, the scattering lengths obtained from both approaches \[[*i.e.*]{}, Eq. (\[boundstatescatteringlength\]) and Eq. (\[scatteringstatescatteringlength\])\] are in excellent agreement with each other. As a side note on terminology, we call attractive (repulsive) regime the one in which the effective interaction between magnons is attractive (repulsive).

We can now generalize the procedure[@Okunishi98] of mapping the antiferromagnetic (unfrustrated) spin-$S$ chain close to saturation onto the low-density limit of the Lieb-Liniger model with a coupling constant $$\label{okunishi}
g_0=-\frac{2}{ma_S} \, . $$

However, in our model the single-magnon dispersion has two minima. The effective theory thus will be a two-component (two species) Lieb-Liniger model. There are two types of low-energy scattering processes, first when momenta of two magnons are in the vicinity of the same dispersion minimum that we have presented above (intraspecies scattering), and second, when the momenta $k_1$ and $k_2$ of two magnons are in the vicinity of different minima of dispersion, [*i.e.*]{}, $k_1=k_{cl}+k$ and $k_2=-k_{cl}-k$ (interspecies scattering). For the latter case we can repeat all steps presented above and extract another coupling constant $\tilde g_0$ from the interspecies scattering length, $\tilde a_S$ in analogy with Eq. (\[okunishi\]). We obtain that $\tilde g_0>0$ (implying that bound states with total momentum $K=0$ are never formed below the scattering continuum), and $\tilde
g_0>g_0$ for any $S$ in the region $-4<J<0$. Since the relation $\tilde g_0>g_0$ always holds, the relevant scattering length at low energies is the intraspecies scattering length $a_S$.

The scattering length, shown in Fig. \[fig:scatt\], can be well described, for small $S>1/2$, by a sum of two terms (resonances): $a_S\simeq \lambda^-_S/(4+J)+\lambda^+_S/(J_{cr}(S)-J)$ \[where $\lambda^{\pm}_S$ are numerical prefactors\]. We emphasize that, for $S\geq 1$, the scattering length is in general much larger than the lattice spacing, as is evident from Fig. \[fig:scatt\]: For $S=1$, $a_S$ takes a minimum at $J\simeq -3.3$ with $ a_S\simeq 80a$. Additionally, the emergence of bound states manifests itself by a diverging scattering length at $J_{cr}(S)$, where $a_S$ changes its sign jumping from $-\infty$ to $+ \infty$. Thus, bound states are typically shallow, with a binding energy given by Eq. (\[bindingenergy\]). In addition, the minima in their dispersion occur at incommensurate momenta $K^*$.

For $S=1/2$, any $J<0$ induces a two-magnon bound state with total momentum $K=\pi$, [@chubukov91a; @hm06a] [*i.e.*]{}, $a_S > 0$ and there is no resonance at $-4 < J<0$ for $S=1/2$. Hence, $S=1/2$ is very different from the $S>1/2$ case where bound states with $K=\pi$ are never below the two-magnon scattering continuum, and, as discussed above, a resonance exists for $-4 < J<0$ and $1\leq S < S_{cr}$.

In order to analyze the finite-size effects with respect to results in the thermodynamic limit, we have numerically diagonalized $H_S$ in the basis given in Eq. (\[eq:mixed-basis\]) (see Ref.  for details of the procedure). We obtained the full spectrum, *i.e.*, the scattering continuum and bound/antibound states if present, for selected values of $J \in[-4,0]$ in systems with up to $L=4000$ sites and several values of $S$. Table \[tab:Num\_Jcr\] shows the numerical determination of $J_{cr}(S)$ for $S=1$ and different system sizes. Although finite-size effects are apparent, a quadratic fit in $1/L$ to numerical data for $J_{cr}(S)$ extrapolates to $J_{cr}\simeq -2.11$ in agreement with the result determined directly in the thermodynamic limit (see the preceding discussion and Table \[tab:Jcr\]).

        $L$          1000     2000     4000    $\infty$
  ---------------- -------- -------- -------- ----------
   $-J_{cr}(S=1)$   $2.25$   $2.16$   $2.13$    $2.11$

  : Finite-size dependence of critical values $J_{cr}(S)$ for the emergence of bound states below the minimum of the two-magnon continuum of scattering states, for $S=1$.[]{data-label="tab:Num_Jcr"}

Mapping of the spin Hamiltonian to a dilute gas of bosons {#sec:dilute}
=========================================================

In this section, we describe our effective theory in the thermodynamic limit, for the case of a finite (though vanishingly small) density of magnons. The mapping to a dilute gas of bosons is motivated by the following observation: For $S>1/2$, we have shown that $a_S$ is large. Hence in the dilute limit, we can safely neglect the hard-core constraint and take the continuum limit. For $S=1/2$, on the contrary, the scattering length $a_S$ is typically of the order of a few lattice constants and only for $-4<J<-3.9$ does $a_S$ become comparable to the smallest value of the scattering length for $S=1$.

For $S\geq 1$, close to saturation, and in the low-energy limit, we therefore map our system onto a dilute two-component gas of bosons interacting with effective short-range interactions. Many-body effects will be incorporated by properly shifting the two-body T matrix off-shell as explained in Ref. . We show that, while the [*inter*]{}species interaction is always repulsive and stronger than the [*intra*]{}species interaction, the latter undergoes a sign change. When the intraspecies interaction becomes negative, the bosons are unstable against a collapse. We show that a $1/S$ expansion captures this physics correctly and, similar to the case of the (unfrustrated) Heisenberg chain, [@Johnson] is applicable to the present problem, albeit its one-dimensional nature.

Effective Hamiltonian
---------------------

Using the Dyson-Maleev transformation [@DysonMaleev] (the Dyson-Maleev representation is used here for convenience. We have checked that the explicitly hermitean Holstein-Primakoff representation, [@HolsteinPrimakoff] to leading order $1/S$, provides equivalent results) $$\begin{aligned}
S_i^z&=&S-a_i^{\dagger}a_i \, ,\quad S_i^+=  \sqrt{2S}a_i \, , \nonumber\\
S_i^-&=&\sqrt{2S}a_i^{\dagger}(1-a_i^{\dagger}a_i/2S)\,,$$ we map Eq.  onto a bosonic problem: $$\label{DysonMaleev}
H=\!\sum_k(2S\epsilon_k-\mu)a_k^{\dagger}a_k
+\!\!\!\sum_{k,k',q} \!\!\! \frac{\Gamma_0(q;k,k')}{2L} a_{k+q}^{\dagger}a_{k'-q}^{\dagger}a_{k}a_{k'} \, ,$$ where $$\epsilon_{k}=J\cos{k}+\cos{2k}-(J\cos{k_{cl}}+\cos{2k_{cl}})\ge 0$$ is the single-magnon dispersion and $k_{cl}=\arccos{(-J/4)}$. Note that in our normalization, the minima of the single-particle dispersion are at zero energy: $\epsilon_{\pm k_{cl}}=0$. The bare interaction vertex $\Gamma_0$ is given by $\Gamma_0(q;k,k') = V_q-\frac{1}{2}(V_k+V_{k'})$ with $V_k= 2J\cos{k}+2\cos{2k}$. The chemical potential is $\mu=h^{cl}_s-h$, where $h_s^{cl}$ is the classical saturation field value $h^{cl}_s=S(J+4)^2/4$. We are interested in the dilute regime $\mu\to 0$.

Concentrating on the low-energy behavior we arrive at, via a Bogoliubov procedure, [@Batyev84; @nikuni95] a two-component Bose gas interacting via a $\delta$-potential with Hamiltonian density $$\mathcal{H}_{\rm eff}=\sum_{\alpha}-\frac{|\nabla \psi_{\alpha}|^2}{2m}+\frac{g_0(S)}{2}(n_1^2+n_2^2)+\tilde g_0(S) \, n_1 n_2\, .
\label{eq:effective}$$ Here $\psi_\alpha$, $\alpha=1,2$ describe bosonic modes with momenta close to $\pm k_{cl}$ ([*i.e.*]{}, the Fourier transforms of $\psi_\alpha$ are $\psi_1(k\!\!\to\!\!0) \approx a_ {k_{cl}+k}$, $\psi_2(k\!\!\to\!\!0) \approx a_ {-k_{cl}+k}$) and $n_{\alpha}=\psi_\alpha^{\dagger}\psi_\alpha$ are the corresponding densities. The bare coupling constants of the effective 1D model of the two-component Bose gas, $g_0(S)$ and $\tilde g_0(S)$, are, in the dilute limit of bosons, related to the renormalized vertices of the microscopic model Eq. (\[DysonMaleev\]) through $$\begin{aligned}
\label{microeffective}
&&  \Gamma(0;k_{cl},k_{cl})=  \frac{g_0(S)}{1+g_0(S){\sqrt{2m}}/({\pi \sqrt{\mu}})} ,\\
&&\Gamma(0;k_{cl},\!-k_{cl})+\Gamma(-2k_{cl};k_{cl},\!-k_{cl})\!=\!\frac{\tilde g_0(S)} {1+\frac{\tilde g_0(S){\sqrt{2m}}}{\pi \sqrt{\mu}}}. \nonumber
 \end{aligned}$$ The relations Eq. (\[microeffective\]) follow from a generalization of the corresponding equation for the case of a one-component Bose gas [@kolomeisky92; @Lee2002] to the two-component case using an RG analysis[@KolezhukRG] (see Appendix \[app:dilute\] for details of the calculation).

$1/S$ expansion {#subsect:1overS}
---------------

Next, we apply a $1/S$ expansion to calculate the interaction vertices and extract the coupling constants $g_0(S)$ and $\tilde g_0(S)$. Using a standard ladder approximation the Bethe-Salpeter equation for the vertices $\Gamma$ reads: $$\begin{aligned}
\label{BetheSalpeterEquation}
\Gamma(q;k,k')&=& \Gamma_0(q;k,k')  \\
&-&\frac{1}{2SL}\sum_{p} \frac{ \Gamma_0(q-p;k+p,k'-p) } {\epsilon_{k+p}+\epsilon_{k'-p}} \Gamma(p;k,k')\nonumber\,.\end{aligned}$$ Setting the transferred momentum $q=0$ in $\Gamma$ and the incoming momenta to $k=k'=k_{cl}$, we get (see Appendix \[app:dilute\] for details): $$\begin{aligned}
\label{nine}
&&\Gamma(0;k_{cl},k_{cl})\left[1+\frac{V_0-V_{k_{cl}}}{2SL}  \sum_{p} \frac{ 1 } {\epsilon_{k_{cl}+p}+\epsilon_{k_{cl}-p}}\right]=\nonumber\\
&&V_0-V_{k_{cl}}+\frac{1}{2SL}\sum_{p}\! \left[1-\frac{  V_{p}-V_0 } {\epsilon_{k_{cl}+p}+\epsilon_{k_{cl}-p}} \right]\Gamma(p;k_{cl},k_{cl})\nonumber\\
&&-\frac{1}{2SL}\sum_{p} \frac{ (V_0-V_{k_{cl}}) \left[\Gamma(p;k_{cl},k_{cl})-\Gamma(0;k_{cl},k_{cl})\right] } {\epsilon_{k_{cl}+p}+\epsilon_{k_{cl}-p}}.\nonumber\\\end{aligned}$$ Now, in the spirit of the $1/S$ expansion, we replace the renormalized vertex with the bare vertex on the right hand side of Eq. (\[nine\]), $\Gamma(p;k_{cl},k_{cl})\to 2\epsilon_p$, which is possible since there are no infrared divergences. Regularizing the left hand side of Eq. (\[nine\]) as in Refs. , and using Eq. (\[microeffective\]) we extract the coupling constants of the effective model. The analytical expression for $g_0(S)$ is $$\label{LLC_wc1}
g_0(S)=\frac{F}{1-\frac{F(J^2-8)}{|J|S(16-J^2)^{3/2}}}$$ (the derivation of Eq.  can be found in Appendix \[app:dilute\] and the constant $F$ is given in Eq. ). To leading order in $J+4$, $$\lim _{J\to -4^+} g_0(S) \simeq \frac{S-S_{cr}}{4S}(J+4)^2+O\left((J+4)^{5/2}\right)\,.$$ To first order in $1/S$, we obtain $S_{cr}=6$, which is not that large a number, hence corrections beyond $1/S$ may affect $S_{cr}$. In the same way we calculate $\tilde g_0(S)$ as $$\label{g0tilde}
\tilde g_0(S)=\frac{\tilde F}{1+\frac{J^2-8}{16S}}$$ (the derivation is presented in detail in Appendix \[app:dilute\]; see Eq.  for the expression for $\tilde F$).

From Eqs.  and , we notice that $g_0(S) < \tilde g_0(S)$ for $-4<J<0$. Thus the state below saturation is a single-component one. Provided the interactions are repulsive, the ground state is a translationally invariant chiral state,[@kolezhuk05] where bosons prefer to ‘condense’ at the same minimum of the single-particle dispersion since they experience a minimal repulsion there.[@jackeli04] We also note that $|g_0|m\ll 1$ for $J<0$, hence interactions between bosons are generically weak. In particular, even though $m\to \infty$ when approaching the ferromagnetic Lifshitz point, $|g_0|m\to 0$.

![(Color online) Effective bare intraspecies interaction $g_0(S)$, for $S=3$ (solid line), representative of the generic behavior for $S<S_{cr}$, $S=10$ (dashed curve), representative of $S>S_{cr}$ \[dot-dashed line: $g_0(S=\infty)$\]. Inset: $g_0$ for $S=1$ and $3/2$. []{data-label="fig:geff"}](figure2.eps){width="0.9\columnwidth"}

The effective bare intraspecies interaction is depicted in Fig. \[fig:geff\] and behaves as $g_0 (S) \sim [J-J_{cr}(S)]$ for $J\to J_{cr}(S)$. The scattering length is related to the effective coupling constant by Eq. , signaling a resonance at $J_{cr}(S)$. Thus, we see that for $S<S_{cr}$ there is a finite region near $J\simeq-4$ where $g_0 (S)<0$ and bosons attract each other, producing a collapsed state.

To corroborate this, using ED for Eq.  and $S=1$ with periodic boundary conditions, we have calculated the ground-state momentum of the states with a small, but finite number of magnons, which is incommensurate, supporting the picture of a uniform chiral state in the repulsive case $J>J_{cr}$, and a collapsed state in the attractive case $J<J_{cr}$ at one of the two minima of the single-particle dispersion. In the attractive case, $\partial^2E_0/\partial n^2<0$, where $E_0$ and $n $ are the bosons’ ground-state energy and density, respectively. In the language of spins, the inverse magnetic susceptibility at saturation becomes negative, and hence, following standard arguments,[@dmitriev01] we conclude that there is a first-order transition at $M=1$, [*i.e.*]{}, a jump in the magnetization curve just below saturation.

As pointed out above, the case of $S=1/2$ is special since the scattering length is typically of the order of the lattice constant here. The mapping of the $S=1/2$ case to a two-component Bose gas (by the procedure presented above for $S>1/2$) can be trusted only for $J \to -4$, where the scattering length becomes much larger than the lattice constant. In that case we can easily incorporate the exact hard-core constraint into our formalism [@nikuni95] and again expect that $S=1/2$ also shows metamagnetic behavior. This conclusion is in agreement with DMRG results for $S=1/2$.[@sudan09] Note that a metamagnetic jump can also be stabilized for spin $1/2$ with suitable anisotropic exchange interactions.[@GMK98; @hirata99] However, with our procedure we cannot account for the formation of stable two-, three-, and four-magnon bound states that is characteristic for most of the region $J>-4$ in the spin-$1/2$ frustrated ferromagnetic Heisenberg chain.[@hm06a; @kecke07; @hikihara08; @sudan09]

Going back to $S>1/2$, at lower $M$, corresponding to higher densities of magnons, the hard-core nature of spins eventually prevails as well, resulting in a uniform ground state at a nonzero momentum. However, as already mentioned, from the finite-size analysis of the two-magnon problem, we observe that bound states disappear with decreasing $L$, suggesting that the attractive effective potential (in the limit of a small magnon density) can become repulsive upon increasing the magnon densities. Thus, the state below the jump ([*i.e.*]{}, $0<M<1-\Delta M_{\mathrm{jump}}$, where $\Delta M_{\mathrm{jump}}$ is the height of the jump) will be similar to the one encountered in the case of $J>J_{cr}$, [*i.e.*]{}, it is a translationally uniform chiral state.

![(Color online) (a) Magnetization curves $M(h)$ for $S=1$ at $J=-2.5,-3,-3.5$. (b) $M(h)$ for $S=1$, $J=-1$ (c) $M(h)$ for $S=3/2$ at $J=-3.5$ (all for $L=128$).[]{data-label="fig:mag_h"}](figure3a.eps "fig:"){width="0.85\columnwidth"}\
![(Color online) (a) Magnetization curves $M(h)$ for $S=1$ at $J=-2.5,-3,-3.5$. (b) $M(h)$ for $S=1$, $J=-1$ (c) $M(h)$ for $S=3/2$ at $J=-3.5$ (all for $L=128$).[]{data-label="fig:mag_h"}](figure3bc.eps "fig:"){width="0.9\columnwidth"}

DMRG results {#sec:dmrg}
============

Next, we turn to numerical results for the case of $S=1$ (unless stated otherwise), solving for the ground state of Eq.  in a finite magnetic field $h$, using ED where possible or DMRG.[@fn-alps] We present data from DMRG simulations using up to $1200$ states, for $L\leq 128$ sites, and for open boundary conditions (OBC), unless stated otherwise.

Magnetization curves and magnetization profiles
-----------------------------------------------

![(Color online) Magnetization profiles for $S=1$, $L=128$ at $J=-1$ (no jump, dashed lines) and $J=-3$ (jump at $S^z=97$) for $S^z=80,90,110,120$ (top to bottom).[]{data-label="fig:mag_prof"}](figure4.eps){width="0.9\columnwidth"}

The main result of this work, namely, the metamagnetic transition from a gapless finite-field phase to full saturation, is clearly seen in the magnetization curves shown in Fig. \[fig:mag\_h\](a). For $S=1$, we observe the appearance of this jump for $-4< J \lesssim -2$, while for $S=3/2$ \[an example is shown in Fig. \[fig:mag\_h\](c)\] the jump exists in a much narrower window $J\lesssim -3 $. This is consistent with our analytical results for $J_{cr}$, listed in Tab. \[tab:Jcr\]. Moreover, for $-2 \lesssim J<0 $, we resolve a plateau in $M(h)$ at $M=0$, which is due to the Haldane gap[@whitehuse] in this $S=1$ system \[see Fig. \[fig:mag\_h\](b)\]. This gap defines the critical field $h_c$ that separates the gapped Haldane phase from the finite-$M$ phase with a smooth $M(h)$-behavior for $h_c<h<h_{\mathrm{sat}}$.

The collapse of magnons manifests itself in the magnetization profiles ($1-\langle S^z_i \rangle$ vs. site $i$) using OBC, displayed in Fig. \[fig:mag\_prof\] for $J=-3$ and $J=-1$. In the former case, there is a jump, but in the latter, there is none. Clearly, in the states that get skipped over ($S^z>97$ for $J=-3$), magnons collapse into the center of the system, whereas for actual ground states below the metamagnetic transition, the magnetization profiles become flat. By contrast, in the case of $J=-1$ where the transition to the fully polarized state is smooth and continuous, [*all*]{} profiles are, apart from boundary effects, flat.

Central charge
--------------

Our effective theory developed in Section \[sec:dilute\] suggests that the gapless phase in the region $h_{c}<h<h_{\mathrm{sat}}$ is a one-component phase (where $h_{\mathrm{sat}}$ is the saturation field). To substantiate this result, one can make use of entanglement measures such as the von-Neumann entropy to extract the central charge, which directly yields the number of components of the gapless state.

The von-Neumann entropy is defined as $$S_{vN}(l)=- \mbox{tr}(\rho_l \ln \rho_{l}) \, ,$$ where $\rho_l$ is the reduced density matrix of a subsystem of length $l$ of our one-dimensional chain of length $L$. In a gapless state that is conformally invariant, the $l$ and $L$ dependence of the von-Neumann entropy is given by[@vidal03; @calabrese] $$S_{vN}(l) =\frac{c}{3} \ln\left( \frac{L}{\pi} \sin(\frac{\pi}{L} l)\right)+g\,,
\label{eq:cch}$$ which is valid for systems with periodic boundary conditions (PBC). PBC are preferable for the calculation of the central charge from Eq.  since for OBC, there may be additional oscillatory terms. $g$ is a non-universal constant that depends on $M$. As DMRG directly accesses the eigenvalues of these reduced density matrices,[@schollwoeck05] it is straightforward to measure $S_{vN}(l)$ with this numerical method.

Some typical DMRG results (squares) for systems of $L=64$ and periodic boundary conditions are presented in Fig. \[fig:svn\]. We have fitted the expression Eq.  to our numerical data (shown as solid lines in the figure) and obtain $c=1.0\pm 0.1$ in all examples. Therefore, we expect the gapless phase to be a (chiral) one-component liquid. Note, though, that at both small $M$ and $|J|$, where the convergence of DMRG is notoriously difficult, we cannot completely rule out the presence of a $c=2$ region, which, however, is irrelevant for the main conclusions of our work.

![(Color online) DMRG results for the von-Neumann entropy $S_{vN}(l)$ in the gapless phase $h_{c}<h<h_{\mathrm{sat}}$ of the $S=1$ system: (a) $J=-1$, $M=1/2$, (b) $J=-2$, $M=1/2$, (c) $J=-3$, $M=5/16$ (symbols). The lines are fits to Eq. , resulting in $c=1.0\pm0.1$ in all cases (we exclude $S_{vN}(l)$ for $l< 10$ and $l>54$ from the fit). In this figure, we display results for periodic boundary conditions and $L=64$ sites. []{data-label="fig:svn"}](figure5.eps){width="0.9\columnwidth"}

Phase diagram for $S=1$
-----------------------

Our results for the $S=1$ chain are summarized in the $h$ vs. $J$ phase diagram Fig. \[fig:phase\]. We identify three phases: (i) a gapped $M=0$ phase at $h<h_c$ (similar to the Double-Haldane phase known for $J>0$, see Ref. ), (ii) a gapless (chiral) finite-field phase for $h_c<h<h_{\mathrm{sat}}$, and (iii) the fully polarized state at $h_{\mathrm{sat}}<h$ (with $h_{\mathrm{sat}}=0$ for $J<-4$).

$\Delta M_{\mathrm{jump}}$ is plotted in the inset of Fig. \[fig:phase\]: the jump sets in at $J\lesssim -2 $ (close to where the zero-field gap becomes small rendering it difficult to resolve it numerically), consistent with our theory.

Since in the limit of $J\to 0$, one has two spin-1 chains with antiferromagnetic interactions which both separately have a Haldane gap at zero field,[@whitehuse] upon coupling the chains, one obtains the so-called Double-Haldane phase (in contrast to the regular Haldane phase that is inherited from a single spin-1 chain with antiferromagnetic interactions). Both phases, Double-Haldane and Haldane phase, are realized in the frustrated, antiferromagnetic spin-1 chain,[@kolezhuk96] yet in our case, only the Double-Haldane phase exists. The determination of the corresponding spin gap $h_c$ is a bit subtle. Namely, a Haldane chain with open boundaries gives rise to spin-1/2 excitations at the open ends.[@kennedy90] Since we have two chains (for small $|J|$), we have a total of four spin-1/2 end spins. Hence, the spin gap in Fig. \[fig:phase\] is determined from $$h_c = E(S^z=3) - E(S^z=2) \, ,
\label{eq:HcOBC}$$ where $E(S^z)$ is the ground-state energy in a sector with a given total $S^z$.

It is worth emphasizing several differences with the phase diagram of the spin-1/2 version of Eq. . First, for $S=1$, there are no multipolar phases, which occupy a large portion of the corresponding spin-1/2 phase diagram.[@hm06a; @vekua07; @kecke07; @hikihara08; @sudan09] Second, the spin-1/2 system features an instability towards nematic order,[@chubukov91a; @vekua07; @kecke07] which can be excluded on general grounds for integer spin-$S$ chains,[@kolezhuk05] even for $0<M\ll 1$.

![(Color online) Phase diagram of Eq.  for $S=1$ (circles: saturation field $h_{\mathrm{sat}}$; squares: spin gap $h_{c}$). Inset: height $\Delta M_{\mathrm{jump}}$ of the metamagnetic jump vs. $J$. The comparison of $L=64$ (open symbols) and $L=128$ (solid symbols) as well as finite-size scaling (not shown here) supports that both $\Delta M_{\mathrm{jump}}$ and $h_c$ are finite in extended regions of $J$.[]{data-label="fig:phase"}](figure6.eps "fig:"){width="0.9\columnwidth"}\[t!\]

Summary {#sec:summary}
=======

In conclusion, we showed that resonances can play a crucial role in determining the low-energy behavior of frustrated quantum spin systems subject to a magnetic field. The proximity of resonances caused by an interplay between frustration ($J>-4$) and quantum fluctuations ($1/2<S 
<S_{cr}$) results in extremely large values of the 1D scattering length that allows to develop an effective theory of a weakly interacting two-component Bose gas. The quasi-collapse of the dilute gas of magnons provides the physical origin of the emergent metamagnetism.

The preditictions of our analytic theory were verified by numerical data. We focused on the case $S=1$ since there the jump in the magnetization curve is the most pronounced. As a by-product, we obtained the phase diagram for the $S=1$ $J$-$J'$ chain with antiferromagnetic $J'$ and ferromagnetic $J$ in a magnetic field. This phase diagram is remarkably simple. In particular, there are no multipolar phases in the case $S=1$, in marked contrast to the spin-1/2 case.[@hm06a; @vekua07; @kecke07; @hikihara08; @sudan09]

We thank H. Frahm, A. Kolezhuk, R. Noack, and D. Petrov for fruitful discussions. A.H. acknowledges financial support from the DFG via a Heisenberg fellowship (HO 2325/4-2). T.V. is supported by the Center of Excellence QUEST. G.R. and M.A. are partially supported by CONICET (PIP 1691) and ANPCyT (PICT 1426).

Details on the mapping to a dilute Bose gas {#app:dilute}
===========================================

In Sec. \[sec:dilute\], we mapped the spin Hamiltonian Eq.  close to saturation to an effective, bosonic field theory Eq. . Here we provide the details of generalizing the procedure of obtaining the coupling constants from the many-body T-matrix for an effective single-component Bose gas model, which is described in Ref. , to the case relevant to us here, namely an effective theory of a two-component Bose gas.

We introduce the ‘mean-field interaction coefficients’, $$\begin{aligned}
\label{effectiveinteraction}
g(S)&=&\Gamma(0;k_{cl},k_{cl}) \\
\tilde g(S)&=&\Gamma(0;k_{cl},-k_{cl})+\Gamma(-2k_{cl};k_{cl},-k_{cl}) \, ,\end{aligned}$$ where $\Gamma(q;k,k')$ is the full interaction vertex of the Hamiltonian Eq. (\[DysonMaleev\]). They satisfy mean-field like relations for $\tilde g_{0}>g_0$, $$\mu=g(S)n$$ and for $\tilde g_{0}<g_0$ $$\mu=(g(S)+\tilde g(S))n/2,$$ where $n=\langle\sum_{\alpha}n_{\alpha} \rangle\to0$ is the total density of bosons. To get the connection between $g(S)$ and $\tilde g(S)$ on the one hand and the bare coupling constants of the effective 1D model of a two-component Bose gas, $g_0(S)$ and $\tilde g_0(S)$ on the other hand, we generalize the corresponding equation for the case of a one-component Bose gas [@Lee2002] to the two-component case, consistent with the $SU(2)$ symmetry of the RG fixed point of a dilute, two-component Bose gas,[@KolezhukRG] $\lim_{\mu\to 0}g(S)= \lim_{\mu\to 0} \tilde g(S)=\pi\sqrt{\mu}/\sqrt{2m}$, $$\begin{aligned}
\label{microeffectiveapp}
&&  g(S)= \frac{g_0(S)}{1+g_0(S){\sqrt{2m}}/({\pi \sqrt{\mu}})},\\
&&\tilde g(S)\!=\! \frac{\tilde g_0(S)}{1+\tilde g_0(S){\sqrt{2m}}/({\pi \sqrt{\mu}})}. \label{microeffective_a}
 \end{aligned}$$ To zeroth order in $1/S$, we have $$g(S=\infty)=V_0-V_{k_{cl}}=\frac{(J+4)^2}{4}$$ and $$\tilde g(S=\infty)= V_{2k_{cl}}+V_0-2V_{k_{cl}}>g(S=\infty) \, .$$ Thus, classically ($S\to \infty, m\to 0$), we have $\tilde g_{0}>g_0$ for $-4<J<0$. We will show below that incorporating quantum fluctuations does not modify this relation. In the following we will calculate the effective intraspecies and interspecies interaction constants, to see how $1/S$ corrections modify $g_0(S)$ and $\tilde g_0(S)$. In order to obtain $g(S)$ we have to set in $\Gamma (q,k,k')$ given by Eq. (\[BetheSalpeterEquation\]), $q=0$ and $k=k'=k_{cl}$, $$\begin{aligned}
g(S)&=& V_0-V_{k_{cl}} +\frac{1}{2S L}\sum_{p} \Gamma_p\nonumber\\
&&-\frac{1}{2S L}\sum_{p} \frac{  V_{p}-V_0+V_0-V_{k_{cl}} } {\epsilon_{k_{cl}+p}+\epsilon_{k_{cl}-p}} \Gamma_p,\end{aligned}$$ where we have denoted $ \Gamma(p;k_{cl},k_{cl})=\Gamma_p$. After straightforward manipulations we obtain, $$\begin{aligned}
\label{nineapp}
&&g(S)\left[1+\frac{V_0-V_{k_{cl}}}{2S L}  \sum_{p} \frac{ 1 } {\epsilon_{k_{cl}+p}+\epsilon_{k_{cl}-p}}\right]=\nonumber\\
&&\quad V_0-V_{k_{cl}} +\frac{1}{2S L}\sum_{p} \Gamma_p \left[1-\frac{  V_{p}-V_0
} {\epsilon_{k_{cl}+p}+\epsilon_{k_{cl}-p}}\right] \qquad \nonumber\\
&&\quad -\frac{1}{2S L}\sum_{p} \frac{ (V_0-V_{k_{cl}}) (\Gamma_p-\Gamma_0) }
{\epsilon_{k_{cl}+p}+\epsilon_{k_{cl}-p}} \, .\end{aligned}$$ Now, on the right hand side of Eq. (\[nineapp\]), we plug in the zeroth order vertex (in the $1/S$ expansion), $\Gamma(p;k_{cl},k_{cl})\to 2\epsilon_p$ (which is possible because there are no infrared divergences anymore), and use that $V_p-V_0=2(\epsilon_p-\epsilon_0)$. Equation (\[nineapp\]), after passing to infinite system size, becomes $$\begin{aligned}
\label{unreg}
g(S)\left[1+\frac{F}{4\pi S} \int\limits_{-\pi}^{\pi} \frac{ 1 }
{\epsilon_{k_{cl}+p}+\epsilon_{k_{cl}-p}}\right]=F \, ,\end{aligned}$$ where $$\begin{aligned}
F&=& V_0- (1+\frac{1}{2S})V_{k_{cl}} - \frac{1}{\pi S}\int\limits_{-\pi}^{\pi}{\mathrm d}p \frac{   \epsilon_p (\epsilon_p-\epsilon_0) } {\epsilon_{k_{cl}+p}+\epsilon_{k_{cl}-p}} \nonumber\\
&&- \frac{ V_{0}-V_{k_{cl}} }{2\pi S}\int\limits_{-\pi}^{\pi}{\mathrm d}p \frac{  \epsilon_p-\epsilon_0 } {\epsilon_{k_{cl}+p}+\epsilon_{k_{cl}-p}}
\label{eq:valueF}\end{aligned}$$ and we have used that $V_0-V_{k_{cl}}=F+O(1/S)$. According to the scheme that we follow the two-body T matrix must be calculated off-shell, thus the denominator in Eq. (\[unreg\]) must be understood as $\epsilon_{p}\to \epsilon_{p}+C\mu/4S$, where the exact value of the numerical constant is $C=\pi^2/8$.[@Lee2002] This leads to (compare Eqs.  and ) $$\label{renormalizedpotential2}
g(S)= \frac{g_0(S)}{1+g_0(S){\sqrt{2m}}/({\pi \sqrt{\mu}})},$$ where we have introduced the intraspecies Lieb-Liniger coupling constant $g_0(S)$ as in Eq. (\[LLC\_wc1\]). Note that all integrals presented in this section are evaluated analytically, which is a nice feature of the $1/S$ treatment of our problem.

Now we outline the calculation of $\tilde g(S)$. We denote $\Gamma(p;k_{cl},-k_{cl})=\tilde \Gamma_p$ and obtain the Bethe-Salpeter equation for $$\tilde \Gamma_0
=V_0 -V_{k_{cl}}   + \frac{1}{4\pi S}  \int\limits_{-\pi}^{\pi} {\mathrm d} p \tilde \Gamma_p
- \frac{1}{4\pi S}  \int\limits_{-\pi}^{\pi} {\mathrm d} p \frac{  V_{-p}-V_{k_{cl}} } {2\epsilon_{k_{cl}+p}}\tilde \Gamma_p
\nonumber$$ and $$\begin{aligned}
\tilde \Gamma_{-2k_{cl}}
&=&V_{-2k_{cl}} -V_{k_{cl}}   + \frac{1}{4\pi S}  \int\limits_{-\pi}^{\pi} {\mathrm d} p  \tilde \Gamma_p\nonumber\\
&&- \frac{1}{4\pi S}  \int\limits_{-\pi}^{\pi} {\mathrm d} p \frac{  V_{-2k_{cl}-p}-V_{k_{cl}} } {2\epsilon_{k_{cl}+p}}\tilde \Gamma_p .\end{aligned}$$ Adding these two equations gives $\tilde g(S)=\tilde \Gamma_0+\tilde \Gamma_{-2k_{cl}}$, $$\begin{aligned}
\label{gamma2}
\tilde g(S)
&=&V_0+V_{-2k_{cl}} -2V_{k_{cl}}   + \frac{2}{4\pi S}  \int\limits_{-\pi}^{\pi}
{\mathrm d} p  \tilde \Gamma_p \qquad \nonumber\\
&&- \frac{1}{4\pi S}  \int\limits_{-\pi}^{\pi} {\mathrm d} p \frac{
V_p+V_{-2k_{cl}-p} -2V_{k_{cl}} } {2\epsilon_{k_{cl}+p}}\tilde \Gamma_p.\end{aligned}$$

We divide the last term in Eq. (\[gamma2\]) into two pieces, $$-\frac{1}{4\pi S}  \int\limits_{-\pi}^{\pi} {\mathrm d} p \frac{
V_p+V_{-2k_{cl}-p}-2V_{k_{cl}} } {2\epsilon_{k_{cl}+p}}\tilde \Gamma_p=I_1+I_2
\, ,$$ where $$\begin{aligned}
I_1&=& - \frac{1}{4\pi S}  \int\limits_{-\pi-k_{cl}}^{-k_{cl}} {\mathrm d} p \frac{  V_p+V_{-2k_{cl}-p}-2V_{k_{cl}} } {2\epsilon_{k_{cl}+p}}\nonumber\\
&& \times (\tilde \Gamma_p- \tilde \Gamma_{-2k_{cl}}+\tilde \Gamma_{-2k_{cl}})
       \end{aligned}$$ and $$\begin{aligned}
       I_2&=& - \frac{1}{4\pi S}  \int\limits_{-k_{cl}}^{\pi-k_{cl}} {\mathrm d} p \frac{  V_p+V_{-2k_{cl}-p}-2V_{k_{cl}} } {2\epsilon_{k_{cl}+p}}\nonumber\\
&& \times ( \tilde \Gamma_p- \tilde \Gamma_{0}+\tilde \Gamma_{0}) \, .
       \end{aligned}$$ Note that, for convenience, we have shifted the first Brillouin zone, $(-\pi,\pi)\to (-\pi-k_{cl},\pi-k_{cl})$. Shifting the T-matrix off-shell, we have (the integral with a dash denotes the principal value) $$\begin{aligned}
I_1&=& - \frac{1}{4\pi S}  \! \dashint\limits_{-\pi-k_{cl}}^{-k_{cl}}
  \!\!\!\! {\mathrm d} p \frac{  V_p+V_{-2k_{cl}-p}-2V_{k_{cl}} } {2\epsilon_{k_{cl}+p}}( \tilde \Gamma_p- \tilde \Gamma_{-2k_{cl}})\nonumber\\
&&- \frac{ \tilde \Gamma_{-2k_{cl}}}{4\pi S}  \int\limits_{-\pi-k_{cl}}^{-k_{cl}} {\mathrm d} p \frac{  V_p+V_{-2k_{cl}-p}-2V_{k_{cl}} } {2\epsilon_{k_{cl}+p}+C\mu/2S}
       \end{aligned}$$ and $$\begin{aligned}
       I_2&=& - \frac{1}{4\pi S}  \int\limits_{-k_{cl}}^{\pi-k_{cl}} {\mathrm d} p \frac{  V_p+V_{-2k_{cl}-p}-2V_{k_{cl}} } {2\epsilon_{k_{cl}+p}}( \tilde \Gamma_p- \tilde \Gamma_{0})\nonumber\\
&&- \frac{ \tilde \Gamma_{0}}{4\pi S}  \int\limits_{-k_{cl}}^{\pi-k_{cl}} {\mathrm d} p \frac{  V_p+V_{-2k_{cl}-p}-2V_{k_{cl}} } {2\epsilon_{k_{cl}+p}+C\mu/2S}.
       \end{aligned}$$ The last terms in $I_1$ can be written as $$\begin{aligned}
&-& \frac{ \tilde \Gamma_{-2k_{cl}}}{4\pi S}  \int\limits_{-\pi-k_{cl}}^{-k_{cl}} {\mathrm d} p \frac{  V_p+V_{-2k_{cl}-p}-2V_{k_{cl}} } {2\epsilon_{k_{cl}+p}+C\mu/2S}=\nonumber\\
&&- \frac{ \tilde \Gamma_{-2k_{cl}}}{4\pi S}
\dashint\limits_{-\pi-k_{cl}}^{-k_{cl}} {\mathrm d} p \frac{
V_p+V_{-2k_{cl}-p}-(V_0+V_{-2k_{cl}}) } {2\epsilon_{k_{cl}+p}}\quad\nonumber\\
&&- \frac{ \tilde \Gamma_{-2k_{cl}}}{4\pi S}  \int\limits_{-\pi-k_{cl}}^{-k_{cl}} {\mathrm d} p \frac{  V_0+V_{-2k_{cl}}-2V_{k_{cl}}} {2\epsilon_{k_{cl}+p}+C\mu/2S}.
       \end{aligned}$$ Similarly, for the last terms in $I_2$ we get, $$\begin{aligned}
&-& \frac{ \tilde \Gamma_{0}}{4\pi S}  \int\limits_{-k_{cl}}^{\pi-k_{cl}} {\mathrm d} p \frac{  V_p+V_{-2k_{cl}-p}-2V_{k_{cl}} } {2\epsilon_{k_{cl}+p}+C\mu/2S}=\nonumber\\
&&- \frac{ \tilde \Gamma_{0}}{4\pi S}  \dashint\limits_{-k_{cl}}^{\pi-k_{cl}}
{\mathrm d} p \frac{  V_p+V_{-2k_{cl}-p}-(V_0+V_{-2k_{cl}}) }
{2\epsilon_{k_{cl}+p}}\quad\nonumber\\
&&- \frac{ \tilde \Gamma_{0}}{4\pi S}  \int\limits_{-k_{cl}}^{\pi-k_{cl}} {\mathrm d} p
\frac{  V_0+V_{-2k_{cl}}-2V_{k_{cl}}} {2\epsilon_{k_{cl}+p}+C\mu/2S} \, .
       \end{aligned}$$ Noting that $$\begin{aligned}
&-& \frac{ 1}{4\pi S}  \dashint\limits_{-k_{cl}}^{\pi-k_{cl}} {\mathrm d} p \frac{  V_p+V_{-2k_{cl}-p}-(V_0+V_{-2k_{cl}}) } {2\epsilon_{k_{cl}+p}}\nonumber\\
&&=- \frac{ 1}{4\pi S}  \dashint\limits_{-\pi-k_{cl}}^{-k_{cl}} {\mathrm d} p
\frac{  V_p+V_{-2k_{cl}-p}-(V_0+V_{-2k_{cl}}) } {2\epsilon_{k_{cl}+p}} \quad \nonumber\\
&&=\frac{J^2-8}{16 S},
       \end{aligned}$$ and gathering all contributions, Eq. (\[gamma2\]) takes the form $$\begin{aligned}
\label{gamma21}
&&\!\!\!\!\!\!\!\!\tilde g(S)\! \left[1\!+\!\frac{J^2-8}{16S}\! +\! \frac{V_0\!+\!V_{-2k_{cl}}-2V_{k_{cl}} }{16\pi S}  \!\int\limits_{-\pi}^{\pi} \!\!\frac{ {\mathrm d} p } {\epsilon_{p}+C\mu/4S}   \right]
\nonumber\\
&&\!\!\!\!=V_0+V_{2k_{cl}}-2V_{k_{cl}}+\frac{8+J^2}{4 S} \nonumber\\
&&\!\!\!\!-\frac{1}{4\pi S}\dashint\limits_{-\pi}^{0}{\mathrm d}p\frac{V_{p-k_{cl}}+V_{p+k_{cl}}-2V_{k_{cl}}}{2\epsilon_{p}}(\tilde \Gamma_{p-k_{cl}}-\tilde\Gamma_{-2k_{cl}})\nonumber\\
&&\!\!\!\!-\frac{1}{4\pi S}\! \int\limits_{0}^{\pi}\!\! {\mathrm d}p\frac{V_{p-k_{cl}}+V_{p+k_{cl}}-2V_{k_{cl}}}{2\epsilon_{p}}(\tilde \Gamma_{p-k_{cl}}-\tilde\Gamma_{0}).\end{aligned}$$

Now we can plug the zeroth-order vertices into the right hand side of Eq. (\[gamma21\]), in the spirit of the $1/S$ expansion: $\tilde\Gamma_{p}\to
2\epsilon_p$. Noting that $V_0+V_{-2k_{cl}} -2V_{k_{cl}} = \tilde F+O(1/S)$, where $$\begin{aligned}
\label{tildeF}
\tilde F&=& V_0+V_{2k_{cl}}-2V_{k_{cl}}+\frac{8+J^2}{4 S} \\
&-&\frac{1}{2\pi S}\dashint\limits_{-\pi}^{0}{\mathrm d}p\frac{V_{p-k_{cl}}+V_{p+k_{cl}}-2V_{k_{cl}}}{2\epsilon_{p}}(\epsilon_{p-k_{cl}}-\epsilon_{-2k_{cl}})\nonumber\\
&-&\frac{1}{2\pi S}\int\limits_{0}^{\pi}{\mathrm d}p\frac{V_{p-k_{cl}}+V_{p+k_{cl}}-2V_{k_{cl}}}{2\epsilon_{p}}(\epsilon_{p-k_{cl}}-\epsilon_{0}),\nonumber\end{aligned}$$ in the same way as for $g(S)$ we obtain $$\label{renormalizedpotential_2}
\tilde g(S)= \frac{\tilde g_0(S)}{1+\tilde g_0(S){\sqrt{2m}}/({\pi \sqrt{\mu}})}$$ with $\tilde g_0(S)$ given in Eq. . We see that in order $1/S$, quantum fluctuations do not modify the relation $\tilde g_0(S)>g_0(S)$ for $J<0$. Therefore, the interspecies interaction is always positive (repulsive) for any $J$ and $S>1/2$, and behaves as $\tilde g_0(S)\sim (J+4)^2$ for $J\to -4$ for all $S$. A hard-core constraint, stemming from the mapping of spins to bosons, is not included in our theory. While it is well-known how to treat the hard-core constraint in an exact way for $S=1/2$,[@nikuni95] this is not the case for $S>1/2$.

DMRG results for periodic boundary conditions {#app:dmrg}
=============================================

The DMRG results presented in the main text were computed for open boundary conditions (OBC), except for those in Fig. \[fig:svn\] where we used periodic boundary conditions (PBC). Alternatively, one may compute all the data, [*e.g.*]{}, the magnetization curves $M(h)$ using PBC. This approach is generally expected to suffer from (i) slower convergence with respect to the number of states kept in the DMRG runs [@schollwoeck05] and (ii) large finite-size effects due to the incommensurability in the problem. We here demonstrate that all main features can be seen with both OBC and PBC, namely the existence of the metamagnetic jump and the Haldane gap. Moreover, the quantitative results are comparable, except for the expected finite-size effects due to the incommensurability.

![(Color online) Magnetization curves $M(h)$ for $S=1$ (a) $J=-3.5$ and (b) $J=-1$ ($L=48$), calculated with DMRG on systems with periodic boundary conditions. In panel (a), we clearly see the metamagnetic jump of height $\Delta M_{\mathrm{jump}}$ and in panel (b), the Haldane gap that defines $h_c$ shows up as a zero-field magnetization plateau.[]{data-label="fig:mag_curve_pbc"}](figure7.eps){width="0.9\columnwidth"}

Magnetization curves from periodic boundary conditions for $S=1$
----------------------------------------------------------------

First, we discuss some examples of magnetization curves obtained for PBC. We kept up to $m=1600$ states for the PBC DMRG computations presented here. In addition, we used exact diagonalization for those sectors which are sufficiently small. Therefore, in particular the data close to the saturation field are free of truncation errors.

Figure \[fig:mag\_curve\_pbc\] shows the magnetization curves for $J=-3.5$ \[panel (a)\] and $J=-1$ \[panel (b)\] for $L=48$ and PBC. As in the results for OBC, we observe the presence of the metamagnetic jump (here in the case of $J=-3.5$) and the Haldane gap (see the $J=-1$ curve), which manifests itself as a plateau in the magnetization curve at $M=0$. The latter defines the critical field $h_c$ that separates the gapped Double-Haldane phase from the gapless finite-$M$ phase (compare Fig. \[fig:svn\]).

Note that at intermediate $M$, the PBC data show spurious small steps with $\Delta S^z>1$. We have checked that these features disappear as one goes to larger system sizes.

![(Color online) Phase diagram of the frustrated ferromagnetic $S=1$ chain, comparing DMRG results from systems with OBC (open symbols) to results from systems with PBC (solid symbols). $L=64$ in all cases. Lines are guides to the eye.[]{data-label="fig:phase_diag_pbc"}](figure8.eps){width="0.9\columnwidth"}

Comparison of OBC vs PBC for $S=1$
----------------------------------

Figure \[fig:phase\_diag\_pbc\] contains the results for the Haldane gap that defines the critical field $h_c$ separating the gapped zero-field phase from the gapless phase at finite magnetizations, the saturation field $h_{\mathrm{sat}}$, and the jump height (inset), comparing data for OBC (open symbols) with data from PBC (solid symbols). Here we choose a chain length of $L=64$ for both PBC and OBC.

For PBC we can determine the spin gap from $$h_c^{\text{PBC}} = E(S^z=1) - E(S^z=0) \, ,
\label{eq:HcPBC}$$ where, as in Eq. , $E(S^z)$ is the ground-state energy in a sector with a given total $S^z$. The good agreement between the OBC and PBC results for $h_c$ in Fig. \[fig:phase\_diag\_pbc\] confirms that it is indeed appropriate to use Eq.  for OBC.

The PBC data for $\Delta M_{\mathrm{jump}}$ in Fig. \[fig:phase\_diag\_pbc\] suffer from the presence of several peaks, resulting in a non-monotonic dependence on $J$. This is due to the incommensurability in the finite-magnetization region, which is incompatible with the lattice vectors of a system with periodic boundary conditions that is translationally invariant. Therefore, in the main text we focus on the discussion of DMRG data from systems with OBC. The results from OBC and PBC data for the saturation field $h_{\mathrm{sat}}$, however, agree very well with each other.

[10]{}

T. Giamarchi, C. Rüegg, and O. Tchernyshyov, Nature Phys. [**4**]{}, 198 (2008).

E. M. Chudnovsky and L. Gunther, Phys. Rev. Lett. [**60**]{}, 661 (1988).

L. Balents, Nature [**464**]{}, 199 (2010).

K. Liu and M. Fisher, J. Low Temp. Phys. [**10**]{}, 655 (1973).

J. Struck, C. Ölschläger, R. Le Targat, P. Soltan-Panahi, A. Eckardt, M. Lewenstein, P. Windpassinger, K. Sengstock, Science [**333**]{}, 996 (2011).

Gyu-Boong Jo, J. Guzman, C. K. Thomas, P. Hosur, A. Vishwanath, and D. M. Stamper-Kurn, arXiv:1109.1591v1 (unpublished).

I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. [**80**]{}, 885 (2008); C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, , 1225 (2010).

J. M. Gerton, D. Strekalov, I. Prodan, and R. G. Hulet, Nature [**408**]{}, 692 (2000).

P. Nozi[è]{}res and D. Saint James, J. de Physique [**43**]{}, 1133 (1982).

E. J. Mueller and G. Baym, Phys. Rev. A [**62**]{}, 053605 (2000).

T. Vekua, A. Honecker, H.-J. Mikeska, and F. Heidrich-Meisner, Phys. Rev. B [**76**]{}, 174420 (2007).

L. Kecke, T. Momoi, and A. Furusaki, Phys. Rev. B [**76**]{}, 060407(R) (2007).

T. Hikihara, L. Kecke, T. Momoi, and A. Furusaki, Phys. Rev. B [**78**]{}, 144404 (2008).

J. Sudan, A. Lüscher, and A. M. Läuchli, Phys. Rev. B [**80**]{}, 140402(R) (2009).

F. Heidrich-Meisner, I. P. McCulloch, and A. K. Kolezhuk, Phys. Rev. B [**80**]{}, 144417 (2009).

H. P. Bader and R. Schilling, Phys. Rev. B [**19**]{}, 3556 (1979).

See, e.g., N. Shannon, T. Momoi, and P. Sindzingre, Phys. Rev. Lett. [**96**]{}, 027213 (2006); M. E. Zhitomirsky and H. Tsunetsugu, Europhys. Lett. [**92**]{}, 37001 (2010).

A. V. Chubukov, Phys. Rev. B [**44**]{}, 4693 (1991).

A. Kolezhuk and T. Vekua, Phys. Rev. B [**72**]{}, 094424 (2005).

F. Heidrich-Meisner, A. Honecker, and T. Vekua, Phys. Rev. B [**74**]{}, 020403(R) (2006).

M. Enderle, C. Mukherjee, B. F[å]{}k, 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. Rheinstädter, and H. M. R[ø]{}nnow, Europhys. Lett. [**70**]{}, 237 (2005).

M. Enderle, B. F[å]{}k, H.-J. Mikeska, R. K. Kremer, A. Prokofiev, and W. Assmus, Phys. Rev. Lett. [**104**]{}, 237207 (2010).

S.-L. Drechsler, O. Volkova, A. N. Vasiliev, N. Tristan, J. Richter, M. Schmitt, H. Rosner, J. Malek, R. Klingeler, A. A. Zvyagin, and B. Büchner, Phys. Rev. Lett. [**98**]{}, 077202 (2007).

F. Heidrich-Meisner, I. A. Sergienko, A. E. Feiguin, and E. R. Dagotto, Phys. Rev. B [**75**]{}, 064413 (2007); I. P. McCulloch, R. Kube, M. Kurz, A. Kleine, U. Schollwöck, and A. K. Kolezhuk, [*ibid.*]{} [**77**]{}, 094404 (2008).

É. G. Batyev and L. S. Braginskiĭ, Zh. Eksp. Teor. Fiz. [**87**]{} 1361 (1984) \[Sov. Phys. JETP [**60**]{}, 781 (1984)\].

T. Nikuni and H. Shiba, J. Phys. Soc. Jpn. [**64**]{}, 3471 (1995).

S. R. White, Phys. Rev. Lett. [**69**]{}, 2863 (1992).

U. Schollwöck, Rev. Mod. Phys. [**77**]{}, 259 (2005).

I. G. Gochev, Theor. Math. Phys. [**15**]{}, 402 (1974).

E. H. Lieb and W. Liniger, Phys. Rev. [**130**]{}, 1605 (1963).

K. Okunishi, Y. Hieida, and Y. Akutsu, Phys. Rev B. [**59**]{}, 6806 (1999).

M. D. Lee, S. A. Morgan, M. J. Davis, and K. Burnett, Phys. Rev. A [**65**]{}, 043617 (2002); M. D. Lee, S. A. Morgan, and K. Burnett, arXiv:cond-mat/0305416.

M. D. Johnson and M. Fowler, Phys. Rev. B [**34**]{}, 1728 (1986).

F. J. Dyson, Phys. Rev. [**102**]{}, 1217 (1956); [*ibid.*]{} [**102**]{}, 1230 (1956); S. V. Maleev, Zh. Eksp. Teor. Fiz. [**33**]{}, 1010 (1957) \[Sov. Phys. JETP [**6**]{}, 776 (1958)\].

T. Holstein and H. Primakoff, Phys. Rev. [**58**]{}, 1098 (1940).

E. B. Kolomeisky and J. P. Straley, Phys. Rev. B [**46**]{}, 11749 (1992).

A. K. Kolezhuk, Low Temp. Phys. [**36**]{}, 752 (2010); Phys. Rev. A [**81**]{}, 013601 (2010).

G. Jackeli and M. E. Zhitomirsky, Phys. Rev. Lett. [**93**]{}, 017201 (2004).

D. V. Dmitriev, V. Ya. Krivnov, and A. A. Ovchinnikov, JETP [**92**]{}, 146 (2001).

C. Gerhardt, K.-H. Mütter, and H. Kröger, Phys. Rev. B [**57**]{}, 11504 (1998).

S. Hirata, arXiv:cond-mat/9912066 (unpublished).

Part of the DMRG runs were done using ALPS 1.3, A. F. Albuquerque, F. Alet, P. Corboz, P. Dayal, A. Feiguin, S. Fuchs, L. Gamper, E. Gull, S. Gürtler, A. Honecker, R. Igarashi, M. Körner, A. Kozhevnikov, A. Läuchli, S. R. Manmana, M. Matsumoto, I. P. McCulloch, F. Michel, R. M. Noack, G. Paw[ł]{}owski, L. Pollet, T. Pruschke, U. Schollwöck, S. Todo, S. Trebst, M. Troyer, P. Werner, and S. Wessel, J. Mag. Mag. Mat. [**310**]{}, 1187 (2007).

See, e.g., S. R. White and D. A. Huse, Phys. Rev. B [**48**]{}, 3844 (1993), and references therein.

G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Phys. Rev. Lett. [**90**]{}, 227902 (2003).

P. Calabrese and J. Cardy, J. Stat. Mech. [**(2004)**]{} P06002.

A. Kolezhuk, R. Roth, and U. Schollwöck, Phys. Rev. Lett. [**77**]{}, 5142 (1996).

T. Kennedy, J. Phys.: Condens. Matter [**2**]{}, 5737 (1990).
