Temperature dependence of the NMR spin-lattice relaxation rate for spin-1/2 chains

# Temperature dependence of the NMR spin-lattice relaxation rate for spin-1/2 chains

E. Coira DQMP, University of Geneva, 24 Quai Ernest Ansermet, 1211 Geneva, Switzerland    P. Barmettler Supercomputing Systems AG, Technoparkstrasse 1, 8005 Zürich, Switzerland    T. Giamarchi DQMP, University of Geneva, 24 Quai Ernest Ansermet, 1211 Geneva, Switzerland    C. Kollath HISKP, University of Bonn, Nussallee 14-16, 53115 Bonn, Germany
###### Abstract

We use recent developments in the framework of time dependent matrix product state method (tMPS) to compute the nuclear magnetic resonance (NMR) relaxation rate for spin-1/2 chains under magnetic field and for different Hamiltonians (XXX, XXZ, isotropically dimerized). We compute numerically the temperature dependence of the . We consider both gapped and gapless phases, and also the proximity of quantum critical points. At temperatures much lower than the typical exchange energy scale our results are in excellent agreement with analytical results, such as the ones derived from the Tomonaga-Luttinger liquid (TLL) theory and bosonization which are valid in this regime. We also cover the regime for which the temperature is comparable to the exchange coupling. In this case analytical theories are not appropriate but this regime is relevant for various new compounds with exchange couplings in the range of tens of Kelvin. For the gapped phases, either the fully polarized phase for spins chains or the low magnetic field phase for the dimerized systems, we find an exponential decrease in of the relaxation time and can compute the gap . Close to the quantum critical point our results are in good agreement with the scaling behavior based on the existence of free excitations.

## I Introduction

Quantum spin systems can exhibit a very rich set of phases, although they are usually described by very simple HamiltoniansAuerbach (1998). These range from phases with long range magnetic order to spin liquids, or even phases for which the order is more complex and rests on non-local order parameters. Understanding the behavior of such systems is thus an extremely challenging task with potential use for quantum computation or quantum simulation of other types of systemsWard et al. (2013). Among all spin systems, low dimensional ones such as spin chains and ladders are particularly interesting since quantum effects are large. In one dimension interactions between the excitations lead to various exotic states, ranging from gapped phases to phases possessing quasi-long range magnetic order known as Tomonaga-Luttinger liquidsGiamarchi (2004) (TLL).

In order to examine the various types of order that can be present, it is important to have a set of probes sensitive to correlation functions of the system. Fortunately a set of such probes such as neutron scattering, electronic spin resonance, Raman scattering and Nuclear Magnetic Resonance (NMR)Abragam (1961); Slichter (1989) exists. The NMR allows for various measurements such as the Knight shift, measuring the local magnetic field, or the so-called relaxation time, essentially sensitive to the decay of the local spin-spin correlation functions.

Although the principle of what is measured by the ratio is simple and relates to local spin-spin correlations, the theoretical determination is far from trivial. Very often various schemes of approximations of the exact formula are used. The first approximation consists in assuming that the NMR frequency is low enough (usually in the hundred of MHz range) compared to the temperature and that it can safely be set to zeroAbragam (1961); Slichter (1989). The second approximation is usually to reduce the local correlation function, which is a sum over all momenta, to a sum taken around special momentum values (e.g. the values and values around either the antiferromagnetic wavevector , where is the lattice spacing, or a corresponding incommensurate one when the magnetization is finite). This last approximation is reasonable when the characteristic scale of excitations (typically the temperature) is low enough compared to the magnetic exchange, so that the excitations around these wavevectors can be well separated. Finally, to compute the correlations, some continuous approximations such as bosonization, exploiting the above points, are usually employed. This set of approximations has allowed a connection between NMR measurements and theoretical predictions for quantum chains and ladder systems.

In the recent years a successful set of magnetic systems, in which the magnetic exchanges are considerably lowerGiamarchi et al. (2008); Zapf et al. (2014) than in previously used materials, typically around 10K, has been developed. These materials have the advantage over previously studied ones that they can be manipulated by the application of experimentally reachable magnetic fields, from zero magnetization to full saturation. This tunability opens the possibility to investigate new physics such as the universality of the TLLKlanjšek et al. (2008); Dóra et al. (2007); Sato et al. (2011), and to use these magnetic systems as quantum simulatorsWard et al. (2013) for other quantum systems such as itinerant bosonsZapf et al. (2006); Bouillot et al. (2011); Mukhopadhyay et al. (2012); Schmidiger et al. (2013); Jeong et al. (2013). Since in these materials the exchange energy scale is now much closer to the typical measurement temperatures, it invalidates partly, or pushes to very low temperatures, the above mentioned approximations. Thus, in connection with this new class of materials, a direct method to compute the NMR relaxation time without having to resort to these approximations is needed.

This is what we undertake in the present paper, by using a time dependent matrix product state method (tMPS) Vidal (2004); White and Feiguin (2004); Daley et al. (2004); Schollwöck (2011) to compute directly the relaxation time at finite temperature Feiguin and White (2005); Verstraete et al. (2004); Zwolak and Vidal (2004); Barthel (2013). An alternative approach would be the direct calculation of the correlation functions in the frequency domain as for example in Ref. Tiegel et al., 2014.

However, here we have chosen to use the time-dependent method to evaluate the time dependence of the local spin-spin correlations, since these can be directly related to the ratio . We calculate the quantity as a function of the temperature, even in regimes for which the temperature is not negligible compared to the magnetic exchange. Previous works have shown numerical results for autocorrelations and relaxation time obtained via DMRG techniquesNaef et al. (1999); Sirker (2006) or exact diagonalization methodsFabricius and McCoy (1998).

The plan of the paper is as follows. In Section II we briefly introduce the models we consider by discussing their Hamiltonians and phase diagrams. Section III defines and discusses the spin-lattice relaxation mechanism, its relation to spin-spin correlation functions, and some analytical results valid in the low temperature limit which we will use to benchmark our numerical results. In Section IV we describe the procedure adopted for the numerical computation. We then move in Section V to the results for the different models considered, namely XXZ spin chains and dimerized ones. We show how the results of the numerical calculations connect with the standard field theoretical approaches for . Section VI presents conclusions and perspectives. In the Appendices we give additional details about the computations and some preliminary tests made to check the robustness of the code.

## Ii Models

We consider spin-1/2 chains characterized by different anisotropies of the coupling, or even dimerization. In addition, a magnetic field is applied along the -direction. The first model we investigate is the antiferromagnetic XXZ chain, whose Hamiltonian is given by

 H=J∑j[12(S+jS−j+1+h.c.)+ΔSzjSzj+1]−h∑jSzj, (1)

where is a spin operator for a spin on site , denotes its direction, and the Pauli matrices. are the spin rising and lowering operators. The parameter gives the spin coupling strength, is dimensionless and measures the anisotropy, is the amplitude of the applied magnetic field along the direction. The factor, the Bohr magneton and have been absorbed into and , which both have here the dimensions of an energy.

For the isotropic case , the model corresponds to the Heisenberg (or XXX) Hamiltonian, while for we have the XX model which can be mapped via a Jordan-Wigner transformation Giamarchi (2004) onto a free-fermion model with a fixed chemical potential. The phase diagram of the XXZ model is given in Fig. 1. The boundary between the XY and ferromagnetic phases is given by . The boundary between the XY and Néel phases is given by the triplet gap, which is a function of Mikeska and Kolezhuk (2004). In this work we will limit ourselves to the case .

Additionally, we consider the dimerized Heisenberg chain which is described by the Hamiltonian

 H=∑j(J+(−1)jδJ)Sj⋅Sj+1−h∑jSzj. (2)

where denotes the vector of the spin at site . Here is the strong exchange coupling on every second bond and the weak coupling of the other bonds. A pictorial representation of this system is given in Fig. 2. At zero magnetic field () such a model has a non-trivial spin-0 ground state (spin liquid) with a gap to the first excitation which is a band of spin-1 excitations (triplons)Giamarchi (2004). This is particularly easy to see in the limit of large values of the dimerization. In this limit, strongly, antiferromagnetically coupled spin dimers are formed due to the strong exchange coupling on every second bond. These dimers are themselves coupled by the weaker interaction . The lowest excitations are the one-triplon excitations which can be accurately described as a single dimer excited from spin 0 () to spin 1 (, ), delocalized on the chain (see Fig. 3).

The role of the magnetic field along is to progressively close the gap to the excitations. At a quantum critical phase arises with gapless excitations. For stronger magnetic fields the spins of the chain become polarized, and above the second critical magnetic field the ground state is fully polarized with a gapped spectrum. For more details on dimerized chains see e.g. Ref. Giamarchi, 2004 and references therein.

## Iii Spin-lattice relaxation time T1

We introduce in this section the so called spin-lattice relaxation time , one of the important time-scales of NMR measurements. In NMR experiments the nuclear spins of the sample, previously polarized by an applied magnetic field, are perturbed using an electromagnetic pulse. The time constant characterizes the process by which the component of the nuclear magnetization along the direction of the applied magnetic field (denoted here by ) reaches thermodynamic equilibrium with its surroundings (the lattice) after the perturbationAbragam (1961); Slichter (1989). The evolution of the nuclear magnetization along is:

 Mz(t)=Mz,eq(1−e−t/T1). (3)

Quite generally in a solid the ratio can be related directly to the spin-spin correlations of the electronic system. Using the fact that the nuclear-electronic coupling (which is the hyperfine one) is very weak, one obtainsSlichter (1989) the relation:

 1T1=γ2n2[A2⊥(Sxx(ω0)+Syy(ω0))+A2∥Szz(ω0)], (4)

where is the nuclear gyromagnetic ratio of the measured nuclear spin, and are the longitudinal and transverse components of the hyperfine tensor, with are the local spin-spin correlation functions at the nuclear Larmor frequency , and

 (5)

where denotes the thermal and quantum average given by

 ⟨⋯⟩=Tr[e−βH⋯]Tr[e−βH]. (6)

Note that in formula (4) we have implicitly assumed that the hyperfine coupling term was essentially independent. This covers a large number of cases, for example the ones in which the relaxation is measured on the site carrying the electronic spin. There are also interesting cases for which the dependence of the hyperfine term can filter some modes, for example the modes at if the relaxation is measured mid-point between two neighboring sites. This leads to different formulas and interesting propertiesSirker et al. (2009, 2011); Karrasch et al. (2015). Note that techniques similar to the ones used here but computing the finite temperature, space and time dependent spin correlations allow to treat this problem as well. We leave this more complicated case for further studies, and focus here to the generic case for which the local spin-spin correlation is sufficient.

The first two terms in Eq. (4) can be conveniently expressed in terms of the and operators

 Sxx(ω0)+Syy(ω0)=12[S+−(ω0)+S−+(ω0)] (7)

The time integral over infinite time is only valid theoretically, since neither in the experiment nor in the simulation one could expect doing the sum over an infinite interval of time. In practice, two time scales compete. One is the typical time above which one can expect the oscillations coming from the frequency to become strong and regularize the integral. The second time scale hidden in the correlation itself is the decay of the correlation linked to the finite temperature. Since typical NMR frequencies are of the order while the typical lowest temperatures at which such experiments are done are of the order of , for all practical purposes we can expect that the decay due to the temperature regularizes the integral. We will thus in the following give the expression by taking this usual limit and keeping in Eq. 5 a finite integration domain up to a maximum time . We will see that this time is important not only from the numerical point of view, but also because in some cases, at high enough temperatures, the approximation of setting the frequency to zero leads to divergences.

Using the approximations discussed above, one obtains that

 Sλμ(ω0→0)≃∫+t0−t0dt⟨Sλj(t)Sμj(0)⟩=2∫+t00dtRe⟨Sλj(t)Sμj(0)⟩. (8)

Now can be or . Since we have set inside the integral in the above expression, and considered that as explained above, the two time integrals of and correlations also become identical (note of course that this is not the case for the correlations themselves at finite time). We can thus compute the one that is the most convenient numerically depending on the specific case.

Since in this work we focus on the parameter dependence of generic Hamiltonians we will omit the factors and , which depend on the specific material. For a specific material they have to be considered and in general both terms might be important. However, in this work we are not focusing on a specific material and have chosen to consider for each example only one of the terms separately. We thus compute numerically

 (1T1)±∓ =2∫+t00dtRe⟨S±j(t)S∓j(0)⟩, (9) (1T1)zz =2∫+t00dtRe[⟨Szj(t)Szj(0)−m2⟩]. (10)

Note that with the definitions in Eq. 9 and Eq. 10 the units of become time and not one over time as for the original definition in Eq. 3.

The NMR relaxation rate in one dimension can be computed in the low energy TLL representationGiamarchi (2004). This calculation is valid when the temperature is low enough compared to the typical spin energy scales. In that case, neglecting the subdominant temperature corrections and the contribution in Eq. 4 (small compared to the term), one finds Bouillot et al. (2011)

 1T1 ≃ limω0→0−2βω0ImχR+−(x=0,ω0) ≃≃ 4Axcos(π4K)u(2πkBTu)12K−1B(14K,1−12K), (11)

where and are the TLL parameters associated to the model, and is the amplitude coefficient relating the microscopic spin operator on the lattice with the operators in the continuous field theory. These coefficients have been computed both analytically and numerically in various contexts ranging from chains to laddersLukyanov and Zamolodchikov (1997); Giamarchi and Tsvelik (1999); Hikihara and Furusaki (2001, 2004); Bouillot et al. (2011). is the retarded, onsite, correlation function at the frequency (for the resolved susceptibility see Refs. Cross and Fisher, 1979; Schulz and Bourbonnais, 1983; Giamarchi, 2004). Note also that in this formula and the lattice spacing has been set to one, thus omitted.

Eq. 11 has provided a quantitative estimation of the NMR in ladder systems for which the relaxation time could be measuredKlanjšek et al. (2008); Jeong et al. (2016). It will thus provide both a benchmark for the numerical evaluation of the relaxation time, as well as an estimation of the deviation from these ideal low energy properties.

## Iv Numerical procedure

As discussed in the previous section we need to compute correlation functions of the form

 ⟨^B(t)^A(0)⟩T=Tr(^ρβ^B(t)^A). (12)

Here and are spin operators with the relation . The expectation values of the operators are taken with the finite temperature density matrix

 ^ρβ=e−βH/Zβ, (13)

where and the inverse temperature . The time evolution of the operators is represented in the Heisenberg picture with .

In order to use the matrix product state (MPS) representation at finite temperatureFeiguin and White (2005); Verstraete et al. (2004); Zwolak and Vidal (2004); Barthel (2013); Schollwöck (2011), the density matrix is encoded by a corresponding purificationNielsen  and  Chuang (2000) which is a pure state in an enlarged Hilbert space:

 ^ρβ⟶∣∣ρβ⟩∈H⊗Haux , (14)

such that

 Traux∣∣ρβ⟩⟨ρβ∣∣=^ρβ. (15)

We choose the auxiliary space and denotes the trace over this space. An MPS representation of can be obtained by applying an imaginary time evolution, starting from the maximally entangled state

 |ρ0⟩∝∑σ|σ⟩⊗|¯σ⟩aux (16)

with denoting the state not equal to . This maximally entangled state corresponds to the physical infinite temperature state , i.e. if one traces out the auxiliary degrees of freedom one obtains the identiy. Further, in each term the state is chosen such that the magnetization is conserved in the following calculations, which enlightens considerably the numerical effort which needs to be spent.

Using the cyclic property of the trace, and expliciting the time dependence of the operator and the density matrix, one can rewrite Eq. 12 as

 ⟨^B(t)^A⟩T=1ZβTr([e−βH/2]^B[e−iHt/ℏ^Ae−βH/2eiHt/ℏ]). (17)

The square brackets indicate which parts of this expression are approximated as an MPS Barthel (2013). The bracketing is not unique and several different approaches exist. However, we found this one tested in Ref. Karrasch et al., 2012 to be often the most efficient for the here considered correlations. The approximation of the bracketed operators is calculated using an imaginary and real time evolution and the application of the local operators and . In Fig. 4 the scheme is sketched. In each step of the real or imaginary time evolution, the evolved operators are approximated by an MPS with bond dimensions as small as possible for a given constraint on the truncated weight. The convergence of our results with the chosen truncated weight is assured. Typical values for maximum bond dimension used here are up to states. Depending on the magnitude of the singular values after each decomposition, we keep those which are bigger than a minimal truncation . This has been chosen of the order of for imaginary time evolution and for real time evolution.

As discussed in the previous section (see Eq. 8), we are especially interested in onsite, dynamical spin-spin correlation functions at finite temperature T:

 ⟨^B(t)^A(0)⟩T⟶⟨Sλj(t)Sμj(0)⟩T, (18)

where is the site index and or . From now on for practical reasons we will denote

 ⟨Sλj(t)Sμj(0)⟩T=SλμT(t). (19)

In order to access the desired results (Eqs. 9 and 10) the time integral of the real part of these correlations from to is required. Numerical results are taken at discrete times which are multiples of the time-step chosen for the real time evolution within tMPS. Typical values of the steps are and respectively for the imaginary and real time evolution, for the XXZ system. For the dimerized system we choose typically and . The convergence with the time step of the time evolution is assured. The amplitude of the time step for real time evolution is chosen to be small enough to guarantee good approximation of the proper integral:

 2∫+tmax0dtReSλμT(t)≈≈N∑l=1δt[ReSλμT((l−1)δt)+ReSλμj(lδt)], (20)

where is the total number of time steps at which the correlations are evaluated, is the amplitude of a single time step and . Depending on the available computational resources and on the constraints on the desired precision, runs are stopped after a certain . In many cases this is large enough such that correlations are practically zero for larger times. In other cases it is not possible due to the numerical complexity to reach such a large . In order to have an idea of the value of the extended integral, we extrapolate its value for and associate to it an error bar. In Appendix A the details of the extrapolation method and the determination of the error bars are given.

The numerical results shown in this work are obtained for the XXZ model for a chain of size and the correlations are evaluated at the central site . For the dimerized model and . The system sizes were chosen such that the perturbations do not yet reach the boundary of the system for times up to . The resulting finite system size effects are small compared to the uncertainties introduced by the finite cut off of the time-integral and are therefore neglected.

To test the accuracy of the described procedure, some calculations have been performed for the XX model and compared with exact analytical results in Appendix B.

## V Results

In the following subsections we present our numerical results for the Heisenberg model, the XXZ model and the dimerized model.

### v.1 Heisenberg model

Let us start by considering the Heisenberg model, i.e. the XXZ model in Eq. 1, with isotropic coupling , of spins 1/2 under a magnetic field, applied along the direction. A pictorial representation of the phase diagram as a function of is given in Fig. 5. At low magnetic field the ground state of the system is a gapless TLL, whereas above the critical magnetic field () a gapped phase develops. In order to explore the different phases and the quantum critical point we focus in the following on the fields and in the gapless phase, at the quantum critical point and in the gapped phase.

In Fig. 6 the results for the relaxation rate of the correlations are shown at the magnetic field values , in the gapless phase, and at the quantum critical point. For the numerical calculations on the Heisenberg model and the XXZ model we have chosen a chain of size , a minimal truncation of , a retained states maximum of 400, and steps for the imaginary time evolution. For the real time evolution we have chosen the minimal truncation , a retained states maximum of 800, a maximal truncated weight of , a time step and .

The behavior of the relaxation time at zero magnetic field has been investigated previously both by analytical Schulz (1986); Sachdev (1994) and numerical methods such as the quantum Monte-Carlo using the maximum entropy method in order to continue to the real time axisSandvik (1995); Starykh et al. (1997) and tDMRG methods Sirker (2006). In the asymptotic low-temperature limit one obtainsCross and Fisher (1979); Schulz and Bourbonnais (1983) if logarithmic corrections are neglected. Since the behavior at has been studied in detail in Sandvik, 1995; Sirker, 2006, we here show the case for mainly for comparison. At temperatures above an almost linear increase of the relaxation time with increasing temperature can be seen. At low temperature shows an almost constant behavior. At temperatures below it even increases again while lowering the temperature. This behavior is consistent with the logarithmic corrections and has been analyzed in Sandvik, 1995. The rise at larger temperatures has been treated in Sirker, 2006 and has been found compatible with an exponential increase with a scale of the order of the magnetic exchange .

In the TLL region (), as described in Sec. III, the low temperature behavior of the relaxation rate should be approximately described by an algebraic decay with the exponent , which is fully determined by the TLL parameter Giamarchi and Tsvelik (1999); Klanjšek et al. (2008). At larger temperature a breakdown of this low energy prediction is expected. In Fig. 6 our numerical results for are quantitatively compared to the TLL predictions. Up to temperatures of about the numerical points agree within the error bars with the prediction made in Eq. 11. This comparison is achieved using previously extracted values for the Tomonaga-Luttinger parameter , the amplitude from Refs. Hikihara and Furusaki, 2001, 2004, and (lattice spacing equal to 1) extracted from separate calculations which we performed using standard finite-size DMRG methods. More details about the determination are given in Appendix C. Thus, all parameters in Eq. 11 are fully determined. For temperatures larger than the numerical results are much higher than the decaying TLL prediction. This is to be expected and clearer in the space for which the local correlation can be seen as a sum over all points. The TLL formula only contains the part coming from one of the low energy points ( in the absence of magnetic field). At higher temperatures other points start to contribute significantly to the sum. The numerical results even seem to show a slight maximum around and then remain more or less constant in value up to the shown maximal temperature.

At the quantum critical point , an algebraic divergence of the with the temperature is also expected in the low T limit. It is predicted to behave as as obtained in Refs. Sachdev et al., 1994; Chitra and Giamarchi, 1997; Orignac et al., 2007. This behavior has been experimentally observed for example in the Heisenberg chain compound copper pyrazine in Ref. Kühne et al., 2003, and discussed in Ref. Jeong and Rønnow, 2015, In this situation the prefactor is not easily extracted and therefore we fit the expected algebraic behavior with a free fit parameter . We obtain very good agreement of our numerical results with the fit using the value in the entire regime of temperatures up to , as shown in Fig. 6. This means that our results are in agreement with the predictions of the quantum critical regime extending up to these temperatures. In addition in Fig. 7 we offer a comparison between our magnetization data computed at finite temperature, and the scaling function close to (or on) a field-induced quantum critical point, which have been calculated and used in the literatureAffleck (1991); Sachdev et al. (1994); Jeong and Rønnow (2015). At , from Eqs. (2) and (3) in Jeong and Rønnow, 2015, we know that the magnetization per site should behave as

 m(T) =mS−(2kBTJ)d/2M(δhc/kBT)= (21) =0.5−0.24312√kBTJ (22)

where is the magnetization per site at saturation, is the dimension of the system, is the distance from the critical field and

 M(δhc/kBT)=1π∫∞0dx 1ex2−(δhc/kBT)+1. (23)

We observe a good agreement between numerical results and the analytical prediction, which is supposed to be valid in the low temperature limit. We see MPS data approaching the analytics as the temperature is lowered.

In Fig. 8 the results for the relaxation rate for the correlations at a field of are reported. The system at this magnetic field exhibits a gapped energy spectrum, and we denote the gap by . Due to the gapped energy spectrum an exponential decay with temperature is expected and indeed observed numerically. In order to validate the exponential form for a different parameter regime we consider also the XX model at the same magnetic field which can be mapped onto free fermions. The corresponding gap in the energy spectrum is given by . We computed the longitudinal part of the (the density-density correlation for the corresponding free fermions) analytically. Error bars come from the extrapolation, since correlations were evaluated up to a finite time. Also in this case a fit with an exponential function of is perfectly compatible with our calculations and validates our procedure.

### v.2 XXZ model

In order to explore more in detail the behavior of the relaxation time in the TLL phase, we move to the spin-1/2 XXZ model, Eq. 1, in absence of magnetic field (). For anisotropies the ground state of this model is a TLL phase. As we discussed for the Heisenberg model at , we expect that at low temperature the relaxation rate corresponding to the correlations shows an algebraic divergence as given in Eq. 11. We consider the cases as shown in Fig. 9. The corresponding values of the Luttinger liquid parameters and are calculated using the Bethe ansatz formulas given e.g. in Ref. Giamarchi, 2004, while the amplitudes are taken from Ref. Hikihara and Furusaki, 2001. Their rounded values are summarized in table 1.

The agreement between the TLL prediction of the algebraic divergences and our numerical results is extremely good at low temperatures. For larger values of the anisotropy the divergence becomes weaker until for one leaves the TLL region and no algebraic divergence is seen. As expected the predictions for the Luttinger liquid behavior disagree above a certain temperature of the order of . Above this scale our numerical results show an upturn and the different curves even cross. Our results clearly show the importance for systems with small exchange constants to be able to go beyond the asymptotic expressions in order to make comparisons with the experiments.

### v.3 Dimerized model

As a final example, we consider the isotropically dimerized spin-1/2 chain in presence of a magnetic field along the direction as defined in Eq. 2. This model can describe very well some interesting compounds like for example the copper nitrate , discussed in Refs. Tennant et al., 2003; Willenberg et al., 2015. For this compound the coupling parameters are determined as K and K and we will focus on these strongly dimerized parameters in the following. The ground state of the system at has zero magnetization. A gap of separates the ground state from the first excited state. In a magnetic field, the system shows a first quantum critical point at a magnetic field . At this point the system undergoes a transition from a gapped phase to a gapless, TLL phase. Here we focus on two cases: the gapped phase for and the TLL phase at .

In our numerical calculations we consider a chain of spins and in the imaginary time evolution a minimal truncation of , a retained states maximum of 500, and a step of . For the real time evolution we choose a minimal truncation of , a time step amplitude and a retained states maximum of 500 for temperatures , 800 for and 2000 for higher temperatures. The maximal truncated weight is in most cases, for the highest temperatures. The reached still decreases from for the lowest temperatures, to for the highest ones, according to the requested precision. We calculate the relaxation time for the onsite correlation (at , ).

Due to the presence of a gap in the absence of a magnetic field, the temperature dependence of the relaxation rate at low temperatures is expected to be exponentially activated. i.e. . In Fig. 10 we show that our results agree very well with this exponential activation.

At larger temperatures a saturation effect seems to set in. In the inset, lower temperature points have been cut because of the difficulties in the extrapolation which led to negative (though pretty close to 0) extrapolated values, as shown in the main panel of Fig. 10. A detailed description of the method used for the extrapolation and the association of an appropriate error bar is given in Appendix A.

In contrast for the case the low energy physics can be described by the TLL theory. Thus, the expected behavior of the relaxation rate as a function of the temperature is an algebraic divergence at low T of the form . The TLL parameter which enters in this formula has been determined by separate calculations of the compressibility using MPS, and the flux dependence of the energy using infinite-size MPS calculations, giving . More details about this method can be found in Appendix C. The numerically obtained relaxation time is shown in Fig. 11 and compared to the TLL predictions.

The black line represents the fit using the separately determined exponent . A constant offset has been added since the behavior is not entirely dominated by the divergence. Deviations between the comparison of the analytical prediction and the numerical calculation are seen. We attribute these deviations to the proximity of the quantum critical point. In this regime, the TLL behavior is valid only for very low temperatures . From our numerical results only the lowest temperature point lies within this region. To verify the influence of the quantum critical point, a fit using the critical power law shifted by a constant offset is performed. This fit leads for the intermediate temperature points to good results (see the green curve). A fit in which also the exponent is a fit parameter leads to an even larger exponent of (red curve).

## Vi Conclusions and outlook

Exploiting recent developments in finite temperature MPS techniques we computed the spin-lattice relaxation rate for a wide range of temperatures, for different Hamiltonians and for different quantum phases. In particular, we have considered the XX, Heisenberg, and XXZ Hamiltonians, plus the isotropically dimerized case. For the non-dimerized cases we have performed a detailed study of the gapless phase, the gapped phase and also of the quantum critical point. Numerical results were in very good agreement with analytical results available in the low-temperature limit. We have shown the deviation from the low-T law at finite temperature and we have swept through quantum critical points, situations in which theoretical results are more difficult to obtain. Our calculations prove that the MPS method can be successfully used to obtain the NMR relaxation time in regimes in which the field theoretical asymptotic values would not be applicable. The overlap between the regimes in which the numerical methods are applicable and the regime covered by the field theoretical asymptotic methods allows essentially a full description of the NMR behavior for the accessible regime of temperatures.

Having a method which can quantitatively compute the NMR relaxation time from a given microscopic Hamiltonian rather than simple asymptotic expressions should allow to test that the microscopic Hamiltonian does not miss an important term, and to fix the various coefficients by comparing the computed temperature dependence with the experimentally measured one. This is similar in spirit to what was achieved by the comparison of the computed neutron spectra with the measured ones for DIMPYSchmidiger et al. (2012). Another interesting direction is the investigation of the behavior of the relaxation mechanism of the spin excitations close to the quantum critical point. Indeed the nature of the relaxation mechanism is potentially different depending on whether one considers the term or the ones. For 3D systems a self energy analysis of the transverse part of was suggesting Orignac et al. (2007) a behavior due to the necessity of making three magnon excitations to be able to scatter a magnon and get a finite lifetime while the part leads, as shown in the present paper, to . Our numerical results which are able to correctly determine the exponential decay in the controlled cases of the longitudinal excitations are thus potentially able to address this issue and potentially make contact on the experiments on that point Mukhopadhyay et al. (2012). Such a study clearly going beyond the scope of the present paper, is thus left for future works.

The present method works efficiently if the systems are one or quasi-one dimensional. One important challenge on the theoretical level is to extend the present analysis to the case of higher dimensional systems. In that case, although other methods such as quantum Monte-Carlo exist, the dynamical correlations in real time are still a challenge for which the MPS methods could bring useful contributions. Indeed the (numerically) rather complete knowledge of the one-dimensional correlation functions allow to incorporate them into approximation schemes such as RPA to capture a large part of the higher dimensional physics. Another route is to solve clusters of one dimensional structures, which allows to at least incorporate part of the transverse fluctuations.

Note added: Just after submitting this work, a related numerical study by M. Dupont, S. Capponi and N. Laflorencie appeared Dupont et al. (2016). Our results are perfectly compatible with each other when comparison can be made.

###### Acknowledgements.
We acknowledge fruitful discussions with P. Bouillot and on NMR with C. Berthier and M. Horvatić. This work was supported in part by the Swiss NSF under Division II and the DFG and the ERC (Grant Number 648166, Phon(t)on).

## Appendix A Extrapolation method

As discussed in Sec. IV, the numerical results for correlations are only available up to a certain time . Since in principle the time integral of these correlations should be performed up to , one needs to find a way to approximate the value of the extended integral and of the associated error bar. In order to do this we study the behavior of the integral as a function of . We perform a linear fit of the value of the integral as a function of at the largest available values of . We use this fit to extrapolate the value of the integral to . If the value of the integral still shows a considerable trend, we associate to the extrapolated value a one-sided error bar corresponding to the difference between the extrapolated value and the value of the integral for the maximum available. An example is shown in Fig. 12 for the dimerized chain in the TLL phase. For the case where the integral oscillates around a certain value and no clear trend is visible, we choose to associate a symmetric error bar with semi-amplitude equal to the distance between the extrapolated value itself and the value of the integral for the maximum available. An example is shown in Fig. 13 for the dimerized chain in the gapped phase. In both cases, the extracted error bars should give a (most probably pessimistic) estimate of the uncertainty on the value of the integral.

## Appendix B Consistency test using the XX model

To test the accuracy of the numerical procedure, we performed calculations for the XX model under a magnetic field along the direction:

 H=J2∑j(S+jS−j+1+h.c.)−h∑jSzj. (24)

For this specific model we focus on two specific cases, (gapless phase) and (gapped phase), and on correlations. In particular, we determine for different temperatures the ratio for correlations, which we define here as:

 (1T1)zz=tmax∫0dtReSzzj,T(t). (25)

We compare our numerical results obtained via the described procedure using MPS, with exact analytically results Giamarchi (2004). In the limit of an infinite-size system, the exact result for the onsite correlations at a temperature and is given by

 Szzj,T(t)=J0(Jt/ℏ)2π⋅+π∫−πdk eiλkt/ℏ⋅fk(β) −−14π2⋅∣∣ ∣∣ +π∫−πdk eiλkt/ℏ⋅fk(β) ∣∣ ∣∣2, (26)

where is the 0th-order Bessel function of the first kind, is the imaginary unit, , where is the dimensionless momentum and

 fk(β)=11+eβλk (27)

is the Fermi function, where is the inverse temperature.

For one obtains

 Szzj,T(t)=J0(Jt/ℏ)2π⋅eiht/ℏ⋅+π∫−πdk eiλ′kt/ℏfk(β) −−14π2⋅∣∣ ∣∣ +π∫−πdk eiλ′kt/ℏfk(β) ∣∣ ∣∣2+ 14 ++14π2⋅⎡⎢⎣+π∫−πdk fk(β)−2π⎤⎥⎦⋅+π∫−πdk fk(β). (28)

Here and .

As for the numerical procedure, simulations are performed for a chain of spins. Onsite correlations are measured in the center of the chain () to avoid boundary effects. Imaginary time evolutions are performed using the following parameter set: minimal truncation , retained states maximum 400 and step . Real time evolutions are performed using: minimal truncation , a retained states maximum of 800, maximal truncated weight , and time step up to . Results of the comparison theory-numerics are reported in Fig. 14. The agreement between the analytical and the numerical results is extremely good at all temperatures, which justifies our procedure.

## Appendix C Determination of the TLL parameters

To get the values of the TLL parameters and we determine first their ratio , which is related to the static TLL susceptibility, and their product , related to the variation of the energy with a flux. Then, the two values of and trivially follow by recombination of the two previous results.

The ratio is determined from the static TLL susceptibility of the system according to the relationGiamarchi (2004)

 Ku=πLd2E0dM2, (29)

where is the size of the system, is the ground state energy, and is the total magnetization. The second derivative has to be discretized since in a spin-1/2 system can only vary by integer steps (thus ):

 Ku(M)=πL[E0(M+1)+E0(M−1)−2E0(M)]. (30)

can be evaluated at fixed values of magnetization via standard finite-size DMRG. The magnetization, defined at the beginning of the simulation by the initial distribution of spins, is set as a conserved quantum number.

The product can be determined by studying the variation of the ground state energy of the system in response to a variation of a flux through the system. To be more precise, for a fixed value M of the magnetizationGiamarchi (2004)

 uK(M)=πLd2E0(Φ,M)dΦ2∣∣∣Φ=0, (31)

The flux is represented by twisted periodic boundary conditions, . This condition can be transferred into the Hamiltonian via the following transformation:

 S+jS−j+1 ⟶S+jS−j+1⋅eiΦL, S−jS+j+1 ⟶S−jS+j+1⋅e−iΦL, (32)

which distributes homogeneously the total flux along the chain. For each fixed value of (and therefore of ), we evaluate within an infinite-size MPS algorithm for symmetric values of around zero. The resulting ground state energy close to can be approximated by a parabola and we fit the points with a second degree polynomial of the form , where and are the fit parameters. According to Eq. 31, the fitting parameter is related to the product

 uK(M)=2πLaM. (33)

## References

• Auerbach (1998) A. Auerbach, Interacting Electrons and Quantum Magnetism (Springer, Berlin, Germany, 1998).
• Ward et al. (2013) S. Ward, P. Bouillot, H. Ryll, K. Kiefer, K. W. Kramer, C. Ruegg, C. Kollath, and T. Giamarchi, J. Phys. Condens. Matter 25, 014004 (2013).
• Giamarchi (2004) T. Giamarchi, Quantum Physics in One Dimension, (Oxford University Press, Oxford, UK, 2004).
• Abragam (1961) A. Abragam, Principles of Nuclear Magnetism (Clarendon Press, Oxford, UK, 1961).
• Slichter (1989) C. Slichter, Principles of Magnetic Resonance, 3rd ed. (Springer, Berlin, Germany, 1989).
• Giamarchi et al. (2008) T. Giamarchi, C. Rüegg, and O. Tchernyshyov, Nat. Phys. 4, 198 (2008).
• Zapf et al. (2014) V. S. Zapf, M. Jaime, and C. D. Batista, Rev. Mod. Phys. 86, 563 (2014).
• Klanjšek et al. (2008) M. Klanjšek, H. Mayaffre, C. Berthier, M. Horvatić, B. Chiari, O. Piovesana, P. Bouillot, C. Kollath, E. Orignac, R. Citro, and T. Giamarchi, Phys. Rev. Lett. 101, 137207 (2008).
• Sato et al. (2011) M. Sato, T. Hikihara, and T. Momoi, Phys. Rev. B 83, 064405 (2011).
• Dóra et al. (2007) B. Dóra, M. Gulácsi, F. Simon, and H. Kuzmany, Phys. Rev. Lett. 99, 166402 (2007).
• Zapf et al. (2006) V. S. Zapf, D. Zocco, B. R. Hansen, M. Jaime, N. Harrison, C. D. Batista, M. Kenzelmann, C. Niedermayer, A. Lacerda, and A. Paduan-Filho, Phys. Rev. Lett. 96, 077204 (2006).
• Bouillot et al. (2011) P. Bouillot, K. Corinna, A. M. Läuchli, M. Zvonarev, B. Thielemann, C. Rüegg, E. Orignac, R. Citro, M. Horvatić, C. Berthier, M. Klanjšek, and T. Giamarchi, Phys. Rev. B 83, 054407 (2011).
• Mukhopadhyay et al. (2012) S. Mukhopadhyay, M. Klanjšek, M. S. Grbić, R. Blinder, H. Mayaffre, C. Berthier, M. Horvatić, M. A. Continentino, A. Paduan-Filho, B. Chiari, and O. Piovesana, Phys. Rev. Lett. 109, 177206 (2012).
• Schmidiger et al. (2013) D. Schmidiger, P. Bouillot, T. Guidi, R. Bewley, C. Kollath, T. Giamarchi, and A. Zheludev, Phys. Rev. Lett. 111, 107202 (2013).
• Jeong et al. (2013) M. Jeong, H. Mayaffre, C. Berthier, D. Schmidiger, A. Zheludev, and M. Horvatić, Phys. Rev. Lett. 111, 106404 (2013).
• Vidal (2004) G. Vidal, Phys. Rev. Lett. 93, 040502 (2004).
• White and Feiguin (2004) S. R. White and A. E. Feiguin, Phys. Rev. Lett. 93, 076401 (2004).
• Daley et al. (2004) A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, J. Stat. Mech.: Theor. Exp. P04005 (2004).
• Schollwöck (2011) U. Schollwöck, Ann. Phys. 326, 96 (2011).
• Feiguin and White (2005) A. E. Feiguin and S. R. White, Phys. Rev. B 72, 220401 (2005).
• Verstraete et al. (2004) F. Verstraete, J. J. García-Ripoll, and J. I. Cirac, Phys. Rev. Lett. 93, 207204 (2004).
• Zwolak and Vidal (2004) M. Zwolak and G. Vidal, Phys. Rev. Lett. 93, 207205 (2004).
• Barthel (2013) T. Barthel, New J. Phys. 15, 073010 (2013).
• Tiegel et al. (2014) A. C. Tiegel, S. R. Manmana, T. Pruschke, and A. Honecker, Phys. Rev. B 90, 060406(R) (2014).
• Naef et al. (1999) F. Naef, X. Wang, X. Zotos, and W. von der Linden, Phys. Rev. B 60, 359 (1999).
• Sirker (2006) J. Sirker, Phys. Rev. B 73, 224424 (2006).
• Fabricius and McCoy (1998) K. Fabricius and B. M. McCoy, Phys. Rev. B 57, 8340 (1998).
• Mikeska and Kolezhuk (2004) H. J. Mikeska and A. Kolezhuk, Quantum magnetism, (Springer, Lecture notes in Physics, 2004).
• Sirker et al. (2009) J. Sirker, R. G. Pereira, and I. Affleck, Phys. Rev. Lett. 103, 216602 (2009).
• Sirker et al. (2011) J. Sirker, R. G. Pereira, and I. Affleck, Phys. Rev. B 83, 035115 (2011).
• Karrasch et al. (2015) C. Karrasch, R. G. Pereira, and J. Sirker, New. J. Phys. 17, 103003 (2015).
• Lukyanov and Zamolodchikov (1997) S. Lukyanov and A. B. Zamolodchikov, Nucl. Phys. B 493, 571 (1997).
• Giamarchi and Tsvelik (1999) T. Giamarchi and A. M. Tsvelik, Phys. Rev. B 59, 11398 (1999).
• Hikihara and Furusaki (2001) T. Hikihara and A. Furusaki, Phys. Rev. B 63, 134438 (2001).
• Hikihara and Furusaki (2004) T. Hikihara and A. Furusaki, Phys. Rev. B 69, 064427 (2004).
• Cross and Fisher (1979) M. C. Cross and D. S. Fisher, Phys. Rev. B 19, 402 (1979).
• Schulz and Bourbonnais (1983) H. J. Schulz and C. Bourbonnais, Phys. Rev. B 27, 5856(R) (1983).
• Jeong et al. (2016) M. Jeong, D. Schmidiger, H. Mayaffre, M. Klanjšek, C. Berthier, W. Knafo, G. Ballon, B. Vignolle, S. Krämer, A. Zheludev, and M. Horvatić,  arXiv:1604.05252  (2016).
• Nielsen  and  Chuang (2000) M. A. Nielsen and I. L. Chuang,  Quantum Computation and Quantum Information (Cambridge University Press, Cambridge, UK, 2000).
• Karrasch et al. (2012) C. Karrasch, J. H. Bardarson, and J. E. Moore, Phys. Rev. Lett. 108, 227206 (2012).
• Schulz (1986) H. J. Schulz, Phys. Rev. B 34, 6372 (1986).
• Sachdev (1994) S. Sachdev, Phys. Rev. B 50, 13006 (1994).
• Sandvik (1995) A. W. Sandvik, Phys. Rev. B 52, 9831 (1995).
• Starykh et al. (1997) O. A. Starykh, R. R. P. Singh, and A. W. Sandvik, Phys. Rev. Lett. 78, 539 (1997).
• Sachdev et al. (1994) S. Sachdev, T. Senthil, and R. Shankar, Phys. Rev. B 50, 258 (1994).
• Chitra and Giamarchi (1997) R. Chitra and T. Giamarchi, Phys. Rev. B 55, 5816 (1997).
• Orignac et al. (2007) E. Orignac, R. Citro, and T. Giamarchi, Phys. Rev. B 75, 140403 (2007).
• Kühne et al. (2003) H. Kühne, H.-H. Klauss, S. Grossjohann, W. Brenig, F. J. Litterst, A. P. Reyes, P. L. Kuhns, M. M. Turnbull, and C. P. Landee, Phys. Rev. B 80, 045110 (2009).
• Jeong and Rønnow (2015) M. Jeong and H. M. Rønnow, Phys. Rev. B 92, 180409(R) (2015).
• Affleck (1991) I. Affleck, Phys. Rev. B 43, 3215 (1991).
• Tennant et al. (2003) D. A. Tennant, C. Broholm, D. H. Reich, S. E. Nagler, G. E. Granroth, T. Barnes, K. Damle, G. Xu, Y. Chen, and B. C. Sales, Phys. Rev. B 67, 054414 (2003).
• Willenberg et al. (2015) B. Willenberg, H. Ryll, K. Kiefer, D. A. Tennant, F. Groitl, K. Rolfs, P. Manuel, D. Khalyavin, K. C. Rule, A. U. B. Wolter, and S. Süllow, Phys. Rev. B 91, 060407(R) (2015).
• Schmidiger et al. (2012) D. Schmidiger, P. Bouillot, S. Mühlbauer, S. Gvasaliya, C. Kollath, T. Giamarchi, and A. Zheludev, Phys. Rev. Lett. 108, 167201 (2012).
• Dupont et al. (2016) M. Dupont, S. Capponi, and N. Laflorencie, Phys. Rev. B  94, 144409 (2016).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters