Bukhvostov-Lipatov model and quantum-classical duality

The Bukhvostov-Lipatov model is an exactly soluble model of two interacting Dirac fermions in 1+1 dimensions. The model describes weakly interacting instantons and anti-instantons in the $O(3)$ non-linear sigma model. In our previous work [arXiv:1607.04839] we have proposed an exact formula for the vacuum energy of the Bukhvostov-Lipatov model in terms of special solutions of the classical sinh-Gordon equation, which can be viewed as an example of a remarkable duality between integrable quantum field theories and integrable classical field theories in two dimensions. Here we present a complete derivation of this duality based on the classical inverse scattering transform method, traditional Bethe ansatz techniques and analytic theory of ordinary differential equations. In particular, we show that the Bethe ansatz equations defining the vacuum state of the quantum theory also define connection coefficients of an auxiliary linear problem for the classical sinh-Gordon equation. Moreover, we also present details of the derivation of the non-linear integral equations determining the vacuum energy and other spectral characteristics of the model in the case when the vacuum state is filled by 2-string solutions of the Bethe ansatz equations.


Introduction
This paper is a sequel to our previous work [2] devoted to the study of the Bukhvostov-Lipatov (BL) model [1]. Originally, this model was introduced to describe weakly interacting instantons and anti-instantons in the O(3) non-linear sigma model in two dimensions, extending the results of [3,4]. It is an exactly soluble (1+1)-dimensional QFT, containing two interacting Dirac fermions Ψ a (a = ±), defined by the renormalized Lagrangian L = a=±Ψ a iγ µ ∂ µ − M Ψ a − g Ψ + γ µ Ψ + Ψ − γ µ Ψ − + counterterms . (1.1) This model provides a remarkable illustration to the idea of the exact instanton counting. Indeed, as explained in [2], the model can be reformulated as a bosonic QFT with two interacting bosons, which upon an analytic continuation into a strong coupling regime, becomes equivalent to the originating O(3) non-linear sigma model. Our interest to the BL model is also motivated by the study of an intriguing correspondence between Integrable Quantum Field Theories (IQFT) and Integrable Classical Field Theories in two dimensions, which cannot be expected from the standard correspondence principle. Over the past two decades this topic has been continuously developed, as can be seen from the works [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. The new correspondence yields an extremely efficient extension of the Quantum Inverse Scattering Method (QISM) [20] by connecting it with well-developed techniques of the spectral theory of ordinary differential equations (ODE) and, most importantly, the Classical Inverse Scattering Method [21]. The resulting approach is usually referred to as "ODE/IM" or "ODE/IQFT" correspondence. Ultimately, for a massive IQFT, it allows one to describe quantum stationary states (and the corresponding eigenvalues of quantum integrals of motion) in terms of special singular solutions of classical integrable partial differential equations. To this moment, the approach has already been applied to a number of IQFT models, such as the sine(h)-Gordon model [14], the Bullough-Dodd model [15] and various regimes [16,17,22] of the Fateev model [23], including the "sausage model" [24] (a one-parameter deformation of the O(3) sigma-model). Despite this progress, the subject, certainly, deserves further studies for a better understanding of the mathematical structure of the ODE/IQFT correspondence, as well as unraveling fundamental underlying reasons of its existence. In this work we continue to address these problems in the case of the BL model. As noted before, this model is integrable. Its coordinate Bethe ansatz solution was presented in the original BL paper [1]. The corresponding factorized scattering theory was proposed in [23]. The non-linear integral equations for calculating the ground state energy were derived in [25].
As in [2] we consider the theory (1.1) in a finite volume imposing twisted (quasiperiodic) boundary conditions on the fundamental fermion fields, Ψ ± (t, x + R) = −e 2πik ± Ψ ± (t, x) ,Ψ ± (t, x + R) = −e −2πik ±Ψ ± (t, x) . (1. 2) 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. The perturbative treatment of (1.1) leads to ultraviolet (UV) divergences. In [2] we used the renormalization procedure which preserves the integrability of the theory. In this case the fermion mass M could only have a finite renormalization and the only UV divergent quantity remaining in the theory is the bulk vacuum energy, defined as E = lim R→∞ E k /R. Therefore, it is convenient to extract the (divergent) extensive part from E k and introduce a scaling function which is simply related to the so-called effective central charge, F = −c eff /6. Notice that it is a dimensionless function of the dimensionless variables r ≡ MR and k, satisfying the normalization condition lim r→+∞ F(r, k) = 0. In [2] the scaling function (1.3) was computed in a variety of different ways. Let us briefly review this work and then describe new considerations contained in this paper. First, we used the renormalized Matsubara perturbation theory (recall that we consider the finite-volume theory) to calculate the first two non-trivial orders of the expansion of the scaling function (1.3) in powers of the coupling constant g. The physical fermion mass M is normalized by the large distance asymptotics of the scaling function. As for the coupling constant g, it is convenient to trade it off the for a physically observable parameter δ, which enters the S-matrix of the model. According to [23] the particle spectrum 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, with a 1 = 2 − a 2 = 1 − δ. Each of the factors S a coincides with the soliton S-matrix in the quantum sine-Gordon theory with the renormalized coupling constant a.
Further, in [2] we also used the conformal perturbation theory and computed the small-r expansion of the scaling function (1.3) to within the forth order terms O(r 4 ), inclusively. Then, applying rather non-trivial integral identities for the hypergeometric function, we discovered that the result of these perturbative calculations can be expressed through a particular solution of the Liouville equation. Guided by this peculiar observation and the previous results [16,17,22] concerning different regimes of the Fateev model we presented an exact formula for the scaling function F(r, k) = −f B 2r cos πδ where k ± = k 1 ± k 2 , in terms of a real valued solutionη of the classical sinh-Gordon equation ∂ w ∂wη − e 2η + e −2η = 0 (1.5) in the domain D BL . This domain is obtained from two inclined half-infinite strips of the complex PSfrag replacements πa 1 πa 2 r/4 w 1 w 2 w Figure 1: Domain D BL has the topology of a thrice-punctured 2D sphere. The horizontal width of D BL is controlled by a length of the segment (w 1 , w 2 ), which coincides with r/4. The boundaries of the strips in the upper and lower half planes are identified, as marked.
plane glued together along their infinite sides, as shown in Fig. 1. It has the topology of a twodimensional sphere with three punctures (w 1 , w 2 , ∞). At the singular points w = w i (i = 1, 2) on the real axis the solution has the following asymptotic behavior: andη → 0 as |w| → ∞ . (1.7) The conical angles at the points w 1 and w 2 (see Fig. 1) are determined by the parameters a 1 = 1−δ and a 2 = 1+δ, while the number k 1,2 are related to the twists parameters k ± = k 1 ±k 2 in (1.2). The quantity f B is the free energy of a free 1D boson, given by where β is the (dimensionless) inverse temperature. In this paper we present a derivation of the formula (1.4). It is based on the fact that the sinh-Gordon equation is a classical integrable equation, which can be treated by the inverse scattering transform method. This allows one to relate the classical conserved charges to the connection coefficients of ODE's arising in the auxiliary linear problem for the sinh-Gordon equation. It turns out that these connection coefficients are determined by the same Bethe ansatz equations that appear in the coordinate Bethe ansatz solution of the BL model. As a result eigenvalues and eigenstates of QFT integrals motion are determined by solutions of the classical sinh-Gordon equation. In addition, using standard approaches of integrable QFT, we derive another expression for the scaling function in terms of solutions of Non-Linear Integral Equations (NLIE) which are equivalent to Bethe ansatz equations, but more efficient for numerical analysis.
The organization of the paper is as follows. The coordinate Bethe ansatz for the BL model is considered in Sect. 2. The lattice-type regularization of the Bethe ansatz equations and associated particle-hole transformations are considered in Sect. 3. The functional relations for the connection coefficients in the auxiliary linear problem for the modified sinh-Gordon equation are derived in Sect. 4. It is shown that these functional relations are equivalent to the Bethe ansatz equation for the BL model. The calculation of vacuum eigenvalues of local and non-local integral of motion from the NLIE and their comparison against small-and large-R asymptotics derived in our previous paper [2] is given in Sect. 5. Appendix A contains details of the derivation of functional relations for the connection coefficients. Their properties in the short distance (CFT) limit, are discussed in Appendix B. Appendix C contains derivations of the NLIE, and Appendix D contains integral representations for connection coefficients. Appendix E contains details of small-and large-R expansions of the vacuum energy as well as eigenvalues of higher integrals of motion.

Bethe ansatz results
In this section we review and extend the original results of [1], concerning the coordinate Bethe ansatz for the BL model.

Coordinate Bethe ansatz
Consider (1+1)-dimensional quantum field theory of two interacting Dirac fermions Ψ a (a = ±), defined by the bare Lagrangian The γ-matrices are defined as γ 0 = σ 1 , γ 1 = i σ 2 , with σ 1 , σ 2 , σ 3 being the usual Pauli matrices. We will consider the theory in a finite volume, imposing twisted (quasiperiodic) boundary conditions on the fermionic fields, (2. 2) The standard perturbation theory leads to ultraviolet divergencies and, therefore, the Lagrangian (2.1) requires renormalization. Typically this procedure is performed over the physical vacuum state of the continuous theory. For instance, the renormalized Matsubara perturbation theory (with the coordinate space cutoff) for one-and two-loop diagrams was considered in our previous paper [2]. However, the traditional Bethe ansatz approach is formulated in terms of unrenormalized quantities and a bare (unphysical) vacuum state. Therefore, in this section we will work with the unrenormalized Lagrangian (2.1). The connection of the bare fermion mass M 0 , coupling constant g 0 and twist parameters p ± to their renormalized physical counterparts will be given below (see (2.15), (2.31), (2.29) (3.35)). We write the bispinors Ψ a (x) in the form where ψ + a (x) and ψ − a (x) stand for their components with the Lorentz spin + 1 2 and − 1 2 , respectively. The Hamiltonian corresponding to (2.1) reads Under the second quantization the quantities ψ s † a (x) and ψ s a (x) become creation and annihilation operators for a quasiparticle with the flavor a and spin s 1 2 (s = ±) at the point x. These operators satisfy the canonical anticommutation relations, (2.5) It is easy to check that the operators of the total number of quasiparticles of each flavor commute among themselves and with the Hamiltonian (2.4). Correspondingly, both these operators N ± are separately conserved quantities. In particular, this means that the stationary Schrödinger equation H |Φ = E |Φ has solutions with a fixed total number N = N + + N − of quasiparticles, Here |0 stands for the bare vacuum state, i.e., the state with no quasiparticles (which should not be confused with the physical vacuum state, obtained after filling the Dirac sea of negative energy states). The variables {a i } and {s i }, taking values a i = ± , s i = ±, label the flavors and the spins of created quasiparticles. Substituting (2.7) into the Schrödinger equation and using the commutation relations (2.5), one obtains a partial differential equation (2.8) where the "wave function" χ is understood as a vector in C 2N with the components given by the coefficients χ s 1 ,...,s N a 1 ,...,a N (x 1 , . . ., x N ). The quantities τ (i) k and σ (i) k (k = 1, 2, 3, i = 1, . . . , N ) denote the Pauli matrices, acting, respectively, on the flavor and spin indices of the i-th quasiparticle. Thus, the diagonalization of the Hamiltonian (2.4) is reduced to a quantum many-body problem (2.8) for N fermions with a δ-function potential. For N = 1 this leads to the free Dirac equations with the solution describing a single quasiparticle with the energy ε and momentum k, where The rapidity variable θ is a free parameter, which can take complex values (e.g., for the positive energy states one should substitute θ → iπ − θ).
For general values of N > 1, the problem (2.8) can be solved by the Bethe ansatz. First, note that if no coordinates {x i } pairwise coincide the interaction potential vanishes and the wave function becomes a linear combination of products of the free quasiparticle solutions (2.9). Using the antisymmetry of the wave function one can reduce the considerations to the case x 1 < x 2 < · · · < x N . Then we use the ansatz where the sum is taken over all permutations Q = {q 1 , q 2 , . . . , q N } of {1, 2, . . . , N }. The last two formulae define a solution of (2.8) with the total energy and momentum which is valid everywhere, except at the hyperplanes x i = x j , separating regions with different orderings of the coordinates. The continuity of the wave function at these boundaries imposes multiple linear relations between the coefficients A q 1 q 2 ...q N a 1 a 2 ...a N . Consider, for instance, the twoparticle case N = 2. Integrating (2.8) over x 2 from x 1 − ǫ to x 1 + ǫ, with ǫ → 0, and using (2.11), (2.12), one obtains exactly four linear relations where the quasiparticle scattering S-matrix is defined as (2.15) The key to the solvability of the problem (2.8) is that this matrix literally coincides with the R-matrix of the 6-vertex model, satisfying the Yang-Baxter equation [26] and the inversion relation where θ ij = θ i − θ j . Note also the "flavor conservation" property, and a special point where the S-matrix becomes proportional to the permutation matrix Eq.(2.14) can be readily generalized for general values of N . Let Q = {q 1 , . . . , q j , q j+1 , . . . , q N } and Q ′ = {q 1 , . . . , q j+1 , q j , . . . , q N } be any two permutations, differing from each other by the order of just two elements q j and q j+1 . Then, repeating the above arguments for the case x q j = x q j+1 , one obtains an overdetermined system of homogeneous linear equations (2.20) Thanks to the Yang-Baxter (2.16) and inversion (2.17) relations this system is consistent and has non-zero solutions for the coefficients A q 1 , ... q N a 1 ... a N . Namely, all of them can be linearly expressed through their subset corresponding to one particular permutation Q, e.g., through the set to which the system (2.20) imposes no restrictions. 1 At this point it is worth remembering that the numbers of quasiparticles of each flavor are separately conserved. This means that the vector (2.21) must be an eigenvector of the operator where, similarly to (2.8), the Pauli matrix τ

Bethe ansatz equations
To find the remaining undetermined coefficients (2.21) and the rapidities {θ i } one needs to use the periodicity of the wave functions implied by (2.2) and (2.7). Note that for this relation one assumes the following ordering x j < x 1 < · · · < x j−1 < x j+1 < · · · < x N of the coordinates on the circle. Then using (2.12) and (2.20) one obtains where the crossed symbol \ j means that the corresponding entry is omitted. With the help of (2.20) and (2.19) the last relation can be viewed as an eigenvalue problem for a set of operators which all have the same eigenvector (2.21). Here T (θ) stands for the transfer matrix of an inhomogeneous six-vertex model with a (horizontal) field, defined as As a consequence of the Yang-Baxter equation (2.16) the transfer matrices with different values of θ commute among themselves [ T (θ), T (θ ′ )] = 0. Moreover, due to (2.18) they also commute with the operator N + , defined by (2.22). Thus, all operators T (θ j ) in (2.25) and the operator N + can be simultaneously diagonalized.
We can now use celebrated results of Lieb [27] and Baxter [28,29] for the lattice 6-vertex model. The eigenvalues of T (θ) in the sector with fixed N + read 2 .
(2.27) This formula contains N + new (complex) parameters {u ℓ }, which are rapidities of auxiliary "magnons". The latter are quasiparticles of yet another kind, required for the diagonalization of the transfer matrix. The new parameters {u ℓ } are determined by the following Bethe Ansatz Equations (BAE) where the indices ℓ, ℓ ′ take N + different values, whereas the index J take N different values. The parameters p 1 , p 2 are defined as Different solutions of (2.28) correspond to different eigenvalues of the transfer matrix (2.26). Next, in order to satisfy the boundary conditions (2.23) the rapidities θ j should satisfy (2.25). Together with (2.27) these conditions take the form  [27][28][29][30], as well as in the original BL paper [1]. We do not reproduce these results here since explicit expressions for the eigenvectors are not used in this paper. It is worth noting again that the parameter M 0 used in this section is the bare mass parameter. Its relationship with the physical fermion mass M in (1.1) follows from the requirement that the scaling function (1.3), calculated from BAE, at large distances should decay as ∝ exp(−MR). As shown in [2] this is achieved if one sets In what follows this relation will always be assumed.

The vacuum state
It is useful to rewrite the above equations (2.28), (2.30) in the logarithmic form For the vacuum state the location of zeroes and the associated phases in (2.32) were analyzed in [1] and [2]. First note that the vacuum is flavor neutral. Thus, we set N + = N − and N = 2N + . We denote N + = N, assuming it is an even number. To make our notations identical to those of [2] we also need to shift the numbering of roots, assuming that the indices J and ℓ run over the values For the relativistic QFT (2.1) the number of quasiparticles filling the bare vacuum state should be infinite. Therefore, the number N plays the rôle of the ultraviolet cutoff. Standard estimates [1] show that for large N the roots behave as This indicates that the products in (2.28) and (2.30) require regularization for N → ∞. Moreover, the total energy (2.13) in this limit diverges quadratically and should be renormalized by subtracting its extensive part, as in (1.3). The most convenient way to do this is to use a "lattice-type" regularization, considered below.
3 Lattice-type regularization and particle-hole duality

Lattice-type regularization
Here we consider a lattice-type regularization of the above BAE (2.28), (2.30). It is achieved by replacing the relativistic phase term r cos(πδ/2) sinh(θ) in (2.30) by a suitable lattice-type expression 3) or, in the logarithmic form, Eqs.(2.28) and (2.32b) remain unchanged. The indices ℓ and J run over the same sets of integers as in (2.35). The resulting equations look like typical BAE for some 1+1 dimensional solvable lattice model, where the parameter N stands for a half of the number of sites for a periodic spin chain; the real valued parameter Θ controls the column inhomogeneity of Boltzmann weights, whereas p 1 , p 2 define twist parameters for quasiperiodic boundary conditions. The energy E N of the corresponding lattice state can be computed by the formula which is quite natural in the context of the light-cone lattice regularization of integrable QFT models (see [31] for the case of the sine-Gordon model). The original BAE (2.30), corresponding to the continuous model (1.1), are recovered in the limit when both N and Θ tend to infinity, while the scaling parameter is kept fixed. The lattice energy expression (3.5) then formally reduces to (2.34), up to a diverging additive constant. Note that our lattice-type regularization (3.1), (3.5) is different from that suggested in [25]. It is not difficult to check that the new BAE possess qualitatively the same patterns of vacuum roots as described in the previous subsection. First consider the case δ > 0. As before, repeating the arguments of [1,2], one concludes that all the roots θ J and u ℓ for the vacuum state are real and the integer phases have the same assignment (2.36). The BAE then imply that these roots form two ordered sets For large N and fixed Θ their distribution, is well approximated by the continuous density, which is determined by the standard lattice model methods [26], .
Similarly for the roots θ J one has For large N and Θ the roots split into two cluster centered at ±Θ. Their asymptotic positions there can be estimated with the formulae where the functions z θ (θ) and z u (u) are given by For large N and fixed Θ the energy (3.5) behaves as where the superscript "(1)" indicates that the vacuum is filled by the real roots {θ J }, which in this context are usually called "1-strings". The constant ε ∞ can be easily calculated by combining (3.5) and (3.11), In particular, for large positive Θ one has Similar considerations apply to the case δ < 0. For practical purposes it is convenient to always assume that the constant δ is positive, and consider different BAE, which are obtained from (2.28) and (3.2) by the negation of δ. In addition it is also convenient to simultaneously interchange p 1 and p 2 . Let {θ J } and {u ℓ } now solve the BAE modified in this way. As explained in [1,2], the θ-roots become complex and form 2-strings, which for large N are asymptotically approaching the values where ℓ runs over the set (2.35). The phase assignment of these complex roots was discussed in details in the first part of this work [2] (see Eqs. (5.10a), (5.10b) therein). The corresponding vacuum energy (3.5), E involves the sum over the 2-strings (3.17). The u-roots remain real. More precisely, as we shall see in the next section, their positions remain unaffected by the transformation (δ, p 1 , p 2 ) → (−δ, p 2 , p 1 ). Then, as follows from (3.17), the asymptotic density of the 2-strings exactly coincides with that of the u-roots given by (3.10), which implies that for N → ∞ Remarkably, the regularized vacuum energies for the above two cases exactly coincide for any finite values of N and Θ. Originally, we have noted this from numerical calculations and then proven analytically (see Appendix D for details). The working is essentially based on the particle-hole duality symmetry, considered below.

Particle-hole duality transformations
The lattice type BAE (2.28) and (3.2) can be brought to various equivalent forms by means of the so-called particle-hole transformations which are well known for sypersymmetric solvable lattice models [32,33]. We use this symmetry to bring the BAE to a form convenient for a derivation of NLIE valid for both positive and negative δ in the interval −1 < δ < 1. Consider the vacuum solution containing real roots (3.8) corresponding to δ > 0. Define two Laurent polynomials Q 1 (λ) and Q 3 (λ) of the degree 2N in the variable λ = e θ , whose zeroes are determined by the vacuum roots, where, as before, the indices ℓ and J run over the values (2.35) (note an iπ shift in the definition of θ (1) J zeroes). Note that Q 3 (λ) only depends on λ 2 , since N is assumed to be even and the products in (3.21) contain even numbers of factors. Eqs. (3.2), (2.28) can now be equivalently rewritten as ℓ }. Assume that these variables are fixed. Then Eq. (3.23a) can be viewed as an algebraic equation for zeroes of some Laurent polynomial P(λ) of the degree 4N. By construction, one half of its zeroes coincides with {λ (1) J }, but there are 2N additional zeroes, therefore, P(λ) factorizes as where Q 2 (λ) is some Laurent polynomial of the degree 2N, It is easy to show that this polynomial satisfies the following BAE . (3.28c) First, note that (3.28b) trivially follows from (3.25) and (3.26). 3 Next, equating the ratios P(λ ℓ q −1 ) obtained from alternative representations of P(λ) from (3.25) and (3.26), and then using the fact that Q 3 (λ) = Q 3 (−λ), one obtains . (3.28d) Together with (3.23b) this leads to (3.28a). Finally, combining (3.28d) and (3.28a) one obtains (3.28c). Thus, we have obtained a set of coupled BAE for the zeroes of the three Laurent polynomials Q 1 (λ), Q 2 (λ) and Q 3 (λ). It is useful to isolate the following three closed subsets of these equations: (I) the pair of original equations (3.23a) and (3.23b) involving Q 1 and Q 3 only; (II) the pair of equations (3.28b) and (3.28c) involving Q 2 and Q 3 only; (III) the set of three Eqs. (3.23a), (3.28a) and (3.28b) involving all three Q 1 , Q 2 , Q 3 .
The arrangement of the vacuum roots is illustrated in Fig. 2.
Note that under the "particle-hole" symmetry transformation the sets (I) and (II) transform into one another, whereas the set (III) remains invariant. Obviously the sets (I) and (II) describe the fillings of the vacuum state for δ > 0 and δ < 0, respectively. The above symmetry proves that the roots {u ℓ }, which determine zeroes of Q 3 (λ), remain unchanged under the transformation (3.29) (this fact has already been mentioned before Eq. (3.19)). Even though the lattice-type regularization is used here as a technical tool, it would be interesting to construct a lattice model, which actually leads to the BAE presented above. This could provide ideas about a lattice regularization for non-linear sigma models, which so far has not been properly understood. From the point of view of quantum groups such a lattice model could be related to the quantized affine superalgebra U q ( D(2, 1|α)) or U q ( sl(2|2)) (see [17] for a hidden quantum group structure of the general Fateev model). As noted in [25] it could also be related to U q (osp(2|2)) algebra R-matrices found in [34][35][36]. Figure 2: Roots of Bethe Ansatz equations with the lattice-type regularization, given in Sect. 3. The (blue) dots show zeroes of Q (1) (e θ ), the (red) asterisks show zeroes of Q (2) (e θ ) and the (green) crosses represent zeroes of Q (3) (e θ ).

Scaling limit
In the scaling limit (3.6) the number of BA roots in the central region |ℜe(θ)| < Θ in Fig. 2 becomes infinite. With the standard Weierstrass regularization the formulae (3.21) and (3.27), with |ℜe(θ)| < Θ, can be easily turned into convergent infinite products, defining analytic functions Q (λ) of the variable λ, having two essential singularities at λ = 0 and λ = ∞: (λ) (which form 2-strings) one needs to follow the reasonings of [37]. For |n| ≫ 1 one obtains where λ is, actually, a function of λ 2 ). For small values of r the above formulae give a very good approximation to the position of zeroes even for small values of |n| (see Appendix B).
The functions (3.30) satisfy the QFT version of BAE (3.23), (3.28). The later are obtained by rewriting the coefficients there in the form and making the substitution which is reverse to (3.1). Finally, following the arguments of [31,38], one can show that in the scaling limit (3.6) an appropriately scaled regularized vacuum energy is expressed in terms of the scaling function (1. A proof of this statement is presented in Appendix D. Note that another variant of the formula (3.34) for a simple momentum cutoff regularization (rather that the lattice-type regularization of this paper) was presented in [2], see Eq. (5.14) therein. Remarkably, Eq. (3.34) gives reasonably accurate results even for finite values of N and Θ, since the main subleading term in the LHS is of the order of O(e −2Θ ).

Connection to classical sinh-Gordon equation
The Bethe ansatz equations derived in the previous section allow one to make a connection of the BL model to the classical inverse scattering problem method for the modified sinh-Gordon equation (see [2,14,16,17]) for a complex-valued function η(z) defined on the Riemann sphere with punctures. Here ρ is an arbitrary constant, while P(z) is a function with three singular points located at z 1 , z 2 and z 3 , where Eq.(1.5), given in Introduction, is connected to (4.1) by a simple change of variables η(z,z) →η(w,w) = η(z,z) − 1 2 log ρ 2 P(z) , dw = ρ P(z) dz (4.4) and the domain D BL , shown in Fig. 1, is the image of the Riemann sphere with three punctures in the coordinates (w,w). However, in this section we prefer to work with Eq. (4.1).
For further references note also that Eq. (4.1) with the choice (4.3) is a particular a 3 → 0 case of a more general equation [16], where which is connected with the symmetric regime of the Fateev model [23].

Symmetries of the auxiliary linear problem
The equation (4.1) is the zero-curvature condition for an sl(2)-valued connection whereP(z) denotes the complex conjugate of P(z), and λ = e θ is the multiplicative spectral parameter.
Recall that the equations (2.28) and (2.30) are the main ingredients of the coordinate Bethe ansatz solution for the Bukhvostov-Lipatov QFT (1.1). Remarkably, as we shell see below, exactly the same equations are also satisfied by connection coefficients of the linear matrix differential equations associated with modified sinh-Gordon equation (4.1). The required calculations are only slightly different from those related to the Fateev model, considered in details in [17]. Therefore, here we only briefly sketch the main steps of working for the case of the vacuum state of the BL model. In this case the function e −η(z) is a smooth, single-valued complex function without zeroes on a Riemann sphere with three punctures at z = z 1 , z 2 , z 3 . The point z = ∞ is assumed to be a regular point on the sphere, where The asymptotic conditions at the punctures are determined by and Note that the above asymptotic conditions are equivalent to (1.6), (1.7) under the substitution (4.4). Below it will be convenient to use parameters p 1 , p 2 related to k ± in (1.2) as, 4 From now on we assume that (i, j, k) is a cyclic permutation of (1, 2, 3). Consider the auxiliary linear problem (4.7). Introduce three matrix solutions normalized by the following asymptotic conditions for z → z 3 . The constants β i are defined by (4.14) The above conditions uniquely determine the solutions provided The connection matrices are defined as and S Further analysis is based on symmetries of the linear differential equations (4.7). Let be a transformation involving a translation of the independent variable z along the contour γ i , going around the point z i anticlockwise (the variablez is translated along complex conjugate contourγ i ), accompanied by the substitution λ → λ q −1 i , where q 1 = e −iπδ , q 2 = e +iπδ , q 3 = 1, q 1 q 2 q 3 = 1 . (4.20) Using (4.2) and (4.6) it is easy to check that the substitutions (4.19) leave the system (4.7) unchanged. Therefore they act as linear transformations in the space of solutions. Namely, in the basis Ψ (i) they read where An important property of the linear system (4.7) is that a combined transformation Ω k • Ω j • Ω i , where (i, j, k) is a cyclic permutation of (1, 2, 3), is equivalent to the identity transformation The proof follows from (4.20) and the fact that γ k • γ j • γ i is a contractible contour which loops around a regular point z = ∞ on the Riemann sphere, see Fig. 3. Combining (4.21) and (4.23) with the definition (4.16) one easily obtains

Functional relations for connection coefficients
The non-trivial functional relation (4.24), complemented with the asymptotic (WKB) analysis of the differential equations (4.7), allows one to completely determine all connection matrices. To within some simple √ λ factors, the connection coefficients are meromorphic function of λ with two essential singularities at λ = 0 and λ = ∞. For the normalization (4.12), (4.13) they could have poles for real negative λ's solving the equations To take this into account it is convenient to introduce a function Note that √ λ Z(λ) is a meromorphic function of λ. It satisfies the functional relation and has the following asymptotics at large real θ where f B (β) is the free-boson free energy, which has already appeared in (1.8). The function Z(λ) was originally introduced in [39]. To parameterize the connection coefficients introduce a set of functions W σσ ′ (λ) of the variable λ, which are analytic everywhere, except the points λ = 0 and λ = ∞. They are defined as where σ, σ ′ , σ ′′ = ±1. Note that W σ ′ σ (λ) is, actually, a function of λ 2 . Specializing Eq.(4.24) for different choices of the indices (i, j, k) and using (4.17) one can obtains a number of important functional relations between W's (see Appendix A for details) σ ′ (λ/q)W (2) σ (λq)) , (4.30d) where q = e i πδ 2 and

Asymptotic expansions
The asymptotic behavior of the W-functions can be found from the WKB analysis of the linear differential equations (4.7). The required calculations are similar to those of the Fateev model (see Sect. 9 of Ref. [17]). For ℜe(θ) → ±∞ and ℑm(θ) < π 2 , one obtains where k 1 and k 2 are the same as in (4.10). The leading coefficient is given by where γ E is the Euler constant and ψ(x) = ∂ x log Γ(x). The constant terms (see Eq. (9.19) of Ref. [17]) are expressed through the regular parts η near the punctures (see (4.9a) and footnote 6 on page 25). Further, taking the limit a 3 → 0 in the results [16] and [17], devoted to the symmetric regime (4.3) of the Fateev model, 5 one can express the coefficient C 1 through the solutionη(w) of (1.5), satisfying the asymptotic conditions (1.6), (1.7), where the integral is taken over domain D BL shown in Fig.1. The asymptotic expansion of W is simply determined by (4.31) and the relation (4.30c).

Connection to the Bethe ansatz
Remarkably, it turns out that the zeroes of the functions W (i) (λ) (introduced above as connection coefficients for the differential operators (4.7)) satisfy exactly the same equations Next, recall that the coefficients W (i) (λ) depend on the parameters p 1 , p 2 defined in (4.10). To indicate this dependence explicitly, we will write these coefficients as W (i) (λ | p 1 , p 2 ). Similarly, the functions Q (BL) i (λ), introduced in Sect. 3.3, depend on the twist parameters p 1 and p 2 , so we will write them as Q (BL) i (λ | p 1 , p 2 ). Below we are going to establish that where σ, σ ′ = ±. Recall that here we assume that p 1 and p 2 are positive, whereas in Sect. 2 and Sect. 3 they were taking both signs. The correspondence is achieved by using the sign variables σ, σ ′ in RHS of (4.37). Moreover, in writing (4.37) we have assumed a particular choice of normalization factors and the coefficients α i , β i in (3.30), which so far were at our disposal since they do not affect the position of zeroes of Q (BL) (λ). With this identification it is easy to check the that the functional equations (4.30) imply all the BAE (3.23) and (3.28) in their scaling limit form (recall that the latter is obtained with the help of substitutions (3.32) and (3.33)). For instance, (4.30c) immediately leads the scaling limit of (3.23a) and (3.28b). Similarly, (4.30d) leads to the scaling limit of (3.28a). Finally, setting λ = q λ n and λ = q −1 λ (3) n in (4.30c) and combining the resulting relations one arrives to (3.23b). Eqs.(3.28c) and (3.28d) are derived in a similar way.
To complete the identification (4.37) one also needs to check that the zeroes of the connection coefficients W (i) (λ), determined by the differential equations (4.7), have precisely the same phase assignments (2.36) as the vacuum state zeroes of the BAE. In particular, that the phases are given by consecutive integers containing no "holes" in their distribution. Alternatively, it is enough to prove that W (i) (λ) and Q (i) (λ) have the same loci of zeroes. We have verified this statement in the small-r limit, when the differential equations (4.7) simplify considerably and become ODE's with explicity known algebraic potentials (see Appendix B). Using asymptotic (WKB) analysis of these ODE's we show that the corresponding distributions of zeroes precisely match (3.31) in the small-r limit. Moreover we have confirmed this coincidence by extensive numerical checks. On this basis we assume (4.37) to hold.

Non-linear integral equations
As is well known [40,41], the BAE can be transformed into NLIE which allow one to accurately calculate vacuum energies as well as eigenvalues of all higher integrals of motion. There are several ways to proceed starting with different sets of BAE discussed after Eq. (3.28). It turns out that the set (III) appears to be most convenient for our purposes. Indeed, in this case the resulting equations admit a regular expansion for small values δ, which is very convenient for the comparison with perturbation theory calculations of our previous paper [2]. The working is presented in Appendix C. The approach is somewhat similar to that of the regime a i > 0 of the Fateev model, considered in details in Ref. [17]. However, the presence of 2-strings among the vacuum roots makes the considerations rather tedious. As a final result, one obtains the following system of two NLIE: Here σ = ±, (χ + , χ − ) = (0, πa 1 /2) and the kernels are given by the relations .
Once the numerical data for ε ± (θ) are available, any (regularized) sum over the Bethe zeroes can be calculated via fastly converging integral representations (see Appendix D). All information about the vacuum eigenvalues of local and nonlocal integrals of motion is contained in the coefficients of the asymptotic expansions, where (ij) = (12) or (21), the parameter ρ is defined in (4.36) and the coefficients C  As demonstrated in [2] the exact formula (5.6) is in a perfect agreement with the results of renormalized perturbation theory and conformal perturbation theory. It is possible to show that (5.6) implies the following large-r asymptotics, where K α (x) is the modified Bessel function, k 1 = k + + k − , k 2 = k + − k − and c(x) ≡ cos(πx). Besides the vacuum energies the non-linear integral equations allows one to determine the vacuum eigenvalues of all higher Integral of Motions (IM). In particular, the vacuum eigenvalues of the local IM, I 2n−1 ,Ī 2n−1 for n > 1, are given by the formula generalizing (5.6), Some details concerning small-and large-R behavior of I 2n−1 , along with numerical data for I 3 and I 5 can be found in Appendix E. Among the nonlocal IM the special rôle is played by the so-called reflection operators [16]. Their vacuum eigenvalues S i (i = 1, 2) are related to the scaling function F(r, k) [16], and can be calculated through the solution of the system (5.1): On the other hand, taking into account (4.37) and comparing the constant terms in the expansions (4.31) and (5.4), one arrives to an alternative expression (4.35) of S i in terms of the regular part of the solution of the modified sinh-Gordon equation (4.1), (4.9). 6 The large-r behavior of S i follows immediately from Eqs.(5.8), (5.10): . (5.14) Using the small-r asymptotic expansion The first term in the RHS reads explicitly [16] S( .

(5.17)
In Fig. 4 the numerical results for 1 2 log(S i ) are compared against the large-and small-r asymptotic formulae (5.13) and (5.16).
Finally, there is an infinite set of non-local IM, whose vacuum eigenvalues H 2), are given by the relations: Some further details can be found in Appendix E. 6 The expression (4.35) can be rewritten in terms of the functionη(w,w), satisfying (1.5) and (1.6),

Summary
This work completes our study of the Bukhvostov-Lipatov (BL) model, started in [2]. Besides its applications to the instanton calculus in the O(3) non-linear sigma model, the BL model provides an ideal testground to developing new and refining existing methods for integrable QFT's. Our attention was mostly devoted to calculation of the scaling function (1.3) (which is simply related to the so-called effective central charge) of the BL model in a finite volume and the quasiperiodic boundary conditions. The scaling function was computed in a variety of different ways. For the readers' convenience let us briefly summarize our main results below: (i) the conformal perturbation theory for a few first terms in the small-r expansion of the scaling function; (ii) the renormalized Matsubara perturbation theory (recall that we consider the finite-volume theory) for two non-trivial orders of the expansion of the scaling function in powers of the coupling constant; (iii) an exact formula for the scaling function in terms of different sets of non-linear integral equations (NLIE) derived from the Bethe ansatz and functional relations; (iv) an exact formula for the scaling function in terms of a special solution of the classical sinh-Gordon equation, based on the ODE/IQFT correspondence and the classical inverse scattering transform; (v) renormalization of the coordinate Bethe ansatz results for the scaling function both with the lattice-type regularization and the simple momentum cut-off.
Remarkably, all the above approaches perfectly agree to each other. This allowed us to establish the ODE/IQFT correspondence between the quantum Bukhvostov-Lipatov model and the classical sinh-Gordon equation. Note that this correspondence has recently been extended to the strong coupling regime of the BL model with δ > 1 [42]. In this case it provides a dual description of the 2D sausage [24], which includes the O(3)-sigma model when δ → ∞. It would be interesting to further develop the "quantum-classical" duality, in particular, to explore its possible connections to the duality between quantum and classical systems with a finite number of degrees of freedom, studied in [43].

Appendix B. Conformal field theory limit
To understand properties of the connection coefficients it is instructive to first consider the small-distance limit ρ → 0. The linear equations (4.7) can be reformulated as second order ordinary differential equations. For example, making the substitution in the first equation of (4.7), one obtains Let ϕ where Wr[f, g] = f ∂ z g − g∂ z f denotes the usual Wronskian. Consider now the limit when ρ → 0, but the product µ = ρλ is kept fixed. In this limit the variablez completely decouples from (B.2) and the potential term there simplifies where z ij = z i − z j . The basis solutions, corresponding to (4.12) and (4.13), become where the constants η (i) reg should be computed at ρ = 0, while µ = ρλ is kept finite. Further, by a change of variables , the differential equation (B.2), (B.4) can be cast into a form more suitable for analytical and numerical analysis, The branch cuts of the power function 1 + e x δ−1 are chosen to be on the lines, where (1 + e x ) is real and negative. Therefore, there are no branch cuts in the x-plane, when ℜe(x) < 0. The symmetry transformations Ω 1 and Ω 3 defined in (4.19) above act as follows:

Zeroes of W
• the zeroes of W (1) σ ′ (µ) are located on the real negative axis of variable µ µ (1) n ≍ − ++ (µ) n µ n (numerical) µ n (asymptotic formula (B.14)) 1 0.58940086 i 0.   • the zeroes of W (2) σ (µ) combine into complex conjugate pairs (µ n , µ * n ) (n = 1, 2, . . . , ∞) in the left half-plane ℜe(µ) < 0, located just outside the wedge with an acute angle πδ centered around the negative real axis of µ (see Fig. 5) σ ′ σ are simple, real and negative, (B.14) The above asymptotic formulae are in a very good agreement with the numerical values of the zeroes obtained from a direct solution of the differential equation (B.5), (B.8) even for small values of n (see Tables 1,2,3). and possess the following symmetry under the complex conjugation. Note that in the strip −π < ℑm(θ) ≤ π the function f 1 (θ) has following poles and zeroes, Taking the logarithms, one obtains where p(θ, δ) and φ α (θ) are defined in ( ℓ } stand for the zeroes of the functions Q 1 (e θ ), Q 2 (e θ ) and Q 3 (e θ ), as defined by (3.21), (3.27). With these notations BAE (3.23a), (3.28a) and (3.28b) take the form Fig . 6 shows the contours C 1 , C 2 and C 3 enclosing zeroes of Q 1 (e θ ), Q 2 (e θ+iπ ) and Q 3 (e θ ), respectively. The horizontal segments of these contours go along the lines ℑm(θ) = ±ω i (mod 2π), i = 1, 2, 3, with the positive constants ω 1 , ω 2 and ω 3 , satisfying the inequalities where ω is the maximal deviation of the roots θ J from the line ℑm(θ) = π. The second inequality in (C.4) ensures that the contour C 2 in Fig. 6 encloses all zeroes of Q 2 (−e θ ). Next, the sum in (C.3a) and the first sum in (C.3b) can be easily written as contour integrals, for instance and similarly The handling of the second sum in (C.3b) is slightly more complicated, since the contour C 2 (which is enclosing the zeroes of Q 2 (e θ+iπ )), also surrounds (unwanted) poles of f 1 (θ) located on the line ℑm(θ) = πδ/2. Taking this into account one obtains Further considerations follow the standard way [40,41] of deriving the non-linear integral equations. Substitute (C.5) back into (C.3). Integrating by parts, stretching the contours towards both infinities along the real axis of θ and introducing the pseudo-energies one obtains , where i, j = 1, 2, 3,

The integral operators R (±)
ij (θ) are conveniently defined by their Fourier transforms. We use the following general convention The matrix elementsR ij (ν) then take the form For completeness we present also the Fourier transform of the source term in (C.8) for which one should use the principal value integral in (C.9). Multiplying both sides of (C.8) by the inverse of the matrix integral operator δ(θ) δ ij − R where the Fourier transform of the kernel G (±) (θ) is given bŷ .

Appendix D. Integral representations for sums over the Bethe roots
Once the non-linear integral equations (C.10) are solved, every sum over the vacuum Bethe roots can be effectively computed via exact integral representations, considered below. The technique works equally well both for finite values of N (when the number of Bethe roots is finite) as well as in the scaling limit (when N → ∞ and (C.10) is replaced by its field theory analog (5.1)).
Let us start with the case of finite N. Consider, for instance, the product formulae (3.21) and (3.27). With an account of (C.10) they can be converted into integral representations The Fourier transform of the kernel is defined by (C.9) (but with the principal value integral) whereF The derivation is standard and very much similar to that of the equation (C.10) itself. The only complication is that one needs to systematically take into account the fact that summing over the zeroes of Q 2 (e θ ) with integrals over the contour C 2 in Fig. 6 (similarly to (C.6)) contains extra contributions determined by the zeroes of Q 3 (e θ q −1 ). Moreover, the function ε 2 (θ) was excluded from the RHS of (D.1) in favor of ε 1 (θ) of by virtue of (C.12 ,F (±) ,F (±) ,F . Similar representations can be obtained for the two expressions for the energy (3.5), when the vacuum state is filled with real roots (1-strings) and complex roots (2-strings), respectively, J , δ) .
The undetermined coefficients α i , β i are fixed by the requirement that the functions (3.30) satisfy the functional relation (4.30c) (upon the identification (4.37)). The integral representations for Q (BL) i obtained in this way, lead to the asymptotic expansions (5.4), with the coefficients given by (5.9) and (5.18).
Finally, note that in the scaling limit (3.6) formula (D.4) leads precisely to (3.34). First, let us split the integration in (D.4) into two regions, where |θ ′ | < Θ and |θ ′ | > Θ. Then, it is not difficult to see that the contribution of the first integral to the LHS of (3.34) exactly equals to F(r, k) (in the form (5.6)), while the second integral gives the term "−F(0, k)" in the RHS of (3.34).

Appendix E. Small and large-R asymptotics of IM
Here we collect some facts concerning the vacuum eigenvalues of integral of motions in the BL model. The local IM are normalized in such a way that Here F n (x, y) are certain polynomials of two variables of degree n, such that where (a) n is the Pochhammer symbol and dots stand for monomials with degree lower than n. The polynomial F 2 and F 3 read explicitly [44]: 7 and The vacuum eigenvalues of local IM (5.9) can be written in the form , (E.4) 7 The polynomials I 2k−1 (P, Q) from Ref. [44] are related to F k (x, y) as I 2k−1 (P, Q) = 2 1−2k F k − P 2 n , Q 2 n + 2 n=δ−1 . where and, as r → ∞, and .
The relations (E.4) are valid for any choice of the sign ± and i = 1, 2. Also notice that for the non interacting case  This gives the corresponding large-r asymptotics (see Fig. 8). In the small-r limit one has lim R→0 R 2π An integral representation for these limiting values where discussed in Ref. [17]. Using the results of this work, it is possible to show that