Temperature dependence of the NMR relaxation rate 1/T 1 for quantum spin chains Maxime Dupont, * Sylvain Capponi, † and Nicolas Laflorencie ‡ Laboratoire de Physique Théorique, Université de Toulouse and CNRS, UPS (IRSAMC), F-31062, Toulouse, France (Dated: October 13, 2016) We present results of numerical simulations performed on one-dimensional spin chains in order to extract the so-called relaxation rate 1/T1 accessible through NMR experiments. Building on numerical tensor network methods using the Matrix Product States (MPS) formalism, we can follow the non-trivial crossover occurring in critical chains between the high-temperature diffusive classical regime 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 on a generic spin-1/2 XXZ chain which is a paradigm of gapless TLL, as well as a more realistic spin-1 anisotropic chain, modelling the DTN material, which can be either in a trivial gapped phase or in a TLL regime induced by an external magnetic field. Thus, by monitoring the finite temperature crossover, we provide quantitative limits on the range of validity of TLL theory, that will be useful when interpreting experiments on quasi one-dimensional materials.

I. INTRODUCTION

One-dimensional (1d) quantum systems are known to be very peculiar due to strong quantum fluctuations that prohibit long-range order and can give rise to unusual phases of matter. In this context, it is remarkable that quantum spin chains fall generically into two classes regarding their low-energy properties [1,[START_REF] Giamarchi | Quantum Physics in One Dimension[END_REF]: (i) critical behavior where gapless low-energy excitations can be described in the framework of Tomonaga-Luttinger liquid (TLL) theory; (ii) gapped behavior.

Nevertheless, condensed-matter experiments are mostly done on quasi -1d materials, hence the role of small inter-chain couplings (as compared to the dominant 1d energy scale J 1d ) may become important at low-enough temperature (eventually leading to magnetic ordering). Conversely, at high temperature (T J 1d ), quantum fluctuations vanish so that a classical picture emerges. As a consequence, for realistic experimental systems the validity of a universal 1d TLL regime is not granted and should be checked in some unbiased way. In particular, understanding the intermediate temperature regime T ∼ J 1d , highly relevant to understand several experimental data, is a great theoretical challenge regarding dynamical observables.

In this paper, we focus on nuclear magnetic resonance (NMR) for quantum spin systems [START_REF] Horvatić | High Magnetic Fields[END_REF], and more specifically on the 1/T 1 spin-lattice relaxation rate. Indeed, this quantity contains lots of information on the dynamical properties of the system since it is directly related to dynamical spin-spin correlations. Moreover, being a local quantity (a crucial property of NMR technique), we will argue that reliable data can be obtained even though we will simulate finite spin chains.

Being of fundamental interest, the low-T behavior of the NMR relaxation rate has been investigated for several 1d or quasi-1d quantum magnets. Spingapped compounds, such as two-leg ladders SrCu 2 O 3 [4], BiCu 2 PO 6 [5], Sr 14-x Ca x Cu 24 O 41 [6], weakly coupled Haldane chains Y 2 BaNiO 5 [7], or dimerized spin chains AgVOAsO 4 [8], exhibit an activated relaxation at low-T . For gapless Heisenberg chain systems, the low-energy critical behavior has been studied [9][10][11][12] for Sr 2 CuO 3 which is an almost ideal realization with a large J 1d ∼ 2000 K and much smaller 3d couplings so that Néel temperature is pushed down to T N 5 K. For such an SU(2) symmetric material, a careful comparison of experimental and numerical NMR data has shown the prominent role of logarithmic corrections [13].

Another route to TLL behavior is to apply an external magnetic field on gapped materials such as spin-1 Haldane gap compound [14] (CH 3 ) 4 NNi(NO 2 ) 3 or dimerized spin-1/2 chains [15]. For such systems, a theoretical analysis of the 1/T 1 behavior has been performed in Refs. 16 and 17.

A useful experimental review on NMR properties of several spin chains can be found in Ref. 18. Note also that 1/T 1 measurements have also been used to characterize one-dimensional metallic phase in carbon nanotube [19] or quasi-1d superconductor [20].

More recently, interesting quasi-1d spin-gapped materials have also been investigated using NMR [21]: an anisotropic spin-1 system NiCl 2 -4SC(NH 2 ) 2 (DTN) and a spin-ladder one (C 5 H 12 N) 2 CuBr 4 (BPCB). In both cases, 1/T 1 measurements could be interpreted either as coming from magnon (respectively spinon) excitations in the gapped (respectively gapless) 1d phase, and the quantum critical regime was also argued to be universal. Most importantly, the whole temperature range, including 1d as well as 3d regimes, was discussed.

Experimentally, when decreasing temperature, the NMR relaxation rate 1/T 1 has been found to diverge in the TLL regime, with power-law governed by a characteristic exponent. Such an analysis is used in experiments to determine the corresponding TLL exponent K [22,23]. For example, it was a smoking-gun signature of attractive TLL in (C 7 H 10 N) 2 CuBr 4 (DIMPY) com-pound [24,25]. However, given that we are generically dealing with quasi-1d materials, critical fluctuations and 3d ordering will limit the low-energy 1d regime, and a genuine TLL critical behavior is observable only within some finite window in temperature. This remains to be analyzed more quantitatively, which is the main purpose of this work.

The rest of the paper is organized as follows. In Sec. II, we present the theoretical models and provide useful definitions. Section III describes the numerical technique based on finite temperature Matrix Product States (MPS) approach. Results are then discussed in Sec. IV. Finally, we present our conclusion in Sec. V.

II. MODELS AND DEFINITIONS

We give in this section the two models that will be studied in this paper and a small discussion on their phase diagram. Both models present a TLL gapless phase and a gapped phase, induced by an external magnetic field. We will also provide definition of the NMR relaxation rate 1/T 1 and discuss its expected behavior with temperature.

A. Theoretical models

The spin-1/2 XXZ chain

We first consider one of the simplest paradigmatic example of TLL liquid, namely the spin-1/2 XXZ chain Hamiltonian:

H XXZ = J L-1 j=1 S x j S x j+1 + S y j S y j+1 + ∆S z j S z j+1 -gµ B h L j=1 S z j (2.1)
where ∆ ∈ (-1, 1] denotes the Ising anisotropy, J the coupling strength and h is an applied magnetic field in the z direction with g the gyromagnetic factor and µ B the Bohr magneton constant. The Hamiltonian is defined with open boundary conditions (OBC), as will be used in our numerical simulations.

In the range ∆ ∈ (-1, 1] the XXZ model can be described by a TLL as long as its spectrum remains gapless [START_REF] Giamarchi | Quantum Physics in One Dimension[END_REF]. As a function of magnetic field, the gapless regimes extends up to a critical field gµ B h c = J(∆ + 1), and the system becomes gapped for h > h c . In the latter regime, the gap increases linearly with the applied magnetic field, ∆ g = gµ B (hh c ), see Fig. . There is a gapless TLL phase below hc = J(∆ + 1)/gµB and a gapped one above when the system is fully polarized. (ii) Lower panel : magnetization curve of the 1d DTN Hamiltonian (2.2) as a function of the magnetic field. Two gapped phases (large-D and polarized) are respectively located below hc1 and above hc2. The intermediate gapless phase can be described by a TLL theory.

The quasi-1d spin-1 compound "DTN"

We also discuss a quasi-1d magnetic insulator compound NiCl 2 -4SC(NH 2 ) 2 , also called DTN, whose relevant 3d structure consists of weakly coupled S = 1 chains in the two other transverse (with respect to the chain axis) directions. Its experimental interest comes from the appearance of a Bose-Einstein condensation (BEC) phase when applying a magnetic field at low temperature [26,27]. More recently, Br-doped (disordered) DTN was suggested to be a good experimental candidate for observing a Bose glass phase [28][29][30].

Although there is 3d magnetic order observed below T N ∼ 1 K in DTN [26] due to weak inter-chain couplings along the two transverse directions, J 3d /J 1d 0.08, one expects 1d physics and a TLL regime at higher T . The effective Hamiltonian to describe this situation reads

H DTN-1d = J L-1 j=1 S j •S j+1 + L j=1 D S z j 2 -gµ B hS z j , (2.2)
where S j = (S x j , S y j , S z j ) are spin-1 operators. In the current literature, [31] J = 2.2 K is the 1d antiferromagnetic coupling and D = 8.9 K is the single-ion anisotropy. The magnetic field h is given in Tesla with g = 2.31 [START_REF]Using this value g = 2.31 for the g-factor for DTN, together[END_REF].

The phase diagram of this 1d Hamiltonian (2.2) is sketched in Fig. 1. In the absence of magnetic field, due to the large on-site anisotropy D, the system is in the so-called large-D phase [START_REF] Schulz | [END_REF]. This is a trivial phase, adiabatically connected to the product state | → . . . → where each state is in a non-magnetic S z = 0 eigenstate. Clearly, this phase has a finite spin-gap, which corresponds to the first critical field h c1 needed to magnetize the system. Its value is known to be, at first order in J/D 1 [34]: h c1 /gµ B = D -2J + O(J 2 /D) 3 T. At finite magnetic field there is a gapless TLL regime for h ∈ [h c1 , h c2 ], with h c2 /gµ B = D + 4J = 11.40 T. Above this critical saturation field, the system becomes gapped again, entering a fully polarized phase. As a side remark, we recall that in the true 3d material DTN, both critical fields are shifted due to interchain couplings, so that h c1 = 2.10(5) T [35], and h c2 = 12.32 T [36].

In the TLL phase and close to the upper critical field h c2 , the DTN Hamiltonian (2.2) can be mapped toward an effective XXZ model of spins S = 1/2 (2.1). Using perturbation theory, effective parameters [21,37] are given by J = 2J and ∆ = 0.5. This result can be refined using 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 effective magnetic field h = h -J -D.

B. Relaxation rate 1/T1

The nuclear spin-lattice relaxation rate T -1 1 measured by NMR [START_REF] Abragam | Principles of Nuclear Magnetism[END_REF] is basically testing the local and dynamical spin correlation function S aa nn (ω 0 ) with a = x, y, z, and ω 0 being the NMR frequency at a given site n,

1 T 1 = γ 2 2 A 2 ⊥ [S xx nn (ω 0 ) + S yy nn (ω 0 )] + A 2 S zz nn (ω 0 ) = 1 T ⊥ 1 + 1 T 1 (2.3)
Here, A ⊥ is the transverse hyperfine coupling constant, A the longitudinal one, γ the gyromagnetic ratio and

S ab nn (ω) = Re 2 ∞ 0 dt e iωt S a n (t) S b n (0) -S a n (t) S b n (0) , (2.4)
with S a = S b † , S a (t) = e iHt S a e -iHt , and is the thermal average defined later in (3.4). For convenience, the x and y spin components can be expressed using the raising and lowering operators,

S xx nn (ω) + S yy nn (ω) = 1 2 S +- nn (ω) + S -+ nn (ω) . (2.5)
It is theoretically justified to take the limit ω 0 → 0 since the NMR frequency is of a few tens or hundreds of MHz, corresponding to temperatures of the order of mK, thus being the smallest energy scale in the problem. Indeed, such temperatures are neither reached in experiments of interest nor in our numerical simulations, in particular for our purpose of probing the finite temperature TLL regime in quasi-1d systems. As a side remark, the two correlations S +- nn (ω 0 ) and S -+ nn (ω 0 ) become equivalent in this limit ω 0 → 0.

The weight of the transverse and longitudinal contributions in the relaxation rate 1/T 1 is experimentally governed by the hyperfine coupling tensors A ⊥ and A . To favor one over the other, a specific nucleus can be targeted for the NMR experiment. In the case of DTN, proton 1 H (nuclear spin I = 1/2) probes both components while nitrogen 14 N (I = 1) probes dominantly the transverse one [36].

Hyperfine coupling tensors put aside or set equal to one, the low temperature behavior of the transverse and longitudinal components of 1/T 1 depends on microscopic parameters of the model. At low temperature, the transverse component is larger than the longitudinal one and it is thus justified to consider 1/T 1 1/T ⊥ 1 . In the following we will fix γ 2 A 2 ⊥, = 1, and compute the relaxation rates

1 T ⊥ 1 = ∞ 0 dt Re S ± n (t) S ∓ n (0) , (2.6) and 1
T 1 = ∞ 0 dt Re [ S z n (t) S z n (0) -S z n (t) S z n (0) ] .
(2.7)

Note that the single operator averages are time independent, S z n (t) = S z n (0) .

C. Tomonaga-Luttinger liquid description in the gapless phase

The 1d gapless phase, experimentally accessible by tuning a control parameter such as the external magnetic field, can be effectively described by the TLL Hamiltonian [START_REF] Giamarchi | Quantum Physics in One Dimension[END_REF],

H TLL = 1 2π dr uK [∂ r θ (r)] 2 + u K [∂ r φ (r)] 2 (2.8)
where u is the velocity of excitations and K is the dimensionless TLL parameter. They both fully characterize the low-energy properties of the system and are thus modeldependent. θ(r) and φ(r) are bosonic fields obeying the commutation relation [φ(x), θ(y)] = iπδ(xy).

For the XXZ model (2.1), the TLL parameters K and u are known from Bethe ansatz equations [START_REF] Luther | [END_REF]. At zero magnetic field, analytical expressions are known as a function of the Ising anisotropy ∆:

K = π 2 arccos (-∆) and u J = π √ 1 -∆ 2 2 arccos ∆ . (2.9)
For generic non-integrable models, the TLL parameters can be obtained numerically using DMRG by fitting static correlation functions [42,43] which has been successfully done in the past years for various quasi-1d compounds [17,44].

In the TLL framework, the dynamical correlation S ab nn (ω) defined in (2.4) can be computed analytically as a 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 defined as,

χ ab ij (t) = -iΘ (t) S a i (t) , S b j (0) , (2.10)
where Θ(t) is the Heaviside function. Since a typical relevant case in experiments is the local quantity i = j = n, we only consider this case in the following. In frequency space, the susceptibility can be related to the dynamical spin correlation function by [START_REF] Sachdev | Quantum Phase Transitions[END_REF] S ab nn (ω) = 2 e -βω -1 Im χ ab nn (ω) .

(2.11)

In the low energy limit βω 1 an analytical expression for Eq. (2.11) can be obtained, leading to [17, 46,

⊥ 1 = 2A x cos π 4K u 2πT u 1 2K -1 B 1 4K , 1 - 1 2K 47] 1 T
, (2.12) and 1

T 1 = A z cos (πK) 2u 2πT u 2K-1 B (K, 1 -2K) + KT 4πu 2 , (2.13)
with B(x, y) the Euler beta function and A x,z prefactors of the static correlation functions. Thus, generically 1/T ⊥ 1 (T ) diverges at zero temperature as a K-dependent power-law, and dominates over 1/T 1 . Note that for finite magnetic field (equivalent to non half-filled case), additional subleading corrections are expected [START_REF] Suzuki | [END_REF].

D. Gapped regime behavior

In contrast to a TLL gapless phase where 1/T 1 has power-law behavior at low temperature, we also consider the gapped regime where fluctuations are exponentially suppressed [16,[49][50][51][52][53], such that for

T < ∆ g 1 T ⊥, 1 ∝ exp -α ⊥, ∆ g /T , (2.14)
where ∆ g is the energy gap of the system, and α ⊥, an O(1) prefactor which depends on the relaxation processes, and also on the temperature range [52]. Below, in section IV B 2, we show numerical results for the highfield gapped regime of the XXZ chain where our data are perfectly described by α ⊥, = 1. At higher temperature above the gap, one may also expect a non-trivial crossover to TLL regime [54].

III. NUMERICAL METHODS

To get the relaxation rate one needs to obtain in the first place the dynamical correlation S a n (t) S b n (0) . We used the TEBD (Time Evolution Block Decimation) algorithm [55] with both real/imaginary time through the MPS formalism adapted for 1d systems [56]. A general one-dimensional system containing L sites with OBC can be represented by the following MPS,

|Ψ = {si} A s1 a1 A s2 a1a2 • • • A s L a L-1 |s 1 |s 2 • • • |s L (3.1)
where the local index s i is the physical index representing an element of the local Hilbert space at site i. Its dimension is d and is equal to 2 (↑ and ↓) for spin-1/2 or 3 (↑, ↓ and →) for spin-1. We note a i the bond index whose dimension is directly related to the "number of states" m to describe the system, meaning that m is a control parameter in the numerical simulations. The first step is to perform an imaginary time evolution on the system to reach the desired temperature. Once the state is at hand, the second step consists of evolving it through a real time evolution. At each time step, the correlation is measured. When all the data in time-space have been obtained, a numerical Fourier Transform can be performed to get the data in frequencyspace.

A. Time evolution with MPS

We will be general and consider the case with a Hamiltonian 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 MPS up to a time t. The operation can be discretized using smaller time steps τ such that t = N τ , leading to e -itH = N e -iτ H . If the time step τ is small enough, a first (or higher) order Trotter decomposition can be performed,

e -itH N e -iτ Heven N e -iτ H odd + O(τ 2 ) (3.2)
where H even and H odd respectively correspond to the even and odd bond Hamiltonians only acting on two nearest-neighbor spins. The decomposition is possible because even -odd -bond Hamiltonians commute with each others. But it is not exact and leads to an error in τ due to the fact that [H even , H odd ] = 0. The advantage is that the bond Hamiltonians can easily be diagonalized and exponentiated since they are only d 2 × d 2 matrices. And so, applying successive evolution gates on the MPS as well as singular value decompositions to restore the MPS original tensor-site dependent form (3.1) will eventually lead to a time-t evolved state.

B. Finite temperature with MPS

It was useful to introduce time-evolution concepts also to discuss finite-temperature with MPS [57]. The main idea is to represent the density matrix ρ β of the physical (mixed) state in an artificially enlarged Hilbert space as a pure state |Ψ β -which is what we can deal with in the MPS formalism. The auxiliary space can simply be constructed as a copy of the original one.

Assuming that we know the purification of the density matrix ρ β=0 as a wave function |Ψ β=0 it can be shown that an imaginary time evolution has to be performed over the infinite temperature state in order to get the finite temperature state,

|Ψ β = e -βH/2 |Ψ β=0 (3.3)
with the Hamiltonian only acting on the physical sites. This imaginary time evolution can be performed using the TEBD algorithm described before in III A. Expectation values can then be measured at inverse temperature,

O β = Tr Oe -βH Tr [e -βH ] = Ψ β |O|Ψ β Ψ β |Ψ β (3.4)
For this procedure to work, the initial state |Ψ β=0 has to be a product state of Bell states between each physical site and its associated auxiliary site,

|Ψ β=0 = 1 √ N L n=1 {s} |p n s a n s (3.5)
with |p corresponding to physical sites, |a to auxiliary ones and N a normalization constant. The summation is over the d possible local states s. Such a state is simple enough to be built exactly in the MPS formalism.

C. Numerical limitations

The main limitations are about the temperature and the final time one can reach using the methods described above. The reason in both cases is directly related to a rapid growth in the entanglement entropy while evolving the state. This implies to keep larger and larger number 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 to reach low temperatures. Indeed, one expects a volumelaw entanglement entropy (i.e. linear with the system size L) due to the auxiliary sites which are used to purify the thermal state. As a consequence, the number of kept states m needed to describe accurately the system will grow exponentially as the temperature T decreases.

On the other hand, the maximal (real) time that can be reached is of the order of few tens of J -1 typically, for similar reasons as discussed above, namely the linear growth of entanglement entropy with time [58]. Thus, fixing a maximum number of kept states m limits simulations to a finite time t max . Note that this limitation applies at all temperatures, even T = 0.

Despite these severe limitations, recent progress in the field has allowed some improvements. For instance, we will make use of the auxiliary degrees of freedom which are used to purify the thermal state by time-evolving them with -H, which is mathematically exactly the same but has been shown to improve substantially the time range [59]. By construction, this trick only applies to finite-temperature simulations though. Last, although its use did not prove to be systematically reliable in our case, we would like to mention the possibility to use so-called linear prediction technique, coming from data analysis [60] which aims at predicting "longer time" behavior from the knowledge of dynamical correlations at "intermediate time".

IV. RESULTS

We provide in this section our numerical results [START_REF]ITensor library[END_REF] about the relaxation rate 1/T 1 using models and techniques presented in the previous sections II and III. First of all we will focus on the XX model (equivalent to free fermions) for which we can compute exactly the dynamical correlations for all temperatures and that will serve as a benchmark for our simulations. Next, we will turn to 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 XX model, can be mapped onto a model of free spinless fermions using a Jordan-Wigner transformation. We restrict ourselves to h = 0. It can be diagonalized in Fourier space with ε k = J cos k and k = nπ L+1 with n = 1, 2, . . . , L considering open boundary conditions,

H XX = J 2 L-1 j=1 c † j c j+1 + h.c. = 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 results presented in this section will be valid for all regimes. We present the details of the calculations in appendix A for the analytical exact expressions of the dynamical correlations for the XX model using precisely the same conditions as in our numerical simulations (i.e. a finite chain length with open boundary conditions).

We show the 'bare' results in Fig. 2 that will be used to obtain the relaxation rate thereafter. Ideally, one is interested in the thermodynamic limit (i.e. L → ∞) but we see that, at finite temperature (hence finite correlation length), working on finite length chains with only a moderate number of sites L allows to get reliable data. Indeed, the MPS estimates agree perfectly with the exact expressions (see appendix A).

First of all, we consider the local dynamical correlation of the site in the middle of the chain reducing de facto boundary effects. Then, as finite size effects are known to be caused by the reflection of the propagating excitations at TLL velocity u on the boundary of the system, one can estimate a time below which the dynamical correlations can 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 the inverse of the NMR frequency ω 0 ) and get reliable values of 1/T ⊥ 1 by integrating over time. Moreover, we have also checked that finite size effects are extremely small since we are computing a local correlation.

The same cannot be said for the longitudinal correlations. They continue to oscillate even for high temperatures and long times, and their amplitude gets (very) slowly smaller with time. This implies severe limitations to 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 -4 at t = 1000 J -1 . This makes the value of 1/T 1 very difficult to estimate. This well-known behavior is related to spin diffusion-like behavior [START_REF] Fabricius | [END_REF]63] which cause a logarithmic divergence at small frequency ω. However, we have to remember that the NMR frequency ω 0 eventually provides a natural cutoff.

For completeness we display the real part of the Fourier transform on the right panels of Fig. 2 for which the 1/T 1 value as defined in section II B corresponds to the ω = 0 value.

B. Spin-1/2 XXZ chain at ∆ = 0

Gapless regime

Building on the perfect agreement observed previously between MPS estimates and the exact analytical solution of the XX model, we are now confident to extend our study of the more generic XXZ case -1 < ∆ ≤ 1, described by a TLL, and compute the relaxation rates. Results are plotted in Fig. 3 for various values of the anisotropy. The simulations were performed on systems of size L = 64 with a cutoff of ε = 10 -10 in the singular values. We kept a maximal number of D = 500 states. A fourth order Trotter decomposition was used with a Trotter step of τ = 0.1.

First, in the gapless regime we do observe an excellent quantitative agreement between numerical estimates and the TLL prediction Eq. ( 2.12) at low-enough temperature. This asymptotic regime with a power-law behavior ∼ T 1 2K -1 occurs only below T /J ∼ 0.1 -0.2 (depending on the anisotropy ∆). Here we stress that there are no free parameters in the analytic expressions. Indeed, the reduced temperature T /J for the spin-1/2 XXZ chain at various ∆ and h = 0 obtained numerically using MPS techniques (circles, from top to bottom: ∆ = -0.8, -0.6, -0.5, -0.4, -0.2, 0.0, 0.2, 0.4, 0.5, 0.6, 0.8, 1). Numerics are compared to TLL theory Eq. ( 2.12) at low temperature (thick lines) for |∆| < 1, and with Eq. (4.2) for the SU(2) Heisenberg point ∆ = 1. The thin lines between the circles are guides to the eyes.

TLL parameters are computed using the exact expressions Eq. (2.9) for u and K, and A x is obtained following Refs. 64 and 65. The isotropic limit ∆ = 1 is a special point where logarithmic corrections appear in several quantities [66][67][68], leading to a very slow divergence of the (isotropic) NMR relaxation rate [12,69] 1

T 1 1 √ 2π 3 ln Λ T + 1 2 ln ln Λ T , (4.2)
where Λ 24.27J. MPS estimates compare well with this parameter-free expression, as visible in Fig. 3.

Interestingly, we notice the non-monotonic behavior of 1/T ⊥ 1 with temperature only when ∆ 0 (which corresponds to repulsive or vanishing interactions in the fermionic language).

As a last comment, we have observed that for infinite temperature (β = 0), the value of 1/T ⊥ 1 does not depend on the sign of ∆, which is expected since the many-body spectrum of H ∆ is an odd function of ∆. Its value is minimum for ∆ = 0 with 1/T ⊥ 1 = √ π/(2J) [70] and increases with |∆|. At the isotropic point |∆| = 1 we expect the relaxation rate to diverge due to the diffusion-like behavior [START_REF] Fabricius | [END_REF]63] of the dynamical correlation function. Our results at infinite-T go beyond Baker-Campbell-Hausdorff expansion developed up to O(t 2 ) in Ref. 71 to compute S ± j (t)S ∓ j (0) at short times, which would suggest

1/T ⊥ 1 ∼ J -1 (1 + ∆ 2 ) -1 2
. This prediction is in contrast to what we found, namely the transverse relaxation rate increasing with |∆|. Indeed, while such an expansion finds 10 -1 plotted against reduced inverse temperature βJ for the spin-1/2 XXZ chain in its gapped phase for the anisotropy ∆ = 0.5. The critical magnetic field is hc = 3J/2gµB and the value of the gap ∆g = gµB (h -hc). Numerical results are obtained using MPS techniques (circles and diamonds) and the exponentially decaying behavior is verified with the straight lines set with the expected gap value 1/T ⊥, 1 = c ⊥, • e -β∆g and c ⊥, a non-universal free parameter.

10 0 1/T ⊥, 1 ∆ g = 0. 5 ∆ g = 1 .0 ∆ g = 1 . 5 h/gµ B = 2.0 J h/gµ B = 2.5 J h/gµ B = 3.0 J 1/T 1 1/T ⊥ 1 1/T 1 1/T ⊥ 1
the correct gaussian behavior for ∆ = 0 (free-fermions), higher-order terms have to be taken into account for |∆| > 0 where the transverse dynamical correlation function at longer times gets larger when increasing |∆|.

Gapped XXZ chain

We then set the anisotropy value to ∆ = 0.5 and apply a magnetic field to move into the gapped phase. Transverse and longitudinal relaxation rates 1/T ⊥, 1 are plotted in Fig. 4 where we observe an excellent agreement with an exponentially activated behavior ∼ exp(-β∆ g ), where ∆ g is the spin gap. We notice that as the gap gets smaller, the lower the temperature has to be to observe the exponential law.

C. DTN

We now move to the DTN compound in its 1d limit described by Eq. (2.2). We compute the relaxation rates for various values of the magnetic field h, mainly close to h c2 which is relevant for NMR experiments [21]. It is a more challenging system to simulate than the XXZ model as it is made of spins S = 1 (enlarged local Hilbert space). The simulations were performed on open chains of size L = 64 with a cutoff of ε = 10 -10 in the singular values. We kept a maximal number of D = 150 states. FIG. 5. (color online) Transverse relaxation rate 1/T ⊥ 1 plotted vs. temperature T for the spin-1 DTN chain obtained numerically using MPS techniques (circles). The low temperature behavior is compared to TLL prediction (straight lines). The magnetic field h is given in Tesla. The inset compares TLL prediction and numerical results for T = 0.4 K and covers the whole TLL phase from hc1 to hc2. The lower panel is a zoom on the low temperature asymptotic power-law regime.

A fourth order Trotter decomposition was used with a Trotter step of τ = 0.02.

Numerical results, shown in Fig. 5, compare extremely well with TLL prediction at low temperature. Interestingly, the TLL power-law behavior starts at slightly higher temperature, as compared to the XXZ model, T 0.5 K (T /J ∼ 0.2), especially as we approach the middle of the TLL phase, away from the critical field h c2 . We point out that there are again no adjustable coefficients, the TLL parameters being computed independently using standard DMRG [72]. The tiny difference that appears at low temperature between numerical data and TLL is due to the limited number of states m kept when performing calculations. Though this does not dispute the TLL prediction, it reveals the challenge in such time-dependent simulations. The inset in Fig. 5 shows the transverse relaxation rate at T = 0.4 K for various values of the magnetic field covering the whole range from h c1 to h c2 . Once more, there is a very good agreement between numerics and TLL theory except when one gets close to the critical fields. Indeed, as we clearly see in the lower panel of Fig. 5 for h = 11.0 T, the power law is not met yet for the lowest temperature we could reach T = 0.2 K.

The non-monotonic behavior of 1/T ⊥ 1 observed in the XXZ model is absent for the DTN and may seem odd at first place since it can be mapped effectively onto a S = 1/2 XXZ chain with ∆ = 0.5 or 0.36 and could thus be compared with Fig. 3. However this non-monotonic variation is observed at high temperature while this mapping is only justified in the low-energy limit as discussed in II A 2.

One can also try to compare the relaxation rates of Fig. 5 with the NMR data for the DTN compound given in Ref 21. What draws our attention is the nonmonotonic regime of 1/T ⊥ 1 observed at high temperature experimentally, which, as we have just discussed, is not theoretically predicted for a single DTN chain. Yet it cannot be attributed to 3d effects as J 3d = 0.18 K is very small compared to the temperature T . We then observed that experiments are performed by proton ( 1 H) NMR which probes both 1/T ⊥ 1 and 1/T 1 . We therefore interpret this effect as due to the parallel contribution of the relaxation rate. We show in Fig. 6 both the transverse and longitudinal 1/T ⊥, 1 as a function of temperature. We cannot precisely estimate the value of 1/T 1 due to its dependence on ω 0 (and therefore on our maximum time in numerical simulations) so that we give a lower bound. Its high temperature contribution to the total relaxation rate clearly dominates over the transverse part and explains well the experimental nonmonotonic regime at high T .

Perhaps more importantly, as displayed in Fig. 6, the 3d BEC ordering observed in DTN [26,36] for m z 0.85 at T N 0.59 K occurs above the asymptotic regime where the genuine TLL power-law behavior is expected. It is therefore impossible to directly extract TLL exponents in DTN, because of interchain effects that eventually lead to an ordering of the coupled TLLs. Ideally we would expect for quasi-1d systems the TLL description of the NMR relaxation to be valid in the following temperature regime: J 1d T J 3d . Concerning the difficulty to obtain reliable data at high temperature for the longitudinal 1/T 1 , it is well known that this is due to spin diffusion-like behavior [START_REF] Fabricius | [END_REF]63]. Therefore, measurements should in principle depend explicitly on the NMR frequency ω 0 .

V. CONCLUSION Performing time-dependent numerical simulations at finite temperature on 1d systems to compute the NMR relaxation rate 1/T 1 , we have discussed the temperature range validity of analytical predictions for two models (i) the paradigmatic example for Tomonaga-Lutinger liquids: the spin-1/2 XXZ chain for various Ising anisotropies, and (ii) a more realistic S = 1 Hamiltonian, relevant for experiments on the DTN compound relaxation rates for the DTN spin-1 chain at h = 11.0 T, corresponding to m z 0.85. As 1/T 1 cannot be estimated for sure, we only provide a lower bound. The non-monotonic behavior observed experimentally at high T in Ref 21 apparently comes from the large contribution of 1/T 1 at high temperature. Experimental data for DTN [21] at the same magnetization are shown for comparison, after a proper rescaling in order to match the low-T regime. The 3d BEC transition temperature TN (m z 0.85) 0.59 K [36] is also shown.

as a function of an external magnetic field.

Both models present in some regime a gapless phase that can be described by TLL "low-energy" theory, with a relaxation rate dominated by its transverse component 1/T ⊥ 1 ∼ T 1 2K -1 algebraically diverging at low temperature, where K is the dimensionless TLL exponent. We observed that the expected power-law behavior occurs only below T /J 1d ∼ 0.1 -0.2, thus defining the lowenergy limit of validity of TLL theory an order of magnitude below the energy scale J 1d of the system. It is important to be able to define this limit as TLL predictions are often used experimentally on quasi-1d compounds to extract the value of K. As a consequence, we believe that it remains experimentally challenging [15], and often impossible, to explore a genuine critical 1d regime in quasi-1d compounds when J 1d is small and 3d ordering prevents a wide TLL regime. For instance, we have shown that for DTN, the BEC ordering temperature is larger than the crossover temperature towards TLL behavior.

We have also studied the transverse relaxation rates of these two models in other regimes than TLL theory. First, we considered high temperatures, with a peculiar non-monotonic behavior in the S = 1/2 XXZ model in the repulsive regime at high T , which does not exist for the 1d S = 1 model of DTN. However, such a nonmonotonic dependence with temperature at high T is experimentally observed in DTN. We showed that this effect comes from the parallel contribution of the relax-ation rate 1/T 1 dominating at high temperature over the transverse part. Finally, we verified that in gapped phases the relaxation rates are exponentially suppressed ∼ exp (-∆ g /T ) and can indeed lead to accurate determinations of the spin gap, at least in a regime of temperature T ∆ g since other relaxation mechanisms can change the activated behavior at higher temperature.

We want to emphasize again the role of 3d ordering at finite temperature, preventing the observation of a 1d TLL regime. As discussed for the particular case of DTN, one needs a hierarchy of energy scale J 1d T J 3d to be able to directly extract the TLL exponent K from the divergence of T -1 1 with T .

Concerning future advances for quasi-1d systems, we can envision trying to simulate imaginary-time dynamics using quantum Monte-Carlo techniques, provided that the model has no minus-sign problem. While it will be necessary to perform a numerical analytic continuation (using for instance Maximum Entropy techniques), we have some hope that this could lead to reliable results to quantitatively capture the influence of interchain effects. For 1d chains, it has been rather successful [10]. As a matter of fact, our extensive 1d results could serve as useful benchmarks for that too.

While we have considered various chains, we are far from being exhaustive. Indeed, there are some TLL models for which elementary excitations may not be simple spin flips, for instance multipolar nematic phase for which 1/T 1 behavior will be different [73,74]. In a similar line of thought, we could imagine simulating more complicated models including charge and spin degrees of freedom to describe NMR relaxation in metallic or superconducting wires.

Note added : While completing this work, a related numerical study by Coira et al. has appeared [75]. Our results are perfectly compatible with each other when comparison can be made, such as the transverse 1/T 1 data for a single spin-1/2 XXZ chain with ∆ ≥ 0. . (A5)

Transverse correlations

The transverse correlations S ± i (t) S ∓ j (0) have a more complicated structure in the fermion representation due to the string of operators in the exponential [77]. We introduce the following identity e iπc † l c l = (-1) c † l c l = A l B l with A l = c † l + c l and B l = c † l -c l , leading to 2 S ± i (t) S ∓ j (0) = i-1 l=1

A l (t)B l (t) A i (t) j-1 l=1

A l (0)B l (0) A j (0) . (A6)

Now thanks to Wick's theorem this product of many fermion operators can be rewritten as elementary expectation values of two operators through the Pfaffian of some skew-symmetric matrix. Its elements above the diagonal (which is purely made of zeros) are And equivalently B i (t)B j (0) = -A i (t)A j (0) as well as B i (t)A j (0) = -A i (t)B j (0) .

A 1 (t)

1

 1

  .

  FIG. 1. (color online) (i) Upper panel : magnetization curve of the XXZ Hamiltonian (2.1) as a function of the magnetic field for ∆ ∈ (-1, 1]. There is a gapless TLL phase below hc = J(∆ + 1)/gµB and a gapped one above when the system is fully polarized. (ii) Lower panel : magnetization curve of the 1d DTN Hamiltonian (2.2) as a function of the magnetic field. Two gapped phases (large-D and polarized) are respectively located below hc1 and above hc2. The intermediate gapless phase can be described by a TLL theory.

FIG. 2 .

 2 FIG. 2. (color online)We compare numerical results (circles) used to determine the 1/T1 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 the longitudinal ( ) 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 for the lowest temperature for the longitudinal correlations. Indeed, this is a priori the hardest to compute and thus the most subject to errors. The right panel shows the real part of the Fourier transform of the real time data. Although we only show data up to t = 40 J -1 the Fourier Transform of the exact zz correlations was performed using data up to t = 1000 J -1 .

FIG. 3 .

 3 FIG. 3. (color online) Transverse relaxation rate 1/T ⊥ 1 vs.reduced temperature T /J for the spin-1/2 XXZ chain at various ∆ and h = 0 obtained numerically using MPS techniques (circles, from top to bottom: ∆ = -0.8, -0.6, -0.5, -0.4, -0.2, 0.0, 0.2, 0.4, 0.5, 0.6, 0.8, 1). Numerics are compared to TLL theory Eq. (2.12) at low temperature (thick lines) for |∆| < 1, and with Eq. (4.2) for the SU(2) Heisenberg point ∆ = 1. The thin lines between the circles are guides to the eyes.

FIG. 4 .

 4 FIG. 4. (color online) Transverse and longitudinal relaxation rates 1/T ⊥, 1

FIG. 6 .

 6 FIG. 6. (color online) Longitudinal 1/T 1 and transverse 1/T ⊥ 1

1 . 1 Zsin ki sin kj 1 + tanh βε k 2 e 4 L + 1 k sin 2

 11242 Longitudinal correlationsThe longitudinal correlations S z i (t) S z j (0) for the XX model are basically density correlations when performing the Jordan-Wigner transformation. Since we are interested in the temperature dependence of the correlations, what we need to compute is actually,S z i (t) S z j (0) = Tr e iHt S z i e -iHt S z j e -βH (A1)with Z = Tre -βH the partition function. Calculations lead to[76], -iε k t (A3) and I 2 ,I 2 = 2 L + 1 k sin ki sin kj 1 -tanh βε k 2 e +iε k t . (A4)Also, the single operator averages are S z i (t) = S z i (0) = 1 -ki e -βε k + 1

  B 1 (t) A 1 (t)A 2 (t) . . . A 1 (t)A j (0) B 1 (t)A 2 (t) . . . B 1 (t)A j (0) Tr e iHt A i e -iHt B j e -βH .

				. . .	. . .	(A7)
				B j-1 (0)A j (0)
	At finite temperature,
	A i (t)B j (0) =	1 Z	(A8)
	The two-body expectation values can be computed going
	in Fourier space,
	A i (t)A j (0) =	2 L + 1 βε k 2	(A9)
	and,
	A i (t)B j (0) =	2 L + 1 k	sin kj sin ki
		-sin ε k t -i cos ε k t tanh	βε k 2	. (A10)

k sin kj sin ki cos ε k ti sin ε k t tanh

ACKNOWLEDGMENT

The authors are grateful to M. Horvatić and M. Klanjšek for their careful reading of this manuscript and thoughtful comments and to R. Blinder, T. Giamarchi and E. Orignac for valuable discussions. This work was performed using HPC resources from GENCI (Grant No. x2015050225 and No. x2016050225), and is supported by the French ANR program BOLODISS and Région Midi-Pyrénées.

Appendix A: Dynamical correlations for the XX model

For completeness, we remind the reader the exact expressions for time-displaced spin correlations in the exactly solvable XX model.