Multipolar phase in frustrated spin-1/2 and 1 chains

TL;DR

This paper investigates the phase diagram of the frustrated $J_1-J_2$ spin chain model, revealing the existence of higher order multipolar phases near the critical point and discussing experimental detection methods.

Contribution

It demonstrates the presence of higher order multipolar phases in the $J_1-J_2$ model and proposes criteria for experimental detection using inelastic neutron scattering data.

Findings

Existence of higher order $p>4$ multipolar phases near the critical point.

Modeling of INS data for LiCuVO$_4$ compound.

Analysis of spin nematic and SDW$_2$ phases in spin-1 systems.

Abstract

The $J_1-J_2$ spin chain model with nearest neighbor $J_1$ and next nearest neighbor anti-ferromagnetic $J_2$ interaction is one of the most popular frustrated magnetic models. This model system has been extensively studied theoretically and applied to explain the magnetic properties of the real low-dimensional materials. However, existence of different phases for the $J_1-J_2$ model in an axial magnetic field $h$ is either not understood or has been controversial. In this paper we show the existence of higher order $p>4$ multipolar phase near the critical point $(J_2/J_1)_c=-0.25$. The criterion to detect the quadrupolar or spin nematic (SN)/spin density wave of type two (SDW$_2$) phase using the inelastic neutron scattering (INS) experiment data is also discussed, and INS data of LiCuVO$_4$ compound is modelled. We discuss the dimerized and degenerate ground state in the quadrupolar phase. The major contribution of binding energy in the spin-1/2 system comes from the longitudinal component of the nearest neighbor bonds. We also study spin nematic/SDW$_2$ phase in spin-1 system in large $J_2/J_1$ limit.

Figures

Tables

Table 1. Table 1: The contribution of different bonds in the E b subscript 𝐸 𝑏 E_{b} are calculated for α = 1.0 𝛼 1.0 \alpha=1.0 and at M = 0.25 𝑀 0.25 M=0.25 and 0.4 for a chain of sizes N = 𝑁 absent N= 16, 20, 24 and 28 with PBC

$\alpha$	$System$	$E_b$	M = 0.25		M = 0.4
			$N\left(24\right)$	$N\left(28\right)$	$N\left(24\right)$	$N\left(28\right)$
1.0	Leg	$E_b^{L,L}$	-0.0408	-0.0713	-0.0332	-0.0237
		$E_b^{T,L}$	0.1935	0.2208	-0.2217	-0.0815
		Total	0.1526	0.1495	-0.2549	-0.1051
	Rung	$E_b^{L,R}$	-0.7103	-0.6600	-0.7150	-0.7493
		$E_b^{T,R}$	0.2070	0.1963	0.6597	0.5756
		Total	-0.5033	-0.4637	-0.0552	-0.1737
	Total binding energy		-0.3507	-0.3142	-0.3102	-0.2788

Table 2. Table 2: B p ​ b ​ c subscript 𝐵 𝑝 𝑏 𝑐 B_{pbc} for three values of α 𝛼 \alpha are shown for N = 24 𝑁 24 N=24 and 28. Eq. 11 is used to calculate the B p ​ b ​ c subscript 𝐵 𝑝 𝑏 𝑐 B_{pbc} in various S z superscript 𝑆 𝑧 S^{z} sectors.

$N$	$S^z$	$B_{pbc}$
		$\alpha =0.8$	$\alpha =1.0$	$\alpha =3.0$
24	1.00	0.10428	0.10820	0.12391
28	1.00	0.08997	0.09343	0.10763
24	3.00	0.10596	0.11176	0.13055
28	3.00	0.09158	0.09672	0.11390
24	5.00	0.09989	0.10933	0.13283
28	5.00	0.08831	0.09633	0.11786
24	7.00	-	0.09673	0.12705
28	7.00	-	0.08963	0.11706

Equations

$$H(J_1,J_2)=J_1\sum _n{S}_n\cdot S_{n+1}+J_2\sum _n{S}_n\cdot S_{n+2}-h\sum _n{S}_n^z,$$

$$\begin{array}{r}\theta =\frac{2\pi }{L},\end{array}$$

$$⟨S_0^{+}\cdot \cdot \cdot S_{p-1}^{+}{S}_r^{-}\cdot \cdot \cdot S_{r+p-1}^{-}⟩=(-1{)}^r⟨b_0{b}_r^{†}⟩$$

$$=\frac{A_m(-1{)}^r}{|r{|}^{1/\eta }}-\frac{{\tilde{A}}_m(-1{)}^r}{|r{|}^{\eta +1/\eta }}\cos (2\pi \rho r)+\cdot \cdot \cdot ,$$

$$\begin{split}C^{L}(r)=\langle S^{z}_{0}S^{z}_{r}\rangle&=\Bigg{\langle}\bigg{(}\frac{1}{2}-pb_{0}^{\dagger}b_{0}\bigg{)}\bigg{(}\frac{1}{2}-pb_{r}^{\dagger}b_{r}\bigg{)}\Bigg{\rangle}\\ &=M^{2}-\frac{p^{2}\eta}{4\pi^{2}r^{2}}+\frac{A_{z}\cos(2\pi\rho r)}{|r|^{\eta}}+\cdot\cdot\cdot,\end{split}$$

$$⟨S_r^z⟩$$

$$z(r;q)$$

$$q$$

$$f_{\nu }(x)$$

$${\rho }_q=⟨{\psi }_{n+2}|S_i^{+}{S}_j^{+}|{\psi }_n⟩,$$

$$E_{b}(n)=\frac{1}{2}\Big{[}E(n+2)+E(n)-2E(n+1)\Big{]},$$

$$E_{b}(n)=\frac{1}{2}\Big{[}E^{x,y}_{b}(n+2)+E^{x,y}_{b}(n)-2E^{x,y}_{b}(n+1)\Big{]}.$$

$$S^{\alpha \alpha }(q,\omega ,M)=\sum _n\frac{\mid ⟨{\psi }_n|S_q^{\alpha }|{\psi }_0⟩{\mid }^2}{E_n-(E_0+\omega )+i\eta },$$

$$|{\psi }_{gs}⟩$$

$$=\prod _{n=2L}({|1⟩}_{n,n\pm 1}+\eta {|-1⟩}_{n,n\pm 1})/\sqrt{1+{\eta }^2},$$

$$B_{pbc}=⟨{\psi }_{+}|(S_i\cdot S_{i+1}-S_{i+1}\cdot S_{i+2})|{\psi }_{-}⟩.$$

Full text

I Introduction

Interaction induced frustration and confinement of electrons in an one dimensional (1D) magnetic system generates many exotic phases Chubukov1991 ; Vekua2007 ; HikiharaMP2008 ; SudanLauchili . Some of these phases can have well defined order parameters, whereas other phases can have hidden order parameter. The 1D spin-1/2 systems with an isotropic $J_{1}-J_{2}$ model MajumdarGhosh ; TonegaHarada ; Chitra1995 ; WhiteAffleck ; ManuScaling2007 ; ManuModfDMRG ; Manu2010BOW ; Sandvik ; ManuDecouple ; ManuAslamJPCM ; ItoiQuin2001 ; Haldane1982 ; HamadaFQCP ; OkamotoNomura ; Bursill ; SirkerPRB ; Mahdavifar in the presence of an axial magnetic field have been extensively studied HikiharaMP2008 ; SudanLauchili ; Dimitriev2006 ; Meisner2006 ; Vekua2007 ; Kecke2007 ; MeisnerMcCulloch ; Furukawa ; SirkerPRB . The $J_{1}-J_{2}$ model in an axial magnetic field $h$ is written as

[TABLE]

where $J_{1}$ and $J_{2}$ are exchange interaction strengths between nearest neighbor (NN) and next nearest neighbor (NNN) spins, respectively.

The model with a ferromagnetic $J_{1}$ shows many interesting phases like spin liquid ManuAslamJPCM ; ItoiQuin2001 ; HamadaFQCP ; SirkerPRB ; Mahdavifar , dimer ManuAslamJPCM ; ItoiQuin2001 ; HamadaFQCP ; SirkerPRB ; Mahdavifar , chiral vector Chubukov1991 ; AslamJMMM , spin multipolar HikiharaMP2008 ; SudanLauchili , decoupled phase ManuAslamJPCM . The spin liquid phase is gapless and possesses quasi-long range order WhiteAffleck ; ItoiQuin2001 . The dimer phase is gapped in nature, and the spin-spin correlation decays exponentially Manu2010BOW ; WhiteAffleck ; ItoiQuin2001 . This model has been extensively used for modelling the magnetization properties of LiCuSbO4 SianManuINS2012 , LiCu2O2 ParkMultifer , (N2H5)CuCl3 N2H5CuCl3 , Rb2Cu2Mo3O12 Rb2Cu2Mo3O12 , Li2CuZrO4 Li2CuZrO4 , Ba3Cu3In4O12, and Ba3Cu3Sc4O12 BaCuIn ; BaCuSc . In the chiral vector phase, both spin parity and inversion symmetry are spontaneously broken ManuSpinParity . This phase has been studied extensively because of its potential application in improper multiferroic systems MathurMultifer ; KatsuraMultifer .

The field theoretical and numerical studies by Hikihara et al. suggest that metamagnetic or spin multipolar phase exist in the presence of the high axial magnetic field $h$ for ferromagnetic $J_{1}$ HikiharaMP2008 . These multipolar phases have hidden order parameters. In this model multipoles of order $p$ depend on the $J_{2}/J_{1}$ ratio HikiharaMP2008 ; SudanLauchili , and the nomenclature of each phase is done based on the number of bound magnons in the systems i.e., the number of paired magnons $p$ in dipolar, quadrupolar, octupolar and hexadecapolar phases are $1,2,3,4$, respectively. The quadrupolar phase is a Tomonaga-Luttinger liquid of hard core bosons HikiharaMP2008 , and each boson is made up of two magnons. In this phase the correlations between bosons and density fluctuations follow a power law. However, the boson propagator is dominant over the density fluctuations in this phase HikiharaMP2008 . In his seminal work Chubukov predicts that this phase has dimerized ground state (gs) Chubukov1991 , but Hikihara et al. show the absence of dimerization HikiharaMP2008 . In the large $J_{2}/J_{1}$ regime, field theoretical calculations show that the SDW2 phase exist in low magnetic field, whereas SN phase exists in the narrow range of magnetic field near the saturation field HikiharaMP2008 . The numerical calculations in $J_{2}/|J_{1}|=\alpha>0.6$ show the finite binding energy of magnon even for a small field AslamJMMM .

The order parameter of the SN phase $\langle S_{i}^{+}S_{j}^{+}\rangle$ is defined in ref. OlegBalents2014 ; Chubukov1991 ; Penc . It is hidden in nature, although the probes like the INS EnderleDynamic ; SianManuINS2012 and the resonant inelastic X-Ray scattering (RIXS) RIXS methods can indirectly measure these phases. The nematic phase in LiCuVO4 compound is confirmed by using the INS data of dynamical structure factor EnderleDynamic , and NMR data of this compound shows a sharp single and solitary line which moves with magnetic field SatoNMR1 ; SatoNMR2 . In this paper we try to show that there is characteristic feature of INS measurement for the SDW2 and SN phase.

In this model there are many unsettled issues such as, the metamagnetic phase in the small $J_{2}/J_{1}$ regime has been completely unexplored, and is difficult to characterize because of very small gaps. We have shown the gs degeneracies in the odd $S^{z}$ sectors AslamJMMM , but dimer order parameter $B$ is vanishingly small in this sector. The existence of quadrupolar phase in spin-1 systems is controversial, as steps of two in magnetization-$h$ curve is absent Arlego2011 ; Kolezhuk2012 , whereas the other studies for general spin show the existence of this phase. We explore this phase for the spin-1 system using the Hamiltonian in Eq. 1.

The rest of the paper goes in the following sequence. In section II the numerical techniques and accuracy of results are discussed. Results are discussed in section III. We start with the higher order multipolar phase and the relation between the pitch angle $\theta$ and magnetization $M$. The quadrupolar phase is discussed thereafter. The dynamical properties in quadrupolar phase of spin-1/2 $J_{1}-J_{2}$ model are discussed in subsection B. The dynamical properties of LiCuVO4 are also discussed in this subsection B. The dimer phase in the SN/SDW2 phase is presented in the subsection C. The results for spin-1 for the same model are discussed in the section IV. The discussion of all the results is done in the next section V.

II Numerical methods

The Density matrix renormalization group method (DMRG) is a state of art numerical technique to calculate accurate gs and a few low lying excited energy states of strongly interacting quantum systems WhiteDMRG ; KarenDMRG . It is based on systematic truncation of irrelevant degrees of freedom. We use modified DMRG algorithm, where four new sites are added to avoid the multiple time of renormalization of operators in the superblock. The modified DMRG has better convergence and also has sparse Hamiltonian matrix of superblock for the model Hamiltonian in Eq. 1, compared to the conventional DMRG where only one site is added in each block at every step ManuModfDMRG . The number of eigenvectors of the density matrix retained up to $m=400$ to maintain the truncation error of density matrix eigenvalues less than $10^{-10}$. In the worst case error in the energy is less than $0.01\%$. The DMRG is used for calculating various properties of large system sizes up to $N=368$ chain with open boundary condition (OBC). The number of finite DMRG sweeps required for an accurate gs and spin correlation function in the different $S^{z}$ sectors is approximately 20. Recently developed PBC algorithm is also employed for calculating the accurate gs and the correlation functions DayaPBC . The dynamical structure factor is calculated using the correction vector method RamseshaCorrection ; JeckelmannCorrection ; OnishiJapan .

III Results

The quantum phase diagram of $J_{1}-J_{2}$ model in an axial magnetic field given in Eq. (1) consists of numerous phases such as the vector chiral (VC) Chubukov1991 ; AslamJMMM , the dimer MajumdarGhosh ; TonegaHarada ; Chitra1995 ; WhiteAffleck ; ManuScaling2007 ; ManuModfDMRG ; Manu2010BOW ; ManuDecouple ; ManuAslamJPCM ; ItoiQuin2001 ; Haldane1982 ; HamadaFQCP ; OkamotoNomura ; Bursill ; SirkerPRB ; Mahdavifar , the decoupled chain ManuDecouple ; ManuAslamJPCM , and multipolar/SDWn phases HikiharaMP2008 ; SudanLauchili . In this paper, the SN/SDW2 phase and other higher order multipolar phases are discussed. This section is divided into three subsections. In subsection A, multipolar phases for spin-1/2 are discussed in the beginning; SN/SDW2 phase is presented in later part of subsection A. The general observations about dynamical property and $M-h$ curve in quadrupolar phase are presented. We model the dynamical structure factor of LiCuVO4 and also compare our results with the experimental data available in literature EnderleDynamic ; MourigalFieldMom in subsection B. The dimer in SN/SDW2 phase is presented in subsection C.

III.1 Multipolar phases in $S=1/2$

The multipolar phase and the spin density wave in the $J_{1}-J_{2}$ model for spin-1/2 chain in the presence of magnetic field $h$ are discussed in this part. We notice that there is a level crossing from ferromagnetic to singlet gs at $\alpha_{c}=0.25$ HamadaFQCP , and near to the critical point $\alpha_{c}$, but $\alpha>0.25$ limit, multiple magnons bind to form multipoles below the saturation magnetic field. It is also noted that number of $p$ changes rapidly with $\alpha$. In this paper, multipolar phase with order $p$ is explored based on the magnetic steps, pitch angle $\theta$ of spin density, and spin correlations in the gs at a finite magnetic field $h$. The angle between two nearest neighbor spin is called pitch angle $\theta$ and is defined as

[TABLE]

where $L$ is the smallest distance between spins whose pitch angle differs by $2\pi$. The field theoretical bosonization calculations HamadaFQCP suggest that for $\alpha>0.25$, the system shows SDWn in low magnetic field, whereas it shows multipolar phase at high magnetic field HikiharaMP2008 . The multipolar correlation of order $p$ or boson propagator $\langle S_{i}^{+}S^{+}_{i+1}\cdot\cdot\cdot S^{+}_{p+i-1}S^{-}_{i+r+1}S^{-}_{i+r+2}\cdot\cdot\cdot S^{-}_{i+r+p-1}\rangle$ is written HikiharaMP2008 as

[TABLE]

where $A_{m}$ and $\widetilde{A}_{m}$ are constants, $\eta$ is twice of the Luttinger liquid parameter, and $r$ represents distance. The density-density correlation is written as

[TABLE]

where $\rho=\frac{1}{p}(1-\frac{M}{M_{0}})$, $M_{0}$ is the saturation magnetization. The pitch angle $\theta=2\pi\rho$ varies with the magnetic field. The spin density $\langle S_{r}^{z}\rangle$ calculated from the field theoretical method is written HikiharaMP2008 as

[TABLE]

The pitch angle $\theta_{T}$ in transverse direction can also be extracted from the transverse correlation function $C^{T}(r)=\langle(S^{x}_{i}S^{x}_{i+r}+S^{y}_{i}S^{y}_{i+r})\rangle$. However, the pitch angle $\theta$ in the longitudinal direction is calculated from $C^{L}(r)$ and spin density $\langle S^{z}_{r}\rangle$.

The $\langle S^{z}_{r}\rangle$ and $C^{L}(r)$ are shown in Fig. 1 (a) for $M=0.35,\alpha=1.0$ and $N=168$. The $C^{L}(r)$ is scaled by 3.25 and $r$ is shifted by 1.0 unit to match the magnitude of $\langle S^{z}_{r}\rangle$. Interestingly, the complex looking equation of $\langle S^{z}_{r}\rangle$ in Eq. (5) has similar variation as that of $C^{L}(r)$. All the $\langle S^{z}_{r}\rangle$ and $C^{L}(r)$ give the same pitch angle. The Friedel oscillation at the edge of chain is seen in both $\langle S^{z}_{r}\rangle$ and $C^{L}(r)$. The spin densities are plotted in Fig. 1 (b) for $M=0.05,0.1,0.25$ and $0.3$. The amplitude of $\langle S^{z}_{0}S^{z}_{r}\rangle$ decreases with the distance, whereas $|\langle S^{z}_{r}\rangle|$ at site $r$ is more or less constant with $r$. Therefore, it is easier to calculate $\theta$ from $\langle S^{z}_{r}\rangle$ than from $\langle S^{z}_{0}S^{z}_{r}\rangle$. We find that $\theta$ decreases with $M$ and reduces to zero at $M=0.5$ for spin-1/2 system. These results are consistent with the Sudan et al. exact diagonalization results SudanLauchili .

As shown in the Fig. 1 (a) $\theta$ calculated from $\langle S^{z}_{r}\rangle$ and $C^{L}(r)$ are same, and it follows a linear relation with $M$. With $M$ the variations in $\theta_{T}$ of the transverse correlation functions $C^{T}(r)$ is less than 5$\%$. The accurate calculation of $\theta$ near $M_{0}\approx 0.5$ requires larger system size, and for these calculations we have used $N=168$ for low magnetization and 368 for higher magnetization. The $\theta$ and $\theta_{T}$ are calculated from the $\langle S^{z}_{r}\rangle$ and the $C^{T}(r)$, respectively, for $\alpha=0.265,0.27,0.3,0.4,$ and 1.0 as a function of $M/M_{0}$ shown in Fig. 2 (a) and 2 (b). The filled symbols are the DMRG calculations for $N=168$ with OBC, and dotted lines are fitted line with $\frac{\theta}{\pi}=\frac{1}{p}(1-\frac{M}{M_{0}})$ where $p$ is the order of the multipole.

In Fig. 2 (a) we notice that the variation of the $\theta$ with $M$ shows linear relation $\frac{\theta}{\pi}=\frac{1}{p}(1-\frac{M}{M_{0}})$ especially at large $M$. For $\alpha>0.4$ and large $M$, $\theta$ varies linearly with $M/M_{0}$ with a slope $2/p=1$. The linear behavior of $\theta$ deviates from straight line at low $M/M_{0}$ for $\alpha\leqslant 0.6$. The deviation point for $\alpha=0.4$ is at $M/M_{0}\approx 0.28$. In the VC phase $\theta$ depends weakly on $M$ as shown in Fig. 2 (a). The phase boundary of the quadrupolar and the VC phase is estimated using the level crossing or magnetic step criterion as in ref. HikiharaMP2008 ; AslamJMMM . For $\alpha\geq 0.4$ results will be discussed in the later part of this section.

The three magnon bound phase or the triatic/SDW3 phase occurs in the vicinity of $\alpha=0.3$ and $\frac{\theta}{\pi}$ is less than $0.26$ at $M<0.21$ as shown in Fig 2 (a). At large $M$ the slope of green line in Fig. 2 (a) is $1/p=1/3$. The phase boundary of the triatic/SDW3 and the VC phase can also be estimated from the deviation of the $\frac{\theta}{\pi}$ from linear relation as shown in Fig. 2 (a). In fact $\theta$ weakly depends on $M$ in the VC phase and remains constant for the given value of $\alpha$, whereas it varies linearly with slope $1/p=1/3$ in the triatic/SDW3 phase. The phase boundary of the triatic/SDW3 and the VC phase calculated with this method is consistent with other calculations SudanLauchili ; HikiharaMP2008 . The maximum value of $\theta$ for a multipole of order $p$ for a given $\alpha$ is $\pi/p$, and it decreases with the number of magnons or $p$. Our DMRG result shows that for $\alpha=0.265$ at large $M$, $p=5$ state shows up for $M/M_{0}>0.4$, and system show the $p=4$ state for intermediate magnetization $0.12<M/M_{0}<0.4$. The vector chiral phase sets in below the $M/M_{0}\leq 0.12$. For $\alpha<0.265$, $\theta$ calculations become difficult with the approximate numerical technique. In this limit energy states are closely spaced, and accurate determination of wavefunctions of the closely spaced energy levels is difficult.

In the quadrupolar phase the binding energy $E_{b}$ of the magnons defined in Eq. 7 below is an important quantity to understand the condensation phenomenon. The $M-h$ curve is analyzed to see the effect of condensation of magnons. Near the critical point $\alpha_{c}=0.25$, energy level spacings are tiny. Therefore, to maintain the accuracy of results, the ED method is used to solve the Hamiltonian for systems with $N=16,20,24$ and 28. In Fig. 3 (a), the finite size effect on the $M-h$ curve is shown for $\alpha=0.254$. The $N=16$ shows jumps of 1 and 7, whereas $N=20$ shows steps in $M$ of size 1 and 8. The gaps between the energy levels decrease with system size $N$, and for $N=24$ and 28 system shows steps of $N/2$. For $\alpha=0.254$ the value of $p$ can be equal or higher than $N/2$. The finite system size effect on the gaps is weak in this parameter regime. The $M-h$ curve for different values of $\alpha=0.254,0,256,0.258$ and 0.26 are shown in Fig. 3 (b) for $N=28$. We notice that the $h$ required for saturation increases with $\alpha$. The step size depends on $\alpha$ for example, at $\alpha=0.256$ and 0.258 system shows jumps of 1 and 12, whereas for $\alpha=0.26$ the jumps are 1, 2 and 6. The VC phase exists in the low $M$ limit, and the phase boundary decrease with $\alpha$. The magnetic steps or the order of multipole $p$ increases rapidly with $1/\alpha$ near the critical point $0.25$, and the magnetic gaps decrease with $\alpha$ as shown in Fig. 3 (a). Unfortunately we need large system size to confirm the large $p>14$, but these results are consistent with the prediction of existence of larger $p$ in ref. OlegBalents2016 .

In case of incommensurate spin density wave there are level crossings or the gs degeneracies, and these two degenerate states have opposite inversion symmetry AslamJMMM . Our ED calculations show that the gs energies are degenerate at large magnetization for $\alpha<0.4$ for both odd and even $S^{z}$ sectors. We calculate the z-component of VC order parameter $\kappa_{i}^{z}=\frac{1}{N}\sum_{i}\bra{\psi_{+}}(S_{i}^{+}S_{i+1}^{-}-S_{i}^{-}S_{i+1}^{+})\ket{\psi_{-}}$ AslamJMMM . In the multipolar phase, the $\kappa_{i}^{z}$ at large $S^{z}$ limit is non-zero for $0.25<\alpha<0.55$ and system sizes up to $N=28$. The $\kappa_{i}^{z}$ for $N=24$ system size is shown for different $M$ for $0.25<\alpha<1.0$ in ref. AslamJMMM .

In the large $\alpha$ limit, the SDW2 and SN phase exist in the presence of the magnetic field $h$. In the SN phase two magnons can condense to form a single boson HikiharaMP2008 ; AslamJMMM , and this phase is determined based on the presence of magnetic step of two $(\Delta S^{z}=2)$ in $M-h$ curve HikiharaMP2008 ; AslamJMMM , order parameter and various correlation functions. The order parameter for the quadrupolar phase is defined as in ref OlegBalents2014 ; Penc ; Chubukov1991

[TABLE]

where $\ket{\psi_{n}}$ and $\ket{\psi_{n+2}}$ are gs of $S^{z}=n$ and $n+2$ spin sector, but both of these are degenerate in the presence of an applied magnetic field.

In this phase the variation of $\theta$, magnitude of $E_{b}$, $\rho_{q}$ and the dynamical structure factor $S(q,\omega)$ as a function of $M$ are calculated in the presence of magnetic field. The variation of magnetization $M$ with $h$ at $\alpha=0.6$ is shown in Fig. 4 (a) for chains with $N=104$ and 168. The magnetic step of $\Delta S^{z}=2$ exists in the full range of $M$. The existing literature shows the SDW2 phase at low magnetic field and SN type at high magnetic field HikiharaMP2008 ; OlegBalents2014 ; OlegBalents2016 . To analyze the quadrupolar phase, $\theta$ is plotted as a function of $M/M_{0}$ in Fig. 4 (b) for two different $N=104$ and 168 chain. The dashed line indicate $\frac{\theta}{\pi}=\frac{1}{p}(1-\frac{M}{M_{0}})$ line with $p=2$. These calculations demonstrate weak size dependence of the pitch angle $\theta$.

The average binding energy of two magnons is defined as

[TABLE]

where $E(n)$ is the energy of the system with even number of magnons $n$. The binding energy of two magnons in the SN/SDW2 phase is shown as a function of $N$ in the inset of Fig. 5. The $|E_{b}|$ has weak finite size dependence in large $M$ limit, whereas it shows significant change with system size in low field limit. The finite size scalings are done for $M=0.05,0.1$ and $M=0.4$ at $\alpha=1.0$ for $N$ up to 200. $|E_{b}|$ increases with $M$ and has finite extrapolated value for $M>0.1$. However, $|E_{b}|$ at low magnetization $M=0.05$ is vanishingly small.

In Fig. 5 the extrapolated values of $|E_{b}|$ as a function of $\alpha$ for different $M=0.05,0.1,0.15,0.2,0.25$ and 0.45 are shown. The error bars reflect the error in extrapolation and inaccuracy in DMRG calculations. We notice that $|E_{b}|$ increases with $\alpha$ and it attains a maximum value around $\alpha_{m}(M)$ for a given $M$, and decreases thereafter. The value of $\alpha_{m}(M)$ increases with $M$. The $|E_{b}|$ increases with $M$ initially and either it saturates or decreases near the saturation magnetic field. This trend of $|E_{b}|$ is consistent with the calculations done by Onishi OnishiJapan .

The bond energies are analyzed to understand the contribution of different bonds in the $E_{b}$. In the large $\alpha$ limit the $J_{1}-J_{2}$ model for a chain behaves like a zigzag chain, and the next nearest neighbor interaction $J_{2}$ of the model act as interaction between the spins along the leg, whereas the nearest neighbor interaction $J_{1}$ becomes the interaction along the rung ManuModfDMRG . The contribution of different bonds in the $E_{b}$ are calculated for $\alpha=1.0$ and at $M=0.25$ and 0.4 for a chain of sizes $N=$16, 20, 24 and 28 with PBC. However, data are shown only for $N=24$ and 28 in the table 1. The binding energy contribution of different bonds $E^{x,y}_{b}$ where $x$ stands for longitudinal $(L)$ or transverse $(T)$ and $y$ stands for leg $(L)$ or rung $(R)$. The $E_{b}$ is defined in terms of $E_{b}^{xy}$ as,

[TABLE]

For $M=0.25$ the major contribution to $E_{b}$ are transverse component $E^{T,R}_{b}$, $E^{T,L}_{b}$ and $E^{L,R}$ ,however, transverse component weakens the $E_{b}$ as shown in table 1. The $E^{L,L}_{b}$ decreases with system size, whereas $E^{T,L}_{b}$ increases with the system size. The magnitude of $E^{T,R}_{b}$ is significantly smaller than the $E^{L,R}_{b}$. The major contribution of $E_{b}$ comes from the $E^{L,R}_{b}$. The magnitude of $E^{T,R}_{b}$ is almost 1/3 of the $E^{L,R}_{b}$, but these two have opposite signs. However, both these quantities increase with $M$. For $M=0.4$ both along leg and rung transverse bonds contributions weaken the total $E_{b}$. The $E^{T,L}_{b}$ also decreases, whereas $E^{T,R}_{b}$ increases. The magnitude of $E^{L,R}_{b}$ and $E^{T,R}_{b}$ are very similar, but opposite to each-other for $M=0.4$. In conclusion rung contributes most of the $E_{b}$ in small $M$, but contribution of leg increases with $M$ of $M<M_{0}$. The $E_{b}$ is still small, however, $E^{T,L}_{b}$ is significantly large.

The quadrupolar phase is directly quantified in terms of the order parameter $\rho_{q}$ defined in Eq. 6. In the inset of Fig. 6 $\rho_{q}$ is extrapolated as a function of $1/N$ for $M=0.1$, 0.2 and 0.4. All the curves of $\rho_{q}$ follows the linear behavior with $1/N$. The extrapolation of $\rho_{q}$ is done with system size upto $N=104$. The extrapolated values of $\rho_{q}$ are shown in the main Fig. 6 for $M=0.1$, 0.2, 0.3, 0.4 for $0.5\geq\alpha\geq 1.0$. The value of the $\rho_{q}$ is within the error limit for $M\leq 0.1$ for smaller $\alpha$. The $\rho_{q}$ increase with $M$, and varies slowly with $\alpha$. The behavior of $\rho_{q}$ is quite consistent with $E_{b}$ as both of these quantity increase with $M$.

III.2 Dynamical structure factor

The dynamical structure factor MullerDynamical is defined as

[TABLE]

where $\ket{\psi_{0}}$ and $\ket{\psi_{n}}$ are the gs wavefunction for fixed $S^{z}=M$ and nth excited states for same $M$ or $M\pm 1$, respectively. $S_{q}^{\alpha}$ is defined as, $S_{q}^{\alpha}=(\sqrt{2\pi/{N}})\sum_{j}S_{j}^{\alpha}e^{iqj}$, where $\alpha=x,y$ and z component. $E_{0}$ and $E_{n}$ are the gs and $n^{th}$ excited state energies, respectively, $\omega$ is the energy transferred to the spin lattice. $\eta$ is broadening factor and is fixed at 0.1 for all the calculations.

The dynamical structure factor $S^{zz}(q_{m},\omega)$ is shown in Fig. 7 for $\alpha=1.0$. The $S^{zz}(q_{m},\omega)$ for the system of size $N=104$ with given $M$ represents structure factor for a given value of momentum $q_{m}$ for which the intensity is the highest. The $S^{zz}(q_{m},\omega)$ for $q_{m}$ is shown as a function of $\omega$ for different $M=0.0,0.1,0.2,0.3,0.4$ and 0.45. As $M$ increases the peak position of $S^{zz}(q_{m},\omega)$ shifts towards lower $q_{m}$ and $\omega_{m}$. However, the longitudinal spin excitation is gapless in the SN/SDW2 in the thermodynamic limit. For $M=0.0$, $q_{m}/\pi$ is at 0.5 and $q_{m}$ decreases with increasing $M$. In the inset of Fig. 7 open circle represents the $q_{m}$ for different values of $M$. The calculated $q_{m}$ is fitted with a function $q_{m}/\pi=(1-{M}/{M_{0}})/2$. These features of SN/SDW2 phase is directly examined by inelastic neutron scattering experiment in the presence of magnetic field $h$.

The existence of SN phase in a real material like LiCuVO4 is confirmed in the presence of high magnetic field $h$. This material consists of planar arrays of spin-1/2 copper chains with a ferromagnetic nearest neighbor $J_{1}$ and antiferromagnetic next nearest neighbor exchange interactions $J_{2}$. The exchange interaction strengths are $J_{1}=-1.6$ meV and $J_{2}=3.8$ meV, found by fitting the data of INS and other experiments EnderleLiCuVO4 .

We use these parameter values for our calculations. The dynamical structure factor $S^{zz}(q,\omega)$ in the absence of the magnetic field is shown in Fig. 8. The intensity is shown by the contour plots. The experimentally observed $S(q,\omega)$ in figure 2 of ref. EnderleDynamic shows as a function of $q$ and $\omega>3$ meV. The random phase approximation (RPA) calculation shows continuous intensity below $\omega<5$ meV, whereas experimental data shows high intensity between $\omega=3$ to 5 meV with momentum between $q/\pi=0.2$ and 0.5. The experimental data is restricted to $\omega\geq 3$ meV and shows only higher level of excitations. For better resolution of intensity we plot the logarithm of $S(q,\omega)$ intensity in Fig. 8. Our DMRG calculations shows that the most intense peak is at $q_{m}=\pi/2$ and $\omega=0.3$ meV. In fact, there are several values of $q$ and $\omega<3$ meV at which this system shows significant intensity of $S^{zz}(q,\omega)$. The $S(q,\omega)$ follows the the sum rule and we notice that major part of the intensity sum is limited to smaller $\omega$, and intensities of $S(q,\omega)$ observed experimentally are only a small fraction of the total intensity. Actually, this is easily justified by the slow variation of intensity for $\omega>3$ meV in experimentally observed $S(q,\omega)$.

The binding energy $|E_{b}|$ and momentum $q_{m}$ in the presence of magnetic field are two important quantities to characterize the SN phase. The INS experiment on LiCuVO4 by Mourigal et al. in ref. MourigalFieldMom shows the linear variation of momentum $q$ with magnetic field in high magnetic field $h$ limit. However, $q$ is independent of field below $h=8T$. Our results for the $|E_{b}|$ and momentum $q_{m}$ are shown in Fig. 9 (a) and (b) for two system sizes $N=104$ and 168. The $q_{m}$ for LiCuVO4 as a function of magnetization is shown in Fig. 9 (b) for $T$=0 K. We notice that $q_{m}$ follows the linear relation with $M$ with a slope of $1/p=0.5$. The linear dependence of momentum is followed in the full range of $M$. The $|E_{b}|$ is shown in Fig. 9 (a) as a function of $M/M_{0}$. For $M=0$, $|E_{b}|$ vanishes and it increases with $M$ up to $M\approx 0.4$ and then remains constant and it increases thereafter.

III.3 Dimers in spin nematic phase

In the paper by Chubukov, he suggested the existence of dimerized uniaxial SN phase which is different from the conventional dimerization where the two nearest spin form singlet pair Chubukov1991 . In this type of dimerization state two neighboring spin forms spin $S=1$ state. The gs wave function is written as

[TABLE]

where $\ket{1}$ and $\ket{-1}$ are $\ket{\uparrow\uparrow}$ and $\ket{\downarrow\downarrow}$ triplet states, respectively. Although, bosonization calculation by Hikihara et al. suggests that dimerization is proportional to $\cos(a{\phi_{-}}+\pi M)$ and their average vanishes to zero HikiharaMP2008 .

We notice that gs is doubly degenerate in odd $S^{z}$ in a finite system with PBC for $\alpha>0.5$ AslamJMMM . These two degenerate gs have opposite inversion symmetry. We notice that these degeneracies are independent of system size. In large $J_{2}$ limit this system is mapped to a zigzag chain with leg A and B. In the odd $S^{z}$ sectors the difference between the total spin densities on each leg A and B differ by 1. Therefore, the extra magnon is confined to either leg A or B depending on the symmetry of the system AslamJMMM . Now the broken symmetry state is defined as $\ket{\psi_{\pm}}=\frac{1}{\sqrt{2}}(\ket{\phi_{+}}\pm\ket{\phi_{-}})$, where $\ket{\phi_{+}}$ and $\ket{\phi_{-}}$ are degenerate states with $+$ and $-$ inversion symmetry. Dimer order parameter $B_{pbc}$ for periodic system is defined ManuModfDMRG as

[TABLE]

The values of $B_{pbc}$ along the rung for odd $S^{z}$ sectors defined in Eq. (11) for PBC systems are listed in table 2 for $N=24$ and 28. We notice that $B_{pbc}$ is approximately constant $0.1\pm 0.02$ with $M$ and decreases with $N$ as shown in the table 2. All the values of $B_{pbc}$ along the leg are zero.

We have also calculated the dimer order parameter $B$ in OBC system in even $S^{z}$ sector. In this sector gs is non degenerate and show spiral behavior. We have followed standard procedure to calculate $B$ in ref. WhiteAffleck ; ManuScaling2007 . The $B$ shows non-monotonic behavior with system size and it is small for large system.

IV Quadrupolar phase in spin-1

In this section, we explore the SN/SDW2 or quadrupolar phase for spin $S=1$ for finite system size with PBC, and assuming the spins interaction follows the model Hamiltonian in Eq. (1).

For the model in Eq. 1 the ferromagnetic to singlet crossover occurs at $\alpha=0.25$ and the singlet state extends for all values of $J_{2}>0.25$. The singlet and the triplet excitation or spin gap near the critical point $\alpha\approx 0.25$ is small compare to the spin gap in anti-ferromagnetic $J_{1}$ model. However, the double Haldane gap is observed in large $\alpha$ limit. The multipolar phase of higher order $p>2$ is observed for $\alpha<0.5$ which is consistent with earlier studies. The gs have spiral arrangement of the spins for $\alpha>0.25$. A detailed study of these properties of the system will be presented somewhere else Spin1New . In this section SN or SDW2 phase is explored for spin-1 chain with PBC in the large $\alpha$ limit. We notice that the energy convergence in DMRG calculation depends on the number of relevant degrees of freedom $m$ kept in the calculation, and energy of odd and even $S^{z}$ sectors follow the linear relation with $m$ but with different slopes. Therefore, we limit our calculations only to ED upto $N=16$.

The $M-h$ plot for three $\alpha=0.97,0.98$ and 0.99 for $N=16$ is shown in Fig. 10. The magnetic steps $\Delta S^{z}$ in chain with OBC is one, however, in PBC chain it is two. This may be because of the edge modes at the end of the chain in OBC case. We notice that there are elementary steps of $\Delta S^{z}=2$ in the magnetization with the magnetic field. The transition of steps $\Delta S^{z}=1$ to $\Delta S^{z}=2$ occurs at high magnetic field, and as $\alpha$ value increases the crossover point shifts to higher magnetic field. We also notice that for $\alpha>1$, all the elementary steps are $\Delta S^{z}=2$. In the main Fig. 10, variation of $M$ with $h$ is shown. The transition from mixed steps of $\Delta S^{z}=1$ and $\Delta S^{z}=2$ to purely $\Delta S^{z}=2$ step occurs at $\alpha=0.98$ for $N=16$ for different $M$ as shown in the main Fig. 10. To see the finite size effect, $M-h$ curve is plotted for $N=8,12$ and 16 for $\alpha=1.0$. We notice that the magnetic gaps decrease with $N$.

V Discussion

In this paper frustrated $J_{1}-J_{2}$ model Hamiltonian in Eq. 1 for spin-1/2 and 1 chains is studied. Our studies are focused on the model with ferromagnetic NN and antiferromagnetic NNN interactions in the presence of a magnetic field $h$. We use the ED and the DMRG numerical techniques to solve the Hamiltonian in Eq. 1. Here we have discussed multipolar phases, and especially, focused on SN phase of this model. The pitch angle $\theta$, the binding energy $E_{b}$, the order parameter $\rho_{q}$ and the steps in the magnetization are used to characterize the SN phase. We modelled the dynamical structure factor $S(q,\omega)$ of the LiCuVO4 compound using the parameter values $J_{1}=-1.6$ meV and $J_{2}=3.8$ meV in the literature EnderleLiCuVO4 . The quadrupolar phase in the spin-1 chain is also discussed in the large $\alpha$ limit.

The multipolar phase is characterized based on the pitch angle $\theta$ calculated from spin density and correlation function. We show that spin density and longitudinal spin-spin correlations are commensurate with each other as shown in Fig. 1. The pitch angle $\theta$ vs. magnetization $M$ plot shows multipolar phase of order up to $p=5$ at $\alpha=0.265$, however, the previous calculations by Sudan et al. are restricted to $p=4$ and all calculations were limited to system size up to $N=28$ SudanLauchili . In this paper the DMRG calculations are done for system size up to $N=368$, especially in the large magnetization limit. We notice that in $M\rightarrow 0$ limit $\theta$ is weakly dependent on $M$, although, in the large magnetization limit, the pitch angle $\theta$ shows a linear behavior in the multipolar phase for $\alpha<0.60$. This result is consistent with the calculations of Sudan et al. SudanLauchili . The junction of flat regime and linear variation of pitch angle $\theta$ is good estimate of the VC and the multipolar phase boundary. The variation of $\theta_{T}$ calculated from transverse correlation is almost independent of magnetization, and it is explained in terms of the finite gap and exponentially decaying correlation function HikiharaMP2008 .

The characterization of multipolar phase of order $p>5$ with approximate numerical technique is a difficult task because of the presence of large number of nearly degenerate states, and in this case it is difficult to get pure gs without using symmetry. To avoid the accuracy problem the ED is used to calculate step in $M-h$ curve. After careful investigation of gaps we show that the multipolar phase of order $p=12$ at $\alpha=0.256$, and $p=N/2$ for $\alpha<0.254$ for $N\geq 24$. Although some of the previous works show that these are metamagnetic phases HikiharaMP2008 ; SudanLauchili , we find these are actually higher order multipolar phases with small binding energy.

The binding energy $|E_{b}|$ in SN/SDW2 phase rapidly increases with $\alpha$ initially, and it has maxima at $\alpha_{m}(M)$. In the large $\alpha$ limit the bond energy contribution of the rung decreases with $\alpha$, therefore $|E_{b}|$ decreases with $\alpha$. The value of $\alpha_{m}(M)$ increases with $M$, and the $|E_{b}|$ have a broad maxima as a function of $\alpha$ as shown in the Fig. 5. The bond energy analysis is done in Table 1. For lower $M$, transverse bond energy for legs and rungs both have contribution to $E_{b}$, whereas longitudinal contribution of rung plays major role in binding of two magnons. The contributions of legs and rungs for higher $M$ have similar trend except that the magnitude of longitudinal contribution decreases in leg, and it increases in rung. The $E^{T,L}_{b}$ decreases for higher $M$, whereas $E^{T,R}_{b}$ increases significantly. The $E^{T,R}_{b}$ actually weakens the binding of the magnons and longitudinal component try to enhance the $E_{b}$. The values of $E_{b}$ for $\alpha=0.5,1,1.5$ and 2 have similar value to the previous calculation by Onishi OnishiJapan .

The earlier studies of $J_{1}-J_{2}$ model of general $S$ chain show the absence of spin nematic phase in $S=1$ chain Arlego2011 ; Kolezhuk2012 . However, the study of Balents et al. shows presence of nematic phase in general spin using the Lifshitz nonlinear sigma model OlegBalents2016 . Our finite size calculations at $\alpha>0.98$ show steps of 2 in $M-h$ curve and this results can be understood using their model OlegBalents2016 . The double Haldane phase agree with ref. Arlego2011 ; Kolezhuk2012 . However, we note doubly degenerate gs in odd $S^{z}$ sectors.

To characterize the SN/SDW2 or quadrupolar phase the dynamical structure factor $S(q,\omega)$ is analyzed, and we notice that the momentum $q_{m}$ of most intense peak of $S(q,\omega)$ for a given $M$ varies linearly with $M$. This result can be directly confirmed by the INS experiments. The LiCuVO4 is the most studied material for SN/SDW2 phase, and we calculate the $S(q,\omega)$ in the absence of magnetic field $h$. The high energy peak is consistent with the earlier results, but the most intense peak is at $q=\pi/2$ and $\omega=0.4$ meV. We also predict the dependence of $q_{m}$ as a function of $M$ and the $M-h$ curve for single crystal of this compound. The linear variation of $q$ with magnetic field $h$ is shown for LiCuVO4 at high magnetic field by Mourigal et al. MourigalFieldMom . However, a more accurate measurement should be performed at low field to verify the theoretical predictions. For this material our calculation shows the linear variation of $q_{m}$ with $M$ for longitudinal $S^{zz}(q,\omega)$ for the given parameter in ref. EnderleLiCuVO4 .

In this paper the order parameter $\rho_{q}$ is also calculated for finite $N$. We notice that the extrapolation of $\rho_{q}$ goes to zero for $M\rightarrow 0$, but it is finite at high value of $M$. This result seems to partially agree with calculation of the exponent $\eta$ with $M$ HikiharaMP2008 . We find difficulties in calculating $\eta$ from spin density and spin correlations. Our order parameter calculation shows that existence of SN phase much below the saturation magnetic field. The bond dimerization in the SN phase at finite $h$ has been under debate. Chubukov claims that there are $S=1$ dimerization and the doubly degenerate gs Chubukov1991 , but the analytical calculation by Hikihara et al. HikiharaMP2008 shows the absence of dimerization. We show that the even $S^{z}$ have non degenerate and spiral gs. The odd $S^{z}$ have doubly degenerate gs for PBC system.

In conclusion, we have studied the $J_{1}-J_{2}$ model in an axial magnetic field $h$ with ferromagnetic $J_{1}$. The multipolar phases with multipole up to $p=14$ are calculated. We have analyzed $E_{b}$ in SN/SDW2 phase, and we show that longitudinal energies of rung have major contribution to the $E_{b}$. We have shown the characterization of the SN phase with INS experiment and also predicted the $q_{m}-M$ relation. We think that the most of intensities of $S(q,\omega)$ in LiCuVO4 is below $3$ meV and most intense peak is at $q=\pi/2$ and $\omega=0.4$ meV. In this paper we have shown that magnitude of dimerization is vanishingly small, and gs is doubly degenerate in the SN phase.

The model Hamiltonian in Eq. 1 supports many quantum phases in 1D system. There are many open questions to be answered, like how to characterize the SN phase and other multipolar phases, what happens to magnon pairing in large $\alpha$ limit, and how to increase the binding energies of magnon pairing. The RIXS is a good experimental tool which may distinguish the SDW2 and SN phase. In the large $\alpha$ limit the $J_{1}-J_{2}$ chain should behave like two decoupled chains and their behavior looks like anti-ferromagnetic Heisenberg chains. We can ask question like what happens to magnon as low lying excitations should be similar to Heisenberg spin-1/2 one dimensional chain.

Acknowledgement

We thank Prof. Z. G. Soos and Prof. G. Baskaran for their valuable comments. M. K thanks DST for Ramanujan fellowship and computation facility provided under the DST project SNB/MK/14-15/137.