Contents

Thermodynamic Bethe ansatz for non-equilibrium steady

[0.1cm] states: exact energy current and fluctuations in integrable QFT

Olalla Castro-Alvaredo, Yixiong Chen, Benjamin Doyon, Marianne Hoogeveen

Department of Mathematics, City University London, Northampton Square EC1V 0HB, UK

Department of Mathematics, King’s College London, Strand WC2R 2LS, UK

We evaluate the exact energy current and scaled cumulant generating function (related to the large-deviation function) in non-equilibrium steady states with energy flow, in any integrable model of relativistic quantum field theory (IQFT) with diagonal scattering. Our derivations are based on various recent results of D. Bernard and B. Doyon. The steady states are built by connecting homogeneously two infinite halves of the system thermalized at different temperatures , , and waiting for a long time. We evaluate the current using the exact QFT density matrix describing these non-equilibrium steady states and using Al.B. Zamolodchikov’s method of the thermodynamic Bethe ansatz (TBA). The scaled cumulant generating function is obtained from the extended fluctuation relations which hold in integrable models. We verify our formula in particular by showing that the conformal field theory (CFT) result is obtained in the high-temperature limit. We analyze numerically our non-equilibrium steady-state TBA equations for three models: the sinh-Gordon model, the roaming trajectories model, and the sine-Gordon model at a particular reflectionless point. Based on the numerics, we conjecture that an infinite family of non-equilibrium -functions, associated to the scaled cumulants, can be defined, which we interpret physically. We study the full scaled distribution function and find that it can be described by a set of independent Poisson processes. Finally, we show that the “additivity” property of the current, which is known to hold in CFT and was proposed to hold more generally, does not hold in general IQFT, that is is not of the form .

August 19, 2019

## 1 Introduction

In recent years, there has been a surge of interest in the study of the thermodynamics of quantum systems out of equilibrium (for reviews, see [1, 2]). In part this is due to the recent advances in experimental techniques, making it possible to drive quantum systems away from equilibrium in a controlled way, and to study their non-equilibrium properties, see for instance [3, 4, 5, 6, 7, 8]. This is also due to theoretical advances which include the discovery of several types of fluctuation theorems, generalizing the fluctuation-dissipation theorem [9] to systems far from equilibrium (see [1, 2] for extensive discussions). These fluctuation theorems express “universal” properties of certain distribution functions. In non-equilibrium steady states, where constant flows of particles, charge or energy exist, one object of interest is the scaled cumulant generating function (SCGF)1. The SCGF fully characterises the statistics of the transferred quantity at large times. Fluctuation theorems are symmetry relations for the SCGF; see for instance [10, 11, 12, 13, 14, 15] concerning the energy-flow SCGF.

In this manuscript we obtain for the first time the exact non-equilibrium current and SCGF for energy transfer in a wide family of integrable models. For this purpose, we propose, develop and test a new approach to the study of quantum integrable models in non-equilibrium steady states, concentrating on the study of non-equilibrium energy flows. We take the setup where two hamiltonian reservoirs at different temperatures are connected homogeneously to each other (local quantum quench), so that at large times an energy current be established. We consider the scaling limit: the quantum system is assumed to be in the universal regime near a quantum critical point (with unit dynamical exponent), with a mass gap and the driving temperatures assumed to be much smaller than any microscopic energy scale. This regime is described by massive quantum field theory (QFT). Our results are the first exact results for the energy SCGF and most general results for the energy current in non-critical models that cannot be described by free particles.

Integrability occurs when enough additional symmetries render the model exactly solvable [16, 17]. Exact solvability means that in principle all energy eigenstates and eigenvalues, and correlation functions of local operators, can be obtained non-perturbatively. Integrability has enjoyed much success over the past four decades. For integrable quantum spin chains several very effective approaches exist such as the algebraic Bethe ansatz (see e.g. [18]). For integrable massive quantum field theories (IQFTs) the problem of computing scattering amplitudes and vacuum correlation functions of local fields has been solved for many models using factorized scattering theory [19, 20, 21, 22, 23, 24]. Factorized scattering and Bethe ansatz ideas led to a very successful approach for the study of the equilibrium thermodynamic properties of integrable models proposed by Al.B. Zamolodchikov [25], known as the thermodynamic Bethe ansatz (TBA) approach. Other exact methods are based on free-fermion techniques (Clifford algebras), and the universal regime of exactly gapless quantum models is described by conformal field theory (CFT), for which a very wide variety of exact techniques are available [26, 27].

Using such techniques, exact results for SCGF in non-equilibrium integrable quantum systems have been obtained in the past, especially for charge flows, in which case the SCGF is referred to as the “full-counting statistics”. The first exact formula in free fermion models was obtained by Lesovik and Levitov in [28, 29] and further studied in other works [30, 31, 32, 33]. Other exact formulas were obtained in Luttinger liquids, which are particular models of CFT, in [34], and in the low-temperature universal regime of quantum critical models in [35] using general CFT. In certain integrable interacting impurity models, the SCGF was obtained in [36, 37] using TBA-like ideas. Exact charge current and shot noise (zero-temperature second cumulant) studies have also been done from similar ideas before in a variety of integrable models and with a variety of techniques, see for instance [38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50].

However for energy flows, exact results for currents in interacting models and for SCGF in general are more recent. As far as we are aware, the first exact SCGF was obtained for a chain of quantum harmonic oscillators in [51]. Exact results for the energy current and SCGF in the low-temperature universal regime of general quantum critical models were obtained very recently in [52, 35] using CFT, and for the quantum Ising chain in a magnetic field in [53] using free fermion techniques (results also exist for energy flows in more general “star-graph” configurations, see for instance [54, 55, 56]).

There are in fact many ways of theoretically constructing a non-equilibrium steady state (see [57] for a discussion); the aforementioned results involve a variety of setups. The hamiltonian-reservoir setup considered here is an old idea going back, in the quantum realm, to Keldysh [58] and Caroli et. al. [59], and studied in detailed classically in [60]. For energy flows, it was rigorously analyzed in the XY model in [61] based on [62], where the resulting non-equilibrium density matrix was constructed2, and likewise in CFT in [52, 35, 56]. The energy-flow non-equilibrium density matrix of the hamiltonian setup in general massive QFT was proposed and derived in [52, 66].

We present here a generalization of the TBA approach to non-equilibrium steady states energy flows with hamiltonian reservoirs: the NESSTBA approach. We base our derivation on the exact massive QFT density matrix derived in [66]. We deduce a set of coupled integral equations which are analog to the standard TBA equations up to some crucial modifications. Solutions to these equations give the exact energy current. All the exact cumulants are then obtained from the extended fluctuation relation that was derived in [15]. We illustrate the working of our approach by numerically evaluating the current and some of the cumulants in three different models: the sinh-Gordon model, the roaming trajectories model, and the sine-Gordon model at the simplest nontrivial reflectionless point. Our work can be seen as a generalization of the CFT results of [52] to IQFTs. In particular, as a consistency check we show that our formulation indeed leads to the known CFT result for the current and SCGF when both right and left temperatures are very high as compared to the gap.

The analysis of our results bring out three interesting observations. First, we find a family of non-equilibrium analogs of Zamolodchikov’s -functions [25], in terms of the current and all higher cumulants. We numerically verify that our non-equilibrium -functions are monotonic along RG trajectories, and equal to the central charge at fixed points. Hence, following the usual wisdom, they are (non-equilibrium) measures of the number of degrees of freedom as functions of the energy scale. Our -functions in turn give rise to exact upper bounds on the current and all higher cumulants. Second, we show numerically that the process of energy transfer can be understood, in its large-time regime, as a family of independent Poisson processes, one for each energy value and direction of transfer. This generalizes nontrivially what was observed in CFT in [52]. It indicates that the generalization of the usual Landauer form of the conductance to the full out of equilibrium statistics is not straightfoward: the weights of the Poisson processes are not simply related to the state density and occupation (contrary to the CFT case). Finally, it is known that in CFT the current is additive [52]: it satisfies (where are inverse temperatures), equivalently for some function . This was suggested to hold more generally from various arguments [67, 66], most importantly from DMRG numerics [67]. We provide a proof that additivity does not hold in a large family of IQFT (including the three models we consider), a fact which is explained via the nontrivial relation between state occupation and state density in interacting models. This result is not in disagreement with the DMRG numerics.

The paper is organized as follows: In section 2, we describe the physical situation and give the density matrix that describes the non-equilibrium steady state. We also give an overview of some existing results that we will use. In section 3 we derive the NESSTBA equations. In section 4 we introduce the integrable models that we will consider. In section 5 we present numerical results for the current and the -functions of these models. In section 6, we introduce a family of -functions that can be obtained by normalizing higher order cumulants, and present numerical evidence of their monotonicity. We discuss the physical implications of this property. In section 7 we provide evidence of the Poissonian nature of the energy transfer process. In section 8 we test the additivity property of the current and show that it is non-additive for three different integrable models. We provide a general argument as well as a mathematical proof as to why this must be the case for IQFTs. Finally, we present our conclusions in section 9. Various derivations, proofs, particular limits and consistency checks of our formulation are presented in appendices A-E.

## 2 Physical situation and overview of relevant previous results

### 2.1 Physical description

Consider a homogeneous open quantum chain of length with local interactions. Suppose it is initially broken into two contiguous halves of lengths with commuting hamiltonians and (for instance, some links are broken near some site, so that the two halves do not interact with each other). The two halves are independently thermalized at temperatures3 and respectively (left and right), so that the system is described by the density matrix . At time the two halves are connected in order to form the original homogeneous quantum chain of length , with hamiltonian . Here represents the connection energy between the two halves, which must be local, i.e. involving only few sites around the connection point, and which we assume to be independent of . The system is then allowed to -evolve unitarily until time . One can view this situation as starting with a profile of temperature on the full quantum chain that has a step change around one site, with on the left and on the right of it, and then evolving unitarily from this profile – this point of view requires the notion of a local temperature, which we do not discuss here.

Averages of observables at time are described by

 ⟨O⟩(L;t)=Tr(e−itHLn[ρ0]eitHLO),ρ0=e−βlHLl−βrHLr (1)

where here and below we use . The steady-state limit is the limit followed by of such an average, with the observation lenght kept fixed, see Figure 1 (below we will denote by , and the limits, when it makes sense, of the operators , and respectively). The observation length is the length of the region over which lie both the support of and the contact point .

If this limit exists, then we say that the system reaches a steady state with respect to the observable :

 ⟨O⟩stat:=limt→∞limL→∞⟨O⟩(L;t). (2)

For local observables, we expect the steady-state limit to exist thanks to the locality of the hamiltonians. Physically, this is because we expect the two infinite () halves of the system to play the role, at large enough times, of far thermal reservoirs, able to absorb and emit independent thermalized excitations unboundedly for all times where is a propagation velocity. Indeed for such times, excitations do not have time to bounce off the far endpoints and thus create non-thermal correlations or be re-emitted. Then, observing only locally, the detailed reservoir information is lost and we only see the steady state. In particular, the energy current observable

 J:=i2[HL,HLr−HLl]=i2[δH,HLr−HLl] (3)

is local (and independent of ) thanks to the locality of , wherefore the system should reach a steady state with respect to the energy current.

The first detailed study of this steady-state limit appeared for classical harmonic oscillators in [60]. In quantum systems, similar steady-state limits were discussed and shown to exist in quantum models under certain conditions in [62] and more particularly in the XY model in [61]; in the context of charge transfer in the Kondo model in [65] and in the resonant-level model in [68, 33]; and in the context of energy and charge transfer in conformal models in [52, 35].

As was discussed in many works, the effect of the connection will propagate from the center outwards as time evolves, so that the effective reservoirs, i.e. the regions which are still thermal, will be situated further and further away from the connection point; specifically, the temperature profile will become flat. This holds true in general local models, not just in CFT (although there are difference in details). Hence, if energy transport is purely diffusive, then the steady-state current should be zero (by Fourier’s law). However, if the initial temperatures are different and the energy propagation in the system has a ballistic component, then we expect that the steady-state limit will support an energy flow, : it is then a nontrivial non-equilibrium steady state. It is expected that any integrable system will display a ballistic component to the energy transport: the conserved charges give enough constraints to counteract diffusion [69, 70, 71]. In particular, the conserved charge associated to the energy current is usually one of the first nontrivial “higher” conserved charges in integrable quantum spin chains.

###### Remark 2.1

In CFT and QFT, the momentum is always conserved; hence, as explicitly shown in [52, 35] in CFT, the steady-state limit exists and supports a nonzero energy flow. This means that in the scaling limit (low-temperature, small-gap universal regime) of any quantum chain (with unit dynamical exponent), integrable or not, the steady-state limit should exist (we refer to [35, 66] for a discussion of the interplay between the steady-state limit and the scaling limit). This was indeed observed using DMRG numerics for a particular non-integrable quantum chain in [67]. Of course, in any quantum chain, the nonzero temperatures and, possibly, nonzero gap mean that the system is never exactly in the scaling limit. Hence, if the microscopic system is not integrable, there may be some large, non-universal time scale where non-integrability “kicks-off” and the current starts to decrease. An interesting open question is the nature of this time scale, or by what other mechanism microscopic non-integrability may appear. In the present paper we discuss integrable QFT, which may be thought of as coming from integrable microscopic systems, where there is true ballistic transport.

###### Remark 2.2

We may attempt to make the steady-state limit a bit more precise as follows. We require that the limit exist when each of the conditions and hold almost surely for all , where are propagation velocities of excitations, and is the observation length scale. Here we have in mind a distribution of propagation velocities associated to the available excited states, and the corresponding measure (discrete / continuous parts of the Hilbert space, etc.); the phrase “almost surely” means that the strong inequalities hold for all velocities except possibly for a set of velocities that is of measure zero. Then, the fact that the steady-state limit exists for any local observable can only be true if there is no stationary excitation with nonzero measure: the Hilbert space shouldn’t have a discrete part containing 0-velocity excitations.

### 2.2 Steady state in massive QFT

Assume now that the quantum chain has a parameter (for instance, an external magnetic field) such that there is a quantum critical point at with unit dynamical exponent. Then, at , the energy gap is zero. At and low temperatures, the chain is described by a CFT. Let us now take the scaling limit where and (which tends to zero). This is the low-temperature, small-gap region near the quantum critical point. This limit is universal, and described by massive QFT, with the mass representing the infinitesimal gap . The main result of [52, 66], which we take as our starting point, is a proposed exact description of the non-equilibrium density matrix in the scaling limit, that represents the steady-state limit (2),

 ⟨O⟩stat=Tr(n[ρstat]O). (4)

It is expressed in terms of massive QFT data, and was derived using general QFT arguments.

The non-equilibrium steady-state density matrix is defined as follows. Let us consider a model of relativistic QFT with a spectrum of particle types, with masses . The basis of asymptotic states (here we take them implicitly as asymptotic states), which span the Hilbert space, is described by

 |vac⟩,|θ1,…,θn⟩i1,…,in:θ1>…>θn,i1,…,in∈{1,2,…,ℓ},n=1,2,… (5)

where is the vacuum, and are the rapidities of the asymptotic particles. The density matrix acts diagonally on the basis of asymptotic states. Its eigenvalues are simply given by and

 ρstat|θ1,…,θn⟩i1,…,in=e−∑kWik(θk)|θ1,…,θn⟩i1,…,in (6)

where

 Wi(θ):=micoshθ⋅{βl(θ>0)βr(θ<0). (7)

If we factorize the Hilbert space into a product of the space formed by positive-rapidity particles, and that formed by negative-rapidity particles, then (6) indicates that the steady-state density matrix factorizes accordingly. Specifically, on each factor, it is thermal with inverse temperature , respectively. Note that this density matrix is stationary and homogeneous (translation invariant).

The density matrix (6) agrees with the mathematical description of the energy-flow steady state in the model in [61] and with that obtained in the Ising model from less rigorous arguments in [53]. Its natural counterpart in CFT, where right-movers and left-movers are thermalized at different temperatures, was proven in [52, 35].

###### Remark 2.3

The physical interpretation of the density matrix is rather simple. The steady-state density matrix is the long-time evolution of the initial density matrix . The initial density matrix represents the left- and right-hand sides of the system separately thermalized at inverse temperatures and , respectively, and the time evolution is generated by the full hamiltonian where both halves are connected. In the scattering language, where is the scattering operator. Also, asymptotic states are, by definition, states on the line with the property that when evolved back in time, they become well-separated, well-defined wave packets behaving like free particles (the incoming asymptotic particles). In the scattering language, an asymptotic state has the form where is a free-particle state representing the asymptotically free particles. Hence, the action of on asymptotic states is isomorphic to the action of on the well-separated, well-defined wave packets: if is an eigenvector of , then is an eigenvector of with the same eigenvalue. Since for positive rapidities these wave packets are far on the left, and for negative rapidities they are far on the right, the expression (7) follows. A more complete argument is presented in [66]. More precisely, this argument shows that for any local operator, the contributions of asympotic states with finite numbers of particles to the steady-state average is given by (6) and (7).

### 2.3 Current and fluctuations

In QFT, the energy current operator is simply the momentum density at the position , so that the average current is

 J=⟨p(0)⟩stat.

Besides the current, its fluctuations are also of great interest in non-equilibrium steady states. In quantum mechanics, one has to define fluctuations of the current via appropriate measurement protocols: see for instance [28, 29] for discussions of an indirect measurement protocol, [30, 14, 72, 33, 52] for results within a two-time measurement protocol and connection with the former protocol, and the review [1]. In order to be specific, let us consider the two-time measurement protocol. We perform a first von Neumann measurement of the quantity

 Q:=12(Hr−Hl) (8)

at time 0, returning the value , and then another von Neumann measurement at time , returning the value . We are interested in the statistics of the difference . By standard quantum mechanics, the probability distribution for is

 Ωt(q)=Tr(PQ=qte−iHtPQ=q0n[ρ0]PQ=q0eiHtPQ=qt) (9)

where is the projector onto the eigenspace of eigenvalue of the operator . One may then evaluate the associated scaled cumulants

 Cn:=limt→∞t−1⟨qn⟩cumulantΩt, (10)

where is the cumulant of with respect to . These may be gathered into a generating function as usual,

 F(z):=∞∑n=1Cnzn/n!,F(z)=limt→∞t−1log(⟨eqz⟩Ωt) (11)

(here is a formal parameter). Note that the first cumulant is the average current.

A result of the discussions cited above is that one expects to have the following expression in terms of averages of quantum-mechanical operators:

 F(z)=limt→∞t−1log(⟨ezQ(t)e−zQ⟩stat). (12)

In fact, we expect that in general, the scaled cumulants can also be expressed in terms of connected correlation functions4 of time-evolved current operators,

 Cn=limϵ→0+∫∞−∞du1⋯dun−1⟨J(un−1+(n−1)iϵ)⋯J(u1+iϵ)J(0)⟩connectedstat (13)

where is an imaginary time shift that regularizes the correlators. This expression is derived, under certain assumptions that hold for integrable models, in Appendix A.

Note that an important symmetry relation satisfied by is the fluctuation relation

 F(βl−βr−z)=F(z). (14)

This energy-flow fluctuation relation was derived in various situations in [10, 11, 12, 13, 14] (but see [73] for limitations in the case of finite systems); [1, 2] present good reviews. It was derived recently in [15] within the present physical setup using a scattering formalism.

A recent result of Bernard and Doyon [15] is that in any system where the energy is purely transmitted (that is, where is the system’s scattering operator, see [15]), the scaled cumulant generating function can be evaluated solely from the current at shifted inverse temperatures5,

 dF(z)dz=J(βl−z,βr+z)⇔F(z)=∫z0dz′J(βl−z′,βr+z′) (15)

(recall that ). The upshot is that the scaled cumulants are simply given by derivatives of the current with respect to inverse temperatures,

 Cn+1=dndznJ(βl−z,βr+z)∣∣z=0 (16)

for . Equation (15) was referred to as an extended fluctuation relation (EFR) in [15], as it actualy implies the usual non-equilibrium fluctuation relations (14) if one assumes parity symmetry. As explained in [15], pure energy transmission, and therefore the EFR, is valid in any homogeneous integrable model and in CFT. Note that it implies a relation between the non-equilibrium differential conductivity (where and ) and the noise (second cumulant), which reads .

### 2.4 Previous exact results

Exact results for energy flow and fluctuations in the present physical setup are available in the XY and Ising models [61, 53, 74], and in conformal field theory (critical models) [52, 35, 56]. To our knowledge, the first exact cumulant generating function for energy transfer in quantum systems was found in [51], within a slightly different physical construction of the non-equilibrium steady-state using Caldeira-Leggett baths attached to a chain of harmonic oscillators. Here we just recall the CFT results, as these are most relevant for our purposes.

In general unitary CFT, the current was first calculated in [52] and found to be

 JCFT=cπ12(T2l−T2r). (17)

This has a clear interpretation via the Stefan-Boltzmann law thanks to [75]. The scaled cumulant generating function was first evaluated exactly in [52], giving the expression , in agreement with (15). This means that the scaled cumulants are

 CnCFT=cπn!12(Tn+1l+(−1)nTn+1r). (18)

The function was evaluated in [52] assuming the fluctuation relations (14); it was proven without this assumption in the CFT context using a local operator formalism in [35], and using the Virasoro algebra in [56] (where a more general star-graph configuration is considered). The expression (17) was verified numerically in the XXZ chain by Karrasch, Iland and Moore [67].

## 3 The non-equilibrium steady-state TBA equations

In the present paper, we obtain for the first time the exact energy current and scaled cumulant generating function in interacting (integrable) models.

In integrable models of relativistic QFTs, two crucial properties hold: the scattering of asymptotic particles is purely elastic (the set of momenta is preserved) and the scattering amplitudes always factorize into two-particle scattering processes. Hence, the only dynamical data necessary in order to fully define the model and its “local structure” (its full scattering and its set of local observables) are the two-particle scattering amplitudes. In particular, two-particle scattering amplitudes are solutions of the Yang-Baxter equation, which is nontrivial whenever there are internal degrees of freedom and the two-particle scattering processes is not diagonal in the internal space. This is the essence of factorized scattering theory.

We now consider a general integrable model of relativistic QFT, under the simplifying assumption that the two-particle scattering matrix is diagonal in the internal space (there is no “back-scattering”: in a two-particle scattering process, only phases may occur).

### 3.1 The integral equations

According to (4), the average energy current is obtained by evaluating the average momentum density at an arbitrary point, say , using the density matrix (6):

 J=Tr(n[ρstat]p(0)). (19)

We have provided the exact density matrix (6) as an operator acting on asymptotic states with finitely many particles. In order to evaluate this trace, we further need the matrix elements of . But we need a bit more: we need to describe how to perform the trace operation. This implies understanding in states with nonzero densities of particles (i.e. infinitely many particles). In such finite-density states, the interaction between particles becomes important, and this affects the way the trace is defined. One could say that the information about the “local structure” is encoded both into the matrix elements of and into the way the trace is performed.

In order to add these two elements of information, we follow Al.B. Zamolodchikov’s TBA argument [25]. This is an argument that mixes the ideas of factorized scattering and of the Bethe ansatz. For this purpose, we consider the model on a finite periodic space of circumference . On the finite space, according to [25], the description of integrable models of QFT can be done similarly to that of Bethe ansatz integrable systems. Each state is characterized by its contents in quasi-particles, which possess momenta and energies whose sum give the total momentum and energy of the state, with relativistic dispersion relation. Then, we construct a density matrix that approximates . The natural choice for is defined by its action on every state as follows:

 ρRstat|v⟩=e−∑kWk|v⟩,Wk=ek⋅{βl(pk>0)βr(pk<0). (20)

The finite- approximation of the average current is the average of the local momentum density with respect to :

 JR=TrR[ρRstatp(0)]TrR[ρRstat],limR→∞JR=J. (21)

Note that is translation invariant. By translation invariance, we may replace by . We then have

 JR=R−1TrR[ρRstatPR]TrR[ρRstat] (22)

where is the total momentum.

The current can then be calculated by conveniently introducing a generating parameter associated to the momentum:

 J=−limR→∞R−1ddalogTrR(ρRstate−aPR)∣∣∣a=0. (23)

We may also define the “free energy” and write the current from it:

 fa:=−limR→∞R−1logTrR(ρRstate−aPR),J=ddafa∣∣∣a=0. (24)

The free energy can be evaluated by extending the TBA arguments of Al. B. Zamolodchikov presented in [25] (a similar extension of TBA arguments beyond Gibb’s equilibrium was first done, to our knowledge, in the context of quantum quenches for the Lieb-Liniger model in [76]). This allows us to obtain the following expressions (see Appendix B):

 J(βl,βr) = ℓ∑i=1∫∞−∞dθ2πmicoshθxi(θ)1+eϵi(θ) xi(θ) = misinhθ+ℓ∑j=1∫∞−∞dγ2πφij(θ−γ)xj(γ)1+eϵj(γ) ϵi(θ) = Wi(θ)−ℓ∑j=1∫∞−∞dγ2πφij(θ−γ)log(1+e−ϵj(γ)) (25)

where as usual and is the two-particle scattering matrix, and where is defined in (7). This is the exact non-equilibrium steady-state TBA (NESSTBA) expression for the energy current in general diagonal-scattering integrable relativistic QFT. It turns out that there are many equivalent expressions for the current. For instance (see Appendix B),

 J(βl,βr) = ℓ∑i=1∫∞−∞dθ2πmisinhθμi(θ)1+eϵi(θ) μi(θ) = micoshθ+ℓ∑j=1∫∞−∞dγφij(θ−γ)μj(γ)1+eϵj(γ). (26)

Note that it is customary to define the -functions

 Li(θ)=log(1+e−ϵi(θ)), (27)

as these functions possess interesting and perhaps clearer features than the pseudo-energies . Note also that the pseudo-energies, hence the -functions, have jumps at contrary to the equilibrium case,

 ϵi(+0)−ϵi(−0)=m(βl−βr). (28)

We present various verifications of (25) in Appendices C and D. We verify that the low-temperature expansion of (25) agrees with an explicit evaluation of the trace (22) in the large- limit, using the finite-volume regularization methods of Pozsgay and Takács [77] (and assuming a single particle spectrum for simplicity).Using and generalizing the analysis performed by Al.B. Zamolodchikov in [25], we also verify that the high-temperature expansion of (25) exactly reproduces the CFT prediction (17) of [52, 35], with the correct central charge.

As mentioned, in [15] it was shown that for integrable models of relativistic QFT, the energy-flow steady state satisfies the extended fluctuation relation (15). This relation immediately allows us to calculate the exact scaled cumulant generating function for such models using the NESSTBA expression above. That is, the generating function defined in (11) for the scaled cumulants defined in (10)-(13), is

 F(z) = ℓ∑i=1∫∞−∞dθ2πmicoshθ∫z0dz′xi(θ,z′)1+eϵi(θ,z′) xi(θ,z) = misinhθ+ℓ∑j=1∫∞−∞dγ2πφij(θ−γ)xj(γ,z)1+eϵj(γ,z) ϵi(θ,z) = Wi(θ)+zsgn(θ)micoshθ−ℓ∑j=1∫∞−∞dγ2πφij(θ−γ)log(1+e−ϵj(γ,z)). (29)
###### Remark 3.1

Our verification that the low-temperature expansion of the trace itself agrees with the low-temperature expansion of (25) indicates that Al.B. Zamolodchikov’s TBA arguments indeed can be extended to the non-equilibrium steady-state density matrix (6) without additional subtleties. This is important, because the derivation presented in [66] of the form (6) of the density matrix is, to be correct, only applicable to states with finite numbers of particles, hence to the low-temperature expansion of the trace.

###### Remark 3.2

It is important to note that the ratio of traces on the right-hand side of (22) is not expected to be related in any way to the average total momentum in the physical non-equilibrium steady state. Indeed, the latter average does not exist because the total momentum is not a local operator (hence does not have a steady-state limit). Expression (22) is just a means to obtain a finite- approximation of the actual steady-state current.

###### Remark 3.3

There are many subtleties involved in the derivation of (25) that we presented. Let us discuss some of them.

First, at equilibrium, the procedure for replacing the infinite-space density matrix by a finite- approximation is clear, since even on finite periodic space, there is a clear notion of thermal states. Out of equilibrium, the procedure is less clear: the physical setup we use to construct the non-equilibrium steady state fundamentally requires the infinite-length limit to be taken. So, we must rely on a “non-physical” density matrix , with the requirement that its infinite- limit reproduces in some “topology” (for instance, for every individual matrix element). This is what we have done in (20) for integrable models using Al.B. Zamolodchikov’s ideas. Of course, the finite- “non-physical” density matrix (20) still describes some non-equilibrium steady state in the periodic system: it is a state where there is an imbalance between the populations of particles moving towards the right and of those moving towards the left. Importantly, however, it is a state whose infinite- limit provides the correct population imbalance of asymptotic states for describing the physical non-equilibrium steady state.

Second, the validity of the statement in (21) relies on the fact that thanks to locality, in the infinite- limit, the average of the local operator does not depend on the actual boundary conditions taken (in particular, there are no topological effects that connect the energy density to boundary conditions). An argument can be made as follows. One can consider a different finite- regularization, which by construction correctly gives the physical steady sate. One evolves the initial density matrix for a finite time , and observes the system on a region of length around the contact point, for the speed associated to some small rapidity . The reduced density matrix on that region is an IR regularization, and its limit gives by construction the correct steady state with respect to any local observable. But the point is that and differ from each other only because of two effects: the boundaries are different (in the former, it’s the transient region; in the latter, it’s a periodic boundary condition), and the presence of excitations at rapidities . Therefore, by locality and assuming no topological effects, averages of operators supported in the bulk are equal with respect to both density matrices in the infinite- limit, except possibly for differences of order . Since can be made arbitrarily small for local operators a finite distance away from the contact point (because the boundary at distance is very far at large for any ), then the averages with respect with give the averages in the physical non-equilibrium steady state.

### 3.2 Physical interpretation

Formula (25) has a simple interpretation: the function is the “dressed” momentum at rapidity , the factor is the “bare” (uninteracting) number of levels around , and is the filling fraction. Formula (26) has a similar and perhaps clearer interpretation: the quantity is the true, interacting number of levels around , and the factor is the momentum at . Formula (25) encodes the steady-state density matrix into the driving term of the pseudo-energy, and the locality information of the QFT, including the state density and the momentum matrix elements, into the dressing operations based on the differential scattering phases .

The specialization of (25) to the Ising model, where there is only one particle type and the differential scattering phase is zero, gives a simple formula for the exact energy current,

 JIsing(βl,βr)=m24π∫∞−∞dθsinh2θ1+eW(θ)=12π∫∞mdEE(11+eβlE−11+eβrE). (30)

This can also be obtained in various other ways, and can be extracted for instance from the exact results in the XY and Ising chains in [61, 53].

The second form of in (30) has the structure of the Landauer formula: with where is the energy and is the momentum, we see that we have the number of available channels , times the current through these channels , weighted by the fermionic filling fractions for particles from the right (with temperature ) and from the left (with temperature ). Such a form was also observed in CFT in [52]: there the filling fraction at energy can be taken as being proportional to the Boltzmann occupation . There is however no immediate Landauer form for the non-equilibrium energy current in the general interacting case, because there is no separation between the state densities and filling fractions of right-movers and left-movers: the presence of particles at some energy affect the state density at other energies.

In fact, the correct interpretation of both the Ising and the interacting energy flow should be obtained following what was done in the CFT context in [52]: one must analyze not only the current, but also the fluctuations, and one must look for a formulation where these appear via independent stochastic processes. It turns out that in CFT, the Boltzmann occupation does provide the correct density for independent Poisson processes describing the current and all fluctuations [52]. But, as we will see in section 7, the fermionic or interacting filling fractions are not the correct densities; the correct one will be calculated there.

## 4 The integrable models considered

In the coming sections we will study the non-equilibrium steady-state current and cumulants for several models. We have chosen here three models as good representatives of increasingly complex families of IQFTs: the sinh-Gordon model which may be regarded as the simplest interacting IQFT and therefore provides an ideal starting point for testing our approach; the roaming trajectories model which is still a very simple theory from a particle content and -matrix point of view, but where the dependence of the -matrix on the free parameter leads to a number of new interesting features of the current and its cumulants; finally, the sine-Gordon model at the simplest nontrivial reflectionless point which gives us the opportunity to study a theory with a more complicated particle content but still a diagonal -matrix. It is also a theory which has a well-known discrete counterpart (e.g. the gapped XXZ chain) and therefore results for this model may prove useful when comparing to numerical simulations of the XXZ quantum spin chain.

### 4.1 The sinh-Gordon and roaming trajectories model

The sinh-Gordon model is one of the simplest interacting QFTs and as such it constitutes an ideal testing ground for the ideas above. It has a single particle spectrum and no bound states. The two-particle -matrix of the model is

 S(θ)=tanh12(θ−iπB2)tanh12(θ+iπB2). (31)

The scattering matrix, hence the kernel, depend on the parameter known as the effective coupling constant. This parameter is related to the coupling constant (not related to the inverse temperatures !) in the sinh-Gordon Lagrangian [78, 79]

 L=12(∂μϕ)2−m2β2cosh(βϕ), (32)

where is a mass scale and is the sinh-Gordon field, as

 B(β)=2β28π+β2, (33)

under CFT normalization [20]. The -matrix is obviously invariant under the transformation , a symmetry which is also referred to as weak-strong coupling duality, as it corresponds to in (33). The point is known as the self-dual point. In our numerics for simplicity we will concentrate on the case for which

 φ(θ)=2sechθ. (34)

The sinh-Gordon model is intimately related to the roaming trajectories model [80]: the latter can be understood as a “deformation” of the sinh-Gordon model. Its -matrix is obtained by setting in (31) where . The resulting -matrix famously satisfies all the required properties (e.g. unitary and crossing symmetry) thus describing a new integrable model. The new kernel is then

 φ(θ)=sech(θ−θ0)+sech(θ+θ0). (35)

The presence of the free parameter turns out to have a profound effect on the features of all TBA quantities (e.g. pseudoenergies, -functions, -functions, etc). Indeed the roaming trajectories model’s name refers to the special features of the effective central charge , which is a particular -function, within the equilibrium TBA approach [25] as calculated in [80]. For massive QFTs it is expected that the function , whose definition we recall in Section 6, “flows” from the value zero in the infrared (large ) to a finite value in the ultraviolet (small ). For many theories, including the sinh-Gordon model, the constant value reached as is the effective central charge of the underlying conformal field theory associated to the model [81, 82, 83] (which equals the CFT central charge in unitary models). In this case, that theory is the free massless boson, a conformal field theory with central charge . Therefore, in the sinh-Gordon model, the function flows from the value zero to the value 1 as decreases.

Crucially, when the same function is computed for the roaming trajectories model it shows a very different behaviour. It still flows from the value 0 to the value 1, but it does so by “visiting” infinitely many intermediate values of giving rise to a staircase (or roaming) pattern. The values of that are visited correspond exactly to the central charges of the unitary minimal models of conformal field theory,

 1−6p(p−1),withp=3,4,5… (36)

As we will discuss later, interesting staircase patterns for several non-equilibrium thermodynamic quantities can also be found for this model.

Another observation made in [80] is that the size of the intermediate plateaux that the function develops at the values (36) is determined by the value of . For there is a single plateau at , thus the usual sinh-Gordon behaviour is recovered, whereas the plateaux at (36) become more prominent as is increased. In the limit a single plateau at remains which reflects the fact that the -matrix (31) becomes in this limit, wherefore the model reduces to the Ising field theory. This interesting limiting behaviour was studied in [84] using the form factor approach.

Generalizations of the roaming trajectories model have been constructed based on more complex theories in which case the associated functions have staircase patterns which visit different families of CFTs [85]. Oher families of theories exist where this staircase patterns also arise naturally, albeit involving a finite number of steps which are in one-to-one correspondence with the on-set of unstable particles in the spectrum. These theories are the homogeneous sine-Gordon models, whose -matrix was constructed in [86]. TBA analysis of these models has yielded a variety of staircase patterns [87, 88, 89] which we would also expect to see replicated in computations of the non-equilibrium normalized current and cumulants.

### 4.2 The sine-Gordon model at reflectionless points

We will now consider the sine-Gordon model. This is a more involved theory including several particle types and bound states. In general, the sine-Gordon model has a non-diagonal scattering matrix, which greatly complicates the TBA equations (in non-diagonal models, there is a different formulation that is more efficient [90, 91, 92, 93, 94]). These -matrices depend however on a continuous parameter which for particular values leads to a diagonal theory. Those values correspond to so-called reflectionless points. At reflectionless points, the theory consists of two fundamental particles of mass usually called soliton () and antisoliton () as well as a “tower” of bound states of the former and of each other known as breathers. The number and masses of the breathers depend on the values of the coupling constant. The scattering matrices involving the soliton, antisoliton and first breather () are [95]:

 Sssss(θ)=S¯s¯s¯s¯s(θ)=exp⎛⎜⎝∫∞0dttsinh(1−ν)t2sinhνt2coshtsinhtθiπ⎞⎟⎠, (37)
 Ss¯ss¯s(θ)=S¯ss¯ss(θ)=sinhθνsinhiπ−θνSssss(θ), (38)
 S¯sss¯s(θ)=Ss¯s¯ss(θ)=sinhiπνsinhiπ−θνSssss(θ), (39)
 Ssbsb(θ)=S¯sb¯sb(θ)=sinhθ+isinπ(1+ν)2sinhθ−isinπ(1+ν)2,Sbbbb(θ)=sinhθ+isinπνsinhθ−isinπν. (40)

As we can see, the -matrix vanishes whenever takes integer values; these are the reflectionless points. In those cases the -matrices become in fact those of a well-known diagonal theory, namely the -minimal Toda theory (see e.g. [96]).

Here we will be interested in the simplest reflectionless point, namely when . In this case only the first breather exists and its mass is given by . The scattering amplitudes above simplify greatly and we obtain:

 Sssss(θ)=S¯s¯s¯s¯s(θ)=Ss¯ss¯s(θ)=S¯ss¯ss(θ)=−sinh12(θ+iπ2)sinh12(θ−iπ2). (41)

As anticipated earlier , and

 Ssbsb(θ)=S¯sb¯sb(θ)=sinhθ+i√2sinhθ−i√2,Sbbbb(θ)=sinhθ+isinhθ−i. (42)

The corresponding TBA kernels are

 φss(θ)=φ¯s¯s(θ)=φs¯s(θ)=φ¯ss(θ)=12φbb(θ)=−% sechθ, (43)

and

 φsb(θ)=φ¯sb(θ)=−2√2coshθsech2θ. (44)

The sine-Gordon model describes the scaling limit of a well-known spin chain system, namely the XXZ quantum spin chain in the gapped regime. We hope therefore that some of our results for this model may be comparable to spin chain results.

## 5 L-functions and current in non-equilibrium TBA

In this section we present results for the -functions and current of several models obtained numerically by solving equations (25). We will start with one of the simplest interacting integrable QFTs: the sinh-Gordon model. Results for this model will provide a benchmark for all the general features of the current and -functions. In addition the subtle differences between the equilibrium and non-equilibrium quantities will become apparent when analysing this model. Throughout this and later secions, we will use the variable

 σ=βrβl=TlTr.

### 5.1 The sinh-Gordon model

Since in this model we only have one particle type we will for now drop the particle indices in (25). Employing the kernel (34) in the equations (25) we can now solve for and by applying a numerical recursive algorithm starting with the free solution . Figure 2 gives the -functions for different values of and , where is the mass of the particle.

The first figure with , that is gives the equilibrium quantities which are known from solving the standard TBA equations. It is a common observation from TBA studies that both the pseudoenergies and -functions develop a plateau at high energies in the region when . For most integrable models the values of and at the plateau can be obtained exactly by solving so-called constant TBA equations which are a high-energy limit of the original equations [25, 99].

From the figures above it appears that we start to see this plateau in the red curves. However, the sinh-Gordon model is rather an exception in this respect as in the limit the value of tends to minus infinity, wherefore tends to infinity. This fact has been understood more generally in [97, 98] where is was noted that it is a feature of all affine Toda field theories (of which the sinh-Gordon model is but the simplest example). We will see later that for other models a clear and finite plateau is reached. Comparing the equilibrium case () to the case when thus we observe two main changes:

• The -functions develop a discontinuity at . This is because of the discontinuity in the function due to the different left and right temperatures. This discontinuity is present for all temperatures and the size of the jump is linear (see (28)). Thus when both temperatures are high no discontinuity can be seen. However it is very apparent for relatively large values of as can be seen in the pink curve of the second figure.

• The -functions cease to be even functions of . Again this can be best appreciated for low energies where it is clear that the maximum of the -functions is not located at the origin. Comparing to the equilibrium case, all functions have experienced a shift towards the right. For high energies we have again a plateau of the same height as in the equilibrium case which now extends in the region .

In summary, for we observe a discontinuity at the origin and a shift towards the right of the -functions. Due to the symmetry of the TBA equations a shift towards the left will occur for .

Let us now turn our attention to the functions and . Numerically solving (25) we obtain as shown in figure 3.

As expected we see that the large scale features of all functions are dominated by the function in the -equation of (25). However, because we find that in general . The oddness of the function is only restored in two limiting cases: 1) at equilibrium (); 2) when both left and right temperatures tend to zero ( large). For this reason all functions in the second figure (except the red-dotted one) are almost perfectly odd.

We have seen in the previous section that the sets of equations (25) and (26) are in fact equivalent to each other. Figure 4 provides several examples of the function for the sinh-Gordon model.

The features of the level density are analogous to those of with the difference that the large scale features are now dictated by the term in (26). As a consequence all figures look like deformed -functions where the evenness is broken for and large temperatures. Recalling that is proportional to the density of energy levels, we see that more energy levels are available at positive rapidities if , which agrees with the presence of a nonzero current from the left to the right (more right-moving energy levels available).

Finally, we present the results for the current in figure 5. As expected the current in the massive model and its CFT counterpart expression agree very closely for high temperatures, and both tend to zero for low temperatures and the disagreement also appears small in that region.

### 5.2 The roaming trajectories model

As indicated above, the roaming trajectories model is an interesting deformation of the sinh-Gordon model (from the -matrix point of view) which displays many new interesting features for all the relevant thermodynamic quantities we study here. Figures 6 and 7 show the new structure displayed by -functions and -functions.

The asymmetry of the -functions with respect to the origin is visible in most of the graphs in figure 6 as all curves appear shifted towards the right with respect to the equilibrium solutions presented in [80]. For lower temperatures the discontinuity at the origin becomes apparent once more. The presence of the parameter leads to intricate staircase patters reminiscent of the equilibrium results.

The function displays a very interesting new structure with an alternation of maxima and zeroes at values of which appear to be multiples of . This structure is surprising at first sight as it is rather what we would expect the function to look like in view of figure 6. However, is a derivative w.r.t. the parameter , not , as expressed in the paragraph after (71) in Appendix B. This can be understood as follows. At equilibrium, when the function as defined by the second equation in (25) can be identified with . In other words, if we differentiate the last equation in (25) with respect to at we would exactly obtain the second equation in (25) up to the identification . This ceases to hold when , however the equilibrium features of are still preserved to a large extent with certain deformations. In particular for the -function now dominates the behaviour of suppressing the negative minima of its equilibrium version. A similar deformation occurs for the function where again the out-of-equilibrium features for are dominated by the -function.

Finally, we take a look at the current in figure 8. Comparting to figure 5 we see that the functions seem hardly to have changed. There are however important differences which only become apparent when studying the normalized current (e.g. the current divided by ) and associated cumulants. We will study these in the next section.

### 5.3 The sine-Gordon model at reflectionless points

As explained earlier we will consider here only a very special case of the sine-Gordon model, namely the situation when the whole spectrum consists of soliton, antisolition and one breather and there is no back-scattering. Once more we will start by examining the general behaviour of the -functions for several values of the energy. These are presented in figure 9. In this case the -functions clearly develop a plateau for high energies in the region , which in the first figure corresponds approximately to . The exact height of the plateu is the same as in the equilibrium case which has been studied in [100] by solving the constant TBA equations. The plateaux are located at and which correspond to the solutions and .

The functions and in figure 10 display the same main features already observed for other models. The same applies to the current which for this reason we will not present here.

## 6 Non-equilibrium c-functions

A -function is a function of the energy scale (some physical quantity like the temperature or the distance in a correlation function), or of the position on a renormalization group (RG) trajectory, which represents well the effective number of degrees of freedom at this scale or at that point on the trajectory. In unitary models, it should satisfy 2 properties: it should be constant and equal to the CFT central charge at fixed points of the renormalization group (critical points), and it should monotonically increase with the energy scale. The original -function was defined by A.B. Zamolodchikov in [101] in terms of two-point correlation functions of the stress-tensor, where he used it to prove the inequality between central charges at UV and IR fixed points of an RG trajectory. Al.B. Zamolodchikov defined a new -function based on (equilibrium) TBA [25], sometimes referred to as the effective central charge or scaling function,

 ceff(r):=3rπ2∫∞−∞dθcoshθlog(1+e−ϵ(θ)), (45)

where is the distance scale or inverse temperature scale entering the TBA equations.

Given the behaviour of the current in CFT (17), we expect that the normalized current

 c1(Tl,Tr):=12J(T−1l,T−1r)π(T2l−T2r)=6π2(T2l−T2r)ℓ∑i=