Time evolution of the Luttinger model with nonuniform temperature profile

TL;DR

This paper analyzes the time evolution of a one-dimensional Luttinger model with a nonuniform temperature profile, deriving exact and approximate expressions for energy, heat current, and correlations, revealing universal and nonuniversal features.

Contribution

It provides analytical formulas for the nonequilibrium dynamics of the Luttinger model with nonuniform temperature, including exact results for local interactions and first-order approximations for nonlocal interactions.

Findings

Exact expressions for energy density and heat current involving Schwarzian derivative.

Universal heat current in the steady state even with broken conformal invariance.

Predicted short-time peaks and dispersion effects consistent with numerical simulations.

Abstract

We study the time evolution of a one-dimensional interacting fermion system described by the Luttinger model starting from a nonequilibrium state defined by a smooth temperature profile $T(x)$. As a specific example we consider the case when $T(x)$ is equal to $T_L$ ($T_R$) far to the left (right). Using a series expansion in $\epsilon = 2(T_{R} - T_{L})/(T_{L}+T_{R})$, we compute the energy density, the heat current density, and the fermion two-point correlation function for all times $t \geq 0$. For local (delta-function) interactions, the first two are computed to all orders, giving simple exact expressions involving the Schwarzian derivative of the integral of $T(x)$. For nonlocal interactions, breaking scale invariance, we compute the nonequilibrium steady state (NESS) to all orders and the evolution to first order in $\epsilon$. The heat current in the NESS is universal even when conformal invariance is broken by the interactions, and its dependence on $T_{L,R}$ agrees with numerical results for the $XXZ$ spin chain. Moreover, our analytical formulas predict peaks at short times in the transition region between different temperatures and show dispersion effects that, even if nonuniversal, are qualitatively similar to ones observed in numerical simulations for related models, such as spin chains and interacting lattice fermions.

Equations

$$\mathcal{G}=\int dx\beta (x)\mathcal{H}(x),$$

$$⟨\mathcal{O}(t)⟩\equiv \mathrm{Tr}\hat{\rho }\mathcal{O}(t),$$

$$H=\sum _r\frac{1}{2}\int dx\left[:{\psi }_r^{+}(x)\left(-ir{v}_F{\partial }_x\right){\psi }_r^{-}(x):+\text{h.c.}\right]$$

$$+\lambda \sum _{r,r^{\prime }}\int dxdyV(x-y):{\psi }_{r^{}}^{+}(x){\psi }_{r^{}}^{-}(x)::{\psi }_{r^{\prime }}^{+}(y){\psi }_{r^{\prime }}^{-}(y):,$$

$$H_r=\frac{1}{2}\int dqv(q):{\tilde{\rho }}_r(-q){\tilde{\rho }}_r(q):,$$

$$\underset{t\to \infty }{\lim }\mathrm{Tr}\hat{\rho }(t)\mathcal{O}=\frac{\mathrm{Tr}{e}^{-{\beta }_{+}{H}_{+}-{\beta }_{-}{H}_{-}}\mathcal{O}}{\mathrm{Tr}{e}^{-{\beta }_{+}{H}_{+}-{\beta }_{-}{H}_{-}}},$$

$$\underset{t\to \infty }{\lim }E(x,t)$$

$$\underset{t\to \infty }{\lim }J(x,t)$$

$$\lim_{t\to\infty}S_{r}(\xi,\tau,x,t)=\frac{i}{2\pi u_{r}}\exp\biggl{(}\int_{\mathbb{R}^{+}}\frac{dq}{q}\,\Bigl{\{}e^{iqu_{r}(q)}-e^{iqu_{r}}\Bigr{\}}\biggr{)}\\ \begin{aligned} &\times\exp\biggl{(}\int_{\mathbb{R}^{+}}\frac{dq}{q}\,\sinh^{2}\varphi(q)\\ &\qquad\quad\times\Bigl{\{}e^{iqu_{r}(q)}+e^{iqu_{-r}(q)}-2e^{-q0^{+}}\Bigr{\}}\biggr{)}\\ &\times\exp\biggl{(}\int_{\mathbb{R}^{+}}\frac{dq}{q}\,\biggl{[}\cosh^{2}\varphi(q)\frac{2\{\cos[qu_{r}(q)]-1\}}{e^{\beta_{r}\omega(q)}-1}\\ &\qquad\quad+\sinh^{2}\varphi(q)\frac{2\{\cos[qu_{-r}(q)]-1\}}{e^{\beta_{-r}\omega(q)}-1}\biggr{]}\biggr{)},\end{aligned}$$

$$\underset{t\to \infty }{\lim }J(x,t)=\sum _{r}r\frac{\pi T_r^2}{12}=\frac{\pi }{12}(T_L^2-T_R^2)\equiv J,$$

$$\underset{t\to \infty }{\lim }E(x,t)=\sum _r\frac{\pi }{12v}{T}_r^2=\frac{\pi }{12v}(T_L^2+T_R^2)$$

$$\lim_{t\to\infty}S_{r}(\xi,\tau,x,t)\\ =\frac{1}{2\pi\tilde{\ell}}\biggl{(}\frac{i\pi T_{r}\tilde{\ell}/v}{\sinh(\pi T_{r}u_{r}/v)}\biggr{)}^{1+\eta/2}\biggl{(}\frac{i\pi T_{-r}\tilde{\ell}/v}{\sinh(\pi T_{-r}u_{-r}/v)}\biggr{)}^{\eta/2},$$

$$E(x,t)$$

$$J(x,t)$$

$$G(x)$$

$$=\frac{\pi }{6v}T(x{)}^2-\frac{v}{12\pi }\left(\frac{T^{\prime \prime }(x)}{T(x)}-\frac{3}{2}{\left(\frac{T^{\prime }(x)}{T(x)}\right)}^2\right).$$

$$G(x)=\frac{\pi T^2}{6v}{g}^{\prime }(x{)}^2-\frac{v}{12\pi }(Sg)(x),$$

$$(Sg)(x)=\frac{g^{\prime \prime \prime }(x)}{g^{\prime }(x)}-\frac{3}{2}{\left(\frac{g^{\prime \prime }(x)}{g^{\prime }(x)}\right)}^2.$$

$$E(z,\bar{z})$$

$$J(z,\bar{z})$$

$$\mathcal{T}(z)={\left(\frac{dw}{dz}\right)}^2\mathcal{T}(w)+\frac{cv}{12}(Sw)(z),$$

$$E(x,t)$$

$$J(x,t)$$

$$E_1(x,t)$$

$$J_1(x,t)$$

$$\begin{gathered} A(p_1,p_2)=e^{i(p_1+p_2)x-i[r_1\omega (p_1)+r_2\omega (p_2)]t} \\ \begin{array}{rl}& \times \frac{[r_{1}v(p_1)+r_{2}v(p_2){]}^2}{4v(p_1)v(p_2)}\frac{[r_1{e}^{2\phi (p_1)}+r_2{e}^{2\phi (p_2)}{]}^2}{4e^{2[\phi (p_1)+\phi (p_2)]}}\\ & \times \frac{e^{\beta [r_1\omega (p_1)+r_2\omega (p_2)]}-1}{r_1\omega (p_1)+r_2\omega (p_2)}\frac{r_1\omega (p_1)}{e^{\beta r_1\omega (p_1)}-1}\frac{r_2\omega (p_2)}{e^{\beta r_2\omega (p_2)}-1}.\end{array} \end{gathered}$$

$$S_r(\xi ,\tau ,x,t)=⟨{\psi }_r^{+}(\xi ,\tau ){\psi }_r^{-}(0,0){⟩}_{\beta }{e}^{ϵB_{1;r}(\xi ,\tau ,x,t)+O({ϵ}^2)},$$

$$B_{1;r}(\xi ,\tau ,x,t)=-\sum _{r_1,r_2}{\text{strokedint}}_{\mathbb{R}}\frac{dp}{2\pi }{\int }_{\mathbb{R}}\frac{dq}{4\pi }\hat{W}(p)C(p-q,q),$$

$$\begin{gathered} C(p_1,p_2)=2\pi e^{i(p_1+p_2)x-i[r_1\omega (p_1)+r_2\omega (p_2)]t} \\ \begin{array}{rl}& \times \frac{r_{1}v(p_1)+r_{2}v(p_2)}{2}\frac{r_1{e}^{\phi (p_2)-\phi (p_1)}+r_2{e}^{\phi (p_1)-\phi (p_2)}}{}\\ & \times \frac{e^{\beta [r_1\omega (p_1)+r_2\omega (p_2)]}-1}{r_1\omega (p_1)+r_2\omega (p_2)}{F}_r^{r_1}(p_1)F_r^{r_2}(p_2)\end{array} \end{gathered}$$

$$F_r^{r^{\prime }}(p^{\prime })=\frac{e^{-\phi (p^{\prime })}+r{r}^{\prime }{e}^{\phi (p^{\prime })}}{2}\frac{e^{i{p}^{\prime }{r}^{\prime }{u}_{r^{\prime }}(p^{\prime })}-1}{e^{\beta r^{\prime }\omega (p^{\prime })}-1}.$$

$${\partial }_{\beta }U(\beta )=-e^{\beta H}\mathcal{W}{e}^{-\beta (H+\mathcal{W})}=-\mathcal{W}(\beta )U(\beta ),$$

$$⟨\mathcal{O}⟩=⟨\mathcal{O}{⟩}_{\beta }-ϵ\left[⟨\mathcal{C}\mathcal{O}{⟩}_{\beta }-⟨\mathcal{C}{⟩}_{\beta }⟨\mathcal{O}{⟩}_{\beta }\right]+O({ϵ}^2),$$

$$\frac{\operatorname{Tr}[e^{-\beta d\hat{\Gamma}(K+W)}d\hat{\Gamma}(O)]}{\operatorname{Tr}[e^{-\beta d\hat{\Gamma}(K+W)}]}-\frac{\operatorname{Tr}[e^{-\beta d\hat{\Gamma}(K)}d\hat{\Gamma}(O)]}{\operatorname{Tr}[e^{-\beta d\hat{\Gamma}(K}]}\\ \begin{aligned} &=\operatorname{tr}\bigl{(}\bigl{\{}[e^{2\beta(K+W)}-1]^{-1}-[e^{2\beta K}-1]^{-1}\bigr{\}}O\bigr{)}\\ &=\frac{1}{\beta}\sum_{\nu}\operatorname{tr}\bigl{(}\bigl{\{}[i\nu-2(K+W)]^{-1}-[i\nu-2K]^{-1}\bigr{\}}O\bigr{)}\\ &=\sum_{n=1}^{\infty}\frac{1}{\beta}\sum_{\nu}\operatorname{tr}\bigl{(}[i\nu-2K]^{-1}(2W[i\nu-2K]^{-1})^{n}O\bigr{)},\end{aligned}$$

$$G_{n}(x)=\int_{\mathbb{R}^{n+1}}\frac{dp_{0}\ldots dp_{n}}{(2\pi)^{n+1}}\,\biggl{(}\prod_{j=0}^{n-1}\hat{W}(p_{j}-p_{j+1})\biggr{)}\\ \times\frac{1}{2}\sum_{r}\frac{1}{\beta}\sum_{\nu}\biggl{(}\prod_{j=0}^{n}\frac{rvp_{j}}{i\nu-rvp_{j}}\biggr{)}e^{i(p_{0}-p_{n})x},$$

Full text

I Introduction

Experiments on ultracold atomic gases have led to renewed interest in nonequilibrium properties of isolated one-dimensional quantum systems EFG ; BDZ ; PSSV ; RDO ; LanEtAl ; EsFa . This field also has roots in a rich history of theoretical works studying both classical RLL ; SpLe ; Dh ; BLRb ; LLP ; BO ; BOS and quantum systems KaFi ; Rue ; JP1 ; JP2 ; AJPP ; BJP ; ARRS ; LaMi ; SaMi ; SPA ; PVC out of equilibrium. One often studied protocol is to join, at time $t=0$, disconnected left and right parts of an infinite system, where each part is in thermal equilibrium with temperatures $T_{L}$ and $T_{R}$, respectively. For $t>0$ the system is evolved with a fully translation invariant Hamiltonian; this produces a heat current and, for long times, the system tends to a nonequilibrium steady state (NESS) if $T_{L}\neq T_{R}$. This is usually referred to as the partitioning protocol.

Using the above protocol, exact results for the NESS were obtained for simple integrable models such as the $XX$ and $XY$ spin chains using $C^{*}$-algebraic methods HoAr ; AsPi ; Og1 ; Og2 and nonequilibrium Green’s functions (Keldysh formalism) ALA . When written in terms of fermions, these models are all noninteracting: they can be mapped to one-dimensional systems of spinless lattice fermions with Hamiltonians that are quadratic in the fermion fields. For general systems of free lattice fermions, results for the NESS were obtained using a generalized Landauer-Büttiker formula in AJPP ; BJP .

For interacting fermions the partitioning protocol was successfully used to obtain exact results for systems described by conformal field theories (CFTs) BeDo1 ; BeDo2 ; BeDo3 ; SoCa ; HoLo . Beyond CFT there are otherwise few exact results for the NESS, and even fewer for the evolution, of interacting fermions; see, e.g., Ca ; IuCa ; RSM ; KRSM ; MaWa ; LLMM ; CaCh . Using the same protocol, the time evolution and properties of the NESS have been studied extensively numerically KIM ; VKM ; BLVRMF and by approximate analytical methods LiAn ; LVBD ; LVMR in various models. Recently, effective hydrodynamic equations for the long-time and large-distance dynamics for Bethe ansatz-solvable models were proposed BCNF ; CBT ; see also Bre ; Spo . We also mention recent studies of the heat current and the thermal Drude weight based on Bethe ansatz Zot , density matrix renormalization group Kar , and hydrodynamics IlNa ; BVKM1 ; BVKM2 .

Most results for systems of interacting fermions, such as those mentioned above, rely on approximate methods or on assumptions, and it is thus interesting to obtain exact results for specific models that can serve as a benchmark. In this paper, we present some exact results for the full time evolution (not just the NESS) of a continuum system of interacting fermions described by the Luttinger model To ; Th ; Lu ; MaLi on the real line starting, at $t=0$, from a nonequilibrium state defined by a smooth temperature profile $T(x)$. This is related to but different from the partitioning protocol. Specifically, if $\mathcal{H}(x)$ is the energy density operator defining the Hamiltonian, $H=\int dx\,\mathcal{H}(x)$, then the initial state is given by $\hat{\rho}=e^{-\mathcal{G}}/\operatorname{Tr}e^{-\mathcal{G}}$, with

[TABLE]

where $\beta(x)\equiv T(x)^{-1}=\beta[1+\epsilon W(x)]$ for some smooth function $W(x)$, with $\beta$ the average inverse temperature and $\epsilon$ the distance from equilibrium. (We use units such that $\hbar=k_{B}=1$.) We will mainly be concerned with the case of a step-like profile $T(x)$ equal to $T_{L}$ ($T_{R}$) far to the left (right), e.g., $W(x)=-(1/2)\tanh(x/\delta)$ with $\delta>0$, where $\beta$ and $\epsilon$ are determined by $\beta(\mp\infty)=T_{L,R}^{-1}$. The evolution of the system is given by $H$, and we are interested in nonequilibrium expectation values ($\epsilon\neq 0$) of local observables $\mathcal{O}$,

[TABLE]

where $\mathcal{O}(t)=e^{iHt}\mathcal{O}e^{-iHt}$. If $\epsilon=0$, then $\langle\mathcal{O}(t)\rangle=\langle\mathcal{O}\rangle_{\beta}$ is an equilibrium expectation value with temperature $T=\beta^{-1}$. For the Luttinger model, such equilibrium properties are well known, since a long time, from the celebrated exact solution in MaLi using bosonization; see also, e.g., Hal ; Gia ; GNT ; Kop ; HSU ; MEJ ; LaMo ; Voit .

We use a series expansion in $\epsilon$ to compute the time evolution and the NESS for the Luttinger model both in the case of local (delta-function) and nonlocal interactions starting from a nonequilibrium state. We show that the NESS is factorized in terms of the eigenmodes of the interacting Hamiltonian (plasmons) MaLi and not in terms of the fermions; the presence of interactions is manifested by interaction-dependent exponents in the fermion two-point correlation function. In contrast, we find that the final heat current is universal even for interactions that break conformal invariance, and the form of its dependence on $T_{L,R}$ confirms previous numerical results for interacting lattice fermion models, such as the $XXZ$ spin chain studied in KIM . For noninteracting lattice models, such results for the temperature dependence were obtained analytically in HoAr ; AJPP ; BJP ; ALA .

For local interactions, our series for the energy and heat current densities can be summed into exact formulas for the time evolution. These results contain the Schwarzian derivative FMS of the integral of $T(x)$, which is very suggestive in view of the conformal invariance in the local case. Its presence produces peaks in the energy and heat current densities at zero time in the transition region between different temperatures. These resemble what is found numerically in related models KIM ; LVBD , even if the shapes of such peaks clearly are nonuniversal.

For nonlocal interactions, breaking conformal invariance, we obtain analytical results for the NESS to all orders and for the time evolution to first order in $\epsilon$. In this case, dispersive effects appear in the evolution, which look qualitatively similar to those seen numerically in lattice models. (Such dispersive effects are absent for local interactions.)

The following two methods are used to compute nonequilibrium expectation values: Method 1, based on the Dyson series, and Method 2, using one-particle operators; see Secs. V.1 and V.2, respectively. Method 1 allows one to compute nonequilibrium results to first order in $\epsilon$ from equilibrium ones, and it can be used even for non-exactly-solvable models. Method 2 allows one to compute results for the Luttinger model to all orders in $\epsilon$, and it is in general applicable only to models that are quasi-free.

We consider the Luttinger model given by

[TABLE]

with fermion fields $\psi_{r}^{-}(x)$ and $\psi^{+}_{r}(x)=\psi^{-}_{r}(x)^{\dagger}$, where $r=+(-)$ denotes right- (left-) moving fermions, $\left.:\!\cdots\!:\right.$ indicates Wick (normal) ordering, $v_{F}>0$ is the Fermi velocity, $V(x)$ is the interaction potential, and $\lambda$ is the coupling constant. We use notation similar to LLMM ; MaLi ; cf. also MaWa ; LaMo and references therein. Let $\hat{V}(p)=\int dx\,V(x)e^{-ipx}$ denote the Fourier transform of the potential. The interactions must satisfy $\lambda\hat{V}(p)>-\pi v_{F}/2$, and $V(x)$ can be local, $V(x)=\pi v_{F}\delta(x)/2$, which requires renormalizations, or nonlocal with an interaction range $a>0$, e.g., $V(x)=\pi v_{F}/[4a\cosh(\pi x/2a)]$. The above examples of potentials are used in Figs. 1 and 2 below to illustrate our analytical results, but we emphasize that these results hold true for a large class of interactions LLMM ; LaMo .

In what follows, we study the evolution of the energy density $E(x,t)\equiv\langle\mathcal{H}(x,t)\rangle$, the heat current density $J(x,t)\equiv\langle\mathcal{J}(x,t)\rangle$, and the fermion two-point correlation function $S_{r}(\xi,\tau,x,t)\equiv\langle\psi^{+}_{r}(x+\xi,t+\tau)\psi^{-}_{r}(x,t)\rangle$, where $\mathcal{J}(x,t)$ is determined by the continuity equation $\partial_{t}\mathcal{H}(x,t)+\partial_{x}\mathcal{J}(x,t)=0$. We start in Sec. II by presenting results for the NESS. This serves as a useful benchmark for the finite-time results presented in Sec. III for local interactions and in Sec. IV for nonlocal interactions. Our methods are described in Sec. V, and concluding remarks are given in Sec. VI. Some computational details are deferred to the Appendix.

II Nonequilibrium steady state

It is well known that the Fourier modes of the fermion densities, $\rho_{r}(p)\equiv\int dx\,\left.:\!\psi^{+}_{r^{\vphantom{\prime}}}(x)\psi^{-}_{r^{\vphantom{\prime}}}(x)\!:\right.e^{-ipx}$, define boson operators MaLi , and that the Luttinger Hamiltonian can be written as $H=H_{+}+H_{-}$, with

[TABLE]

using Bogoliubov transformed fermion densities $\tilde{\rho}_{r}(p)=\rho_{r}(p)\cosh\varphi(p)-\rho_{-r}(p)\sinh\varphi(p)$, where $\tanh 2\varphi(p)=-\lambda\hat{V}(p)/[\pi v_{F}+\lambda\hat{V}(p)]$, and the renormalized Fermi velocity $v(p)=v_{F}\sqrt{1+2\lambda\hat{V}(p)/\pi v_{F}}$ MaLi ; LLMM ; Voit . The $\tilde{\rho}_{r}(p)$ are commonly referred to as plasmons, and the Luttinger Hamiltonian is diagonal in terms of these MaLi . To find the NESS, we write $\hat{\rho}(t)=e^{-iHt}\hat{\rho}e^{iHt}=e^{-\mathcal{G}(-t)}/\operatorname{Tr}(e^{-\mathcal{G}(-t)})$ with $\mathcal{G}(t)=\int dx\,\beta(x)\mathcal{H}(x,t)$ and express $\mathcal{H}(x,t)$ in terms of $\tilde{\rho}_{r}(p,t)=\tilde{\rho}_{r}(p)e^{-ir\omega(p)t}$, where $\omega(p)\equiv v(p)p$. Taking $t\to\infty$ in $\hat{\rho}(t)$ by making use of the Riemann-Lebesgue lemma (cf., e.g., LLMM ), which can be justified for expectation values using Method 2, we find

[TABLE]

with $\beta_{+}=T_{L}^{-1}$ and $\beta_{-}=T_{R}^{-1}$. This NESS describes a translation invariant state factorized into right- and left-moving plasmons at equilibrium with temperatures $T_{\pm}=\beta_{\pm}^{-1}$. A similar NESS was obtained in BeDo1 ; BeDo2 ; BeDo3 for CFTs and in HoAr ; AsPi ; Og1 ; Og2 for the $XX$ chain; in the latter case the same factorization of the NESS is valid also in terms of right- and left-moving fermions, whereas in our case only the plasmons factorize in such a way but not the fermions.

The long-time limit of expectation values for all local observables can be computed using (5) by straightforward generalizations of well-known equilibrium computations. By recalling that $\int dx\,\mathcal{H}(x)=\sum_{r}H_{r}$ with $H_{r}$ in (4) and using the continuity equation to show that $\int dx\,\mathcal{J}(x)=\sum_{r}({r}/{2})\int dq\,({d\omega(q)}/{dq})v(q)\!\left.:\!\tilde{\rho}_{r}(-q)\tilde{\rho}_{r}(q)\!:\right.$, we obtain

[TABLE]

where $w_{\lambda}$ is the ground state energy density MaLi ; LLMM , using that the NESS is translation invariant. Similarly, for the fermion two-point correlation function, using the well-known bosonization formula expressing fermions as exponentials of plasmons (see, e.g., HSU ; MEJ ; LLMM ; LaMo and references therein), we find

[TABLE]

where $u_{r}(p)\equiv r[\xi-rv(p)\tau]+i0^{+}$ and $u_{r}\equiv u_{r}(0)$.

The second integral in (6) gives the final energy flow and appears to depend on the interactions. However, by the change of variables $u=\beta_{r}\omega(q)$, we obtain

[TABLE]

due to the presence of the group velocity $d\omega(q)/dq$ in the integrand [assuming $d\omega(q)/dq>0$, which is true for a large class of interaction potentials LLMM ]. It follows that the final heat current only depends on $T_{L,R}$ and is independent of microscopic details. Such universal behavior, previously observed in CFTs BeDo1 ; BeDo2 ; BeDo3 , thus remains true for the Luttinger model even when scale invariance is broken by the interactions. This result supports the conjecture, based on numerical simulations of the $XXZ$ chain KIM , that for interacting fermions, $J=f(T_{L})-f(T_{R})$, where, in general, $f$ is a nonuniversal function tending to the universal CFT result BeDo1 in the low-temperature limit.

For noninteracting fermions, the temperature dependence $J=f(T_{L})-f(T_{R})$ corresponds to the above-mentioned factorization of the NESS and was previously obtained analytically by different methods HoAr ; AJPP ; BJP ; ALA . In fact, using these analytical results, the function $f$ for the $XX$ chain can be computed analytically, $f(T)=({\pi}/{12})T^{2}\left[1-R(b_{+})-R(b_{-})\right]$, with nonuniversal corrections $R(b_{\pm})=(6/\pi^{2})\int_{b_{\pm}}^{\infty}dx\,x/(e^{x}+1)$ governed by $b_{\pm}=({v_{F}}/{Ta_{0}})[{1\pm\cos(\nu\pi)}]/{\sin(\nu\pi)}$, where $a_{0}$ is the lattice spacing and $0<\nu<1$ is the filling factor (specifying the Fermi momentum $k_{F}=\nu\pi/a_{0}$). If $Ta_{0}/v_{F}$ is small, the corrections are exponentially suppressed, and the universal result becomes exact in the scaling limit $Ta_{0}/v_{F}\to 0$.

The first integral in (6) expresses the energy density in the NESS as a sum of energy densities at equilibrium with temperatures $T_{L,R}$ and is nonuniversal. Indeed, it depends on the interactions, and only in the local case, when $v(p)=v\equiv v(0)$ and $\varphi(p)=\varphi$ are constant, does it simplify to

[TABLE]

after an additive renormalization corresponding to subtracting the (diverging) constant $w_{\lambda}$. Similarly, the two-point correlation function in the local case, after a multiplicative renormalization of the fermion fields (not needed in the nonlocal case), becomes

[TABLE]

where $u_{r}=r[\xi-rv\tau]+i0^{+}$, with the equilibrium anomalous exponent $\eta=2\sinh^{2}\varphi$ MaLi and a length parameter $\tilde{\ell}$ due to the renormalization; cf. also LLMM ; LaMo . This exponent depends on the interactions and is nonzero if the interactions are nonzero. Clearly, unless $\eta=0$, the NESS does not factorize into right- (left-) moving fermions with temperatures $T_{L}$ ($T_{R}$) as for the $XX$ chain.

III Finite-time results: Local interactions

The Luttinger model with local interactions is conformally invariant, implying that $\mathcal{H}(x,t)$ and $\mathcal{J}(x,t)$ satisfy the wave equation, and thus

[TABLE]

for some function $G(x)$. Using Method 2 $G(x)$ can be computed as a series expansion in $\epsilon$ to all orders (see the Appendix), and, after summation, we obtain the following remarkably simple result:

[TABLE]

The term $({\pi}/{6v})T(x)^{2}$ is expected from the equilibrium result for a uniform temperature profile, but the presence of the derivative terms has apparently been overlooked in the previous literature. Thus, in the case of a nonuniform temperature profile, (11) and (12) show that $E(x,t)$ and $J(x,t)$ also depend on the first and second derivatives of $T(x)$, but not on higher-order ones. This is true even at $t=0$.

The evolution of the energy flow can be easily understood using (11) and (12). For a step-like $\beta(x)=\beta[1+\epsilon W(x)]$ with $W(x)=-(1/2)\tanh(x/\delta)$, as in the Introduction, the energy profile at $t=0$ away from $x=0$ is essentially proportional to the local temperature, i.e., $E(x,0)$ equals $({\pi}/{12v})T_{L,R}^{2}$ far to the left and right. However, in the transition region, for small $\delta>0$ and $\epsilon\neq 0$, the derivative terms in (12) produce peaks; see Fig. 1(a). As $t$ increases, a region develops around the origin with a uniform energy density bounded by two rigid fronts (their shape does not change with time) that move ballistically to the right (left) with constant velocity $v$ ($-v$); see Fig. 1(b). In the same region the current has a nonvanishing constant value. For large times we recover the results for the NESS in (8) and (9).

As we discuss in Sec. IV, peaks qualitatively similar to those described above are seen in other related models, including interacting lattice models, such as the $XXZ$ chain, and noninteracting models, such as the $XX$ chain. It is important to stress that the shape of the peaks is nonuniversal and depends on short-distance details: this is clear already from the interaction dependence of the derivative terms that appear in (11) due to (12).

It is interesting to note that $G(x)$ can be written as

[TABLE]

using the function $g(x)=\int_{0}^{x}dx^{\prime}\,T(x^{\prime})/T$ and the so-called Schwarzian derivative FMS

[TABLE]

By recalling that the Luttinger model with local interactions is a CFT with central charge $c=1$, this result has a simple interpretation as follows. In a CFT, the energy and heat current densities are given by expectation values of the renormalized energy-momentum tensor $\mathcal{T}(z)$ and $\bar{\mathcal{T}}(\bar{z})$,

[TABLE]

using $z=x+iv\tau$ and $\bar{z}=x-iv\tau$, with $\tau$ denoting imaginary time FMS ; SoCa . Moreover, under a conformal transformation $z\to w(z)$, the renormalized energy-momentum tensor in a CFT transforms as $\mathcal{T}(z)\to\mathcal{T}(w)$, with

[TABLE]

using the Schwarzian derivative $S$ FMS . From the above one obtains (11) by a Wick rotation $\tau\to it$ and the identification $G(x)=-\pi^{-1}\langle\mathcal{T}(z)\rangle|_{z=x}$, using that $E(x,t)=E(z,\bar{z})$ and $J(x,t)=-iJ(z,\bar{z})$. Our results in (11) and (13) are therefore equivalent to what one would obtain by a conformal transformation determined by the function $g(x)=\int_{0}^{x}dx^{\prime}\,T(x^{\prime})/T$ from the equilibrium result $\langle\mathcal{T}(w)\rangle_{\beta}=\langle\bar{\mathcal{T}}(\bar{w})\rangle_{\beta}=-c\pi^{2}T^{2}/6v$ (for the latter see, e.g., BeDo3 ). As we discuss in Sec. VI, it would be interesting to check if this is true also for other observables.

IV Finite-time results: Nonlocal interactions

We now consider the Luttinger model with nonlocal interactions. Such interactions break conformal invariance and give rise to dispersion effects since the renormalized Fermi velocity $v(p)$ depends on momenta. These effects are qualitatively similar to ones observed in lattice models. (The interaction range introduces a scale similar to the lattice spacing.) We compute quantities only to first order in $\epsilon$ using Method 1. Comparison with our all-order results for the NESS and for finite times in the local case suggests that such a first-order approximation works well for small $\epsilon$: for example, for $\epsilon=-0.01$, used below in Figs. 1 and 2, first- and all-order results are practically indistinguishable, and thus the deviations seen in these figures between the plots for local and nonlocal interactions can be fully attributed to dispersive effects.

For the energy and heat current densities, we obtain

[TABLE]

where $E_{0}$ is equal to $\lim_{t\to\infty}E(x,t)$ in (6) for $\beta_{+}=\beta_{-}=\beta$,

[TABLE]

with

[TABLE]

Similarly, for the two-point correlation function, we obtain

[TABLE]

where $\langle\psi^{+}_{r}(\xi,\tau)\psi^{-}_{r}(0,0)\rangle_{\beta}$ is equal to $\lim_{t\to\infty}S_{r}(\xi,\tau,x,t)$ in (7) for $\beta_{+}=\beta_{-}=\beta$,

[TABLE]

with

[TABLE]

and

[TABLE]

The above results agree, to first order in $\epsilon$, with (6) and (7) as $t\to\infty$.

Similar to the discussion for the local case in Sec. III for a step-like $\beta(x)$, our analytical results in (17) and (18) show, for small $\delta>0$ and $\epsilon\neq 0$, that peaks are produced in the transition region between different temperatures; see Fig. 1(a). As $t$ increases, a region develops around the origin with a uniform energy density bounded by two ballistically moving nonrigid fronts (their shape changes with time); see Fig. 1(b).

In Fig. 2 we plot the current through $x=0$ as a function of time. The plotted results contain an initial peak. As seen from the dotted line in Fig. 2, such a peak is absent in the local case if the second term in (13) is omitted. A qualitatively similar peak is present in numerical results for the $XXZ$ chain; see, e.g., Fig. 1(a) in KIM and Fig. 3 in LVMR showing the heat current through the contact point in the partitioning protocol. As emphasized in Sec. III and also in KIM , the shape of such peaks is nonuniversal: in the Luttinger model the shape depends on the interactions and in the $XXZ$ chain on the anisotropy and the dispersion relation. However, the presence of the peaks seems to be a generic feature.

To further support our claim about the peaks, we also present, as an example for noninteracting lattice fermions, plots of the corresponding results for the $XX$ chain computed to first order in $\epsilon$ using Method 1; see Figs. 3 and 4. Peaks and dispersion effects that are qualitatively similar to the ones in Figs. 1 and 2 are clearly visible. Moreover, we checked numerically and analytically, to first order in $\epsilon$, that the results for the $XX$ chain approach those of the noninteracting Luttinger model in the scaling limit; plots of the latter are given by the red (plain) line in Figs. 3 and 4. This is true even at finite times.

V Methods

Our results are based on rigorous bosonization methods well known from studies of the Luttinger model in equilibrium; see, e.g., MaLi ; HSU ; LLMM ; LaMo . We work on the circle $-L/2\leq x\leq L/2$ of length $L>0$ with the fermion fields $\psi_{r}^{\pm}(x)$ satisfying anti-periodic boundary conditions and take the thermodynamic limit $L\to\infty$ only after computing expectation values for finite $t\geq 0$. The order, first $L\to\infty$ and then $t\to\infty$, is important for computing results in the long-time limit LLMM ; BeDo1 .

V.1 Method 1

To compute $\langle\mathcal{O}\rangle$, we write $\mathcal{G}$ in (1) as $\beta(H+\mathcal{W})$ with $\mathcal{W}=\epsilon\int dx\,W(x)\mathcal{H}(x)$ and use the fact that $U(\beta)\equiv e^{\beta H}e^{-\beta(H+\mathcal{W})}$ satisfies

[TABLE]

with $\mathcal{W}(\beta)=e^{\beta H}\mathcal{W}e^{-\beta H}$. Solving this by iteration we obtain a series expansion in $\epsilon$ (the Dyson series),

[TABLE]

with $\mathcal{C}=\int_{0}^{\beta}d\beta^{\prime}\int dx\,W(x)\mathcal{H}(x,-i\beta^{\prime})$. It follows that nonequilibrium expectation values are expressed as sums of equilibrium ones. This method can be used for computing nonequilibrium results to first order in $\epsilon$ for any model where equilibrium results are computable. Computations of the energy and heat current densities and the two-point correlation functions for the Luttinger model are straightforward but tedious using Wick’s theorem; the details will be presented elsewhere.

V.2 Method 2

Higher-order terms can be computed using general mathematical results for quasi-free models; see, e.g., GrLa . For the bosonized Luttinger Hamiltonian we write $H=d\hat{\Gamma}(K)$ to mean boson second quantization of the one-particle operator $K$, and similarly $\mathcal{W}=d\hat{\Gamma}(W)$ for some $W$. (We note that the second quantization map $d\hat{\Gamma}$ is in a nontrivial representation of the boson field algebra and that there are certain technical requirements on the one-particle operators GrLa ; NNS that are fulfilled in the cases of interest to us.) For $\mathcal{O}=d\hat{\Gamma}(O)$ with some one-particle operator $O$, one can show (e.g., using results in GrLa ) that $\langle\mathcal{O}\rangle-\langle\mathcal{O}\rangle_{\beta}$ can be written as

[TABLE]

where $\operatorname{tr}$ is the one-particle trace and the $\nu$ sum is over all boson Matsubara frequencies $\nu\in(2\pi/\beta)\mathbb{Z}$. [Note that the second and third identities in (23) are standard expansions.] The computation of $G(x)$ in (12) using (23) for the Luttinger model with local interactions is explained in the Appendix.

VI Conclusions

We derived analytical results for the NESS and for the full time evolution of the Luttinger model with both local and nonlocal interactions starting from a nonequilibrium state defined by a smooth nonuniform temperature profile. These results were computed using methods based on a series expansion in the distance $\epsilon$ from equilibrium in the initial state. We showed that the NESS is factorized in terms of the eigenmodes of the interacting Hamiltonian and that its fermion two-point correlation function contains interaction-dependent exponents. On the contrary, the final heat current is equal to the universal CFT result BeDo1 even if conformal invariance is broken by the interactions. Moreover, the form of the temperature dependence of the heat current agrees with the one found numerically in KIM for interacting fermions and analytically in HoAr ; AJPP ; BJP ; ALA for noninteracting fermions.

For local interactions (and thus a priori for the noninteracting case), the series for the energy and heat current densities were computed to all orders in $\epsilon$ and summed into simple exact formulas valid at all times. These formulas contain a Schwarzian-derivative term [cf. (11) and (13)], which captures a qualitative feature that appears rather generically, namely, the presence of nonuniversal peaks at short times in the transition region between different temperatures. We also showed that these formulas coincide with the result obtained by a particular conformal transformation from the corresponding equilibrium result. It would be interesting to find an explanation for this and to check if this is true also for other observables and in other CFT models; if true, this would be similar in spirit to results in SoCa but for a different physical situation. Also, it would be interesting to investigate if this can be used to gain some insight into nonequilibrium properties of interacting lattice models, such as the $XXZ$ chain.

For nonlocal interactions, we computed the time evolution of the energy and heat current densities and of the fermion two-point correlation function to first order in $\epsilon$. This truncated expansion can be seen as a linear-response approach (cf., e.g., KIM ) and can, in principle, be used even for models that are not exactly solvable.

Acknowledgements.

We would like to thank Natan Andrei, Jens H. Bardarson, Spyros Sotiriadis, and Herbert Spohn for valuable discussions. E.L. acknowledges support by VR Grant No. 2016-05167 and thanks the Erwin Schrödinger Institute (ESI) in Vienna for their hospitality during the workshop “Synergies between Mathematical and Computational Approaches to Quantum Many-Body Physics,” which provided useful input to this work. The work of J.L.L. is supported by AFOSR Grant No. FA-9550-16-1-0037 and NSF Grant No. DMR1104501. V.M. thanks Rutger’s University and the Institute for Advanced Study in Princeton, where parts of this work were done, for their kind hospitality. P.M. is grateful to the organizers of the 2016 EMS-IAMP Summer School in Mathematical Physics in Rome and thanks L’Università degli Studi di Roma “La Sapienza” for support during a visit that was important for the completion of this work.

Appendix: Computational details

For local interactions, the one-particle operators $K$ and $W$ in (23) are given by $K_{r,r^{\prime}}(p,p^{\prime})=(rvp/2)\delta_{r,r^{\prime}}\delta_{p,p^{\prime}}$ and $W_{r,r^{\prime}}(p,p^{\prime})=\epsilon[rv\operatorname{sgn}(p)\sqrt{|pp^{\prime}|}/2L]\delta_{r,r^{\prime}}\hat{W}(p-p^{\prime})$, respectively. Since $G(x)=E(x,0)$ [cf. (11)] it follows that $G(x)=\langle\mathcal{H}(x)\rangle=\langle d\hat{\Gamma}(O)\rangle$, with $O$ given by $O_{r,r^{\prime}}(p,p^{\prime})=[rv\operatorname{sgn}(p)\sqrt{|pp^{\prime}|}/2L]\delta_{r,r^{\prime}}e^{i(p^{\prime}-p)x}$. Using (23) we obtain $G(x)=\sum_{n=0}^{\infty}\epsilon^{n}G_{n}(x)$, where $G_{0}(x)=\pi/6v\beta^{2}$ is the equilibrium result and

[TABLE]

for $n=1,2,\ldots$. While this formula can be generalized to nonlocal interactions, the local case is special in that it is possible to compute the integrals exactly: changing variables to $q_{j}=p_{j-1}-p_{j}$ for $j=1,\ldots,n$ and $p=p_{n}$, and renaming $\nu\to r\nu$, we can write

[TABLE]

with

[TABLE]

where ${\bf q}=(q_{1},\ldots,q_{n})$ and $Q_{j}=\sum_{k=j+1}^{n}q_{k}$. The integral in (A3) can be computed exactly, and, after a lengthy computation, we obtained the following remarkably simple result,

[TABLE]

where $\simeq$ is defined by $q_{j}q_{k}\simeq q_{1}^{2}$ if $j=k$ and $q_{j}q_{k}\simeq q_{1}q_{2}$ if $j\neq k$. Inserting (A4) into (A2) yields

[TABLE]

Using this, the series $G(x)=\sum_{n=0}^{\infty}\epsilon^{n}G_{n}(x)$ can be summed, giving the result in (12).