On Factorizable S-matrices, Generalized TTbar, and the Hagedorn Transition

We study solutions of the Thermodynamic Bethe Ansatz equations for relativistic theories defined by the factorizable $S$-matrix of an integrable QFT deformed by CDD factors. Such $S$-matrices appear under generalized TTbar deformations of integrable QFT by special irrelevant operators. The TBA equations, of course, determine the ground state energy $E(R)$ of the finite-size system, with the spatial coordinate compactified on a circle of circumference $R$. We limit attention to theories involving just one kind of stable particles, and consider deformations of the trivial (free fermion or boson) $S$-matrix by CDD factors with two elementary poles and regular high energy asymptotics -- the"2CDD model". We find that for all values of the parameters (positions of the CDD poles) the TBA equations exhibit two real solutions at $R$ greater than a certain parameter-dependent value $R_*$, which we refer to as the primary and secondary branches. The primary branch is identified with the standard iterative solution, while the secondary one is unstable against iterations and needs to be accessed through an alternative numerical method known as pseudo-arc-length continuation. The two branches merge at the"turning point"$R_*$ (a square-root branching point). The singularity signals a Hagedorn behavior of the density of high energy states of the deformed theories, a feature incompatible with the Wilsonian notion of a local QFT originating from a UV fixed point, but typical for string theories. This behavior of $E(R)$ is qualitatively the same as the one for standard TTbar deformations of local QFT.


Introduction
The so-called TTbar deformations [1,2] of two-dimensional quantum field theories (QFTs) has brought about a renewed interest to UV properties of Renormalization Group (RG) flows generated by higher dimensional (a.k.a. "irrelevant") operators. The TTbar deformation is defined as the one-parameter family of formal "actions" A α , determined by the flow where (TT ) α (x) is a special composite operator built from the components of the energy-momentum tensor associated with the theory A α [3]. The deformation (1.1) has a number of notable properties. The theory A α is "solvable", in the sense that certain characteristics can be found exactly in terms of the corresponding ones in the undeformed theory A α=0 . This is remarkable, because the deformation operator (TT ) α has exact dimension 4, meaning the perturbation in (1.1) is "irrelevant" in the RG sense. Normally, such deformations are expected to break the short-distance structure of the quantum field theory, generally rendering the theory UV incomplete, and possibly violating causality at short scales. The abnormal UV properties of the theory A α are manifest already in the short-scale behavior of its finite-size ground-state energy. If the spatial coordinate of the 2D space-time is compactified on a circle of circumference R, its ground-state energy E α (R) is determined exactly, via the equation [1,2] E α (R) = E 0 (R − αE α (R)) , (1.2) in terms of the ground state energy E 0 (R) of the undeformed theory, at α = 0. The equation (1.2) shows that, depending on the sign of the deformation parameter α, the ground state energy either develops a square root singularity at some R * ∼ 1/ |α|, or has no short-distance singularity at all. Neither of these types of behavior is compatible with the idea of QFT as the RG flow stemming out of a UV fixed point. The theory defined by (1.1) therefore is not a local QFT in the Wilsonian sense [4]. Moreover, at negative α, the singularity at finite R signals a very fast growth of the density of states at high energies, a common hallmark of string theories, leading to the Hagedorn transition [5]. The behavior of E α (R) at positive α is possibly even more puzzling, as it suggests a finite number of states per unit volume, an unlikely feature if one thinks of a QFT as a system of continuously many interacting degrees of freedom, unless quantum gravity is involved 1 . Therefore, the deformed theory determined by (1.1) cannot be considered a conventional UV complete local QFT. At the same time, however, the TTbar deformation has a number of robust features which makes one reluctant to simply dismiss it as "pathological". It is instead tempting to think that the deformation (1.1) exemplifies some meaningful extension of the notion of local QFT.
In particular, an interesting interpretation of the theory A α in terms of its gravitational dual was proposed in [8], where a relation to the state of the bulk gravity in the dual theory was suggested. Several questions about 2D physics of the deformed theory need to be elucidated in order to put such suggestions on a solid ground. For example, does the deformation preserve any part of the local structure of QFT? Notice how the very definition (1.1) depends on the notion of the energy-momentum tensor, conventionally a part of such a local structure. Another important question concerns the macro-causality in 2D space-time. While the deformation (1.1) with positive α is suspected to display super-luminal propagation [6,8], the case of negative α is most likely free from this problem. We will not dwell on this question, as it is the negative-α deformation which will be of interest to the present discussion. In any case, we believe it is important to understand the physical origin of the above abnormal short-distance properties.
Another exact result about the theory A α concerns the deformation of its S-matrix, whose elements differ from the corresponding undeformed ones by a universal phase factor, available in closed form [6]. In particular, the 2 → 2 elastic scattering amplitude has the form S α (θ) = S 0 (θ) exp −iαM 2 sinh θ , (1.3) where S 0 (θ) = S α=0 (θ) is the 2 → 2 scattering amplitude of the undeformed theory.
Here θ = θ 1 − θ 2 is the difference between rapidities of the two particles involvedassumed for simplicity to be identical -and M denotes their mass; in what follows we set the units so that M = 1. A notable feature of the additional phases acquired under the deformation is their abnormally fast high-energy growth, which is evident already in the form (1.3) 2 . The scattering phase in (1.3) determines the density of two-particle states, suppressing it when α > 0 but greatly enhancing it at negative α. In the latter case, one might be led to believe that the Hagedorn behavior is directly related to this rapid growth of the 2 → 2 scattering phase. One of the results of the present work is to show that the situation is more subtle: the growth of the two-particle scattering phase in (1.3) is not a necessary condition for the formation of the singularity of the finite-size energy at finite real R. We will study certain generalizations of the TTbar deformation which can be defined whenever the original QFT is integrable [1]. In most of such deformations, the scattering phases present a less exotic high-energy behavior -i.e., they have finite limit at θ → ∞ -while, at the same time, the overall density of states grows nonetheless exponentially with the energy, leading to the Hagedorn singularity. The generalizations of the TTbar deformations we will be interested in are based on the integrability of the original QFT. This assumes that the theory possesses infinitely many conserved local currents of higher Lorentz spins s + 1, with s taking values in the set {s} of odd natural numbers: s = 1, 3, 5, 7, ... 3 . The deforming operators TT (s) (x) are constructed from these currents in the exact same way as the operator TT (x) is built from the energy-momentum tensor, see [1] for details. It can be then shown that the theory deformed by adding such operators retains its integrability, preserving the same set of conserved local currents. Therefore the deformations of an Integrable QFT (IQFT) by the operators TT (s) generate an infinite-dimensional family of flows generalizing (1.1),  1). The latter corresponds to the special case α s = 0 for s > 1, and α 1 = α. To distinguish them, below we often refer to (1.1) as the "TTbar proper", or simply TTbar, reserving the term "Generalized TTbar" to the generic deformation (1.4). It was argued that the deformation (1.4) leads to the following deformation of the elastic two-particle S-matrix with the same notations as in (1.3) and (1.4) 5 . The phase factor Φ {α} (θ) is known with the name of CDD factor [13]. Generally, it is an energy-dependent phase factor Φ(θ) 3 Generally, the set of spins {s} of local Integrals of motion may be different in different integrable theories. Here we assume, again for simplicity, the most common situation -represented e.g. by sinh-Gordon or sigma models -where {s} involves all odd natural numbers. In different models the CDD factor discussed below may be constrained by additional conditions, which however do not change the overall conclusions below. 4 In [10], a different family of generalizations of the TTbar flow, in which the deforming operators TTs are asymmetrically constructed from the energy-momentum tensor and a higher-conserved current, was explored.  {α} (x) is chosen, otherwise the terms in the sum in (1.3) would have additional normalization-dependent numerical coefficients. The form (1.3) was explicitly derived in [1] for the deformed sine-Gordon model, to leading order in the deformation parameters. However, this form of the S-matrix deformation under the flow (1.4) can be proven in the general case, using the methods of [11] or the approach developed in [12]. We will elaborate this point elsewhere.
which can be added to the 2 → 2 scattering amplitude without violating the analyticity, unitarity and crossing symmetry conditions. The unitarity and crossing demand that Φ(θ) satisfies the functional relations which Φ {α} (θ) in (1.5) obviously does term by term in the sum over s. Moreover, it is easy to see that (once the overall sign ambiguity is ignored) any solution of (1.6) can be represented by the form (1.5), with the series in the exponential converging in some vicinity of the point θ = 0. However, the series does not need to converge at all θ. The S-matrix analyticity forces Φ(θ) to be a meromorphic function of θ, with the locations of the poles constrained by the condition of macro-causality (more on this momentarily). Therefore, for (1.5) to represent a physically sensible S-matrix, the sum over s is allowed to have a finite domain of convergence, while its analytic continuation must admit the representation where the first factor absorbs all the poles located at finite θ, whose number N is in general arbitrary (possibly infinite), In this last factor, the series in the exponential is assumed to converge at all θ, so that Φ entire (θ) represents an entire function of θ. Macro-causality restricts possible positions of the poles θ p to either the imaginary axis Re θ p = 0, or to the strips Im θ p ∈ [−π, 0] mod 2π since, in virtue of (1.6), Φ(θ) is a periodic function, Φ(2πi + θ) = Φ(θ). Let us stress here that the representation (1.7-1.9) of the generic CDD factor Φ {α} (θ) differs from the one given in (1.5) only in the parameterization: any factor (1.7-1.9) can be written in the form (1.5), with the parameters α s expressed in terms of a s and θ p , and conversely any factor Φ {α} (θ) defined in (1.5), being analytically continued to the whole θ-plane, can be written in the form (1.7).
In the present work we focus our attention on the class of S-matrices (1.5) having CDD factors (1.7) for which the entire part (1.9) is absent 6 , 10) and the product in (1.8) involves finitely many factors, i.e. N < ∞. Note that, unlike (1.3), such CDD factors have regular limits at θ → ±∞. Therefore, if the undeformed S-matrix S 0 (θ) behaves regularly -presenting no abnormal growth of the scattering phase -at large θ, so does the deformed S-matrix S 0 (θ)Φ(θ). We now raise the following question: how does an S-matrix deformation such as the one just described affect the short-distance behavior of the theory? Unfortunately, for the general TTbar deformation (1.4) no closed form of the finite-size energy levels similar to (1.2) is available with which one could analyze their dependence on the size R of the system. However, having an exact expression for the deformed IQFT S-matrix, the finite-size groundstate energy E(R) can be obtained by solving the associated Thermodynamic Bethe Ansatz (TBA) equation [15,16]. In general, the form of the TBA equations depends on the particle spectrum of the theory. Here we consider, for simplicity, the case of a factorizable S-matrix involving only one kind of particles, having mass M = 1. In this case the two-particle S-matrix consists of a single amplitude S(θ), which itself satisfies the equations (1.6). Therefore we can limit attention to the functions S(θ) of the form (1.8) 7 . There are two substantially different cases, depending on the sign of S(0) = σ = ±1. Following [16], we refer to these cases as the "bosonic TBA" when σ = +1 and "fermionic TBA" for σ = −1. Given S(θ), let ϕ(θ) be the derivative of the scattering phase, Then the TBA equation takes the form of a non-linear integral equation for a single function (θ), the pseudo-energy, The ground state energy can then be recovered from the pseudo-energy via the following expression (1.14) In most cases the TBA equations are impervious to the exact analytic derivation of their solutions but are amenable to numerical approaches. These can yield important insight into high energy, viz. short distance, properties of the deformed theories (1.4). A numerical solution can be obtained, with practically arbitrary accuracy, by numerical integration of (1.12). This approach was employed to obtain E(R) in a number of IQFT's with known S-matrices, see e.g. [16,17]. Usually, the numerical solution is obtained by iterations, starting from a seed function, conventionally taken to be (θ) = R cosh θ, and successively substituting the result of the previous iteration in the righthand-side in (1.12). We will review this approach in §3.1. If one considers the S-matrix associated with a UV complete local IQFT -such as a conformal field theory (CFT) perturbed by a relevant operator, the sine-Gordon model, or an integrable sigma-model -the iterations turn out to converge for all R > 0, and the resulting ground-state-energy E(R) happens to be analytic at all positive real R, developing a Casimir singularity at R = 0. But how adding a CDD factor to the S-matrix will affect the TBA solution? This question was addressed in the early 90's by Al. Zamolodchikov, who has considered the modification of the trivial fermionic S-matrix S(θ) = −1 by the simplest possible rational CDD factor, namely (1.8) with N = 1. In the resulting theory, the celebrated "staircase model" [17], the iterative solution of the TBA still converges at all positive R, producing a ground-state-energy E(R) analytic for R > 0. He also observed that when adding more general CDD factors the situation changes qualitatively. Typically, the convergence of the iterative solution breaks down at R below a certain critical value R * , and the form of the numerical solution at R > R * , where the iterations converge, strongly indicates the existence of a square-root singularity of E(R) at R * [18]. A similar observation was made in [19], where a particular CDD deformation of the trivial bosonic S-matrix S(θ) = 1 was studied and the numerical solution of the associated TBA equation was found consistent with the existence of a singularity at finite R * > 0. We wish to stress that the presence of the singularity at finite R * and, moreover, its square-root character, are features very similar to the ones displayed by E(R) in the TTbar deformed QFTs, as shown in Fig 2 below.
In this work we study a few simple cases of CDD deformed TBA equations, using a refined numerical routine based on the so-called "pseudo-arc-length continuation" (PALC) method. This allows one to recover solutions to the TBA equation (1.12) which are unstable under the standard iterative approach. This method is explained in detail in §4. The object of our attention will be trivial S-matrices S(θ) = σ = ±1 deformed by CDD factors (1.8) with N = 1, 2. The case N = 1 with σ = −1 corresponds to either the sinh-Gordon or the staircase model, depending on the position of the pole. As mentioned just above, these models do not display any abnormal short-distance behavior and were extensively studied in the literature. The bosonic TBA with N = 1 was considered in [19] and we will comment on it in §5, along with the N = 2 case. Of these, we mostly address the fermionic cases, although some results for the bosonic TBA are also presented. We find that for all allowed values of the parameters θ p (p = 1, 2) the fermionic TBA equation (1.12) with sufficiently large R possesses two real solutions, or "branches", which merge at some finite R = R * . For R < R * these branches are likely to continue as a pair of conjugated complex-valued solutions. Of these two real solutions at R > R * , one reproduces the iterative solution of the TBA equations (1.12). We will call this solution the "primary branch", while referring to the other one as the "secondary branch". Let us stress here that it is the primary branch which directly corresponds to the deformed theory: E(R) on the primary branch represents the finite-size vacuum energy of the deformed theory (in particular, at R → ∞ the effect of the deformation disappears, as expected); it also gives the specific free energy of the deformed theory at temperature T = 1/R (in particular, it is the primary branch solution which correctly sums up the virial expansion associated with the input particle theory). In this sense, one could call the primary branch the "physical" one, although we will not use such a term 8 . The secondary branch always has lower energy E(R) than the primary one, which is qualitatively similar to the behavior observed in the TTbar deformations with negative α, see Fig 2. Since the two branches merge at some finite R = R * , this can be regarded as a "turning point", where the continuation along the graph of E(R) turns backward into the secondary branch. This is precisely the kind of situation the PALC method is designed to deal with. The secondary branch remains real for all R > R * and, moreover, develops a linear asymptotic ∼ e ∞ R as R → ∞. This, again, is in full qualitative agreement with the TTbar deformations, together with the important fact that the singularity of the pseudo-energy (θ|R), viewed as a function of R, occurs at R = R * that is independent of θ. Of the above features, the existence of primary and secondary branches with the turning point at finite R * , independent from θ, repeat verbatim in the bosonic 2CDD model. On the other hand, we still cannot check the large R behavior of the secondary branch with sufficient accuracy, due to some instability in the numerical procedure. We will return on this problem in a future work.
It is likely that the general situation displayed in the models studied here, i.e. the solution of the TBA equation developing a square-root singularity at finite R * , which signals the presence of a Hagedorn transition, remains qualitatively the same when more CDD poles are added in (1.8) -with the possible exceptions of special domains, hypersurfaces of lower dimension, in the parameter space 9 . This of course will have to be carefully verified. We regard the present work as a first step in the program of systematically studying the short-distance behavior of the generalized TTbar deformations (1.4) of IQFTs. The qualitative similarity to the TTbar-deformed QFTs, with 8 The reason is that this would imply that the secondary branch is "unphysical", which we are reluctant to claim. Although the secondary branch definitely does not have direct interpretation in terms of "physics" of the input S-matrix, it might very well have some physical content of its own. In fact, understanding physical mechanism behind the secondary branch is one of the outstanding problems which remains open both for the generalized TTbar deformations and for the TTbar proper. 9 Examples of such cases can be found in [20,21,22].
negative α, suggests that the same mechanism behind the formation of the Hagedorn singularities is at play in all of these models. Understanding the physics underlying this phenomenon remains the most important open problem in this context, as well as the main motivation for the present work.
2 From TBA to Hagedorn: the TTbar case Henceforth we will assume that the theory under consideration is integrable, with a factorizable S-matrix. Let us briefly remind how, in this case, equation (1.2) can be derived from the S-matrix deformation (1.3) via the TBA equations. We will present a somewhat simplified version of the much more general arguments of [2] (for related work see [23,24] and the more recent [25,26]). Whereas the analysis in [2] applies to all the energy eigenvalues of the TTbar deformed theory (1.1), we limit our considerations to the ground-state energy, which we denote as E(R). The advantage is that the simple arguments presented below apply to the deformation (1.3) of an essentially generic integrable theory. The only assumptions, made for simplicity, are that the particle scattering theory associated with A 0 involves only one kind of neutral particles, with the factorizable scattering of fermionic type 10 , i.e. S 0 (0) = −1. The goal is to emphasize some important properties of the solution which, as we will see, are shared by the TBA solutions by more general CDD deformations. The TBA equation (1.12) associated with the deformed S-matrix (1.3) has the following kernel Recall that the ground state energy E α (R) is given by (1.14), which in our case reads where L α (θ|R) := log 1 + e − α(θ|R) satisfies the deformed TBA equation (1.12), Due to the fact that the pseudo-energy is even, as is easily shown, we can separate the dependence on θ and θ in the rightmost term in the kernel (2.1) so that the TBA equation can be written as follows where we used the definition (2.2). For reasons that will become clear shortly we have made explicit the fact that (θ|R) and L(θ|R) are functions of R as well as of the rapidity θ. This last form (2.4) shows that α (θ|R) satisfies the same TBA equation as . It then follows that which immediately implies the equation (1.2) for the deformed energy. It is also worth reminding here how the singularity of E α (R), signifying the Hagedorn density of states, follows from (1.2). This takes a particularly simple form in terms of the function R α (E), inverse to the function E α (R), where α is regarded as a fixed parameter, This expression shows that the graph of the deformed function E α (R) differs from the graph of E 0 (R) just by an affine transformation (R, E) → (R + αE, E) of the (R, E) plane. If we assume, as we do, that the undeformed theory A 0 is a conventional QFT, definedà la Wilson as the RG flow from some UV fixed point down to an IR one (see [4]), then the graph of E 0 (R) looks qualitatively as shown in  Its R → 0 behavior −πc/6R is controlled by the UV fixed point. At large R, E 0 (R) shows the linear behavior ε 0 R, with the slope ε 0 representing the bulk vacuum energy density. We have to stress that the TBA equations actually compute the difference E vac (R) − ε 0 R, and in our subsequent analysis E(R) stands for this difference. (That is why in all plots below the R → ∞ slope of the primary branch is always set to zero.) At large R the function E 0 (R) approaches a linear asymptotic ε 0 R, where ε 0 is the vacuum energy density of the infinite system, with the rate of the approach controlled by the IR fixed point, which, typically, is a non-critical one. On the other hand, at R → 0 it diverges as the Casimir energy determined by the UV fixed point, E 0 (R) → −πc/6R, where c is the Virasoro central charge of the UV fixed point CFT. Then, according to (2.6), the plot of E α (R) will look as either one of the panels a) or b) in Fig 2, depending on the sign of α. In what follows we will concentrate our attention to the case of negative α, shown in panel a). Note that the curve E α (R) has two branches, each of which having real values for R above a certain critical value R * . It is the upper "primary" branch that corresponds to the ground state energy of the TTbar-deformed theory (1.1). The two branches merge at R = R * , where the function E α (R) develops a squareroot branch point, i.e. the derivative dE α (R)/dR diverges as (R − R * ) −1/2 . At R < R * , the analytic continuation of E α (R) returns complex values and the two branches are complex conjugate.
It is the singularity at R * that signals the Hagedorn phenomenon in the deformed theory, which can be inferred as follows. When the Euclidean theory is considered in the geometry of a very long cylinder of circumference R, as shown in Fig 3, its partition function Z is saturated by the finite-size ground state where L → ∞ is the length of the cylinder. This corresponds to the picture in which the coordinate y along the cylinder is taken as the Euclidean time. Alternatively, if one uses the picture where x plays the role of Matsubara time, the same partition function The coordinate x is compactified on a circle of circumference R, while the length L of the cylinder is assumed to be asymptotically large. In the picture where y is regarded as the Euclidean "time" the partition function (2.7) is dominated by the finite-size ground state contribution. In the complementary picture, where y is interpreted as spatial coordinate while x plays the role of the Matsubara "time", the same partition function is given by the thermal trace (2.8).
is represented as the trace where dEΓ(E) ∼ dE e S(E) denotes the density of states, i.e. the number of states in the energy interval dE. While in a local QFT whose high-energy limit is governed by the UV fixed point the entropy S grows as S(E) 2πc/3 √ LE as E → ∞ -this is known as Cardy formula [27] -the singularity of F (R) at finite positive R = R * is formed when the entropy S(E) grows much faster, so that the partition sum diverges at R < R * . We will discuss the density of states in the TTbar deformed theories in more details in §6.
In the above discussion we have denoted by R * the position of the singularity of E α (R) as a function of R. It is important to observe that the solution α (θ|R) displays a singularity at the same position R = R * , independent on the value of the rapidity θ. In other words, in the two-dimensional space spanned by the variables (θ, R) the singularity of α (θ|R) occurs along the line (θ, R = R * ). We will see that this feature of the singularity associated with the Hagedorn transition will be reproduced in the generalized TTbar flows studied below.
As already mentioned, the enhancement of the density of states in the deformed theory is well expected. The scattering phase in (1.3) grows fast with the center of mass energy, leading to the increase of the density of two-particle states, implying a yet greater increase of the density of all multi-particle states. The calculation presented above demonstrates that the resulting entropy displays the Hagedorn behavior (2.9). It is then tempting to assume that the formation of the Hagedorn density (2.9) is directly related to the fast growth of the scattering phase. In the next section we will show that the Hagedorn singularity develops just as well in models whose associated CDD factor has a finite behavior at high energies, as in (1.8) with finite N , which indicates that the physical origin of the Hagedorn transition in the deformed theories is substantially more intricate.

The models
Here we study the CDD deformations of the trivial (fermionic or bosonic) S-matrix by the pole factor (1.8), which we write as where, as before, σ = − (resp. σ = +) corresponds to the fermionic (resp. bosonic) case. The parameters u p may be taken to be complex and, in view of the obvious periodicity of S(θ), we may limit our attention to the strip −π < Re(u p ) < π. The standard analytic requirements for the physical S-matrix, however, impose restrictions on the possible locations of the poles θ p = iu p . Taking these restrictions into consideration, the parameters u p are allowed to be either real or complex with negative real parts. The poles θ p = iu p , with real positive u p signal the existence of bound states -new stable particles of mass 2M cos(u p /2). Since the presence of such particles violates our working assumption that the mass spectrum of the theory only involves a single kind of stable particle with mass M , henceforth we will assume that all parameters u p in (3.1) possess a negative real part 11 : This leaves us with poles θ p = iu p lying in the unphysical region, i.e. the region of the complex center-of-mass energy s-plane reached by analytically continuing the scattering amplitude through the two-particle branch cut. When u p has nonzero imaginary parts, such poles are associated to unstable particles, having complex masses M p = 2M cos(u p /2), with the real and imaginary parts identified as usual with the mean center of mass energy and the width of the resonances. The poles with real negative u p do not have clear particle interpretation, but the number of such poles signify the increment of the scattering phase as the function of θ at low energies; these poles are often referred to as the virtual states (see e.g. [28]). A final requirement is that of unitarity of the physical S-matrix, which demands that S(−θ) = S * (θ) at all real θ, or, equivalently, that S(θ) takes real values at pure imaginary θ. It follows that any non-real parameter u p in (3.1) either has fixed real part Re(u p ) = −π/2 or appears together with its conjugate u * p . We can then refine the range (3.2) to the following three cases Thus, each subfamily (σ, N ) of (3.1) contains a number of, in principle, different models, determined by a given combination of the ranges (3.3) for each of the parameters {u p } N p=1 . Some simple combinatorics 12 tells us that this number is Since for any model determined by (3.1), with parameters in the ranges (3.2), the mass spectrum contains a single stable excitation, the resulting single-particle TBA equation takes the simple form (1.12), with the kernel ϕ(θ) being the derivative of the scattering phase which, in the case of (3.1), explicitly reads An equivalent, sometimes more useful, expression of this kernel is its partial fractions expansion .
(3.6) 12 Given the number N of poles one needs to partition it into three non-negative integers na, n b and nc with the constraint that na + n b + 2nc = N . Once a value of nc = 0, 1, . . . , N/2 is chosen, one is obviously left with N − 2nc + 1 non-equivalent arrangements of poles between the cases a) and b). Thus the number of different models is given by In what follows, we are going to concentrate our attention on two particular subfamilies: the "1CDD models" where N = 1 and the "2CDD models" with N = 2.
The 1CDD models When N = 1 the S-matrix (3.1) consists of a single factor According to the breakdown of cases (3.3), for each choice of the TBA statistics we only have two possible models, corresponding to the following ranges of the parameter u 1 : (a) u 1 ∈ R and −π < u 1 < 0, Considering at first the fermionic case σ = −1, one recognizes in (3.7), for the case (a), the well-known S-matrix of the sinh-Gordon model On the other hand, the case (b) corresponds to the S-matrix of the "staircase model", introduced in [17] S stair (θ) = sinh θ − i cosh θ 0 sinh θ + i cosh θ 0 , θ 0 ∈ R . For what concerns the two bosonic 1CDD models, the solution of the TBA equation has a considerably more intricate behavior. The case (a) of u 1 real was first addressed in [19], where it was observed that the iterative solution of the TBA equation only converges for sufficiently large radius R > R * > 0. The authors also noticed that the function E(R) appears to develop some sort of singularity at R = R * . Below in §5 we will show that the solution to the TBA equation, and, consequently, the ground state energy E(R), possesses, as a function of R, two branches. These merge at R = R * , meaning that R * is a square-root branching point. We also show that this behavior extends to the case (b) of complex parameter u 1 = −π/2 + iθ 0 .
The 2CDD model In the N = 2 subfamily, a pair of CDD factors is present in (3.1): Following the breakdown (3.3), we see that there are 4 possibly distinct models, corresponding to the following ranges of the parameters u 1 and u 2 (a) u 1 ∈ R and −π < u 1 < 0, u 2 ∈ R and −π < u 2 < 0, The model (a) can be considered as a special instance of the more general case d). On the other hand the models (c) and (b) -equivalent to (b') -are genuinely distinct. All the models above display, both for the bosonic and fermionic statistics, the same type of behavior observed in the bosonic 1CDD models mentioned above: the iterative procedure for solving the TBA equation (1.12) only converges for R larger than some positive value R * > 0 and the ground state energy E(R) apparently develops a singularity at R = R * . While we are going to present some data for all the various 2CDD cases, we devoted most of our attention to the case (d), that we will call, with some definitional abuse, the "2CDD model". Its S-matrix and TBA kernels explicitly read as follows . (3.12)

Iterative solution
The chances of a non-linear integral equation of the form (1.12) to be amenable to an explicit analytic solution are considerably slim. For this reason the main investigation approach to the TBA equations is of numerical nature 13 . In most situations, a simple iterative procedure of the following type appropriately discretized, is shown to converge to the actual solution lim n→∞ n (θ) = (θ) , (3.14) when the seed function 0 (θ) is chosen as the driving term 14 The existence and uniqueness of the limit (3.14) has been proven rigorously in [30] for the fermionic single-particle 15 TBA equation (1.12) with a kernel satisfying the requirement The fermionic 1CDD models do satisfy this condition and, as such, the iteration procedure is guaranteed to converge nicely in the whole range R ∈ R >0 , a fact which is easily verified numerically. All the other models we considered above, on the other hand, violate one or more of the hypotheses of the existence and uniqueness theorem in [30] -being either of bosonic statistic, or having a kernel with L 1 measure ||ϕ|| 1 = 2, or both -and are not guaranteed to possess a convergent iterative solution. Notice that the L 1 measure of the TBA kernel (3.6) counts the number of CDD factors meaning that, in the class of models described by the S-matrix (3.1), only the subfamily with (σ, N ) = (−1, 1) is guaranteed to have a convergent iterative solution. 13 In some limiting cases, it is possible to derive exact expressions, e.g. for the ground-state energy in the conformal limit, via the so-called "dilogarithm trick", as explained nicely in [29].
14 In the case in which the iterative procedure does converge, there is actually a vast freedom in the choice of the seed function. However the standard choice indicated in the main text is the most natural one. 15 See also [31] for an extension to fermionic multi-particle TBA equations.  Figure 4: Ground-state energies for the various models discussed above, along with that of the TT -deformed free fermion (black dots). The empty (resp. filled) markers correspond to models with bosonic (resp. fermionic) statistics. The fermionic sinh-Gordon and staircase models can be solved iteratively all the way to the R → 0 limit, while the rest fail to converge below a certain model-specific scale R * . The parameters of the models were chosen as to allow a comfortable visual comparison between the curves and are the same for both bosonic and fermionic versions of the same model. Insets: inverse square of the (numerical) derivative. As shown by the fits (dotted lines), the fermionic sinh-Gordon and staircase models show the conventional UV behavior ∝ R 4 , while the other models develop a ∝ R behavior reminiscent of the square-root branching singularity of the ground state energy.
We investigated numerically the 1CDD models (a) and (b) and the 2CDD models (b) to (d) 16 , for both the bosonic and fermionic statistic, using the iterative procedure (3.13). As already mentioned above we observed that only for the 1CDD fermionic models this procedure converges for all positive values of the radius R. In every other case, there exists a positive "critical radius" R * > 0 such that for R ≤ R * the iterative routine stops converging. As R approaches R * from larger values, we noticed that the rate of convergence of the iterative numerical routine slows down dramatically, a telltale sign of the existence of some kind of singularity nearby 17 . In Figure 4 we collected the plots of the ground-state energy E(R) for one representative point in the parameter space for each of the models we mentioned above along with one for the TTdeformed free fermion. The shape of the curves suggests that all the cases, apart from the fermionic 1CDD models, behave qualitatively in the same way as the TT -deformed free fermion, that is to say they develop a square-root type singularity at some critical value of the radius R = R * > 0: In order to further confirm this suspicion we plotted the derivative of the ground-state energy to the power −2 in the vicinity of the supposed critical point. As we can see in the insets of Figure 4, the numerical results are in good accord with the hypothesis that R * is a singular point of square root type, as expressed by (3.18).

Two branches
Having our expectation confirmed leaves us with the question of how to deal numerically with such a square root critical point. In particular, the behavior (3.18) implies the existence of a secondary branch of the ground-state energy, behaving as 19) in the vicinity of the critical point. Here and below we are going to use the notatioñ E(R) for the secondary branch. We would like to be able to access numerically to this secondary branch and to explore its properties, e.g. its large R behavior and the possible existence of further critical points. The iterative routine (3.13) is ill suited for this job and we need to employ a more refined method, the PALC mentioned in the introduction and described in §4. Deferring a more thorough analysis of the properties of E(R) to §5, let us present here its main qualitative features, concentrating on a single point in the parameter space of the fermionic 2CDD model (d) as a representative case. More specifically let us set θ 0 = 1/2 and γ = 3π/20 and compute numerically the ground-state energy of the model defined by the S-matrix (3.11). The result is displayed in Figure 5. We see that the function E(R) does indeed possess two branches with distinctly different IR behavior. The primary branch is characterized by the universal IR behavior where K 1 stands for the modified Bessel function while the secondary branch approaches a linear behavior at large RẼ with a rate of approach likely to be some negative power of R. For the specific case depicted in Figure 5 the coefficient of the linear term is found to be while the constant term is vanishing up to the precision we used for our numerical routines. We will see in §5 that this is the asymptotic behavior predicted by analytical considerations. In the zoomed box in Figure 5 we also plotted a fit of the function E(R) in the vicinity of the critical point R * . As expected the behavior in this region is best described by the square-root function ( Another notable fact is that we see no trace of additional singular points: the PALC method can, apparently, reach arbitrarily large values of R on the secondary branch and the resulting ground-state energy quickly approaches the expected asymptotic linear behavior. We note again that the behavior of E(R) depicted in Figure 5 is qualitatively identical to the one exhibited by the ground-state energy of TT -deformed models for negative values of the deformation parameter α, as described in §2 (see e.g. Figure 2).
Finally, we stress that the features of E(R) described here for a point in the parameter space of a specific model really are representative of the general behavior of the ground-state energy in the family of models defined by the S-matrices (3.1), at least for what concerns the case of fermionic statistics. As we will discuss in §5 the status of the models with bosonic statistics is still not completely settled. In particular it is still unclear whether the secondary branch of E(R) displays additional critical points or continues undisturbed in the deep IR and, if this was the case, what type of behavior it follows.

Numerical Method
The results displayed in the previous section suggest that the solution to the TBA equation (1.12), for S-matrices of the form (3.1), may generally possesses a singular dependence on the parameter R. In particular the slope of the tangent to the graph of E(R) apparently diverges at some R = R * . Such critical points are known as turning points. Their presence in the dependence of the ground-state energy E on the system size R evokes the case of the TTbar deformed models, in which all the quantities obtainable from the TBA display a square-root singularity at the same value R = R * .
The iterative procedure described in §3.1 becomes unstable at R → R * , therefore it is not particularly suitable for analyzing the vicinity of the singular point. Fortunately, many powerful methods exist that are capable of handling numerically critical points in non-linear equations. We refer to the nice monograph by Allgower and Georg [32] for an introduction, paired with an extensive literature, to the subject. The simplest of these numerical routines is the already mentioned PALC method which, in spite of the simplicity of its implementation, will be entirely sufficient to handle the situations of interest for us. In this section we will quickly review this method and its main features.

The pseudo-arc-length continuation method
Before starting let us point out a trivial fact: the TBA equation (1.12) is non-linear. It is then not at all surprising that its solutions can develop a highly non-trivial dependence on the parameters. Conversely, what is remarkable is that in the vast majority of instances known in the literature, the solution to the TBA equations display a simple behavior as functions of R. In full generality, we should expect a solution (θ|R) to potentially present, as a function of R 18 , any type of critical point imaginable. As we will see later, in the cases of the 1CDD and 2CDD models we are concerned with here, only turning points appear. We will thus restrict our attention to the simple cases in which every critical point is a turning point. This considerably simplifies both the discussion and the actual implementation of the PALC method, although, if needed, it is entirely possible -and not exceedingly difficult -to include the existence of bifurcations in the game.
Since our goal is to analyze the TBA equation (1.12) numerically, we are going to describe the principles of the PALC for maps between finite-dimensional spaces. Let us then truncate and discretize the real θ-line on a N -point lattice {θ k | k = 1, 2, . . . , N } which, for the moment, we are not going to specify further. Now, consider a parametrized map H which takes as input a parameter R ∈ R together with the values k = (θ k ) ∈ R of some real function on the lattice, and yields N real numbers: We call this the solution curve for the map H.
Our goal is to follow the solution curve from a given starting point C i = ( i , R i ) to a final one C f = ( f , R f ). The most straightforward way to achieve this is to simply parametrize the curve by R and employ some numerical iterative routine, such as the one reviewed in §3.1, to move from C i = C(R i ) to C f = C(R f ). However this simple-minded approach fails at any point in the parameter space where the rank of the Jacobian is not maximal. There we can no longer rely on the implicit function theorem to solve (4.2) for in terms of R. More geometrically, what happens is that the curve C(R) displays a turning point, where d dR C(R) diverges. Fortunately there exists a very simple cure for this problem: instead of parameterizing the curve C by the parameter R, we can use an auxiliary quantity s, traditionally chosen to be the arc-length of C or a suitable numerical equivalent, whence the name pseudo-arc-length given to this approach. The condition (4.2) then becomes In order to proceed, let us take a derivative of this condition with respect to the parameter s. We immediately obtain where the extended Jacobian is a N × (N + 1) block matrix, whilė is an (N + 1) column vector. At this point we seem to be short of 1 condition, since we introduced an additional parameter. However, remember that we decided to choose s as the (pseudo-)arc-length of C, which means ||Ċ(s)|| = 1 . capable of dealing with the presence of turning points. Still, this formulation is somewhat unnatural as it completely disregards the fact that the curve C is the fixed point of the map H, and, as such, should enjoy powerful local contractive properties with respect to iterative solution methods -such as Newton's method. We are then led to an integrated approach in which we numerically integrate (4.11) very coarsely and subsequently employ some kind of iterative method to solve (4.6) locally. This is the general strategy behind the approaches known as predictor-corrector routines. In Appendix A we are going to describe the one that we employed in this work and present a pseudo-code of its implementation.

Results for the 2CDD model
Here we present some results obtained using the numerical techniques of the previous Section. We first concentrate on the fermionic 2CDD models and then discuss some facts about the bosonic models.

Fermionic case
The numerical data we collected, of which we have shown some example in §3.2, strongly indicate the following properties of the ground-state energy E(R) as a function of R: We could not find a convincing analytic argument proving the first three properties and we regard them as experimental observations. On the other hand, the last property (3.21) can be verified analytically, as we are now going to show.

The large R behavior
Let us analyze the possible behaviors of the TBA equation (1.12) at large R. To this end, we write the equation as follows where d(θ) is the driving term and χ(θ) the convolution: As R → ∞, the driving term becomes large, ∼ R, and, in order for the equation (5.1) to be satisfied, it has to be balanced by a similar behavior in either (θ), χ(θ) or both. The standard assumption is that which turns out to be consistent, since, as one easily verifies, However this is not, in general, the only possibility. It might be the case that the convolution term χ(θ) is diverging as R → ∞ and becomes comparable with either (θ), d(θ) or both. It is then not difficult to check that only two possibilities are consistent: 1. (θ) −→ R→∞ 0 and the kernel ϕ(θ) is not integrable on the real line; is positive only in some finite 19 subset Θ ⊂ R of the real line and negative everywhere else. 19 The subset Θ cannot be infinite, since the equation (5.1) forces (θ) to behave as d(θ) for θ → ±∞.
The scenario 1 cannot arise for the class of models we are dealing with 20 , since the kernels (3.6) are obviously bounded functions of θ ∈ R. The situation 2 is, on the other hand, a possible one. Let us explore its consequences.
In the hypothesis that the convolution can be approximated as follows Discarding the second term in the right-hand side, we arrive at the linear equation Due to our hypothesis on the function f (θ), we see that the integrand in the right-hand side above is positive for any (θ, θ ) ∈ R 2 , which implies the following bound , then the following inequalities are Rearranging the right inequality above, we find that which we can interpret as a constraint on the class of models which allow for this scenario. In fact, remember that the integral of the kernel on the whole real line, (3.17), counts the number N of CDD factors appearing in the S-matrix (3.1). But, since we assumed that Θ is a finite subset of R, we find that Thus we have found that the fermionic 1CDD models, namely sinh-Gordon and the staicase models, can only display the standard large R behavior (5.3), (5.4). We stress that this result should not be read as a proof of the absence of turning points in these models, but rather as a sanity check for the correctness of our computations, since the ground-state energy for fermionic 1CDD models is well known to be a smooth and monotonously increasing function of the radius in the whole range R > 0. Conversely, all fermionic N CDD models with N > 1 allow for both the standard large R behavior (5.3), (5.4) and the non-standard one (5.5). Consequently, their ground-state energy will possibly display both the asymptotic behavior (3.20) and (3.21), where (5.12) in accordance with the numerical data we have obtained.

Analysis of the numerical data
The fermionic 2CDD models were classified in §3 into cases (a) to (d). We have performed numerical analysis for all the different cases and the results show that the behaviors are qualitatively the same. Thus, we are going to show here the details of the numerical analysis only for the representative case (d). We begin by analyzing the numerical solution obtained through the PALC method for large values of R. It was argued in the previous section that the pseudoenergy should behave as in (5.5), assuming negative values in a finite subset of the real line and positive values elsewhere. This is indeed checked to be true for all the 2CDD models under consideration, as illustrated for a particular member of this family in Figure 6, and to be contrasted with the standard iterative solution (the primary branch) which is positive everywhere. The numerics indicate that the negativity region is always a single interval centered at the origin of the form Θ = {θ ∈ R | − Λ ≤ θ ≤ Λ}. They also indicate that the interval size Λ is model-dependent. In particular, it seems to grow with θ 0 and decreases with γ. Nevertheless, the precise dependence of Λ on the parameters deserves further investigation.
We then proceed to analyze the secondary branch solution in the opposite extremum of R, i.e., as R approaches the critical value R * . For some of the plots it will be that alleviates the exponential dependence (with x * = log(R * /2) for the corresponding critical point). Here we find it more instructive to display L(θ) instead of the pseudoenergy itself in order to ease the comparison with the primary branch solution. The situation is illustrated in Figure 7. The two branches approach each other as the value of R decreases, eventually merging at R = R * after which they become complex-valued. For each R, the function L(θ) for the secondary branch is everywhere larger than the corresponding primary branch counterpart, which is compatible with the previously mentioned fact that it has lower energy (recall the overall minus sign in (1.14)). The critical value R * could in principle have a dependence on θ. We ran an extensive numerical test exploring this possibility, but all the numerical results indicate θ-independence to high accuracy, even though at this moment we do not have an analytic proof of this property. The analyses were as follows. We first ran the iterative numerical routine and computed the pseudoenergy (θ) for at least ten different values of x differing from each other and from x * by 10 −8 . Then, we selected several values of θ and for each value we performed a square root fit of the form a(θ)+b(θ) −x * (θ) + x. The fits were done using Mathematica's NonlinearModelFit function by giving an initial estimate for x * (θ). By comparing all the obtained x * (θ), we verified that they agree up to errors greater than 10 −8 which was our minimal working precision. The analysis was performed for several values of θ 0 and for γ in the range 0 ≤ γ ≤ (99/200)π. In many cases, when the number of necessary points in the discretized θ grid was not very high it was possible to work with even higher precision. In those cases, another way of getting x * with high precision is by assuming a square root behavior for the pseudoenergy and solving the resulting equations using Mathematica's FindRoot function.
In addition, we also verified that R * depends smoothly on the model parameters θ 0 and γ, as shown in Figure 8 for both the fermionic and bosonic models. In particular, for large θ 0 we have the asymptotic behavior x * = log(R * /2) ≈ −θ 0 +x (0) * (see §5.3 for a derivation in the special limit where γ is close to π/2, for which x (0) * = log log(2 + 2 √ 2); for other values of γ the linear term remains the same, though x (0) * is different).

Bosonic case
We have also repeated the analysis described above using the PALC method to the case of bosonic systems. The numerical routine used in this case only differs from the fermionic case by a few signs. As already mentioned in §3, the solutions to the TBA equation for the bosonic models have intricate behavior already for the 1CDD cases. It was first noticed in [19] (for the case of real u 1 , in the notation of (3.7)) that the numerical iterative routine stops converging for some R * , signaling the presence of a singularity. In fact, we have verified numerically that all the bosonic models up to two Figure 8: Dependence of the critical x * in the model parameters. Black lines correspond to fermionic 2CDD models, red lines correspond to bosonic ones. On (a), we demonstrate the validity of the narrow resonance limit approximation for x * (red and black bullets/boxes), see in 5.3.
CDD factors behave similarly to the fermionic 2CDD models of previous section, i.e., they have a "primary branch" and a "secondary branch" which merge at a critical scale R * , where the energies E(R) have square-root singularities in R * , and the value of R * is independent of θ.
There is a simple argument based on the well-known relation between bosonic and fermionic TBA which makes this behavior of the bosonic 1CDD model rather natural. (see (3.6)) 21 . Recalling the arguments presented in §5.1.1, we conclude that bosonic N CDD models admit two different types of large R behaviors whenever N > 0. The large R regime of the pseudoenergy (θ) for the primary branch is as expected and it is easily accessed numerically, however for the secondary branch it is more involved to compute it. By increasing the value of R, eventually we reach a value R where the PALC method suddenly ceases to provide a real solution and reverts back to the primary branch solution. Analyzing the behavior of (θ) for complex values of θ, we verified that a pair of complex conjugate zeros of z(θ) = 1 − e − (θ) is approaching the real axis and causing the numerical instability. In principle it is possible to refine the numerical methods so as to obtain solutions for R > R . However, it is not clear at the moment whether or not those singularities of L(θ) ever cross the real axis. In case they do, an analysis similar to the one performed in [33] for the excited state TBA could be carried out. We leave the analysis of the large R behavior of the secondary branch in bosonic models for a future study.
The behavior of the models for R close to R * is illustrated in Figure 9 by the L(θ) function for a 2CDD model of type (d). The qualitative picture is similar to the fermionic case, i.e. the function L(θ) for the secondary branch solution is greater everywhere than the one for the primary branch and the two merge as the critical point is approached. We conclude this subsection by showing in Figure 8 the smooth dependence of x * on the model parameters and in particular in the limit γ → π/2. In addition, notice that the bosonic curve is always above of the fermionic curve for the same parameters. This can be understood by analyzing the map (5.15) and the fact that the additional delta function term always give a positive contribution to the convolution term of the TBA equations.

Narrow resonance limit
Here we consider the special limit γ → π 2 of the kernel (3.12). In this limit the poles of the kernel get closer to the real line, finally forming two Dirac δ functions. We shall refer to this as the Narrow Resonance (NR) limit. After integration of the delta functions and exponentiation, TBA (1.12) becomes the difference equation where we introduced the notation Y (θ|R) = e − (θ|R) . Note that this can be seen as an infinite set of equations relating the values of Y on the grid points θ ∈ (−θ 0 , θ 0 ) + θ 0 Z. Let us focus on the fermionic case (σ = −1). Introducing y k = Y (θ + kθ 0 ) and g k = e −R cosh(θ+kθ 0 ) we can write (5.16) as and look for a solution for different grids specified by a choice of θ. This is an infinite set of equations, however starting with k = 0 one can obtain an approximate solution by truncating the system for some |k| ≤ m, since g k and y k decay very rapidly with R and θ 0 , and hence with k. Truncating to m = 1 leads to the quadratic equation One can now choose the integer lattice (i.e., θ = 0), to get The solution develops a square root singularity at x * ≈ −θ 0 + log log(2(1 + √ 2)), which is compatible with our findings in §5.1.2. This point is shown as a red bullet in Figure  8a. In contrast to the general case, it is clear that here the branching point depends on the choice of θ lattice. Let us also comment that a similar analysis of the truncated system in the bosonic case (σ = +1) using the half-integer lattice leads to the black box shown in Figure 8a Note that the truncation to m = 1 is only valid for sufficiently large R and θ 0 . Increasing the truncation order leads to more coupled equations, which in turn can be recast as an (more complicated) algebraic equation for y 0 , with parameters depending on θ. The number of solutions increases accordingly. However, for any θ ∈ (−θ 0 , θ 0 ) there is always a single pair of solutions which collide and form a branching point at some x * (θ) ≈ −θ 0 + const., corresponding to real, positive R * (θ), a feature that is not altered by increasing the truncation order.
Finally, we remark that in the further special limit θ 0 = 0, the difference equations (5.16) become simple algebraic equations for Y (θ) that can be exactly solved both in the fermionic and the bosonic case, leading to exact expressions x * = log log 2 and x * = log log 3 2 √ 3 for the critical points, respectively. These points are also shown in Figure 8a, emphasizing the smooth nature of the limit γ → π 2 . In Figure 10 we present as an example a solution with m = 8 truncation together with the iterative solution of the integral equation (1.12) for θ 0 = 2 and γ approaching π/2, just before reaching the (first) critical R * (θ) of the NR limit. The transition seems to be smooth, however we do not yet have a complete understanding of the nature of this limit. We plan to revisit the narrow resonance model in a more sistematic way in the future.

Discussion
There are two general questions which we believe our results shed some light upon. One concerns the short-distance behavior of the theory under the generalized TTbar deformation (1.4). Our results supports the expectation that, at least in the cases when the CDD factor in the associated S-matrix deformation has the form (1.10), (1.8) with finite N , the theory develops the Hagedorn singularity corresponding to a density of high-energy states much greater than what is allowed in a Wilsonian QFT. Although we demonstrated this in a limited set of examples -the 2CDD deformations of the free S-matrix with both fermionic and bosonic statistics and the 1CDD deformations of the free boson S-matrix -this result likely extends to more general N CDD deformations, at least for massive theories involving only one kind of particles. In fact the case N = ∞, a model known as Elliptic sinh-Gordon, is shown to display the same behavior as the ones studied here [34]. We note that this behavior is qualitatively the same as the one encountered under the "TTbar proper" deformation (1.1) of a generic local QFT. Moreover, the singularity of E(R) at the Hagedorn point R * is a square-root branching point, exactly as in the TTbar deformations with negative α. From a formal point of view, this nature of the singularity is not entirely unexpected. Indeed, the character of the singularity relates to the rate of approach of the Hagedorn asymptotic (2.9) at high energy E → ∞. Assume that the approach is power-like 23 where κ is some positive number, L is the spatial size of the system which is assumed to be asymptotically large, and the dots represent yet higher negative powers of E. The dependence on L of the subleading term reflects the extensive nature of the entropy, which must behave as L σ(E/L) in the limit L → ∞, with the intensive quantity -the entropy density σ -depending on the energy density E/L. Inspection of (6.1) reveals the mass dimension of the coefficient a to be a ∼ [mass] 2κ+1 . Having in mind that all the deformation parameters α s in (1.4) have even integer dimensions, one could expect that the exponent 2κ+1 is an integer. The lowest positive κ consistent with this assumption is κ = 1, and then (6.1) leads exactly to the square-root singularity of E(R). Still, the physics behind this simple character of the singularity appears mysterious. Analytic continuation of E(R) below R * returns complex values of E. This likely signals an instability of the ground state at R < R * against some sort of decay. If so, what is the product(s) of the decay? Usually in a theory with finite range of interaction the decay of the unstable ground state goes through the process of nucleation, as in the "false vacuum" decay studied in [36,37]. However such a decay would imply a much weakerand analytically more complicated -singularity at R * . Therefore the simple algebraic character of the actual singularity appears puzzling. A different, but possibly related, question is the physical interpretation of the secondary branch of E(R) discovered in §5.
An even more general question concerns the relation between the S-matrix and the underlying local structure. Suppose we are given an S-matrix, i.e. a collection of masses of stable particles as well as the full set of scattering amplitudes, satisfying all the standard requirements of the S-matrix theory -unitarity, analyticity, crossing and bootstrap conditions (see e.g. [38,39]), with the singularity structure consistent with the macro-causality [40]. Is there a local QFT generating such a scattering theory? The answer is generally no. There are consistent S-matrices which cannot be derived from Wilsonian QFT, and indeed do not have an underlying local structure, meaning a complete algebra of local operators. This possibility is famously realized in string theories. The results presented here support the expectation that the overwhelming majority of self-consistent S-matrices are not derivable from local QFT. Although this expectation arise from a general analysis of the RG flows [4], we substantiate it by providing concrete examples in 1+1 dimensions with factorizable S-matrices consisting of pure CDD factors. We studied a number of examples of such S-matrices and verified that they lead to the Hagedorn density of high-energy states (2.9), familiar to the string theories. What's more, it looks likely that this situation is rather general: with the exception of a small subset of "local field-theoretic" S-matrices, the bulk part of the space of consistent, factorizable S-matrices in 1+1 dimensions, leads to a Hagedorn transition. This statement of course needs to be verified on a more systematic level, but it is tempting to conjecture that this is a general situation, not limited to integrable theories and to low space-time dimensions. If so, would it mean that the majority of consistent S-matrices correspond to some kind of string theories? Or maybe there is a more general class of theories, besides the strings, which break the standard local structure of QFT while preserving macro-causality and exhibit the stringy density of states?
The present work represents a first step of a project having as a goal the systematic analysis of the TBA equations for completely general CDD-deformed factorizable Smatrices, with arbitrarily complicated CDD factors (1.7), possibly including the factors (1.9) with singular behavior at high energies. Clearly, also CDD deformations of more complicated S-matrices, involving more than one kind of particles -possibly having mass degeneracies, a situation leading to off-diagonal scattering -have to be studied. Such S-matrices lead to systems of TBA equations more complicated than the simple equation (1.12). Nonetheless, we believe that the numerical methods adopted here, in particular the PALC routine, can be adopted in full generality. Finally, a similar analysis can be extended to the CDD deformed "massless TBA systems" (see e.g. [41,42,29]). Although the physical foundation here is less firm -since the notion of S-matrix is ambiguous for massless theories in 1+1 dimensions -these cases might yield welcome surprises.
• Predictor. This part of the routine takes as input a point c(s j ) = ( j , R j ) on the solution curve and uses the initial value problem form (4.11), which we recall here H (c(s))ċ(s) = 0 , ||ċ(s)|| = 1 , c(s j ) = ( j , R j ) , (A.5) to yield a reasonable guess for a new point c (0) (s j+1 ) = ( ). The simplest way to obtain such a point is to employ the so-called Euler predictor, which implements the equation where the N + 1 vector t j is tangent to the extended Jacobian H (c(s)) at the point ( j , R j ): • Corrector. This second part of the routine engages in the problem of adjusting the predictor's output ( j+1 ) to a point actually lying on the solution curve. It does so by some iterative method for solving the equation H = 0 starting from an initial, reasonably close, guess. The fastest and least expensive of these methods is the Newton's one, which in our case would take the following form if only we were not worried to encounter a point where J is not invertible. In fact we are concerned precisely with such an eventuality, it being the very reason that led us to consider the PALC method and the associated predictor-corrector routine. Hence, we need to appropriately modify Newton's method in order to accommodate the possibility of a singular J , with H of maximal rank N . The way to handle such a situation is to consider the concept of quasi-inverse (also called Moore-Penrose inverse) A + of a matrix A, defined as where a superscript T denotes standard matrix transposition. Notice that, if A is a square matrix, the above definition is equivalent to the standard inverse. Now, if A is instead an N × (N + 1) matrix of maximal rank N and t is its tangent vector At = 0, then the following statements are equivalent 1. Ax = b and t T x = 0, 3. x = min v ||v|| Av = b which, in plain words, means that x is the vector of minimal norm which solves the equation Ax = b.
Without going too much in the details (see chapter 3 of [32]), the takeaway is that we can implement Newton's method in the usual way, as long as we trade the inverse of the Jacobian for the quasi-inverse of the extended Jacobian: Here follows a pseudo-code summarizing the procedure expounded above. As we can immediately see, the algorithm requires an initial point solving the TBA equation. This can be provided by using the standard iterative procedure of §3.1 to solve the equation at some value of R > R * . This will yield a solution ( 0 , R 0 ) on the first branch, from which to start the PALC.