# Thermodynamic Bethe Ansatz past turning points: the (elliptic) sinh-Gordon model

Lucía Córdova, Stefano Negro, Fidel I. Schaposnik Massolo

► To cite this version:

Lucía Córdova, Stefano Negro, Fidel I. Schaposnik Massolo. Thermodynamic Bethe Ansatz past turning points: the (elliptic) sinh-Gordon model. Journal of High Energy Physics, 2022, 2022 (01), pp.035. 10.1007/JHEP01(2022)035. hal-03424946

HAL Id: hal-03424946

https://hal.science/hal-03424946v1

Submitted on 10 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Thermodynamic Bethe Ansatz past turning points: the (elliptic) sinh-Gordon model

Lucía Córdova,ª Stefano Negroᵇ and Fidel I. Schaposnik Massoloᶜ

E-mail: lucialgc@gmail.com, stefano.negro@nyu.edu, fidels@ihes.fr

**ABSTRACT:** We analyze the Thermodynamic Bethe Ansatz (TBA) for various integrable S-matrices in the context of generalized $T\bar{T}$ deformations. We focus on the sinh-Gordon model and its elliptic deformation in both its fermionic and bosonic realizations. We confirm that the determining factor for a turning point in the TBA, interpreted as a finite Hagedorn temperature, is the difference between the number of bound states and resonances in the theory. Implementing the numerical pseudo-arclength continuation method, we are able to follow the solutions to the TBA equations past the turning point all the way to the ultraviolet regime. We find that for any number $k$ of resonances the pair of complex conjugate solutions below the turning point is such that the effective central charge is minimized. As $k \to \infty$ the UV effective central charge goes to zero as in the elliptic sinh-Gordon model. Finally we uncover a new family of UV complete integrable theories defined by the bosonic counterparts of the $S$-matrices describing the $\Phi_{1,3}$ integrable deformation of non-unitary minimal models $M_{2,2n+3}$.

**KEYWORDS:** Integrable Field Theories, Conformal Field Theory

ARXIV ePRINT: 2110.14666

<footnote>ªInstitut de Physique Théorique Philippe Meyer,
Laboratoire de Physique de l'École Normale Supérieure, PSL University,
CNRS, Sorbonne Universités, UPMC Univ. Paris 06,
24 rue Lhomond, 75231 Paris Cedex 05, France</footnote>

<footnote>ᵇCenter for Cosmology and Particle Physics, New York University,
New York, NY 10003, U.S.A.</footnote>

<footnote>ᶜInstitut des Hautes Études Scientifiques,
35 route de Chartres, 91440 Bures-sur-Yvette, France</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Elliptic sinh-Gordon models and the thermodynamic Bethe ansatz</td><td>5</td></tr><tr><td>3</td><td>Numerical methods</td><td>7</td></tr><tr><td>4</td><td>Models with resonances: turning points in the TBA</td><td>10</td></tr><tr><td>4.1</td><td>Real branches and asymptotic properties</td><td>10</td></tr><tr><td>4.2</td><td>Fermionic/bosonic sinh-Gordon</td><td>12</td></tr><tr><td>4.3</td><td>Fermionic/bosonic elliptic sinh-Gordon</td><td>16</td></tr><tr><td>5</td><td>Bound state explorations: deformed minimal models</td><td>21</td></tr><tr><td>6</td><td>Discussion</td><td>24</td></tr><tr><td>A</td><td>Some properties of elliptic functions</td><td>26</td></tr><tr><td>B</td><td>Approximate algebraic relation between kernel and convolution</td><td>27</td></tr><tr><td>C</td><td>Pseudo-arclength continuation of the TBA equation</td><td>28</td></tr></table>

## 1 Introduction

Quantum Field Theory (QFT) provides a framework capable of describing a wealth of different phenomena with astonishing accuracy. As such, it constitutes one of the pillars sustaining our modern understanding of Nature. It is hard to overstate the importance of obtaining a complete and satisfactory understanding of QFTs. These can be thought as populating a vast landscape, the *theory space*, which can be explored with the help of methods such as the Renormalization Group (RG). Conformal Field Theories (CFTs) play a pre-eminent role in this description as fixed points<fnref n="1" /> of the RG flows. From these, one can then explore the theory space by perturbing the corresponding CFT with a given set of operators. Depending on their scaling dimension, the deforming operators and the associated RG flows are arranged in three distinctive classes: *relevant*, *marginal* and *irrelevant*. While considerable work has been devoted to the first two, the irrelevant deformations have received much less attention and to this date their high-energy physics constitute a largely unexplored territory. There is a very good reason for this. General arguments of renormalization show that the perturbative analysis of irrelevant deformations leads to an

<footnote n="1">¹More precisely, the fixed points of the RG flows are scale invariant theories. For a moderately recent account on the distinction between scale and conformally invariant theories, see [1].</footnote>
accumulation of UV divergences which necessarily require an infinite number of counter-terms, destroying the predictive power of the theory. Ultimately the high-energy behavior of such theories is completely ambiguous and therefore confined to the realm of effective field theory.

The grim destiny looming over the UV behavior of irrelevant deformations is not inescapable, at least not for all theories. In the last few years, certain special classes of irrelevant deformations in two space-time dimensions were shown to be under exceptionally good control, even deep in the UV. The poster child of these is the so-called T$\bar{T}$ deformation [2, 3], triggered by the composite irrelevant operator built from the components of the energy-momentum tensor [4]. This deformation possesses several remarkable properties. It was shown to be intimately related to two-dimensional gravity [5–8], random geometry [9], string theory [3, 10–14] and holography [15–17]. The T$\bar{T}$ deformed theories are solvable, in the sense that physical observables of interest, such as the finite-volume spectrum, the S-matrix [2, 3, 5] and the partition functions [6, 9], can be determined exactly in terms of the corresponding undeformed quantities. This high degree of control means that it is possible to follow the irrelevant flow stemming from the undeformed CFT all the way to the UV, effectively reversing the renormalisation group trajectory, and to obtain exact results on the UV physics. The latter are remarkable: the finite-size density of states grows exponentially at high energies, in a Hagedorn fashion [18] reminiscent of string theories [19–21]. This exceptional feature is already manifest in the short-scale behavior of the finite-size groundstate energy. If we let R denote the circumference of the circle on which the space component of two-dimensional space-time is compactified, then the groundstate energy $E_α(R)$ of a T$\bar{T}$ deformed system obeys the functional equation [2, 3]

$$E_{\alpha}(R) = E_{\alpha=0}(R - \alpha E_{\alpha}(R)), \qquad (1.1)$$

where $α$ is the deformation parameter. Depending on the sign of $α$, the function $E_α(R)$ determined by (1.1) either presents a square-root singularity at a positive radius $R_* \sim |α|^{-1/2}$ (for $α < 0$) or possesses no short-scale singularity at all (for $α > 0$). Neither of these is compatible with Wilson’s paradigm of local QFTs [22], leading us to the conclusion that T$\bar{T}$ deformed theories cannot be considered conventional UV-complete theories and that, thanks to their robust features, they represent a sensible extension of the Wilsonian notion of a local QFT.

The compelling features of T$\bar{T}$ deformations are in large part shared by larger families of two-dimensional theories that we will refer to collectively as *solvable irrelevant deformations*. Important examples are the Lorentz-breaking J$\bar{T}$ deformations [23], triggered by the composite operator J$\bar{T}$ that can be defined whenever the undeformed theory possesses a conserved holomorphic U(1) current $J_μ$. Another very large and important family of solvable irrelevant deformations is available if one considers integrable QFTs [2]. In this case we have access to an infinite space of conserved currents $T_{\mu\nu}^{(s)}$ labeled by a half-integer spin<fnref n="2" /> $s \in S \subset \mathbb{Z}/2$, and each of these can be used to construct a two-parameter family of

<footnote n="2">²To be more precise, in complex coordinates $(z, \bar{z})$ the components $T_{zz}^{(s)}$ and $T_{\bar{z}\bar{z}}^{(s)}$ have spin $s+1$, while $T_{z\bar{z}}^{(s)}$ and $T_{\bar{z}z}^{(s)}$ have spin $s-1$. The requirement of locality restricts the spins to be either integers or</footnote>
bi-linear irrelevant operators $T\bar{T}^{(s,s')} \equiv \varepsilon^{\mu\rho}\varepsilon^{\nu\sigma}T_{\mu\nu}^{(s)}T_{\rho\sigma}^{(s')}$ (here $\varepsilon^{\mu\nu}$ is the totally antisymmetric Levi-Civita symbol in two dimensions). Integrable QFTs deformed by these operators were shown to preserve their integrability and the corresponding set of local conserved charges [2, 24, 25], allowing the use of powerful exact, non-perturbative methods in their analysis. Hence, this vast family of theories represents the perfect playground on which to study the ultraviolet properties of theories not complying with the Wilsonian QFT paradigm. Of primary importance is the question concerning the nature of the short-scale physics of these models, with particular focus on the status of their local structures and the mechanisms responsible for the appearance of a Hagedorn temperature. Understanding the physics underlying these matters remains the most important open problem in this context, as well as one of the main motivations for the present work.

Amongst the vast set of irrelevant deformations triggered by the operators $T\bar{T}^{(s,s')}$, a sub-family of Lorentz-breaking ones obtained by fixing $s' = 1$ was introduced and studied in [26], while [27] and [28] initiated the analysis of integrable QFTs deformed by linear combinations of the operators $T\bar{T}^{(s)} \equiv T\bar{T}^{(s,s)}$. In this work we will focus our attention on particular instances of the latter type of deformations that we will refer to as generalized $T\bar{T}$ deformations. The most convenient setting to analyze these is provided by the factorised scattering framework [29]. As explained thoroughly in [28], given an integrable QFT with elastic two-body $S$-matrix $S_0(\theta)$, the space of its generalized $T\bar{T}$ deformations stands in one-to-one correspondence with the space of CDD factor deformations of its $S$-matrix,

$$S(\theta) = S_0(\theta)\Phi(\theta). \qquad (1.2)$$

Here the CDD factor [30] $\Phi(\theta)$ is a scalar phase factor, arbitrary up to the request that its presence does not violate the unitarity, crossing symmetry, analyticity and macro-causality conditions of the $S$-matrix [31–33]. Ultimately these requirements force the CDD factor to be a meromorphic function of the rapidity $\theta$ of the form

$$\Phi(\theta) = \Phi_{\text{prod}}(\theta)\Phi_{\text{entire}}(\theta), \qquad (1.3)$$

where the first factor incorporates all the poles at finite $\theta$, whose number $N$ is arbitrary (and possibly infinite)

$$\Phi_{\text{prod}}(\theta) = \prod_{j=1}^{N} \frac{\sinh \theta + \sinh \theta_j}{\sinh \theta - \sinh \theta_j}. \qquad (1.4)$$

The second factor represents an entire function of $\theta$ of the form

$$\Phi_{\text{entire}}(\theta) = \exp \left[ -i \sum_{s \in \mathcal{S}} a_s \sinh(s\theta) \right], \qquad (1.5)$$

with $\mathcal{S}$ being the characteristic set of spins mentioned above. The series in the exponent of (1.5) is assumed to be convergent for any $\theta$. Finally, the condition of macro-causality [33] restricts possible positions of the poles in (1.4) to either the imaginary axis $\text{Re}\theta_j = 0$, or

<footnote>half-integers. The specific subset $\mathcal{S}$ in which $s$ takes values is a characteristic of the integrable QFT under consideration.</footnote>
to the strips Im $\theta_j \in [-\pi, 0]$ mod $2\pi$. Each pole has a standard interpretation in terms of the spectrum of the theory. The purely imaginary ones signal, when their imaginary part is positive, the presence of bound states of mass $2m \cosh(\theta_j/2)$, where $m$ stands for the mass of the lightest excitation in the theory. The poles $\theta_j$ with Re $\theta_j \neq 0$ and $-\pi \le \text{Im} \theta_j < 0$ lie in the so-called *unphysical strip*<fnref n="3" /> and are associated with unstable particles, also known as *resonances*. The real and imaginary parts of their complex mass $2m \cosh(\theta_j/2)$ are identified, respectively, with the center-of-mass energy and the width of the resonance. The last possibility is to have poles lying on the negative imaginary axis. These are usually referred to as *virtual states* [34], and while they do not have a clear interpretation in terms of particles their presence effects an increment of the scattering phase as a function of $\theta$ at low energies.

In this work, we continue the program laid down in [28] by studying a number of models obtained by deforming a trivial $S$-matrix $S_0(\theta) = \pm 1$ with $\Phi_{\text{entire}} = 1$ and a product CDD factor (1.4),

$$S(\theta) = \pm \prod_{j=1}^{N} \frac{\sinh \theta + \sinh \theta_j}{\sinh \theta - \sinh \theta_j}, \qquad (1.6)$$

where the poles $\{\theta_j\}_{j=1}^N$ all lie in the unphysical strip, thus corresponding to resonances. Here the upper (lower) sign corresponds to the fermionic (bosonic) TBA statistics (see §2 of [35] for an explanation of this point). We will be especially interested in a particular exponent of this class with $N = \infty$ resonances: the *elliptic sinh-Gordon model* [36].

Our interest in this particular model is two-fold. On one hand we can view it as the effect of deforming a free theory with a particular combination of generalized $T\bar{T}$ operators as discussed above. On the other hand its periodic nature allows us to make direct contact with the periodic $S$-matrices arising in the recently revived $S$-matrix bootstrap program [37, 38], and more generally to theories with a large number of resonances.

In the $S$-matrix bootstrap program one explores the space of consistent $S$-matrices by imposing general principles like unitarity, crossing symmetry and analyticity. When restricting to two spacetime dimensions, integrable theories naturally appear at the boundary of the allowed space of theories [37–46]. In this way one can rediscover known integrable models, but also obtain integrable $S$-matrices with no Lagrangian realization. A paradigmatic example of the second situation is the $O(N)$ *periodic Yang-Baxter model*. Its integrable $S$-matrix was proposed in [47] and later found sitting at a special point of the *monolith* parametrizing the space of theories with a global $O(N)$ symmetry and no bound states [44]. This model, as the elliptic sinh-Gordon one, has infinite resonances spaced periodically in the rapidity plane (hence its name). The global symmetry makes it however more complicated than the single-particle models analyzed in this article, so we restrict our attention to the simplest periodic $S$-matrices and leave the study of global symmetry analogues for future work.

The plan of this paper is as follows. In section 2 we review the elliptic sinh-Gordon models and the Thermodynamic Bethe Ansatz (TBA), the primary tool we will use to

<footnote n="3">³This corresponds to the region of the complex center-of-mass energy s-plane reached by analytically continuing the scattering amplitude through the two-particle branch cut.</footnote>
study them. Section 3 contains a review of the numerical approaches to the solution of the TBA equation, with a special emphasis on the application of the pseudo-arclength continuation method required to overcome the limitations of the more traditional iterative methods typically used in the literature. Most of our new results are presented in section 4. After considering the standard sinh-Gordon model as a warm-up we tackle its elliptic deformation, in both its fermionic and bosonic flavors, devoting special attention to the UV and IR behavior of the effective central charge. Section 5 initiates the analysis of CDD deformations of theories with more than a single stable particle by looking at the simple examples provided by the S-matrices of the $\Phi_{1,3}$ integrable deformation of the non-unitary minimal models and their bosonic counterparts. We conclude in section 6 with a discussion of our results and an outlook for future research directions. Three appendices provide some background on elliptic functions, an interesting albeit somewhat mysterious relation between the TBA kernel and the convolution term in the elliptic models, and a more detailed description of our implementation of the pseudo-arclength continuation method.

## 2 Elliptic sinh-Gordon models and the thermodynamic Bethe ansatz

The main model of interest in the present work is an elliptic deformation of the sinh-Gordon model. Its S-matrix was proposed in [36] and reads

$$S_{a,l}^{\text{f,b}}(\theta) = \pm \frac{\text{sn}_l(2iK_l \theta/\pi) + \text{sn}_l(2K_l a)}{\text{sn}_l(2iK_l \theta/\pi) - \text{sn}_l(2K_l a)}, \quad (2.1)$$

where a is the coupling constant, the modulus l parametrizes the periodicity in real rapidity, sn is the Jacobi sine function and $K_l$ is the complete elliptic integral

$$K_l = \int_0^{\pi/2} \frac{d\phi}{\sqrt{1 - l^2 \sin^2 \phi}}. \quad (2.2)$$

The upper (lower) sign in (2.1) corresponds to the fermionic (bosonic) case. This S-matrix has two periods, the usual one along the imaginary direction arising from crossing symmetry and unitarity $S(\theta) = S(\theta + 2i\pi)$ and an unusual one in the real direction $S(\theta) = S(\theta + T_l)$. The latter period is given in terms of elliptic integrals,

$$T_l = \pi \frac{K \sqrt{1-l^2}}{K_l}. \quad (2.3)$$

The elliptic sinh-Gordon model contains an infinite tower of resonances spaced periodically in the rapidity plane. We see them as complex zeros<fnref n="4" /> in the physical rapidity strip defined by $\text{Im}(\theta) \in (0, \pi)$. The location of the “fundamental” zeros is $\theta_n = ia\pi + nT_l$ with $n \in \mathbb{Z}$ and the arrangement of the remaining zeros and poles outside the physical strip follows from the crossing and unitarity conditions of the S-matrix. Figure 1 shows the analytic structure in the physical strip.

<footnote n="4">⁴Due to the unitarity condition $S(\theta)S(-\theta) = 1$, complex zeroes in the physical strip are in one-to-one correspondence with poles in the unphysical strip. Thus their existence signals the presence of a resonance in the spectrum of the theory.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Analytic structure of the elliptic sinh-Gordon S-matrix in the complex rapidity plane. Simple zeros (poles) are depicted by crosses (dots). The shaded region is the physical rapidity strip $0 \le \operatorname{Im}(\theta) \le \pi$ which contains two crossing symmetric towers of infinite resonances separated in the real direction by the period $T_l$.</figcaption></figure>

In order to ensure we have only zeros inside the physical strip — and, consequently, that the model possesses no bound-state — we need to limit the coupling constant $a$ to take values in the interval $[0, 1]$. Moreover, since the model is invariant under $a \leftrightarrow 1-a$, we will consider $a \in [0, 1/2]$ without any loss of generality. The limit $a \to 0$ leads to free theory, $S_{0,l}^{\text{f,b}} = \mp 1$. At zero modulus $l \to 0$ — equivalent to infinite period $T_l \to \infty$ — we recover the usual sinh-Gordon S-matrix with a single resonance

$$S_{a,0}^{\text{f,b}}(\theta) = S_{\text{shG}}^{\text{f,b}} = \pm \frac{\sinh \theta - i \sin(\pi a)}{\sinh \theta + i \sin(\pi a)}. \qquad (2.4)$$

Even though the S-matrix (2.1) approaches well-known theories in the above limits, its physical nature is still mysterious, e.g. there is no known Lagrangian formulation for the theory. To gain more information about these models and, in particular, about their ultraviolet behaviour, we will use the Thermodynamic Bethe Ansatz (TBA) to compute the groundstate energy and the effective central charge of the theory compactified in a circle of circumference $R$. In the following we give a brief review of the method, for a more comprehensive treatment the reader is referred to [48, 49] and references therein.

For the simplest scenario where we have only one type of stable particle of mass $m$, the TBA equation is

$$\epsilon(\theta) = r \cosh \theta \mp \int_{-\infty}^{\infty} \varphi(\theta - \theta') \ln[1 \pm e^{-\epsilon(\theta')}] \frac{d\theta'}{2\pi}, \qquad (2.5)$$

where $r = mR$, the S-matrix enters through the kernel $\varphi(\theta) = -i \partial_{\theta} \ln S(\theta)$ and one solves for the pseudo-energy $\epsilon(\theta)$. The upper and lower signs refer to the fermionic and bosonic cases, respectively, and in what follows we will sometimes refer for brevity to the L-functions $L(\theta) = \pm \ln[1 \pm e^{-\epsilon(\theta)}]$.

Except for very rare cases such as the standard $T\bar{T}$ deformation, one needs to solve the TBA equations (2.5) numerically. As we will see, the S-matrices studied in this paper require a refined numerical method which we review in section 3. In any case, once a
solution $\epsilon(\theta)$ to (2.5) has been found, the effective central charge and groundstate energy are computed from the expressions

$$ \tilde{c}^{\text{f,b}}(r) = \pm \frac{6}{\pi} \int_{-\infty}^{\infty} r \cosh(\theta) \ln \left[ 1 \pm e^{-\epsilon(\theta)} \right] \frac{d\theta}{2\pi}, \quad E_0(R) = -\frac{\pi}{6} \frac{\tilde{c}(r)}{R}. \quad (2.6) $$

Due to the absence of massless excitations, the infrared limit $r \to \infty$ consists of a trivial theory. Accordingly, the effective central charge tends to zero in a precise and universal fashion, controlled by the modified Bessel function $K_1$,

$$ \tilde{c}(r) \stackrel{r \to \infty}{\sim} \frac{6r}{\pi^2} K_1(r) \stackrel{r \to \infty}{\sim} 3 \sqrt{\frac{2r}{\pi^2}} e^{-r}. \quad (2.7) $$

In the usual situation in which the theory under consideration is UV complete, one can relate the central charge $c_{\text{UV}}$ of the CFT controlling the UV behavior to the effective central charge in the $r \to 0$ limit

$$ \lim_{r \to 0} \tilde{c}(r) = c_{\text{UV}} - 12(\Delta_{\min} + \bar{\Delta}_{\min}), \quad (2.8) $$

where $\Delta_{\min}$ is the lowest dimension in the UV theory. The exact value for $\tilde{c}(r = 0)$ is model-dependent but can be obtained fairly easily with the help of the so-called *dilogarithm trick*<fnref n="5" /> [51]. The key point in this calculation is the observation that as $r \to 0$ the pseudo-energy becomes constant in a region $|\theta| < \ln(2/r)$. One therefore proceeds to solve the TBA equation for $r = 0$ and a constant pseudo-energy $\epsilon$,

$$ \epsilon = \mp \Omega \ln (1 \pm e^{-\epsilon}) \quad \text{with} \quad \Omega = \int_{-\infty}^{\infty} \varphi(\theta) \frac{d\theta}{2\pi}. \quad (2.9) $$

The UV effective central charge is finally computed in terms of Rogers' dilogarithm function $L(x) = \text{Li}_2(x) + \frac{1}{2} \ln(x) \ln(1-x)$, as

$$ \tilde{c}(0) = \frac{6}{\pi^2} \times \begin{cases} L\left(\frac{1}{1+e^{\epsilon}}\right) & \text{fermionic} \\ L(e^{-\epsilon}) & \text{bosonic} \end{cases} \quad (2.10) $$

In the following section we will describe the different algorithms that can be used to numerically solve the TBA equation (2.5). These methods, specifically the pseudo-arclength continuation, were implemented to extract the results that will be presented in section 4 and section 5.

## 3 Numerical methods

We would like to solve the TBA equation (2.5) numerically for the elliptic sinh-Gordon models with kernel $\varphi_{a,l}(\theta)$ given by

$$ \varphi_{a,l}(\theta) = -i \partial_{\theta} \ln S_{a,l}(\theta) = \frac{4K_l \operatorname{cn}_l(2iK_l \theta/\pi) \operatorname{dn}_l(2iK_l \theta/\pi) \operatorname{sn}_l(2aK_l)}{\pi \left(\operatorname{sn}_l(2aK_l)^2 - \operatorname{sn}_l(2iK_l \theta/\pi)^2\right)}, \quad (3.1) $$

<footnote n="5">⁵To the best of our knowledge, this first appeared in [50] where it was employed to derive the heat capacity of the Anderson model. For an application of this trick to the TBA of relativistic theories, see [48].</footnote>
where sn$_l$, cn$_l$ and dn$_l$ are respectively the basic Jacobi elliptic sine, cosine and delta amplitude. Note the kernel is identical for the fermionic and bosonic cases, the only difference between these resides in the TBA equation itself.

Already in [48] it was noted that this equation can be solved numerically by an iterative process, whereby we start with an initial approximation $\epsilon^{(0)}(\theta)$ for the pseudo-energy, and produce increasingly more accurate ones through

$$ \epsilon^{(k+1)}(\theta) = r \cosh \theta \mp \int_{-\infty}^{\infty} \varphi_{a,l}(\theta - \theta') \ln (1 \pm e^{-\epsilon^{(k)}(\theta')}) \frac{d\theta'}{2\pi} \quad \text{for } k = 0, 1, \dots \quad (3.2) $$

This process is not guaranteed to converge, but if it does the iterations can be stopped once the desired level of accuracy is reached for the solution.

The only peculiarity in the numerical formulation of our problem is that the kernel $\varphi_{a,l}(\theta)$ has a real period, and therefore the result of the convolution is also periodic with period $T_l$. The pseudo-energy being an even function of $\theta$, for numerical purposes it is thus convenient to discretize an interval in the real-$\theta$ line centered at $\theta = 0$, whose size is an even number of $T_l$ periods. More concretely, we take the following discretization of the interval $[-MT_l, MT_l)$,

$$ \theta_i = -MT_l + \frac{i-1}{N}T_l \quad \text{for} \quad i = 1, 2, \dots, 2NM \quad \text{and} \quad N, M \in \mathbb{N}, \quad (3.3) $$

which then allows us to immediately apply the circular convolution theorem to improve the efficiency of the computation of convolutions. The resulting algorithm is then $O(N \log N + NM)$, since the number of iterations required to converge to a solution is typically $O(1)$.

The procedure described above is well-suited to study families of solutions depending on continuous parameters. Indeed, one may use a solution found at a given point in parameter space as an initial approximation to start the iterative process at another nearby point. For our purposes, this is true both for the coupling constant $a$ and the radius $r$, since the discretization (3.3) does not depend on their values. We will however concentrate mostly on the latter, i.e. we will consider families of solutions related by a continuous variation of the radius $r$.

For large values of $r$ we can neglect the convolution term in the TBA equation and start from the initial approximation $\epsilon^{(0)}(\theta) = r \cosh \theta$. We may then slowly decrease the radius, but in some cases it turns out that the iterative process takes longer to converge, and eventually stops converging altogether, for small enough values of $r$. However, failure to converge is only indicative of a numerical instability and should not be confused with an indication that solutions to (2.5) do not exist beyond what can be found by iteration. Indeed, this numerical instability could in principle be the result of the family of solutions approaching a singular point, but this is not necessarily always the case as will be the case e.g. in section 5.

To address the issue above, we follow$^6$ [28] and implement a pseudo-arclength continuation method [52, 53] on the parameter $r$. In other words, keeping the discretization (3.3)

<footnote>$^6$Our implementation in Mathematica 12.3 is a simplified version of the algorithm described there, relying at its core on the FindRoot function to numerically solve the resulting systems of coupled nonlinear equations.</footnote>
we directly write (2.5) as a set of 2NM coupled non-linear equations on as many variables, $\epsilon_i = \epsilon(\theta_i)$ for $i = 1, 2, \dots, 2NM$. We then introduce an arclength parameter $s$ such that a family of solutions is described by a curve $c(s) = (\epsilon_i(s), r(s))$, and impose the normalization condition $||\frac{dc}{ds}|| = 1$. The whole process can be viewed as a numerical stabilization mechanism for the iterative solution of (2.5), allowing us to follow the solution curve past turning points and bifurcations which may constitute the sources of instability for the more traditional iterative method. We spell out further the pseudo-arclength continuation procedure in appendix C, while the reader is referred to [28] for more details and to [52, 53] and references therein for the theoretical background on this procedure.

**Complex pseudo-energies.** While the TBA equation (2.5) involves a real-valued pseudo-energy $\epsilon : \mathbb{R} \to \mathbb{R}$, we can at least in principle consider the natural generalization to complex-valued pseudo-energies $\epsilon : \mathbb{R} \to \mathbb{C}$. From a purely mathematical point of view this is a well-posed problem, so it may therefore be interesting to explore the space of complex solutions to (2.5), of which the real-valued physical solutions turn out to be a non-trivial subspace.

For this purpose, the numerical algorithms described above require only minimal modifications. Indeed, the iterative method remains unchanged, whereas in the pseudo-arclength continuation method we simply need to double the number of variables to account for both the real and imaginary parts of $\epsilon(\theta_i)$, updating normalization conditions accordingly.

**Bifurcation detection and branch-switching.** When studying families of solutions to the TBA equation related by the variation of the radius $r$ it is important to identify bifurcation points, i.e. points whose immediate vicinity, no matter how small, includes multiple branches of solutions. There is a considerable literature dedicated to this issue in the context of numerical continuation methods [52, 53]. In particular, when performing pseudo-arclength continuation we can detect bifurcation points by keeping track of the augmented Jacobian determinant, which in our case reads

$$J = \begin{vmatrix} \frac{\partial F(\epsilon_i, r)}{\partial \epsilon_i} & \frac{\partial F(\epsilon_i, r)}{\partial r} \\ \frac{d\epsilon_i}{ds} & \frac{dr}{ds} \end{vmatrix}, \qquad (3.4)$$

where $F(\epsilon_i, r) = 0$ is the discretized TBA equation. As we construct a given family of solutions, the quantity $J$ defined above changes signs whenever a bifurcation point is crossed. Therefore, as long as we take small enough steps while performing the continuation we can make sure no bifurcations have inadvertently been skipped. Moreover, once a bifurcation point has been detected we can pinpoint its precise location using binary search to increasingly tighten the bounds on the parameters where the sign-change occurs.

Finally, various techniques have been developed in order to “switch branches” for the purpose of exploring the intersecting family of solutions once a bifurcation point has been identified. In our case, to accomplish this it is sufficient to introduce small perturbations to the initial starting point provided to the underlying Newton-like solver. This allows us to exhaustively explore the space of complex solutions to (2.5) which are connected to the $r \to \infty$ asymptotic solution $\epsilon(\theta) = r \cosh \theta$.
# 4 Models with resonances: turning points in the TBA

## 4.1 Real branches and asymptotic properties

Before presenting the numerical results for the different models, we review some asymptotic properties of the TBA equation (2.5).

The large radius analysis of the fermionic TBA equation in [28] shows that the pseudo-energy admits two possible behaviors for $r \gg 1$: the first one is the standard asymptotic $\epsilon(\theta) \sim r \cosh(\theta)$, whereas the second one is $\epsilon(\theta) \sim -rf(\theta)$, where $f(\theta) > 0$ for $\theta \in \Theta$ and $\Theta$ is a finite subset of the real line.<fnref n="7" /> The latter case is only allowed when the kernel $\varphi(\theta)$ is positive on the real line and satisfies the inequality

$$ \Omega = \int_{-\infty}^{\infty} \varphi(\theta) \frac{d\theta}{2\pi} > 1. \qquad (4.1) $$

The integral above counts the number of resonances minus the number of bound states, so having two more resonances than bound states in the theory will generally allow for a second branch of solutions. For example, a single CDD-zero S-matrix (i.e. having one resonance and zero bound states), such as the sinh-Gordon model, has a single branch of solutions to the TBA equation interpolating between the expected UV and IR behaviours. On the other hand the two CDD-zeros models studied in [28] have both real branches.<fnref n="8" />

Now let us turn to the bosonic TBA equation. The asymptotic analysis follows from the fermionic one through an appropriate redefinition of pseudo-energies and kernel. The bosonic TBA equation with pseudo-energy $\epsilon^+(\theta)$ and kernel $\varphi^+(\theta)$ can be recast into a fermionic TBA equation with pseudo-energy $\epsilon^-(\theta)$ and kernel $\varphi^-(\theta)$ provided we define

$$ \begin{aligned} \epsilon^{-}(\theta) &= \ln[e^{\epsilon^{+}(\theta)} - 1], \\ \varphi^{-}(\theta) &= \varphi^{+}(\theta) + 2\pi\delta(\theta). \end{aligned} \qquad (4.2) $$

The delta function in the kernel above effectively adds a resonance to the counting in (4.1). This interpretation of the bosonic TBA as a fermionic one with an added resonance is supported by the identity

$$ \lim_{u \to 0^{-}} \log \left[ \frac{i \sin u + \sinh \theta}{i \sin u - \sinh \theta} \right] = i\pi \operatorname{sign}(\theta), \qquad (4.3) $$

showing that the $\delta$-function in (4.2) arises as a resonance CDD factor in the limit of the resonance mass going to zero. Consequently, we expect the single-particle bosonic TBA to allow for a second branch of solutions already in the presence of a single CDD-zero in the S-matrix. The leading asymptotic behaviour in the second branch for various quantities of interest is summarized in table 1.

<footnote n="7">⁷A third possibility with $\epsilon(\theta) \sim r^{-1} \cosh(\theta)$ is only viable for $T\bar{T}$ deformations, as opposed to generalized $T\bar{T}$ deformations with general integrable S-matrices.</footnote>

<footnote n="8">⁸Note that this argument does not constrain the local properties of the solutions, only their $r \to \infty$ asymptotic behavior. One might, in theory, observe any number of branch points, of any order, at finite $r$, so long as only two of the branches flow to $r \to \infty$ (or one, in case (4.1) is violated). Numerical evidence coming from various models, however, shows that the simplest possible case of a single square-root branch point (or no branch point if (4.1) is false) is the one presenting itself.</footnote>
<table><caption>**Table 1.** Large radius behaviour of the pseudo-energy, L-function and effective central charge in the second branch of real solutions, to first order in r. The function f(θ) is non-negative for rapidities θ ∈ Θ in a finite subset of the real line, and negative elsewhere.</caption><thead><tr><th></th><th>Fermionic</th><th>Bosonic</th></tr></thead><tbody><tr><td>ε(θ) ~</td><td>-rf(θ)</td><td>e<sup>-rf(θ)</sup></td></tr><tr><td>L(θ) ~</td><td>rf(θ)</td><td>rf(θ)</td></tr><tr><td>c̄(r) ~</td><td>r<sup>2</sup></td><td>r<sup>2</sup></td></tr></tbody></table>

**Higher order corrections.** We are interested in finding the sub-leading corrections to the asymptotic behavior $\epsilon(\theta) \sim -rf(\theta)$. We can write in full generality<fnref n="9" />

$$ \epsilon(\theta|r) = -rf(\theta) + g(\theta|r), \quad \lim_{r \to \infty} \frac{g(\theta|r)}{r} = 0. \qquad (4.4) $$

As mentioned, the function $f(\theta)$ is such that

$$ \left\{ \begin{array}{ll} f(\theta) \ge 0 & \theta \in \Theta, \\ f(\theta) = 0 & \Leftrightarrow \theta \in \partial\Theta, \\ f(\theta) < 0 & \theta \in \Theta^{\perp}, \end{array} \quad \Theta \subset \mathbb{R}. \right. \qquad (4.5) $$

A careful analysis of the TBA equation in the large r limit, with the above conditions on the solution, reveals that $f(\theta)$ obeys the equation

$$ f(\theta) = -\cosh \theta + \int_{\Theta} \varphi(\theta - \theta') f(\theta') \frac{d\theta'}{2\pi}, \qquad (4.6) $$

and that the function $g(\theta|r)$ is expanded in negative odd powers of r,

$$ g(\theta|r) = \sum_{\ell=0}^{\infty} \frac{g_{2\ell+1}(\theta)}{r^{2\ell+1}}. \qquad (4.7) $$

Using this expansion one can then systematically extract from (2.5) the equations for the coefficient functions $g_{2\ell+1}(\theta)$. The first of these reads

$$ g_1(\theta) = \frac{\pi}{12} \left( \frac{\varphi(\theta - \theta')}{f'(\theta')} \right) \Big|_{\theta' \in \partial\Theta} + \int_{\Theta} \varphi(\theta - \theta') g_1(\theta') \frac{d\theta'}{2\pi}. \qquad (4.8) $$

An interesting fact that follows from the above equations is that<fnref n="10" />

$$ \int_{\Theta} \cosh \theta g_1(\theta) d\theta = -\frac{\pi^2}{6} \left( \frac{\cosh \theta}{f'(\theta)} \right) \Big|_{\theta \in \partial \Theta}. \qquad (4.9) $$

<footnote n="9">⁹In the next few formulae we temporarily display explicitly the dependence of all functions on r.</footnote>

<footnote n="10">¹⁰One can see this by introducing the resolvent $K(\theta - \theta')$ of the kernel $\frac{1}{2\pi}\varphi(\theta - \theta')$. Then $f$ and $g_1$ are given by $f(\theta) = -\cosh\theta - \int_{\Theta} d\theta' K(\theta - \theta')\cosh\theta'$ and $g_1(\theta) = \frac{\pi^2}{6} \left( \frac{K(\theta - \theta')}{f'(\theta')} \right) \Big|_{\theta' \in \partial\Theta}$. Finally, since $f(\theta)\Big|_{\theta \in \partial\Theta} = 0$ by definition and $K(-\theta) = K(\theta)$, the identity (4.9) follows.</footnote>
The left-hand side of the above contributes to the constant term in the large $r$ expansion of the effective central charge (2.6). A careful analysis of this expansion shows that another term contributes to this order in $r$,

$$ \int_{\mathbb{R}} d\theta \cosh\theta \log[1 + e^{-|f(\theta)|}] \sim \frac{\pi^2}{6} \left(\frac{\cosh\theta}{f'(\theta)}\right) \Big|_{\theta \in \partial\Theta} \quad (4.10) $$

where we used the Laplace method to compute the right-hand side above. This term exactly cancels the one provided by (4.9), so that the effective central charge displays no constant term in its large-$r$ expansion. Sparing the reader of the uninteresting details, one can compute the contributions to order $r^{-2}$, witnessing a similar, albeit not total, cancellation of terms and arriving at the following expression

$$ \tilde{c}(r) = \frac{3}{\pi^2} r^2 \int_{\Theta} \cosh \theta f(\theta) d\theta + \frac{1}{r^2} \left( \frac{g_1(\theta)}{f'(\theta)} - \frac{7\pi^2}{30} \frac{f''(\theta)}{f'(\theta)^3} \right) \Big|_{\theta \in \partial\Theta} + \mathcal{O}(r^{-4}). \quad (4.11) $$

In principle it should be possible to observe the $\mathcal{O}(r^{-2})$ correction above in our numerical results, but comparison to the analytic expression is hampered by the fact that we do not have an explicit form for $f(\theta)$. We leave this point to be addressed in future work.

## 4.2 Fermionic/bosonic sinh-Gordon

Before considering the more interesting bosonic sinh-Gordon model, we first review some aspects of the well-known fermionic theory defined by the action

$$ \mathcal{A} = \int d^2x \left[ \frac{1}{4\pi} (\partial\phi)^2 + 2\mu \cosh(2b\phi) \right]. \quad (4.12) $$

Here $\mu$ is a coupling constant of dimensions $[mass]^{2+2b^2}$, which determines the scale of the model, while $b$ is a dimensionless parameter that we take, for reasons that will be clarified momentarily, to lie in the interval $b \in [0, 1]$. The action (4.12) is amenable to various interpretations. First, we can understand it as defining a free Gaussian theory deformed by the relevant operator $\cosh(b\phi)$ with negative dimension$^{11}$ $\Delta = -b^2$. Another possibility is to consider (4.12) as a perturbation of Liouville CFT by the operator $e^{-2b\phi}$, again of negative dimension $\Delta = -b^2$. A further possible interpretation, presented in [54], sees the sinh-Gordon model as a conformal affine $sl_2$ Toda field theory with spontaneously broken conformal symmetry. Here we will adopt the first point of view and take sinh-Gordon theory to be a relevant RG flow from a free bosonic theory in the UV to a trivial IR fixed point. As is well-known, this model is integrable, its spectrum consisting of a single neutral particle of mass$^{12}$ $m$ subject to a factorised two-body scattering with amplitude [56]

$$ S_{\text{shG}}^{f}(\theta) = \frac{\sinh \theta - i \sin(\pi a)}{\sinh \theta + i \sin(\pi a)}, \quad (4.13) $$

$$ \pi\mu \frac{\Gamma(b^2)}{\Gamma(1-b^2)} = \left[ \frac{m}{8\sqrt{\pi}} a^a (1-a)^{1-a} \Gamma\left(\frac{a}{2}\right) \Gamma\left(\frac{1-a}{2}\right) \right]^{2+2b^2}. $$

<footnote>$^{11}$The reader might object that a perturbation by an operator of negative dimension destroys unitarity. This is actually not the case for sinh-Gordon theory, due to the presence of fields with non-vanishing vacuum expectation values which alter the two-point function’s UV behaviours from their expected CFT ones. A discussion of this point can be found in [35].</footnote>

<footnote>$^{12}$The particle mass is related to the parameter $\mu$ in the action (4.12) by the famous formula [55]</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2. Effective central charge of the fermionic sinh-Gordon model at various values of the coupling $a$. Numerical solutions were computed using the iterative algorithm described in section 3 with $10^5$ points in the interval $\theta \in [-10, 10)$.</figcaption></figure>

where $a$ and $b$ are related by

$$a = \frac{b^2}{1 + b^2}. \qquad (4.14)$$

Notice that the S-matrix (4.13) is invariant under the weak-strong duality $b \to 1/b$, or equivalently $a \to 1-a$. This is the reason why we restricted the parameter $b$ to the interval $b \in [0, 1]$.

Given the S-matrix (4.13) we can compute the effective central charge via the TBA from (2.6). In this case, the TBA equation can be solved by a simple iterative routine as described in section 3, and figure 2 shows the resulting curves $\tilde{c}(r)$ for different values of the coupling $a$. In particular, we see that there is a single branch of solutions displaying the expected interpolation between the IR and UV behaviours,

$$\lim_{r \to \infty} \tilde{c}(r) = 0, \quad \lim_{r \to 0} \tilde{c}(r) = 1. \qquad (4.15)$$

The above result for the ultraviolet limit of the effective central charge can be obtained following the lines illustrated in section 2. The peculiarity with this case is that the solution is such that $L \to \infty$ and $\epsilon \to -\infty$. The dilogarithm trick can nonetheless still be applied and at the end of the day one sees that

$$\lim_{r \to 0} \tilde{c}(r) = \lim_{\epsilon \to -\infty} \frac{6}{\pi^2} L \left( \frac{1}{1 + e^{\epsilon}} \right) = 1. \qquad (4.16)$$

Sub-leading effects were computed in [57] and [58], and include soft logarithmic corrections of the form

$$\lim_{r \to 0} \tilde{c}(r) = 1 - \frac{3\pi^2 a(1-a)}{2 (\log r)^2} + \mathcal{O}(\log r)^{-3}. \qquad (4.17)$$
<figure><img src="image_4.png" /><figcaption>Figure 3. Effective central charge for the bosonic sinh-Gordon model at various values of the coupling $a$. Above a critical radius $r_c$ there are two branches of real solutions, whereas below this radius there is a pair of complex-conjugate solutions with non-vanishing imaginary parts (dashed lines). The complex effective central charge in the ultraviolet $\tilde{c}(r=0) \approx 0.5 \pm 0.617i$ can be computed from the plateau solutions as explained in the main text. Numerical solutions were computed using the pseudo-arclength continuation method described in section 3 with 500 points in the interval $\theta \in [-6, 6)$.</figcaption></figure>

We now turn to the bosonic sinh-Gordon model defined through the S-matrix

$$S_{\text{shG}}^{\text{b}}(\theta) = - \frac{\sinh \theta - i \sin(\pi a)}{\sinh \theta + i \sin(\pi a)}, \qquad (4.18)$$

which differs from the usual sinh-Gordon one by an overall sign. We call this the bosonic sinh-Gordon model since (4.18) is of bosonic type, i.e. $S_{\text{shG}}^{\text{b}}(0) = 1$. Even though there is no associated Lagrangian or known physical realization for this S-matrix, we take the point of view in which we explore the space of quantum field theories by analyzing consistent S-matrices. In this sense, the S-matrix (4.18) along with the assumption of elastic scattering provides the definition of the theory.

The seemingly harmless change of overall sign in the S-matrix (4.18) has profound consequences for the TBA. As figure 3 shows, there is now a critical radius $r_c$ where we encounter a bifurcation. Above this critical radius we find two branches of real solutions, whereas below it there is a pair of complex conjugate solutions. From the point of view of the real solutions, the critical radius is a turning point in which the two possible branches, i.e. the one with free asymptotic behaviour and the second one with asymptotics described in section 4.1, merge. Let us focus on these solutions first.

The existence at large radius of two, and only two, real branches is compatible with the analysis of [28] for generalized $T\bar{T}$ deformations. That is, for any coupling $a \in (0, 1/2]$ the $L_1$ norm of the kernel is equal to 1, and for the bosonic TBA equation this is enough to
<figure><img src="image_5.png" /><figcaption>**Figure 4.** $L(\theta) = -\ln[1 - e^{-\epsilon(\theta)}]$ functions for the bosonic sinh-Gordon model with $a = 1/2$ at different points along the curve of figure 3. The top figure shows the solutions for $r > r_c$ where we have two real branches ($r_c \approx 0.93$ for $a = 1/2$). In blue we have the usual branch with free asymptotic behaviour $\epsilon(\theta) \sim r \cosh(\theta)$; and in purple the second branch with asymptotic behaviour as in table 1. In the bottom figure we have the complex solutions below $r_c$, which form a purely imaginary plateau (dashed lines) as $r \to 0$.</figcaption></figure>

allow for a second real branch. The branch with lower effective central charge can be found numerically with the usual iterative method starting from the free solution in the IR. To get very close and eventually pass the critical radius, we need however to use the pseudo-arclength continuation method described in section 3. Close to the turning point the curve is well approximated by a square root, $\tilde{c} \approx \alpha\sqrt{r - r_c}$ for some constant $\alpha$, reminiscent of the situation in $T\bar{T}$ deformed theories. The second branch behaves like $\tilde{c} \sim r^2$ at large radius, as expected.
Note that the free theory limit $a \to 0$ is not smooth, since the $L_1$ norm of the kernel changes discontinuously from 1 for any $a \in (0, 1)$ to 0 for $a = 0, 1$. We find numerically that as the coupling decreases the critical radius grows, see right panel of figure 7, while the second real branch approaches the first one at small scales, see figure 3. Of course, both branches later separate as is required by their corresponding $\tilde{c} \sim 0$ and $\tilde{c} \sim r^2$ asymptotics. The $L$-functions are plotted for different radii in figure 4 (top) and are again consistent with the expected asymptotic behaviour summarized in table 1.

Below $r_c$ solutions to the TBA equation become complex. As mentioned in section 3, to find these solutions it suffices to apply the pseudo-arclength continuation method allowing pseudo-energies to have real and imaginary parts, introducing small perturbations close to the critical radius to switch from the real to the complex branch. Note that once an imaginary part has been acquired by the pseudo-energy the perturbations can be turned off as the corresponding family of solutions is continued further.

As seen in figure 3, the effective central charge reaches in the ultraviolet limit $r \to 0$ a pair of complex-conjugate values that is independent of the coupling $a$. We can verify numerically that the $L$-functions form a plateau at small radii, as shown in figure 4 (bottom), so that the dilogarithm trick outlined in section 2 is applicable. Since the integral of the kernel is $\Omega = 1$ we find two complex solutions,

$$ \epsilon = +\ln(1 - e^{-\epsilon}) \quad \implies \quad \epsilon = \pm i\pi/3. \qquad (4.19) $$

resulting in $\tilde{c}_{\text{shG}}^b(0) \approx 0.5 \pm 0.617i$ which matches the pair of complex conjugate values found numerically, see figure 3. The fact that these values are the same for any coupling can thus be traced back to the $L_1$ norm of the kernel being independent of the coupling $a$.

## 4.3 Fermionic/bosonic elliptic sinh-Gordon

We now turn to the elliptic sinh-Gordon model. As reviewed in section 2 the model has an $S$-matrix with periodicity $S(\theta+T_l) = S(\theta)$, controlled by the modulus $l$. There is an infinite set of resonances for any $l \in (0, 1)$ and no bound states. We can therefore immediately infer that two real branches are allowed in the TBA, and thus expect qualitatively similar behavior to the bosonic sinh-Gordon model, including the existence of a critical radius representing a turning point. However, this needs to be checked numerically since the arguments we have access to at the moment only constrain the asymptotic behavior of the solutions and tells us nothing about their local properties.

As it turns out, the elliptic sinh-Gordon model has two real branches in both its fermionic and bosonic flavors, see figures 5 and 6.<fnref n="13" /> We find that the critical radius where the turning point is found is always greater in the bosonic than in the fermionic case, see figure 7. Moreover, the critical radius decreases as we decrease the modulus $l$ (left panel). In the fermionic case this was expected because in the limit $l \to 0$ we should recover the results for the fermionic sinh-Gordon model, which does not have two branches. In the bosonic case, the critical radius approaches in the $l \to 0$ limit the corresponding value

<footnote n="13">¹³In a previous attempt to solve the TBA equation for this model [59], the authors mention poor convergence of their iterative method at small radii.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Effective central charge for the fermionic elliptic sinh-Gordon model. Plots for various values of $l$ with fixed $a = 1/2$ (top), and various values of $a$ with fixed $l = 1/2$ (bottom). Above a critical radius $r_c$ there are two branches of real solutions, whereas below this radius there is a pair of complex-conjugate solutions with non-vanishing imaginary parts (dashed lines). The effective central charge in the ultraviolet $\tilde{c}(r=0)=0$ can be computed from the plateau solutions as explained in the main text. Numerical solutions were computed using the pseudo-arclength continuation method described in section 3 with $N=200$ and $M=2$.</figcaption></figure>

in the bosonic sinh-Gordon model, cf. rightmost point of the dashed black curve on the right panel.

The dependence of the critical radius on the coupling $a$ is milder than that on the modulus $l$, as can be seen in the right panel of figure 7. We find $r_c$ decreases with $a$, but this is only noticeable when $a \approx 0$. Similarly to the bosonic sinh-Gordon case, the $a \to 0$ limit is not smooth, which again traces back to the fact that the norm of the kernel changes discontinuously at $a=0$.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Effective central charge for the bosonic elliptic sinh-Gordon model. Plots for various values of $l$ with fixed $a = 1/2$ (top), and various values of $a$ with fixed $l = 1/2$ (bottom). Above a critical radius $r_c$ there are two branches of real solutions, whereas below this radius there is a pair of complex-conjugate solutions with non-vanishing imaginary parts (dashed lines). The effective central charge in the ultraviolet $\tilde{c}(r=0) = 0$ can be computed from the plateau solutions as explained in the main text. Numerical solutions were computed using the pseudo-arclength continuation method described in section 3 with $N = 200$ and $M = 2$.</figcaption></figure>

There is a surprising property of the family of real solutions, namely the appearance of an algebraic relation between the kernel $\phi_{a,l}(\theta)$ and the convolution term in the TBA equation. We leave this discussion to appendix B, and now turn instead our attention to the complex solutions.

Evaluating the augmented Jacobian determinant (3.4) along the family of real solutions, we are able to see that the turning point is actually a simple bifurcation point for the complexified theory. Indeed, following the real branch from the standard large $r$ regime
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Critical radius as a function of the modulus $l$ (left, $a = 1/2$) and of the coupling $a$ (right, $l = 1/2$). Results for the standard bosonic, elliptic fermionic and elliptic bosonic sinh-Gordon models are shown in dashed black and solid blue and orange, respectively. On the left panel the colored dashed lines show extrapolations for $l \to 0$: the fermionic case linearly approaches $r_c(l=0) = 0$, whereas the bosonic case has $r_c(l) = r_c(l=0) + \frac{l^2}{4(1-l)^{1/5}}$ with $r_c(l=0) \approx 0.995$, which fits all the datapoints and matches at $l=0$ the critical radius obtained for the bosonic sinh-Gordon model. On the right panel, the critical radii of the bosonic standard and elliptic sinh-Gordon models are very similar for small $a$, displaying more significant differences as $a \approx 1/2$, while staying always above the critical radius of the fermionic elliptic model.</figcaption></figure>

with $\epsilon(\theta) \approx r \cosh(\theta)$ to the non-standard regime with $\tilde{c}(r) \sim r^2$, we have $J = -1$ in the lower branch, i.e. before the critical point, and $J = 1$ afterwards. The precise location of the critical point can thus be pinpointed using binary search on the pseudo-arclength step size $\Delta s$. Below the critical radius $r_c$ we have a pair of complex conjugate solutions, which in this case have $\tilde{c} \to 0$ as $r \to 0$. Interestingly, the UV behavior is only observed very close to $r=0$.

To better understand the $\tilde{c} \to 0$ behaviour, we consider $k$-CDD factors approximating the elliptic sinh-Gordon $S$-matrix, and take the $k \to \infty$ limit. That is, we take the $S$-matrix

$$S_{k\text{CDD}}^{\text{f,b}}(\theta) = \pm \prod_{j=-\lfloor k/2 \rfloor}^{\lfloor k/2 \rfloor} \frac{\sinh(\theta) - \sinh(i\pi a + j T_l)}{\sinh(\theta) + \sinh(i\pi a + j T_l)}, \quad (4.20)$$

and analyze it as before with the help of the TBA. In practice, the differences between the elliptic and the $k$-CDD models above is that the corresponding kernels start to differ for $|\theta| \sim \lfloor k/2 \rfloor T_l$. Since we are truncating all numerical integrations, we need to make sure to include the tails of the $k$-CDD kernel at large $\theta$. Consequently, deviations are only noticeable for small $r$.

Guided by experimental numerical evidence for several $k$'s, let us assume the associated TBA has a turning point and the complex solutions reach a plateau as $r \to 0$. Since the $L_1$ norm of the kernel counts the number of resonances, these constant solutions satisfy equation (2.9) with $\Omega = k$. There is a simple geometric depiction in terms of the $Y$-functions $Y(\theta) = e^{-\epsilon(\theta)}$, for which the equations read

$$Y = (1 \pm Y)^{\pm k}, \quad (4.21)$$

where the upper (lower) signs refer as usual to the fermionic (bosonic) TBA.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Constant solutions to equations (4.21) in the complex plane, where $k$ labels the number of resonances and $Y$ are the constant $Y$-functions $Y = e^{-\epsilon}$ in the ultraviolet. There are $k$ solutions to the fermionic equations (left) and $k+1$ to the bosonic ones (right). The filled circles identify the pair of solutions chosen by the TBA. In the $k \to \infty$ limit, the solutions condense into the circles centered at $Y = \mp 1$, depicted in black. In this limit the TBA gives $Y \to 0$ for both bosonic and fermionic equations, resulting in $\tilde{c} = 0$ in the ultraviolet.</figcaption></figure>

In figure 8 we plot solutions to the above equation for different values of $k$. In general, there are $k$ solutions to the fermionic equation and $k+1$ to the bosonic one.<fnref n="14" /> As $k$ increases, the solutions condense into circles of unit radii centered around $Y = \mp 1$, see figure 8. The TBA, however, picks only two of these $O(k)$ solutions. Which ones? The compelling answer is: the solutions which minimize $|\tilde{c}|$ as computed with Rogers’ dilogarithm (2.10). In terms of the plots in figure 8, the TBA chooses the pair of solutions closer to $Y = 0$. The analysis above is confirmed by our numerical solution to the TBA equations for the first few non-trivial values of $k$, namely $k \le 7$ both in bosonic and fermionic equations. Note that the norm of the kernel is independent from the parameters entering the CDD factors. Moreover, since we can change these parameters continuously and a different solution to (4.21) would involve a discontinuous change, we expect the TBA to choose the solutions closer to $Y = 0$ no matter the resonance locations.<fnref n="15" /> This intuition should be nonetheless confirmed with a detailed numerical analysis.

In figure 9 we show the ultraviolet effective central charge in the complex plane as we increase the number of resonances in the theory. The curve starts at the bosonic 1-CDD model as the bosonic sinh-Gordon model in section 4.2 and goes to zero as the number of resonances goes to infinity as in the elliptic sinh-Gordon models.

<footnote n="14">¹⁴This counting is in agreement with the map between bosonic and fermionic TBA in (4.2).</footnote>

<footnote n="15">¹⁵For instance, for the 2-CDD model it was shown in [28] that the TBA has a turning point for any values of the parameters entering the $S$-matrix, so that the associated complex solutions would agree in the ultraviolet with figure 8.</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 9. Effective central charge $\tilde{c}(r = 0)$ in the complex plane as we increase the number $k$ of resonances in the theory, computed for the bosonic TBA with the constant solutions of figure 8 and Rogers’ dilogarithm formula (2.10). We show only one of the complex conjugate pairs for simplicity. The corresponding curve for fermionic theories is the same provided we map $k \to k + 1$. The first point $k = 1$ in red matches the bosonic sinh-Gordon model value, cf. figure 3. In the limit with infinite resonances $k \to \infty$ the complex effective central charge goes to zero as we find in the elliptic sinh-Gordon models, see figures 5 and 6.</figcaption></figure>

## 5 Bound state explorations: deformed minimal models

In this section we initiate the study of the CDD deformations of models with more than one stable particle. We choose as the undeformed theory the well-known $\Phi_{1,3}$ integrable deformation of the non-unitary minimal models $M_{2,2n+3}$ — the first of which is the Lee-Yang model — and investigate its simplest possible CDD deformation,

$$ \Phi(\theta) = \frac{i \sin u + \sinh \theta}{i \sin u - \sinh \theta}, \qquad (5.1) $$

in the limit $u \to 0^-$. In other words, recalling (4.3) we will effectively focus our attention on the bosonic counterparts of the $\Phi_{1,3}$ integrable deformations. Most curiously, at least to our knowledge these models have not been constructed previously even though they give a standard (single branch) TBA.

These $S$-matrices provide a paradigmatic example of how the iterative method to solve TBA equations can give rise to non-physical instabilities that can be circumvented applying the pseudo-arclength continuation method described in section 3. Furthermore they offer a very different testing ground for the condition on the $L_1$ norm of the kernel (4.1) being the determining factor for the appearance of a critical radius.
Let us now briefly review some properties of the models. The integrable theories resulting from the $\Phi_{1,3}$ deformation have a spectrum of $n$ particles of mass $m_a = \sin(\frac{a\pi}{2n+1})$, $a = 1...n$. The scattering amplitudes defining the interactions between them can be obtained with the usual integrable bootstrap and — contrary to the $S$-matrices in the previous section — have bound-state poles. The scattering matrices read

$$S_{11}(\theta) = \mathrm{th}_{\frac{2}{2n+1}}(\theta), \qquad (5.2)$$

$$S_{ab}(\theta) = \mathrm{th}_{\frac{|a-b|}{2n+1}}(\theta) \mathrm{th}_{\frac{a+b}{2n+1}}(\theta) \prod_{k=1}^{\min(a,b)-1} \left[ \mathrm{th}_{\frac{|a-b|+2k}{2n+1}}(\theta) \right]^2, \qquad (5.3)$$

where the building blocks $\mathrm{th}_x(\theta)$ are simple CDD poles

$$\mathrm{th}_x(\theta) = \frac{\tanh[(\theta + i\pi x)/2]}{\tanh[(\theta - i\pi x)/2]} = \frac{\sinh(\theta) + i \sin(\pi x)}{\sinh(\theta) - i \sin(\pi x)},$$

see [60] and references therein.

Following the logic of the previous section, we will study these fermionic models and their bosonic counterparts defined purely from their $S$-matrix $S_{ab}^{b}(\theta) = -S_{ab}(\theta)$. The TBA equations are given by

$$\epsilon_a(\theta) = m_a r \cosh \theta \mp \sum_{b=1}^{n} \int_{-\infty}^{\infty} \varphi_{ab}(\theta - \theta') \ln[1 \pm e^{-\epsilon_b(\theta')}] \frac{d\theta'}{2\pi} \quad \text{for } a=1, \dots, n, \quad (5.4)$$

and while this is now a coupled system of $n$ equations over the unknown pseudo-energies $\epsilon_a(\theta)$ for $a = 1, \dots, n$, all of the numerical techniques described in section 3 can be straightforwardly generalized to account for this change.

If we solve the TBA equations (5.4) iteratively, we find that for $n > 2$ (or $n > 1$ in the bosonic case) there is an instability at a finite radius $r_* > 0$. For all values of $n$ we notice that the bosonic case displays this instability before the fermionic one. We could suspect that this indicates a true physical critical radius, but switching to the pseudo-arclength continuation method proves this to be false. Indeed, there is nothing strange about these models and we can continue the solutions all the way to $r = 0$, reproducing the well-known effective central charge in the UV for the fermionic models.

Looking closer, it turns out that below $r_*$ in the iterative method even and odd iterations converge to different values.<fnref n="16" /> This is a purely numerical artifact which we can expect to avoid using a variation of the iterative algorithm with weighted-updates, i.e. replacing (3.2) with

$$\epsilon_a^{(k+1)}(\theta) = (1 - \alpha)\epsilon_a^{(k)}(\theta) + \alpha \left( m_a r \cosh \theta \mp \sum_{b=1}^{n} \int_{-\infty}^{\infty} \varphi_{ab}(\theta - \theta') \ln[1 \pm e^{-\epsilon_b^{(k)}(\theta')}] \frac{d\theta'}{2\pi} \right), \qquad (5.5)$$

$$\epsilon_a^{\text{even,odd}}(\theta) = m_a r \cosh \theta \mp \sum_{b=1}^{n} \int_{-\infty}^{\infty} \varphi_{ab}(\theta - \theta') \ln[1 \pm e^{-\epsilon_b^{\text{odd,even}}(\theta')}] \frac{d\theta'}{2\pi}.$$

<footnote n="16">¹⁶Note that none of these values are true solutions to the TBA equations, since we have</footnote>
<figure><img src="image_11.png" /><figcaption>Figure 10. Effective central charge (black) for the $\Phi_{1,3}$ deformed minimal model $M_{2,2n+3}$ with $n=3$, in its fermionic (top) and bosonic (bottom) flavours. The blue curves refer to the single particle contributions $\tilde{c}_a(r)$ as in equation (5.6). The vertical gray line shows the radius $r_*$ where the iterative method becomes unstable. For $r < r_*$ we use the pseudo arclength numerical method. The effective central charges converge in the ultraviolet to 2/3 in the fermionic model and to 0.724253 in the bosonic one.</figcaption></figure>

for $a = 1, \dots, n$, $k = 0, 1, \dots$ and some $\alpha \in (0, 1)$. In practice, $\alpha = 1/2$ works well, and we could skip the pseudo-arclength continuation altogether. Note however the latter is a more robust stabilization method for the TBA equation, as seen in the previous section.

Of course, the fact that there is no critical radius is also expected from the observation that in this theory there are no resonances and only bound states. Indeed, a condition on the norm of the kernel in the spirit of (4.1) would not be satisfied, so that we should not find more than one branch in this theory.
The effective central charge in these models receives contributions from each particle in the following way

$$ \tilde{c}(r) = \sum_{a=1}^{n} \tilde{c}_a(r) = \sum_{a=1}^{n} \pm \frac{6}{\pi} \frac{m_a}{m_1} \int_{-\infty}^{\infty} r \cosh(\theta) \ln[1 \pm e^{-\epsilon_a(\theta)}] \frac{d\theta}{2\pi}. \quad (5.6) $$

In figure 10 we show the fermionic (top) and bosonic (bottom) effective central charges for $n = 3$. In the fermionic TBA, the effective central charge at $r = 0$ converges to the expected values for $M_{2,2n+3}$ minimal models, given by the formula $\tilde{c}^{(n)} = 2n/(2n + 3)$. In the bosonic case the UV effective central charge is given by larger values in the range $[1/2, 1]$, which for $n > 1$ appear to be irrational and for which we have not found an analytic expression.<fnref n="17" /> It would be very interesting to see if there exists a simple realization of these bosonic models.

## 6 Discussion

In this work we have analyzed the Thermodynamic Bethe Ansatz for various integrable models defined through consistent $2 \to 2$ scattering matrices. With the exception of the results discussed in section 5, we concentrated on theories with a single stable particle and no bound states. We studied the sinh-Gordon model, which has one resonance, its elliptic deformation, with infinite resonances, and in-between models with $k$ resonances described by $k$ CDD-zeros. We have taken a bottom-up approach and studied both the fermionic and the lesser studied bosonic realizations of the TBA equations for all of the above models.

Following [28], we have confirmed that the determining factor for the TBA to have a turning point at a critical length scale $r_c$ is the $L_1$ norm of the kernel, which measures the difference between the number of bound states and resonances in the theory. Namely, all the models with a single stable particle and any number of resonances — with the exception of the single resonance fermionic model, i.e. sinh-Gordon — have a bifurcation point. The exact length scale at which the bifurcation happens depends on the details of the model.

An important result in the present work is that we were able to continue the TBA solutions past the turning point and follow the pair of complex conjugate solutions all the way to the UV regime. This required an efficient implementation of the pseudo-arclength continuation method to solve the (complexified) TBA equations. We found that the complex solutions to the TBA equations are such that the UV effective central charge is minimized. Moreover, we showed that the latter decreases with the number of resonances and finally goes to zero as we consider an infinite number of them, as in the elliptic sinh-Gordon models.

The existence of a turning point in the TBA — interpreted as the presence of a finite Hagedorn temperature — and of a complex effective central charge at short distances implies a non-Wilsonian UV behaviour for these theories, just as it is the case in the $T\bar{T}$-deformed theories and in their generalizations [28]. Our results shed light on the possible nature of the theories with infinite resonances ubiquitous in the $S$-matrix bootstrap, and

<footnote n="17">¹⁷For the curious reader we report $\tilde{c}(r = 0) \approx 0.5, 0.641304, 0.724253, 0.778979, 0.817083$ for $n = 1, 2, 3, 4, 5$.</footnote>
add to the body of experimental evidence substantiating the expectation that the majority of self-consistent $S$-matrices are not derivable from local QFT. Indeed, there are further hints that some of the consistent $S$-matrices found saturating bounds cannot describe conventional UV-complete theories. In particular, for the $O(N)$ periodic Yang-Baxter (describing a global symmetry generalization of our infinite resonance model) there are two concrete pieces of evidence. First, in [61] it was argued with an $S$-matrix/Form Factor bootstrap involving c-minimization that this theory seems incompatible with a standard UV fixed point. Secondly, the authors of [62] observed walking behaviour [63, 64] in this model for $N > 2$ and proposed a relation to complex CFTs.

With the above hints and the fact that we are obtaining complex effective central charges, it is very tempting to think of the underlying UV theories as complex CFTs. This is an intriguing idea that warrants further investigation. A first step would be to extract the scaling dimensions of primary operators by considering the excited state TBA equations obtained by analytically continuing the groundstate ones [65]. Regardless of the details of the underlying UV theories, we believe that analyzing theories outside the paradigm of a standard UV fixed point is a worthy endeavour. Moreover, while our work focused on integrable theories we expect their qualitative properties to be of a universal nature, applying equally well to non-integrable and non-Wilsonian QFTs, just as it is the case for the $\bar{T}\bar{T}$ deformation.

Finally, the results discussed in section 5 constitute a first incursion into the realm of CDD deformations of theories with an arbitrary number of stable particles. Specifically we studied the bosonic counterparts of the $S$-matrices describing the $\Phi_{1,3}$ deformed non-unitary minimal models $M_{2,2n+3}$. We found that these theories possess a standard Wilsonian UV limit, corresponding to CFTs whose central charges appear to be irrational numbers between 1/2 and 1. We find it remarkable that these models have never, to the extent of our knowledge, appeared before in the literature. We believe it would be interesting to look at the bosonic counterparts of other well-known theories — e.g. the unitary minimal models — as well as to search for their QFT construction. The next natural step in this direction is to analyze a general family of theories having $n$ bound states and $N$ resonances. According to the argument we have already put forward before, we expect these to have a standard Wilsonian UV behavior for $n \ge N$, and conversely to develop a Hagedorn temperature for $n < N$. We hope to return to this point in a future publication.

## Acknowledgments

We would like to thank Yifei He, Alexandre Homrich, Slava Rychkov, Roberto Tateo, Alexander Zamolodchikov and Bernardo Zan for useful and insightful discussions. Work of SN is supported by NSF grant PHY-1915219.

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (QUASIFT grant agreement 677368).
## A Some properties of elliptic functions

Elliptic functions are meromorphic functions in the complex plane with two periods $ω₁, ω₂$ satisfying Im($ω₁/ω₂$) > 0. In this appendix we review some of their nice properties.<fnref n="18" />

**Poles and zeros.** Any non-constant elliptic function must have at least one pole in a period parallelogram. The sum of all residues at the poles inside this parallelogram is zero. Therefore, any non-constant elliptic function has at least order two: inside the period parallelogram we have either two simple poles with opposite residues or a double pole with zero residue. The former case gives rise to Jacobi elliptic functions while the latter corresponds to Weierstrass elliptic functions. Any elliptic function can be expressed in terms of Jacobi or Weierstrass elliptic functions.

**Algebraic relations.** Two elliptic functions $X(θ)$, $Y(θ)$ with the same periods $ω₁, ω₂$ satisfy an algebraic relation of the form

$$F(X, Y) = \sum_{\substack{i,j=0 \\ i+j \le n}}^{n} c_{ij} X^i Y^j = 0. \qquad (A.1)$$

In the following we describe how these polynomial relations arise and in particular how the total degree of the polynomial is related to the order of the poles in $X(θ)$, $Y(θ)$.

Let $a₁, a₂, ..., a_m$ be the set of distinct poles that belong to at least one of the functions. The functions might have a pole at the same position with different orders,

$$X(\theta) \sim \frac{r_k}{(\theta - a_k)^{m_k}}, \quad Y(\theta) \sim \frac{r'_k}{(\theta - a_k)^{n_k}}. \qquad (A.2)$$

We define the total pole order $K$ as the sum of the largest order for each pole,

$$K = \sum_{k=1}^{m} \max(m_k, n_k). \qquad (A.3)$$

Now let us consider an elliptic function $Φ(θ)$ defined by a polynomial with no constant term

$$\Phi(\theta) = \sum_{\substack{i,j=0 \\ 0<i+j \le n}}^{n} c_{ij} X^i Y^j. \qquad (A.4)$$

The idea is to fix the coefficients $c_{ij}$ such that $Φ(θ)$ is an entire function and therefore equals a constant by Liouville's theorem.<fnref n="19" /> The function $Φ(θ)$ can have poles only at $a₁, a₂, ..., a_m$, so for it to reduce to a constant the principal parts near these points must all vanish. For example, close to $a_k$ the function $Φ(θ)$ has a Laurent expansion

$$\Phi(\theta) = \sum_{l=0}^{\infty} A_l(\theta - a_k)^l + \sum_{l=1}^{n \times \max(m_k, n_k)} B_l(\theta - a_k)^{-l}. \qquad (A.5)$$

<footnote n="18">¹⁸There are many textbooks on the theory of elliptic functions, here we follow [66].</footnote>

<footnote n="19">¹⁹This theorem holds because the function would also be bounded in the period parallelogram with vertices {0, $ω₁$, $ω₁ + ω₂$, $ω₂$}.</footnote>
We require that the $B_l$ coefficients in the second sum vanish for all $l$. In this way we get $n \times \max(m_k, n_k)$ equations $B_l[c_{ij}, X(a_k), Y(a_k)] = 0$. Repeating this procedure for all poles gives at most $n \sum_k \max(m_k, n_k) = nK$ linear homogeneous equations with variables $c_{ij}$. For a non-trivial solution to exist we need to have more unknowns than equations, i.e. $n(n+3)/2 > nK$. This leads to a bound on the total degree of the polynomial in terms of the total pole order, $n > 2K - 3$.

Having determined the coefficients $c_{ij}$ we can now write $\Phi(\theta) = -c_{00}$, from which we recover a non-trivial algebraic equation as (A.1).

## B Approximate algebraic relation between kernel and convolution

A curious fact stems from the otherwise trivial observation that in the elliptic sinh-Gordon models the convolution term in the TBA equation (2.5) has the same real period as the kernel $\varphi_{a,l}(\theta)$. If these functions were to share their imaginary period as well it would imply an algebraic relation between them, as explained in appendix A. Given the analytic structure of the kernel, we would expect this relation to be of at least second order, i.e. to be of the form

$$c_{00} + c_{10} C(\theta) + c_{01} \varphi_{a,l}(\theta) + c_{11} C(\theta) \varphi_{a,l}(\theta) + c_{20} C(\theta)^2 + c_{02} \varphi_{a,l}(\theta)^2 = 0, \quad (B.1)$$

where $C(\theta) = \epsilon(\theta) - r \cosh(\theta)$. Note that if such a relation holds, it would not only require some physical interpretation but also be of practical importance, since it would permit a much more efficient analysis and solution of the TBA equation.

Having at our disposal numerical solutions for the convolution, determining the coefficients $c_{ij}$ reduces to performing a linear regression on the $2NM$ datapoints (3.3) to minimize the mean squared deviation to (B.1). Restricting ourselves to the real branches, we find in practice that $c_{11} = c_{02} = 0$, so that fixing the normalization through $c_{20} = 1$ we are left with only three undetermined coefficients, $c_{00}$, $c_{10}$ and $c_{01}$. In our experiments with $N = 200$ and $M = 2$, fitting these three parameters produces residuals of relative magnitude $O(10^{-5})$ for all real solutions. This is very close to the accuracy of the solutions themselves, therefore we can assert that the algebraic relation (B.1) is indeed satisfied to the numerical precision we are working with. For illustration purposes, in figure 11 we display plots of the non-zero coefficients as a function of $r$, for the fermionic (left panels) and bosonic (right panels) models at $a = 1/2$ and $l = 1/2$.

While we can verify the algebraic relation (B.1) holds for the real branch of solutions to the TBA equations, we were unable to do the same for the solutions belonging to the complex branch. In fact, in order to prove this relation, it would be necessary to establish the analytic structure of the convolution in the complex plane, as well as the non-trivial periodicity in the imaginary direction. This involves a non-trivial analytic continuation of the TBA solutions from the real rapidity line into the complex plane, which we are not able to do reliably all the way to $\theta \in \mathbb{R} + 2\pi i$ as needed to confirm the validity the algebraic relation. We leave a deeper exploration of this interesting curiosity to future work.
<figure><img src="image_12.png" /><figcaption>Figure 11. Non-vanishing coefficients in the algebraic relation (B.1) between the kernel $\varphi_{a,l}(\theta)$ and the convolution $C(\theta) = \epsilon(\theta) - r \cosh \theta$ of the elliptic sinh-Gordon models with $a = 1/2$ and $l = 1/2$. Panels on the l.h.s. correspond to the fermionic model, whereas panels on the r.h.s. correspond to bosonic model.</figcaption></figure>

## C Pseudo-arclength continuation of the TBA equation

In this appendix we explain in more detail the application of the pseudo-arclength continuation method to the TBA equation (2.5). After introducing the discretization (3.3), we have to solve

$$ \epsilon_i = r \cosh \theta_i \mp \frac{\Delta\theta}{2\pi} \sum_{j=1}^{2NM} (\varphi_{a,l})_{ij} \ln (1 \pm e^{-\epsilon_j}) \quad \text{for} \quad i = 1, 2, \dots, 2NM, \quad (C.1) $$

where $\Delta\theta = T_l/N$, $(\varphi_{a,l})_{ij} = \varphi_{a,l}(\theta_i - \theta_j)$ and $\epsilon_i = \epsilon(\theta_i)$ are $2NM$ unknowns. We can accomplish this task using standard Newton-like iterative methods broadly applicable to
the solution of nonlinear equations, as long as we can provide a good-enough initial starting point for the discretized pseudo-energy $\epsilon_i$.

We consider the situation where the parameters $a$ and $l$ in the kernel above are fixed beforehand, and we want to solve this problem repeatedly for various values of $r$. In this setup, small changes in $r$ and the $\epsilon_i$ are related by

$$ \delta\epsilon_i = \delta r \cosh \theta_i + \frac{\Delta\theta}{2\pi} \sum_{j=1}^{2NM} (\varphi_{a,l})_{ij} \frac{\delta\epsilon_j}{e^{\epsilon_j} \pm 1} \quad \text{for} \quad i=1,2,\dots,2NM. \quad (C.2) $$

If we now take $(\delta\epsilon_i, \delta r)$ to be unknown, we need to introduce an additional equation through the normalization condition

$$ \frac{1}{2NM} \sum_{i=1}^{2NM} \delta\epsilon_i^2 + \delta r^2 = 1, \quad (C.3) $$

which leaves us with $2NM + 1$ equations as required. Once again, providing a starting point for $(\delta\epsilon_i, \delta r)$ that is close enough to the actual solution is all that is necessary to solve these equations applying standard Netwon-like methods.

Assuming we already have a sequence of $n$ solutions to (C.1) we wish to extend,<fnref n="20" /> say $\{\epsilon_i^{(1)}, \epsilon_i^{(2)}, \dots, \epsilon_i^{(n)}\}$ corresponding to $\{r^{(1)}, r^{(2)}, \dots, r^{(n)}\}$, we can use its last two elements to provide an initial guess for $(\delta\epsilon_i, \delta r)$,

$$ \delta\epsilon_i \approx \mathcal{N}(\epsilon_i^{(n)} - \epsilon_i^{(n-1)}) \quad \text{and} \quad \delta r \approx \mathcal{N}(r^{(n)} - r^{(n-1)}), \quad (C.4) $$

where the normalization constant $\mathcal{N}$ is given by

$$ \mathcal{N} = \left( \frac{1}{2NM} \sum_{i=1}^{2NM} (\epsilon_i^{(n)} - \epsilon_i^{(n-1)})^2 + (r^{(n)} - r^{(n-1)})^2 \right)^{-1/2} . \quad (C.5) $$

Solving the derivative equation (C.2) and normalization condition (C.3) for $(\delta\epsilon_i, \delta r)$ then allows us to construct an approximate solution to the original TBA equation. Indeed, introducing a pseudo-arclength parameter $\Delta s$ controlling how far we want to move away from the last available solution $\epsilon_i^{(n)}$, we can use $\epsilon_i \approx \epsilon_i^{(n)} + \delta\epsilon_i \Delta s$ as an initial guess for the solution of (C.1). This corresponds to $r \approx r^{(n)} + \delta r \Delta s$, therefore we need to promote $r$ to become an $(2NM + 1)$-th unknown, while introducing an additional constraint

$$ \frac{1}{2NM} \sum_{i=1}^{2NM} \delta\epsilon_i(\epsilon_i - \epsilon_i^{(n)}) + \delta r(r - r^{(n)}) = \Delta s. \quad (C.6) $$

Thus, the solution of (C.1) and (C.6) provides a new term $\epsilon_i^{(n+1)}$ in our sequence, and the process can be repeated ad arbitrium.

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

<footnote n="20">²⁰These may come, for example, from the application of the iterative method described in section 3.</footnote>
References

[1] Y. Nakayama, Scale invariance vs conformal invariance, Phys. Rept. 569 (2015) 1
[arXiv:1302.0884] [INSPIRE].

[2] F.A. Smirnov and A.B. Zamolodchikov, On space of integrable quantum field theories, Nucl. Phys. B 915 (2017) 363 [arXiv:1608.05499] [INSPIRE].

[3] A. Cavaglià, S. Negro, I.M. Szécsényi and R. Tateo, $T\bar{T}$-deformed 2D quantum field theories, JHEP 10 (2016) 112 [arXiv:1608.05534] [INSPIRE].

[4] A.B. Zamolodchikov, Expectation value of composite field $T\bar{T}$ in two-dimensional quantum field theory, hep-th/0401146 [INSPIRE].

[5] S. Dubovsky, V. Gorbenko and M. Mirbabayi, Asymptotic fragility, near $AdS_2$ holography and $T\bar{T}$, JHEP 09 (2017) 136 [arXiv:1706.06604] [INSPIRE].

[6] S. Dubovsky, V. Gorbenko and G. Hernández-Chifflet, $T\bar{T}$ partition function from topological gravity, JHEP 09 (2018) 158 [arXiv:1805.07386] [INSPIRE].

[7] R. Conti, S. Negro and R. Tateo, The $T\bar{T}$ perturbation and its geometric interpretation, JHEP 02 (2019) 085 [arXiv:1809.09593] [INSPIRE].

[8] T. Ishii, S. Okumura, J.-I. Sakamoto and K. Yoshida, Gravitational perturbations as $T\bar{T}$-deformations in 2D dilaton gravity systems, Nucl. Phys. B 951 (2020) 114901 [arXiv:1906.03865] [INSPIRE].

[9] J. Cardy, The $T\bar{T}$ deformation of quantum field theory as random geometry, JHEP 10 (2018) 186 [arXiv:1801.06895] [INSPIRE].

[10] M. Baggio and A. Sfondrini, Strings on NS-Ns backgrounds as integrable deformations, Phys. Rev. D 98 (2018) 021902 [arXiv:1804.01998] [INSPIRE].

[11] A. Hashimoto and D. Kutasov, $T\bar{T}$, $J\bar{T}$, $T\bar{J}$ partition sums from string theory, JHEP 02 (2020) 080 [arXiv:1907.07221] [INSPIRE].

[12] A. Sfondrini and S.J. van Tongeren, $T\bar{T}$ deformations as $T\bar{s}T$ transformations, Phys. Rev. D 101 (2020) 066022 [arXiv:1908.09299] [INSPIRE].

[13] N. Callebaut, J. Kruthoff and H. Verlinde, $T\bar{T}$ deformed CFT as a non-critical string, JHEP 04 (2020) 084 [arXiv:1910.13578] [INSPIRE].

[14] A.J. Tolley, $T\bar{T}$ deformations, massive gravity and non-critical strings, JHEP 06 (2020) 050 [arXiv:1911.06142] [INSPIRE].

[15] L. McGough, M. Mezei and H. Verlinde, Moving the CFT into the bulk with $T\bar{T}$, JHEP 04 (2018) 010 [arXiv:1611.03470] [INSPIRE].

[16] A. Giveon, N. Itzhaki and D. Kutasov, $T\bar{T}$ and LST, JHEP 07 (2017) 122 [arXiv:1701.05576] [INSPIRE].

[17] S. Chakraborty, A. Giveon and D. Kutasov, $T\bar{T}$, $J\bar{T}$, $T\bar{J}$ and String Theory, J. Phys. A 52 (2019) 384003 [arXiv:1905.00051] [INSPIRE].

[18] R. Hagedorn, Statistical thermodynamics of strong interactions at high-energies, Nuovo Cim. Suppl. 3 (1965) 147 [INSPIRE].

[19] S. Dubovsky, R. Flauger and V. Gorbenko, Effective string theory revisited, JHEP 09 (2012) 044 [arXiv:1203.1054] [INSPIRE].
[20] S. Dubovsky, R. Flauger and V. Gorbenko, *Solving the simplest theory of quantum gravity*, JHEP 09 (2012) 133 [arXiv:1205.6805] [INSPIRE].

[21] M. Caselle, D. Fioravanti, F. Gliozzi and R. Tateo, *Quantisation of the effective string with TBA*, JHEP 07 (2013) 071 [arXiv:1305.1278] [INSPIRE].

[22] K.G. Wilson and J.B. Kogut, *The renormalization group and the $\epsilon$-expansion*, Phys. Rept. 12 (1974) 75 [INSPIRE].

[23] M. Guica, *An integrable Lorentz-breaking deformation of two-dimensional CFTs*, SciPost Phys. 5 (2018) 048 [arXiv:1710.08415] [INSPIRE].

[24] R. Conti, L. Iannella, S. Negro and R. Tateo, *Generalised Born-Infeld models, Lax operators and the $T\bar{T}$ perturbation*, JHEP 11 (2018) 007 [arXiv:1806.11515] [INSPIRE].

[25] B. Le Floch and M. Mezei, *$KdV$ charges in $T\bar{T}$ theories and new models with super-Hagedorn behavior*, SciPost Phys. 7 (2019) 043 [arXiv:1907.02516] [INSPIRE].

[26] R. Conti, S. Negro and R. Tateo, *Conserved currents and $T\bar{T}_s$ irrelevant deformations of 2D integrable field theories*, JHEP 11 (2019) 120 [arXiv:1904.09141] [INSPIRE].

[27] G. Hernández-Chifflet, S. Negro and A. Sfondrini, *Flow equations for generalized $T\bar{T}$ deformations*, Phys. Rev. Lett. 124 (2020) 200601 [arXiv:1911.12233] [INSPIRE].

[28] G. Camilo, T. Fleury, M. Lencsés, S. Negro and A. Zamolodchikov, *On factorizable S-matrices, generalized TTbar, and the Hagedorn transition*, JHEP 10 (2021) 062 [arXiv:2106.11999] [INSPIRE].

[29] A.B. Zamolodchikov, *Integrable field theory from conformal field theory*, Adv. Stud. Pure Math. 19 (1989) 641.

[30] L. Castillejo, R.H. Dalitz and F.J. Dyson, *Low's scattering equation for the charged and neutral scalar theories*, Phys. Rev. 101 (1956) 453 [INSPIRE].

[31] R.J. Eden, P.V. Landshoff, D.I. Olive and J.C. Polkinghorne, *The analytic S-matrix*, Cambridge University Press, Cambridge U.K. (1966).

[32] D. Iagolnitzer, *Macrocausality, physical region analyticity and independence property in S matrix theory*, Lect. Notes Math. 449 (1975) 102.

[33] D. Iagolnitzer, *Scattering in quantum field theories: the axiomatic and constructive approaches*, Princeton University Press, Princeton U.S.A. (1994).

[34] A.M. Perelomov and Y.B. Zeldovich, *Quantum mechanics*, Selected Topics, World Scientific, Singapore (1998).

[35] A.B. Zamolodchikov, *Two point correlation function in scaling Lee-Yang model*, Nucl. Phys. B 348 (1991) 619 [INSPIRE].

[36] G. Mussardo and S. Penati, *A quantum field theory with infinite resonance states*, Nucl. Phys. B 567 (2000) 454 [hep-th/9907039] [INSPIRE].

[37] M.F. Paulos, J. Penedones, J. Toledo, B.C. van Rees and P. Vieira, *The S-matrix bootstrap. Part I: QFT in AdS*, JHEP 11 (2017) 133 [arXiv:1607.06109] [INSPIRE].

[38] M.F. Paulos, J. Penedones, J. Toledo, B.C. van Rees and P. Vieira, *The S-matrix bootstrap II: two dimensional amplitudes*, JHEP 11 (2017) 143 [arXiv:1607.06110] [INSPIRE].

[39] Y. He, A. Irrgang and M. Kruczenski, *A note on the S-matrix bootstrap for the 2d O(N) bosonic model*, JHEP 11 (2018) 093 [arXiv:1805.02812] [INSPIRE].
[40] L. Córdova and P. Vieira, *Adding flavour to the S-matrix bootstrap*, JHEP 12 (2018) 063
[arXiv:1805.11143] [INSPIRE].

[41] M.F. Paulos and Z. Zheng, *Bounding scattering of charged particles in 1+1 dimensions*,
JHEP 05 (2020) 145 [arXiv:1805.11429] [INSPIRE].

[42] J. Elias Miró, A.L. Guerrieri, A. Hebbar, J. Penedones and P. Vieira, *Flux tube S-matrix
bootstrap*, Phys. Rev. Lett. 123 (2019) 221602 [arXiv:1906.08098] [INSPIRE].

[43] A. Homrich, J. Penedones, J. Toledo, B.C. van Rees and P. Vieira, *The S-matrix Bootstrap.
Part IV. Multiple amplitudes*, JHEP 11 (2019) 076 [arXiv:1905.06905] [INSPIRE].

[44] L. Córdova, Y. He, M. Kruczenski and P. Vieira, *The O(N) S-matrix Monolith*, JHEP 04
(2020) 142 [arXiv:1909.06495] [INSPIRE].

[45] C. Bercini, M. Fabri, A. Homrich and P. Vieira, *S-matrix bootstrap: supersymmetry, $Z_2$, and
$Z_4$ symmetry*, Phys. Rev. D 101 (2020) 045022 [arXiv:1909.06453] [INSPIRE].

[46] M. Kruczenski and H. Murali, *The R-matrix bootstrap for the 2d O(N) bosonic model with a
boundary*, JHEP 04 (2021) 097 [arXiv:2012.15576].

[47] M. Hortacsu, B. Schroer and H.J. Thun, *A two-dimensional $\sigma$ model with particle
production*, Nucl. Phys. B 154 (1979) 120 [INSPIRE].

[48] A.B. Zamolodchikov, *Thermodynamic Bethe ansatz in relativistic models. Scaling three state
Potts and Lee-Yang models*, Nucl. Phys. B 342 (1990) 695 [INSPIRE].

[49] S.J. van Tongeren, *Introduction to the thermodynamic Bethe ansatz*, J. Phys. A 49 (2016)
323005 [arXiv:1606.02951] [INSPIRE].

[50] A. Tsvelick and P. Wiegmann, *Exact results in the theory of magnetic alloys*, Adv. Phys. 32
(1983) 453.

[51] T.R. Klassen and E. Melzer, *Purely elastic scattering theories and their ultraviolet limits*,
Nucl. Phys. B 338 (1990) 485 [INSPIRE].

[52] E.L. Allgower and K. Georg, *Numerical continuation methods: an introduction*, Springer,
Germany (2012).

[53] E.L. Allgower and K. Georg, *Numerical path following, Handbook of numerical analysis 5*
(1997) 3.

[54] O. Babelon and L. Bonora, *Conformal affine sl(2) Toda field theory*, Phys. Lett. B 244
(1990) 220 [INSPIRE].

[55] A.B. Zamolodchikov, *Mass scale in the sine-Gordon model and its reductions*, Int. J. Mod.
Phys. A 10 (1995) 1125 [INSPIRE].

[56] S.N. Vergeles and V.M. Gryanik, *Two-dimensional quantum field theories having exact
solutions*, Sov. J. Nucl. Phys. 23 (1976) 704 [Yad. Fiz. 23 (1976) 1324] [INSPIRE].

[57] A.B. Zamolodchikov, *Resonance factorized scattering and roaming trajectories*, J. Phys. A 39
(2006) 12847 [INSPIRE].

[58] A.B. Zamolodchikov, *On the thermodynamic Bethe ansatz equation in sinh-Gordon model*, J.
Phys. A 39 (2006) 12863 [hep-th/0005181] [INSPIRE].

[59] O.A. Castro-Alvaredo and A. Fring, *Constructing infinite particle spectra*, Phys. Rev. D 64
(2001) 085005 [hep-th/0103252] [INSPIRE].
[60] G. Mussardo, Statistical field theory, Oxford Graduate Texts, Oxford University Press, Oxford U.K. (2020).

[61] D. Karateev, S. Kuhn and J. Penedones, Bootstrapping massive quantum field theories, JHEP 07 (2020) 035 [arXiv:1912.08940] [INSPIRE].

[62] V. Gorbenko and B. Zan, Two-dimensional $O(n)$ models and logarithmic CFTs, JHEP 10 (2020) 099 [arXiv:2005.07708] [INSPIRE].

[63] V. Gorbenko, S. Rychkov and B. Zan, Walking, weak first-order transitions, and complex CFTs, JHEP 10 (2018) 108 [arXiv:1807.11512] [INSPIRE].

[64] V. Gorbenko, S. Rychkov and B. Zan, Walking, weak first-order transitions, and complex CFTs II. Two-dimensional Potts model at $Q > 4$, SciPost Phys. 5 (2018) 050 [arXiv:1808.04380] [INSPIRE].

[65] P. Dorey and R. Tateo, Excited states by analytic continuation of TBA equations, Nucl. Phys. B 482 (1996) 639 [hep-th/9607167] [INSPIRE].

[66] K. Chandrasekharan, Elliptic functions, Springer, Germany (2012).