Vacuum energy of the Bukhvostov-Lipatov model

Bukhvostov and Lipatov have shown that weakly interacting instantons and anti-instantons in the $O(3)$ non-linear sigma model in two dimensions are described by an exactly soluble model containing two coupled Dirac fermions. We propose an exact formula for the vacuum energy of the model for twisted boundary conditions, expressing it through a special solution of the classical sinh-Gordon equation. The formula perfectly matches predictions of the standard renormalized perturbation theory at weak couplings as well as the conformal perturbation theory at short distances. Our results also agree with the Bethe ansatz solution of the model. A complete proof the proposed expression for the vacuum energy based on a combination of the Bethe ansatz techniques and the classical inverse scattering transform method is presented in the second part of this work [40].


Introduction
The "instanton calculus" is a common approach for studying the non-perturbative semiclassical effects in gauge theories and sigma models. One of the first and perhaps the best known illustration of this approach is the O(3) Non-Linear Sigma Model (NLSM) in two dimensions, where multi-instanton configurations admit a simple analytic form [1]. It is less known that the O(3) NLSM provides an opportunity to explore a mechanism of exact summation of the instanton configurations in the path integral. In order to explain the purpose of this paper, we start with a brief overview of the main ideas behind this summation.
The instanton contributions in the O(3) NLSM were calculated in a semiclassical approximation in the paper [2]. It was shown that the effect of instantons with positive topological charge can be described in terms of the non-interacting theory of Dirac fermions. Moreover, every instanton has its anti-instanton counterpart with the same action and opposite topological charge. Thus, neglecting the instanton-anti-instanton interaction, one arrives to the theory with two non-interacting fermions. Although the classical equation has no solutions containing both instanton-anti-instanton configurations, such configurations must still be taken into account. In ref. [3] Bukhvostov and Lipatov (BL) have found that the weak instanton-anti-instanton interaction is described by means of a theory of two Dirac fermions, ψ σ (σ = ±), with the Lagrangian (1.1) The perturbative treatment of (1.1) leads to ultraviolet (UV) divergences and requires renormalization. The renormalization can be performed by adding the following counterterms to the Lagrangian which preserve the invariance w.r.t. two independent U(1) rotations ψ ± → e iα ± ψ ± , as well as the permutation ψ + ↔ ψ − : In fact the cancellation of the UV divergences leaves undetermined one of the counterterm couplings. It is possible to use the renormalization scheme where the renormalized mass M, the bare mass M 0 = M + δM and UV cut-off energy scale Λ UV obey the relation where the exponent ν is a renormalization group invariant parameter as well as dimensionless coupling g. For ν = 0 the fermion mass does not require renormalization and the only divergent quantity is the zero point energy. The theory, in a sense, turns out to be UV finite in this case. Then the specific logarithmic divergence of the zero point energy can be interpreted as a "smallinstanton" divergence in the context of O(3) NLSM. Recall, that the standard lattice description of the O(3) sigma model has problems -for example, the lattice topological susceptibility does not obey naive scaling laws. Lüscher has shown [4] that this is because of the so-called "small instantons" -field configurations such as the winding of the O(3)-field around plaquettes of lattice size, giving rise to spurious contribution to quantities related to the zero point energy.
To the best of our knowledge, there is no any indication that the fermionic QFT is integrable for general values of the parameters (g, ν) [5]. However, it is expected to be an integrable theory for ν = 0, which is of prime interest for the problem of instanton summation. The corresponding factorizable scattering theory was proposed in [6], by extending previous results of [7][8][9]. According to the work [6] the spectrum of the model contains a fundamental quadruplet of mass M whose two-particle S-matrix is given by the direct product (−S a 1 ⊗S a 2 ) of two U(1)symmetric solutions of the S-matrix bootstrap. Each of the factors S a coincides with the soliton S-matrix in the quantum sine-Gordon theory with the renormalized coupling constant a. The couplings are not independent but satisfy the condition a 1 + a 2 = 2, so that, without loss of generality, one can set a 1 = 1 − δ and a 2 = 1 + δ with δ ≥ 0. A relation between δ and the fourfermion coupling g is not universal, i.e., depends on regularization procedure involved in the perturbative calculations. Nevertheless, g = πδ 1−δ 2 if one uses the regularization that preserves the underlying U(1)⊗U(1) symmetry of the BL model. Together with the fundamental particles of mass M, there are also bound states whose masses are given by M n = 2M sin πn 2 (1 − δ) , where the integer n run from 1 to an integer part of 1 1−δ . As δ → 1 − , the fermion coupling g approaches infinity g → ∞, and an increasing number of particles with vanishing mass occur in the theory. The theory can also be continued into the strong coupling regime with δ > 1 by means of the bosonization technique. Namely, the fermionic BL model can be equivalently formulated as a theory of two Bose scalars ϕ i governed by the Lagrangian [3] (1.4) The interacting term here can be written as 4µ cosh 2 ϕ 2 , and, hence, the bosonic description is still applicable as δ > 1. As it was pointed out by Al.B. Zamolodchikov (unpublished, see also [6]), the Lagrangian (1.4) with a 1 = 2−a 2 < 0 provides a dual description of the so-called sausage model [10], which is a NLSM whose target space has a geometry of a deformed 2-sphere. As a 1 → −∞ the sausage metric gains the O(3)-invariance and we come back to the O(3) NLSM. Notice that the formal substitution δ ≡ 1 − a 1 = ∞ into the relation g = πδ 1−δ 2 leads to the vanishing fermionic coupling in the initial Lagrangian (1.1).
Putting the theory on a finite segment x 1 ∈ [0, R], one should impose boundary conditions on the fundamental fermion fields. We shall consider the twisted (quasiperiodic) boundary conditions, which preserve the U(1) ⊗ U(1) invariance of the bulk Lagrangian, The pair of real numbers (k + , k − ) labels different sectors of the theory and, therefore, one can address the problem of computing of vacuum energy E k in each sector. Notice that twisted boundary conditions is of special interest for application of resurgence theory to the problem of instanton summation [11]. There is no doubt to say that the above scenario of the instanton summation deserves a detailed quantitative study. Perhaps the simplest question in this respect concerns an exact description of finite volume energy spectrum for the theory (1.4) in both regimes 0 < δ < 1 and δ > 1. In this work we will focus on the perturbative regime 0 < δ < 1, where the fermionic description (1.2) can be applied. We propose an exact formula which expresses the vacuum energies in terms of certain solutions of the classical sinh-Gordon equation. The formula is perfectly matching both the conformal perturbation theory as well as the standard renormalized perturbation theory for the Lagrangian (1.2). The result also agrees with the original coordinate Bethe ansatz solution of ref. [3] and the associated non-linear integral equations derived in [12]. The aim of this paper is to review and further develop all these approaches to facilitate future considerations of the NLSM regime of the theory with δ > 1.
The paper is organized as follows. In the first two sections we discuss the perturbative approaches for calculating E k . In Sec. 2, the small-R behavior of E k is studied by means of the conformal perturbation theory for the bosonic Lagrangian (1.4). Then, in Sec. 3, using the fermionic Lagrangian (1.2), the vacuum energies are calculated within the second order of standard renormalized perturbation theory. The exact formula for the vacuum energies expressed through solutions of the classical sinh-Gordon equation is presented in Sec. 4. Our considerations there are essentially based on the previous works [13][14][15]. These connections allows one to derive a system non-linear integral equations which is well suited for perturbative analysis around δ = 0. Finally, Sec. 5 contains a summary of the original coordinate Bethe ansatz results [3] and the corresponding non-linear integral equations [12], as well as their numerical comparison with our calculation.

Small-R expansion
In this paper we shall mainly focus on the BL model with the vanishing exponent ν (1.3). Nevertheless it is useful to start with the theory characterized by a general set (g, ν). In the bosonic formulation, the model is still described by the Lagrangian (1.4), where the couplings (a 1 , a 2 ) substitute the pair (g, ν). These two pairs of renormalization group invariants are related as follows [3,6] Due to the periodicity of the potential term in ϕ i , the space of states splits into the orthogonal subspaces H k characterized by two "quasimomenta" k = (k 1 , k 2 ), As usual in the bosonization, the quasimomenta are related to the fermionic twists (1.5): The neutral (w.r.t. U(1) ⊗ U(1)) sector of the theory is described by the Bose fields with periodic boundary conditions: In the Euclidean version of (1.4), the periodic boundary corresponds to the geometry of infinite (or very long in the"time" direction x 0 ) flat cylinder Then the ratio E k /R would correspond to the specific (per unit length of the cylinder) free energy with the scalar operator exp i(k 1 ϕ 1 + k 2 ϕ 2 ) "flowing" along the cylinder. The UV conformal dimension of this operator is ∆ = 1 4 2 i=1 a i k 2 i . Therefore, we expect that at R → 0 The conformal perturbation theory for E k is constructed in the usual way [16] and yields an expansion in the dimensionless variable λ = 2πµ R 2π 1−ν , Here the first coefficient e (ν) 0 coincides with − c k 6 , while the subsequent ones are given by the perturbative integrals. In particular e (ν) 1 = I(p + ) + I(p − ), where p ± = 1 2 (a 1 k 1 ± a 2 k 2 ) and . (2.8) In the opposite large-R limit, the vacuum energy is composed of an extensive part which is proportional to the spatial size of the system and does not depends on the quasimomenta. The specific bulk energy, E ≡ lim R→∞ E k /R, has dimension [ mass ] 2 , i.e., E/M 2 is a certain function of the dimensionless couplings (g, ν). This universal ratio, along with another dimensionless combinations µ/M 1−ν , are fundamental characteristic of the theory, which allows one to glue together the small-and large-R asymptotic expansions. It is convenient to extract the extensive part from E k and introduce the scaling function Notice that it is a dimensionless function of the dimensionless variables r ≡ MR and k (and, of course, the couplings), satisfying the normalization condition lim r→+∞ F(r, k) = 0 . which can be resolved as We will assume that 0 < a 1 ≤ 1 ≤ a 2 , i.e., 0 ≤ δ < 1. A formal substitution of ν = 0 in (2.8) leads to a divergent expression. In order to regularize I(p), we cut a small disk |x| < ǫ in the integration domain D. As ǫ → 0, the regularized integral diverges logarithmically: where ψ stands for the logarithmic derivative of the Γ-function and γ E is the Euler constant.
In the case ν = 0, the general small-R expansion is substituted by the asymptotic series of the form e n (δ) (µR) 2n , (2.14) where explicitly and In ref. [6] Fateev presented strong arguments supporting the integrability of the BL model with ν = 0 and found an exact µ − M relation, (2.17) Using his results it is straightforward to obtain (see Sec. 4 bellow) the following expression for the bulk specific energy (2.18) One can see now that F, defined by eq. (2.9), does not contain any UV divergences, i.e., it is an universal scaling function of the dimensionless variable r = MR. Its small-R expansion can be written in the form where k ± = 1 2 (k 1 ± k 2 ), ρ = r 4π cos πδ (2.20) A few comments are in order here. As it was already mentioned in the introduction, the logarithmic divergence of E is well expected in the context of application of the BL model to the problem of instanton summation in the O(3) sigma model. The integration over the instanton moduli space leads to the divergent contribution of the small-size instantons [4]. So that ǫ can be interpreted as a cut-off parameter which allows one to exclude the divergent contribution of the small-instantons. Another comment concerns to the symbol ≍, which is used in eqs. (2.14) and (2.19) to emphasize the asymptotic nature of these power series expansions. To see that they have zero radius of convergence, it is sufficient to consider the case δ = 0. Returning to the fermionic description, the model (1.1) with g = 0 constitutes a pair of non-interaction Dirac fermions, so that there exists a closed analytic expression for the scaling function where πf/R 2 coincides with the specific free energy of the free Dirac fermion at the temperature 1/R and (imaginary) chemical potential 2πik/R, i.e., It is now straightforward to see that the power series (2.19) for δ = 0 is indeed an asymptotic expansion and where the superscript stands for derivative of (2n − 2)-order w.r.t. the argument. For nonzero δ the asymptotic coefficients e n (δ) with n ≥ 2 can be expressed in terms of the multiple integrals. Unfortunately such representation can not be used for any practical purposes. The only exclusion is e 2 (δ), whose integral representation can be simplified dramatically. For future references we describe here major steps in this calculation. First of all, using the complex coordinate z = exp(2π(x 0 + ix 1 )/R), the asymptotic coefficient e 2 (δ) can be represented as a 6-fold integral, Now let us substitute the integration variable , and integrate over z 1 by means of the identity 2 F 1 stands for the conventional hypergeometric function, and Finally, the integral over z 3 can be performed using a remarkable relation As a final result one obtains the following integral representation This formula allows one to achieve a reliable accuracy in the numerical calculation of e 2 (δ). For illustration, we present in Fig.1 the numerical results for k 1 = k 2 = 0. Notice that in this case the corresponding function τ 00 (ζ) in (2.28) can be expressed in terms of the complete elliptic integral of the first order K(ζ) = π 2 2 F 1 ( 1 2 , 1 2 , 1; ζ): Note that τ p 1 p 2 (ζ) given in (2.25) is a particular case of a more general function τ p 1 p 2 p 3 (ζ) defined by (4.13), namely τ p 1 p 2 (ζ) ≡ τ p 1 p 2 p 3 (ζ)| p 3 =0 . This function defines a real solution (4.16) of the Liouville equation (4.17), satisfying the asymptotic conditions (4.18).

Weak coupling expansion
We now consider a weak coupling expansion of the scaling function F. Since g = πδ 1−δ 2 = πδ + O(δ 3 ), no needs to distinguish g π and δ within the first two perturbative orders. It is convenient to define the perturbative coefficients through the relation: ). The results obtained in the previous section allows one to predict the leading small-R behavior of F i . Generally speaking the coefficients in the power series (2.19) admit the Taylor expansion e n (δ) = e n (0) + ∂en(0) In particular, as it follows from eq.(2.20), Also, using the original integral representation (2.23) for e 2 (δ), one can show that In the case p 1 = p 2 = 0 the weak coupling expansion includes only even powers of δ (see Fig. 1) and All of these can be used to study the short distance expansion of F i in eq.(3.1). 2 In particular, it is possible to show that Recall that the relations (2.16) between k i and p i involve the perturbative coupling. This should be taken into account since it is assumed that the expansion (3.1) is performed for fixed values of k i rather then p i . − + Figure 2: A diagrammatic representation of F 1 in eq. (3.1). The signs ± label the fermion "colors" ψ ± propagating along the loops (see Lagrangian (1.2)). and also where In the case k 1 = k 2 = 0, For finite values of r the perturbative coefficients F i can be calculated within the renormalized perturbation theory based on Lagrangian be the fermionic Matsubara propagator with the temperature 1/R and chemical potential 2πik σ /R. It can be expressed in terms of the modified Bessel function of the where γ a are Euclidean γ-matrices, {γ a , γ b } = 2 δ ab , and At the first perturbative order one has (see Fig. 2) Here Ψ σ andΨ σ stand for the components of the Dirac bispinors ψ σ with the Lorentz spin + 1 2 and − 1 2 , respectively. In zero-temperature limit the Lorentz invariance is restored and hence one observes that F 1 takes the form (3.5). It is also easy to see that where f = f(r, k) is given by (2.21). This is in a complete agreement with the short distance prediction (3.6). The second-order diagrams are depicted in Fig. 3. The type I diagram gives the contribution Because of the UV divergence at x = 0, the integration domain D ǫ here is chosen to be the cylinder (2.5) without an infinitesimal hole |x| < ǫ. One can show that, as ǫ tends to zero, In fact, since E k = R E + π R F, the quadratic divergence ∝ 1/ǫ 2 should be relocated to the specific bulk energy. Generally speaking, the specific bulk energy has a form where Λ ≫ M is some lattice energy scale and w is some (nonuniversal) function of the coupling g. Notice that, in writing eq. (2.18), the quadratic divergence was omitted (as usual in QFT). The type II diagrams from Fig. 3 leads to the UV finite integral over the whole cylinder D: Finally, the counterterm ∝ g 1 in (1.2) contributes through the type III diagrams, schematically visualized in Fig. 3. This can be written in the form Contrary to the one point functions (3.13), the condensate Ψ σ Ψ † σ (0) diverges logarithmically: where C is some constant. Since the UV divergence ∝ log 2 (ǫ) is canceled from the sum of types I and III diagrams if we choose g 1 = g 2 2π + O(g 3 ). As well as the quadratic divergence, the remaining logarithmic divergence should be relocated to the specific bulk energy. Expanding cos 2 πδ 2 in (3.18) one can find the value of the constant C: This way the second order correction takes the form where the finite constant should be adjusted to satisfy the normalization condition (2.10). It reads explicitly as Further calculations show that Here the shortcut notation c(k) = cos(πk) is used and K s (z) denotes the modified Bessel function of the second kind: Also in eq. (3.26) and bellow, the symbol o e −2r denotes a remaining term that decays faster than r −N e −2r for any positive N as r → +∞. Notice that the normalization condition F 2 = o e −r implies an absence of the finite renormalization of the fermion mass. It can be used for fixing the constant C in (3.21) and hence avoid any reference to the exact relation (3.18). For k 1 = k 2 = 0, the result of perturbative calculation is presented in Fig. 4. Using eqs. (3.5) and (3.14) it is easy to show that Thus, at least at the first-two perturbative orders, the leading large-r behavior of the scaling function F is defined by F 0 only and therefore . counted with the phase factor e iπ(σ 1 k 1 +σ 2 k 2 ) and, therefore, the summation over four possible sign combinations with σ 1,2 = ±1 gives rise eq.(3.29). Thus we may expect that the asymptotic formula (3.29) holds true as the mass of the first bound state M 1 = 2M cos( πδ 2 ) is greater than M, i.e., for δ ∈ [0, 2 3 ). Before concluding this section let us make a few remarks about the (non-integrable) case with a non-zero value of ν. Instead of adjusting the counterterm coupling g 1 , the logarithmic divergences can be absorbed by the mass counterterm with δM = ν log(M/Λ UV ), where . (This is an infinitesimal version of eq. (1.3) where M 0 = M + δM.) As it was mentioned in the introduction, the exponent ν and the four-fermion coupling g can be thought of as independent parameters for the family of BL models. Using eq. (3.20), it is easy to see that Finally we note that for ν = 0, (an universal part of) the specific bulk energy has a valid Laurent expansion of the form where h n (g) admit power series expansions in g 2 .

Exact formula for F(r, k)
The BL model with non-vanishing ν = 1 2 (a 1 + a 2 − 2) can be thought as a sort of analytical regularization of the model with ν = 0 -the integrals appearing in the conformal perturbation theory converge for negative values of ν, but become singular at ν → 0 − . A brief inspection of eq. (2.8) shows that a simple pole 1 ν replaces the logarithmic divergence 2 log(ǫµ) + const in (2.13) which occurs when the integral is regularized by excluding a neighborhood of the singular point from the integration domain. The BL with non-vanishing ν is a well defined QFT and it is interesting in itself in a context of applications in condensed matter physics [17]. However, as it was already mentioned in the introduction, the "ν-deformation" spoils the integrability. Remarkably that there exists an integrable deformation of the BL model with ν = 0. The corresponding model was introduced by Fateev in the works [6,18] and it will be referred to bellow as the Fateev model.
Contrary to the BL model, the Fateev (F) model involves three Bose fields governed by the LagrangianL Here α i = 1 2 √ a i and the coupling constants a i satisfy a single constraint which implies that the parameter µ has a dimension of mass. As α 3 → 0, the field ϕ 3 decouples in (4.1) and the interacting part coincides with the bosonic version of the BL Lagrangian (1.4) with a 1 + a 2 = 2. In fact, this observation requires a more careful assessment. Performing the limit α 3 → 0, one should expand the exponentials e ±i α 3 ϕ 3 in (4.1) to the terms ∝ a 3 = 4α 2 3 . The mass of the decoupled field is given by the relation m 2 = 8πµa 3 cos(α 1 ϕ 1 ) cos(α 2 ϕ 2 ) , where the vacuum expectation value is taken for the BL model with ν = 0. This expectation value is simply related to the corresponding specific bulk energy, cos(α 1 ϕ 1 ) cos(α 2 ϕ 2 ) = − 1 4 ∂E ∂µ , and hence Eq. (2.18) shows that E = 4πµ 2 log(µǫ) + . . . and, as has been argued above, should be replaced by E = 4πµ 2 a 1 +a 2 −2 + . . . within the analytical regularization. This, combined with (4.3) and the constraint a 3 = 2 − a 1 − a 2 , means that the field ϕ 3 has the mass m = 4πµ in the decoupling limit. Taking  where we use δ = 1 − a 1 = a 2 − 1.
One of Fateev's important results concerning the theory (4.1) is an elegant analytical expression for the specific bulk energy [6]: .

(4.5)
The linear constraint imposed on parameters a i , can be resolved by setting a 1 = 1 − δ − a 3 2 and a 2 = 1 + δ − a 3 2 , and, therefore, as a 3 → 0 one has Keeping in mind that 1 a 3 can be substituted by (− log(µǫ)) one find the relation where E is the specific bulk energy for the BL model (2.18), whereas the term ∝ m 2 is a contribution of the free massive field. Notice that const does not depend on the coupling δ, and it can be always set to zero. We can consider now the Fateev model in finite volume with the periodic boundary conditions ϕ i (x 0 , x 1 +R) = ϕ i (x 0 , x 1 ) imposed on all three fields ϕ i . Similar to the definition (2.9) for the BL model, let us introduce F F = R (E k − RE F )/π. Then the above consideration suggests that where the second term in the r.h.s. with corresponds to a contribution of the free boson of mass 2Mc δ 2 with c(x) ≡ cos(πx). Notice that the limit in (4.8) should be taken from negative values of a 3 , so that the Lagrangian (4.1) is real. For a i > 0 (i = 1, 2, 3) the Lagrangian is complex, but the QFT is still well defined. In this case the potential term in (4.1) is periodic w.r.t. all fields ϕ i and the space of states splits on the orthogonal subspaces characterized by a triple of quasimomenta k = (k 1 , k 2 , k 3 ). For a 3 < 0 different sectors of the theory are labeled by a pair of quasimomenta, similar to the case of the BL model, so that eq. (4.8) can be understood literally as a relation between the vacuum energies in the Fateev and BL models characterized by the same k = (k 1 , k 2 ).
A major advantage of the case with all positive a i is that the general structure of the small-R expansion in this regime is considerably simple compared to the case a 3 < 0. For a i > 0 the potential term in the Lagrangian (4.1) is a uniformly bounded perturbation for any finite value of the dimensionless product µR. Therefore the conformal perturbation theory yields an expansion of the form where p i = 1 2 a i k i . Notice that the integral diverges at a 3 → 0 + and formula (2.23) for the asymptotic coefficient e 2 (δ) in the BL model is a regularized version of e 2 + p 1 + p 2 + p 3 , 1 2 + p 1 + p 2 − p 3 , 1 + 2p 1 ; ζ . (4.14) The derivation follows the same steps outlined in Sec. 2; Fist of all, one should substitute the integration variables z 2 by ζ = (1−z 1 )(z 2 −z 3 ) (1−z 2 )(z 1 −z 3 ) . Then the integral over z 1 is performed using the identity (2.24) where p 1 is substituted by p 1 +p 3 . Finally one should use the identity generalizing (2.27): An important observation is that τ p 1 p 2 p 3 (ζ), considered as a function on the Riemann sphere, is regular except for three points ζ = 0, 1, ∞ and is a real solution of the Liouville equation for |p i | < 1 2 , i |p i | < 1 2 (for details, see e.g. ref. [19]). Notice that τ p 1 p 2 p 3 (ζ) = τ p 1 p 2 p 3 (1 − ζ) = |ζ| 2 τ p 3 p 2 p 1 (ζ −1 ) and therefore η L satisfy the following asymptotic conditions at the punctures: This way the result of conformal perturbation theory can be expressed in terms of solution of the Liouville equation on the three-punctured sphere S 2 /{0, 1, ∞}: where ρ = 1 2 µR and In ref. [13] it was conjectured that where η is a real solution of the so-called modified sinh-Gordon equation satisfying the the same asymptotic conditions as (4.18) (i.e., η L should be substituted by η in (4.18)). The last term in (4.22) ∝ ρ 4 and can be treated perturbatively for |p i | < a i 4 . Therefore the small-R behavior (4.19) follows immediately from the exact formula (4.21). One can show that the leading large-R asymptotic of (4.21) correctly reproduces the specific bulk energy (4.5) (see ref. [13] for details). Additional arguments in support of eq. (4.21) were presented in the work [14].
Eq. (4.21) can be transformed to a formula for the scaling function F F ≡ R (E F − RE F )/π. For this purpose, one should consider the Schwarz-Christoffel mapping which maps the upper half plane ℑm(ζ) ≥ 0 to the triangle (w 1 , w 2 , w 3 ) in the complex wplane (see Fig. 5). The lower half plane ℑm(ζ) ≤ 0 is mapped into the congruent triangle (w 1 , w 2 ,w 3 ). It is straightforward to show that the real functionη = η − 1 4 log(PP ) is a solution of the sinh-Gordon equation ∂ w ∂wη − e 2η + e −2η = 0 (4.24) in the open domain D (+) F , which is obtained by gluing together the triangles along their sides, as it shown in Fig. 5. At the singular points w = w i (i = 1, 2, 3) the solution has the following asymptotic behavior: In ref. [13] it was shown that formula (4.21) implies the relation Then, in the consequent paper [15], it was argued that (4.26), with some minor modifications, also applies to the case a 1,2 > 0, a 3 < 0. Namely, where nowη is a solution of the sinh-Gordon equation (4.24) in the domain shown in Fig. 6, satisfying the asymptotic conditions (4.25) at the vertices w 1 and w 2 , and η → 0 as |w| → ∞ .   Fig. 7. With the relation (4.8), this leads to the following exact formula for the scaling function F(r, k) in the BL model, As we shall see below this formula is in a perfect agreement with all perturbation theory calculations, considered in Sec. 2 and Sec. 3 of this paper, as well as with all other known results on the BL model, including the Bethe ansatz results of [3,12]. The sinh-Gordon equation (4.24) is a classical integrable equation which can be treated by the inverse scattering transform method. Thus the relation (4.29) allows one to apply this powerful method to the problem of determining the vacuum energies. The working is very similar to that for the Fateev model, considered in [14], where all a 1 , a 2 , a 3 > 0, though contains a few original details. We postpone these derivations to our future publication [40] but present the final result here. The scaling function (4.29) is expressed through the solution of a system of two Non-Linear Integral Equations (NLIE): Here σ = ±, (χ + , χ − ) = (0, πa 1 /2) and the kernels are given by the relations .
Once the numerical data for ε ± (θ) are available, F(r, k) (4.29) can be computed by means of the relation where Notice that (4.33) is valid for both choices of the sign ±. Eq.(4.33) can be compared against the predictions of renormalized perturbation theory in several ways. First, note that the integral equation (4.30) have a smooth limit for δ → 0 (its kernel vanishes linearly in δ). Using this property we have verified that the function F 2 in eq. (3.1), extracted from the numerical solution of (4.30)-(4.34) for k 1 = k 2 = 0 and 0.1 ≤ r ≤ 5, within nine significant digits coincides with the result of the perturbative calculations, shown with the solid line in the right panel of Fig. 4. Second, one can show that the exact formula (4.33) implies the following large-R asymptotics where k 1 = k + + k − , k 2 = k + − k − and c(x) ≡ cos(πx). Expanding this relation to the second order in δ = 1 − a 1 = a 2 − 1, one finds that the result is consistent with eqs. (3.1), (3.28) and (3.26) from Sec. 3. Third, the numerical values for F(r, k) obtained from (4.33) and presented in Fig. 8 and Tab. 1 on page 28, show an excellent agreement with the large-R asymptotic formula (4.35) and also with the predictions of the conformal perturbation theory, given by (2.19), (2.20) and (2.28). Finally, the exact expressions (4.29) and (4.33) perfectly agree with the Bethe ansatz results, considered in the next section.

Bethe ansatz results
As shown already in the original BL paper [3] the fermionic model (1.1) could be solved by the coordinate Bethe ansatz. In this section we review and extend their results. Within the Bethe ansatz approach the eigenvectors and eigenvalues of the Hamiltonian are parameterized through rapidities of pseudoparticles filling the bare vacuum state. These rapidities are determined by the Bethe Ansatz Equations (BAE). In the context of relativistic QFT models the number of pseudoparticles is infinite and, therefore, the related BAE require some regularization which makes that number finite. Following the BL paper [3] here we will impose a straightforward cutoff to the number of pseudoparticles. An alternative and in many respects more efficient lattice-type regularization is considered in our next paper [40].
Let N ≥ 2 be an even integer. The BAE of ref. [3] involve two sets of unknown rapidities Throughout this section we will assume that the indices ℓ and J always run over the above sets of values, respectively. With a slight change of notations and some minor corrections 3 the Bethe ansatz equations of ref. [3] (generalized for the twisted boundary conditions (1.5)) can be written as The parameter M is the bare mass parameter entering the coordinate Bethe ansatz calculation of [3] (denoted as "m" therein). Its relationship with the physical fermion mass M used in the previous sections follows from the requirement that the scaling function, determined by the BAE, at large distances should vanish as ∝ exp(−MR), i.e., exactly as the one in (3.29). As we shall see below this is achieved if one sets (see remarks after eq. (5.14)) This relation will be assumed in what follows. For practical purposes it is useful to rewrite BAE (5.2) in the logarithmic form As usual, the most difficult question in the analysis of BAE is to determine patterns of zeroes and the corresponding phase assignment in (5.4) for different states, in particular for the vacuum state. For the untwisted boundary conditions, p 1 = p 2 = 0, this question was studied in [3]. It was shown that for small values of |δ| ≪ 1 the vacuum roots {u ℓ } are real and their positions are given by an asymptotic formula whereas the roots {θ J } split into pairs The θ-roots become complex and form the so-called 2-strings with a more subtle phase assignment. Near the origin ℜe(θ J ) < 2/(π 2 δ) the phases are still consecutive, as stated in however for larger ℜe(θ J ) this is no longer true and the consecutive phase segments are divided by regions of "holes", where the RHS of (5.10b) jumps over several integers. A general description of this pattern is unknown. The arguments of [3] are based on the perturbation theory around the free fermion case with the untwisted boundary conditions (corresponding to δ = 0 and p 1 = p 2 = 0) and expected to work well for sufficiently small δ's and vanishing p's. We have verified this picture numerically. The arrangement of the vacuum roots for N = 16 and |δ| = 0.05 is illustrated in Fig. 9, where only a part of complex plane, containing a half of the roots is shown. For δ < 0 the formula (5.8) is valid for |ℓ| < 2/(π 2 δ). For larger values of ℓ the θ-roots form almost perfect 2-strings Our numerical analysis shows that essentially the same picture of zeroes 5 holds also for small non-zero values of p 1 and p 2 . In particular, the integer phases (5.9) and (5.10) remains the same, as they cannot change under continuous deformations of the boundary conditions.
Using BAE (5.4) one can show [40] that the vacuum energy (5.6) diverges quadratically for large N (cf. eq. (3.18)) where The phases of complex roots are not uniquely defined. Here we adopt the convention that the functions (5.5) entering (5.4) should not have jumps under small variation of roots near their exact positions. For that reason for δ < 0 we replace φ 2δ (θ) in (5.4) withφ 2δ (θ), wherẽ φ α (θ) = 1 2πi log sinh θ − 1 4 iπα sinh θ + 1 4 iπα , differs from (5.5) by the sign of the argument of the logarithm. As a result our 2-strings phases assignment in (5.10b) looks different, but nevertheless equivalent to the corresponding eq. (92) in [3]. 5 When p 1 , p 2 = 0, eqs. (5.7) and (5.8) should be modified, but (5.9), (5.10) and (5.11) remain intact. Then from the finite-size scaling arguments (applied in the context of the Bethe ansatz regularization of massive field theory models [20,21]) one expects that for N → ∞ the regularized expression for the energy where c k = 2 i=1 1−6a i k 2 i , reduces to the scaling function F(r, k) for the integrable case of the QFT model (1.1), (1.2). The constant C is non-universal, it is determined by the requirement F(r, k) → 0 as r → ∞. The relation (5.3) follows from the requirement that (5.14) has the same large distance decay exponent as in (3.29). For δ > 0 the formula (5.14) has been verified numerically. The values F(r, k) obtained from (5.6) with the solution of (5.4a), (5.4b) and (5.9) for N = 500 display a good agreement (to within at least three decimal places) with the more accurate results obtained from the NLIE (4.30), see Fig. 10.
Finally note that, as shown by Saleur [12], the BAE (5.4a), (5.4b), describing the δ > 0 vacuum state (5.9) filled by the real θ-roots, can be converted to a set of NLIE. After some minor corrections 6 these NLIE (generalized for the twisted boundary conditions (1.5)) can be written as where j = 1, 2, r 1 = 2r cos πδ 2 ,r 2 = r ,k 1 = k 2 ,k 2 = k + , k ± = 1 2 (k 1 ± k 2 ) . (5.16) The kernel G jl reads , (5.17) and . (5.18) Note that G 22 (θ) coincides with G ++ (θ) defined in (4.31). With these notations the scaling function (5.14) can be written as Note, that even though the equations (5.15) look totally different from (4.30) the resulting expression (5.19) for the scaling function is, in fact, exactly equivalent to (4.33). A complete proof of this equivalence is presented in our next paper [40]. It is also worth noting that from the point of view of numerical analysis the system (4.30) displays a much faster convergence than (5.15) and, therefore, requires lesser computational resources. Moreover, the system (4.30) is well suited for small δ analysis, whereas the eq. (5.15) becomes singular for δ → 0 (the latter fact has already been noted in [12], where the NLIE (5.15) for the untwisted case k ± = 0 were originally derived).

Conclusion
The Bukhvostov-Lipatov (BL) model [3] describes weakly interacting instantons and antiinstantons in the O(3) non-linear sigma model in two dimensions. In this paper we have studied various aspects of the BL model with twisted boundary conditions, using all wellestablished approaches to 2D massive integrable QFT, including the conformal perturbation theory (Sec. 2), the standard renormalized perturbation theory (Sec. 3) and the Bethe ansatz (Sec. 5). Moreover, in Sec. 4 we have proposed an exact formula (4.29) for the vacuum energy of the model, expressing it via a special solution of the sinh-Gordon equation (4.24) in the domain D BL (see Fig. 7). The required solutionη(w) decays at |w| → ∞ and obey the boundary conditions (4.25) at the singular points w 1 and w 2 . The connection to the classically integrable sinh-Gordon equation is rather powerful, since it allows one to obtain the non-linear integral equations (4.30), determining the vacuum energy in the form (4.33). We have shown that this formula perfectly matches all our perturbation theory calculations as well as the previously known coordinate Bethe ansatz results of Bukhvostov and Lipatov [3], and Saleur [12]. The comparisons were done both analytically (where possible) and numerically. Complete proofs and derivations of our exact results are postponed into the forthcoming publication [40]. The main idea of that work is to connect the functional equations for connection coefficients for the auxiliary linear problem for the sinh-Gordon equation (4.24) to the Bethe ansatz equations (5.2), arising from the coordinate Bethe ansatz [3]. This requires rather substantial works involving the particle-hole transformation and lattice-type regularization of the BAE, as well as some generalization of arguments of ref. [14], devoted to the Fateev model. Clearly, further study of the BL model is desirable. Indeed, almost all the considerations in this paper concerns the weak coupling regime 0 < δ < 1. However, the most interesting regime is the strong coupling regime δ > 1, where the BL model admits a dual description as the so-called sausage model [10]. Interestingly, this model turns into the O(3) NLSM, in the limit δ → ∞. This suggests that the instanton counting becomes exact in the strong coupling limit of the BL model. We intend to address this problem in the future.
The description of the vacuum state energy of the BL model in terms of the classical sinh-Gordon equation can be viewed as an instance of a remarkable, albeit unusual correspondence between integrable quantum field theories and integrable classical field theories in two dimensions, which cannot be expected from the standard quantum-classical correspondence principle. In the past two decades this topic has undergone various conceptual developments, which can be traced through the works [13,14,[22][23][24][25][26][27][28][29][30][31][32]. The commonly accepted mystery of this correspondence is slightly unveiled by our conformal perturbation theory calculations in Sec. 2 and Sec. 4. Indeed eqs. (2.28) and (4.19), expressing the vacuum energy in terms of the solutions (4.13), (4.16) of the Liouville equation (4.17) arise as a direct result of calculations, without any additional assumptions. It would be interesting to check whether these calculations can be generalized to other integrable QFTs where the correspondence to classical integrable equations is known. More generally, it would be very important to better understand connections of the above correspondence to mathematical structures arising in 4D gauge theories [33][34][35], calculations of amplitudes of high energy scattering [36][37][38] and dualities in finite dimensional quantum-mechanical systems [39].