Temperature dependence of the NMR relaxation rate 1/T_1 for quantum spin chains
TTemperature dependence of the NMR relaxation rate 1 / T for quantum spin chains Maxime Dupont, ∗ Sylvain Capponi, † and Nicolas Laflorencie ‡ Laboratoire de Physique Th´eorique, Universit´e de Toulouse and CNRS, UPS (IRSAMC), F-31062, Toulouse, France (Dated: August 31, 2018)We present results of numerical simulations performed on one-dimensional spin chains in orderto extract the so-called relaxation rate 1 /T accessible through NMR experiments. Building onnumerical tensor network methods using the Matrix Product States (MPS) formalism, we can followthe non-trivial crossover occurring in critical chains between the high-temperature diffusive classicalregime and the low-temperature response described by the Tomonaga-Luttinger liquid (TLL) theory,for which analytical expressions are known. In order to compare analytics and numerics, we focus ona generic spin-1 / I. INTRODUCTION
One-dimensional (1d) quantum systems are known tobe very peculiar due to strong quantum fluctuations thatprohibit long-range order and can give rise to unusualphases of matter. In this context, it is remarkable thatquantum spin chains fall generically into two classes re-garding their low-energy properties [1, 2]: (i) critical be-havior where gapless low-energy excitations can be de-scribed in the framework of Tomonaga-Luttinger liquid(TLL) theory; (ii) gapped behavior.Nevertheless, condensed-matter experiments aremostly done on quasi -1d materials, hence the role ofsmall inter-chain couplings (as compared to the dom-inant 1d energy scale J ) may become important atlow-enough temperature (eventually leading to magneticordering). Conversely, at high temperature ( T (cid:29) J ),quantum fluctuations vanish so that a classical pictureemerges. As a consequence, for realistic experimentalsystems the validity of a universal 1d TLL regime is notgranted and should be checked in some unbiased way. Inparticular, understanding the intermediate temperatureregime T ∼ J , highly relevant to understand severalexperimental data, is a great theoretical challengeregarding dynamical observables.In this paper, we focus on nuclear magnetic resonance(NMR) for quantum spin systems [3], and more specif-ically on the 1 /T spin-lattice relaxation rate. Indeed,this quantity contains lots of information on the dynam-ical properties of the system since it is directly relatedto dynamical spin-spin correlations. Moreover, being alocal quantity (a crucial property of NMR technique), wewill argue that reliable data can be obtained even thoughwe will simulate finite spin chains.Being of fundamental interest, the low- T behav-ior of the NMR relaxation rate has been investigated ∗ [email protected] † [email protected] ‡ nicolas.lafl[email protected] for several 1d or quasi-1d quantum magnets. Spin-gapped compounds, such as two-leg ladders SrCu O [4],BiCu PO [5], Sr − x Ca x Cu O [6], weakly coupledHaldane chains Y BaNiO [7], or dimerized spin chainsAgVOAsO [8], exhibit an activated relaxation at low- T . For gapless Heisenberg chain systems, the low-energycritical behavior has been studied [9–12] for Sr CuO which is an almost ideal realization with a large J ∼ T N (cid:39) ) NNi(NO ) or dimerizedspin-1 / /T behavior has been performed in Refs. 16and 17.A useful experimental review on NMR properties ofseveral spin chains can be found in Ref. 18. Note also that1 /T measurements have also been used to characterizeone-dimensional metallic phase in carbon nanotube [19]or quasi-1d superconductor [20].More recently, interesting quasi-1d spin-gapped ma-terials have also been investigated using NMR [21]: ananisotropic spin-1 system NiCl -4SC(NH ) (DTN) anda spin-ladder one (C H N) CuBr (BPCB). In bothcases, 1 /T measurements could be interpreted either ascoming from magnon (respectively spinon) excitations inthe gapped (respectively gapless) 1d phase, and the quan-tum critical regime was also argued to be universal. Mostimportantly, the whole temperature range, including 1das well as 3d regimes, was discussed.Experimentally, when decreasing temperature, theNMR relaxation rate 1 /T has been found to divergein the TLL regime, with power-law governed by a char-acteristic exponent. Such an analysis is used in exper-iments to determine the corresponding TLL exponent K [22, 23]. For example, it was a smoking-gun signatureof attractive TLL in (C H N) CuBr (DIMPY) com- a r X i v : . [ c ond - m a t . s t r- e l ] O c t pound [24, 25]. However, given that we are genericallydealing with quasi-1d materials, critical fluctuations and3d ordering will limit the low-energy 1d regime, and agenuine TLL critical behavior is observable only withinsome finite window in temperature. This remains to beanalyzed more quantitatively, which is the main purposeof this work.The rest of the paper is organized as follows. InSec. II, we present the theoretical models and provideuseful definitions. Section III describes the numericaltechnique based on finite temperature Matrix ProductStates (MPS) approach. Results are then discussed inSec. IV. Finally, we present our conclusion in Sec. V. II. MODELS AND DEFINITIONS
We give in this section the two models that will bestudied in this paper and a small discussion on theirphase diagram. Both models present a TLL gapless phaseand a gapped phase, induced by an external magneticfield. We will also provide definition of the NMR relax-ation rate 1 /T and discuss its expected behavior withtemperature. A. Theoretical models
1. The spin- / XXZ chain
We first consider one of the simplest paradigmatic ex-ample of TLL liquid, namely the spin-1 / H XXZ = J L − (cid:88) j =1 (cid:0) S xj S xj +1 + S yj S yj +1 + ∆ S zj S zj +1 (cid:1) − gµ B h L (cid:88) j =1 S zj (2.1)where ∆ ∈ ( − ,
1] denotes the Ising anisotropy, J thecoupling strength and h is an applied magnetic field inthe z direction with g the gyromagnetic factor and µ B the Bohr magneton constant. The Hamiltonian is definedwith open boundary conditions (OBC), as will be usedin our numerical simulations.In the range ∆ ∈ ( − ,
1] the XXZ model can be de-scribed by a TLL as long as its spectrum remains gap-less [2]. As a function of magnetic field, the gaplessregimes extends up to a critical field gµ B h c = J (∆ + 1),and the system becomes gapped for h > h c . In the lat-ter regime, the gap increases linearly with the appliedmagnetic field, ∆ g = gµ B ( h − h c ), see Fig. 1. h c . . . m z TLL phase Polarized phase h c h c . . . m z TLL phase P o l a r i z e dp h a s e L a r g e - D p h a s e FIG. 1. (color online) (i) Upper panel : magnetization curveof the XXZ Hamiltonian (2.1) as a function of the magneticfield for ∆ ∈ ( − , h c = J (∆ + 1) /gµ B and a gapped one above when the system isfully polarized. (ii) Lower panel : magnetization curve of the1d DTN Hamiltonian (2.2) as a function of the magnetic field.Two gapped phases (large- D and polarized) are respectivelylocated below h c and above h c . The intermediate gaplessphase can be described by a TLL theory.
2. The quasi- d spin- compound “DTN” We also discuss a quasi-1d magnetic insulator com-pound NiCl -4SC(NH ) , also called DTN, whose rele-vant 3d structure consists of weakly coupled S = 1 chainsin the two other transverse (with respect to the chainaxis) directions. Its experimental interest comes fromthe appearance of a Bose-Einstein condensation (BEC)phase when applying a magnetic field at low tempera-ture [26, 27]. More recently, Br-doped (disordered) DTNwas suggested to be a good experimental candidate forobserving a Bose glass phase [28–30].Although there is 3d magnetic order observed below T N ∼ J /J (cid:39) .
08, oneexpects 1d physics and a TLL regime at higher T . Theeffective Hamiltonian to describe this situation reads H DTN − = J L − (cid:88) j =1 S j · S j +1 + L (cid:88) j =1 (cid:104) D (cid:0) S zj (cid:1) − gµ B hS zj (cid:105) , (2.2)where S j = ( S xj , S yj , S zj ) are spin-1 operators. In the cur-rent literature, [31] J = 2 . D = 8 . h is given in Tesla with g = 2 . D , the system is in theso-called large- D phase [33]. This is a trivial phase, adi-abatically connected to the product state | → . . . →(cid:105) where each state is in a non-magnetic S z = 0 eigenstate.Clearly, this phase has a finite spin-gap, which corre-sponds to the first critical field h c needed to magnetizethe system. Its value is known to be, at first order in J/D (cid:28) h c /gµ B = D − J + O ( J /D ) (cid:39) h ∈ [ h c , h c ], with h c /gµ B = D + 4 J = 11 .
40 T. Abovethis critical saturation field, the system becomes gappedagain, entering a fully polarized phase. As a side remark,we recall that in the true 3d material DTN, both criti-cal fields are shifted due to interchain couplings, so that h c = 2 . h c = 12 .
32 T [36].In the TLL phase and close to the upper critical field h c , the DTN Hamiltonian (2.2) can be mapped towardan effective XXZ model of spins S = 1 / J = 2 J and ˜∆ = 0 .
5. This result can be refined us-ing contractor renormalization (CORE) method [38, 39]leading to the same value of ˜ J but a slightly reduced˜∆ = 0 .
36. Both mappings lead to a value of the effectivemagnetic field ˜ h = h − J − D . B. Relaxation rate 1 / T The nuclear spin-lattice relaxation rate T − measuredby NMR [40] is basically testing the local and dynamical spin correlation function S aann ( ω ) with a = x, y, z , and ω being the NMR frequency at a given site n ,1 T = γ (cid:110) A ⊥ [ S xxnn ( ω ) + S yynn ( ω )] + A (cid:107) S zznn ( ω ) (cid:111) = 1 T ⊥ + 1 T (cid:107) (2.3)Here, A ⊥ is the transverse hyperfine coupling constant, A (cid:107) the longitudinal one, γ the gyromagnetic ratio and S abnn ( ω ) = Re (cid:26) (cid:90) ∞ d t e iωt (cid:2) (cid:104) S an ( t ) S bn (0) (cid:105)− (cid:104) S an ( t ) (cid:105)(cid:104) S bn (0) (cid:105) (cid:3) (cid:27) , (2.4)with S a = S b † , S a ( t ) = e i H t S a e − i H t , and (cid:104)(cid:105) is the ther-mal average defined later in (3.4). For convenience, the x and y spin components can be expressed using the raisingand lowering operators, S xxnn ( ω ) + S yynn ( ω ) = 12 (cid:2) S + − nn ( ω ) + S − + nn ( ω ) (cid:3) . (2.5) It is theoretically justified to take the limit ω → S + − nn ( ω ) and S − + nn ( ω ) become equiva-lent in this limit ω → /T is experimentally gov-erned by the hyperfine coupling tensors A ⊥ and A (cid:107) . Tofavor one over the other, a specific nucleus can be tar-geted for the NMR experiment. In the case of DTN,proton H (nuclear spin I = 1 /
2) probes both compo-nents while nitrogen N ( I = 1) probes dominantly thetransverse one [36].Hyperfine coupling tensors put aside or set equal toone, the low temperature behavior of the transverse andlongitudinal components of 1 /T depends on microscopicparameters of the model. At low temperature, the trans-verse component is larger than the longitudinal one and itis thus justified to consider 1 /T (cid:39) /T ⊥ . In the follow-ing we will fix γ A ⊥ , (cid:107) = 1, and compute the relaxationrates 1 T ⊥ = (cid:90) ∞ d t Re (cid:2) (cid:104) S ± n ( t ) S ∓ n (0) (cid:105) (cid:3) , (2.6)and1 T (cid:107) = (cid:90) ∞ d t Re [ (cid:104) S zn ( t ) S zn (0) (cid:105) − (cid:104) S zn ( t ) (cid:105)(cid:104) S zn (0) (cid:105) ] . (2.7)Note that the single operator averages are time inde-pendent, (cid:104) S zn ( t ) (cid:105) = (cid:104) S zn (0) (cid:105) . C. Tomonaga-Luttinger liquid description in thegapless phase
The 1d gapless phase, experimentally accessible bytuning a control parameter such as the external mag-netic field, can be effectively described by the TLLHamiltonian[2], H TLL = 12 π (cid:90) d r (cid:110) uK [ ∂ r θ ( r )] + uK [ ∂ r φ ( r )] (cid:111) (2.8)where u is the velocity of excitations and K is the dimen-sionless TLL parameter. They both fully characterize thelow-energy properties of the system and are thus model-dependent. θ ( r ) and φ ( r ) are bosonic fields obeying thecommutation relation [ φ ( x ) , θ ( y )] = iπδ ( x − y ).For the XXZ model (2.1), the TLL parameters K and u are known from Bethe ansatz equations [41]. At zeromagnetic field, analytical expressions are known as afunction of the Ising anisotropy ∆: K = π − ∆) and uJ = π √ − ∆ . (2.9)For generic non-integrable models, the TLL parameterscan be obtained numerically using DMRG by fittingstatic correlation functions [42, 43] which has been suc-cessfully done in the past years for various quasi-1d com-pounds [17, 44].In the TLL framework, the dynamical correlation S abnn ( ω ) defined in (2.4) can be computed analytically asa function of the temperature in the “ low energy limit ”which we will try to define more precisely in this paper.Let us recall that the dynamical spin susceptibility is de-fined as, χ abij ( t ) = − i Θ ( t ) (cid:104) (cid:2) S ai ( t ) , S bj (0) (cid:3) (cid:105) , (2.10)where Θ( t ) is the Heaviside function. Since a typicalrelevant case in experiments is the local quantity i = j = n , we only consider this case in the following. Infrequency space, the susceptibility can be related to thedynamical spin correlation function by [45] S abnn ( ω ) = 2 e − βω − (cid:2) χ abnn ( ω ) (cid:3) . (2.11)In the low energy limit βω (cid:28) T ⊥ =2 A x cos (cid:0) π K (cid:1) u (cid:18) πTu (cid:19) K − B (cid:18) K , − K (cid:19) , (2.12)and1 T (cid:107) = A z cos ( πK )2 u (cid:18) πTu (cid:19) K − B ( K, − K )+ KT πu , (2.13)with B ( x, y ) the Euler beta function and A x,z prefac-tors of the static correlation functions. Thus, generically1 /T ⊥ ( T ) diverges at zero temperature as a K -dependentpower-law, and dominates over 1 /T (cid:107) . Note that for fi-nite magnetic field (equivalent to non half-filled case),additional subleading corrections are expected [48]. D. Gapped regime behavior
In contrast to a TLL gapless phase where 1 /T haspower-law behavior at low temperature, we also considerthe gapped regime where fluctuations are exponentiallysuppressed [16, 49–53], such that for T < ∆ g T ⊥ , (cid:107) ∝ exp (cid:0) − α ⊥ , (cid:107) ∆ g /T (cid:1) , (2.14) where ∆ g is the energy gap of the system, and α ⊥ , (cid:107) an O (1) prefactor which depends on the relaxation pro-cesses, and also on the temperature range [52]. Below,in section IV B 2, we show numerical results for the high-field gapped regime of the XXZ chain where our dataare perfectly described by α ⊥ , (cid:107) = 1. At higher temper-ature above the gap, one may also expect a non-trivialcrossover to TLL regime [54]. III. NUMERICAL METHODS
To get the relaxation rate one needs to obtain in thefirst place the dynamical correlation (cid:10) S an ( t ) S bn (0) (cid:11) . Weused the TEBD (Time Evolution Block Decimation) al-gorithm [55] with both real/imaginary time through theMPS formalism adapted for 1d systems [56]. A generalone-dimensional system containing L sites with OBC canbe represented by the following MPS, | Ψ (cid:105) = (cid:88) { s i } A s a A s a a · · · A s L a L − | s (cid:105)| s (cid:105) · · · | s L (cid:105) (3.1)where the local index s i is the physical index represent-ing an element of the local Hilbert space at site i . Itsdimension is d and is equal to 2 ( ↑ and ↓ ) for spin-1 / ↑ , ↓ and → ) for spin-1. We note a i the bond in-dex whose dimension is directly related to the “numberof states” m to describe the system, meaning that m isa control parameter in the numerical simulations.The first step is to perform an imaginary time evo-lution on the system to reach the desired temperature.Once the state is at hand, the second step consists ofevolving it through a real time evolution. At each timestep, the correlation is measured. When all the datain time-space have been obtained, a numerical FourierTransform can be performed to get the data in frequency-space. A. Time evolution with MPS
We will be general and consider the case with a Hamil-tonian H consisting of nearest-neighbor interactions only– it is the case for the XXZ (2.1) or DTN Hamiltonian(2.2) introduced before. Now, we need to evolve our MPSup to a time t . The operation can be discretized us-ing smaller time steps τ such that t = N τ , leading to e − it H = (cid:81) N e − iτ H . If the time step τ is small enough,a first (or higher) order Trotter decomposition can beperformed, e − it H (cid:39) N (cid:89) e − iτ H even N (cid:89) e − iτ H odd + O ( τ ) (3.2)where H even and H odd respectively correspond to theeven and odd bond Hamiltonians only acting on twonearest-neighbor spins. The decomposition is possiblebecause even – odd – bond Hamiltonians commute witheach others. But it is not exact and leads to an error in τ due to the fact that [ H even , H odd ] (cid:54) = 0. The advantageis that the bond Hamiltonians can easily be diagonalizedand exponentiated since they are only d × d matrices.And so, applying successive evolution gates on theMPS as well as singular value decompositions to restorethe MPS original tensor-site dependent form (3.1) willeventually lead to a time- t evolved state. B. Finite temperature with MPS
It was useful to introduce time-evolution concepts alsoto discuss finite-temperature with MPS [57]. The mainidea is to represent the density matrix ρ β of the physical(mixed) state in an artificially enlarged Hilbert space asa pure state | Ψ β (cid:105) – which is what we can deal with inthe MPS formalism. The auxiliary space can simply beconstructed as a copy of the original one.Assuming that we know the purification of the densitymatrix ρ β =0 as a wave function | Ψ β =0 (cid:105) it can be shownthat an imaginary time evolution has to be performedover the infinite temperature state in order to get thefinite temperature state, | Ψ β (cid:105) = e − β H / | Ψ β =0 (cid:105) (3.3)with the Hamiltonian only acting on the physical sites.This imaginary time evolution can be performed usingthe TEBD algorithm described before in III A. Expecta-tion values can then be measured at inverse temperature, (cid:104)O(cid:105) β = Tr (cid:2) O e − β H (cid:3) Tr [ e − β H ] = (cid:104) Ψ β |O| Ψ β (cid:105)(cid:104) Ψ β | Ψ β (cid:105) (3.4)For this procedure to work, the initial state | Ψ β =0 (cid:105) hasto be a product state of Bell states between each physicalsite and its associated auxiliary site, | Ψ β =0 (cid:105) = 1 √N L (cid:89) n =1 (cid:88) { s } | p ns a ns (cid:105) (3.5)with | p (cid:105) corresponding to physical sites, | a (cid:105) to auxiliaryones and N a normalization constant. The summation isover the d possible local states s . Such a state is simpleenough to be built exactly in the MPS formalism. C. Numerical limitations
The main limitations are about the temperature andthe final time one can reach using the methods describedabove. The reason in both cases is directly related to arapid growth in the entanglement entropy while evolvingthe state. This implies to keep larger and larger num-ber of states m in the MPS if one wants to be accurate,strongly limiting numerical simulations in practice. On the one hand, it becomes increasingly difficult toreach low temperatures. Indeed, one expects a volume-law entanglement entropy ( i.e. linear with the systemsize L ) due to the auxiliary sites which are used to purifythe thermal state. As a consequence, the number of keptstates m needed to describe accurately the system willgrow exponentially as the temperature T decreases.On the other hand, the maximal (real) time that canbe reached is of the order of few tens of J − typically,for similar reasons as discussed above, namely the lineargrowth of entanglement entropy with time [58]. Thus,fixing a maximum number of kept states m limits sim-ulations to a finite time t m ax . Note that this limitationapplies at all temperatures, even T = 0.Despite these severe limitations, recent progress in thefield has allowed some improvements. For instance, wewill make use of the auxiliary degrees of freedom whichare used to purify the thermal state by time-evolvingthem with −H , which is mathematically exactly the samebut has been shown to improve substantially the timerange [59]. By construction, this trick only applies tofinite-temperature simulations though. Last, althoughits use did not prove to be systematically reliable in ourcase, we would like to mention the possibility to useso-called linear prediction technique, coming from dataanalysis [60] which aims at predicting “longer time” be-havior from the knowledge of dynamical correlations at“intermediate time”. IV. RESULTS
We provide in this section our numerical results [61]about the relaxation rate 1 /T using models and tech-niques presented in the previous sections II and III. Firstof all we will focus on the XX model (equivalent to freefermions) for which we can compute exactly the dynam-ical correlations for all temperatures and that will serveas a benchmark for our simulations. Next, we will turnto the interacting XXZ case for S = 1 /
2, and then to a S = 1 chain model relevant to the DTN material. A. Case study : XX point ( ∆ = 0 ) The XXZ Hamiltonian (2.1) at ∆ = 0, known as XXmodel, can be mapped onto a model of free spinlessfermions using a Jordan-Wigner transformation. We re-strict ourselves to h = 0. It can be diagonalized in Fourierspace with ε k = J cos k and k = nπL +1 with n = 1 , , . . . , L considering open boundary conditions, H XX = J L − (cid:88) j =1 (cid:16) c † j c j +1 + h . c . (cid:17) = (cid:88) k ε k c † k c k . (4.1)Unlike the bosonization expressions (2.12) and (2.13)which are only valid in the low-energy limit, the resultspresented in this section will be valid for all regimes. We tJ − . . . . . R e (cid:20) (cid:28) S ± , z L ( t ) S ∓ , z L ( ) (cid:29) (cid:21) k correlations
10 20 30 40 tJ − − − − − − − ⊥ correlations − − − ω/J . . . . . . R e (cid:20) F T t → ω (cid:28) S ± , z L ( t ) S ∓ , z L ( ) (cid:29) (cid:21) k correlations − − ω/J . . . . . ⊥ correlations β = 0 . β = 1 . β = 2 . β = 3 . β = 4 . β = 5 . β = 6 . β = 8 . β = 10 . ExactNumeric
FIG. 2. (color online) We compare numerical results (circles) used to determine the 1 /T with analytical results (straight lines)for the XX model. For the transverse ( ⊥ ) case, exact results are computed on a chain of size L = 64 (OBC). As for thelongitudinal ( (cid:107) ) case, the chain size is L = 1000 (OBC). Numerics on their side are performed on a chain of L = 64 (OBC)sites. The left panel shows the real value of the dynamical correlations. For readability, we only display numerical results forthe lowest temperature for the longitudinal correlations. Indeed, this is a priori the hardest to compute and thus the mostsubject to errors. The right panel shows the real part of the Fourier transform of the real time data. Although we only showdata up to t = 40 J − the Fourier Transform of the exact zz correlations was performed using data up to t = 1000 J − . present the details of the calculations in appendix A forthe analytical exact expressions of the dynamical corre-lations for the XX model using precisely the same condi-tions as in our numerical simulations ( i.e. a finite chainlength with open boundary conditions).We show the ‘bare’ results in Fig. 2 that will be usedto obtain the relaxation rate thereafter. Ideally, one isinterested in the thermodynamic limit ( i.e. L → ∞ ) butwe see that, at finite temperature (hence finite correla-tion length), working on finite length chains with only amoderate number of sites L allows to get reliable data.Indeed, the MPS estimates agree perfectly with the exactexpressions (see appendix A).First of all, we consider the local dynamical correlationof the site in the middle of the chain reducing de factoboundary effects. Then, as finite size effects are known tobe caused by the reflection of the propagating excitationsat TLL velocity u on the boundary of the system, one canestimate a time below which the dynamical correlationscan be considered as free of finite size effects (basically, ut ∼ L ).We first discuss the transverse correlations, see Fig. 2.For all temperatures, they decay rather quickly to zero,so that we can safely truncate data to a maximum time t max (which is anyway a natural cutoff provided by theinverse of the NMR frequency ω ) and get reliable valuesof 1 /T ⊥ by integrating over time. Moreover, we have alsochecked that finite size effects are extremely small sincewe are computing a local correlation.The same cannot be said for the longitudinal correla-tions. They continue to oscillate even for high temper-atures and long times, and their amplitude gets (very)slowly smaller with time. This implies severe limitationsto get data in the thermodynamic limit. For instance,exact computations using (A2) were done on L = 1000 and still displayed oscillations of amplitude around 10 − at t = 1000 J − . This makes the value of 1 /T (cid:107) verydifficult to estimate. This well-known behavior is relatedto spin diffusion-like behavior [62, 63] which cause a log-arithmic divergence at small frequency ω . However, wehave to remember that the NMR frequency ω eventuallyprovides a natural cutoff.For completeness we display the real part of the Fouriertransform on the right panels of Fig. 2 for which the 1 /T value as defined in section II B corresponds to the ω = 0value. B. Spin-1/2 XXZ chain at ∆ (cid:54) = 0
1. Gapless regime
Building on the perfect agreement observed previouslybetween MPS estimates and the exact analytical solu-tion of the XX model, we are now confident to extendour study of the more generic XXZ case − < ∆ ≤ L = 64 with a cutoff of ε = 10 − in the singularvalues. We kept a maximal number of D = 500 states.A fourth order Trotter decomposition was used with aTrotter step of τ = 0 . ∼ T K − occurs only below T /J ∼ . − . there are nofree parameters in the analytic expressions . Indeed, the .
04 0 . T /J . / T ⊥ ∆ = − . ∆ = − . ∆ = − . ∆ = − . ∆ = − . ∆ = . ∆ = . ∆ = . ∆ = . ∆ = . ∆ = . ∆ = . TLL TheoryNumericTLL TheoryNumeric ∞ FIG. 3. (color online) Transverse relaxation rate 1 /T ⊥ vs. reduced temperature T /J for the spin-1 / h = 0 obtained numerically us-ing MPS techniques (circles, from top to bottom: ∆ = − . , − . , − . , − . , − . , . , . , . , . , . , . , | ∆ | <
1, and with Eq. (4.2) forthe SU(2) Heisenberg point ∆ = 1. The thin lines betweenthe circles are guides to the eyes.
TLL parameters are computed using the exact expres-sions Eq. (2.9) for u and K , and A x is obtained followingRefs. 64 and 65. The isotropic limit ∆ = 1 is a spe-cial point where logarithmic corrections appear in sev-eral quantities [66–68], leading to a very slow divergenceof the (isotropic) NMR relaxation rate [12, 69]1 T (cid:39) √ π (cid:115) ln Λ T + 12 ln (cid:18) ln Λ T (cid:19) , (4.2)where Λ (cid:39) . J . MPS estimates compare well withthis parameter-free expression, as visible in Fig. 3.Interestingly, we notice the non-monotonic behavior of1 /T ⊥ with temperature only when ∆ (cid:38) β = 0), the value of 1 /T ⊥ does not dependon the sign of ∆, which is expected since the many-bodyspectrum of H ∆ is an odd function of ∆. Its value is min-imum for ∆ = 0 with 1 /T ⊥ = √ π/ (2 J ) [70] and increaseswith | ∆ | . At the isotropic point | ∆ | = 1 we expect therelaxation rate to diverge due to the diffusion-like behav-ior [62, 63] of the dynamical correlation function. Our re-sults at infinite- T go beyond Baker-Campbell-Hausdorffexpansion developed up to O ( t ) in Ref. 71 to com-pute (cid:104) S ± j ( t ) S ∓ j (0) (cid:105) at short times, which would suggest1 /T ⊥ ∼ J − (1 + ∆ ) − . This prediction is in contrast towhat we found, namely the transverse relaxation rate in-creasing with | ∆ | . Indeed, while such an expansion finds βJ − − − − − − − − − / T ⊥ , k ∆ g = . ∆ g = . ∆ g = . h/gµB = 2 . Jh/gµB = 2 . Jh/gµB = 3 . J /T k /T ⊥ /T k /T ⊥ FIG. 4. (color online) Transverse and longitudinal relax-ation rates 1 /T ⊥ , (cid:107) plotted against reduced inverse temper-ature βJ for the spin-1 / .
5. The critical magnetic field is h c = 3 J/ gµ B and the value of the gap ∆ g = gµ B ( h − h c ).Numerical results are obtained using MPS techniques (cir-cles and diamonds) and the exponentially decaying behavioris verified with the straight lines set with the expected gapvalue 1 /T ⊥ , (cid:107) = c ⊥ , (cid:107) · e − β ∆ g and c ⊥ , (cid:107) a non-universal freeparameter. the correct gaussian behavior for ∆ = 0 (free-fermions),higher-order terms have to be taken into account for | ∆ | > | ∆ | .
2. Gapped XXZ chain
We then set the anisotropy value to ∆ = 0 . /T ⊥ , (cid:107) are plot-ted in Fig. 4 where we observe an excellent agreementwith an exponentially activated behavior ∼ exp( − β ∆ g ),where ∆ g is the spin gap. We notice that as the gap getssmaller, the lower the temperature has to be to observethe exponential law. C. DTN
We now move to the DTN compound in its 1d limitdescribed by Eq. (2.2). We compute the relaxation ratesfor various values of the magnetic field h , mainly closeto h c which is relevant for NMR experiments [21]. Itis a more challenging system to simulate than the XXZmodel as it is made of spins S = 1 (enlarged local Hilbertspace). The simulations were performed on open chainsof size L = 64 with a cutoff of ε = 10 − in the singularvalues. We kept a maximal number of D = 150 states. . . . / T ? h = 9 . h = 9 . h = 10 . h = 10 . h = 11 . TLL TheoryNumericTLL TheoryNumeric . . T [ K ]0 . . / T ? Zoom at low T . . . . . . m z . . . . / T ? T = 0 . . . . . . h [T] FIG. 5. (color online) Transverse relaxation rate 1 /T ⊥ plot-ted vs. temperature T for the spin-1 DTN chain obtainednumerically using MPS techniques (circles). The low temper-ature behavior is compared to TLL prediction (straight lines).The magnetic field h is given in Tesla. The inset comparesTLL prediction and numerical results for T = 0 . h c to h c . The lower panel isa zoom on the low temperature asymptotic power-law regime. A fourth order Trotter decomposition was used with aTrotter step of τ = 0 . T (cid:39) . T /J ∼ . h c . We point out that there are again no adjustablecoefficients, the TLL parameters being computed inde-pendently using standard DMRG [72]. The tiny differ-ence that appears at low temperature between numericaldata and TLL is due to the limited number of states m kept when performing calculations. Though this doesnot dispute the TLL prediction, it reveals the challengein such time-dependent simulations. The inset in Fig. 5shows the transverse relaxation rate at T = 0 . h c to h c . Once more, there is a very good agree-ment between numerics and TLL theory except when onegets close to the critical fields. Indeed, as we clearly seein the lower panel of Fig. 5 for h = 11 . T = 0 . /T ⊥ observed in theXXZ model is absent for the DTN and may seem oddat first place since it can be mapped effectively onto a S = 1 / . .
36 and could thusbe compared with Fig. 3. However this non-monotonicvariation is observed at high temperature while this map-ping is only justified in the low-energy limit as discussedin II A 2.One can also try to compare the relaxation rates ofFig. 5 with the NMR data for the DTN compoundgiven in Ref 21. What draws our attention is the non-monotonic regime of 1 /T ⊥ observed at high temperatureexperimentally, which, as we have just discussed, is nottheoretically predicted for a single DTN chain. Yet itcannot be attributed to 3d effects as J = 0 .
18 K isvery small compared to the temperature T . We then ob-served that experiments are performed by proton ( H )NMR which probes both /T ⊥ and 1 /T (cid:107) .We therefore interpret this effect as due to the parallelcontribution of the relaxation rate. We show in Fig. 6both the transverse and longitudinal 1 /T ⊥ , (cid:107) as a functionof temperature. We cannot precisely estimate the valueof 1 /T (cid:107) due to its dependence on ω (and therefore onour maximum time in numerical simulations) so that wegive a lower bound. Its high temperature contributionto the total relaxation rate clearly dominates over thetransverse part and explains well the experimental non-monotonic regime at high T .Perhaps more importantly, as displayed in Fig. 6, the3d BEC ordering observed in DTN [26, 36] for m z (cid:39) . T N (cid:39) .
59 K occurs above the asymptotic regimewhere the genuine TLL power-law behavior is expected.It is therefore impossible to directly extract TLL expo-nents in DTN, because of interchain effects that even-tually lead to an ordering of the coupled TLLs. Ideallywe would expect for quasi-1d systems the TLL descrip-tion of the NMR relaxation to be valid in the followingtemperature regime: J (cid:29) T (cid:29) J .Concerning the difficulty to obtain reliable data at hightemperature for the longitudinal 1 /T (cid:107) , it is well knownthat this is due to spin diffusion-like behavior [62, 63].Therefore, measurements should in principle depend ex-plicitly on the NMR frequency ω . V. CONCLUSION
Performing time-dependent numerical simulations atfinite temperature on 1d systems to compute the NMRrelaxation rate 1 /T , we have discussed the temper-ature range validity of analytical predictions for twomodels (i) the paradigmatic example for Tomonaga-Lutinger liquids: the spin-1 / S = 1 Hamil-tonian, relevant for experiments on the DTN compound . T [ K ]0 . . / T ⊥ , k /T ⊥ /T k mz ’ . h
1d = 11 . T h
3d = 11 . TTN ’ . K exp. ∞ FIG. 6. (color online) Longitudinal 1 /T (cid:107) and transverse 1 /T ⊥ relaxation rates for the DTN spin-1 chain at h = 11 . m z (cid:39) .
85. As 1 /T (cid:107) cannot be estimated forsure, we only provide a lower bound. The non-monotonic be-havior observed experimentally at high T in Ref 21 apparentlycomes from the large contribution of 1 /T (cid:107) at high tempera-ture. Experimental data for DTN [21] at the same magne-tization are shown for comparison, after a proper rescalingin order to match the low- T regime. The 3d BEC transitiontemperature T N ( m z (cid:39) . (cid:39) .
59 K [36] is also shown. as a function of an external magnetic field.Both models present in some regime a gapless phasethat can be described by TLL “ low-energy ” theory, witha relaxation rate dominated by its transverse component1 /T ⊥ ∼ T K − algebraically diverging at low tempera-ture, where K is the dimensionless TLL exponent. Weobserved that the expected power-law behavior occursonly below T /J ∼ . − .
2, thus defining the low-energy limit of validity of TLL theory an order of magni-tude below the energy scale J of the system. It is im-portant to be able to define this limit as TLL predictionsare often used experimentally on quasi-1d compounds toextract the value of K . As a consequence, we believethat it remains experimentally challenging [15], and of-ten impossible, to explore a genuine critical 1d regimein quasi-1d compounds when J is small and 3d order-ing prevents a wide TLL regime. For instance, we haveshown that for DTN, the BEC ordering temperature islarger than the crossover temperature towards TLL be-havior.We have also studied the transverse relaxation ratesof these two models in other regimes than TLL theory.First, we considered high temperatures, with a peculiarnon-monotonic behavior in the S = 1 / T , which does not exist forthe 1d S = 1 model of DTN. However, such a non-monotonic dependence with temperature at high T isexperimentally observed in DTN. We showed that thiseffect comes from the parallel contribution of the relax- ation rate 1 /T (cid:107) dominating at high temperature overthe transverse part. Finally, we verified that in gappedphases the relaxation rates are exponentially suppressed ∼ exp ( − ∆ g /T ) and can indeed lead to accurate deter-minations of the spin gap, at least in a regime of tem-perature T (cid:28) ∆ g since other relaxation mechanisms canchange the activated behavior at higher temperature.We want to emphasize again the role of 3d orderingat finite temperature, preventing the observation of a 1dTLL regime. As discussed for the particular case of DTN,one needs a hierarchy of energy scale J (cid:29) T (cid:29) J tobe able to directly extract the TLL exponent K from thedivergence of T − with T .Concerning future advances for quasi-1d systems, wecan envision trying to simulate imaginary-time dynamicsusing quantum Monte-Carlo techniques, provided thatthe model has no minus-sign problem. While it will benecessary to perform a numerical analytic continuation(using for instance Maximum Entropy techniques), wehave some hope that this could lead to reliable results toquantitatively capture the influence of interchain effects.For 1d chains, it has been rather successful [10]. As amatter of fact, our extensive 1d results could serve asuseful benchmarks for that too.While we have considered various chains, we are farfrom being exhaustive. Indeed, there are some TLLmodels for which elementary excitations may not besimple spin flips, for instance multipolar nematic phasefor which 1 /T behavior will be different [73, 74]. Ina similar line of thought, we could imagine simulatingmore complicated models including charge and spindegrees of freedom to describe NMR relaxation inmetallic or superconducting wires. Note added :
While completing this work, a relatednumerical study by Coira et al. has appeared [75]. Ourresults are perfectly compatible with each other whencomparison can be made, such as the transverse 1 /T data for a single spin-1 / ≥ ACKNOWLEDGMENT
The authors are grateful to M. Horvati´c and M.Klanjˇsek for their careful reading of this manuscript andthoughtful comments and to R. Blinder, T. Giamarchiand E. Orignac for valuable discussions. This work wasperformed using HPC resources from GENCI (Grant No.x2015050225 and No. x2016050225), and is supported bythe French ANR program BOLODISS and R´egion Midi-Pyr´en´ees.0
Appendix A: Dynamical correlations for the XXmodel
For completeness, we remind the reader the exact ex-pressions for time-displaced spin correlations in the ex-actly solvable XX model.
1. Longitudinal correlations
The longitudinal correlations (cid:104) S zi ( t ) S zj (0) (cid:105) for the XXmodel are basically density correlations when performingthe Jordan-Wigner transformation. Since we are inter-ested in the temperature dependence of the correlations,what we need to compute is actually, (cid:104) S zi ( t ) S zj (0) (cid:105) = 1 Z Tr (cid:2) e i H t S zi e − i H t S zj e − β H (cid:3) (A1)with Z = Tr e − β H the partition function. Calculationslead to[76], (cid:104) S zi ( t ) S zj (0) (cid:105) = 14 I I (A2)with I , I = 2 L + 1 (cid:88) k sin ki sin kj (cid:20) βε k (cid:21) e − iε k t (A3)and I , I = 2 L + 1 (cid:88) k sin ki sin kj (cid:20) − tanh βε k (cid:21) e + iε k t . (A4)Also, the single operator averages are (cid:104) S zi ( t ) (cid:105) = (cid:104) S zi (0) (cid:105) = 1 − L + 1 (cid:88) k sin kie − βε k + 1 . (A5)
2. Transverse correlations
The transverse correlations (cid:104) S ± i ( t ) S ∓ j (0) (cid:105) have a morecomplicated structure in the fermion representation due to the string of operators in the exponential[77]. We in-troduce the following identity e iπc † l c l = ( − c † l c l = A l B l with A l = c † l + c l and B l = c † l − c l , leading to2 (cid:104) S ± i ( t ) S ∓ j (0) (cid:105) = (cid:42)(cid:34) i − (cid:89) l =1 A l ( t ) B l ( t ) (cid:35) A i ( t ) (cid:34) j − (cid:89) l =1 A l (0) B l (0) (cid:35) A j (0) (cid:43) . (A6)Now thanks to Wick’s theorem this product of manyfermion operators can be rewritten as elementary expec-tation values of two operators through the Pfaffian ofsome skew-symmetric matrix. Its elements above the di-agonal (which is purely made of zeros) are (cid:104) A ( t ) B ( t ) (cid:105) (cid:104) A ( t ) A ( t ) (cid:105) . . . (cid:104) A ( t ) A j (0) (cid:105)(cid:104) B ( t ) A ( t ) (cid:105) . . . (cid:104) B ( t ) A j (0) (cid:105) . . . . . . (cid:104) B j − (0) A j (0) (cid:105) (A7)At finite temperature, (cid:104) A i ( t ) B j (0) (cid:105) = 1 Z Tr (cid:2) e i H t A i e − i H t B j e − β H (cid:3) . (A8)The two-body expectation values can be computed goingin Fourier space, (cid:104) A i ( t ) A j (0) (cid:105) = 2 L + 1 (cid:88) k sin kj sin ki (cid:20) cos ε k t − i sin ε k t tanh βε k (cid:21) (A9)and, (cid:104) A i ( t ) B j (0) (cid:105) = 2 L + 1 (cid:88) k sin kj sin ki (cid:20) − sin ε k t − i cos ε k t tanh βε k (cid:21) . (A10)And equivalently (cid:104) B i ( t ) B j (0) (cid:105) = −(cid:104) A i ( t ) A j (0) (cid:105) as wellas (cid:104) B i ( t ) A j (0) (cid:105) = −(cid:104) A i ( t ) B j (0) (cid:105) . [1] F. D. M. Haldane, Phys. Rev. Lett. , 1153 (1983).[2] T. Giamarchi, Quantum Physics in One Dimension (Clarendon Press, Oxford, UK, 2004).[3] M. Horvati´c and C. Berthier, in
High Magnetic Fields ,Lecture Notes in Physics No. 595, edited by C. Berthier,L. P. L´evy, and G. Martinez (Springer Berlin Heidelberg,2002) pp. 191–210.[4] M. Azuma, Z. Hiroi, M. Takano, K. Ishida, and Y. Ki-taoka, Phys. Rev. Lett. , 3463 (1994). [5] L. K. Alexander, J. Bobroff, A. V. Mahajan,B. Koteswararao, N. Laflorencie, and F. Alet, Phys. Rev.B , 054438 (2010).[6] K. Magishi, S. Matsumoto, Y. Kitaoka, K. Ishida,K. Asayama, M. Uehara, T. Nagata, and J. Akimitsu,Phys. Rev. B , 11533 (1998).[7] T. Shimizu, D. E. MacLaughlin, P. C. Hammel, J. D.Thompson, and S.-W. Cheong, Phys. Rev. B , R9835(1995). [8] N. Ahmed, P. Khuntia, H. Rosner, M. Baenitz, A. A.Tsirlin, and R. Nath, arXiv:1512.01650 (2015).[9] S. Sachdev, Phys. Rev. B , 13006 (1994).[10] A. W. Sandvik, Phys. Rev. B , R9831 (1995).[11] O. A. Starykh, A. W. Sandvik, and R. R. P. Singh, Phys.Rev. B , 14953 (1997).[12] V. Barzykin, Phys. Rev. B , 140412 (2001).[13] M. Takigawa, O. A. Starykh, A. W. Sandvik, andR. R. P. Singh, Phys. Rev. B , 13681 (1997).[14] T. Goto, T. Ishikawa, Y. Shimaoka, and Y. Fujii, Phys.Rev. B , 214406 (2006).[15] K. Izumi, T. Goto, Y. Hosokoshi, and J.-P. Boucher,Physica B , 1191 (2003).[16] E. Orignac, R. Citro, and T. Giamarchi, Phys. Rev. B , 140403 (2007).[17] M. Klanjˇsek, H. Mayaffre, C. Berthier, M. Horvati´c,B. Chiari, O. Piovesana, P. Bouillot, C. Kollath,E. Orignac, R. Citro, and T. Giamarchi, Phys. Rev.Lett. , 137207 (2008).[18] H.-A. Krug von Nidda, N. B¨uttgen, and A. Loidl,The European Physical Journal Special Topics , 161(2009).[19] Y. Ihara, P. Wzietek, H. Alloul, M. H. R¨ummeli, T. Pich-ler, and F. Simon, EPL (Europhysics Letters) , 17004(2010).[20] H. Z. Zhi, T. Imai, F. L. Ning, J.-K. Bao, and G.-H.Cao, Phys. Rev. Lett. , 147004 (2015).[21] S. Mukhopadhyay, M. Klanjˇsek, M. S. Grbi´c, R. Blin-der, H. Mayaffre, C. Berthier, M. Horvati´c, M. A. Conti-nentino, A. Paduan-Filho, B. Chiari, and O. Piovesana,Phys. Rev. Lett. , 177206 (2012).[22] M. Klanjˇsek, D. Arˇcon, A. Sans, P. Adler, M. Jansen,and C. Felser, Phys. Rev. Lett. , 057205 (2015).[23] M. Klanjˇsek, M. Horvati´c, S. Kr¨amer, S. Mukhopadhyay,H. Mayaffre, C. Berthier, E. Can´evet, B. Grenier, P. Le-jay, and E. Orignac, Phys. Rev. B , 060408 (2015).[24] M. Jeong, H. Mayaffre, C. Berthier, D. Schmidiger,A. Zheludev, and M. Horvati´c, Phys. Rev. Lett. ,106404 (2013).[25] M. Jeong, D. Schmidiger, H. Mayaffre, M. Klanjˇsek,C. Berthier, W. Knafo, G. Ballon, B. Vignolle,S. Kr¨amer, A. Zheludev, and M. Horvati´c, Phys. Rev.Lett. , 106402 (2016).[26] V. S. Zapf, D. Zocco, B. R. Hansen, M. Jaime, N. Har-rison, C. D. Batista, M. Kenzelmann, C. Niedermayer,A. Lacerda, and A. Paduan-Filho, Phys. Rev. Lett. ,077204 (2006).[27] A. Paduan-Filho, Brazilian Journal of Physics , 292(2012).[28] T. Giamarchi and H. J. Schulz, Phys. Rev. B , 325(1988).[29] M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S.Fisher, Phys. Rev. B , 546 (1989).[30] R. Yu, L. Yin, N. S. Sullivan, J. Xia, C. Huan, A. Paduan-Filho, N. F. Oliveira Jr, S. Haas, A. Steppke, C. F. Mi-clea, et al. , Nature , 379 (2012).[31] S. A. Zvyagin, J. Wosnitza, C. D. Batista, M. Tsukamoto,N. Kawashima, J. Krzystek, V. S. Zapf, M. Jaime, N. F.Oliveira, and A. Paduan-Filho, Phys. Rev. Lett. ,047205 (2007).[32] Using this value g = 2 .
31 for the g-factor for DTN, to-gether with the most frequently used set of couplings(
D, J , J ) = (8 . , . , .
18) K yields a second criticalfield in perfect agreement with the most precise estimates from NMR at h c = 12 .
32 T [21, 36].[33] H. J. Schulz, Phys. Rev. B , 6372 (1986).[34] N. Papanicolaou, A. Orend´aˇcov´a, and M. Orend´aˇc, Phys.Rev. B , 8786 (1997).[35] E. Wulf, D. H¨uvonen, R. Sch¨onemann, H. K¨uhne, T. Her-rmannsd¨orfer, I. Glavatskyy, S. Gerischer, K. Kiefer,S. Gvasaliya, and A. Zheludev, Phys. Rev. B , 014406(2015).[36] R. Blinder, M. Dupont, S. Mukhopadhyay, M. S. Grbi´c,N. Laflorencie, S. Capponi, H. Mayaffre, C. Berthier,A. Paduan-Filho, and M. Horvati´c, arXiv:1610.03312(2016).[37] C. Psaroudaki, J. Herbrych, J. Karadamoglou,P. Prelovˇsek, X. Zotos, and N. Papanicolaou, Phys.Rev. B , 224418 (2014).[38] C. J. Morningstar and M. Weinstein, Phys. Rev. D ,4131 (1996).[39] S. Capponi, A. L¨auchli, and M. Mambrini, Phys. Rev.B , 104424 (2004).[40] A. Abragam, Principles of Nuclear Magnetism (Claren-don Press, Oxford, UK, 1961).[41] A. Luther and I. Peschel, Phys. Rev. B , 3908 (1975).[42] T. Hikihara and A. Furusaki, Phys. Rev. B , 134438(2001).[43] T. Hikihara and A. Furusaki, Phys. Rev. B , 064427(2004).[44] D. Schmidiger, P. Bouillot, S. M¨uhlbauer, S. Gvasaliya,C. Kollath, T. Giamarchi, and A. Zheludev, Phys. Rev.Lett. , 167201 (2012).[45] S. Sachdev, Quantum Phase Transitions , 2nd ed. (Cam-bridge University Press, 2011).[46] P. Bouillot, C. Kollath, A. M. L¨auchli, M. Zvonarev,B. Thielemann, C. R¨uegg, E. Orignac, R. Citro,M. Klanjˇsek, C. Berthier, M. Horvati´c, and T. Gia-marchi, Phys. Rev. B , 054407 (2011).[47] P. Bouillot, Statics and dynamics of weakly coupled anti-ferromagnetic spin-1/2 ladders in a magnetic field , Ph.D.thesis, University of Geneva (2011).[48] T. Suzuki and S.-i. Suga, Phys. Rev. B , 172410 (2006).[49] T. Jolicœur and O. Golinelli, Phys. Rev. B , 9265(1994).[50] M. Troyer, H. Tsunetsugu, and D. W¨urtz, Phys. Rev. B , 13515 (1994).[51] J. Sagi and I. Affleck, Phys. Rev. B , 9188 (1996).[52] J. Kishine and H. Fukuyama, J. Phys. Soc. Jpn. , 26(1997).[53] S. Sachdev and K. Damle, Phys. Rev. Lett. , 943(1997).[54] B. D´ora, M. Gul´acsi, F. Simon, and H. Kuzmany, Phys.Rev. Lett. , 166402 (2007).[55] G. Vidal, Phys. Rev. Lett. , 040502 (2004).[56] U. Schollw¨ock, Annals of Physics , 96 (2011).[57] F. Verstraete, J. J. Garc´ıa-Ripoll, and J. I. Cirac, Phys.Rev. Lett. , 207204 (2004).[58] N. Laflorencie, Physics Reports , 1 (2016), quantumentanglement in condensed matter systems.[59] C. Karrasch, J. H. Bardarson, and J. E. Moore, NewJournal of Physics , 083031 (2013).[60] T. Barthel, New Journal of Physics , 073010 (2013).[61] ITensor library, http://itensor.org .[62] K. Fabricius and B. M. McCoy, Phys. Rev. B , 8340(1998).[63] J. Sirker, Phys. Rev. B , 224424 (2006). [64] S. Lukyanov and A. Zamolodchikov, Nuclear Physics B , 571 (1997).[65] S. Lukyanov, Phys. Rev. B , 11163 (1999).[66] I. Affleck, D. Gepner, H. J. Schulz, and T. Ziman, J.Phys. A: Math. Gen. , 511 (1989).[67] S. Eggert, I. Affleck, and M. Takahashi, Phys. Rev. Lett. , 332 (1994).[68] I. Affleck, J. Phys. A: Math. Gen. , 4573 (1998).[69] V. Barzykin, J. Phys.: Condens. Matter , 2053 (2000).[70] K. Fabricius, U. L¨ow, and J. Stolze, Phys. Rev. B ,5833 (1997).[71] T. Moriya, Progress of Theoretical Physics , 23 (1956). [72] See Ref. 36 for TLL parameters dependence on the mag-netic field.[73] M. Sato, T. Momoi, and A. Furusaki, Phys. Rev. B ,060406 (2009).[74] M. Sato, T. Hikihara, and T. Momoi, Phys. Rev. B ,064405 (2011).[75] E. Coira, P. Barmettler, T. Giamarchi, and C. Kollath,Phys. Rev. B , 144408 (2016).[76] S. Katsura, T. Horiguchi, and M. Suzuki, Physica ,67 (1970).[77] J. Stolze, A. N¨oppert, and G. M¨uller, Phys. Rev. B52