Low temperature dynamics of nonlinear Luttinger liquids

We generalize nonlinear Luttinger liquid theory to describe the dynamics of one-dimensional quantum critical systems at low temperatures. Analyzing density-matrix renormalization group results for the spin autocorrelation function in the XXZ chain we provide, in particular, direct evidence for spin diffusion in sharp contrast to the exponential decay in time predicted by conventional Luttinger liquid theory. Furthermore, we discuss how the frequencies and exponents of the oscillatory contributions from the band edges are renormalized by irrelevant interactions and obtain excellent agreement between our finite temperature nonlinear Luttinger liquid theory and the numerical data.

Introduction-The dynamics of quantum critical systems at finite temperatures is an outstanding challenge in many-body physics [1]. Understanding the combined effects of thermal fluctuations and interactions is crucial for the interpretation of experiments that probe frequencydependent responses and scattering cross sections. Addressing this problem has become even more pressing since fascinating new experiments in cold atomic gases [2][3][4][5] and condensed matter [6] have given us access to correlation functions directly in the time domain.
Much of the interest in real time evolution of manybody states has focused on one-dimensional (1D) models. In particular, the difference in the nonequilibrium dynamics of integrable versus nonintegrable systems is currently a hotly debated topic [7][8][9][10][11]. In parallel, studies of ground state dynamical correlations have recently forced us to revise our understanding of quasiparticles in critical 1D systems [12][13][14][15][16][17][18][19], culminating with the development of the nonlinear Luttinger liquid (NLL) theory [20]. This theory provides a framework for calculating exponents of edge singularities in dynamical response functions based on a picture of fermionic quasiparticles with nonlinear dispersion, reminiscent of elementary excitations in Bethe ansatz (BA) solvable models [21,22]. NLL theory predicts that the long-time decay of correlation functions is dominated by excitations involving particles or holes near band edges, explaining the high-frequency oscillations observed numerically [21,23] and confirmed by an exact form factor approach [24].
The purpose of this Letter is to extend NLL theory to describe the low temperature, long time decay of correlation functions of 1D quantum fluids. Determining the precise effects of temperature is certainly relevant for experiments aimed at testing the phenomenology of NLLs. Another motivation for this work is to provide analytic expressions to examine numerical results for time dependent correlation functions, accessible by finitetemperature versions [25][26][27][28][29] of time-dependent [30][31][32][33][34] density matrix renormalization group (tDMRG) methods [35,36], and for the thermal broadening of edge singularities in the frequency domain [37]. In the following we will generalize NLL theory to T > 0 and compare the predictions with state-of-the-art tDMRG results. We will concentrate on the spin autocorrelation function G(t) of the XXZ model at zero magnetic field, but our approach can easily be generalized to other correlations and 1D models. We stress that the time decay of correlation functions in the XXZ model at finite T is still an open problem despite the integrability of the model [38,39]. Our main results are: (i) G(t) at intermediate times t is well described by a generalization of NLL theory that takes into account the effects of irrelevant operators on the dispersion of high-energy quasiparticles and on their coupling to low energy modes; (ii) at long times, G(t) contains a non-oscillating, ∼ 1/ √ t decaying diffusive term, in sharp contrast to the exponential decay predicted by Luttinger liquid theory. Diffusive behavior in spins chains has been invoked to explain the spin-lattice relaxation rate [40] and muon-spin relaxation [41] measured in quasi-1D antiferromagnets and is attributed to inelastic umklapp scattering [42]. However, fingerprints of spin diffusion have so far only been found indirectly in the current-current correlation function [27,42,43] and in the imaginary time dependence of the dynamical susceptibility [44]. Here we provide the first direct numerical evidence of spin diffusion in the autocorrelation function.
Model-Consider the 1D XXZ Hamiltonian where S a j are spin-1/2 operators, ∆ is the anisotropy parameter, and periodic boundary conditions are assumed. This model can be realized, for instance, in the Mott insulating phase of two-component Bose mixtures in 1D optical lattices [2,45]. The spin autocorrelation function at temperature T is defined by where Z = Tr e −H/T is the partition function. The XXZ model is equivalent to spinless fermions, c ( †) j , via a Jordan-Wigner transformation with S z j → c † j c j − 1/2, so that ∆ plays the role of a nearest-neighbor densitydensity interaction [46].
Noninteracting case-For ∆ = 0, the XXZ model reduces to a free fermion model. The autocorrelation factorizes into a product of free particle and hole Green's functions and is given exactly by [47] with dispersion k = − cos k and Fermi distribution f k = 1/(e k /T + 1) at half-filling.
Let us discuss an approximation that captures the asymptotic long time decay of G (0) (t). First note that G (0) (t) oscillates at arbitrary temperatures due to saddle points of the integrand where d k /dk = 0, which occur at k = 0 and k = π. In the low-temperature regime T 1, we employ a mode expansion of the fermion field which keeps only states within sub-bands near the smeared Fermi surface and near the saddle points, in the form Here ψ R,L are the low-energy right-and left-moving components, whiled and d are high-energy modes:d † creates a hole at the bottom of the band (k = 0), and d annihilates a particle at the top of the band (k = π).
The mode expansion reduces the problem to the calculation of free propagators for ψ R,L and d,d. The low energy modes have linear dispersion, thus their propagator at T 1 is given by the standard conformal field theory result ψ R,L (x, t)ψ † R,L (x, 0) ∼ πT / sinh(πT t). On the other hand, the high-energy modes do not feel the temperature, up to an exponentially small correction in the Fermi distribution. Since the dispersion is parabolic near the band edges, k ≈ 1 − k 2 /2, the propagator of d,d de- Substituting the contributions from low and high energy modes into Eq.
The first term stems from the contribution in which both particle and hole are high-energy modes. The second term is due to a high-energy particle (hole) plus a low-energy hole (particle) at either one of the Fermi points. The latter decays exponentially for t 1/T . The third contribution is due only to the low-energy modes ψ R,L , and decays more rapidly than the oscillating terms. Therefore, G (0) (t) for t 1/T is dominated by the term oscillating with frequency ω = 2 and decaying as 1/t.
Interacting case: NLL theory-We focus on the regime 0 < ∆ < 1 corresponding to a critical phase of fermions with repulsive interactions. The factorization of G(t) as in Eq. (3) is then a priori lost. However, we can investigate the long time decay of G(t) using the framework of NLL theory [20]. The idea is to keep the mode expansion with both low and high energy modes, which are now identified as quasiparticles in a renormalized band. The ground state is a vacuum of d andd, and the spin operator S z j ∼ c † j c j creates at most one d particle and/or oned hole. At low temperatures, we neglect an exponentially small thermal population of the band edge modes. The d,d quasiparticles can then be treated as mobile "impurities" distinguishable from the low energy modes.
At T = 0 the dynamics is described by an effective field theory with Hamiltonian density [21] Here ε is the energy and m the absolute value of the effective mass of the band edge modes, V is the impurityimpurity interaction, φ(x) and θ(x) are dual fields representing the bosonized low-energy modes [46] and obey . v is the spin velocity, K the Luttinger parameter, and α the dimensionless coupling constant of the impurity-boson interaction. For the integrable XXZ model one finds Eq. (5) must be regarded as a fixed point Hamiltonian which includes all marginal interactions allowed by symmetry. This model can be solved exactly by performing a unitary transformation that decouples the impurities from the bosonic modes, but attaches "string" operators to the d,d fields, [14]. This allows one to predict the long time decay of G(t) at T = 0, up to non-universal amplitudes [21].
We now extend NLL theory to 0 < T ε (low temperatures compared to the renormalized bandwidth), using Eq. (5) as a starting point. We obtain the leading T dependence by analyzing the effects of irrelevant interactions. The leading corrections to the oscillating terms in G(t) allowed by symmetry [48] stem from the irrelevant dimension-three operators Substituting the mode expansion into Hamiltonian (1) and bosonizing the low-energy modes, we find for ∆ 1: g ≈ −∆ , µ − ≈ −∆/ √ π, while g and µ + are not generated to first order in ∆. The g interaction is particularly important: it can be identified with a three-body scattering process [16,21] and gives rise to a nonzero impurity decay rate for T > 0 [50,51]. However, by imposing nontrivial conservation laws in the XXZ model we can show that g = 0 exactly [48,49], as expected from the lack of three-body scattering in integrable models.
We calculate the impurity self-energy Σ and effective impurity-boson interactionα by perturbation theory in the irrelevant operators. The leading T dependence is determined by loop diagrams which contain only one irrelevant coupling constant but arbitrary factors of the bare α. In the calculation of loop diagrams, it is convenient to treat the quadratic term in the impurity dispersion (which is also a dimension-three operator) as a perturbation, expanding the internal impurity propagators in powers of 1/m. We find The prefactors c Σ and c α are linear functions of g, µ + , µ − .
, we obtain the weak coupling approximation Omitted terms are O(∆ 3 ) or higher. To first order in ∆, the self-energy and vertex correction are both governed by g ≈ −∆. For 0 < ∆ 1, Eq. (7) thus implies that the effective impurity energyε(T ) = ε + Σ(T ) decreases ∼ T 2 , whereas the effective couplingα increases linearly with T . Phenomenologically, we find g = πv 2 2K ∂ 2 ε ∂h 2 h=0 relating the coupling constant g in Eq. (6) to a change in the impurity energy when applying a magnetic field h. This allows one to obtain g(∆) exactly using the BA and Wiener-Hopf techniques. Unfortunately, the corrections of higher order in α in Eq. (7) quickly become of the same order as the O(g) term making it impossible to fix c Σ , c α beyond the lowest order in ∆. Importantly however, the T dependence does hold for arbitrary 0 < ∆ < 1, as long as the temperature is small enough.
The renormalization of the impurity-impurity interaction appears at order g 2 . The basic process involves a two-boson loop connecting two impurity lines. The correction depends on the momentum and frequency exchange between the impurities:Ṽ (q, ω, T ) ≈ V + 2πg 2 T 2 ω 3v 4 q , where we simplified the result in the physically relevant regime ω vq for impurities with parabolic dispersion. We can now describe the decay of G(t) using the methods of NLL theory with renormalized parameters at T > 0. First, consider the contribution from the excitation with a single impurity, equivalent to the second term in Eq. (4). The unitary transformation introduces a "string" operator whose scaling dimension depends on temperature throughα(T ). The result is is the unknown prefactor. Note thatα(T ) > α for 0 < ∆ 1 implies thatη(T ) decreases with temperature, slightly slowing down the decay of G 1 (t).
The two-impurity contribution to G(t), analogous to the first term in Eq. (4), is strongly modified by the V interaction. At T = 0, the 1/t decay for ∆ = 0 changes to 1/t 2 for ∆ > 0 and t 1/V 2 [21]. This asymptotic behavior is associated with two impurities scattering in a ladder series with small energy and momentum transfer, ω ∼ q 2 /m ε. Neglecting the renormalization of V for ω vq, we obtain the two-impurity contribution where B(T ) is the unknown prefactor. For generic models, Eq. (10) must be modified to include the exponential decay due to relaxation by the three-body process g .
Diffusive decay-At temperatures T 1, conventional Luttinger liquid theory predicts that the nonoscillating terms in G(t), associated only with the lowenergy modes ψ R,L , are given in the interacting case by A [πT / sinh(πT t)] 2 + B [πT / sinh(πT t)] 2K , where the amplitudes A , B are known [52]. However, we must also consider how irrelevant operators affect the low-energy contributions. In [42] it was shown that the formally irrelevant umklapp scattering qualitatively changes the long time decay of the lowenergy, long-wavelength contribution to G(t). There appears a new time scale set by the decay rate γ ∝ λ 2 T 8K−3 . For the XXZ model, γ can be calculated exactly [42,52]. At low T and long times t 1/γ 1/T , we find a diffusive contribution [58] Strikingly, this diffusive (∼ 1/ √ t) decay in the autocorrelation function coexists with ballistic transport [42,53]. Eqs. (9,12) are the two contributions to G(t) which are dominant at intermediate and long times, respectively.
Numerical results-In order to test our theory, we now turn to a comparison with tDMRG results for finite system size. Non-zero temperatures are incorporated via a purification of the density matrix. By using the disentangler introduced in Ref. [27] and exploiting time translation invariance S z j (t)S z j (0) = S z j (t/2)S z j (−t/2) [29]    we can substantially extend the accessible time scale. Details of the algorithm are described in Refs. [28,36]. By varying both the system size and the bound for the maximally discarded weight we ensure that the finite size and truncation errors are smaller than the symbol size for all data presented in the following. We find a striking difference between the weakly interacting case, Fig. 1, and the strongly interacting case, Fig. 2. While the data for ∆ = 0.8 can be very well fitted by a sum of the single impurity contribution (9) and the diffusive part (12), it is necessary to also include the twoimpurity contribution (10) for ∆ = 0.3. In the latter case, we also allow for a decay rate ρ by multiplying Eq. (10) by e −ρt . The fits yield very small decay rates which seem to be of order ∼ e −1/T and could possibly be related to thermal excitations at the band edges which we have neglected in our analysis [48]. At T = 0.25 we find for both ∆ values clear evidence for spin diffusion with diffusion constants Γ fit close to the theoretically predicted values, Γ th , see Eq. (12). In agreement with theory we also find numerically that the diffusive contribution seems to vanish for magnetic fields h T (data not shown) where the Umklapp term (11) is oscillating and should be dropped from the effective theory [42]. The oscillation frequencỹ ε(T ) and exponentη(T ), obtained from fits of tDMRG data at ∆ = 0.3, are shown in Fig. 3. The data confirm thatε ∼ −T 2 andη ∼ −T while the prefactors, even for ∆ = 0.3, already seem to deviate significantly from the lowest order result, Eq. (8).
Exact parameters-The integrability of the XXZ model raises the question whether parameters such as the frequencyε(T ) can be determined exactly, beyond the lowest order in Eq. (7). At T = 0 the exact ε is identified with the half bandwidth of the elementary excitations computable by BA [21]. The natural extension to T > 0 should be based on the energies of excitations on top of an equilibrium state in the thermodynamic Bethe ansatz (TBA) approach [54,55]. However, by solving TBA nonlinear integral equations numerically [48] we find that the bandwidth of the dressed energy for a single-hole excitation increases with T , whereas the perturbative expression (7) and the tDMRG results in Fig. 3 show thatε(T ) decreases with T . This disagreement is rather puzzling given that a generalized TBA approach has been shown to be applicable even to nonequilibrium dynamics [56].
Conclusion-Finite temperatures and quantum fluctuations lead to large non-perturbative effects on dynamical correlations in NLL's. The exponents of oscillating contributions are renormalized by temperature while umklapp scattering leads to spin diffusion dominating the long-time asymptotics. These predictions are in excellent agreement with tDMRG calculations which provide, in particular, the first direct evidence for spin diffusion in the XXZ model and show striking changes in the oneand two-impurity contributions as a function of temper-ature and interaction strength. The theory is easily extended to models such as the Bose gas and can also be used to study the propagation of an impurity through a 1D quantum fluid at finite temperatures. We expect, in particular, that these results will help to set up and interpret experiments in ultracold atoms aiming at an observation of spin diffusion [3,5]. Supplemental Material: "Low temperature dynamics of nonlinear Luttinger liquids" by C. Karrasch, R. G. Pereira, and J. Sirker

Irrelevant impurity-boson interactions
In this section we discuss how discrete symmetries and integrability of the XXZ model constrain the coupling constants of irrelevant operators in the effective impurity model.
Particle-hole (C) and parity (P ) transformations act on the bosonic fields and impurity fields as follows [52] : C : . (1) The Hamiltonian can only contain operators that are invariant under C and P . Let us denote by H (n) a All the operators in Eq. (2) conserve the total number of particles and holes in the high-energy subbands. In the above list we have omitted operators which couple the impurity modes to umklapp type operators. The lowest dimension operator in this family is (d † d +d †d ) cos(4 √ πKφ), whose scaling dimension varies continuously from 5 at ∆ = 0 to 3 at ∆ = 1. The set of constraints on the coupling constants that we shall derive in the following is not affected by this family of operators.
The integrability of the XXZ model affects the effective impurity model by constraining the coupling constants of irrelevant interactions. Following [15], we shall examine the consequences of integrability by imposing the existence of nontrivial conservation laws. For the XXZ model, it is fortunate that the first nontrivial conserved quantity can be identified with the energy current operator J E , which is defined from the continuity equation for the Hamiltonian density [49]. In the continuum limit, the energy current density is given by where H = dx H(x) = a dx H a is the total Hamiltonian. The energy current operator is J E = dx j(x). Let us denote by j a,b , with a < b, the contribution to the energy current density obtained by taking the commutator of terms H a and H b in the Hamiltonian as follows: The notation implies that when we take the commutator of a dimension-n operator with another dimension-m operator, the corresponding contribution to J E (when nonvanishing) has dimension n + m − 2.
To check that J E = a<b J a,b is conserved, we need to take the commutator with all the terms in H again: We organize the expansion by operator dimension. In order to find nontrivial relations between the coupling constants in Eq. (2), it suffices to compute [J E , H] to the level of dimension-four operators. For this we need to consider up to dimension-three operators in J E . We find the following list of operators The calculation of the commutators in Eq. (5) is tedious but straightforward. To simplify the result, we use the known relations for the XXZ model ε = 1/m = v. We find that the conservation law [J E , H] = 0 imposes the constraints Most importantly, integrabiliy rules out the g interaction. This is precisely the operator considered in [50] which accounts for a finite decay rate of a mobile impurity in a Luttinger liquid at finite temperatures.

Excitation energies from the thermodynamic Bethe ansatz
In this section we describe the calculation of dressed energies using the thermodynamic Bethe ansatz (TBA) approach to the XXZ model [54,57].
Bethe ansatz states are parametrized by a set of rapidities {x α j } which satisfy the Bethe equations. Here j labels the type of string and α specifies a particular rapidity. More precisely, x α j refers to the real part of the rapidity, since the imaginary part is fixed by the string hypothesis [54,57]. Strings with length n j > 1 are interpreted as bound states of n j particles. For simplicity, we choose the anisotropy parameter to be ∆ = cos(π/ν), with ν ∈ Z. In this case, we can restrict ourselves to a finite number of strings j = 1, 2, . . . , ν. The strings with j = 1, 2, . . . , ν − 1 have length n j = j and parity υ j = +1; the string with j = ν has length n ν = 1 and parity υ ν = −1.
The Bethe equations (in logarithmic form) for a chain of length N and periodic boundary conditions read where M j is the number of strings of type j in the Bethe ansatz state and I j α are integers (for M j odd) or half-integers (for M j even). A particular Bethe ansatz wave function is determined by the set of I j α 's. The functions t j (x) and the scattering phase shifts Θ jk (x) are given by where we define the function In the TBA approach, we take the limit N → ∞ and characterize the macroscopic state by the density of particles ρ(x) and density of holes ρ h (x) in rapidity space. The equilibrium state is obtained by minimizing the free energy as a functional of ρ(x) and ρ h (x). This leads to a set of coupled nonlinear integral equations for the dressed energies where β = 1/T is the inverse temperature and a j (x) = 1 2π Eq. (17) can be solved numerically by iteration. In the limit T → 0, the dressed energy for j = 1 (the even-parity one-string) reduces to the dispersion of the single-hole excitation over the ground state. The dressed energies can be used to calculate the free energy and other thermodynamic properties [57]. They also show up as the energies of elementary excitations over the equilibrium state [38]. The thermal excitation spectrum for the gapless phase of the XXZ model was calculated by Puga [55]. The energy required to create a single hole with rapidity x in the density of type-j strings is Notice that the excitation energies ∆E j (x) differ from the dressed energies in Eq. (17) by a constant term that involves all strings. Fig. 4 shows the excitation energy for the j = 1 string for three different values of temperature. We are particularly interested in the bandwidth, which is given by ∆E 1 (0). At T = 0, the bandwidth is known analytically, lim T →0 ∆E 1 (0) = π

Fits of the DMRG data
We have fitted the tDMRG data using the fit function where A, B,η,ε,φ, ρ, Γ are real fitting parameters. The first, constant term Γ is the diffusive contribution. The second term is the single impurity contribution, Eq. (9) in the main text, while the third term represents the two-impurity contribution, Eq. (10), where we have allowed for a small decay rate which seems to be of order ∼ e −1/T and might possibly be related to thermal excitations at the band edges. In table I we present the fit parameters for the fits shown in Fig. 2 and Fig. 3   One of our main findings based on the analysis of the numerical data is that the two-impurity contribution becomes very small for large interaction strengths leading to a fit parameter B for ∆ = 0.8 which is essentially zero. To further support that B for the considered temperatures is strongly reduced with interaction, we present in Fig. 5 below tDMRG data and fits for intermediate interaction strength ∆ = 0.5. Here a two-impurity contribution is still visible but the amplitude B is already very small, see table II. In order to illustrate the sensitivity of the fit parameters on the fit interval we concentrate on the case ∆ = 0.5, T = 0.167 and show in table II parameters for fits using three different time intervals. Except for the phase shiftφ, and, to a lesser extent, the amplitudeÃ, all fit parameters show little variation implying, in particular, that it is possible to extract the temperature dependence ofε(T ) andη(T ) with reasonably accuracy from the tDMRG data. On the other hand, we want to emphasize that the phase shiftφ cannot be fixed reliably from numerical data even at zero temperature, see Refs. [21,23].    Table II: Parameters for various fits of the tDMRG data at ∆ = 0.5.