Thermodynamic Bethe ansatz for nonequilibrium steady
[0.1cm] states: exact energy current and fluctuations in integrable QFT
Olalla CastroAlvaredo, 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 largedeviation function) in nonequilibrium 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 nonequilibrium 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 hightemperature limit. We analyze numerically our nonequilibrium steadystate TBA equations for three models: the sinhGordon model, the roaming trajectories model, and the sineGordon model at a particular reflectionless point. Based on the numerics, we conjecture that an infinite family of nonequilibrium 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
Contents
 1 Introduction
 2 Physical situation and overview of relevant previous results
 3 The nonequilibrium steadystate TBA equations
 4 The integrable models considered
 5 functions and current in nonequilibrium TBA
 6 Nonequilibrium functions
 7 The probability distribution and its Poisson process interpretation
 8 (Non)additivity of the current
 9 Conclusion
 A Cumulants in terms of correlators
 B Derivation of the NESSTBA equations
 C Hightemperature limit
 D Lowtemperature expansion
 E Proof of nonadditivity of the current
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 nonequilibrium 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 fluctuationdissipation theorem [9] to systems far from equilibrium (see [1, 2] for extensive discussions). These fluctuation theorems express “universal” properties of certain distribution functions. In nonequilibrium steady states, where constant flows of particles, charge or energy exist, one object of interest is the scaled cumulant generating function (SCGF)
In this manuscript we obtain for the first time the exact nonequilibrium 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 nonequilibrium steady states, concentrating on the study of nonequilibrium 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 noncritical 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 nonperturbatively. 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 freefermion 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 nonequilibrium integrable quantum systems have been obtained in the past, especially for charge flows, in which case the SCGF is referred to as the “fullcounting 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 lowtemperature 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 TBAlike ideas. Exact charge current and shot noise (zerotemperature 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 lowtemperature 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 “stargraph” configurations, see for instance [54, 55, 56]).
There are in fact many ways of theoretically constructing a nonequilibrium steady state (see [57] for a discussion); the aforementioned results involve a variety of setups. The hamiltonianreservoir 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 nonequilibrium density matrix was constructed
We present here a generalization of the TBA approach to nonequilibrium 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 sinhGordon model, the roaming trajectories model, and the sineGordon 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 nonequilibrium analogs of Zamolodchikov’s functions [25], in terms of the current and all higher cumulants. We numerically verify that our nonequilibrium functions are monotonic along RG trajectories, and equal to the central charge at fixed points. Hence, following the usual wisdom, they are (nonequilibrium) 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 largetime 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 nonequilibrium 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 nonadditive 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 AE.
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 temperatures
Averages of observables at time are described by
(1) 
where here and below we use . The steadystate 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 :
(2) 
For local observables, we expect the steadystate 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 nonthermal correlations or be reemitted. Then, observing only locally, the detailed reservoir information is lost and we only see the steady state. In particular, the energy current observable
(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 steadystate limit appeared for classical harmonic oscillators in [60]. In quantum systems, similar steadystate 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 resonantlevel 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 steadystate 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 steadystate limit will support an energy flow, : it is then a nontrivial nonequilibrium 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 steadystate limit exists and supports a nonzero energy flow. This means that in the scaling limit (lowtemperature, smallgap universal regime) of any quantum chain (with unit dynamical exponent), integrable or not, the steadystate limit should exist (we refer to [35, 66] for a discussion of the interplay between the steadystate limit and the scaling limit). This was indeed observed using DMRG numerics for a particular nonintegrable 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, nonuniversal time scale where nonintegrability “kicksoff” and the current starts to decrease. An interesting open question is the nature of this time scale, or by what other mechanism microscopic nonintegrability 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 steadystate 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 steadystate 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 0velocity 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 lowtemperature, smallgap 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 nonequilibrium density matrix in the scaling limit, that represents the steadystate limit (2),
(4) 
It is expressed in terms of massive QFT data, and was derived using general QFT arguments.
The nonequilibrium steadystate 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
(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
(6) 
where
(7) 
If we factorize the Hilbert space into a product of the space formed by positiverapidity particles, and that formed by negativerapidity particles, then (6) indicates that the steadystate 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 energyflow 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 rightmovers and leftmovers are thermalized at different temperatures, was proven in [52, 35].
Remark 2.3
The physical interpretation of the density matrix is rather simple. The steadystate density matrix is the longtime evolution of the initial density matrix . The initial density matrix represents the left and righthand 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 wellseparated, welldefined wave packets behaving like free particles (the incoming asymptotic particles). In the scattering language, an asymptotic state has the form where is a freeparticle state representing the asymptotically free particles. Hence, the action of on asymptotic states is isomorphic to the action of on the wellseparated, welldefined 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 steadystate 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
Besides the current, its fluctuations are also of great interest in nonequilibrium 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 twotime measurement protocol and connection with the former protocol, and the review [1]. In order to be specific, let us consider the twotime measurement protocol. We perform a first von Neumann measurement of the quantity
(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
(9) 
where is the projector onto the eigenspace of eigenvalue of the operator . One may then evaluate the associated scaled cumulants
(10) 
where is the cumulant of with respect to . These may be gathered into a generating function as usual,
(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 quantummechanical operators:
(12) 
In fact, we expect that in general, the scaled cumulants can also be expressed in terms of connected correlation functions
(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
(14) 
This energyflow 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 temperatures
(15) 
(recall that ). The upshot is that the scaled cumulants are simply given by derivatives of the current with respect to inverse temperatures,
(16) 
for . Equation (15) was referred to as an extended fluctuation relation (EFR) in [15], as it actualy implies the usual nonequilibrium 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 nonequilibrium 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 nonequilibrium steadystate using CaldeiraLeggett 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
(17) 
This has a clear interpretation via the StefanBoltzmann 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
(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 stargraph configuration is considered). The expression (17) was verified numerically in the XXZ chain by Karrasch, Iland and Moore [67].
3 The nonequilibrium steadystate 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 twoparticle 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 twoparticle scattering amplitudes. In particular, twoparticle scattering amplitudes are solutions of the YangBaxter equation, which is nontrivial whenever there are internal degrees of freedom and the twoparticle 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 twoparticle scattering matrix is diagonal in the internal space (there is no “backscattering”: in a twoparticle 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):
(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 finitedensity 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 quasiparticles, 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:
(20) 
The finite approximation of the average current is the average of the local momentum density with respect to :
(21) 
Note that is translation invariant. By translation invariance, we may replace by . We then have
(22) 
where is the total momentum.
The current can then be calculated by conveniently introducing a generating parameter associated to the momentum:
(23) 
We may also define the “free energy” and write the current from it:
(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 LiebLiniger model in [76]). This allows us to obtain the following expressions (see Appendix B):
(25) 
where as usual and is the twoparticle scattering matrix, and where is defined in (7). This is the exact nonequilibrium steadystate TBA (NESSTBA) expression for the energy current in general diagonalscattering integrable relativistic QFT. It turns out that there are many equivalent expressions for the current. For instance (see Appendix B),
(26) 
Note that it is customary to define the functions
(27) 
as these functions possess interesting and perhaps clearer features than the pseudoenergies . Note also that the pseudoenergies, hence the functions, have jumps at contrary to the equilibrium case,
(28) 
We present various verifications of (25) in Appendices C and D. We verify that the lowtemperature expansion of (25) agrees with an explicit evaluation of the trace (22) in the large limit, using the finitevolume 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 hightemperature 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 energyflow 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
(29) 
Remark 3.1
Our verification that the lowtemperature expansion of the trace itself agrees with the lowtemperature expansion of (25) indicates that Al.B. Zamolodchikov’s TBA arguments indeed can be extended to the nonequilibrium steadystate 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 lowtemperature expansion of the trace.
Remark 3.2
It is important to note that the ratio of traces on the righthand side of (22) is not expected to be related in any way to the average total momentum in the physical nonequilibrium steady state. Indeed, the latter average does not exist because the total momentum is not a local operator (hence does not have a steadystate limit). Expression (22) is just a means to obtain a finite approximation of the actual steadystate 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 infinitespace 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 nonequilibrium steady state fundamentally requires the infinitelength limit to be taken. So, we must rely on a “nonphysical” 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 “nonphysical” density matrix (20) still describes some nonequilibrium 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 nonequilibrium 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 nonequilibrium 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 steadystate density matrix into the driving term of the pseudoenergy, 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,
(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 nonequilibrium energy current in the general interacting case, because there is no separation between the state densities and filling fractions of rightmovers and leftmovers: 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 nonequilibrium steadystate current and cumulants for several models. We have chosen here three models as good representatives of increasingly complex families of IQFTs: the sinhGordon 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 sineGordon 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 wellknown 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 sinhGordon and roaming trajectories model
The sinhGordon 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 twoparticle matrix of the model is
(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 sinhGordon Lagrangian [78, 79]
(32) 
where is a mass scale and is the sinhGordon field, as
(33) 
under CFT normalization [20]. The matrix is obviously invariant under the transformation , a symmetry which is also referred to as weakstrong coupling duality, as it corresponds to in (33). The point is known as the selfdual point. In our numerics for simplicity we will concentrate on the case for which
(34) 
The sinhGordon model is intimately related to the roaming trajectories model [80]: the latter can be understood as a “deformation” of the sinhGordon 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
(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 sinhGordon 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 sinhGordon 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,
(36) 
As we will discuss later, interesting staircase patterns for several nonequilibrium 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 sinhGordon 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 onetoone correspondence with the onset of unstable particles in the spectrum. These theories are the homogeneous sineGordon 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 nonequilibrium normalized current and cumulants.
4.2 The sineGordon model at reflectionless points
We will now consider the sineGordon model. This is a more involved theory including several particle types and bound states. In general, the sineGordon model has a nondiagonal scattering matrix, which greatly complicates the TBA equations (in nondiagonal 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 socalled 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]:
(37) 
(38) 
(39) 
(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 wellknown 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:
(41) 
As anticipated earlier , and
(42) 
The corresponding TBA kernels are
(43) 
and
(44) 
The sineGordon model describes the scaling limit of a wellknown 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 functions and current in nonequilibrium 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 sinhGordon 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 nonequilibrium quantities will become apparent when analysing this model. Throughout this and later secions, we will use the variable
5.1 The sinhGordon 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 socalled constant TBA equations which are a highenergy 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 sinhGordon 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 sinhGordon 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 reddotted 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 sinhGordon 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 rightmoving 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 sinhGordon 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 outofequilibrium 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 sineGordon model at reflectionless points
As explained earlier we will consider here only a very special case of the sineGordon model, namely the situation when the whole spectrum consists of soliton, antisolition and one breather and there is no backscattering. 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 Nonequilibrium 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 twopoint correlation functions of the stresstensor, 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,
(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