Approaching the self-dual point of the sinh-Gordon model

One of the most striking but mysterious properties of the sinh-Gordon model (ShG) is the b → 1/b self-duality of its S-matrix, of which there is no trace in its Lagrangian formulation. Here b is the coupling appearing in the model’s eponymous hyperbolic cosine present in its Lagrangian, cosh(bϕ). In this paper we develop truncated spectrum methods (TSMs) for studying the sinh-Gordon model at a finite volume as we vary the coupling constant. We obtain the expected results for b ≪ 1 and intermediate values of b, but as the self-dual point b = 1 is approached, the basic application of the TSM to the ShG breaks down. We find that the TSM gives results with a strong cutoff Ec dependence, which disappears according only to a very slow power law in Ec. Standard renormalization group strategies — whether they be numerical or analytic — also fail to improve upon matters here. We thus explore three strategies to address the basic limitations of the TSM in the vicinity of b = 1. In the first, we focus on the small-volume spectrum. We attempt to understand how much of the physics of the ShG is encoded in the zero mode part of its Hamiltonian, in essence how ‘quantum mechanical’ vs ‘quantum field theoretic’ the problem is. In the second, we identify the divergencies present in perturbation theory and perform their resummation using a supra-Borel approximate. In the third approach, we use the exact form factors of the model to treat the ShG at one value of b as a perturbation of a ShG at a different coupling. In the light of this work, we argue that the strong coupling phase b > 1 of the Lagrangian formulation of model may be different from what is naïvely inferred from its S-matrix. In particular, we present an argument that the theory is massless for b > 1.


Introduction
The sinh-Gordon model (ShG) is a canonical quantum integrable field theory. It has a number of different descriptions, but in this work we are going to take as a starting point the Lagrangian formulation of the model given by (∂ µ φ) 2 − 2µ cosh(bφ) .

JHEP01(2021)014
Here φ(x, t) is a real non-compact scalar field, µ is some dimensionful mass scale and b is a dimensionless coupling constant. Upon quantization, µ is replaced by a renormalized coupling constant depending on the chosen quantization scheme. The spectrum of the model is exceedingly simple, consisting of a single massive particle of mass M ShG . This makes the ShG the simplest of interacting integrable field theories. Much is known about its properties. Its elastic S-matrix was found in [1]. The form factors of local operators were obtained in [2,3], while the vacuum expectation values of exponential operators were found in [4]. The exact relationship between the physical mass M ShG and the renormalized coupling in the perturbed gaussian CFT scheme, here denoted later as µ ShG , as a function of the coupling b, was derived in [5]. Its thermodynamic Bethe ansatz for the ground state and for the excited states was studied in [6,7], while the thermal correlation functions of the model were discussed in [8][9][10]. A suggestive connection between the ShG model and roaming renormalization group trajectories among the minimal models of CFT was studied in [11], while a direct mapping between the ShG and the Ising model was established in [12]. Furthermore, beyond simply being a model that is amenable to analytic manipulation, the ShG finds applications in a wide range of areas of physics running from toy models of quantum gravity [13], to cold atomic gases [14,15], studies of thermalization in classical field theories [16,17], and lattice models with non-compact quantum group symmetries [18]. It is also worth stressing that the ShG model is the simplest example of Toda field theories, a large class of models with exponential interactions based on root systems of Lie algebras, see for instance [19] and references therein. The main difference between the ShG model and the rest of the Toda field theories is that the ShG does not have bound states.
One of the most striking but mysterious aspects of the ShG model 1 is its apparent weak-strong duality: In the presence of such a symmetry, the self-dual point b = 1 clearly emerges as a special value of the ShG model, for it divides the weak-coupling regime, b < 1, from the strong coupling regime, b > 1. It is important to underline that this duality is not at all manifest in the Lagrangian of the theory but is apparent, as discussed later, in its S-matrix formulation. It is the primary aim of this paper to develop truncated spectrum methods (TSMs) in order to study the model at finite volume by varying its coupling constant in the vicinity of this self-dual point. For reasons which will become clear later, the obvious regime in which these methods can be implemented is the weak-coupling regime b < 1, but we can extend them to approach the b = 1 self-dual point. As we shall see, at b = 1 the physical mass vanishes and the theory appears to be (at least naively) critical. Ultimately we aim to explore the ShG at b = 1 and understand what theory is described by the Lagrangian given in eq. (1.1). The theory's duality, as expressed in eq. (1.2), is built on results established at b < 1 which are then subsequently analytically continued to regimes beyond their nominal validity. It is then an important question to understand whether the Lagrangian corresponding to these analytic continuations is the same as given in eq. (1.1) with b > 1. 1 The Toda field theories have a similar duality.

JHEP01(2021)014
As we explore in this paper, the development of TSM for the ShG nearby the self-dual point b = 1 proves to be surprisingly challenging. Truncated spectrum methods treat a model by firstly defining it on a finite volume (typically an infinitely long cylinder of width R) and then, secondly, introducing a hard UV cutoff, E c , in the number of energy levels which are included in the computation [20][21][22][23][24][25][26][27][28][29][30]. Under these two conditions, numerics can be performed (either exact diagonalization or Lánczos based approaches) and the low lying energy spectrum, together with vacuum expectation values and matrix elements of several operators, can be computed. Of course, in this treatment it is crucial to understand the effect of E c on the computed results. For certain models, even small values of E c lead to results that are, in effect, independent of the cutoff (i.e. the c = 1/2 Ising model perturbed by the presence of a magnetic field being a classic example [21]). In other cases, as for instance those analysed in refs. [22,[31][32][33][34][35][36], the results are instead noticeably affected by E c and, to ameliorate the effects of the introduction of E c , various renormalization group (RG) strategies have been employed, both analytic and numerical [23,24,[27][28][29][30][37][38][39].
The premise of all these RG strategies is that cutoff dependent effects are in some sense small. However, in the case of the ShG model, we shall see that such cutoff effects near b = 1 can be on the contrary extremely large and therefore the traditional RG strategies do not work. In order to deal with this new situation, we propose herein three different approaches to tackle the problem: 1. In the first, we explore more carefully the small-volume regime and its 'quantum mechanical nature'. In particular, one may expect that the UV behaviour of the spectrum is dominated by the quantum mechanics of the zero mode of the field. By using either this quantum mechanical picture or the TBA equations combined with the mass-coupling relation, it is possible to derive a systematic expansion for certain energy levels (more precisely, their corresponding scaling functions) in terms of (ln(µ ShG R)) −1 . The two expansions are however different in subleading orders. As the energies contain an additional R −1 factor relative to the scaling function, the difference between TBA and zero mode energy levels eventually diverge for R → 0 for all b > 0. We derive an effective potential, partially taking into account the effect of oscillators. At the one hand, we analytically reproduce the exact expansion up to O(b 12 ), confirming that the oscillators are able to explain the differences in the logexpansion. On the other hand, we show that TSM numerics significantly outperforms even the numerical solution of the complete effective potential. We then use this fact to provide a more precise measurement of the IR parameters from TSM, combining UV numerics with the small-volume expansion of TBA.
2. In the second strategy, we recast the analytic RG strategy used to remove the effect of the cutoff. Typically this RG strategy is pursued by initially performing low order perturbation theory in the conformal coupling, here µ ShG . However this fails for the ShG near b = 1 as the perturbation theory of this model is divergent term by term. These divergences, we show, actually appear for any value of b, although for small values of b their appearance is delayed until higher orders. Facing this, we argue that the diverging perturbative series can be resummed. However, this is not -3 -

JHEP01(2021)014
a Borel resummation per se, as the series is diverging more rapidly than n!, but it does nonetheless admit a supra-Borel resummation.
3. In the third and last strategy, we abandon the use of the non-compact boson Hilbert space as a computational basis. One way to understand the difficulties in using TSM's about b = 1 is to think of them as arising due to a poor choice of computational basis.
We have already said that the theory becomes critical at b = 1 (i.e. the mass scale M ShG vanishes at fixed µ ShG ). Hence, using a non-interacting field to describe the vicinity of what it could be a non-trivial conformal field theory (presumably strongly interacting) may then simply be inappropriate. Thus we explore the possibility of using an interacting basis of states as a computational basis. The natural choice here is to use, as a computational basis, the basis of exact eigenstates of the ShG at one value of b to study the theory at a different (relatively close) value of b.
The paper is organized as follows. In section 2, we review basic information on the ShG model, pointing out its origin and possible limitations. Although this section reviews previously known results, it is crucial for understanding the rest of the paper. In section 3 we discuss truncated spectrum methods, in particular the key role played by the choice of computational basis. In section 4 we then present our particular choice of computational basis. In section 5, we discuss our numerical results for various quantities including the finite volume spectrum, the S-matrix, and the vacuum expectation values of various exponential operators of the model. In section 5 we further demonstrate how the standard renormalization group techniques used to improve TSM results fail to do so for the ShG model close to the self-dual point, thus setting up the rationale for the next three sections. In section 6 we explore in detail the information carried by the zero modes of the theory and the 'quantum-mechanical' nature of the ShG model in certain regimes of the coupling and volume. In section 7 we analyse the nature of perturbation theory which defines the ShG model as a massive deformation of a Gaussian theory and we argue that the perturbative series is badly behaved and is non-Borel resummable. This leads us to consider a supra-Borel resummation in order to give meaning to these divergent sums. In section 8 we come back to the issue of a proper choice of the basis for the TSM and we explore the possibility to study the ShG model at a given coupling b in terms of states and matrix elements of a ShG model defined at a different value of b. As we will see, this approach admits of a series of sanity checks. In section 9 we finally discuss our conclusions and future directions.

Basic features of the sinh-Gordon model
In this section we briefly review the basic properties of the ShG necessary to understand the TSM results and their interpretation presented in the main body of the paper. The scale dimension of the renormalized counterpart of the coupling µ appearing in the ShG Lagrangian depends on the quantization scheme of the model. Hereafter we are going to discuss three such schemes: (i) a perturbative scheme based on Feynman diagrams; (ii) -4 -JHEP01(2021)014 T < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > δ m 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > + < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > = 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > Figure 1. Diagrams for the first order tadpole diagram T and the corresponding relative mass counterterm δm 2 = −T (m 0 ). treating the theory on the same grounds as its analytically continued cousin, the sine-Gordon model, namely as a perturbed Gaussian CFT; and (iii) as a perturbation of a Liouville quantum field theory. The model's self-duality is often encoded in the parameter which we record here for the reader to emphasize its importance.

Feynman diagrammatic analysis
In the first scheme, the ShG model is considered by employing perturbation theory in the coupling constant b and evaluating all quantities in terms of Feynman diagrams. This can be done by introducing a momentum cutoff Λ and expanding the potential of the theory in terms of b: Here µ(μ, Λ, b) is a bare parameter of dimension mass squared. Using the cutoff Λ, we have introduced a renormalized couplingμ, which we aim to keep fixed as we tune µ such that the physical quantities are finite: Above the unique ground state of the theory, there is a massive excitation, whose mass at the lowest order in b is m 0 given by Of course the actual mass of the particle will get corrections by all the higher order interactions. However, the perturbative series contains divergences. Fortunately, in 1+1d theories with local interactions all divergences come from the tadpole diagrams. These divergences can be cured by introducing a mass counter-term δm 2 and imposing that, order by order, δm 2 cancels the infinities coming from the tadpole diagrams. At the lowest order in b 2 , for instance, we have the condition expressed in figure 1, where the tadpole is regularized in terms of the momentum cutoff Λ as The counterterm δm 2 = −T (m IR ), involving in general an arbitrary mass scale m IR , is absorbed by the bare parameter µ such that This prescription is equivalent to defining a normal ordering for the Lagrangian (eq. (1.1)). The quantization scheme is fixed by the choice of m IR . In particular, setting m IR = m 0 leads to the usual scheme of a perturbed massive boson, where normal ordering is with respect to the free mass m 0 , eliminating altogether the tadpole diagrams at each order. The exact relation between µ andμ in the normal ordering scheme m IR is easily obtained by means of the Baker-Campbell-Hausdorff formula. It reads (2.7) In this way, all n-point correlation functions of the theory are finite to all orders in perturbation theory. In particular, one can compute the physical mass M ShG of the theory, as a function of m 0 and b 2 , by looking at the pole of the 2-point correlation function. In the scheme m IR = m 0 ≡ m, we obtain where g = 1 8π and each term comes from the Feynman diagrams of figure 2. We will point out later that this perturbative analysis is consistent with an exact formula for the mass that we present in section 2.2.4.
We note that at O(b 4 ), the ShG coincides with a φ 4 Landau-Ginzburg model. Given the repulsive nature of this latter theory, the ShG is expected to have no bound states, its spectrum consisting of multi-particle states of the same particle. As we will see shortly, this conclusion is in agreement with the exact S-matrix of the model.

Relation with the sine-Gordon model
In the second approach, properties of the ShG model are extracted from a closely related model, the sine-Gordon (SG) model. The SG model has a Lagrangian given by This can be obtained from the ShG Lagrangian (1.1) by making the substitutions It is important to stress that, although the two theories are related by this simple transformation, their underlying nature is rather different and there are indeed a series of hidden subtleties behind the innocent looking analytic continuation (2.10), some of which are discussed below.

SG and ShG models as deformations of a Gaussian theory
Both theories may be regarded as deformations of the Gaussian fixed point action given by the kinetic term in eq. (1.1) With respect to this CFT of central charge c = 1, the chiral conformal dimension of a vertex operator, V (a) = e iaφ , is ∆(a) = a 2 . The sine-Gordon model involves the vertex operators V (±b) = e ±i bφ which are compact and bounded, while the sinh-Gordon model employs the vertex operators V (∓ib) = e ±bφ which are instead non-compact and unbounded. Moreover, while in the sine-Gordon model the conformal dimensions of the vertex operators are positive and given by ∆ ±b = b 2 , (2.12) in the sinh-Gordon model they are instead negative and given by How the sinh-Gordon model turns out to be a unitary quantum field theory, despite the negative conformal dimension of its basic vertex operators, is one of the remarkable aspects of this model. The way the theory restores its unitarity is through the existence of non-zero vacuum expectation values (VEV), whose exact values are provided in eq. (2.33) below. With x 12 = x 1 − x 2 , consider for instance the operator product expansion (OPE) with respect to the Gaussian fixed point: (2.14) and taking the vacuum expectation value of both terms of this equation, we have cosh(bφ(x 1 )) cosh(bφ(x 2 )) cosh(2bφ(0)) Hence, if cosh(2bφ(0)) = 0, we see that the two-point function of the vertex operators e ±bφ has effectively the same leading short-distance singularity as it would have in the case of a positive ∆ = b 2 conformal dimension.

Coleman bound in SG and its formal absence in ShG
From a renormalization group point of view, the vertex operators which give rise to the sine-Gordon model are relevant operators for b 2 ≤ 1, where the upper value b 2 = 1 is known in the literature as Coleman's bound [40]. The values 0 ≤ b 2 ≤ 1 are those for which the SG is ultraviolet stable (i.e. we do not need extra non-trivial counter-terms in its Lagrangian to cure its ultraviolet divergencies). As we already know, the only divergences come from the tadpoles, which can be absorbed by a normal ordering prescription under which the vertex operators get renormalized multiplicatively. Defining m IR as the mass scale by which normal ordering is defined and using a −1 as the UV cutoff, the multiplicative renormalization appears as When b 2 > 1, the vertex operators are irrelevant: hence the SG model becomes essentially a massless theory [41]. In the ShG model, the renormalization of the operators (or the coupling) occurs with the inverse factor of the SG model (2.17) Typically in this multiplicative renormalization the power of a that arises is absorbed into the bare coupling µ so defining a renormalized dimensionful parameter µ SG/ShG : where + is for the SG theory and − is for the ShG model. µ SG/ShG then has engineering dimension in the two theories of 2 ∓ 2b 2 . The scale m IR that appears in this multiplicative renormalization is then typically absorbed into the definition of the normal ordered vertex operator so that the OPE has the conventions expressed in eq. (2.14). Henceforth it is understood as part of the definition of µ SG/ShG that m IR is chosen this way. The relation between the free mass m appearing in eq. (2.8) and the coupling µ ShG is [42] µ ShG = m 2+2b 2 2 4+2b 2 πb 2 e 2b 2 γ E , (2.19) as derived in appendix A. Because of the negative conformal dimension of its vertex operators (which makes them relevant operators), at least formally the ShG model does not have a Coleman bound. However, according to the argument given above, the singularity structure of the OPE for the ShG interaction (eq. (2.15)) is the same as for the SG model. One thus may suspect that there is in fact a Coleman bound for the ShG model, namely that the theory is properly defined only for b 2 < 1, has a singularity at b 2 = 1 and a massless phase for b 2 > 1. This is the scenario we will actually present later in the paper.

The spectrum of SG model
Let us now turn our attention to the spectrum of the SG model. This quantity is key as the spectrum and S-matrix of the SG model will be connected to that of the ShG model by analytic continuation. Reproducing this spectrum will be one of the major targets of our TSM studies.
We note that this analytic continuation is subtle. While the sinh-Gordon model has only one vacuum state, the sine-Gordon model has instead an infinite number of vacuum states, |n , which are associated to the minima of the potential, φ n = 2πn/b. These multiple vacua give rise to solitons and anti-solitons, excitation which interpolate between two neighboring vacua, |n and |n ± 1 . For the integrability of the theory, scattering among solitons and anti-solitons is elastic and the relative amplitudes can be computed exactly [43]. Here it is sufficient to remind the reader of the main results of this analysis. It is convenient to define 20) as this parameter controls the spectrum of the SG theory. The number of neutral solitonanti-soliton bound states (breathers) is given by where [x] denotes the integer part of x. Denoting by M s the mass of the soliton, the breather masses are given by Hence, the first breather exists provided ξ ≤ 1, namely only in the range b 2 ≤ 1/2. Ignoring this restriction, we plot m 1 (b) vs b for the entire interval (0, 1) (see figure 3.b). Notice that even though the breather does not exist for b 2 > 1/2, its mass remains positive until b 2 = 2/3. After this, its value turns negative and begins to rapidly oscillate, reflecting its possession of an essential singularity at b 2 = 1.
The mass scale M s can be related to the renormalized coupling of the theory µ SG (as defined in eq. (2.18)). In the SG model the ground state energy in finite volume and in the presence of an external field coupled to the topological charge of the model can be computed in two different ways: using the thermodynamic Bethe ansatz (TBA) and using conformal perturbation theory. The former approach employs the physical mass M s while the latter, the renormalized mass scale µ SG . Comparing the results coming from the two different approaches, Al. Zamolodchikov [5] was able to obtain an exact formula encoding the µ SG − M s relation: (2.23) We see that this formula is consistent with µ SG in the SG model having engineering and scaling dimension 2 − 2b 2 . It is also important to stress that this formula assumes the vertex operators are normalized with the convention of eq. (2.14).  This formula is physical in the interval 0 ≤ b 2 ≤ 1. It has an essential singularity when Its behaviour in the interval 0 ≤ b 2 ≤ 1 is shown on figure 3a.

ShG model as analytic continuation of SG model
As we have stated, the ShG model can be thought of as the analytic continuation of the SG. How then to connect the rich spectrum of SG containing topological excitations and their bound states to the much simpler spectrum of ShG consisting of a single parity odd excitation? The choice typically made is to identify the first breather of the SG model with the massive excitation of ShG.
To obtain the mass, M ShG , of the fundamental excitation in ShG, we then take M ShG (b) = m 1 (ib). Using eqs. (2.23) and (2.22) we arrive at [4] where we have replaced µ SG with −µ ShG -necessary as only µ ShG has the correct engineering dimension. The plot of this quantity can be found on figure 4. A few remarks are in order: 1. Keeping µ ShG fixed, the mass formula (eq. (2.26)) is not invariant under the weakstrong duality b → 1/b of the ShG model. Moreover, its analytic continuation for b 2 > 1 gives generally complex values for M ShG .  2. For any finite value of µ ShG , the mass M ShG vanishes both at b 2 = 0 and b 2 = 1. The nature of these zeros is however very different. Indeed, the zero at b 2 = 0 disappears if we rescale µ ShG → µ ShG /(8πb 2 ) adopting the more conventional definition of the coupling constant of the model used in the Feynman diagram expansion. However on the approach to b 2 = 1, we instead find a singular point: 3. In order to compare with the perturbative expansion (eq. (2.8)), it is necessary to connect the renormalization scheme used to define µ ShG with that used to define m (eq. (1.1)). Substituting eq. (2.19) into (eq. (2.26)) and expanding in b 2 , we have agreeing with the series expansion of the square root of expression eq. (2.8).

S-matrix of the ShG model and its duality
Having argued the fundamental excitation of the ShG model is to be identified with the analytically continued breather of SG, we are now in a position to derive the S-matrix of the ShG's excitation. This S-matrix is nothing but the analytically continued breatherbreather S-matrix, S B 1 B 1 (θ), of SG which is given by where θ = θ 1 − θ 2 and θ i (i = 1, 2) is the rapidity of each of the breathers involved in the scattering, with energy and momentum given by E i = m 1 cosh θ i and p i = m 1 sinh θ i . This amplitude has a pole at θ = iπξ, and, for ξ < 1, its residue is positive, i.e. this pole signals a further bound state, while for ξ > 1 the residue changes sign, which can be interpreted as another signal of the absence in the spectrum of the first breather for ξ > 1.
If we now continue this expression analytically by substituting −b 2 for b 2 in eq. (2.29), we obtain for the exact 2-body S-matrix of the ShG model (2.32) This expression coincides with the S-matrix of the ShG model proposed in [1]. Although this argument is amazingly simple, the final result is nonetheless surprising because a duality has appeared. The S-matrix now is invariant under the weak/strong duality b ↔ 1/b or B ↔ 1 − B. However this duality is nowhere apparent in the Lagrangian (1.1) of the model.

Vacuum expectation values in the SG and ShG model
Similar to the mass and S-matrix, the vacuum expectation values (VEVs) of the vertex operators of the ShG model can be obtained as analytic continuations from the corresponding expressions for the SG model, the eponymous FLZZ formula [4,44]: (2.33) with G(α) = G(−α). Notice that the integral above converges for a bound conceived for physical operators by N. Seiberg in his study of the allied Liouville problem [45]. For values of α beyond the Seiberg bound, one can exploit an analytic continuation of G(α). Obtaining this continuation is facilitated by the expression [46]: From it one can see that the VEV, as a function of α, does not have poles but only zeros.
Besides the zero at α = Q/2, there is an infinite set of generically simple zeros located at: with m ≥ n ≥ 2 n ∈ 2 Z or n > m ≥ 2 m ∈ 2 Z.

JHEP01(2021)014
Formula (2.35) is positive for − Q 2 < α < Q 2 , but it changes sign at its zeroes. The vertex operators, being the exponentials of Hermitian operators, are positive (semi-)definite. This means that, outside the above domain, the analytic continuation cannot directly correspond to the expectation value. We will therefore consider these values "unphysical". The function G(α) itself is self-dual, i.e. invariant under b → 1/b -for the proof one may benefit from the identity but the VEV is not itself self-dual because of the presence of M −2α 2 ShG . From the dependence on α of this term, we can infer that the scaling dimension of the vertex operator V (α) is ∆(α) = −α 2 , a value which coincides with its conformal dimension with respect to the Gaussian fixed point. Using the VEV (2.33) we can compute the expectation value of the trace of the stress-energy tensor, an operator that, on general terms, is defined as where β(µ ShG ) is the β-function of the coupling µ ShG which perturbs a critical point and O its conjugate field. For the case at hand, we have Using the mass formula (eq. (2.26)) and the simplified expression of the VEV at α = b

Questions arising from the analytic continuation b ↔ ib
In this section we have presented a number of results for the ShG model (its spectrum, its S-matrix, and the VEVs of its exponential operators) that are arrived at by analytically continuing results from SG. The question of the validity of these analytic continuations has to be raised. While the S-matrix of the ShG model is physically sensible for all b, the expressions for the mass and VEVs are not. Given that the fundamental excitation of the ShG is identified with the breather of SG and the SG breather ceases to exist for b < 1/ √ 2, what exactly can be said for b > 1/ √ 2 is not entirely clear. And certainly the mass formula for M ShG for b > 1 breaks down entirely giving complex-valued results.
That we are able to match perturbative computations of the mass formula with the exact expression is thus important. This gives us some confidence that the results for M ShG are valid for b < 1. This confidence will be increased in the following section where we discuss the ShG model as a perturbation of a Liouville theory. However the validity and interpretation of formulae at b > 1 including the S-matrix arising from the analytic continuation remains, in our opinion, an open question.

ShG and Liouville models
We have now considered the ShG model from the perspective of perturbation theory and an analytically continued SG model. We now present a third way to look at the ShG model: as a deformation of a Liouville field theory [4,47,48]. This third way will be essential for us in what follows and results presented here will be used in our discussion of quantum mechanical reductions of the ShG model in section 6.
The Liouville conformal field theory is defined by the action [4,47] Here the operator e bφ has conformal dimension 1 and µ L is dimensionless. This vertex operator has this dimension because we have coupled the field to a background charge of In general, in the presence of the charge at infinity, the conformal dimension of the vertex operator e αφ becomes In order to ensure the theory is IR finite, the theory can be placed on a Riemann sphere (of area A) with a metric The ShG model then can be obtained as a perturbation of the Liouville action: Note that we have made a scale d that comes from normal ordering the vertex operator explicit instead of absorbing it into µ L . We use d here to distinguish this scale from ones (i.e. a) previously introduced in normal ordering vertex operators in different (Gaussian) schemes. µ L here is the same dimensionless constant that appears in the original Liouville action, eq. (2.41).
A key property of the Liouville field theory is that the exponential operators are pairwise identified as [47] e αφ( where R(α) is related to the Liouville reflection amplitude S L (P ) as This identification of the operators implies that their VEV must satisfy the reflection relations

JHEP01(2021)014
Here the L-subscripts indicate that we are taking these VEVs as defined in the perturbed Liouville formulation of ShG and are not assuming these relations are the same as in eq. (2.35). A solution of these equations can be obtained by an infinite iteration With the further assumption of minimality, we can find a result equivalent to combining (2.26) and (2.33). This provides below an alternate way of understanding these formulae without resorting to analytic continuation from results derived for SG. In order to present this argument, we need to trace carefully the dimensions of the quantities involved. This is an issue because in the Liouville approach the exponential operator, e aφ , has dimension (2.42) whereas the dimension of this operator in the perturbed Gaussian formulation of the ShG model is instead −α 2 . To understand this dimensional transmutation, we follow along with ref. [4] and interpret properly the results.
To begin, we use the action (eq. (2.44) to write the following perturbative expansion of the VEV of e aφ : Since in Liouville field theory the coupling µ L can be absorbed into the field φ via a redefinition of this field by an additive shift, it is easy to obtain the explicit µ L dependence of all correlators appearing in eq. (2.49). With the IR regulator in place, this expression can be written as the following series: whereG Ln (α) is independent of µ L . The prefactor µ −α/b L also ensures the satisfiability of the reflection relations (2.47), given the dependence on µ L of the refection amplitude R(α), see eq. (2.46). We see that our IR regulator appears in a dimensionless combination with the scale d in this perturbative expansion. If we assume that this expression has a sensible large A limit, the above series must behave asymptotically as .
(2.51) Thus in the large A limit, we obtain LG L (α). (2.52) We thus see the VEV has dimension −2α(α−Q) as set by the UV cutoff d -the dimension set by the Liouville CFT. At the same time its dependence upon µ L is exactly what would be expected from thinking of the ShG model as a perturbation of a free non-compact boson.
If we identify the couplings in the two formulations via we can identify the expressions for the VEVs in the two formulations via The appearance of the factor d in this expression is a reflection of the different normal ordering schemes in the Gaussian vs Liouville pictures.

Generalized TBA equations for ground and excited state energies
One of the tools that we will use extensively in characterizing our TSM data is the thermodynamic Bethe ansatz (TBA). The exact ground state energy, E 0 (R), at a finite volume R from the TBA by using the 'finite volume-finite temperature' equivalence of the partition function Z: is the 'vacuum' energy density equal to the VEV of eq. (2.33): The pseudo-energy (θ) in eq. (2.56) is defined as the solution of the nonlinear integral equation where the kernel is related to the S-matrix as For excited states, the exact finite volume energies can be obtained either by careful analytic continuations of the ground state TBA (following [49]), or by examining the continuum -16 -JHEP01(2021)014 limit of an integrable lattice regularization [7]. A finite volume n-particle state can thus be described by the multiparticle pseudo energy (θ|{ϑ j } n j=1 ) ≡ (θ|{ϑ}) and a set of quantization numbers {I j } n j=1 satisfying the non-linear integral equation together with the additional quantization conditions Given quantization numbers I j ∈ Z, the rapidities {ϑ} and the pseudo energy (θ|{ϑ}) can be determined self-consistently by efficient numerical methods. Note that the above phase convention is 'bosonic' in the sense that it is permitted to have I j = I k , j = k (see also [50]). Nevertheless, the resulting rapidities are always different, reflecting the fermionic nature of the particles. The above ingredients provide the finite volume energy of the multiparticle state as Notice that, neglecting all terms exponentially suppressed in M ShG R in the finite volume expression of the energies, these are nothing else but the well-known Bethe ansatz equations of the ShG model. Moreover, regarding the mass M ShG just as a parameter of these equations (i.e. as a quantity independent on b), these equations are invariant under the duality b ↔ 1/b present in the S-matrix.

Summary
Let us summarise the main points of this section: • As a Lagrangian theory, the ShG model has three equivalent descriptions: (i) one arising from perturbation theory in b; (ii) one as a deformation of a Gaussian c = 1 CFT; and (iii) finally one as a deformation of a Liouville theory.
• In the Gaussian picture, the ShG model can be thought of as an analytic continuation of the SG model. In this way, results can be derived for the mass spectrum, the Smatrix, and the VEVs of the vertex operators.
• The S-matrix so obtained suggests the theory has a weak-strong duality: b ↔ 1/b. However the mass formulas (and so the VEVs) are not invariant under this duality.
• The validity of this analytic continuation for values of b > 1 is unclear and there are doubts about its validity even for b > 1/ √ 2.
• There is however reason to believe the expressions for the mass and the VEV up to b = 1 because of the availability of an alternate derivation of the VEVs using the Liouville picture as well as consistency with perturbation theory.

JHEP01(2021)014
• For the ShG model, we have a formalism that allows us to compute in a numerically exact fashion the excited state energies at any volume R. These equations are invariant under duality if we consider M ShG to be merely a parameter In the next two sections we will present the results coming from the truncated spectrum method and testing the various expressions for the mass M ShG , the S-matrix, S(θ), and the VEVs e αφ presented here. This analysis will give us some insight, even if not definitive, on the nature of the theory for b > 1.

Truncated spectrum methods
The purpose of this section is to give the reader an overview of truncated spectrum methods (TSMs) and their application to the sinh-Gordon model. TSMs were introduced by Yurov and Zamolodchikov [20] to study the low-energy spectrum of 2D perturbed conformal field theories. However the method is able to study the spectrum and matrix elements of any theory whose Hamiltonian can be conveniently written as a sum of two terms where H 0 is a base theory of which we assume to have a complete control of its energy eigenvalues and eigenstates |E n 0 . In particular, we assume that we are able to write down the matrix elements of the second term in the full Hamiltonian, V pert (x), in the basis {|E n 0 } of eigenvectors of H 0 . From a computational point of view, the actual implementation of the method requires both a denumerable set of energy states and finite-dimensional subspaces of the Hilbert space of the model: the former condition is typically achieved by putting the model onto a cylinder of finite circumference R in the spatial direction; the latter condition is satisfied by restricting the set of eigenstates of H 0 to those whose energies fall below a cutoff E c . Once a finite basis is obtained in this way, the truncated Hamiltonian is constructed. This operator possesses the same matrix elements as the original Hamiltonian in the truncated subspace, but acts trivially in the orthogonal subspace. Having this in hand, one then solves, numerically, the eigenproblem of the truncated Hamiltonian.
Assuming for the moment that the dependence of the data on the cutoff E c is smooth, once the truncated Hamiltonian is diagonalized for different volumes, the infinite volume quantities can be obtained by extrapolation via Lüscher's principles [31,[51][52][53]. TSMs were first applied to the scaling Lee-Yang model [20] and the Ising model [21]. In both cases the numerical results reported therein were strongly convergent in E c . In the study of the perturbed tri-critical Ising theory, though, it was argued that the convergence in E c of the TSM results depends on the scaling dimension of the perturbing operator [22,31]. Various renormalization group approaches have been advocated to treat cases where convergence in E c is suboptimal. These strategies are both numerical [23] and analytical [37][38][39]54] (For a comprehensive review of such strategies see [24].) We will demonstrate later that these strategies require modification (at the very least) for the case of the ShG.

JHEP01(2021)014
The performance of TSMs is dependent on the choice of the computational basis used to perform the calculations (or, in other words, how we split the Hamiltonian H into H 0 and dxV pert ). As with any variational method, we want to use a computational basis that captures at the start features of the physics of the model at hand. In this paper we study the ShG model by means of two different choices of H 0 or computational bases: 1. In the first, discussed in section 4, we consider the ShG model as a deformation of a Gaussian CFT and the corresponding non compact bosonic field expanded in terms of an infinite number of oscillators and a single zero mode. When H 0 is a compact CFT on a cylinder and V is a relevant operator, one typically has control on the magnitude of the interaction between different energy scales in the theory. The ShG model is, however, different, and the low and high energy scales in the problem become strongly coupled on the approach towards b = 1.
2. This leads us to our second basis choice, outlined in section 8, where we use the basis of the ShG model itself as the computational basis. In particular, we use the basis of the ShG model at one value of b to compute the properties of the model at a different value of b. In this scheme, V pert is the difference between two hyperbolic cosines. This approach does immediately raise questions of circularity. We are, after all, using conjectured information about the model at a point b 0 as input, to obtain results at point b 1 = b 0 . We will address this question in section 8, and attempt to ameliorate this concern.

TSM for the ShG using a non-compact massless bosonic basis
In our first attempt to study the ShG model using TSMs, we employ a computational basis based on a non-compact massless basis. In this section we review the details surrounding this choice of basis.

Non-compact massless boson
In describing this basis the starting point is the mode expansion of the massless noncompact bosonic field on an infinite cylinder of radius R: where the oscillators are subject to the usual Fock commutator relations, while the zero mode defines an effective 1D quantum mechanical system with the canonical commutator [ϕ 0 , Π 0 ] = i.
The computational basis of states follows from the specification of H 0 in eq. (3.1). Here we will divide H 0 into a zero mode H ZM and non-zero mode H N ZM part: : e bϕ 0 : + : e −bϕ 0 : ; where the Virasoro generators L 0 andL 0 appearing in H N ZM are related to the Fock mode operators as Notice that H ZM , unlike H N ZM , is an interacting Hamiltonian. With this writing of H 0 , our computational eigenbasis has a tensor product structure composed of a zero mode and an oscillator sector: Here H osc is decomposed into a chiral subspace, H R , is spanned by right-moving particles and an anti-chiral subspace, H L , is spanned by left-moving particles

Zero modes
Unlike the non-interacting H N ZM , we have chosen a form for H ZM that is non-trivial. We do so following [26,28,42] so that H ZM consists of a countable (i.e. discrete) basis of states. We will denote this basis of states as follows Unlike the states in H N ZM , in our implementation of the TSM the eigenstates |m will be found numerically. To do so, we need to choose a computational basis to represent H ZM in eq. (4.3) and, for this aim, we choose the position basis in the zero mode coordinate where 2L is the length of the truncated zero mode space (rather than having a noncompact zero mode, we assume it lies between −L and L) and a is our spatial discretization parameter. In performing our computations here, we have always taken L large enough and a small enough so that the eigenvalues and eigenstates (or at least their matrix elements) of H ZM have converged completely.

Truncated Hilbert spaces
Having determined H ZM , we are now in a position to define the truncated basis, H T , for the problem as a whole. We truncate each part of the Hilbert space separately. In particular we write The cutoff is then implemented in terms of two separate parameters, N c , the level of which we cut off the chiral oscillator mode part of the Hilbert space, and N ZM , the number of zero mode eigenstates of smallest energy (w.r.t. to H ZM ) that we keep. We typically work not in this full space, but its zero-momentum counterpart composed of tensored states from H R,T and H L,T that satisfy The last step of forming the Hamiltonian matrix consists of specifying the interaction part of the Hamiltonian and its corresponding matrix elements. In the zero-momentum subspace, P = 0, we can write where δ P,0 reflects the projection onto the zero momentum subspace and we have separated out the zero mode from the field: In the above, normal ordering is defined as where The matrix elements of the chiral parts of : e bφ(x) : admit a closed analytic expression where the chiral state vector |n 1 , . . . , n k is a normalized state having n q right-moving particles with momentum 2πq R for each q ∈ {1, . . . , k}. Using this expression together with the knowledge of the (numerical) zero mode matrix elements m|e bϕ 0 |n , we can construct the full matrix elements of H int .

Methods of diagonalization
Once we have collected all the matrix elements of both H 0 and H int , for any truncated space we have a finite dimensional matrix to diagonalise. To find its eigenstates and eigenvalues we can proceed in two ways: • We can use exact diagonalization perhaps augmented with a numerical renormalization group. This latterprocedure will be discussed further in section 5.2.
• We can also use iterative methods that, thanks to their reverse communication protocols, do not require us to store the full Hamiltonian in memory. The only cost that we need to pay is that we are restricted here to computing the low-lying eigenvalues. However for our purposes here this is not a limitation. Using the Jacobi-Davidson method and exploiting the tensor product structure (i.e. eq. (4.4)) of the Hilbert space, we can treat matrices arising from truncation parameters of up to (N c = 20, N ZM = 24), corresponding to a truncated Hilbert space of size approximately 2.4 × 10 7 . We elaborate on the usage of the tensor product structure in appendix E. We specifically use the JDQMR_ETOL algorithm 2 provided in the package PRIMME [55,56].

TSM results for ShG model
In this section we present our TSM results based on the non-compact bosonic computational basis for various quantities in the ShG model. In the first part, we show how the TSM results for the ShG are robust at small b (b 1), but begin to have strong cutoff effects in N c (and, to a lesser extent, N zm ), deviating noticeably from exact results predicted by the thermodynamic Bethe ansatz (TBA), as b exceeds 1/2. We then discuss our first strategy in dealing with these cutoffs: a power law extrapolation in N c (and N zm ). We show that this procedure in fact produces robust results -albeit still imperfect as b → 1 is approached. We then turn to other quantities that we are able to measure using TSM methods, such as the VEVs of the exponential operators and the S-matrix.
In the second part of this section we consider standard renormalization group strategies for alleviating the effects of the cutoff. We show that while these strategies work at small values of b, they lead to sub-optimal (or even unphysical) results at larger values of b which JHEP01(2021)014 are closer to the self-dual point. However this failure provides the motivating drive to consider other strategies for treating the sensitivity of TSM results to the cutoff that form the next three sections that follow this one. It will also set the scene for understanding why the power law extrapolation used in the first part of this section is robust.

Results
We present our numerical results with the aim to answer the following specific questions: 1. What is the performance of the TSM applied to the ShG model? What level of precision can be achieved below the self-dual point (compared directly to finite volume theoretical quantities), and how does it depend on the coupling b and the volume R?
2. Is there a simple extrapolation that robustly improves the accuracy of the numerics? How much does it improve?
3. Not assuming any special properties of the ShG model (in other words, relying only on standard TSM analysis), to what extent can the conjectured infinite volume parameters (mass, vacuum energy density, S-matrix) be reproduced?

How effectively does TSM reproduce the one-point functions of vertex operators?
In particular, what happens when we probe them outside the region of validity of the FLZZ formula, eq. (2.33)?
The following results are organized according to the four points listed above.

Finite volume spectrum
Let us begin with the finite volume spectrum. In figure 5 we present results for the ground state energy and the first excited state energy for different values of b at fixed volume M ShG R = 1, where M ShG is the physical mass. The computations were done using different chiral cut-offs N c , at a fixed zero mode cutoff N ZM . On the other hand, we have also calculated the corresponding quantities by numerically integrating the (excited state) TBA equations eq. (2.62). We consider the difference between TSM and TBA (taking into account the vacuum energy density (2.57)) to be the error of the former. The energy levels are normalized with respect to the free mass m defined in eq. (2.19) and we plot the differences between the TBA computations and the TSM data on a log scale. The largest cutoff, N c = 12, N ZM = 24, at which data are presented has been obtained using a truncated Hilbert space of size 3 × 10 5 . It is apparent that the errors are slightly different for the ground state and the excited state, but the overall pattern is very similar: for small b, even a raw cutoff can produce precise results, and a reasonable increase in the cutoff N c actually has a strong positive effect on the precision. On the other hand, the error increases exponentially in increasing the coupling constant b and, at the same time, the precision becomes less sensitive to the cutoff. In the immediate vicinity of the self-dual point, the error essentially becomes O(1), indicating that the naive TSM is limited to a region below the self-dual point. As a first JHEP01(2021)014   Here we present the two extrapolation schemes for the cutoffs N c and N ZM . We perform the extrapolation for data of the ground state energy at b = 0.95, M ShG R = 1). In a) we first extrapolate in N ZM , plotting the results at each N c , then performing a further extrapolation in N c . The results of this extrapolation is shown in the legend. In b) we reverse the order of extrapolations. Note that the data has a much stronger dependence on N c than on N ZM . step to improve the results, we propose an (at this point 'empirical') extrapolation scheme, which involves fitting the numerical results with power laws in N c and N ZM .
In figure 6, we present two implementations of this power law fitting at b = 0.95, M ShG R = 1. In the first, presented in figure 6(a), we first extrapolate in zero mode number N ZM and then perform a further extrapolation in N c . Figure 6(a) then shows the result of this first extrapolation in N ZM at each N c . We see that the results of this extrapolation differ little from the raw data determined at N ZM = 42. Having done this, we then perform a separate extrapolation with respect to the chiral cutoff N c . The result, reported in the legend of figure 6(a) is about a 3% error in units of the free mass m. We  also perform the same extrapolation in N c for the raw data obtained at N ZM = 42. The result is essentially identical. In figure 6(b), we consider the second implementation. This is obtained by performing the two extrapolations in the opposite order but for the same set of input data. Shown in figure 6(b) is the result of the first extrapolation in N c at fixed N ZM . The second extrapolation in N ZM lead to results essentially the same as the results reported with the first scheme. We now consider TSM data over a range of values of b. Having seen that the data is essentially converged at N ZM sufficiently large, we work at fixed N ZM and only consider extrapolations in N c . In particular we use the fitting function:

JHEP01(2021)014
In figure 7 we report our results for the ground state and first excited state energies at M ShG R = 1. We show the raw TSM data at different N c together with the extrapolated values (red dots). In most cases, the extrapolation improves the numerical data by at least an order of magnitude. We note that the cusp-like feature appearing in the extrapolated data is due to a sign change of the extrapolated error, of which we take the absolute value to produce the log-scale plot. The precise position of the cusp is also volume-dependent.
In figure 8 we present the fitting exponent ν(b) (see eq. (5.1)) coming from these extrap-   Here we present TSM data (dots) for the first 8 energy levels after N c extrapolation as functions of the dimensionless volume M ShG R in the zero momentum sector. The vacuum energy E 0 is subtracted and energies are normalized with respect to the mass M ShG . The color coding of the TBA (solid) curves is as follows: the ground state is depicted in blue, the one-particle state in green, two-particle states in red, three-particle states in orange, four-particle states in brown, five-particle states in teal, and six-particle states in pink.
olations. We see that at small b the exponent is large indicating that the data is rapidly converging in N c while at values of b approaching the self-dual point, the exponent becomes much smaller. We will provide a partial explanation for the behavior of the power law in section 7.
Having presented results as a function of b, we now consider the spectrum as a function of R. In figure 9 we present data for the low-lying finite volume spectrum (after extrapolation in N c ) after subtraction of the exact vacuum energy density (2.57) for two different couplings, b = 0.4 and b = 0.8. The TSM data is plotted against the numerical solution of the exact TBA equations (shown in the plots with continuous curves). It is apparent that TSM follows very closely the theoretical excited-TBA data for b = 0.4, while small discrepancies become visible at b = 0.8, especially at larger volumes. Contrary to the previous plots, here we opted for normalizing the energies with respect to the physical mass M ShG , owing to the emphasis of finite volume corrections presented in these plots.
We close the first part of this subsection with a contour plot which shows the order of magnitude of errors as a function of both b and the dimensionless volume M ShG R at the same time. The error can be smaller than 10 −8 in the small-R region of the perturbative sector, and remains below 10 −4 over a wide range of couplings and volumes. On the other hand, the error increases exponentially as either the volume or the coupling is increased. We note that the apparent 'islands' on the top and the 'valley' around M ShG R = 9 are due to the same sign-changing phenomenon that causes the cusps in figure 7.

Determination of mass, bulk energy and S-matrix
In this subsection, we present and discuss the TSM numerical results for the particle mass, M ShG , the bulk energy density, and the S-matrix of the ShG model. Physical mass. We have measured the physical mass M ShG of the ShG model through taking the difference of the two lowest energy levels. Ideally, this difference converges to the physical mass in the R → ∞ limit. In practice, for large volumes, truncation effects produce an overestimate for the mass. On the other hand, small-volume effects also produce an overestimate. As a consequence, we determine the mass as the minimum of the volume-dependent energy difference. The results are shown as the function of b in figure 11(a) and (c). In (a) we plot the TSM data (extrapolated) against the theoretical curve expected from eq. (2.26). In (c) we plot, on a logarithm scale, the absolute value of the differences between the measured mass and the same theoretical value. For comparison here, we show the differences of the perturbative expansions of the exact mass formula, truncated at various orders with their exact counterpart. This gives an idea of the relative precision of TSM with respect to a perturbative expansion. We see that for intermediate couplings, the TSM outperforms a 5-loop (up to and including O(b 10 )) perturbative expansion of the mass, coming close to 6-loop accuracy. Approaching the self-dual point, the region of viable TSM data shrinks to a region in M ShG R where the exponential corrections become relevant and therefore the standard TSM methods are not available. (Of course, in the sinh-Gordon model everything is supposedly known about these exponential corrections, but for now we intentionally opt for neglecting any a priori knowledge on the integrability of the model.) Vacuum energy density. Let us now turn our attention to the vacuum or bulk energy density, E 0 . The measurement of this quantity proceeds by measuring the slope of the ground state energy E 0 (R) ≈ E 0 R. For small R, this function enjoys a conformal R −1 dependence up to logarithms, and is monotonically increasing.  for larger volumes, truncation errors are expected to dominate. Therefore the best first approximation to E 0 is the minimum of the numerical derivative of E 0 (R). In a general field theory, the leading exponential (Lüscher) correction to the ground state energy is of the form Substituting the mass measured previously and subtracting this correction from the numerical ground state energy improves the precision. The results as a function of b are shown in figure 11(b) and (d) in the same fashion as presented in (a) and (c) of this same figure for the physical mass. Like with figure 11(c), we show the differences between a finite order Feynman perturbative computation and the exact value for orders b 4 through b 14 . Note that the convergence radius of the series for E 0 is only b max = 1/ √ 2 as opposed to 1 for the mass, M ShG . 3 We thus only show the perturbative curves up to this point in b.

JHEP01(2021)014
Conventional TSM is able to measure the bulk energy with an error of 10 −3 even in this strongly coupled region.

S-matrix.
Finally, let us consider the measurement of the S-matrix from TSM data. For asymptotically large volumes, two-particle states with zero overall momentum (each particle having rapidity ±θ) are quantized by the requirement that the multi-particle wavefunction be one-valued on the cylinder which, after taking the logarithm, provides the Bethe-Yang quantization condition where we have introduced the phase shift S(θ) = e iδ(θ) . Eq. (5.4) is a quantization condition which determines the rapidity θ. In fact, it is the large-volume limit of the TBA quantization condition (2.61) for the state {I 1 = −n, I 2 = n}. Once this quantity is known, we have access to the energy of the two-particle state since, up to exponential corrections, the energy is a sum of one-particle terms We focus on the lowest energy two-particle states in the zero-momentum sector by taking n = 0, 1 in (5.4). Numerically, for large enough volumes, this corresponds to the fourth lowest energy level. In this domain, we express the rapidity θ in terms of the energy difference between the two-particle state and the vacuum: Thus we can directly measure the phase shift appearing in eq. (5.4). The result extracted from the n = 0 and n = 1 two-particle states is shown on figure 12(a). In the following, we have assumed that the S-matrix indeed consists of a single CDDlike factor, namely but we have treated the quantity B in the S-matrix amplitude as a parameter to be fitted to the numerical data. To perform this fitting, we have only utilized our n = 1 TSM data in region of R where no level crossings occur. The numerical phase shift obtained in this way is more robust than that coming from the n = 0 state. We estimate the phase shift by two methods. In the first method, we increase the parameter B until one of the numerically determined phases coincide with the theoretical curve. The value of B at this point is the estimate. In the second method, we instead look at the two largest rapidities available from the n = 1 data. We then find the value of B for which eq. (5.7) best approximates the values of δ(θ) at these rapidities. The difference of the results of the two methods is considered to be the error of the measurement. We will return to the measurement of the above quantities by an alternative method in section 6.5.

One-point functions
TSMs can be used to measure one-point functions (either on the vacuum or on exited states) by sandwiching Schrödinger-picture operators between the numerically obtained eigenstates. Since the method directly uses the eigenvectors, the resulting precision is inevitably more limited as compared to the energy spectrum. Nevertheless, it is possible to compare the numerical estimate of the VEVs of the exponential operators to their theoretical FLZZ formula eq. (2.33) in a wide region.
Beyond the usual b and R-dependence, the convergence of the VEV of the exponential operators, e aϕ , heavily depends on the exponent a of the operator. We have again tried to control the cutoff dependence of the one-point functions by power-law fits on the chiral cutoff N c . The functional form eq. (5.8) is helpful as long as b is small enough. Generally the cutoff extrapolation becomes less stable as the coupling or the volume is increased. For larger couplings, we have found it advantageous to use a sum of two power laws as a fitting function: We show the dimensionless quantity G(a) as a function of a in figure 13, for two different couplings and at R = 6. For small to moderate a, the extrapolated numerics agrees with the FLZZ formula to at least to 1%. Regardless of the coupling, the error (and the cutoff dependence) always becomes significant before reaching the first pair of zeros (located at the Seiberg bounds ±Q/2) of the analytic VEV formula. For larger couplings, it is possible to achieve higher precision by performing the computations at a small volume. However, in this case one needs to take into account the finite volume corrections of the VEV, which are available through the LeClair-Mussardo formula [8] in the form of an infinite series:   where C a n is the connected evaluation of the n-particle diagonal form factor for the operator e aϕ : In figure 14, we show the numerical results for R = 2. Exploiting finite volume corrections extends the availability of TSM estimates of G(a) up to b ≈ 0.8, as long as a Q/2.

Renormalization group improvements
In the previous section, we presented a series of results for various quantities in the ShG model as measured with TSM. We showed in general that as one moves towards the selfdual point, the quality of the results found using TSM deteriorates in comparison to the available exact predictions. This deterioration is largely due to the presence of finite cutoff -31 -

JHEP01(2021)014
effects. There are a standard set of renormalization group-like techniques that are employed in ameliorating the effects of a finite cutoff. We show in this section that these strategies are not practical for the sinh-Gordon model close to the self-dual point. However this failure is instructional and point the way to a better understanding of some of the peculiarities of the model and new strategies to tackle them. The strategies take two forms: analytical and numerical. We discuss the analytic form first.

Analytic renormalization group
In presenting how one can analytically take into account the effects of states above the cutoff, we follow the discussion in ref. [54]. The first step is to divide the Hilbert space, H into two parts: Here H l , the low energy Hilbert space, consists of all states of the form Here note we have expressed H l and H h in a way reflective of the tensor nature of the computational Hilbert space (at least as conceived as that of a free massless non-compact boson). We are also working at a fixed number of zero mode states, N ZM , assuming in effect, that this number of zero mode states leads to completely convergent results (an assumption borne out by our numerics reported in the previous section). In particular the high and low energy parts of the Hilbert space have the same zero mode content. We can thus write our Hamiltonian in the following manner:

JHEP01(2021)014
By eliminating c h from the above set of equations, we have In doing so, we have reformulated the eigenvalue problem in terms of coefficients of states that live in the low energy Hilbert space alone. Now we are studying a Hamiltonian of the form H = H 0 + µ ShG V with V given by We can then expand δH in powers of µ ShG , giving Introducing the (imaginary) time dependence of operators in the interaction picture, we can rewrite eq. (5.18) as From here on we introduce another chiral cutoff Λ > N c . We are going to focus upon the most singular (in Λ) contribution to δH 2 and so drop from V in eq. (5.17) the term proportional to cosh(bφ 0 ). Of course if we were interested in using δH 2 in a quantitative fashion, we would need to include this term. We can readily analyze δH 2 through the use of OPEs. OPEs allow us to take into account the insertion of the partial resolution of identity in eq. (5.20) that involves only the states from the high energy part of Hilbert space, H h . Following the procedure outlined in [24], the matrix elements of δH 2 satisfy where we have defined and Λ → ∞. In eq. (5.21), the states a, b are drawn from H l , δ Pa,P b enforces momenta conservation, the sum ϕ runs over all fields, φ (of chiral dimension ∆ φ ), that appear in -33 -
: e bφ(x 1 ,τ 1 ) :: where σ = ±1. As discussed in ref. [24], we obtain the sum n>Nc appearing in eq. (5.21) by expanding the term |z 1 − z 2 | −4σb 2 that appears in the OPE of the oscillator part of the fields into a Taylor series in |z 1 /z 2 |. We then only keep the terms in the series at order N c + 1 and above, up to Λ. We now focus on the part of δH 2 involving the operator cosh(2bφ): We can see that this term diverges as b 2 → 1/2 and that furthermore for b 2 > 1/2, the correction tends to ∞ as the chiral cutoff, Λ, tends to ∞. This means any strategy to compute corrections to TSM results perturbatively due to states coming from above the cutoff fails for values of b close to the self-dual point. This result is actually worse than the second order result implies. At the third order, we can again use OPEs and find that the most singular third order contribution to δH 3 goes as Here we see the third order term has a pathological dependence on Λ when b 2 > 1/3, even further away from the self-dual point. We can continue this to n-th order, finding Here we see that the situation becomes worse and worse as we go to higher and higher perturbative order: at n-th order, the correction diverges as Λ → ∞ for b 2 > 1/n. From this analysis we can see that the perturbative series developed here is essentially a small-volume expansion in the parameter R 2+2b 2 . This implies that the ground state energy does not have a proper expansion in powers of R 2+2b 2 around the CFT limit R → 0.

JHEP01(2021)014
We also want to remark that the pathologies identified here for the ShG do not apply to its analytically continued cousin, the sine-Gordon model. In the sine-Gordon, this perturbative analysis will give rise to divergences for b 2 > 1/2. However these divergences occur in the identity channel in terms of the OPE of eq. (5.23). This means the most singular part of δH n is proportional to 1 and so leads only to corrections to the energies that are state independent, i.e. energies measured relative to the ground state energy are unaffected. Alternatively one can add a single counterterm to the sine-Gordon Hamiltonian to remove this divergent behavior.
As we have stated, in this section we have focused on the leading order singularity at each in perturbative theory as determined by the OPEs. However it is not difficult to determine the subleading divergences. At order n, one can write (δH n ) ab as is some numerical coefficient. There are further singular contributions from descendants. We thus see that not only the leading terms are problematic. We see that the operators : cosh(kbφ): have diverging coefficients (as Λ → ∞) once k 2 > n + (n − 1)/b 2 .

Numerical renormalization group
In section 5.2.1 we demonstrated that a perturbative analytic renormalization group is not a tool that can be used to take into account the states above the truncation N c . In this section we show that the non-perturbative numerical renormalization group [23], while not beset by pathological divergences, also is challenged for values of b close the self-dual point. The basic idea of the numerical renormalization group (NRG) for the TSM is to adapt the Wilsonian renormalization group invented to attack the Kondo problem to the case at hand. Normally in applying TSM, one introduces a cutoff, here N c , and either does a single exact diagonalization or uses Jacobi-Davidson (JD) methods to obtain the energies. With the NRG, one trades a large single diagonalization for a sequence of smaller exact diagonalizations. This sequence is determined by two parameters N s and ∆. The size of matrices that one diagonalizes in the sequence is (N s + ∆) × (N s + ∆). These parameters should be thought of as variational in nature. In general the larger these parameters are, the closer one gets to reproducing the exact diagonalization result. We will not describe this procedure in further detail here but refer the reader to refs. [23,24].
In figure 15 we present results for the computation of the ground state energy at two different b's. We do so at cutoffs of N c = 14 and N ZM = 24. The Hilbert space size at such cutoffs is 492888. In figure 15 we show the results of the NRG computation for different values of N s = ∆ ranging from 2500 to 10000 (i.e. we are diagonalizing sequences of matrices with size from 5000 × 5000 to 20000 × 20000). We see that at even the largest value of N s = ∆ = 10000 considered, the results are not converged. We thus fit a power law to the evolution of E gs as a function of N s = ∆ and extrapolate the power law to where it would correspond to solving the problem exactly (i.e. finding the low lying eigenenergies of a 492888 × 492888 matrix, corresponding to evaluating the power law fit function at N s = 492888/2). The result is reported in figure 15. We see that for b = 1/ The performance of the NRG close to the self-dual point is considerably worse than for a model like the sine-Gordon model where we see agreement at the 5 significant digit level for NRG block sizes far smaller than those considered here (N s = 1500, ∆ = 500) and without extrapolation (for example compare the results here with table IV of section VI of ref. [24]). This we believe is a manifestation of the slow convergence as a function of N c in the ShG model that we have observed elsewhere in this paper. While obtaining 4 significant digits is usually sufficient (say at b = 1/ √ 2), we are of course interested in further extrapolating our results in N c . These N c -extrapolations turn out to be sensitive to the errors on the order of 10 −3 . This makes the use of NRG-based data, at least for values of b close to 1, problematic.

Quantum mechanical reductions of the sinh-Gordon model
In section 5 we argued that a straightforward implementation of analytical RG improvements is hindered because the small volume expansion of energy levels is not perturbative with respect to the parameter µ ShG R 2+2b 2 . In this section we take a closer look at the small-R UV spectrum. Since the energy of oscillators behave as R −1 in the R → 0 limit,

JHEP01(2021)014
one might expect that the UV behaviour of the spectrum is dominated by the quantum mechanics of the zero mode of the field. Using this quantum mechanical picture, a systematic expansion for certain energy levels (more precisely, their corresponding scaling functions) can be developed in terms of 1 ln(µ ShG R) . Alternatively, one can expand the TBA equations, yielding a similar expansion, but involving IR parameters (the physical mass, M ShG , and the S-matrix parameter B). Using the mass-coupling relation and expressing B in terms of the coupling b, we get another expansion in 1 ln(µ ShG R) , which is however different from the zero mode expansion in subleading orders. As the energies contain an additional 2πR −1 factor relative to the scaling function, the difference between TBA and zero mode energy levels eventually diverge for R → 0 for all b > 0.
In subsection 6.3 we derive an effective potential, partially taking into account the effect of oscillators. On one side, we analytically reproduce the exact expansion up to O(b 12 ), confirming that the oscillators are able to explain the differences in the log-expansion. On the other side, we show that the TSM numerics significantly outperforms even the numerical solution of the complete effective potential. We then use this fact to provide an alternative measurement of the IR parameters from TSM, combining UV numerics with the small-volume expansion of the TBA.

Semiclassical reflection amplitude
In the semiclassical limit b → 0 the small volume behavior of energy levels is dominated by the contribution of the zero mode [47]. For R → 0 the potential walls that the zero mode sees (i.e. the points where the µ ShG R 2+2b 2 cosh(bφ 0 ) potential exceeds 1) are far from one another. Then it is sensible to consider first the quantum mechanical problem of a particle reflecting from a single wall: Introducing the coordinate representation ϕ 0 ≡ x, Π 0 = −i∂ x , it is possible to solve the Schrödinger equation Its general solution is given by modified Bessel functions. Requiring that the wave function vanishes at x → ∞ and evaluating the x → −∞ asymptotics, we can write the relative phases of the left-moving and right-moving wave as where the semi-classical reflection amplitude is defined as This expression is the semi-classical b → 0 limit of the Liouville reflection amplitude (2.46). As the other exponential term is turned on, we can get an approximate quantization condition for the energy levels of the full potential through the quantization condition of the Denoting S sc (P ) = −e iδ(P ) and taking the logarithm, the quantization condition (6.5) reads δ (P ) = n π, n ≥ 1 (6.6) and the branch cuts of δ(P ) are to be chosen such that it is continuous for real P and δ(0) = 0. The equation for the ground state wave number P 0 is then given by δ (P 0 ) = π . (6.7) Making a formal expansion of δ (P ) δ (P ) = δ 1 P + δ 3 P 3 + δ 5 P 5 + . . . (6.8) we can expand the ground state momentum as a function of z = δ −1 1 : P 0 = πz − π 3 δ 3 z 4 − π 5 δ 5 z 6 + . . . , (6.9) where the parameter z reads explicitly Hence, the (semi-classical) ground state energy admits the expansion

JHEP01(2021)014
It is important to notice that the terms of this series contain log (µ ShG ) −1 factors. Therefore, it is not so surprising if a power expansion in µ ShG around µ ShG = 0 turns out to be pathological.

Quantization condition from the UV limit of TBA
A similar quantization condition exists for the exact energy levels and can be obtained from the small-volume expansion of the TBA system [7]. In this section we denote the coupling appearing in TBA byb, to emphasize that this parameter is directly related to the S-matrix parameter B =b 2 1+b 2 , and not (immediately) to the parameter b appearing in the Hamiltonian.
In the small-R limit, the TBA equations decouple into a right-and a left-moving part. To obtain these equations, one first performs a shift in the rapidity variables {θ, ϑ} → {θ ± |ln M R| , ϑ ± |ln M R|}, leading to a pair of volume-independent equations up to O (M R) corrections. One then neglects the O (M R) corrections and reverses the previous rapidity shift. Let us introduce the notation Y ± (θ | {ϑ}) = e − ± (θ|{ϑ}) , where the ± denotes the right-and left-moving solutions. In the following it is advantageous to construct the so-called Q-functions, defined through the pair of functional relations ln (1 + Y ± (θ | {ϑ})) cosh (θ − θ ) + (source terms) . (6.14) Note that the source terms and quantization conditions of the TBA system can be understood as a prescription for the zeros of the function 1 + Y . Correspondingly, Q ± needs to have an analogous set of zeroes to be compatible with eq. (6.13). It was shown in [57] that the Q-functions obtained from the decoupled TBA equations possess the asymptotic form where P is a real parameter and Θ (P ) is an antisymmetric phase (to be obtained below). For small volumes, the asymptotic eq. (6.15) is expected to dominate the θ-dependence over a wide region. The parameter P is quantized by the requirement that Q + and Q − corresponds to the same Y -function Y + = Y − , which is nothing else but the UV limit of the Y (θ | {ϑ}) = e − (θ|{ϑ}) .

JHEP01(2021)014
Let us focus on the zero mode sector defined by restricting the Bethe quantum numbers to be I j = 0, ∀j. Taking into account that ±Q leads to the same Y , we arrive at the condition 2 Θ (P ) = nπ, n ∈ Z . (6.16) Due to the antisymmetry of Θ, it is sufficient to restrict to the cases n ≥ 0 (subsequently we will see that the ground state corresponds to n = 1). In the zero mode sector, energy levels behave in the UV as as can be seen by direct integration (for details, see appendix C of ref. [7]). It is hard to extract the phase Θ (P ) directly from eq. (6.14). However, as shown in [58], there exists an ingenious trick to obtain its explicit expression. To this aim, consider the second-order ODE (6.18) and the pair of solutions defined by their asymptotical behaviour Hence, the Wronskian, W (p), constructed in terms of these solutions satisfies the same set of functional equations as the Q-system, provided that the parameters κ and p are tuned appropriately. Let us focus on the right-moving part. W (p) can be evaluated in a small-κ expansion (using the reflection quantization) and a large-κ expansion (by means of the WKB approximation). It is convenient to parametrize κ as κ = c e θ . Then, comparing the form of the pseudo-energy obtained from the small-κ expansion (θ → −∞) to the asymptotic formula eq. (6.15), we can fix while, from the large-κ expansion, we obtain (6.23) Finally, the phase Θ (P ) is obtained by comparing the small-κ expansion to eq. (6.15) and is given by Notice that if we take advantage of the mass-coupling relation and equateb ≡ b, as it was observed earlier [47,59], the quantization condition eq. (6.16) can be expressed as where S L (P ) is the Liouville reflection amplitude given in eq. (2.46). This is analogous to the semiclassical formula eq. (6.5) but with the miraculously appearing Liouville reflection amplitude which replaces the quantum mechanical amplitude S sc (P ) introduced in subsection 6.1. At the same time, S L (P ) reduces to the semiclassical S sc (P ) in the b → 0 limit. Comparing to the quantization condition (6.6) and assuming continuity of energy eigenvalues as functions of b, it is then natural to exclude n = 0 from the quantization condition eq. (6.16). Repeating the UV expansion of the energy levels using the condition eq. (6.26), we get a modified expansion for the ground state energy Even though the leading term of this expansion coincides with the semiclassical expansion eq. (6.11), the R −1 factor in the front ensures that a difference in any term of the u-expansion leads to a singular discrepancy in the R → 0 limit. This small-volume discrepancy is not trivially accounted for perturbatively. The µ ShG -expansion introduces terms that vanish in the R → 0 limit. On the other hand, in the massive scheme of eq. (A.28), where an expansion in b is natural, the coefficient of the perturbing operator diverges in the R → 0 limit. In the next subsection, we overcome these obstacles by deriving an effective potential which partially takes into account the corrections appearing in eq. (6.27).

An effective quantum mechanical potential
In the previous subsections, we compared the small-volume expansion of the ground-state energy obtained from the zero mode quantum mechanics to the UV expansion of TBA, and found that they differ by R −1 (ln R) −k type terms. It is then an important question whether the oscillator states neglected in the zero mode calculation can account for these inverse logarithmic differences, or is this a sign that we are missing additional terms from the Lagrangian. In any case it is not straightforward to reproduce these terms perturbatively. In the following we derive an effective quantum mechanical potential from the Lagrangian, which provides a more precise description of the UV spectrum of the zero mode subspace, by partially taking into account the effect of oscillators. This is done by JHEP01(2021)014 means of a Bogoliubov transformation applied to the Hamiltonian, involving the oscillators only (keeping the zero mode intact). Following the notations of appendix A, we can start by adding and subtracting an auxiliary quadratic term to the Hamiltonian eq. (3.1), but considering the zero mode ϕ 0 as a free parameter with m eff ≡ m eff (ϕ 0 ) a function that depends on it. Upon transforming to the oscillator eigenbasis of the Hamiltonian given by the first two terms of eq. (6.28), the normalization of the cosh term changes Here :: m eff implies that we are normal ordering w.r.t. to massive (of mass m eff ) oscillator modes). The value of m eff is then given by the requirement that the explicit quadratic term precisely cancels that of the cosh function. This leads to the transcendental equation The effective Hamiltonian H eff ZM is then obtained by dropping all higher order oscillator terms of the cosh interaction. This leads to consisting of the contribution of a Bogoliubov ground state energy plus a correction due to the change of normalization of the cosh term.
Let us now focus on the small volume limit of the spectrum of H eff ZM . The effective mass of this Hamiltonian is then approximated by

JHEP01(2021)014
while the sums S 1 andS 2 admit the small-volume behavior Inserting the expansions eq. (6.31) and eq. (6.32) into the effective Hamiltonian eq. (6.30) and using the coordinate representation, we obtain a Schrödinger equation with the asymptotic form where y is a zero mode wave-function. In writing (6.33) we have introduced the notations The potential V (x → ±∞) = re ±x − 2 e ±2x that we obtain in this way does not possess normalizable eigenfunctions. This is an artefact of expanding the effective mass according to eq. (6.31). However, for small R, the potential builds up a flat plateau around x = 0, which is bounded by a large peak on either side. WKB analysis predicts that tunneling through the peaks can be neglected as long as P In this domain, the amplitude of reflection from the peaks can be approximated as . (6.35) For b 1, the absolute value of S 2 (P ) is essentially 1 over a wide interval of P . To now determine the energy levels, we can once more follow the tactic outlined in subsection 6.1, writing a quantization condition for the phase shift δ eff ≡ −i ln (−S 2 (P )) = nπ. Expanding this in powers of P and using (6.3), we find the small-volume expansion: 36) where Notice that while the difference between the zero mode expansion eq. (6.11) and the smallvolume expansion eq. (6.27) from the exact UV quantization condition is order O b 8 , the difference between eq. (6.36) and eq. (6.27) is instead only order O b 12 .

Numerical performance in the UV
We have now at hand six different ways to estimate the spectrum of the ShG model in the UV small-R limit: i) raw TSM data; ii) extrapolated TSM data; iii) diagonalization of the zero mode Hamiltonian, H ZM (i.e. eq. (4.3)); iv) reflection quantization using the semiclassical reflection amplitude; v) reflection quantization using the full Liouville reflection amplitude; and finally vi) diagonalization of the effective zero mode Hamiltonian, H eff ZM , derived in section 6.3. In this subsection we make an effort to compare the numerical accuracy of these different approaches. We present our results in figure 17 for four different values of b. Here the ground state energy corresponding to the Hamiltonian in eq. (6.30) was obtained numerically and compared to the raw zero mode potential using a real-space basis of size 16000 -see eq. (4.6). Let us summarise what we learn from these plots.
• The UV limit of TBA is indeed different from that calculated from the zero mode Hamiltonian H ZM alone, apart in the b → 0 limit.
• In precisely the b → 0 limit, the validity of the reflection quantization method shrinks -44 -

JHEP01(2021)014
to extremely small volumes. This is intuitive at least in the semiclassical case as the potential increases relatively mildly around the minimum, so neglecting the overlap of the two exponentials in the middle is not well-founded.
• For small couplings, the effective potential eq. (6.30) efficiently accounts for the discrepancy between the eigenvalues of H ZM and the exact TBA energies.
• As the coupling b is increased, higher order oscillator terms in µ ShG not included in eq. (6.30) become significant for arbitrarily small volumes and our expression for the effective potential leads to inaccurate results.
• In this same regime of large b (but always b < 1), the accuracy of the ground state energy as obtained from reflection quantization extends to ever larger volumes well into the IR regime of the model.
• However, TSMs (especially after power-law extrapolation) are able to reproduce the TBA result surprisingly well, even for larger b. This provides an important confirmation of the efficiency of TSMs as applied to the sinh-Gordon model.

Extracting M ShG , E 0 and B from UV spectrum
We have shown that the TSM is able to reproduce the UV limit of TBA equations remarkably well. It is thus reasonable to take advantage of this fact and to try to extract infinite volume parameters combining TSM numerics with the quantization condition eq. (6.16).
As the exact quantization condition is expressed in terms of the IR parameters, M ShG , E 0 , and B, it is possible to fit these parameters using TSM data. To do so we use the lowest two energy levels (the vacuum and the one-particle state) and two values of the volume (mR = 0.1 and mR = 0.2), and then we minimize the function as function of M ShG , E 0 and B. In figure 18 we compare the mass and bulk energy obtained in this way to the standard TSM methods discussed in section 5. It is worth stressing that this method produces an estimate for the mass that is an order of magnitude better than the extrapolated TSM data reported in figure 11 (see panel (a) of figure 18). The determination of E 0 as shown on figure 18(b) provides a precision comparable to the previous analysis. The real power of the UV method is revealed when we consider the measurement of the S-matrix parameter B. Here the improvement of the error generally exceeds two orders of magnitude (see figure 19).

What have we learned?
In this section we have put the UV behavior of sinh-Gordon theory under scrutiny. After pointing out the incompleteness of the zero mode in describing the small-volume limit of the ground state energy for any finite b, we have derived an effective quantum mechanical   potential. This effective potential accounts for the above discrepancies up to an O(b 12 ) error, by partially taking into account the oscillators neglected from the zero mode.
Comparing the zero mode, the effective zero mode, and the TBA-Liouville reflection quantization results to TSM and exact TBA numerics, we realized that the effective potential significantly outperforms the zero mode. On the other hand, numerical TSM (especially after extrapolation) significantly outperforms the effective potential, providing an even more effective incorporation of the oscillators.
At the same time we noticed that the validity of the reflection quantization approximation extends to larger volumes as the coupling b is increased. This lead to a more precise, UV-based measurement of the mass and the S-matrix parameter B, and even provided a consistency check for the energy density E 0 .
We emphasize that both the S-matrix parameter B and (up to the mass scale) the Liouville reflection amplitude are self-dual quantities, which are successfully reproduced -46 -

JHEP01(2021)014
directly from the Lagrangian, which is not in any sense manifestly self-dual. In this sense, the 'self-duality' of the Lagrangian ShG model (i.e. the dependence of the model on coupling constant as encoded in the expression as b 2 /(1+b 2 )) is confirmed up to the precision of TSM. This result, together with the partial analytical reproduction of the reflection amplitude, indicate that the discrepancy in the ground state energy as derived in the UV expansion coming from the semi-classical and Liouville quantizations is solely due to the effect of oscillators. Importantly, no extra terms (like cosh(b −1 φ)) in the Lagrangian are necessary to account for it.
However, note that all the above measurements were done below the self-dual point. As the self-dual point is traversed, the problems arising from the mass-coupling relation are inherited by the Liouville quantization condition as well. In the next section, we develop a supra-Borel resummation technique to motivate the power-law fit and to reach a better understanding of the cutoff dependence in applying TSMs to the ShG model. In section 8, we actually provide an argument that if the Lagrangian with a finite positive µ ShG parameter defines a meaningful theory at all for b > 1, this theory should actually be massless.

Supra-Borel resummation
In section 5.2.1 we demonstrated that the leading order corrections, δE(N c , Λ), to the energies, E(N c ), obtained at a given cutoff, N c , that arise from taking into account states above the cutoff up to a new cutoff Λ > N c can be expressed as a perturbative series in µ ShG of the form where a n is constant independent of N c , Λ.
The first thing to note about this series is that no matter how small b is, terms in the series for sufficiently large n will diverge as Λ → ∞. The second thing to stress about this quantity is that δE(N c , Λ) is ultimately a finite quantity with a well-defined limit as n → ∞ for fixed Λ. This follows as all UV singularities in the theory, at least for b < 1, have been removed through the normal ordering of vertex operators. The divergences that one sees in the perturbative representation of δE(N c , Λ) found in eq. (7.1) are not reflective of a need to add additional counterterms to the theory in order to make it UV finite. The last point is evident by simply choosing a Λ slightly larger than N c and performing TSM numerics up to the cutoff Λ. In this case we get the exact (and finite) δE Now in other models, such counterterms can be a necessity. For example, in the sine-Gordon model, when β 2 > 1/2 (β BKT = 1), a counterterm proportional to the identity operator is required for the finiteness as Λ → ∞. Here it is important to emphasize that the need for this counterterm is readily observable in TSM results for sine-Gordon. When one studies the sine-Gordon model for β 2 > 1/2, one finds that the ground state energy computed by TSM has a marked cutoff dependence in N c . This is not the case for sinh-Gordon. In all of our TSM studies for the b < 1 ShG model, we observe no such -47 -JHEP01(2021)014 cutoff dependence, giving numerical 'proof' to the fact that indeed no new counterterms are required in this theory.
So then, what are we to make of divergences that appear in each of the terms of eq. (7.1)? These divergences are, in fact, an artefact. They are a result of computing δE(N c , Λ) by doing perturbation theory in µ ShG . Because the µ ShG cosh(bφ) potential is (strongly) unbounded from below if µ ShG < 0, the radius of convergence of any perturbation theory in µ ShG is 0. This is reflected in the TBA energy in µ ShG being a series, not in µ ShG but in 1 log µ ShG . It is thus not surprising to find an ill-defined perturbative series around µ ShG = 0, even though the TBA energies are finite. Thus for any N c , this series must be treated as asymptotic. The question then becomes whether it can be resummed. We know E(N c ) itself is finite. But can we resum the series and explicitly demonstrate this? 4 Now this series is not Borel resummable in the standard sense, namely the coefficients at order n are diverging more quickly than n! (in fact, they are diverging as e 2b 2 log(Λ)(n 2 −n) ). We will show in this section that nonetheless it is possible to resum this series in a meaningful way for any N c < Λ and obtain a finite result. From this resummation we will obtain some understanding of the slow convergence of TSM results in N c as well as a partial justification of the power law scaling that we observed in section 5.
Before doing so we consider a toy example. It is a bit unusual from a standard field theory perspective to perform a resummation of a perturbative series at fixed cutoff, N c . Normally one would take the cutoff to ∞ first and then perform the resummation. But, again, here the theory as is is UV finite and, moreover, our purpose here is to connect to our TSM results which are always computed at fixed N c . The toy example with which we will discuss in the next section shows that there are simple functions which are finite, but whose perturbative expansions suffer from the same behavior in the presence of a cutoff as our quantity of interest, δE(N c , Λ).

Toy example
Here we show that it is easy to find functions for which when one performs a perturbative expansion about a non-analytic point, one can obtain pathological behavior in a cutoff used to regulate each term in the expansion. Consider the function, f (µ, ), defined as Here is a regulator which if set to zero (analogous to taking N c → ∞), leaves f finite, i.e.
This function is clearly well behaved for any µ ≥ 0, b > 0. However if we now naively expand about µ = 0, we have

JHEP01(2021)014
We see that while the original form of the function f (µ, ) was finite as → 0, the expansion in the coupling µ leads to a series where terms beyond a certain order in k are all divergent as the cutoff goes to 0. This series is thus asymptotic. The k-th order term here goes as e k log(1/ ) (so more slowly than k!) and here at least the series at any fixed should be resummable via standard Borel resummation techniques. Now, with this toy example in mind, let us turn to the real problem at hand.

Minimal resummation
In the following, we are going to take a slightly different approach. Instead of considering eq. (7.1) (with N c fixed and Λ → ∞), we are going to employ a perturbative expansion in µ ShG of the ground state energy, E(N c ) (with N c → ∞), that takes into account all contributions up to cutoff N c . 5 Here the advantage will be that we will have a specific form for the a n 's of eq. (7.1), the coefficients at each order of the expansion that are independent of N c . This expansion, as we demonstrate in appendix C using certain assumptions, has the following form: is the dimensionless ShG coupling and a n , c are dimensionless constants independent of R, N c .
Despite this series not being Borel resummable, it admits what we call a supra-Borel resummability, something G. Hardy in his classical treatise on divergent series [60], termed resummation via moment constant methods. In this procedure, one supposes that one has the asymptotic series S(z) = k a k z k . In order to resum it, one chooses a function ρ and then defines its moments on R + as The series, S(z), is then said to be r − ρ resummable if converges in some neighbourhood of t = 0 and B(t) has an analytic continuation to a neighbourhood of the positive real axis. If the integral, is convergent for a z = 0, then the function g(z) exists, is analytic in some domain −∞ < Re(log(z)) < c 0 , and has the asymptotic Taylor series S(z). S(z) can then be identified as the resummation of S(z) in this same domain. 5 Of course δE(Nc, Λ) can be written here as E(Λ) − E(Nc).

JHEP01(2021)014
For the case at hand we choose the function ρ(t) as 16γ . (7.8) As discussed in ref. [61] this choice arises generically in problems with exponential potentials. Ref. [61] discusses the technical conditions needed for resummation for this choice of ρ(t). Assuming for the moment that they are met, we can resum our expression for the ground state energy and rewrite it in the form: In appendix C, we show that in fact a n = (−α) n with α > 0. Assuming this, we can then write B(t) as (7.10) and so is analytic along all of the positive real axis. Now that we have resummed our asymptotic series, let us investigate its properties as a function of N c . It is not a priori obvious that this resummation will necessarily lead to a sensible result, i.e. something that is finite as N c → ∞. But nonetheless it does. To investigate the asymptotics, we write the integral expression for E in a suggestive form: We can now do a Sommerfeld-type expansion and obtain asymptotics in large N c of the form: where C, C are constants. We note that in the large-N c limit, the Borel 'chemical potential', µ B , can be written as and so |µ B | in this limit is invariant under b → 1/b. We also note that as N c → ∞, the resummed part (i.e. involving terms n ≥ 2) of E vanishes if b < 1, but tends to a finite constant if b > 1.
We thus expect that for large N c that We compare in figure 20 the exponents predicted by this result to those determined from extrapolating TSM data (see figure 8 in section 5). We see that as the self-dual point is approached, the exponent predicted by the asymptotics describes approximately the trend observed from power law fits of the data. It however fails to describe the fitted exponents at small b. Instead these are described well by the exponent coming from the first perturbative correction to the TSM data (see eq. (5.26)).

General supra-Borel approximates
In the previous section, we considered the resummation of a perturbative expansion of the ground state energy. This expansion, discussed in appendix C, and given in eq. (7.5), makes certain approximations to the evaluation of the constants c and a n . How important are these approximations to the final result? We show here that in fact the final answer is sensitive to the exact form of these constants. We begin by developing a simple constraint on the resummation arising from the behaviour of the energy as a function of N c . Let us call the energy of a state without truncation E ∞ . Then we can write i.e. if δE(N c , ∞) is, again, the correction to the energy due to including all states above the truncation N c . This means that the fitting form for the TSM energies is

JHEP01(2021)014
In general we expect E(N c ) approaches E ∞ from above, a consequence of the variational principle. So this means δE(N c , ∞) must be negative. For b > 1 we see from asymptotics, i.e. eq. (7.12), that this is indeed the case. But for b < 1 it is not, and we thus conclude that the simple form computed for E(N c ) in eqs. (7.5) cannot strictly be correct. For the purpose of resummation, we require more precise forms for a n 's and the constant c.
In general it is difficult to compute to high precision the a n 's and so derive an exact form to the supra-Borel transform that we presented in eq. (7.10). One way to resolve this problem heuristically is to suppose the Borel transform has a more general form than in eq. (7.10). The simplest possibility is a Padé form: We will only look at forms of B(t) whose large t asymptotics are B(t → ∞) ∝ t as this form leads to a naive cancellation of the powers of N c in eq. (7.9). With this more general form, we can arrange as we will see in the next section, through an appropriate choice of couplings, that δE(N c , ∞) < 0, (7.18) despite the overall minus sign. To ensure this, we allow the r i 's to be of either sign and so that we have no poles on the positive real axis, we assume the s i 's to be positive. So, for example, with M = 2, we take r 1 > 0, r 2 < 0, and s 1,2 > 0. In principle these free parameters can be fixed by carefully computing the a n 's to finite order.
In our discussion here we have only treated Borel approximates that have no poles along the positive real axis. Whether this is actually correct is open to question. We know that the path integral describing the sinh-Gordon model will have complex saddle points. This suggests that a correct Borel approximate, B(t), might then have poles in R + . In order to make sense of such an occurence (i.e. in order to interpret the integral of eq. (7.9)) would require the use of resurgence theory and the associated identification of Lefshetz thimbles related to complex saddle points [62][63][64][65]. While this possibility is very interesting, it is beyond the scope of the current work. We however plan to pursue it in the future.

Fits and discussion of limitations
In this section we look at select examples where we fit as a function of N c our TSM data to the functional forms suggested by the Borel resummations. We compare these fits to a simpler power law fit. We look at two fits suggested by the Borel resummations. In the first we look at a fit that incorporates only the asymptotics suggested by the minimal resummation discussed in section 7.2:

JHEP01(2021)014
In the second we look at a fit to a form involving an M = 2 supra-Borel transform (7.20) Finally we consider a power law fit: The latter two forms are generically applicable a priori. However the form in eq. (7.19) is not. The Borel resummation 'chemical potential', µ B , has a zero, N c0 as a function of N c . Thus if we are to use this form to fit our numerical TSM data we need to be sure our data runs over values of N c > N c0 . N c0 is given by Becauseμ ShG = R 2+2b 2 µ ShG , we see that at large R the regime where the Borel asymptotic form is expected to apply recedes to ever larger N c . We note in general that if we are in a regime where the Borel asymptotics applies then the exponent from power law fit should approximate that predicted by eq. (7.13). We will see in the examples below for larger values of b that this does not happen.
In figure 21 we show the results of these fits for two different values of b: b = 0.778 and b = 0.919. On the left hand side of the figure (panels a and c), we show the fits over the existing range of the data (for b = 0.778, 5 ≤ N c ≤ 20, and for b = 0.919, 6 ≤ N c ≤ 14) on a linear scale. We see that in all three cases (Borel asymptotic, eq. (7.19), M = 2 Borel, eq. (7.20), and power law, eq. (7.21), the fits appear identical. What is different however is the value of E ∞ inferred from the fits. These differ markedly between the fits. We see in comparison to the exact value (as determined from the TBA), the power law extrapolation, E ∞,power law , is greater than the exact value (i.e. there is an under extrapolation) and the Borel asymptotic fit leads to a value E ∞,Borel asymp. that is less than the exact value (i.e. there is an over extrapolation). The M = 2 Borel fit for these examples seems to perform best (see below, however), leading to an asymptotic value, E ∞,Borel M=2 , that is closest to the exact value. One conclusion that we draw from this is that the values of N c at which our data is computed are decidedly not in the large-N c asymptotic regime. Our data does obey a power law form (as our fits indicate), but this power law is different from the one derived from the Borel asymptotics (i.e. eq. (7.13)). In particular the power law being obeyed by the data at N c in the range (10,20) is greater than the one indicated by eq. (7.13). This is consistent with our finding that the power law fit under extrapolates and leads to a value of E ∞ that is too large while the Borel asymptotics do the opposite. That the M = 2 Borel form does better is thus not surprising. However this should not be taken (at all) as the final word of which fitting form to use. Certainly we have made no real attempt to explore the space of supra-Borel transforms, B(t), and leave this to later -53 -JHEP01(2021)014    work. As we have already stated, a complete analysis will require the incorporation of any putative complex saddle points of the path integral for the ShG model and the attendant use of resurgence theory [64,65].
To give the reader a better picture of the extrapolations that are being performed by the fits, we have also plotted the fits on a semi-log scale (the right hand side of figure 21). Here we explicitly show the extrapolations to large N c , N c ∼ e 20 . We see that the different fits only begin to diverge from one another for N c far in excess of our ability to perform numerical TSM computations. We also see the necessity for the extrapolations. Even the most rapidly converging of the fits, the power law fits, only see convergence for N c e 10 at the values of b that we are considering in figure 21. We do note that the convergence at b = 0.778 is much more rapid than at b = 0.919, something we already expected from the power law inferred from the Borel asymptotics eq. (7.13).  Table 1. Fits including an uncertainty analysis. All data is for µ ShG = 0.1.
In the fits described in figure 21, we made no attempt to assign an uncertainty to our fits. In table 1 we attempt to do this. For each of the three different forms, we perform four fits. Each of these four fits is over different subsets of our numerical data. The different subsets are obtained by dropping data points for the lowest values of N c at which we have numerical data. For the b = 0.778 fit, the subsets are formed from data at values of N c corresponding to the sets {5, · · · , 20}, {6, · · · , 20}, {7, · · · , 20}, and {8, · · · , 20}, while for the b = 0.919 fits the N c -subsets are given by {6, · · · , 14}, {7, · · · , 14}, {8, · · · , 14}, and {9, · · · , 14}. Having done these four fits, we take the average and standard deviation of the obtained values of E ∞ and report these in table 1. In this table we provide an analysis for one more pair of (b, R), (0.919, 3.0), beyond those considered in figure 21.

What have we learned?
We have shown in this section that the divergences observed in the perturbative corrections to the TSM can be made sense of. In particular, we have shown that these corrections can be both computed at all orders in the coupling µ ShG and resummed with a resulting finite value via a supra-Borel procedure. We note that the ShG model is not the only theory where resummation is needed. The ability to do so, for example, might be useful for studies of the Φ 4 (and other polynomial) theories. To date at least partial (useful) resummations of the expansion eq. (5.18) are available and tested in the context of 1 + 1d Φ 4 theory in refs. [29,30].
We have shown the functional form of the resummed corrections gives us important insights into the large N c asymptotics of the TSM. It provides a partial explanation at values of b close to the self-dual point for the power law fit that we have observed (see figure 20). It also provides a guide to how to extrapolate TSM results to N c = ∞. In the section we considered two fitting functions suggested by the supra-Borel resummation to use with our TSM data. In the first, we employed only the asymptotic form (valid for large N c ) suggested by the resummation. In the second we used a Padé form for the Borel approximant. Unsurprisingly the Padé form outperformed the asymptotic form. The asymptotic form is likely only really appropriate for values of N c beyond that we have performed TSM computations.
In presenting this resummation, we have really only scratched the surface. We are fairly certain that there are more appropriate Borel approximates than used here that are better tailored to the properties of the ShG model. We view this an important open problem presented by this work.

JHEP01(2021)014 8 Form factor truncated spectrum method
In the previous sections we have focused on numerical TSM data arrived at from a computational basis based on a non-compact free massless boson and its zero mode. We have discussed in detail the strong dependence of the data from the cutoff, E c , in particular approaching the self-dual point b = 1. In light of these results, it is natural to ask the following question: Is there a starting point H 0 for the TSM much 'closer' to the model of interest?
In this context, 'closer' means having a set of eigenstates such that TSMs is better suited to approximate the spectrum of the perturbed theory. For example, instead of separating the zero mode, we might have started with the Hamiltonian (A.28) and so used a computation basis built on a standard massive Fock basis. Having considered this, a more radical possibility suggests itself. This, however, comes at the cost of abandoning the comfort provided by the free nature of H 0 .

Using the ShG basis as a computational basis
The main idea is as follows. Our target is a ShG model at coupling b 1 and we choose as our starting point a ShG model at a different coupling b 0 , setting . While such a starting point does suffer from the critique that we are using as input information that is itself non-trivial and is meant to be, in part, verified by our TSM studies, we believe that there are good reasons to explore this path: • As we will see soon, this method provides a set of a highly nontrivial self-consistency checks.
• It provides useful insights regarding the validity of the exact VEV formulae in certain domains.
• It is a non-traditional starting point for the application of the TSM to interacting massive field theories 6 and insight gained here will be useful for the study of deformations of other massive integrable field theories.
In the approach taking in here, we write the Hamiltonian as: where the solvable piece is given by where H 0 is the Hamiltonian of the non-compact c = 1 boson, while the perturbation is given by Certainly the TSM has a long history of employing an H0 that is interacting. However H0 in such cases is an interacting massless CFT [20,22,31,34,35,53,[66][67][68]].

JHEP01(2021)014
In order to proceed with the diagonalization of the Hamiltonian (8.1) on a cylinder of finite circumference R we need two highly nontrivial sets of data: • Finite volume energy levels of H Fortunately, both pieces of data are available in the ShG model (as well as a number of other integrable models), at least up to exponentially small corrections in the volume. Indeed, according to Lüscher's principle, all information about finite volume quantities are encoded in infinite volume data, like the physical mass M and the S-matrix. Indeed, the exact energies of all finite volume eigenstates can be computed efficiently from the excitedstate thermodynamical Bethe ansatz (see section 2.4). Formulae for the matrix elements, at least up to exponential small corrections [42,69], are also available. In the special case of diagonal matrix elements, one can use a generalization of the Leclair-Mussardo formula [70] to systematically compute exponential corrections. In the following, we are going work mostly in the regime 5 ≤ M ShG R ≤ 15, where exponential corrections should be small and can be neglected.

Finite volume energies
We have already discussed in section 2.4 the excited state TBA that provides the finite volume energies of H 0 . However we revisit this discussion to establish notation that is needed to discuss the finite volume matrix elements. The finite energy eigenstates can be described by a finite number n of particles with momenta p j = M sinh θ j quantized due to finite volume R. We can assign an integer quantum number I j to each momenta. The Bethe-Yang equations essentially impose the one-valuedness of the quantum mechanical multi-particle wave-function in presence of the purely elastic two-particle scattering S (θ) = e iδ(θ) , is given by the set I j , j ∈ {1, . . . , n}. The system eq. (8.4) is solved for the set of rapidities θ j giving the total energy and momentum of the state as:

Finite volume matrix elements
What we have not considered so far are matrix elements in finite volume. Here the Pozsgay-Takács formulae [71,72] come to our aid. For completely non-diagonal form factors, where all rapidities are pairwise different between the bra and the ket states, the formula expresses -57 -

JHEP01(2021)014
the finite volume form factors in terms of infinite volume form factors: [2,3] (see appendix D). Here N k is defined as while ρ k is given in terms of a Gaudin determinant: The normalization coefficients of (8.7) make all matrix elements (8.6) real (apart from the e iM x factors).
We still need to treat the case of when rapidities in matrix elements happen to coincide. In such a case, the infinite volume form factors have singularities that require regulation. This requires then significantly more complicated formulae than eq. (8.6). Coinciding rapidities happen in two different circumstances. In the first the bra and ket states are the same state -so called diagonal form factor). The second possibility is that the two different states have an odd number of particles 2n + 1 and 2m + 1, respectively, with quantization numbers In the latter case, a particle of exactly zero momentum is present in both states. There is, however, a simple trick [73] to circumvent the complications in numerical computations by exploiting the factorization property of exponential operators. Introducing an auxiliary particle with rapidity θ Λ → ∞, and where Λ is the momentum quantization number of the auxiliary particle. Choosing instead a large but finite Λ, the equality is not exact, but the error can be made arbitrarily small -58 -

JHEP01(2021)014
by increasing Λ. In practice we need to solve the set of BYE (8.4) twice: once for k + 1 particles involving the extra quantization number Λ, and once for the quantization number set {Ĩ j } l j=1 . Thus we obtain rapidities of the form θ i = ϑ i + i , where i are small corrections depending on all the I i and Λ. These are then put back into the left hand side of eq. (8.11), yielding the sought for matrix element.

Sanity checks
The setup of the TSM in the basis of the ShG model itself allows us to make some nontrivial consistency checks. The first such check is offered by an infinitesimal change in the coupling constant b 1 = b 0 + δb.

Comparison with form factor perturbation theory
Denote µ 0 = µ(b) and µ 1 = µ(b + δb). Let us calculate the first order correction of the mass of the particle in form factor perturbation theory [74]. It is convenient to keep the system in a finite (very large) volume R. Here the volume factors from spatial integration and the Hilbert state norms cancel. The first perturbative correction to the mass is then: (8.12) where |M ShG denotes the one-particle state of zero momentum and |0 is the vacuum. The diagonal form factor contains a disconnected term which be cancelled by a similar term arising from the VEV of eq. (8.12): Inserting the explicit form of the two-particle form factor calculated using eq. (D.4), we arrive at the relation: where the vacuum expectation values are given by the FLZZ formula (2.33). Therefore, for an infinitesimal δb, we can write down a differential equation for the mass shift equation where

JHEP01(2021)014
The physical mass M ShG , given by the mass-formula (2.26), numerically agrees with the result of integrating equation (8.14). In the special case lim b→0 µ ∝ b −2 , we also have to explicitly specify M ShG (b = 0). The integral appearing in the r.h.s. of (8.15) can be evaluated explicitly in terms of a digamma function Ψ(z): (8.16) A somewhat tedious but straightforward calculation shows that substituting eq. (8.15) into the r.h.s. of eq. (8.14) indeed agrees with the derivative of the exact mass formula (eq. (2.26)), once squared.

Massive boson limit
In the second check, we fix µ(b i ) = m 2 2b 2 i and take the limit b 0 → 0. Here we want to show that we then obtain the Hamiltonian in the massive boson basis. In particular in this limit all matrix elements of the perturbation vertex operators need to approach the limit: with a = ±b 1 and the quantities on the l.h.s. are understood with respect to theory b 0 . In formula eq. (8.17) we denoted the number of particles with quantum number k as n k in the bra vector and m k in the ket vector. Furthermore, matrix elements of vertex operators with a = ±b 0 need to approach those of the operator − m 2 2 : ϕ 2 :. These relations are highly nontrivial, especially for diagonal matrix elements, but we have confirmed numerically their validity on a large number of matrix elements. Note that the b 0 → 0 limit of the form factors is numerically unstable and so the precision often needs to be increased from the usual machine (double) precision.

Implementation and results
Due to the high numerical precision needed, especially in taking diagonal limits, we opted for an implementation in Mathematica. The truncated basis was selected with two cutoffs, a momentum cutoff k max and a limit on the number of particles N max . Figure 22 summarizes the numerical results obtained starting from theories with four different values of b 0 .
To facilitate comparison of different bases, we introduce the following double cutoff. We limit the maximum number of particles N max as well as the sum of Bethe-Yang quantization numbers of each sign: j:I j >0 I j <= k max . Let us remark that enlarging the basis size through the increase of k max only is computationally easier than increasing N max . This is because the latter results in having to evaluate much more complicated expressions   for the matrix elements. However, from numerical tests we deduce that the majority of cutoff dependence is found when changing N max and k max simultaneously. Therefore in the following we fix N max = k max .
In figure 22 we present sample results for the mass gap, starting from four different basis theories (denoted by unfilled circles). Computations were done with a raw cutoff N max = 6. This resulted in Hilbert spaces of dimension 203 in the even particle number sector and 124 in the odd. Numerics were performed at M ShG R = 6 with M ShG being the physical mass of the unperturbed (basis) theory.
In figure 23 we follow the precision of results as the TSM basis is changed from the free -61 -

JHEP01(2021)014
massive boson towards the self-dual point. Here we show explicitly the cutoff dependence of the ground state energy at R = 6. The b 0 = 0 points result from using the free massive Hamiltonian (A.28). The largest cutoff considered for the form factor approach is N max = 8, corresponding to an 1171 dimensional (even sector) truncated Hamiltonian. In the free massive basis, we also show the point corresponding to N max = 10, 12, with basis sizes of 5830 and 25488, respectively. Not surprisingly, the result corresponding to any fixed cutoff improves as we start closer to b 1 . However, it is striking that as we approximate b 0 to b 1 , the power-law exponent of the cutoff dependence apparently worsens. As we have seen previously, when b 1 is infinitesimally close to b 0 , the first order term of the form factor perturbation series dominates. However, as we move b 0 slightly away, corrections from all energy scales play a role.

What have we learned: massless regime for b > 1?
While one can see from the right panel of figure 23 that it is difficult to obtain accurate results (again because of slow convergence in the cutoff) for b 1 close to the self-dual point (even when b 0 itself is close to b = 1), we can use the properties of this massive interacting basis to draw conclusions on the nature of the theory beyond b > 1. We begin, as before, with a starting point b 0 < 1 and aim to describe the theory for a fixed b 1 > 1. We are completely free to choose b 0 < 1 and we so chose it such that it satisfies However for this particular choice of b 1 , the VEV of the vertex operator e bϕ (relative to the b 0 theory) is exactly zero! This in turn means that all matrix elements of e bϕ vanish. Hence according to eq. (8.1), the perturbation consists in simply subtracting the original vertex operators from the original Hamiltonian, which (formally at least) leads to a massless Gaussian model. This argument suggests then that for b > 1, the ShG model, like its SG counterpart is trivial. Because Q(b 0 ) 2 ≥ 1, ∀b 0 , this argument only works when the target theory satisfies b 1 ≥ 1. We note that perturbing a theory with a vertex operator with b > Q 2 is meaningless as the vacuum expectation value becomes negative, thus (formally) making the spectrum unbounded. A plausible explanation is that the VEV e aϕ in the Lagrangian formulation is only identical to the FFLZ formula in the domain − Q 2 ≤ a ≤ Q 2 and vanishes identically outside this interval. We note that this argument then requires us to use the FFLZ formula at the boundary of its validity. These arguments are certainly heuristic but give nevertheless a hint on the phase of the ShG model in the strong coupling regime b > 1.

Conclusions and future directions
In this paper we have studied in detail the ShG model in finite volume as a function of its coupling constant b. This analysis has explored several important conceptual and numerical features of this model. At the conceptual level, we have seen that the model admits two different formulations, given in terms of (i) a Lagrangian and (ii) a S-matrix. The S-matrix formulation is inherently an infrared theory, giving rise in particular to a thermodynamic -62 -

JHEP01(2021)014
Bethe ansatz equation describing in finite volume the energy of the ground and excited states. Moreover, the S-matrix formulation is ecumenical as to the actual value M ShG of the mass of the physical particle. The Lagrangian formulation, on the other hand, gives us the opportunity to study the theory from several directions. Here the theory can be conceived as a perturbation of a massive free boson or as a perturbation of a conformal field theory, either Gaussian or Liouvillian. While the S-matrix formulation is manifestly invariant under the weak/strong duality transformation, b → 1/b, the Lagrangian formulation has no sign of such a symmetry.
Herein we have used truncated spectrum methods to study in finite volume the ShG model. This proved to be particularly challenging. Indeed, while at small values of the coupling constant b, dynamical quantities of the model -such as physical mass, vacuum expectation values of vertex operators, finite volume energies of the ground state and excited states, the two-body S-matrix -could be accurately determined and were found to coincide with their theoretical predictions, this ceased to be true upon approach to b = 1. There we observed that the values of these quantities started to deviate noticeably from their predicted exact results. Indeed, in the vicinity of the self-dual point b = 1, the TSM data showed a marked sensitivity to the cut-off N c related to the number of states employed by the TSM. This sensitivity is of a different nature from other quantum field theories studied so far by means of TSM, a consequence of the unbounded exponential nature of its interaction.
Understanding and attempting to ameliorate these difficulties have had a series of positive by-products. In the first, we were able to come to a detailed understanding of the small volume region of the theory through exploring the quantum mechanics of the zero mode of the field. This analysis led us, in particular, to derive an effective potential which took into account the effect of the oscillator modes. This permitted a more precise measurement of the infrared parameters from the TSM through combining UV numerics with the small-volume expansion of the TBA.
In a second happy after effect, we were able to extend the usual RG scheme for improving TSM numerical data. Normally these improvements are perturbative in the coupling of the theory, here µ ShG . This series proved however to be pathological in µ ShG in the sense that its coefficients diverge absent a UV cutoff and the coefficients at order n increase more rapidly than n!, making the series not even Borel resummable. We were able to overcome these pathologies using generalized resummation techniques for asymptotic series -here dubbed a supra-Borel resummation. In particular we were able to show that summation of this series led to a finite result upon taking the UV cutoff to infinity and that we could partially explain the observed power law dependence in N c observed in our TSM data. This is the first time that the leading corrections at all orders in this perturbative RG treatment for the TSM have been explicitly summed. This perhaps might be useful for studies of other theories. In particular let us mention [29,30] in the context of the φ 4 model, wherein efficient partial resummations of the series eq. (5.18) are achieved by employing a clever variational Ansatz. We believe that here there is considerable work to be done in terms of better undertanding the appropriate supra-Borel approximates to use and the potential role resurgence theory [62][63][64][65] has to play.

JHEP01(2021)014
It is important to stress that all results obtained to this point indicate that the Lagrangian of the theory provide a faithful definition of the theory for b < 1, even though the Lagrangian does not share the self-duality of the S-matrix. A hint about the ShG model for b > 1 may however come from the third spinoff explored in this paper: treating the ShG at a given coupling b 1 as perturbation of a ShG model at another coupling, b 0 . As explained in section 8.4, for b > 1 there is always the possibility of identifying a point b 0 for which the perturbation is given in terms of a vertex operator V a with a = (b + b −1 )/2. Since the VEV of this operator vanishes, it simultaneously ensures the vanishing of all of the matrix elements involving cosh(bφ) and therefore leads to the remarkable conclusion that the ShG model in the strong coupling regime b > 1 is a free massless theory!
We have already noted that in the strong coupling regime of the ShG model that all exact and analytic formulas for the physical mass and vacuum expectation values (see eqs. (2.26) and (2.33)) have a singularity at b = 1. If one analytically continues these expressions beyond b = 1, they give, in general, complex values which make their physical interpretation challenging. One can take the point of view that for b > 1 the Lagrangian should be ignored and one should rely only on the S-matrix (and its explicit duality) to define the theory. In this way one uses explicitly results for b < 1 to define the theory for b > 1. However this is a tautological way of defining the duality as it leads to no predictions with regards to the mass formula, the VEVs, and the energy levels.
It is worth stressing that the same conclusions may apply as well to all Toda field theories which share with the ShG model all of its basic features, namely an apparent duality of their S-matrix which is absent in their Lagrangian formulation. Using the exact formulae reported in [75] of the Toda field theories for the physical mass, vacuum expectation values, reflection amplitudes of the underlying (generalized) Liouville field theory, one can repeat indeed the analysis done in this paper, applying in particular the form factor basis to argue that the Toda field theories for b > 1 are massless models. Indeed, establishing whether the ShG model and all Toda field theories are not in fact selfdual models as their S-matrix suggests, but on the contrary are massless, as conjectured in [76], is one of the most interesting and important open problems for the future.

A Relating the couplings between the Gaussian and massive formulations of the ShG model
In this appendix we derive the relation shown in eq. (2.19), connecting the coupling µ ShG in the Gaussian formulation of the ShG model with the mass m that appears in using Feynman perturbation theory in b. We do so by recasting the finite volume Hamiltonian of the ShG theory where normal ordering is done with respect to a massless basis to a Hamiltonian where the normal ordering is done with respect to a massive oscillator basis. We use the form of the massive basis form (A.28) in section 8.3 where we compare it to the limit of form factors computed in a massive interacting basis.

A.1 Massless c = 1 boson
We first recall how the Hamiltonian of the ShG model on the cylinder as a perturbed Gaussian theory is derived from the theory defined on the plane. This will inform how we derive the ShG model on the cylinder using massive oscillators and ultimate how we arrive at eq. (2.19). Consider the c = 1 boson on the plane. The stress-energy tensor is given by In the main body of the paper we fixed g = 1 8π . In the following we keep this normalization parameter explicit. The generator of dilatations is defined in terms of the stress-energy tensor as where the contour C may be chosen as the positively directed unit circle around the origin.
If we now map the theory to a cylinder of circumference R, where τ is (Euclidean) time and x is the space coordinate, the dilaton operator becomes the Hamiltonian on the cylinder: If we use the expansion of the field in terms of modes satisfying canonical commutation relations 1 |n| a n e iknx + a † n e −iknx , (A.6) -65 -

JHEP01(2021)014
where k n = 2πnR −1 , the Hamiltonian eq. (A.5) can be written as 2π |n| R a † n a n . (A.7) While the − π 6R term can be thought of as arising from the anomalous properties of the stress-energy tensor under conformal transformation, it can also be seen as a result of bringing the oscillator modes into normal order. This will be important going forward.
We end this subsection by noting that under the mapping eq. (A.4), normal ordered vertex operators on the plane and on the cylinder are related by

A.2 Free massive boson as a perturbed c = 1 boson
We now consider the massive boson as a perturbation of the massless one. To this end we define the perturbed dilatation on the plane as

A.3 Expressing the sinh-Gordon model in terms of a massive oscillator basis
Here we finally provide the derivation of eq. (2.19). Our starting point again is the perturbed dilatation operator that defines the sinh-Gordon theory on the plane to be: Mapped to the cylinder, this gives rise to the Hamiltonian The spatial integration can easily be carried out at the cost of imposing momentum conservation explicitly, symbolized by the presence of δ P . Because we are interested in making contact with the analysis of the sinh-Gordon as a massive perturbed boson, we add and subtract at the same time an auxiliary mass term We now perform the same Bogoliubov transformation on the remaining terms of eq. (A.23): Here the normal ordering :: m indicates the normal ordering is being done w.r.t. the massive oscillator modes. The sum appearing in the above is evaluated in appendix B.1 to be where ρ (x) is defined as .

JHEP01(2021)014
After transformation, our entire Hamiltonian appears as We are now in a position to establish eq. (2.19). Using the large R limit of ρ, lim x→∞ ρ (x) = 0, we choose m so that the quadratic term in the second line of eq. (A.28) vanishes, resulting in . (A.29)

B Calculating S 1 and S 2
In this appendix we provide integral representations for the sums S 1 and S 2 introduced in appendix A.
To evaluate this sum, let us first introduce a cutoff Λ = 2πN R , N 1, and then add and subtract an auxiliary integral term: (B.1) In the first term on the r.h.s. we have rewritten the sum so that it includes the n = 0 term. Using the definition of the Euler-Mascheroni constant and the explicit expression of the integral we can immediately evaluate the second term on the r.h.s. of eq. (B.1): The first term of eq. (B.1) can be evaluated in the same way that Matsubara sums are: (B.5) -69 -

JHEP01(2021)014
Here the contour C consists of distinct, small circles around all poles between p = −Λ and p = +Λ on the real axis. We deform this disjoint contour to form two straight vertical line sections running slightly above and below the real axis. The contour in the lower-half plane is then combined with the first integral in eq. (B.1). At this point the limit Λ → ∞ can be taken and the two contours tightened around the branch cuts of the square roots.
In this way we obtain the representation, .
Collecting terms, we see that where now the contour C encloses all poles (in the positive direction) on the real axis except for the one at the origin. The function s 2 (p) has two overlapping pairs of branch cuts all along the imaginary axis of the complex p plane. The contour C can be deformed in the first step into two connected contours at each side of the branch cut, starting and ending at infinity and enclosing all poles lying on one half of the real axis. In the second step, this pair of contours is unbent to form two straight vertical lines aligned in the immediate vicinity of the branch cuts. The contour integrals can then be written as a sum of integrals -70 -

JHEP01(2021)014
the usual Rayleigh-Schrödinger series are generated by using the formula iteratively and expanding in µ. Now we are interested in the divergences appearing in the sinh-Gordon model. It can easily seen that these disconnected terms are less divergent. We thus focus on the connected contributions of this series: where V ij denotes the matrix elements of the perturbation w.r.t. the unperturbed theory. The n-th order term in the above can be rewritten as where, if we introduce τ i = t i − t i+1 , we can rewrite the nth order term as an integral over an Euclidean correlator of the unperturbed theory: From now on, the time dependence of any operator is understood in the sense of eq. (C.5).
We now turn to the Hamiltonian of direct concern: H = H 0 + µ ShG R 2π with Here we have added and subtracted a term quadratic in the zero mode allowing us to formulate the perturbation theory about a massive zero mode. This choice allows for some additional explicit steps in the following calculations. Let us begin by evaluating the second order term in this perturbative expansion. The vertex functions of the oscillator part of bosonic field,φ, can be written as
(C.11) Note that a factor 2 arises due to an analogous term with all exponents negative. The zero mode quantum mechanics (C.7) is that of a harmonic oscillator, so the correlator on the left is easily evaluated. We introduce The Euclidean correlator takes the form For b < 1 √ 2 , eq. (C.15) evaluates to a finite value. For larger couplings, we need to introduce a chiral cutoff N c . In this parameter domain, we can approximate the sum by an integral, and using the asymptotic eq. (C.10), we obtain δE (2) Rωϕ 0 The ω −1/2 ϕ 0 terms make an expansion singular around ω ϕ 0 = 0. On the other hand, having in mind b > 1 √ 2 and a fixed ω ϕ 0 > 0, we argue that the corrections due to the double exponential are subleading for large N c . Let us expand the outer exponential into a Taylor series. Then we need to perform simple exponential time integrals, which yield 1 m -like terms. Rescaling m → N cm and factoring out N c leaves an explicit kω ϕ 0 N −1 c term in the denominator (k is the order of the Taylor expansion term). Even though k eventually -73 -

JHEP01(2021)014
becomes comparable to N c , the corresponding term is multiplied by an overall 1/k! factor, which renders it negligible. What remains is an effective e b 2 4π Rωϕ 0 constant multiplier δE (2) Note that the spurious singularity at b = 1 √ 2 is due to the asymptotic approximation and is unphysical.
In eq. (C.22) the time dependence of the zero mode correlators were neglected by an analogous argument to that of the δE (2) 0 case. The cutoff N c can be scaled out by transforming to new variables p ij = N cpij , q ij = N cqij , resulting in the final form of the leading order asymptotic This is the form that we use as a starting point in section 7.2.

D Form factors of the sinh-Gordon model
In this appendix we provide the explicit form of infinite volume form factors needed for the implementation of the form factor TSM in section 8. Using the S-matrix (eq. (2.31)) and the integrability of the model, one can compute the exact form factors of the local operators O on the multi-particle states of the theory [2,3]: where e kbφ is the VEV given in eq. (2.33) and F n (θ 1 , . . . , θ n ) is given by with σ (n) s the elementary symmetric polynomials in n variables of total degree s. Form factors of the ShG model were also studied in [77,78].

E Reverse communication protocols and chiral factorization
There are several iterative numerical methods (Lánczos, Arnoldi, Jacobi-Davidson, etc.) to obtain the smallest eigenvalues (and the corresponding eigenvectors) of a Hamiltonian. A common feature of these algorithms is that they do not need all of the Hamiltonian's matrix elements to be stored in memory. What they do require is a routine by which arbitrary vectors are acted upon by the Hamiltonian matrix. We will describe in this section how the tensor product nature of the Hilbert space for the sinh-Gordon model makes this matrix-vector multiplication remarkably efficient.

E.1 Separating zero and oscillator modes
As a demonstration of the numerical efficiency that can be obtained, we begin by factorizing the Hilbert space into (merely) two pieces, one encompassing the zero mode H ZM and one the oscillator modes, H osc = H L ⊗ H R . Let N = dimH, N ZM = dimH ZM , N osc = dimH osc be the sizes of these various Hilbert spaces. We want to compute the action of our Hamiltonian upon a vector |v . Normally this would require N 2 multiplications. We show that this can be reduced by a factor of N ZM . To do so, we write |v in terms of our tensored Hilbert space: with |j 0 a basis vector for H ZM and |j a basis vector of H osc . The Hamiltonian consists of a sum of matrices Φ whose matrix elements respect the tensor product structure: where Φ 0 /Φ act in H 0 /H osc . Applying Φ to |v then leads to the need to evaluate:

JHEP01(2021)014
To do so we first perform the matrix-vector multiplication involving the oscillator modes: This involves N ZM N 2 osc multiplications. We now follow this by the matrix-vector product over the zero mode space: This involves another N 2 ZM N osc multiplications. If we write |u = Φ|v , its components, u I are defined as As is typical in our computations, N ZM N osc , we see that we have reduced the total number of multiplications by an approximate factor of N ZM . The above algorithm is a variant of the Shuffle Algorithm where 'shuffle' refers to the reshuffling of elements of the vector |v into a multi-index tensor v j 0j .

E.2 Exploiting the structure of H osc
The above multiplication algorithm can be further optimized by taking into account the chiral structure of the oscillator Hilbert space. Here we will account for momentum conservation. Up to a chiral cutoff N c , the oscillator Hilbert space in the momentum zero sector takes the form: We can proceed just as in the previous section. We write the vector |v that we want to multiply as |v = × |i 0 0 j 0 | 0 ⊗ |i, n i L j, n j | L ⊗ |i,n i R j,n j | R . (E.13) Applying Φ to |v yields (E.14) We now determine the number of multiplications necessary to evaluate this expression. The first step is to evaluate the action of Φ R . This costs N osc N L N ZM multiplications. Similarly the application of the Φ L matrix involves N osc N R N ZM multiplications. Finally the action of Φ 0 costs N L N R N 2 ZM multiplications. The total number of required multiplications is thus For larger values of N c , the number of needed multiplications can be several hundred times smaller than using the naive matrix-vector multiplication involving (N osc N ZM ) 2 multiplications.
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.