Thermodynamic bootstrap program for integrable QFT's: Form factors and correlation functions at finite energy density

We study the form factors of local operators of integrable QFT's between states with finite energy density. These states arise, for example, at finite temperature, or from a generalized Gibbs ensemble. We generalize Smirnov's form factor axioms, formulating them for a set of particle/hole excitations on top of the thermodynamic background, instead of the vacuum. We show that exact form factors can be found as minimal solutions of these new axioms. The thermodynamic form factors can be used to construct correlation functions on thermodynamic states. The expression found for the two-point function is similar to the conjectured LeClair-Mussardo formula, but using the new form factors dressed by the thermodynamic background, and with all singularities properly regularized. We study the different infrared asymptotics of the thermal two-point function, and show there generally exist two different regimes, manifesting massive exponential decay, or effectively gapless behavior at long distances, respectively. As an example, we compute the few-excitations form factors of vertex operators for the sinh-Gordon model.


Introduction
One of the ultimate goals in the study of quantum many-body systems and quantum field theory (QFT) is the computation of correlation functions of local (physical) operators. Correlation functions are directly related to observable quantities such as scattering amplitudes in high energy physics or response functions in condensed matter settings. Correlation functions are also the key quantity to understand the emergence of macroscopic behavior from the underlying quantum microscopic description.
One long-standing challenge has been the determination of correlation functions in strongly correlated quantum field theories at finite energy densities. Such situations arise naturally in condensed matter settings, where systems at finite temperataure or finite chemical potential (e.g. biased systems) are commonly studied [1]. Similar problems are also of interest in high energy physics and cosmology. Some prominent examples are the study of hot and dense quark matter inside neutron stars [2], and the process of quark-gluon plasma formation in heavy ion collisions [3,4].
Recent years have also seen a renewed interest in the non-equilibrium physics, especially in low dimensions, of QFT's and other quantum many-body systems [5,6,7]. Driven by questions of thermalization in quantum physics, this line of research led to a wider realization that, at least in principle, it is feasible to realize many exotic states of matter. They are exotic because they are very far from the vacuum, and cannot be described by traditional equilibrium thermodynamics. Examples of such states are those characterized by the generalized Gibbs ensembles [8] (the precise structure of the GGE in QFT's has been studied in [9,10,7,11]), inhomogeneous quasi-stationary steady states [12,13], or even states capturing parts of the dynamics (Floquet dynamics) in periodically driven systems [14,15,16].
The natural question then arises: how can we compute correlation functions in such far from the vacuum, highly-excited states, characterized by a finite energy density. Especially, that the standard QFT techniques were mostly developed to tackle computation of the vacuum correlation functions.
In this work we address this problem for a specific subset of QFT's: integrable QFT's in 1+1 dimensions. To achieve the aim of computing the correlation functions we develop a generalization of the form factor bootstrap approach [17]. In this approach computation of the correlation function is divided in two steps. In the first step one has to determine form factors -matrix elements of the operator under the consideration. The second step is the evaluation of the spectral sum whose main ingredients are the form factors. The main result of our work is the generalization of the form-factor bootstrap program for states with finite energy density. This bootstrap program allows us to determine the form factors and realize the first step of the outlined above procedure. We then present an expression for the two-point correlation functions as an expansion in terms of these new form factors. The actual explicit evaluation of leading terms of the correlation functions is still a demanding analytic problem, as the form factors have singularities that need to be handled with care, however, we present an unambiguous prescription for regularizing these singularities.

Integrable QFT's, bootstrap program and correlation functions
Integrable field theories (IQFT) in 2d form an interesting subset of QFT's. They are characterized by elastic and factorizable scattering processes. Elasticity means that the set of particles momenta in a given state is fully conserved after any scattering process, and factorizability means in this context that any multiparticle scattering process can be written as a product over 2-particles scattering matrices. In turn, these features are related to the existance of an infinite family of local conserved charges. This rich structure allows often for an exact determiantion of the S-matrix [18], just from a few physical principles.
Besides the intrinsic mathematical interest in IQFT's as exactly solvable toy models, they have also found numerous applications in a wide variety of physical problems. To give a small and not exhaustive list of examples, IQFT's can be used to describe the low energy dynamics of different quantum spin chains [19] and systems of cold atomic gases [20] which are relevant to condensed matter experiments. Integrable QFT's can also be viewed as specific perturbation of 2d conformal field theories [21], and provide a nonperturbative approach to study near-critical systems. Applications in high energy physics include a recent proposal, showing that the worldsheet dynamics of the confing string in quantum chromodynamics can be approximately described by a massless 2d IQFT [22], as well as numerous applications of integrability in the study of planar N = 4 supersymmetric Yang-Mills theory and the Ads/CFT correspondence [23].
One of the cornerstones in the development of the IQFT was construction of the form-factor bootstrap program [24,17]. Following from a few physical principles, the bootstrap program allows to significantly constrain the form factors of local operators between states involving a finite number of particls on top of the vaccuum. In many circumstances the constraints are so strong that they allow the complete determination of the form factors, without relying on a Lagrangian formulation, or perturbation theory. Once these form factors are known, it is easy to write expressions for correlation functions as a form factor expansion. As a result, the vacuum correlation functions of integrable QFT's are, at least in principle, exactly computable. In practice, for models with rich spectrum and complicated scattering determinantion of the vaccumm correlation function might still be technically difficult.
There have been several previous approaches to computing correlation functions of integrable QFT's in finite-energy-density states. In particular, in Ref. [25], Leclair and Mussardo conjectured a general expression for one-point functions of local operators at finite temperature, in terms of a sum of standard, zero-temperature form factors. A similar expression in terms of zero-temperature form factors was then also conjectured in [25], attempting to compute the two-point function at finite temperature.
Later independent checks have proven more rigorously that the Leclair-Mussardo conjecture is indeed correct for computing one-point functions [26] (though some of concepts required for this proof were discovered as far back as in [27]). The validity of the two-point function proposal, however, has been much more contested, with some counter examples first explored in [28] and [29]. A more fundamental problem with the two-point Leclair-Mussardo proposal was first addressed in [30]. Namely, the form factor expansion proposed in [25] contains some divergent contributions, arising from kinematical poles of the form factors, and no clear unambiguous proposal is given on how to regularize this expression. It has been suggested that the Leclair-Mussardo two-point formula is only valid for cases where the kinematic poles do not present a problem, and the thermal background has no effect of dressing the particles of the theory, the only known examples being free theories [31,29] and field theories in 'tHooft's planar large-N limit [32].
While the Leclair-Mussardo proposal for the two-point function, which in principle is an expression applicable to all values of the temperature, does not seem to be valid, some attempts have been successful at developing low temperature expansions [33,30]. These approaches consist simply on expanding the Gibbsensemble expressions at low temperatures, and carefully handling the kinematical singularities order by order in this expansion. These expansions naturally are expressed in terms of standard zero-temperature form factors, since the low-temperature expansion is centered around the zero-temperature ground state. These approaches do not seem to provide any viable pathway towards computing correlators at higher temperature values, except if one is somehow able to resum an infinite number of terms.
A new approach towards two-point functions in highly excited states was recently developed in [34], by showing a new way to regularize the kinematical singularities of form factors at all orders in the lowtemperature expansion. This approach consisted on treating the two operators in the two-point function as a single quasilocal operator, whose form factors have similar kinematical pole structure as those of local operators. One then can rely on the results of the one-point Leclair-Mussardo formula, which provide the prescription to regularize kinematical poles of the quasilocal operator. While this approach is conceptually groundbreaking, solving the puzzle of how the kinematic poles can be regularized to all orders, pragmatically, it can still only be used as a low-temperature expansion centered around the zero-temperature ground state, unless one can find another way to resum higher order terms.
We finally point out that there has been recent progress in the computation of thermodynamic correlation functions numerically through the truncated conformal spectrum approach (TCSA) [35]. This approach consists on treating the interacting QFT as a perturbation of a conformal field theory, whose finite-volume spectrum, and matrix elements are known. At finite volume, the spectrum becomes discrete, an can be truncated at some high-energy cutoff. The interacting Hamiltonian is then a finite, discrete matrix that can be diagonalized numerically. Only recently [36], the TCSA has been applied towards the computation of correlation functions in highly excited states, for the sine-Gordon model. We will concentrate only on simpler QFT's with diagonal scattering, such that the computation of sine-Gordon correlators is beyond the scope of this paper, so we will not be able, at present, to make a quantitative comparison with this new numerical approach.
The approach we will develop in this paper has the novel feature that it is not a low-temperature expansion, but the expansion can be centered around any given finite energy density state.

Summary of the results
The starting point of our approach is the idea that thermal expectation values of local operators (also expectation values in generalized Gibbs ensembles) are equivalent in the thermodynamic limit (where the system size, L, goes to infinity) to the expectation value of the operator in just one single representative eigenstate. The representative eigenstate, which we denote as |ϑ is chosen such that the particle occupation number is the same when evaluated in the ensemble average, or in the single eigenstate, where A † (θ), A(θ) are creation and annihilation operators, respectively, of a particle with rapidity θ. The idea of the representative states orginated during the early studies of integrable models [37]. Since then it was applied succesfully in numerous situations. These include studies of thermal correlation functions [38,39] and also studies of non-equilibrium QFT's, in an approach known as the "quench action" method [40,41,42,43].
To specify the representative state one then simply needs to specify the distribution ϑ(θ).
Integrable field theories are characterized by a stable particle content (due to the property of elastic scattering). The multi-particle eigenstates in such theories are constructed simply by specifying rapidity θ of each particle. States of finite energy density are states in which there is an extensive number of particles such that the total energy of the state scales with the size of the system where H is the QFT Hamiltonian. For simplicity of presentation, for the majority of this paper we assume that we are working with an IQFT with only one species of particle in the spectrum, though it is also possible to generalize our results for other models with a richer spectrum. To each state |ϑ we associate a smooth function ϑ(θ) that plays a role of the filling function. The factor ϑ(θ)dθ determines how densely the interval [θ, θ + dθ] is filled with particles. We can then consider a general two-point correlation function where for simplicity of the presentation we take both operators to be the same, which doesn't need to be the case in general. One way of handling the computation of such two point functions is by utilizing the Lehmann (spectral) representation. To this end, a crucial step is to assume that we can construct all the eigenstates |ϑ such that lim In other words we assume that we can span a subspace of the Hilbert space relevant for the computation of a given correlation function. Given that, we have The matrix elements ϑ |O(0)|ϑ , the form factors, are building blocks of the correlation function in the spectral representation. The aim of our work is to present an axiomatic way of computing such form factors in 2d integrable field theories, in the case where the states, |ϑ and |ϑ are states with finite energy density. The form factors ϑ|O(0)|ϑ are just matrix elements between two finite energy density states. For local physical operators, the form factor between two arbitrary states, |ϑ and |ϑ , is in most cases equal to zero because these two states will, in general, differ macroscopically: there will be an infinite number of particles that are present in one of the states and not in the other. Local physical operators cannot connect such macroscopically distant states. Therefore, the only kind of non-zero matrix element are those where the state |ϑ is not too far from |ϑ , and involves only microscopic deviations from this state. Therefore, |ϑ must be given by the same filling function ϑ(θ) plus some countable number of modifications: like adding or removing a particle with certain rapidity θ. We write |ϑ → |ϑ, θ 1 , . . . , θ n .
Note that the removal of a particle from the state |ϑ can be considered as the insertion of a "hole" excitation. Given the relativistic invariance of IQFT's, a hole-excitation in the incoming state, is equivalent by crossing symmetry to a particle in the outgoing state. The crossing of a particle from incoming to outgoing states amounts to the transformation of its rapidity, θ → θ + πi, such that the notation (6) includes the insertion of both, particles and holes, as long as we allow πi shifts in rapidity.
With this notation the form factors of local operators are defined as The filling function ϑ(θ) plays a role of a background over which we construct additional excitations. We call such form factors "finite background form factors". Sending the background to zero means considering form factors and excitations over the vacuum.
To anticipate some of the results we will present now the most important axioms for the finite background form factors. Additional axioms are to be found in the text. The main axioms are: • the scattering axiom • the periodicity axiom • the annihilation-pole axiom These three axioms give information about how the form factors behave under the exchange of two particles in the state, under the application of crossing symmetry, and under the annihilation of particle-hole excitation, respectively. The axioms depend on two functions: the scattering matrix S(θ 1 , θ 2 ) and R ϑ (θ|θ 1 , . . . , θ n ) which describes the collective effect of the background particles on the scattering of the particle with rapidity θ due to the presence of the additional particles θ 1 , . . . , θ n . This function describes the back-reaction of the background upon inserting extra particles. In the low density limit the R ϑ function tends to 1 and the axioms turn into the standard vacuum form-factors axioms. We finish the introduction with a few important remarks. First, we point out that the expansion (5) in terms of form factors, is not a low-temperature (or lowenergy-density) expansion. The expansion is centered around the representative state, |ϑ , and the leading contributions to (5) will come from the cases where the states |ϑ are small deviations around |ϑ . The state |ϑ is far from the vacuum, and there is no restriction for the energy density to be low.
Second, an approach for the computation of correlation functions similar to the one we are describing here has already proven to be successful in the study of correlation functions of the non-relativistic Lieb Liniger model [44,45,46]. In this model, form factors with an arbitrary number of particle excitations are known from using algebraic Bethe ansatz techniques [37]. The form factors of the form ϑ |O(0)|ϑ can then be computed by taking the thermodynamic limit of the many-particle form factors. Considering QFT's instead of the Lieb Liniger model will give us a set of additional analytic tools, since relativistic invariance places very strong restrictions on such form factors, and the form factors will be computed by finding functions that satisfy all these restrictions. Our approach will then be very different than that used for the Lieb Liniger model in [44]. Our approach will not be to take the thermodynamic limit of a previously known many-particle form factor (which are extremely difficult to calculate exactly in QFT). Instead, we will derive the finite-energy-density form factors by finding functions which solve the large number of constraints which arise from integrability and relativistic invariance in the presence of the finite-energy-density background, without ever having to compute complicated many-particle form factors.
Third, we point out that the idea of computing form factors in the presence of a finite-energy density background in general is not new. In particular in [47], an approach was developed for Ising field theory for computing form factors with excitations on top of the thermal Gibbs ensemble. The authors dervied a set of axioms for local operators (for which the new axioms are simple in this non-interacting theory), as well as non-trivial semilocal operators (which we will not discuss in this paper). It was however, not easy at all to generalize this approach to interacting field theories. The reason this approach is so much more difficult than ours, is that in [47], it was attempted to compute form factors including excitations on top of the genuine Gibbs ensemble, which involves a sum over an infinite number of states. Our approach is made much simpler by using the concept of the representative state (1). It is conceptually much simpler to understand what it means to have an excitation over the representative state, while adding an excitation on top of an averaging ensemble is a much more elusive concept, making it much more difficult to handle the form factors and derive axioms for them.
The organization of the manuscript is the following. In Sections 2 and 3 we introduce the neccessary background on the vacuum form factor bootrap program and on integrable QFT's at finite volume. Some standard computations of the thermodynamic limit of the finite volume QFT are recalled in Appendix A. Those two sections (and Appendix) serve as a basis for Section 4, where we formulate the axioms satisfied by the finite-background form factors. The logic behind the derivation of these axioms is the following. We start with knowledge of the vacuum form-factor axioms in the infinite volume. We then place the system in a finite volume, and use previously known relations between the finite-and infinite-volume form factors [48]. We use these relations to derive constraints satisfied by finite-volume form factors, with a large, but finite number of particles. Then we go back to the infinite volume limit but this time we keep the density of particles fixed. In taking this thermodynamic limit, the relations we found for the finite-volume form factors yield the new axioms for the finite-background form factors, now at infinite volume.
After deriving the axioms, in Section 5 we analyze the consequences they impose on the form factors. Specifically we introduce a notion, familiar from the vacuum form factors axioms, of minimal form-factors. In Section 6 we conjecture a formula for the correlation function which resembles the well-known LeClair-Mussardo formula. The main difference lies in the form factors used, which in our case are now dressed by the finite background. Our two-point function is also already regularized, since the problems with kinematical poles are already absorbed into the definition of the finite-background form factors. In the Section 7 we show that there are generally two distinct infrared limits of the correlation functions, which show either exponential or polynomial decay at long distances. The closing section 8 is devoted to Sinh-Gordon theory, a simple example of an integrable field theory. We show there, in a concrete example, how the form factor axioms can be used to compute the one-and two-particle form factors of vertex operators in the presence of a thermal background.

The form factor bootstrap program
We first review the standard bootstrap approach to form factors in integrable QFT. We assume diagonal scattering for simplicity, and consider theories with one species of particle. We parametrize the energy and momentum of a particle of mass m, as E = m cosh θ, p = m sinh θ, respectively, where θ is the particle's rapidity.
We define the n-particle form factor of some local operator, O(x), as ordered as θ 1 > θ 2 > · · · > θ n . We can also consider matrix elements with particles in the outgoing state, using the crossing relation [17] f O (θ |θ 1 , . . . , θ n ) ≡ θ |O(0)|θ 1 , . . . , θ n where S(θ i − θ j ) is the two-particle S-matrix, which depends only on the difference of the two rapidities, as a consequence of Lorentz invariance. The last term in the right hand side of (11) is the connected piece of the form factor. The form factors satisfy a set of very restricting axioms, which are a consequence of elasticity and factorization in the integrable theory, as well as other physical considerations, such as Lorentz invariance, and unitarity. An exact expression for the form factors can often be computed by finding a solution to the set of axioms. We now list all these restrictions.

Lorentz symmetry
If O(x) is an operator with Lorentz spin, s, then under a boost, where all the rapidities are shifted as θ i → θ i + α, the form factor transforms as In particular, if O(x) is a scalar operator, with s = 0, then it is invariant under a Lorentz boost. The Lorentz symmetry implies, that besides the s-dependent pre-factor, the form factors depend only on the differences of rapidities θ ij = θ i − θ j . Therefore the form factor depends independently only on n − 1 of the n rapidities.

Scattering axiom
The order of two particles in the incoming state may be exchanged by scattering the two particles, as where S(θ i − θ j ) is the two-particle S-matrix. The S-matrix can typically be derived from a set of axioms which follow from integrability, we refer the reader to [18,49] for more information. Here, we assume the S-matrix of the IQFT is known, and use it as an input in the form factor bootstrap program.

Periodicity axiom
Following the generalized crossing relation, (11), a particle in the incoming state may be turned into an outgoing particle by shifting its rapidity by πi. By applying crossing symmetry twice, we find that we can turn the last particle in the incoming state, into the first one, as This axiom may be modified when the particle θ n and the operator O(x) are not local with respect to each other (for example when the particle is a topological soliton), but for simplicity, we will not consider such cases for now. The three axioms presented so far are enough to compute the so-called "minimal form factors", which are independent of the particular operator, O. These minimal form factors are defined as the maximally analytic solutions of the scattering and periodicity axioms, i.e. satisfying the so-called "monodromy properties", with no poles or zeroes in the "physical region", θ ij ∈ [0, 2πi].
The minimal solution of (14) can be found [24] by expressing the S-matrix as Then the minimal form factor is The general form factors can be expessed as where K O is necessarily a periodic function in all the rapidities, for any i = 1, . . . , n. These functions K O must contain all necessary poles in the physical region, and any dependence on the particular operator, O(x). These function can be expressed generally as where Q O (θ 1 , . . . , θ n ) and D(θ 1 , . . . , θ n ) can be chosen to be polynomials involving powers of cosh(θ ij ). The denominator, D(θ 1 , . . . , θ n ) is chosen as the minimal function that includes all the necessary poles in the physical region, and does not depend on the operator, O(x). All the operator dependence is then incorporated in the numerator function, Q O (θ 1 , . . . , θ n ).

Annihilation-pole axiom
In a relativistic QFT, it is possible for a particle in the incoming state to annihilate with its antiparticle on the outgoing state (in the simplest cases we consider, the particle is its own antiparticle). As a consequence, the form factor has a kinematical annihilation pole, for example for the particles with rapidities θ 1 and θ 2 , at the point θ 1 = θ 2 + πi. The annihilation-pole axiom relates the residue of the n-particle form factor at the annihilation pole, to the n−2 particle form factor, where the two particles have annihilated, explicitly, This axiom then provides a recursive relation, linking the n and n − 2 particle form factors. This recursion can often be used to normalize higher-particle form factors, based on the normalization of lower form factors.

Cluster properties
If one boosts the rapidities of only a subset of the particles in the incoming state, the form factors are expected to satisfy the so-called "cluster properties" [50], where N O 1 cluster is some normalization constant, and O 2,3 are some other local operators. In particular, this relation is expected to simplify in simple theories with no additional internal symmetries [51], where some of the form factors should "self cluster", such that

Bounds on growth
One can further restrict the function Q O (θ 1 , . . . , θ n ) by placing bounds on how fast this function can grow as we take |θ i | → ∞. An upper bound on the growth can be derived (for unitary theories) by considering the two-point function of the operator O, which at very short distances can be described by an underlying CFT, using an operator product expansion. The derivation can be found in [52,53], and here we only cite the result, with y O ≤ ∆ O , where ∆ O is the conformal weight of the operator. This constraint (17) places an upper bound on the degree of the polynomial, Q O (θ 1 , . . . , θ n ), and this bound is specific to the particular operator.

Translation
Once the form factor of the operator at the origin, x = 0, is known, it is easy to obtain the form factor for the operator at any space-time point, x. Using Lorentz symmetry, one can write where P n = n i=1 m sinh θ, and E n = n i=1 m cosh θ are the total momentum and energy, respectively, of the n-particle state.
In addition to these axioms, for theories whose spectrum includes particle bound states, there exist additional axioms, concerning the additional bound-state pole structure. If the the interactions between particles are such that they may form bound states, the form factors will have additional poles in the physical region, corresponding to the fusion angle between the two elementary particles. It is also possible for the form factor to have higher-order poles corresponding to processess with bound states as intermediate states [17]. For simplicity we will not consider theories with bound states here. For simplicity, we consider only QFT's with no bound states, and only one particle species, so we will not list any bound-state axioms, or their modification at finite background.
We note also that we are only considering operators, O(x), which are local, relative the particle excitations. This means we are excluding the computation of semi-local operators, which can also be studied in the bootstrap approach by modifying some of the form factors with some semi-locality factor. Semi-local operators commonly arise, for example, in cases where the particle excitations are topological kinks (such as in the Ising field theory, or the sine-Gordon model), and are therefore not purely local excitations of the fundamental field. The derivation of finite-background form factor axioms for semi-local operators is a much harder problem, since there are fundamental problems with the definition of semi-local operators at finite volume (which we need to formulate as an intermediate state before taking the thermodynamic limit), which induce non-trivial effects, such as operator mixing. These difficulties were encountered and discussed for the non-interacting Ising field theory in [47], but they are beyond the scope of this paper.

Integrable QFT in a finite volume
Now we want to consider form factors of local operators with incoming and outgoing states with a finite energy density background, and finite excitations on top of this background. This background can be defined by first placing the system at finite volume, L, and then taking the thermodynamic limit, L, n → ∞, with the particle density, n/L, fixed. We then first need to understand how to work with form factors at finite volume. Most of the background knowledge introduced in this section can be found in more detail in the original references [54] (the first use of the thermodynamic Bethe ansatz for QFT), and [48] (derivation of expressions for finite-volume form factors).
We place system at finite volume and assume periodic (anti-periodic) boundary conditions for bosons (fermions). This introduces a quantization of allowed particle rapidities. We shall be mainly concerned with the integrable theories of the fermionic type 1 in which The quantization conditions are given by the Bethe equations, It is convenient to consider the logarithm of this expression, to find the allowed rapidities for an n-particle state satisfy where δ 0 (θ) ≡ −i log(−S(θ)), I k is a quantum number, and there is one such equation for each k = 1, . . . , n.
For fermionic (antiperiodic) boundary conditions, the quantum numbers are integers, while for bosonic (periodic) boundary conditions they are integers when n is odd and half-odd integers for n even. The tildes on the rapidity variables denote that these are finite-volume solutions of the quantization equations (20). A state is then fixed by choosing a set of integers, {I k }, which are an input for finding the solutions {θ k } of (20). The total energy of the n-particle state is given by At finite volume, it is then convenient to switch to an "integer" basis of states, such that where ρ(θ 1 , . . . ,θ n ) is the determinant of the Jacobian of the transformation between the two bases, We point out that these determinants, ρ({θ}) are periodic under the transformationsθ i →θ i + 2πi. They are also symmetric functions ofθ j . The general form factors with particles in the incoming and outgoing states, without any disconnected pieces is given by [48] The expression (21) is exact at finite volume to all orders of powers of 1/L corrections, and the only corrections are exponentially suppressed. These corrections, however interesting for the structure of the finite volume IQFT [55], do not play any role in the thermodynamic limit. It is important to note that the only way the form factor (21) can have any disconnected contributions, is if the two sets of rapidities in the incoming and outgoing states are exactly equal to each other, or {I } = {I}, this is because changing one of the integers, will produce an O(1/L) shift in all of the other rapidities,θ i , which means the rapidities in the incoming and outgoing states will no longer match. The thermodynamic limit is taken by increasing L while leaving the energy density constant. The total energy of the finite-density state in the thermodynamic limit can be written as [54] where ρ p (θ) is a distribution describing the density of particles. We can also define the density of available states, ρ s (θ), with which we define the particle occupation number ϑ(θ) = ρ p (θ)/ρ s (θ). The quantization conditions in the thermodynamic limit become an integral equation for the densites (c.f. Appendix A) so that only one of the distributions ρ p (θ), ρ s (θ), ϑ(θ) is independent and contains all macroscopic information. We choose ϑ(θ) to play this role. In the formula above we have defined ϕ(θ) = ∂δ(θ)/∂θ. We also note that function ρ s (θ) is 2πi periodic. It will later be convenient to also define a density of holes, ρ h (θ), or density of available states which are not occupied by particles, such that ρ s (θ) = ρ p (θ) + ρ h (θ).
Acquired now with the tools to describe the thermodynamic limit we consider this limit for a number of quantities. We define, in the thermodynamic limit, a state corresponding to some macroscopic distribution, ϑ(θ), as where one sums over all the microscopic states with particle content {θ}, leading to the same macroscopic distribution, ϑ. S[ϑ] is the Yang-Yang entropy associated with this state, The norm of the state |ϑ is given that the finite-volume states are normalized as {I}|{I} = 1. The thermodynamic limit on the right hand side is not well-defined because functions ρ({θ}) contain some irregular terms [56]. If we consider a normalized form factor, the irregular terms cancel with corresponding terms in the form factor. Therefore, only the normalized form factor has a well-defined thermodynamic limit. For this reason the right hand side of this formula should be understood as an intermediate expression. We will formulate the axioms for normalized form factors and therefore we shall not be concerned with giving precise description of this thermodynamic limit. Introducing the additional particles induces a change in the rapidities of the background. This backreaction of the background particlesθ is captured by the back-flow function F (θ|{θ i }). We assume that the back-flow function is separable. The full back-flow can be written as a sum of separate contributions coming from each particle {θ i }, The shift in the background rapidities has two origins. First, it is the effect of interactions between the particles. Second, in bosonic theory with periodic boundary conditions, there is an additional shift. This additional shift occurs to keep the boundary conditions periodic. When adding or removing a particle, the quantum numbers must be shifted by 1/2. This results in an extra shift of rapidities. We shall call the back-flow function which takes into account only the shift due to the interactions, the fermionic back-flow function and denote it F F (θ|{θ i }). The bosonic back-flow function F B (θ|{θ i }) takes into account both effects. The single particle fermionic back-flow is given by (c.f. Appendix A) We find the hole back-flow to be minus back-flow of the particle, We will introduce the bosonic back-flow function later. Let us now turn to the problem of defining a finite-volume representation of a state with excitations over the background.
We denote the state specified by a particle filling function ϑ(θ) and a set of additional excitation {θ i } by |ϑ, θ 1 , . . . , θ n . This state is understood as thermodynamic limit of the finite-volume state where again one sums over all the microscopic states leading to the same macroscopic state. The prime in {θ } indicates that the rapidities in the background state are modified by the back-flow function. At this stage it is convenient to split the quantum numbers in two sets: {I} is a set of integers which leads to the distribution ϑ(θ) in the thermodynamic limit, and I 1 , . . . , I n are the integers corresponding to the additional excitations. The choice of {I} is not unique, indeed there are many choices of a quantum number that in the thermodynamic limit yield the same excitation θ, which is why we average by summing over all possible such sets {I}.
We are interested in form factors of the form ϑ|O|ϑ, θ 1 , . . . , θ n . We relate them to the thermodynamic limit of the finite system form factors in the following way. From the definition (23) and (28) it follows that where function δS[ϑ, {θ i }] is the differential entropy [57] δS[ϑ, with s[ϑ, θ ] defined in (24). The differential entropy is antiperiodic in θ 1 , To evaluate the sum in (29) we follow the prescription introduced in [58]. The main point is that the form factors {θ}|O|{θ },θ 1 , . . . ,θ n in the thermodynamic limit are almost diagonal in sets {θ} and {θ }. Therefore, choice of {θ} fixes {θ } up to small (of order 1/L) displacements. These displacements correspond to tiny non-dispersive excitations, referred to as soft modes. The summation over the soft modes leads to a renormalized, dressed, form factor The dressed form factor is, by the definition, independent of the discretization {θ} and simply depends on the choice of ϑ. Therefore At the present state of the development of the Bethe Ansatz techniques the evaluation of the dressed form factors is a difficult task. However certain progress in this direction was achieved in [56,59]. Here, following the logic presented in [44] and further developed in [45,46], we approximate the sum over soft modes by simply picking one specific, reference state. The choice of this state is operator dependent and it also might depend on the excitationsθ 1 , . . . ,θ n . Therefore we replace the dressed form factor with a single form factor times the phase space of soft modes, which is exp(δS[ϑ]), The thermodynamic form factor then becomes At this stage it is worth noting that the axioms we will derive will not really rely on approximating the dressed form factor by a single reference state form factor. Instead the assumption is that the new factors coming from exchaning particles or from periodicity axiom are the same for the dressed form factor and for the form factor of the reference state. Finally, we write form factors as functions of quantum numbers The formula, relating the finite volume form factors with their thermodynamic limit, is the main result of this section. In the next section we will use it to conjecture the axioms for the finite background form factors from the axioms for the vacuum ones. The derivation presented here follows closely conceptually similar derivations in non-relativistic integrable models shown for example in [44,58].
Note that we assume we know, as an input, the one-point expectation value, ϑ|O|ϑ / ϑ|ϑ , which we can think of as a zero-particle form factor in our approach. In general, this value can be computed from the one-point Leclair-Mussardo formula [25], which is an expansion in terms of standard zero-background form factors. For some specific operators, more convenient closed-form expressions may be found, which essentially resum the Leclair-Mussardo series. One example is when the operator is given by the trace of the energy momentum tensor, O = T µ µ , in which case an exact thermal expectation value can be computed directly from the thermodynamic Bethe ansatz formalism [54]. Another more recent example is the conjectured Negro-Smirnov formula [60], which gives a closed-form expression for the expectation value of vertex operators in sinh-Gordon on a given thermodynamic background.
The analysis presented so far works for fermions with antiperiodic boundary conditions. However, for bosons with periodic boundary conditions there is a small caveat. The problem is that in such theories quantum numbers change from integers to half-odd integers when we add or remove a particle. Consider a form factor with a single particle in a bosonic theory ϑ|O|ϑ, θ 1 .
Following the procedure outlined above we can discretize state |ϑ to find a set of integers {I} defining a finite system reference state |{I} . However what quantum numbers {Ĩ} should we choose to represent background of |ϑ, θ 1 ? They cannot be the same quantum numbers as {I} because there is now one more particle. Therefore {Ĩ} must differ from {I} by ±1/2. In principle any choice of ±1/2 is valid. However, guided by the physical intuition, we choose the quantum numbers {Ĩ} to spread from the newly inserted particle,Ĩ j = where J is the quantum number corresponding to the extra particle with rapidity θ 1 in state |ϑ, θ 1 . The shift of quantum numbers by a factor 1/2 causes an additional change to the rapidities of the state |{Ĩ}; J comparing with the state |{I} . In the thermodynamic limit this leads to a modification of a back-flow function. The back-flow functions obeys the modified integral equation where Θ H (θ) is the Heaviside function defined to be This is the back-flow function that we will use for bosonic theories. The back-flow function for the hole is again a minus back-flow function for a particle In the case of a multiparticle state, the total back-flow is still a sum of the single excitation back-flows In the rest of the manuscript we adopt the following shorthand notation. Given a set of quantum numbers {I} of the background state |ϑ we denote the finite volume counterpart of the excited state |ϑ, θ 1 , . . . , θ n as |{Ī}, I 1 , . . . , I n .
The bar reminds us that the set of quantum numbers {Ī} is shifted, according to the above prescription, due the presence of I 1 , . . . , I n . We will use the same notation also for theories with antiperiodic boundary conditions remembering that in that case there is no shift. Also, in what follows, we will write just F (θ|θ) and δ(θ) keeping in mind that both functions have to be choosen appropriately as well.

Form factors on a finite energy density background
We now want to find axioms for the (normalized) form factor, which we define as the function There are several reasons to choose this particular normalization. First, this normalization ensures that in the low density limit where ϑ(θ) → 0, the form factor reduces to the standard form factor defined on the vacuum, lim ϑ(θ)→0 f ϑ (θ 1 , . . . , θ n ) = f (θ 1 , . . . , θ n ).
The normalization (43) also ensures that the form factor axioms we will derive in this section look very natural and similar to the axioms obeyed by the standard vacuum form factors. Finally, we will show this normalization leads to a simple and natural expression for correlation functions. Our strategy is now the following. We start with axioms for form factors with a finite set of excitations at finite volume. The finite-volume form factors can be expressed in terms of infinite volume form factors, via Eq. (21). We then take the thermodynamic limit, taking the large-volume limit, but keeping the energy density constant. This leads to relations that we propose to be axioms for (thermodynamic) finite energy density form factors.
The finite size form factors we need can be written as The form factor {I}|O|{Ī}, {I i } contains no disconnected pieces, since the rapidities corresponding to {I} and {Ī} in the incoming and outgoing state do not match, as the rapidities in the incoming state have been displaced by the additional particles, {I i }.

Lorentz symmetry
It is no longer true that f O ϑ (θ 1 , . . . , θ n ), for a scalar operator, O, should be invariant under a Lorentz boost of the form θ i → θ i + α. This is because the new rapidities can be compared with those of the particles in the background. This implies that the one-particle form factor, f O ϑ (θ 1 ), is not necessarily a constant, yet the dependence on θ 1 is restricted by the periodicity axiom, as we shall see later.

Scattering axiom
We now investigate what happens when we exchange two adjacent particles in the form factor, At finite volume, the corresponding form factor is The form factor in the numerator in the right-hand side of (45) is the standard, infinite-volume form factor. Scattering in this form factor is a local process, between two neighboring particles, and it is not affected by the presence of background particles. The i-th and (i + 1)-th particles can be exchanged by multiplying with a factor of the S-matrix, S(θ i − θ i+1 ). The Jacobian in the denominator in (45) is symmetric under exchange of two rapidities. In the finite-volume form factor, we can therefore exchange two particles as This S-matrix exchange factor is not affected by the presence of the background. As we increase the system size, and introduce more background particles accordingly, this exchange factor remains the Figure 6: The periodicity axiom in the presence of the background is modified by an extra phase factor which depends on the background, the particle content and on the particle which is moved around. This extra phase factor comes from scattering of the particle going around with the background particles of the incoming and outgoing states.
same. In the thermodynamic limit, we can then write the scattering axiom as, see also Fig. 5, Therefore, it takes the same form as in the vacuum

Periodicity axiom
We now find what happens if we shift the last excited particle to become the first particle by applying crossing symmetry twice, We recall that the presence of the n-additional particles causes a shift in the background particles θ =θ + F (θ|{θ n })/Lρ s (θ). The finite-volume form factor can then be written as With this expression, we are now able to apply twice crossing symmetry on the n-th particle, and bring it to the position before the first excited particle. After doing this process of double crossing, the n-th particle needs to scatter with all the particles in the background in the outgoing and incoming states, such that we can write .
We notice that the two infinite products of S-matrices in (47) would completely cancel each other, if the set of background rapidities in the incoming state were not shifted by F (θ|{θ i })/Lρ s (θ). Taking the thermodynamic limit of expression (47), we then find the new periodicity axiom where and R ϑ (θ n |θ j ) is defined as This a good moment to recall our convention according to which functions appearing in the exponent have to be choosen appropriately depending whether the theory under the consideration is fermionic or bosonic. In any case function R ϑ (θ i |θ j ) has the following asymptotic behavior 4. Annihilation-pole axiom 2 In the presence of the finite-density background, two excitations in the form factor may still annihilate with each other, requiring an annihilation pole, but the residue of the form factor at this pole will now involve the scattering terms with the particles on the background. Taking the thermodynamic limit of the finite-volume expression, as was done for the previous axioms, it is easy to find that the annihilation pole axiom now becomes . . , θ n ) This axiom provides a recursive relation, connecting the n-particle and n − 2 particle form factors.

Cluster properties
Again starting with the finite-volume expression (21), we can study the cluster properties that follow from shifting the rapidities of some subset of particles, by the same amount, α → ∞.
It can now be easily derived that . Now taking the thermodynamic limit, we find where we point out that the second form factor on the right hand side of (53) is the standard, zerobackground form factor.

Bounds on growth
The same bounds apply on the form factors, as in the zero-background case. This can be seen from the finite volume expression for the form factors (21). The form factor on the numerator of the righthand-side of (21) has known bounds as we take the the rapidities of any of the excitations on top of the background to be large. It is clear from the properties of the determinant of the Jacobian, (52), that if one of the rapidities is taken to be very large, its contribution to this determinant simply factorizes. The behavior of the denominator then does not affect this bound.
The same bound then holds for the form factor even in the presence of the background,

Translation
To obtain the form factors of the operator at some space-time point, x, t, we can use Lorentz symmetry to write where E(ϑ), P (ϑ) are the energy and momentum of the background state, respectively, and E(ϑ; θ 1 , . . . , θ n ) P (ϑ; θ 1 , . . . , θ n ) are the energy and momentum of the background state with the particle excitations on top of it, respectively.
These total energies and momenta can be expressed in terms of the dressed energy and momentum of particles, which can be obtained from the TBA formalism [54], such that, where the dressed energy and momentum are given by (see Appendix A), respectively,

Analyzing the consequences of the form factors axioms
We will now explore the restrictions that the new form factor axioms place on the functions f ϑ (θ 1 , . . . , θ n ). We point out that the procedure outlined here yields the minimal solution to all the form factor axioms. Without any further evidence to the contrary, we assume the form factors are described by the simplest solutions which satisfy all the axioms. For instance, we assume the form factors have maximally analytic structure, with no zeros or poles which are not kinematically or otherwise required.

One-particle form factors
We begin by examining the one-particle form factor, f ϑ (θ), which as we have established, now may depend on the particle's rapidity, in contrast to the zero-background case where the one-particle form factor is a constant. The periodicity axiom for this form factor gives the condition We can search for a minimal solution to this equation, which we denote as the function f min ϑ (θ), which satisfies the condition while being maximally analytic. We notice that Equation (58) is very similar to the monodromy equation (14) satisfied by the minimal two-particle form factors at zero density background. We can therefore easily find the minimal solution to (58) in the same way. If we can find an expression then the minimal solution to (58) is Function R ϑ (θ|θ ) is a pure phase and therefore From eq. (59) and definition (50) of R ϑ (θ|θ) it follows that C ϑ (θ) is related to the Fourier transform of the back-flow function F (θ|θ). This implies that representation (59) is indeed possible. The general one-particle form factor can then be expressed as where K O ϑ (θ) is a periodic function, invariant under θ → θ + 2πi, which may depend on the properties of the operator, and should not have any poles (since the one particle state should not have any kinematical annihilation or bound state pole properties). These properties imply that we can express this function, in general, as the polynomial expansion where all of the information about the background, ϑ, and the operator, O, is contained in the coefficients, a O ϑ n , and the series is truncated at some operator and background dependent level, N O ϑ . From the cluster properties of the dressed form factors (53), we expect the one-particle form factor to satisfy where the right-hand-side involves the zero-background one-particle form factor, which is a constant and does not depend on the particle's rapidity, and the finite-background expectation value, which could be determined using the one-point Leclair-Mussardo formula. This provides a normalization for dressed form factors in terms of the zero-background form factors. For self clustering operators, such that The condition (62) provides an upper bound for the level, N O ϑ , of the expansion (61). The level N O ϑ is constrained by the fact that the function f O ϑ (θ) must go to a constant value at large rapidities, while the powers of hyperbolic cosine grow as lim θ→∞ cosh n θ → e nθ /2. If the minimal function f min ϑ (θ) does not exponentially decay at large rapidities, then only the constant term in (61) will survive (and the constant a O ϑ 0 can be fixed using (62)), or N O ϑ = 0.

Form factors with two or more particles
We can also attempt to find a minimal solution for the two-particle form factor by imposing the restrictions from the scattering and periodicity axioms. Combining these two axioms, we find a modified monodromy equation One can similarly derive a monodromy equation by applying the periodicity axiom on the θ 2 particle instead of the θ 1 particle. A minimal solution for the two-particle form factor satisfying the monodromy equation can be found of the form where f min (θ) is the standard (background independent) minimal form factor (15). We define the function A min ϑ (θ; θ 1 , . . . , θ n ) to satisfy the monodromy conditions for i ∈ [1, n], and R ϑ (θ|θ)R(θ|θ 1 . . . θ n )A min ϑ (θ + 2πi; θ 1 , . . . , θ n ) = A min ϑ (θ; θ 1 , . . . , θ n ).
A minimal solution to the monodromy equations (65) and (66) can again be found in terms of the function C ϑ (t) defined in (59), and a functionC ϑ (t|θ 1 , . . . , θ n ) defined by The minimal solution of (65) and (66) is then given by Again, representation (67) This procedure to find minimal solutions can be easily generalized to n-particle form factors, which can then be generally expressed as where K O ϑ (θ 1 , . . . , θ n ) is a periodic function in all the rapidities, symmetric under exchange of two rapidities, and contains all the information about the necessary annihilation and bound state poles, and any information identifying the operator, O.
Functions K O ϑ (θ 1 , . . . , θ n ) can be furher restricted. First of all, the annihilation pole axiom can be applied to derive recursive relations between the functions K ϑ (θ 1 , . . . , θ n ) and K ϑ (θ 3 , . . . , θ n ). Second, it is worth noting, that the location of annihilation poles is at the same values of rapidities as in the zero-background form factors. We can therefore assume the n-particle form factors are given by the expression (70), with the function where the denominator function contains all the necessary pole structure, and is not modified by the background density. The dependence on the background and on the operator is only in the numerator, Q O ϑ (θ 1 , . . . , θ n ). Unlike the zero-background case, the numerator Q O ϑ (θ 1 , . . . , θ n ) is not necessarily a function of only rapidity diferences, θ ij , but it can depend on each rapidity individually.
With this we conclude our general analysis of the form factors axioms. In section 8 we will compute the minimal form factors in the Sinh-Gordon model. We will also determine full form factors for the special class of the vertex operators. Before doing so, in the following two sections we consider the structure of the correlation functions in a state with a finite energy density.

Correlation function
In this section we find an expression for the two-point correlation function, in the presence of a finite density background Such correlation functions, evaluated on finite-energy-density states have many applications in strongly correlated systems. For instance, by properly choosing the function ϑ(θ), these correlators coincide in the thermodynamic limit with finite temperature [54], or generalized-Gibbs-ensemble averages [61]. The procedure, that we follow, is to write down first a correlation function in the infinite volume over a vaccum. Then we transform the formula to the finite volume. Finally we take the thermodynamic limit keeping the density of particles fixed. At infinite volume, the vacuum correlation function can be computed by inserting a complete set of intermediate states between operators In finite volume, rapidities are quantized and are intertwined with each other as solutions of the Bethe equations. Therefore to span the Hilbert space it is more convenient to use integer quantum numbers. The correlation functions is where E J and P J are the total energy and momentum, respectively, of the background set of particles {J}.
We are interested now in situations in which the state {J} in the thermodynamic limit follows the particle occupation number ϑ(θ) and define (75)  3 The additional labels on them, σ i = ±1, distinguish particles from holes. The positive sign stands for the introducion of a particle and a new quantum number, I i , which is not included in the set of background numbers, {J}. The negative sign, σ i = −1 stands for the introduction of a "hole", that is, removing a particle with quantum number I i , from the background, {J}. As mentioned above, for reasonable local physical operators we can assume that in the thermodynamic limit the matrix operators {J}|O|{J}, I where E J and P J are the total energy and momentum, respectively, of the particles {J}, in the presence of the additional n particles and holes. The sum {I σ } is over available integers, that is, for particles, σ k = +1 we sum only over integer values which are not included in the set, {J}, and for holes, σ k = −1, we sum only over the values of integers included in {J}.
We want now to take the thermodynamic limit of the correlation function. In that case, the background energy and momentum differences, E J − E J and P J − P J can be absorbed into the energy and momentum contributions of the excitations, by substituting with their dressed energy and momentum, ε(θ) and k(θ), respectively. Furthermore, the sum over integers, {I σ } can be replaced by a smooth integral over rapidities. This can be explicitly done using the technique introduced in [30], where the sum over some function f (θ k ) is expressed as where the contour of the integrals are around the possible integer solutions of the Bethe equations, where there is a simple pole in the denominator (these poles lie on the real axis of θ k ), and ρ σ k is a distribution ensuring that we only sum over available integer values for the given excitation. In the thermodynamic limit, these become ρ −1 (θ) = ρ p (θ), and ρ +1 (θ) = ρ h (θ). We have also defined the shorthand notation Q(θ k ) = Q k (ϑ, θ 1 , . . . , θ n ). Focusing on one type of excitations, particles or holes, the sum over the different contours can be replaced by a single contour encircling the real line, such that where we introduced for n being the order of the pole of f (z) atz. where the second term, and the fact that rapidity integrals are slightly shifted off the real axis, ensure the removal of any poles in the neighborhood of the real line (Γ ), which are not the poles corresponding to the integer solutions, I σ k . In our case, the function f (z) has such poles, which are interpreted as annihilation poles of the form factors, every time the rapidity of a given particle and a hole coincide. Since the terms in the correlation function involve the form factors squared, these will contribute double poles in Γ . In the thermodynamic limit, the integration below the real axis, R − i vanishes, such that We point out that the expression (81) is finite, and the possibly divergent contributions from the (double) annihilation poles are systematically regularized. It may still be a challenging task in practice to compute explicitly the contributions from the residues in the right-hand side of (81) (see similar calculations for the few-particle contributions in [30]), but the prescription is nonetheless finite and unambiguous. We can define the integration symbol, ffl , to denote this precise regularization where the integration is shifted from the real axis by i and the residues from the double poles are subtracted as prescribed.
The two-point function is then written as an expansion in terms of the dressed form factors, with particle and hole excitations on top of the background, ϑ.
The expression (83) can be contrasted with the LeClair-Mussardo conjecture for the finite-temperature two-point function: where only the zero-background form factors are used, and therefore the conjecture is not expected to be correct [25]. Furthermore, the original LeClair-Mussardo proposal did not show how to properly regularize singularities arising from annihilation poles, while our expression (83) is fully regularized.

Infrared regimes of the correlation function
Form factor expansions of correlation functions are known to be rapidly convergent in the infrared limit, in the zero-background case. In this case infrared refers to the limit of large separations between the operators (compared to the natural length scale defined by the mass gap, m). Defining the Euclidean distance, r = √ x 2 + t 2 , the zero-background two-point function may be written as At large values of mr, it is easy to see that this correlator decays exponentially as e −mr , and the insertion of intermediate states with higher number of particles leads to corrections in powers of e −mr . The exponential supression leads to this expansion being very quickly convergent. The leading exponentially supressed contribution comes from the one-particle form factor, and can be written explicitly as for large mr, where K α (x) is a modified Bessel function of the second kind. We now consider the correlation function in the presence of a background density, ϑ(θ), (83). In this case, there are two distinct IR limits we can consider. The first IR limit corresponds again to considering large separations between operators, or large values of mr. The second IR limit is the limit of small particle density, or small ϑ(θ) for all θ.
Let us discuss first in some details what happens for large separtion between the operators. In such situations the factors become highly oscillatory. We then face different scenarios. First, there might be a saddle point of the exponent which then localizes contribution to the correlation function to the vicinity of the saddle point configuration. This happens when we consider asymptotes of the autocorrelation (x = 0, large t) or large distance and time asymptotics with fixed ratio x/t. This leading behaviour is then due to particle-hole configuration which minimizes the exponent. However, there might be no saddle-point, as is the case for large x asymptote of an equal-time correlation function. Then still the particle-hole contribution arises then as the leading one, but now this is caused by the form factor itself. The form-factors have a kinematic pole whenever θ i → θ j + iπ. The contribution from the pole is excluded by the regularization procedure. However, in the vicinity of the pole the derivative of the form-factor squared is not continous. Therefore, the integration picks up contribution whenever particle and hole are close to each other. The leading term is therefore again controlled by the particle-hole contributions.
There might another source of non-smoothness of the integrand: a discontinuity of the filling function ϑ(θ). This is not the case at finite temperature, where the filling function is smooth, but e.g. at finite chemical potential. We don't consider this situation in the detail focusing solely on the finite temperature. In summary, for the thermal like states, the first infrared regime controlled by (88) is dominated by the particle-hole contributions.
Let us now consider the second infrared regime, when the filling function ϑ(θ) is small. To be concrete we concentrate on one simple and physically interesting example: the thermal equilibrium, with a finite temperature, T = 1/β. In this case, the function ϑ(θ) can be derived by minimizing the entropy of the system, subject to a constraint fixing the total energy density. The particle occupation number function is then given by [54] where ε(θ) is the dressed particle energy, which is determined as the solution of the integral equation At very low temperatures, or large mβ, the function (89) becomes small, and exponentially suppressed for all values of θ. One can then express the low temperature two-point function as an expansion in powers of e −mβ . The suppresion factor arises then through in the integrand in (83). It is important to notice that the two different factors that lead to suppressions in the IR limits, (90) and (88), have different behaviors, depending on whether the set of excitations are particles or holes, that is, depending on the signs σ i . At large mβ, the factor (90) leads to exponential supression of holes, while particle excitations are not supressed. On the other the space-time asymptotes are controlled by the particle-hole pairs and states with only particles or holes are supressed by oscillating factors (88).
Determining what are the leading contributions in the IR limit of the expansion (83) depends on whether mβ or mx is larger. In the doubly IR limit, with mβ → ∞, and mx → ∞, the two-point function (83) becomes simply lim β,x→∞ that is, only the disconnected piece of the correlator remains. We can then examine what are the leading corrections to (91) in the two IR regimes.
In the IR regime where mβ mx, the thermal corrections, e −mβ are much smaller than the spatial corrections, e −mx , therefore we can neglect contributions from hole excitations. The leading correction to (91) is then given by the one-particle intermediate state, In the IR regime where mβ mx, contributions from states with only particles or only holes are supressed exponentially at long distances, while contributions from mixed particle-hole states are not. The leading correction to (91) is then given by one particle-hole pair, Correlation functions in this regime, exhibit effectively gapless behavior, since there is no energy threshold needed to create a particle-hole pair. Depending on the properties of the background, ϑ(θ), the correlator may exhibit polynomial or exponential decay at long-distances. The regime, where particle-hole pairs are the dominant corrections, is analogous to the so-called hydrodynamic regime, that has been previously discussed in the non-relativistic Lieb Lininger model [44]. This is the only IR regime in the non-relativistic theory, since the number of particles is constant, so only particlehole pairs can be created. The IR regime described in (92) is not present in the non-relativistic theory, since it involves the creation of new particles. A proper non-relativistic limit of the QFT correlator in a thermal background, is then given by the limit of very large mass, m, particularly such that the supression factor (88) is much larger than the supression factor (90).

Comparison with previous low-temperature expansions
We now compare the infrared limits of our two-point function to a previously known low-temperature expansion. The first terms in a low temperature expansion of the thermal two point function, in an integrable QFT were rigorously computed in [30], by using a finite-volume regularization to take care of kinematical singularities, order by order.
This approach consists of a straightforward expansion of the thermal correlator in terms of the finitevolume form factors, with where |N L and |M L are N -and M -particle states, respectively, at finite volume. The main idea is that the kinematical singularities in the form factors in the numerator, C N M can be cancelled exactly, order by order, with analogous singularities arising from an expansion of ( N Z N ) −1 . The terms from the numerator and denominator can be combined into terms that contain no singularities, and are properly regularized in the L → ∞ limit, a procedure known as the "linked cluster expansion", where the details of these regularized terms D N M can be found in the original reference [30], and will not be discussed here. We only point out that the term D N M is suppressed at finite temperature by a factor of e −N βm cosh(θ) . This low-temperature expansion consists on computing the first few orders in N , while the label M controls the exponential suppression at long distances, x. At low temperatures, the leading terms are where the first sum, M D 0M is simply the vacuum two-point function, and M D 1M is the two-point function averaged on a one-particle state, and so on.
The two-point function formula we have presented, Eq. (83) involves a complete reorganization of terms, compared to the expansion (96). For example, it is easy to see that even the zeroth order term in our expansion, (91), involves a complete resummation of infinite terms of the expansion (96). To fully explore the differences between the two expansions, we can introduce the short-hand notation for Eq. (83), where G HP is defined as the term in the expansion (83), containing H holes and P particles. Each term in our expansion (98) involves an infinite resummation of terms with all powers of e −βm , so that (98) is not a low-temperature expansion, and the indices H, P are not in direct correspondence with N, M in (96). The reason for the discrepancy is that the expansion (98) is not centered around the vacuum state, but around the excited state, |ϑ , so it already starts from an infinite resumation of finite-temperature terms, even at zeroth order The way we can make a connection between the expansions (96) and (98) is by expanding each of the terms G HP in powers of e −βm , as where G N HP is the term in the expansion proportional to e −N βm . The correspondence between our expansion (98) and the linked cluster expansion (96) is then given by the relation which clearly indicates a complete reorganization of leading terms, between both expansions. This being said, the main advantage of our expansion, is indeed the fact that it is centered around a highly excited state, while (96) is centered around the vacuum. This means that even just computing a handful of leading terms in our expansion can give nonperturbative insight that is inaccessible to the expansion (96) unless one is somehow able to resum an infinite number of terms. One simple explicit example of this is the clear separation of the two regimes discussed in Eqs. (92) and (93). Additionally, it is completely trivial from our expansion to verify the clustering property, since this is just the zeroth order in our expansion. For Linked-Cluster and low-temperature expansions, the proof of this clustering property was shown in [34], but it was a highly non-trivial task.

Sinh-Gordon model
In this section we materialize the ideas introduced above in a concrete example of an Integrable QFT, namely the Sinh-Gordon (ShG) model. The action of the ShG model is where g is the interaction parameter. The particle content of the ShG theory is especially simple. There is only one type of particle, and these particles are also their own antiparticles: the field φ(x) is real. Moreover there are no bound states in the spectrum. Therefore the scattering matrix is just a single function, given by [62] S(θ) = tanh 1 The renormalized coupling constant B(g) is To formulate the axioms we need to know functions R ϑ (θ n |θ 1 , . . . , θ n ). The builiding block for them is the back-flow function which can be computed once the filling function ϑ(θ) is known. We formulate the finite volume theory with periodic boundary conditions. Therfore we will use the bosonic back-flow function F B (θ|θ). We will focus on the case of thermal equilibrium. In this case the filling function follows from the standard Thermodynamic Bethe Ansatz (TBA) which we now briefly recall.
The thermal equilibrium at inverse temperature β is described by the filling function with the dressed energy ε(θ) determined from the integral equation Knowing the filling function, the particle density function ρ p (θ) can be computed from eq. (22). The crucial ingredient for the TBA is the kernel ϕ(θ) (the derivative of the phase shift). From eq. (103) it follows that .
The TBA equations can be solved numerically. In Fig. 7 we show the resulting particles density function ρ p (θ) and the R ϑ (θ 1 , θ 2 ) factor. Recall that the R ϑ (θ 1 , θ 2 ) factor is a pure phase function that describes the contribution of the background particles to the effective scattering between particles with rapidities θ 1 and θ 2 . Equipped with the back-flow function we can now find the minimal solutions f min ϑ (θ 1 , . . . , θ n ) to the form-factor bootstrap. We follow the procedure outlined in section 5. We first analyze the one-particle form factors. We numerically find the representation (59) of R ϑ (θ|θ) in terms of C ϑ (t), that is Function C ϑ (t) enters then the expression for the minimal form-factor Figure 8: The one particle minimal form factor in the Sinh-Gordon theory for B = 0.5 and β = 1. Contrary to the vacuum case, the one-particle form factor depends on the rapidity θ. The asymptote lim |θ|→∞ f min ϑ (θ) = 1.
The results for the one-particle minimal form factor are presented in fig. 8. The general one-particle form factor takes the form Function K O ϑ (θ) is a polynomial in cosh θ. The order of this polynomial can be bounded from above by considering the θ → −∞ limit. The minimal form factor f min ϑ (θ) goes to 1 in this limit. This implies that the polynomial must actually be a constant. This constant can still depend on the operator and on the state |ϑ . Therefore The constant can be found from the cluster properties. For an operator O 1 that clusters into O 2 × O 3 the results is Finally, the general form of background dependent one particle form factor in the Sinh-Gordon theory is given that We can now consider the particular case of vertex operators, O = exp (iαφ), for some constant, α. In this case it is particularly simple to fix the one-particle form factor for two reasons: vertex operators are self-clustering, and one-point functions of vertex operators are computable exactly through the so-called Negro-Smirnov formula [60].
The one-particle form factors of the vertex operator are then The three constants involved in expression (115), namely 0|e iαφ |0 , ϑ|e iαφ |ϑ ϑ|ϑ , and f e iαφ 1 are exaclty known, and we include their expressions in Appendix B, for completeness. This completes the determination of  Figure 9: We plot the minimal two particle form factor normalized by the vacuum minimal form factor, defined as g min ϑ (θ 1 , θ 2 ) ≡ f min ϑ (θ 1 , θ 2 )/f min (θ 1 − θ 2 ). Both plots are again for the Sinh-Gordon model with interaction parameter B = 1/2 and at temperature β = 1.
one-particle form factors of vertex operators on a thermal background. We point out that it is easy to see that in the zero background limit, as is expected.
Computation of form factors of two or more particles follows in a similar way. We first compute the functionC ϑ (t|θ). It is given byC The building block A min ϑ (θ; θ 1 , . . . , θ n ) of the minimal form factor is then with C ϑ (t|θ 1 , . . . , θ n ) = n j=1 C ϑ (t|θ j ). Finally, the expression for the minimal form factors comes from combining the minimal vacuum form factors with buliding block A min ϑ (θ; θ 1 , . . . , θ n ). In the two particle case it is f min The standard vacuum minimal form factor f min (θ) in the Sinh-Gordon theory is with the normalization given by This minimal form factor satisfies lim θ→±∞ f min (θ) = 1.
Let us now consider the asymptotic of the background minimal form factor when θ 1 → −∞. The problem boils down to determining the asymptotes of A min ϑ (θ 1 , θ 2 ) in both variables. We have  Figure 10: We plot the minimal two particle form factor over the background (of temperature β = 1) and over the vacuum. Both plots are again for the Sinh-Gordon model with interaction parameter B = 1/2.
The limit with the second variable approaching the inifinity is lim θ 1 →−∞ This expression (127) can be further simplified by applying the cluster property again, to find We can then fix the normalization constant as which completely fixes the dressed two-particle form factors of vertex operators.

Conclusions and Outlook
In this paper we have introduced a set of axioms for finite background form factors of local operators in Integrable Quantum Field Theory in 2d. These form factors are building blocks for the computation of correlation functions in states of finite energy density, such as finite temperature states. We have also found the minimal form factors which are the maximal analytic solutions to the axioms. These results open a path towards an explicit computation of the correlation functions in states of finite energy density. Our work includes a number of simplifying assumptions. We focused on integrable QFT's with diagonal scattering, with a single species of particles and without bound states. The last two assumptions are easier to generalize, since it should not be too difficult to introduce an additional bound-state bootstrap axiom in the presence of a thermodynamic background, which takes care of new bound-state poles appearing in the S-matrix. The generalization to non-diagonal scattering theories seems much more challenging at present, since the thermodynamic Bethe ansatz formulation of non-diagonal theories is much more complicated, and requires the introduction of additional magnonic and string-like pseudoparticles [63]. While non-diagonal theories are too difficult for this introductory work, it is definitely an interesting direction to study in the future, since some physically relevant QFT's, such as the sine-Gordon model, and the O(N ) sigma model, have non-diagonal scattering [18]. Developing further in this direction will allow us to compare our approach to the numerical results available from TCSA computations [36].
Also for the sake of simplicity, we have chosen to study only operators which are local, relative to the particle excitations. This excludes theories which have soliton, or kink-like particles. While it is easy to generalize the standard form factor bootstrap on the vacuum to include semi-local operators, the generalization in the presence of a thermodynamical background is much more involved. This difficulty can be seen even in the simple case of Ising field theory [47], which is a theory of free Majorana fermions. The fermionic particles are not local relative to the spin operator, which makes the computation of finite-temperature form factors highly non-trivial, even for this non-interacting case.
As a simplifying step in our calculation, in equation (33) we assumed that the form factors on top of the background can be defined while ignoring the averaging over soft modes that lead to the same particle distribution. At very high energy densities, however, these soft modes may have a non-negligible effect, which can be interpreted as a renormalization of the form factor. At high energy densities, many soft particle-hole pairs (with momenta close to each other) can be created with no energy gap. Therefore the physical renormalized form factors would become those where one sums over soft particle-hole pairs leading to the same physical state. The computation of such form factor renormalization factors, which will be relevant at high temperatures, will be left for a future.
Formulating the axioms is of course only the first step. The next step is to actually determine form factors of interesting operators in a variety of theories. Here, we have shown, as a proof of principle, how to determine the one-and two-particle form factors of vertex operators in the Sinh-Gordon model. There exist many other integrable QFT's of physical and mathematical interest, which have been, over the course of many years, studied through the standard form factor bootstrap program (there are too many examples to refer individually to all here, but for an introductory overview of the applications of the bootstrap program in a variety of field theories, see Ref. [53]). After establishing our new set of axioms in the presence of a thermodynamic background, a large amount of work is set for the future, to re-derive the catalog of known exact form factors for the large number of IQFT's, but now in the presence of a thermodynamic background.
The main concrete applications we discussed here concern the simplest kind of finite-energy-density background, namely, the case of thermal equilibrium. An important task for the future will also be to apply the program developed here to other interesting cases, such as computing correlation functions after a quantum quench, perhaps by taking advantage of the quench action approach [40], which formulates non-equilibrium correlators in terms of representative eigenstates like the ones we considered here.
Another physical problem to address in the future is that of inhomogeneous quantum quenches, and in particular the problem of quantum transport. It has been recently understood that properties of quantities such as the Drude weight [64], which can be computed in terms of correlation fucntions in the presence of a spatially inhomogeneous background, can determine qualitative aspects of quantum transport in integrable systems, in particular, whether energy is transported ballistically or diffusively [65]. It would be interesting to study in the future these kind of correlations in the presence of an inhomogenous background, and give some insight into the transport properties of IQFT's.
Another interesting direction to consider in the future is the non-relativistic limit of the results derived in this paper. In particular, it is known that the Lieb-Liniger model can be obtained as a non-relativistic limit of the sinh-Gordon QFT, and several useful results have been recently derived using this correspondence [66,67,68]. It would be interesting to consider the non-relativistic limit of form factors derived for the sinh-Gordon model in the presence of a background, and see how (and if) we can recover thermodynamic form factors, for example, those that have been previously computed for the Lieb-Liniger model [44].
We will conclude by remarking that the standard form factor bootstrap on a vacuum state has seen too many applications in a variety of fields to provide a complete list here. The program we have outlined in this paper, which generalizes the bootstrap program to arbitrary highly excited states, should then also be of a similarly wide applicability, as we can now begin to explore the same kind of applications of QFT correlation functions, but in a variety of physically relevant thermodynamic states.

A Finite Volume Field Theory
In this appendix we derive the integral equations for the densities and the back-flow functions. These results are not new. They are included for the convenience to the reader.

Density equation
To derive the equation for the densities it is convenient to start with a finite size version of the density defined in terms of the rapidities solving the Bethe equations (20).
From the definition it follows that x(θ k ) = I j /L. Computing dx/dθ we find On the other hand dx dθ is the Jacobian of a mapping from the x-space to the θ space. Therefore it expresses the density of accessible rapidities in the vicinity of θ, that is the ρ s (θ). Taking the thermodynamic limit of the sum we find 2πρ s (θ) = m cosh θ +ˆdθ ρ p (θ )ϕ(θ − θ ).

Back-flow
The back-flow function describes a shift in rapidities due to an introduction of an extra particle or a hole. The equation for the back-flow can be derived by comparing the Bethe equations (20) for the background state |ϑ and an excited state |ϑ, θ 1 . Let {θ} be a set of rapidities of a finite version of the |ϑ state and {θ } together with θ be a set of rapidities of a finite version of |ϑ, θ . Compare now the Bethe equations for the corresponding rapidities in the two states. Because the quantum numbers are the same we get 0 = mL sinhθ k − mL sinhθ k + l =k δ 0 (θ k −θ l ) − δ 0 (θ k −θ l ) + δ(θ k − θ 1 ).
For the equality to be fulfilled the differenceθ k −θ k must be of order 1/L. We parametrize it, anticipating the result, in the following way,θ Substituting and expanding in inverse powers of L, in the leading order, we find In the thermodynamic limit, the bracket becomes equal 2πρ s (θ k ) and the second sum turns into an integral. Reorganizing we find This is the equation for the fermionic back-flow function. The bosonic back-flow function F B (θ|θ 1 ) can be found in analogous way. The back-flow can be used to express the energy (θ) and momentum k(θ) of the excited state |ϑ, θ with respect to the background state |ϑ . From the finite-size expressions for the energy and momentum, in the thermodynamic limit, we find k(θ) = m sinh θ −ˆdθ ϑ(θ ) cosh θ F (θ |θ), (θ) = m cosh θ −ˆdθ ϑ(θ ) sinh θ F (θ |θ).
The energy and momentum of a state with more excitations is the sum of the energy and momentum carried by each single excitation. At the thermal equilibrium and for fermions, these functions can be expressed in a standard TBA [54] form ε(θ) = m cosh θ +ˆd θ 2π ϕ(θ − θ )ϑ(θ )ε(θ ),

B Normalization constants in Sinh Gordon
In this appendix, for the sake of completeness, we list several previously derived formulas that we make use of, to normalize the one-and two-particle form factors of vertex operators in the Sinh-Gordon model. In both, the one-and two-particle cases, the normalization constants were found in terms of three values, namely, the vacuum expectation value of the vertex operator, the one-particle form factor (on top of a vacuum state), and the one-point function at finite background, of the same operator.
An integral expression for the vacuum expectation value was found in [69], and it is given by where q = g + g −1 .
An integral expression has also been found for the one-particle form factor in [70], and it is given by .
The computation of the expectation value in the presence of the background can be performed in two ways. The first approach is through the one-point Leclair-Mussardo formula [25], which has the advantage of being easily generalizable to any local operator of any field theory with diagonal scattering. This Leclair-Mussardo formula, for a given operator, O(x), and a background given by the distribution, ϑ(θ), is given by The disadvantage of the Leclair-Mussardo formula is that, in practice, it can only be evaluated as a lowenergy density expansion, by computing one by one the few-particle form factor terms. It is not clear if a complete resummation of all terms is possible, to yield any closed form expression valid at higher energy densities. A more convenient expression is available, proposed by Negro and Smirnov [71,72], exclusively for the finite-background expectation values of vertex operators in the sinh-Gordon model. While being specific to one kind of operator in a particular QFT, the Negro-Smirnov formula has the advantage of being a fully resummed integral expression, that does not need to be computed order by order like the Leclair-Mussardo formula. The one-point functions are given by the integral expression ϑ|e (k+1)gφ |ϑ ϑ|e kgφ |ϑ = 1 + 2 sin(πB(k + 1/2)) πˆ∞ −∞ dθ ϑ(θ) e θ p k (θ), where p k (θ) = e −θ +ˆ∞ −∞ dµ ϑ(µ) χ k (θ − µ)p k (µ), and χ k (θ) = i 2π e −ikBπ sinh(θ + iπα) − e ikBπ sinh(θ − iπα) .
and B(g) is defined in Eq. (104). The formula (144) only computes ratios of expectation values of vertex operators. Nevertheless, as was pointed out in [43], if B(g) is an irrational number, this formula can be used to compute the expectation value of the vertex operator for arbitrary k, in a procedure outlined in [43].