Strong coupling expansion in $\mathbf{\mathcal N=2}$ superconformal theories and the Bessel kernel

We consider strong 't Hooft coupling expansion in special four-dimensional $\mathcal N=2$ superconformal models that are planar-equivalent to $\mathcal N=4$ super Yang-Mills theory. Various observables in these models that admit localization matrix model representation can be expressed at large $N$ in terms of a Fredholm determinant of a Bessel operator. The latter previously appeared in the study of level spacing distributions in matrix models and, more recently, in four-point correlation functions of infinitely heavy half-BPS operators in planar $\mathcal N=4$ SYM. We use this relation and a suitably generalized Szego-Akhiezer-Kac formula to derive the strong 't Hooft coupling expansion of the leading corrections to free energy, half-BPS circular Wilson loop, and certain correlators of chiral primaries operators in the $\mathcal N=2$ models. This substantially generalizes partial results in the literature and represents a challenge for dual string theory calculations in AdS/CFT context. We also demonstrate that the resulting strong-coupling expansions suffer from Borel singularities and require adding non-perturbative, exponentially suppressed corrections. As a byproduct of our analysis, we determine the non-perturbative correction to the above mentioned four-point correlator in planar $\mathcal N=4$ SYM.

Abstract: We consider strong 't Hooft coupling expansion in special four-dimensional N = 2 superconformal models that are planar-equivalent to N = 4 super Yang-Mills theory. Various observables in these models that admit localization matrix model representation can be expressed at large N in terms of a Fredholm determinant of a Bessel operator. The latter previously appeared in the study of level spacing distributions in matrix models and, more recently, in four-point correlation functions of infinitely heavy half-BPS operators in planar N = 4 SYM. We use this relation and a suitably generalized Szegő-Akhiezer-Kac formula to derive the strong 't Hooft coupling expansion of the leading corrections to free energy, half-BPS circular Wilson loop, and certain correlators of chiral primaries operators in the N = 2 models. This substantially generalizes partial results in the literature and represents a challenge for dual string theory calculations in AdS/CFT context. We also demonstrate that the resulting strong-coupling expansions suffer from Borel singularities and require adding non-perturbative, exponentially suppressed corrections. As a byproduct of our analysis, we determine the non-perturbative correction to the above mentioned four-point correlator in planar N = 4 SYM.

Introduction
In this paper, we address the problem of computing a special class of observables in fourdimensional N = 2 and N = 4 superconformal Yang-Mills theories (SYM) for an arbitrary 't Hooft coupling λ = g 2 YM N. A distinguished feature of these observables (denoted as F ℓ (g)) is that they can be expressed in terms of determinants of a certain semi-infinite matrix exp (F ℓ (g)) = det δ nm − K nm (g, ℓ) 1≤n,m<∞ , g = √ λ 4π . (1.1) The semi-infinite matrix K nm carries information about the dynamics of the underlying gauge theory. In addition to the coupling g it may depend on other ("kinematical") parameters such as a non-negative integer ℓ to be specified below. The derivation of the representation (1.1) is rather non-trivial and relies on different techniques -integrability [1] (in the case of four-point correlation function in planar N = 4 SYM) and localization [2,3] (in the case of leading non-planar correction to N = 2 free energy on S 4 ).
In N = 2 superconformal models that are planar-equivalent to N = 4 SYM (in particular, in some SU(N) models with matter in the fundamental, rank-two symmetric or antisymmetric representations) F ℓ (g) stands for the leading non-planar correction to free energy on S 4 [18][19][20][21]. The localization technique allows one to express it in terms of a matrix model integral whose evaluation (to leading non-planar order) leads to the determinant representation (1.1) (see Appendix A for details).
There is a priori no reason why the semi-infinite matrices in (1.1) corresponding to these different observables in the two different theories should be related to each other in a simple way. It is therefore surprising that in both cases the matrix K nm turns out to have essentially the same, universal form. Its matrix elements are given by integrals involving the product of two Bessel functions with indices of the same parity (1. 2) The dependence of K nm on the coupling constant g enters through a function χ, conventionally called the symbol of the matrix K. As we will explain below, the value of a non-negative integer number ℓ and the explicit form of the function χ(x) depend on a choice of a particular observable.
It is remarkable that the same matrix (1.2) has previously appeared in a completely different area of mathematical physics and mathematical analysis. Namely, for ℓ = 0 and for the special choice of the function χ(x) = θ(1 − x), the determinant (1.1) describes the level spacing distributions in the Laguerre ensemble in random matrix theory [22,23]. More precisely, the function exp(F ℓ=0 (g)) gives the probability that there are no eigenvalues on the interval [0, g 2 ]. This function has a number of interesting properties and can be found exactly in terms of a Painlevé transcendent [24].
It is convenient to think about the semi-infinite matrix K nm as representing a certain integral operator B ℓ on a space spanned by basis functions, B ℓ ψ m (t) = K nm ψ n (t) (see Appendix B for details). In mathematical literature, this operator is known as a truncated (or finite temperature) Bessel operator [25]. The function (1.1) is then a Fredholm determinant of this operator, det(1 − B ℓ ). Such an identification proves to be very useful because it allows us to study the observable (1.1) using powerful methods developed in the literature (for reviews see [26][27][28][29]).

Applications
In the application to superconformal gauge theories, the relations (1.1) and (1.2) provide a concise representation of the observable F ℓ (g) that is valid for an arbitrary 't Hooft coupling. In this paper, we exploit this representation to evaluate F ℓ (g) at both weak and strong coupling and, then, study a transition between the two regimes. A detailed discussion of F ℓ (g) in (1.1) in N = 2 superconformal theories is presented below in Section 4. Each observable (1.1) corresponds to a specific choice of the non-negative integer ℓ and the function χ(x). We will encounter three different examples of the choices of ℓ and χ(x).
The first example corresponds to χ BES (x) = 2 1 − e x = 1 − coth (x/2) . (1.3) For ℓ = 0 the resulting semi-infinite matrix (1.2) governs the cusp anomalous dimension in planar N = 4 SYM through the BES equation [30]. 1 The second example comes from the study of four-point correlation functions of half-BPS operators in planar N = 4 SYM. In the limit, when the R-charge of the operators becomes infinitely large, this correlation function factorizes into a product of two building blocks, the octagons [4,5]. They depend on the kinematical variables y and ξ which are expressed in terms of the two cross-ratios built out of the coordinates of the four operators 2 as well nonnegative integer ℓ, the so-called bridge length. This parameter defines the smallest number of scalar propagators stretched between the four operators at zero coupling. Its value depends on the choice of polarizations of the operators. For the so-called "simplest" correlation function, the scalar propagators connect the four operators in a sequential (and not pairwise) manner. In this case, the octagon is given by (1.1) and (1.2) with ℓ = 0 and χ oct (x) = cosh y + cosh ξ cosh y + cosh x 2 + ξ 2 . (1.4) The properties of the octagon for arbitrary y and ξ have been studied in Refs. [7,[9][10][11][12]. For our purposes it will be sufficient to consider the special kinematical configuration y = ξ = 0, in which case χ oct (x) = 2 1 + cosh x . (1.5) The choice of y = ξ = 0 corresponds to the so-called bulk point singularity of four-point correlation function [31]. One reason why this kinematical point is interesting is that, as we show in Appendix A, the corresponding expression for the octagon admits a matrix model representation similar to that coming from the localization. The third example is related to a class of special N = 2 superconformal gauge theories that are planar-equivalent to N = 4 SYM. In these theories, the localization technique can be used to compute some observables (free energy, circular Wilson loop, two-point correlation functions of chiral operators) as functions of λ and N in terms of a non-gaussian matrix model (see [2,32]). It turns out that the leading 1/N 2 non-planar corrections in the corresponding matrix integrals can be expressed in terms of the determinant (1.1) and (1.2) with ℓ = 1, 2 3 and . (1.6) Below we first study the observable (1.1) for a generic function χ(x) and, then, specify the resulting expressions to the cases of the symbols χ(x) in (1.3), (1.5) and (1.6) that are relevant for physics applications.

Weak and strong-coupling expansions
Computing F ℓ (g) in (1.1) and (1.2) for an arbitrary value of 't Hooft coupling is a challenging problem. A significant simplification occurs in the limits of weak and strong coupling g.
At weak coupling, it is straightforward to expand the determinant (1.1) in powers of g 2 . Changing the integration variable in (1.2) as x → xg, we find that the matrix elements scale as K nm = O(g n+m+ℓ+1 ). As a consequence, to any finite order in g 2 , an infinite-dimensional matrix K nm can be replaced by its finite-dimensional minor and the determinant (1.1) can be expanded in powers of tr(K r ) (with r = 1, 2, . . . ).
At strong coupling, the study of asymptotic behaviour of the determinants like (1.1) has a long history in mathematical analysis, see, e.g., [27][28][29]. Relying on the strong Szegő limit theorems [33,34], we expect that, for sufficiently smooth function χ(x), the determinant (1.1) should admit a semiclassical expansion in the effectiveh ∼ 1/g. This leads to the following asymptotic behaviour where ∆F ℓ (g) vanishes for g → ∞. Each term on the right-hand side of (1.7) depends on the symbol χ(x) and has a different origin. The first three terms in (1.7) give the expression for F ℓ (g) which is known as the Szegő-Akhiezer-Kac (SAK) formula [33][34][35][36]. It involves the coefficients A 0 , A 1 and B, where B is called the Widom-Dyson constant. For the matrix (1.2), these coefficients are known in mathematical literature only for the unphysical values ℓ = ± 1 2 , see [37]. For arbitrary ℓ, the expressions for these coefficients were conjectured in [12].
The appearance of O(log g) term in (1.7) can be attributed to Fisher-Hartwig singularity of the symbol χ(x), see [38]. It corresponds to the behaviour where the parameter β defines the strength of the singularity. The coefficient A 1 in (1.7) depends on β and vanishes for β = 0. For the symbols (1.3), (1.5) and (1.6) the parameter β is different from zero and takes (half) integer values. Notice that 1 − χ(x) is an even (odd) function of x for even (odd) 2β. The last term in (1.7) can be split into the sum of two terms of different kinds ∆F ℓ (g) = f ℓ (g) + ∆f ℓ (g) . (1.9) Here f ℓ (g) is given by a perturbative series in 1/g A k+1 2k(k + 1) g −k , (1.10) where the expansion coefficients can be expressed in terms of the symbol χ(x) [11,12]. The function ∆f ℓ (g) in (1.9) describes non-perturbative exponentially small (inh −1 ∼ g) corrections to F ℓ (g). We show below that, for the symbols (1.3)-(1.6), the properties of the perturbative series (1.10) depend on the value of the parameter β, i.e. the strength of the Fisher-Hartwig singularity in (1.8) or, equivalently, on the parity of the function 1 − χ(x). For half-integer β, the series (1.10) can be resummed to all orders in 1/g to yield a well-defined function of g. 4 For integer β, the expansion coefficients in (1.10) grow factorially. Performing the Borel transform 11) we find that the function B f (σ/g) has poles at positive σ. As a consequence, the function f ℓ (g) is ill-defined and requires a regularization of the Borel singularities.

Non-perturbative corrections
The non-perturbative function ∆f ℓ (g) in (1.9) takes the form ∆f ℓ (g) = c 1 g n 1 e −8πg x 1 (1 + O(1/g)) , (1.12) where the O(1/g) term denotes a series in 1/g. We show below that the parameters c 1 , n 1 and x 1 depend on the choice of the function χ(x) and have a simple interpretation, e.g., x = 2πx 1 is the root of 1 − χ(x) of degree n 1 closest to the origin. In general, the expression on the right-hand side of (1.12) contains an infinite sum of terms of the same form but with different parameters c i , n i and x i . The leading contribution to (1.12) comes from the term with the minimal value of x i . Despite the fact that the non-perturbative corrections (1.12) are exponentially small at large g, they play an important role in understanding the properties of the function (1.7) in the transition region from strong to weak coupling. Determining the non-perturbative function ∆f ℓ (g) is an important open problem that we address in this paper.
For the symbol χ(x) given by (1.3), a closely related question of finding non-perturbative corrections to the cusp anomalous dimension in planar N = 4 SYM was studied in [39,40].
In what follows we employ the approach developed in these papers to determine ∆f ℓ (g) for a generic symbol χ(x) including the cases of (1.5) and (1.6).
We show below that, similarly to the perturbative series (1.10), the non-perturbative function ∆f ℓ (g) has different properties for integer and half-integer β in (1.8). For halfinteger β, the coefficients in (1.12) can be determined unambiguously in terms of the function χ(x). For integer β, due to the presence of Borel singularities in the perturbative series (1.10), the functions f ℓ (g) and ∆f ℓ (g) are not well-defined separately. However, all ambiguities related to a freedom in regularizing these singularities cancel in their sum (1.9). We determine the leading non-perturbative correction (1.12) by specifying a regularization procedure of the perturbative series (1.10). It amounts to a deformation of the integration contour in the Borel transform (1.11) in the vicinity of poles of the integrand.
As a relevant example, illustrating different properties of (1.1) for integer and halfinteger β, one can consider an asymptotic expansion of a (properly normalized) modified Bessel function (1.13) At large λ, the second term is exponentially small. For half-integer β, P (x) is a polynomial in x of degree 2β + 1. For integer β, P (x) is given by the product of √ x and a Borel non-summable series in x. We show below that F ℓ (g) has similar properties for the symbols (1.3)-(1.6). For β = 1, the relation (1.13) yields the strong-coupling expansion of the circular Wilson-Maldacena loop in planar N = 4 SYM [41,42] (see (4.13) below). In this case, the first and the second terms on the right-hand side of (1.13) define, respectively, the perturbative and the non-perturbative contribution to the circular Wilson loop at strong coupling.
The rest of this paper is organized as follows. In Section 2 we present the relation of the observable (1.1) to the Fredholm determinant of the truncated Bessel operator and discuss its expansion at weak and strong coupling. We show that the strong-coupling expansion of (1.1) has different properties for odd and even functions 1 − χ(x). In the latter case, the perturbative series in 1/g is not Borel summable and requires a regularization. In Section 3, we analyse the resulting non-perturbative, exponentially small corrections that appear in (1.1). We apply the approach developed in [40] to establish the relation (1.12) and present the explicit expressions for the coefficients there. In Section 4 we use the general results of Sections 2 and 3 to compute the strong-coupling expansion of the leading non-planar corrections to observables in special N = 2 superconformal SU(N) models. We find, in particular, the analytic expressions for the leading strong-coupling coefficients that were previously estimated only numerically (cf. [19][20][21]). We summarize our results and make some concluding comments (in particular, on possible dual string theory interpretation) in Section 5. Some technical details are presented in Appendices A-F.

Szegő-Akhiezer-Kac formula
In this Section, we summarize the properties of the function F ℓ (g) defined in (1.1) and (1.2) and present the expressions for the coefficients in its strong-coupling expansion (1.7).

Truncated Bessel operator
Discussing the properties of (1.1), it is convenient to switch from the semi-infinite matrix K nm to an integral operator B ℓ (χ) defined as where f (t) is a test function and the kernel is given by In mathematical literature, it is called truncated Bessel operator, see, e.g., [25]. The relation between the semi-infinite matrix K nm and the Bessel operator B ℓ (χ) is discussed in Appendix B. The function (1.1) admits a representation as a Fredholm determinant of the Bessel operator 3) follows from the following identity for the trace of the product of n copies of the matrices (1.2) (see [12]) The representation (2.3) is very useful for deriving an expansion of the function F ℓ (g) at small g. As mentioned in the Introduction, it can be obtained by expanding the determinant (2.3) in terms of traces of powers of the operator B ℓ According to its definition (2.4), tr(B n ℓ ) is given by the n-fold integral (2.4). Changing the integration variables in (2.4) as t i → gt i and taking into account (2.2), we find that tr(B n ℓ ) = O(g 2n(ℓ+1) ) and, therefore, the expansion in (2.5) runs in powers of g 2(ℓ+1) . The leading term of the expansion looks as where we introduced The subleading terms in (2.6) are given by multi-linear combinations of the coefficients q k (their expressions can be found in [11,12] where ζ(n) is the Riemann zeta-function. Notice that the function F ℓ (g) receives corrections only starting at order O(g 2(ℓ+1) ). This property does not depend on the explicit form of the function χ(x). In particular, it holds in both cases (1.5) and (1.6) mentioned above but its physical interpretation is different. For the octagon, the leading correction to F ℓ is associated with a scattering of elementary excitations (magnons) off a heavy state built of ℓ scalar particles in N = 4 SYM. At weak coupling, the corresponding amplitude behaves as g 2(ℓ+1) . In N = 2 superconformal models planar-equivalent to N = 4 SYM, the partition function does not depend on the matter content of the theory to first few orders in g 2 and, therefore, it does not get corrections (coinciding with the partition function in N = 4 SYM which is protected).

Specifying the symbol
The determinants (1.1) and (2.3) depend on the function χ(x) in a non-trivial way. In mathematical literature, this function is called a symbol of the Bessel operator (2.1). It is assumed to be a smooth function on real semi-axis.
For the strong-coupling expansion (1.7) to be well-defined, the symbol χ(x) has to verify additional conditions. The expansion coefficients in (1.7) and (1.10) are given by multilinear combinations of the integrals (see (2.25) For these integrals to be finite, the function 1 − χ(x) should be positive definite for x > 0 and decrease at infinity faster than 1/x. It is easy to see that the symbols (1.3), (1.5) and (1.6) satisfy these conditions. Notice that these symbols are expressed in terms of the same function and, as a consequence, they are related to each other as Substituting these expressions into (2.9) we find in particular, Generalizing the relations (1.3), (1.5) and (1.6), let us choose the following ansatz for the function χ(x) where the parameters b, x n and y n are real positive numbers. For this choice of χ(x), the integrals (2.9) become functions of x n and y n , e.g. (2.14) We also assume that χ(x) is analytical at the origin, so that 2β takes integer values. For integer (half-integer) β, the symbol (2.13) is an even (odd) function of x.
As was mentioned in the Introduction, the parameter β defined in (1.8) plays an important role in our analysis. According to (2.13), it controls behaviour of the symbol around the origin Substituting this relation into (1.2) and requiring the matrix elements K nm to be finite for n, m ≥ 0, we find that ℓ and β have to satisfy (2.16) The relation (2.15) implies that the symbol (2.13) possesses the Fisher-Hartwig singularity [38]. It is responsible for the appearance of the O(log g) term in the exponent of (1.7). The corresponding coefficient is given by [12] It only depends on ℓ and β and is insensitive to the values of x n and y n in (2.13). As follows from (2.13), the function 1 − χ(x) has an infinite number of poles and zeros in the complex x-plane. They are located along the imaginary axis at x = ±2πiy n and x = ±2πix n , respectively. It proves convenient to decompose the function 1−χ(x) into a product of functions analytical in upper and lower half-planes (the Wiener-Hopf decomposition) where x n and y n are positive. The function Φ(x) has poles and zeros located in the lower half-plane. By definition, it satisfies the normalization condition Φ(0) = 1. Some of the parameters x n and y n may take infinite value so that the number of poles and zeros can be different. In addition, some of x's and y's can coincide so that the roots and poles can be double, triple, etc. Recall that the symbol χ(x) has to vanish at infinity. Then, it follows from (2.18) that for x → ∞ This relation imposes non-trivial conditions on large n behaviour of x n and y n in (2.18). Let us examine the Wiener-Hopf decomposition (2.18) of the symbols (2.10). For the (2.20) Matching this relation to (2.18), we identify the values of roots x n = n − 1/2 and poles y n = n (with n ≥ 1). For the two remaining symbols in (2.10) we get It is straightforward to check that the functions Φ(x) in (2.20) and (2.21) verify the relation (2.19). Notice that the poles and roots of Φ oct (x) and Φ loc (x) are double degenerate. We will show below that this has important consequences for the properties of non-perturbative corrections at strong coupling. Recall that ℓ and β have to satisfy the condition (2.16). For the symbols (2.10) this leads to ℓ BES ≥ 0, ℓ oct ≥ 0 and ℓ loc ≥ 1.

Widom-Dyson constant
For the symbol (2.13) the Widom-Dyson constant is given by [12] , (2.22) where the subscript on the left-hand side was introduced to indicate its dependence on ℓ.
Here G(x) is the Barnes function satisfying G(x + 1) = G(x) Γ(x) and ψ(k) is given by a Fourier transform of log(1 − χ(x)) The relation (2.15) translates to ψ(k) ∼ −β/k at large k. The second term inside the brackets in the first term in (2.22) ensures that the integral converges at large k.

Perturbative corrections at strong coupling
With A 1 given by (2.17) the expressions for the remaining expansion coefficients in (1.7) and (1.10) are [12] A 0 = 2I 0 , where ℓ β = ℓ + β and I n are defined in (2.9) and (2.11). Replacing the coefficients in (1.10) with their explicit expressions (2.25), we obtain the function f ℓ (g) in (1.9) which describes the subleading corrections to (1.7) suppressed by powers of 1/g. This function has the following interesting properties.
According to (2.25), the expansion coefficients involve powers of I 1 in (2.9). All such terms can be eliminated at once by shifting the coupling constant as thus getting (2.27) Because the coefficients I n are independent of ℓ, it is obvious from this relation that f ℓ (g) depends on ℓ only through ℓ β = ℓ + β. We show below that the same is true for the nonperturbative function ∆f ℓ (g) in (1.7). We can use this observation to show that the function f ℓ (g) ≡ f ℓ β (g) has different properties for integer and half-integer β.
For half-integer β, or, equivalently, half-integer ℓ β = ℓ + β, the series (2.27) simplifies dramatically. For instance, for ℓ β = 1 2 and ℓ β = 3 2 all but the first few terms of the expansion (2.27) vanish. For an arbitrary half-integer ℓ β = n + 1 2 , or equivalently β = n + 1 2 − ℓ, the series (2.27) can be resummed to all orders in 1/g, see [11,12]. The resulting expression for f ℓ (g) looks as f ℓ (g) = n(n + 1) 2 log(g ′ /g) + log P n(n+1) where n = β + ℓ − 1 2 is non-negative integer and g ′ = g − 1 2 I 1 . Here P n(n+1) with the expansion coefficients given by multilinear combination of I k with k ≥ 2. For instance, for n = 0, 1, 2, 3 we have For integer β, or equivalently integer ℓ β , the situation is different. For a generic symbol χ(x) the expansion coefficients (2.25) are different from zero and the function f ℓ (g) is given by an asymptotic series in 1/g with factorially growing coefficients. Moreover, as we will show below, this series is not Borel summable and its Borel transform (1.11) develops a pole at positive σ where x 1 is the smallest root of the symbol (2.13) of degree n 1 . As a consequence, the series (2.27) approximates the function f ℓ (g) up to an exponentially small correction proportional to the residue of (1.11) at the pole σ = 8πgx 1 . The latter takes the same form as the non-perturbative correction (1.12). To define f ℓ (g) unambiguously, one has to specify the prescription for deforming the integration contour in (1.11) in the vicinity of the Borel pole.
We return to this question in Section 3.4 below.

Physics applications
In this subsection, we combine together the above relations and present the results for the strong-coupling expansion (1. where dots stand only for non-perturbative corrections of the form (1.12). Indeed, according to (2.20), the parameter β is half-integer in this case and, therefore, the perturbative 1/g n corrections to F BES ℓ (g) are expected to be very simple. Namely, the first two relations in (2.31) do not receive perturbative corrections in 1/g at all. In the last relation, they all come only from the expansion of log(1 − log 2 2πg ). These properties are in a perfect agreement with the relations (2.28) and (2.29) for n = ℓ − 1 and ℓ = 0, 1, 2.
Yet another remarkable feature of the symbol (1.3) is that the functions (2.31) can be found exactly for arbitrary g, see [9,43,44]  It is easy to check that, at strong coupling, these relations reproduce (2.31). In addition, they allow us to identify non-perturbative corrections to (2.31). We discuss them in Section 3.5 below.
We can use (2.32) to define the following functions As was shown in [9,45], these "octagon anomalous dimensions" determine the asymptotic behaviour of the four-point correlation function of infinitely heavy half-BPS operators in the limit when the four operators are located in the vertices of a light-like rectangle. 6 In the cases of the octagon (1.5) and localization (1.6) symbols we similarly find (cf. (1.7) and (1.9)) where f ℓ (g) stand for 1/g n corrections (1.10) and dots denote non-perturbative corrections (1.12). Notice that the leading O(g) term has an opposite sign in (2.34) and (2.35). The Widom-Dyson constants B oct ℓ and B loc ℓ are given by (2.24). The O(log g) term in both expressions is generated by the Fisher-Hartwig singularity of the symbol. According to (2.21), the corresponding parameters β oct and β loc are integer and, therefore, the perturbative functions f oct ℓ (g) and f loc ℓ (g) are given by Borel non-summable series. For instance, for ℓ = 2 and β = −1 we find from (1.10), (2.25) and (2.11) One can verify that, in agreement with (2.27), (2.26) and (2.12), all terms in this expression involving log 2 π can be eliminated by changing the expansion parameter to g ′ = g + log 2 π . There exists an interesting relation between the two different perturbative functions f oct ℓ (g) and f loc ℓ (g). It follows from the identity I loc n = −I oct n in (2.11). A close examination of (2.27) shows that the function f ℓ (g) is formally invariant under transformation g → −g and I n → −I n (with n = 0, 1, 2, . . . ). This leads to We recall that both functions depend on ℓ β = ℓ + β and the shift ℓ → ℓ + 2 on the right-hand side is needed to compensate the difference β oct − β loc = 2. Because the functions on both sides of (2.37) suffer from Borel singularities, this relation is rather formal and it should be understood as an equality between the expansion coefficients in the two series. Equivalently, upon the Borel transform (1.11), Eq. (2.37) leads to the relation which maps the Borel singularities of the two functions into each other. In particular, the leading Borel pole of B oct ℓ (σ) for σ > 0 is in one-to-one correspondence with the pole of B loc ℓ+2 (σ) for σ < 0 closest to the origin.

Non-perturbative corrections to SAK formula
Let us now employ the method developed in [40] to compute non-perturbative corrections to F ℓ (g) in (1.7) at strong coupling. They are described by the function ∆f ℓ (g) which is expected to have a general form (1.12).
To find the dependence of ∆f ℓ (g) on the coupling constant g and non-negative integer ℓ, we shall examine two functions related to F ℓ (g). The first one is where we used the determinant representation (1.1) of F ℓ (g) in terms of matrix K. The second one is F ℓ+2 (g) − F ℓ (g), i.e. the difference of functions with shifted indices. It follows from the definition (1.1) and (1.2) that exp(F ℓ+2 ) is given by the determinant of the matrix (1 − K) with the first row and column removed, exp(F ℓ+2 ) = det(1 − K)| n,m≥2 . Then, the ratio of the determinants exp(F ℓ+2 ) and exp(F ℓ ) can be evaluated using Cramer's rule as Having computed (3.1) and (3.2) at strong coupling, we can obtain relations for g∂ g ∆f ℓ (g) and ∆f ℓ+2 (g) − ∆f ℓ (g). It proves convenient to introduce an auxiliary function Γ(x, y) [40,43] Γ(x, y) = 1 The rationale for defining this function is that both quantities (3.1) and (3.2) can be expressed in terms of Γ(x, y) in a simple way.
Indeed, taking into account (1.2), we can rewrite (3.1) as To find the ratio (3.2), we examine asymptotic behaviour of the function Γ(x, y) for small x and/or y. In both cases, the leading contribution to (3.3) comes from the Bessel functions with the minimal index where in the first relation we applied (2.15) and replaced J ℓ+1 (x) ∼ x ℓ+1 . Combining together the two limits, x → 0 and y → 0, we find from (3.3) where dots denote terms suppressed by powers of x and y. Thus, the ratio (3.2) can be obtained from the leading behaviour of Γ(2gx, 2gy) at small x and y.
In what follows, we first determine the function Γ(x, y) and, then, apply (3.4) and (3.6) to compute the two quantities defined in (3.1) and (3.2). We use (1.2) together with (3.3) to find that Γ(x, y) verifies the (infinite) system of integral equations The function Γ(x, y) has to satisfy the additional conditions that follow from its definition. According to (3.3), it is given by the product of the entire function γ(x, y) and the meromorphic function (1 − χ (x/(2g))) /y. It follows from (2.13) that the latter function vanishes at x = ±4πigx n and has poles at x = ±4πigy n (with n = 1, 2, . . . ). The function Γ(x, y) inherits these properties, e.g., We show below that the integral equation (3.7), supplemented with the information about analytical properties of Γ(x, y), is sufficient to construct its solution.
As was explained in Section 2, the perturbative part of the strong-coupling expansion (1.7) takes a different form for integer and half-integer β. This property can also be seen from (3.7). We verify using (3.3) and (2.18) that the function Γ(x, y) satisfies where we took into account the Bessel function identity J k (−x) = (−1) k J k (x). As a result, the integrand on the left-hand side of (3.7) is an odd/even function of x for even/odd 2β. For half-integer β, this allows us to extend the integration in (3.7) to the whole real x-axis and evaluate the integral by residues at the poles of Γ(x, y). For integer β, we show below that the relation (3.7) leads to a Riemann-Hilbert problem for Γ(x, y). Let us consider separately the cases of half-integer and integer β.

'Easy' case: half-integer β
It follows from (3.9) that Γ(−x, y) = (−1) ℓ Γ(x, y) and, therefore, Γ(x, y) is an even (odd) function of x for ℓ even (odd). This allows us to rewrite the relation (3.7) as Solving the integral equation (3.10), we follow the same steps as in [40]. Details of the calculation can be found in Appendix C. We start with performing a Fourier transform of Substituting (3.11) into (3.10) and exchanging the order of integration, we find that the integral over x vanishes for k 2 > 1. Then, solving the relation (3.10) we can determine the functionΓ(k, y) for k 2 ≤ 1. To find the functionΓ(k, y) for k 2 > 1, we invert the relation (3.11) and replace Γ(x, y) with its representation (3.3) as a product of 1 − χ(x/(2g)) and an entire function γ(x, y)/y. Computing the integral over x by residues, we obtain the representation forΓ(k, y) at k 2 > 1 as a sum over poles of the function 1 − χ(x/(2g)).
Finally, splitting the integration region in (3.11) into k 2 ≤ 1 and k 2 > 1, we replacẽ Γ(k, y) with its expressions in each of these regions and arrive at the following result for Γ(x, y) (see Appendix C) for even ℓ, and for odd ℓ.
The relations (3.12) and (3.13) are valid for β ≥ −(ℓ + 1)/2 in which case the function Γ(x, y) remains regular at small x, see (3.5). 7 The first two lines in (3.12) and (3.13) come from integration over k 2 ≤ 1 in (3.11) and the last line from k 2 > 1. Both expressions involve the polynomials as well as some functions a n (y) and c j (y). A distinguished property of p n (x) is that the x-dependent coefficients in front of a n (y) in (3.12) and (3.13) are regular for x → 0. The functions a n (y) describe the contribution of zero modes of the integral equation (3.10) whereas c j (y) define the residue of Γ(x, y) at the poles x = ±4πigy j . Both sets of functions can be found by requiring the functions (3.12) and (3.13) to satisfy the relations (3.5) and (3.8). To illustrate this, we take β = −1/2 and ℓ = 0, 1.

Special solutions
For β = −1/2 and ℓ = 0 the relation (3.12) simplifies as It does not involve the zero modes contribution and satisfies (3.5).
To find the functions c j (y) we impose the condition (3.8). Because (3.15) is an even function of x, it is sufficient to require that Γ(x, y) vanishes for x = 4πigx n . This leads to the following (infinite) system of equations where n ≥ 1 and x n are positive. At strong coupling, the expression on the second line is exponentially suppressed. This suggests to look for a solution to (3.16) in the form of expansion in powers of e −8πgxn . Arranging the roots x n in ascending order, 0 < x 1 < x 2 < . . . , we find that the leading contribution to c j (y) comes from the smallest root where dots denote subleading corrections of the form e −8πg(m 1 x 1 +m 2 x 2 +... ) with m i nonnegative integer. If all roots are multiples of the smallest one, e.g., x n /x 1 are positive integer, the expansion in (3.17) runs in powers of e −8πgx 1 . It is easy to see that this is indeed the case for the symbols in (2.10). Combining together (3.15) and (3.17) we obtain where Γ (0) (x, y) is given by (3.15) with c j (y) replaced by the leading result c (0) j (y). In a similar manner, Γ 1) (x, y) is given by the sum in (3.15) with c j (y) replaced by c (1) j (y). The dots in (3.18) denote subleading exponentially small corrections.
Substituting (3.17) into (3.16) and matching the terms on both sides of the relation we find , (3.19) where n ≥ 1. The left-hand side of both relations involve a Cauchy matrix 1/(x n − y j ).
Inverting this matrix, we can determine the functions c The above relations were obtained for β = −1/2 and ℓ = 0. The same analysis can be carried out for β = −1/2 and ℓ = 1. We can start with (3.13) and go through the same steps to find that the function Γ(x, y) takes the same form (3.18). The only difference as compared to the previous case is that the coefficients c (0) j (y) and c (1) j (y) satisfy the system of equations that differ from (3.19) by signs in front of various terms on the right-hand side.
Using the obtained expressions for c (0) j (y), it is possible to express Γ (0) (x, y) for β = −1/2 and ℓ = 0, 1 in terms of the function Φ(x) defined in (2.18), see Appendix D. This leads to We verify that this function satisfies the relation (3.8) for g ≫ 1. For x = ±2πix n each term on the right-hand side of (3.20) is either exponentially small at strong coupling, or proportional to Φ(−2πix n ) = 0. Note that the function (3.20) is regular for x = ±y. Similarly, the function Γ (1) (2gx, 2gy) admits the following representation where we introduced the notation .
The additional factor in the denominator ensures that F (−2πix 1 ) is different from zero. As above, we verify that Γ (1) (2gx, 2gy) vanishes for x = ±2πix n and n ≥ 2. For x = ±2πix 1 the function (3.21) is different from zero and scales as e 4πgx 1 . Going back to (3.18) we observe that its contribution to Γ(±4πigx 1 , 2gy) cancels against the O(e −4πgx 1 ) contribution coming from Γ (0) (±4πigx 1 , 2gy). We would like to emphasize that the relations (3.20) and (3.21) hold only for β = −1/2 and ℓ = 0, 1. In particular, they automatically satisfy the condition (3.5). For β > −ℓ/2 the situation is different. For arbitrary a n (y) and c j (y) the functions (3.12) and (3.13) scale at small x as O(x 0 ) and O(x), respectively. The relation (3.5) implies that the first 2β + ℓ + 1 terms of the small x expansion of both functions have to vanish. Imposing this condition on (3.12) and (3.13) allows us to express the zero modes a n (y) in terms of functions c j (y) (with j ≥ 1) and, in addition, obtain non-trivial relations for infinite sums j c j (y)/y m j with m = 1, . . . , β + 1/2. Requiring the resulting expression for Γ(x, y) to verity (3.8) and going through the same steps as above, we find that at strong coupling c j (y) and Γ(x, y) have the same general form as before, Eqs. (3.17) and (3.18). Important difference as compared with the previous case is that for β = −1/2 + p the expansion coefficients c (0) j (y) and c (1) j (y) in (3.17) scale at strong coupling as O(g p ). They satisfy a system of linear relations analogous to (3.19) supplemented with the additional relations for the sums j c j (y)/y m j (with m = 1, . . . , p) mentioned above. Solving these relations we can obtain the expressions for Γ (0) (2gx, 2gy) and Γ (1) (2gx, 2gy) that are valid for arbitrary half-integer β ≥ −1/2 and non-negative ℓ.

Strong coupling expansion
Taking into account (1.7), (1.12), (2.17) and (2.25), we expect that the strong-coupling expansion of g∂ g F ℓ looks like where dots denote terms suppressed by powers of 1/g as well as subleading exponentially small corrections. Using the obtained expressions for the function Γ(x, y) in (3.18), we can apply (3.4) to compute (3.23). In this way, we should be able to reproduce the first two terms on the right-hand side of (3.23) and, most importantly, identify the values of the parameters c 1 , n 1 and x 1 defining the leading non-perturbative correction.
Let us start again with β = −1/2 and ℓ = 0, 1. We use (3.20) to get Substituting this expression into (3.4) and taking into account (2.9), we find that the first term on the right-hand side of (3.24) gives rise to (−2gI 0 ) term in (3.23). The second term in (3.24) yields the O(g 0 ) correction to (3.23) Here in the first relation we replaced the symbol (1 − χ(x)) with its expression (2.18), introduced notation for φ(x) = ∂ x log Φ(x) and extended the integration to the whole real axis. In the second relation, we took into account that φ(x) has poles in the lower half-plane and deformed the integration contour to the upper half-plane to become an arc of infinite radius, x = R e iα with R → ∞ and 0 ≤ α ≤ π. It follows from (2.19) that φ(x) ∼ −β/x on this arc and the integral (3.25) can be easily evaluated at β = −1/2. Similarly, the contribution of the last term in (3.24) to (3.4) can be written as The last factor arises because the sum of the two terms inside the brackets on the second line of (3.24) is regular for x → 0 but each term separately has a pole 1/x. The integral (3.26) can be evaluated by closing the integration contour to the lower half-plane and by picking up residue at the poles. Replacing the functions χ(x) and Φ(x) with their expressions (2.13) and (2.18), we find that the integrand has simple poles at x = 0, x = −4πiy n and double poles at x = −4πix n . The residue at the pole x = 0 yields a constant (−1) ℓ /4, whereas the contribution of the two remaining sets of poles is exponentially small at strong coupling. The residue at the double pole is enhanced by the factor of g and the leading contribution comes from the double pole closest to the origin, x = −4πix 1 . As a result, the integral (3.26) is given by where the notation was introduced for Here the function F (x) is given by (3.22).
The first term inside the brackets (3.27) contributes to the O(g 0 ) term in (3.23) whereas the second term defines a non-perturbative correction to g∂ g F ℓ . Similar contribution also comes from the second term in (3.18) defined in (3.21). Substituting e −8πgx 1 Γ (1) (2gx, 2gx) into (3.4) and carrying out the integration, we find that its contribution to g∂ g F ℓ scales at strong coupling as e −8πgx 1 and, therefore, is subleading compared to (3.27).
Finally, combining together (3.25) and (3.27) we arrive at where dots stand for subleading exponentially suppressed corrections. We recall that this relation holds for β = −1/2 and ℓ = 0, 1. It is easy to see that for these values of the parameters the first two terms on the right-hand side of (3.23) coincide with the analogous terms in (3.29).
Matching (3.29) to (1.7) we identify the leading non-perturbative correction to F ℓ (g) for β = −1/2 and ℓ = 0, 1 ( 3.30) It is straightforward to generalize the relation (3.30) to arbitrary non-negative ℓ and halfinteger β. As we will see in a moment, the non-perturbative correction to the difference of functions ∆f ℓ+2 −∆f ℓ scales as O(g −1 e −8πgx 1 ) and, therefore, it is suppressed by the factor of g as compared with (3.30). This means that the leading non-perturbative correction cancels in the difference ∆f ℓ+2 −∆f ℓ and, as a consequence, the relation (3.30) holds for an arbitrary ℓ.
To restore the β-dependence of (3.30), we recall (see footnote 7) that sending one of the roots of the function (2.18) to zero, say x i → 0, generates the shift β → β + 1. It is easy to see from (3.28) that under this transformation Λ(x 1 ) → −Λ(x 1 ). Thus, in order to restore the β-dependence of ∆f ℓ , it is sufficient to insert (−1) β+1/2 on the right-hand side of (3.30) (3.31) This relation defines the leading non-perturbative correction to F ℓ (g) for half-integer β.
To find a general expression for ∆f ℓ+2 − ∆f ℓ , we repeated the calculation of F ℓ+2 − F ℓ for ℓ = 2, 3 and β = p − 1/2 with p = 1, 2, 3. In this way, we reproduced the first two terms inside the brackets on the right-hand side of (3.32) and identified the leading nonperturbative correction to ∆f ℓ+2 − ∆f ℓ . We found that this correction can be obtained by applying the following transformation to the O(1/g) perturbative term in (3.32) The resulting relation for the leading non-perturbative correction is then Notice that the expression on the right-hand side depends on the sum ℓ + β. Recall that the perturbative function (2.27) has the same property.

'Hard' case: integer β
The main difference as compared to the previous case is that, as follows from (3.9), Γ(−x, y) = −(−1) ℓ Γ(x, y) and, therefore, Γ(x, y) is an odd (even) function of x for ℓ even (odd). As will see below, this entails a change of the properties of both perturbative and non-perturbative expansions of F ℓ . As we have seen in the previous subsection, for half-integer β the non-perturbative corrections depend on the sum ℓ + β. We assume below that the same is true also for integer β. Indeed, this property holds for the perturbative function f ℓ (g) that is given by a Borel non-summable series for integer β. The Borel singularities induce an exponentially small ambiguous contribution to the perturbative function f ℓ (g) which is canceled in the sum of f ℓ (g) with ∆f ℓ (g). Because the former depends on ℓ + β, the same should be true for ∆f ℓ (g). Taking the advantage of this property, we can restrict our analysis to ℓ = 0, 1 and an arbitrary integer β.
Let us start with β = 0 and then extend consideration to arbitrary integer β > 0. For ℓ = 0, 1 an infinite system of equations (3.7) can be cast into a compact form To recover the relations (3.7), it is sufficient to replace the trigonometric functions with their Bessel series expansion and match the coefficients on both sides. It is important to emphasize that the relations (3.38) only hold for −1 < u < 1.
Applying the Fourier transform (3.11) we get from (3.38) where the integral is defined using the principal value prescription. The functionΓ(k, y) has a definite parity,Γ and its analytical properties are in one-to-one correspondence with analogous properties of the function Γ(x, y) described at the beginning of this Section. In particular, we can show, following [40], that because Γ(x, y) has an infinite sequence of simple poles at x = ±4πiy j , the functionΓ(k, y) has the following form for k > 1 Here the functions c j (y) define the residue Γ(x, y) at x = ±4πiy j . The resulting Riemann-Hilbert problem, Eqs. (3.39) - (3.41), is similar to that discussed in [40]. Going along the same lines as in that paper we find These relations are valid for |k|≤ 1. They allow us to express the functionΓ(k, y) for k 2 < 1 in terms of the same functionΓ(p, y) defined for p 2 > 1. ReplacingΓ(p, y) in the last relation with (3.41), we obtain a representation forΓ(k, y) in terms of an infinite set of functions c j (y) (with j ≥ 1). As in the previous subsection, we can determine these functions from the requirement of Γ(x, y) to satisfy the additional conditions (3.6) and (3.8). Substituting (3.41) and (3.42) into (3.11) we find after some algebra where J n , I n and K n are Bessel functions. We verify that at small x and y these expressions satisfy the relation (3.5) for β = 0. For positive integer β, the condition Γ(x, y) ∼ x 2β+ℓ+1 as x → 0, leads to the additional relations for the coefficient functions c j (y). For β = 0 and ℓ = 0, we substitute (3.43) into (3.8) to obtain the system of linear relations for the coefficient functions c j (y) 0 = 4πgx n J 0 (y) + r n yJ 1 (y) 16π 2 g 2 x 2 n + y 2 where n ≥ 1 and we defined For β = 0 and ℓ = 1, combining together (3.44) and (3.8), we find thatc j (y) satisfy similar relations. The relation (3.45) can be further simplified at strong coupling. For g ≫ 1, the ratio K 0 (4πgy j )/K 1 (4πgy j ) is given by an asymptotic sign-alternating series in 1/g. Replacing it by the leading term we find (for n ≥ 1) It is instructive to compare the last relation with (3.16). We note that the exponential functions are replaced by a linear combination of the Bessel functions and an exponentially small factor of exp(−8πgx n ) with Here both terms are accompanied by series in 1/g and dots denote terms suppressed by powers of e −8πgxn .
Because the left-hand side of (3.48) is a well-defined real function of the coupling constant g, the appearance of the factor 'i' on the right-hand side of (3.48) may be surprising. It is related to the fact that the series entering the second term in (3.48) suffer from Borel singularities. To make it well-defined, we can move the coupling constant slightly above the real axis by g → g + i0. This transformation amounts to deforming the integration contour in the Borel plane in the vicinity of the singularities (see (1.11)). It induces an imaginary contribution to the second term in (3.48) which cancels against a similar contribution from the first term in (3.48) in such a way that their sum remains real and well-defined (see Appendix E for details).
In a close analogy with (3.17), we look for solutions to (3.47) in the form where dots denote corrections suppressed by powers of 1/g. Substituting (3.49) into (3.47) and taking into account (3.48), we find that the coefficient functionsc . (3.50) These relations are similar to (3.19). As mentioned above, they can be obtained from (3.19) by replacing exponential functions with a linear combination of the Bessel functions e ±iy → J ± (y) ≡ J 0 (y) ± iJ 1 (y) . It is therefore not unexpected that the solutions to (3.19) and (3.50) are related to each other by the same transformation.
Solving the system of equations (3.50) as well as the analogous system for β = 0 and ℓ = 1, we can find the functions (3.43) and (3.44) at strong coupling where dots stand for terms suppressed by powers of 1/g and by exponentially small factors e −8πgxn (with n ≥ 1). Here the functions Γ (0) (x, y) and Γ (1) (x, y) are obtained from (3.20) and (3.21) by applying the transformation (3.51) where the functions Φ, F and J ± were introduced in (2.18), (3.22) and (3.51), respectively. The important difference between (3.18) and (3.52) is that, in virtue of (3.48), the perturbative series in 1/g on the right-hand side of (3.52) contains Borel singularities. As discussed above, these singularities are regularized by replacing g → g + i0.
We can now apply (3.53) to determine the strong-coupling expansion of the observable F ℓ . Following the same steps as in Section 3.1, we examine (3.53) at the coinciding point y = x. Using the properties of Bessel functions (3.51), we find from (3.53) for g ≫ 1 (3.54) The analogous expression for Γ (1) (2gx, 2gx) scales as O(1/g 2 ) and produces a subleading contribution to (3.4). Compared to (3.24), the two terms inside the brackets in (3.54) have an opposite relative sign. As we will see in a moment, this has important consequences for the properties of F ℓ . Substituting (3.54) into (3.4) we find that the first term on the right-hand side of (3.54) yields (−2gI 0 ), whereas the second one leads to the integral that is similar to (3.25). Because ∂ x log Φ(x) ∼ −β/x at large x (see (2.19)), this integral vanishes for β = 0.
The contribution of the last term in (3.54) to (3.4) can be written as where the integration contours C + and C − start at the origin and go to +∞ slightly above (C + ) and below (C − ) the real axis. Here we took into account that the two terms inside the brackets on the last line of (3.54) have poles located along the imaginary axis and rotated the integration contour as x → −ix and x → ix, respectively. The integrals on the left-hand side of (3.55) define two branches of the function ϕ(g). After the change of variable x = σ/(4g), they take the form of the Borel transform (1.11) Using (2.18) we find that B ϕ (σ) is a meromorphic function with poles located at real σ. As a consequence, the function ϕ(g) has a cut running along the real axis in the complex g-plane. 8 The discontinuity of this function across the cut is given by an (infinite) sum over the residues of B ϕ (σ/g) at the poles located at positive σ. At large positive g the leading contribution comes from the double pole at σ = 8πgx 1 where Λ(x 1 ) was defined in (3.28). 9 Combining together the above relations, we obtain from (3.4) where dots denote subleading corrections including those coming from Γ (1) (2gx, 2gy). The sum of two terms on the right-hand side of (3.58) is a real function of g, whereas each term separately develops an exponentially small imaginary part (3.57). Applying (3.58) and (3.57) we can obtain two other representations of the same quantity The three representations, Eqs. (3.58) and (3.59), are equivalent. This illustrates a universal feature of the strong-coupling expansion that has been mentioned previously. Namely, the definition of non-perturbative, exponentially small corrections to g∂ g F ℓ (g) depends on a regularization that one employs to integrate through the Borel singularities in (3.56). According to (3.56), the functions ϕ(g+i0) and ϕ(g−i0) are obtained by shifting the integration contour in the Borel plane slightly above and below the real axis, respectively. In the representation (3.58), the Borel poles are integrated using the principal value prescription. Comparing (3.59) with (3.29) we notice that non-perturbative corrections to both expressions involve the same quantity 8(−1) ℓ gπx 1 e −8πgx 1 Λ(x 1 ). For half-integer β, the perturbative series f ℓ (g) is well-defined and the coefficient in front of the non-perturbative term in (3.29) can be determined unambiguously. In contrast, for integer β this coefficient depends on the regularization of the Borel singularities of the perturbative function ϕ(g).
Integrating the second relation in (3.59) and matching F ℓ to (1.7), we can find the leading non-perturbative correction (3.60) We would like to emphasize that this relation should be supplemented with the analogous relation for the perturbative function where the '+i0' prescription on the right-hand side specifies a deformation of the integration contour in the vicinity of the Borel singularities in (3.56). If we used '−i0' prescription, the expression on the right-hand side of (3.60) would have an opposite sign. Note that the relation (3.60) and its complex conjugate coincide with (3.31) evaluated at β = 0 and √ −1 = ±i. Finally, we can find the difference ∆f ℓ+2 −∆f ℓ by evaluating non-perturbative corrections to the ratio (3.2). To obtain D ℓ using (3.6), we examine the leading asymptotic behaviour of the function Γ(2gx, 2gy) given by (3.52) and (3.53) for x, y → 0 (3.62) The first term inside the parentheses in both relations is real whereas the second one is pure imaginary. The situation here is similar to the one we encountered in (3.59). The non-perturbative correction to (3.62) corresponds to a particular regularization of the perturbative function f ℓ (g) → f ℓ (g +i0), see Eq. (3.61). It amounts to deforming the integration contour in the Borel transform (1.11) slightly below the real axis, thus avoiding the Borel poles. Matching (3.62) to (3.6) we find where Λ(x 1 ) is given by (3.28). This relation is valid for β = 0 and ℓ = 0, 1. As in the previous case of half-integer β, we repeated the calculation of D ℓ for β = 1, 2, 3, 4 and found that the leading non-perturbative correction to D ℓ can be obtained through the transformation (3.36). This leads to up to corrections suppressed by 1/g. Notice that the leading correction (3.60) cancels on the right-hand side of (3.64).

Degenerate symbols
The obtained relations for the non-perturbative corrections (3.60) and (3.64) are valid for the symbol of the form (2.18). A distinguished feature of (2.18) is that (1 − χ(x)) has simple roots at x = ±2πix n . Examining different expressions of the symbols in (2.10) we observe that this condition is satisfied for 1−χ BES (x) whereas the two remaining symbols, 1−χ oct (x) and 1 − χ loc (x) have double roots.
To define the symbols with double roots it is sufficient to take the limit of (2.18) as x 2 → x 1 . We recall that deriving the non-perturbative corrections (3.60) and (3.64), we neglected subleading corrections to ∆f ℓ (g) suppressed by powers of exp(−8πx n ) with n ≥ 2. For x 2 → x 1 we also have to retain corrections proportional to e −8πx 2 . They are given by the same expressions (3.60) and (3.64) with x 1 replaced by x 2 . For instance, we get from (3.60) where Λ(x 2 ) is obtained from (3.28) by replacing x 1 ↔ x 2 . In particular, for where the notation was introduced for (3.67) Combining the above relations we get from (3.65) in the limit As compared with (3.60), this expression contains the additional factor of g. In the similar manner, we get from (3.64) The relations (3.68) and (3.69) hold for integer β. They can be obtained from the analogous relations (3.60) and (3.64) by replacing For half-integer β, we can apply the same transformation to (3.31) and (3.37) to get the corresponding expressions for the function ∆f ℓ (g) for the symbol with double roots.

Non-perturbative corrections at strong coupling
Let us summarize the obtained results for the non-perturbative corrections to (1.7) at strong coupling.
For half-integer β, perturbative f ℓ (g) and non-perturbative ∆f ℓ (g) terms in (1.7) are well-defined functions of the coupling constant g. The function f ℓ (g) takes into account perturbative corrections in 1/g and it takes the form (2.28). The function ∆f ℓ (g) describes non-perturbative, exponentially small corrections to (1.7). It satisfies the relations (3.31) and (3.37).
For integer β, the function f ℓ (g) is given by a Borel non-summable series (1.10) and requires a regularization. Combining together the relations (1.11), (3.61) and (3.56), we find that its Borel transform satisfies (3.71) According to (3.56), the function B ϕ (σ) has double poles at σ = 8πx n and σ = −8πy n (with n ≥ 1). As a consequence, B f (σ) has the following schematic form where we displayed the contribution of the two poles closest to the origin. Substituting (3.72) into (1.11) and expanding the integral in powers of 1/g, we find As expected, the expansion coefficients grow factorially. For x 1 < y 1 and x 1 > y 1 the large order behaviour is controlled, respectively, by the first and second term inside the brackets. In both cases, the series (3.73) suffers from Borel singularities. Due to the presence of the pole in (3.72) at x = 8πx 1 , the integral (1.11) is not welldefined and requires a regularization. Different ways of deforming the integration contour in (1.11) in the vicinity of the pole lead to different results for the perturbative term f ℓ (g). They differ from each other by exponentially small terms O(e −8πgx 1 ) which are proportional to the residue at the pole. The dependence on the regularization disappears in the sum of perturbative and non-perturbative terms (1.9). We demonstrated that this sum can be written in three equivalent ways Here f ℓ (g +i0) and f ℓ (g −i0) are given by the Borel transform (1.11) in which the integration contour is shifted, respectively, slightly below and above the Borel poles. The additional exponentially small term ∆f ℓ (g) = 1 2 (f ℓ (g − i0) − f ℓ (g + i0)) satisfies the relations (3.60) and (3.64).

Physics applications
Let us now specify the non-perturbative function ∆f ℓ (g) for three physically relevant cases of the symbol χ(x) in (1.3), (1.5) and (1.6).

Easy case: BES
Let us first consider the symbol (1.3) with β = −1/2. According to terminology adopted in the previous subsection, this corresponds to an 'easy' case.
The relations (3.31) and (3.37) for the non-perturbative function ∆f ℓ (g) depend on x 1 and Λ(x 1 ). We recall that x 1 is the solution to Φ BES (−2iπx 1 ) = 0 closest to the origin. We apply (2.20), (3.22) and (3.28) to find that x 1 = 1 2 and We then use (3.31) and (3.37) to find the leading non-perturbative corrections as It follows from the first relation that the leading non-perturbative correction to F BES ℓ is 1 2 (−1) ℓ e − √ λ . It is easy to verify that for ℓ = 0, 1, 2 this result is in an agreement with the exact relations (2.32).
In a similar manner, we apply the second relation in (3.76) together with (2.31) to obtain the ratio (3.2) for ℓ = 0 where the last term denotes subleading non-perturbative corrections. This relation should be compared with the exact expression of D 0 that follows from (2.32) We observe a perfect agreement.

Hard cases: octagon and localization
Let us now consider the symbols (1.5) and (1.6). As compared to the previous case, there are two important differences. First, because the strength of the Fisher-Hartwig singularity is integer in this case, β oct = −β loc = 1 (see Eq. (2.21)), the perturbative correction to F ℓ has a more complicated form (3.73). Matching (2.21) to (2.18), we identity the leading root (x 1 ) and pole (y 1 ) of the two symbols as It follows from (3.73) that the perturbative series f oct ℓ (g) is sign alternating, 10 whereas f loc ℓ (g) has the expansion coefficients of the same sign. This property is in agreement with the explicit expressions (2.36) and (2.37).
Secondly, the leading root of the symbols (2.21) is double degenerate and, as a consequence, the non-perturbative corrections to F ℓ satisfy the relations (3.68) and (3.69). We apply (3.67) and (2.21) and identify the corresponding non-perturbative parameters as Then, it follows from (3.68) that the leading non-perturbative correction satisfies Note that the leading large g terms in (3.81) cancel out in the difference. It is important to emphasize that the relations (3.81) and (3.82) hold for large positive g. In fact, the functions ∆f oct ℓ (g) and ∆f loc ℓ (g) have different asymptotic behaviour at large positive and negative g due to the Stokes phenomenon.
Recall that the perturbative series for f loc ℓ (g) and f oct ℓ (g) are related to each other as (2.37). Due to the presence of Borel singularities in both series, this relation is formal. As discussed above, we can regularize these singularities by deforming the integration contour in (1.11) in the vicinity of the Borel poles. In this way, we get . This discontinuity yields the non-perturbative correction ∆f ℓ (g) = 1 2 (f ℓ (g − i0) − f ℓ (g + i0)), i.e. (3.83) leads to ∆f oct ℓ (g) = −∆f loc ℓ+2 (−g) .

(3.84)
Combined together with (3.81) and (3.82), this relation allows us to determine the asymptotic behaviour of the non-perturbative functions at large negative g. Moreover, substituting (3.83) into the first relation in (3.74) we find that the sum of the perturbative and non-perturbative contributions to F oct ℓ (g) and F loc ℓ (g) are related to each other as Thus, the functions ∆F oct ℓ (g) and ∆F loc ℓ+2 (g) can be identified as two branches of the same function defined for negative and positive g, respectively.

Applications to N = 2 superconformal models
Our aim in this section is to compute the leading non-trivial corrections to special observables in N = 2 four-dimensional superconformal models that are planar-equivalent to N = 4 SYM. These observables are controlled by the localization matrix model and, as a result, can be expressed in terms of the semi-infinite matrix K nm defined in (1.2) (see [18][19][20][21] and Appendix A).
The relevant N = 2 models may be split into two classes I and II: I. Models with gauge group SU(N) or Sp(2N) and matter content summarized in Ta  rank-2 symmetric and antisymmetric representations. Guided by geometrical engineering, these models are expected to be holographically dual to IIB superstring theory on orientifolds/orbifolds of the type AdS 5 × S 5 /G [46,47]. They correspond to a combination of 2N D3-branes together with an orientifold O7 plane and, in models with n F = 0, also with several D7-branes.
II. Quiver theories (that we call Q L ) with gauge group SU(N) ⊗L , bi-fundamental matter, and equal gauge couplings. They are orbifold projections of the N = 4 SYM theory and are dual to type IIB superstring on AdS 5 × S 5 /Z L orbifold [48][49][50][51].

Observables in terms of Bessel operator
The three simplest observables in these models are the free energy on S 4 , half-BPS circular Wilson loop, and some correlators of chiral primary operators. They can be computed using localization matrix model techniques. It turns out that the leading non-trivial corrections to them can be expressed in terms of the Bessel operator (1.2). The interaction potential of the localization matrix model is given by an (infinite) sum of terms weighted by powers of the 't Hooft coupling. The evaluation of the observables is straightforward in weak coupling expansion, but the strong-coupling expansion poses a non-trivial problem. The latter limit is of main interest from the point of view of establishing correspondence with dual string theory, i.e. for tests of the AdS/CFT correspondence.
Explicitly, semi-infinite matrix in (1.1) appearing in the context of the above N = 2 superconformal models is given by [18] It is easy to see that, upon change of the variable t → t/ √ λ, the matrix X nm coincides with the Bessel matrix (1.2) with the symbol χ = χ loc (x) given in (1.6) (4.2)

Free energy
The free energy is defined as F = − log Z, where Z is the partition function of gauge theory on S 4 . It is a function of λ and N. Planar equivalence implies that at large N the free energy of the above N = 2 models is equal to the free energy of N = 4 SYM theory 11 The leading O(N 2 ) term cancels in the difference ∆F = F N =2 − F N =4 which is thus a genuine N = 2 quantity. More precisely, let us define (4.4) We will be interested in the leading N → ∞ limit of ∆F M . This quantity has been studied previously in the SA [20], FA and FA [21], and Q 2 models [19].
In the SA and Q 2 models, one finds the following explicit representations for the leading correction to the free energy in terms of the Bessel operator observable F ℓ defined as in (1.1), (1.2), (1.6) and (4.2) The relation (4.5) was proved in [20], and (4.6) is derived below in Appendix F. From a string theory argument for the non-planar correction to the half-BPS Wilson loop [53] and its relation to the free energy [19], one expects that in both cases (4.5) and (4.6) the leading term at strong coupling should take the form The constant C Q 2 in the L = 2 quiver model was estimated in [19] by a Monte Carlo numerical simulation (in a moderate range λ < 450) with the result In the SA model, an analytical determination of C SA was attempted in [20] by considering the leading-order (LO) large λ contribution to the matrix elements of (4.1) at ℓ = 2 before computing the determinant. As a result, the expected λ 1/2 scaling of ∆F SA was reproduced with the coefficient being A numerical high-precision resummation of the weak-coupling expansion of ∆F SA was performed in [20] using an improved conformal mapping Pade' analysis. While the large λ scaling exponent 1/2 in (4.7) was confirmed, the Pade' estimate for its coefficient was substantially smaller C SA Pade ≃ 0.12 than in (4.9), so the final picture was not totally satisfactory. Turning to the FA and FA models with fundamental hypermultiplets, one finds that the large N expansion of the free energy contains both odd and even powers of 1/N. In the FA model with gauge group SU(N) one gets [21] (4.10) Here the leading O(N) term F 1 (λ) can be found in a closed form [21] (4.11) It is straightforward to work out its strong coupling expansion [21] F 1 (λ) = log 2 Note that the second line of (4.12) contains only one perturbative term and that the coefficient of the non-perturbative correction 12 is real (cf. (4.32) and (4.34) below). According to (4.10), the function F 1 (λ) determines the part F 2 (λ) of the subleading O(N 0 ) correction. The remaining part of F 2 (λ) is the same as in free energy in the SA theory given by (4.6).
In the FA model with the Sp (2N) gauge group, the large N expansion of the free energy ∆F FA is much simpler (it is essentially determined by F 1 (λ)) and thus does not involve the Bessel operator. 13 In this case one may directly identify the non-perturbative corrections to ∆F FA that are exponentially suppressed at large λ and are partially controlled by resurgence properties [21].

Non-planar correction to half-BPS circular Wilson loop
In the above N = 2 models the localization yields a matrix model representation for the expectation value of the (suitably defined) half-BPS circular Wilson loop W . The planar equivalence implies that W M in these models has a universal large N limit where I 1 is the modified Bessel function. The subscript '0' stands for the planar limit, i.e. the leading term in the large N expansion. The leading non-planar correction to W M defines the model-dependent function q M (λ) 14 (4.14) Remarkably, one can establish simple relations between q M and the leading correction ∆F M to the free energy [19,20] ∆q 16) 12 We omit terms with higher odd powers of e − √ λ , see [21]. 13 The technical reason for this simplification is that the localization matrix model for FA model has the interaction potential which is free from double-trace terms, cf. (A.4). As a consequence, ∆F FA can be found analytically and an exact resummed expression is available for the leading strong coupling terms at each order in the 1/N expansion [21].
14 For a discussion of the Wilson loop in the models FA and FA with the matter in the fundamental representation see [21].
. . at strong coupling. Combining these relations together with (4.5) and (4.6), one can express the function q M (λ) in terms of the Fredholm determinant F ℓ of the Bessel operator.
The string theory argument suggests [53] that q M (λ) should scale at strong coupling as λ 3/2 . The relations (4.15) and (4.16) then imply that the coefficient of λ 3/2 should be proportional to C M in (4.7).

Correlation functions of chiral operators
In contrast to the free energy and circular Wilson loop discussed above, the correlation functions of some chiral operators in N = 2 models of class I and II differ from their counterparts in N = 4 SYM already at the leading large N order.
In the class I N = 2 models this happens for correlators of chiral primary operators O k (x) = tr ϕ k (x) involving odd power k of the complex scalar ϕ(x) from the N = 2 vector multiplet. Defining the ratio of two-point functions of anti-chiral/chiral operators in N = 2 model and N = 4 SYM one finds that, in the planar limit, it is 1 for even k = 2n but a nontrivial function of λ for odd k = 2n + 1. This property can be understood as follows. Considering an orbifold projection of a conformal gauge theory with respect to some discrete subgroup of a global symmetry group one can split the observables into "untwisted" (even under the action of orbifold group) and "twisted" (odd under the orbifold group) [49,[54][55][56]. In the planar limit, the former are the same as in the original theory (but they differ at subleading order in 1/N), while the latter deviate from the original theory ones already at the leading large N order. The free energy, the circular Wilson loop and the ratio of the correlators R M 2n are examples of "untwisted" observables whereas the ratio R M 2n+1 and analogous ratio of the correlators in class II models (see Eq. (4.25) below) belong to the "twisted" sector.
In the Q 2 orbifold model with the SU(N) × SU(N) gauge group the operators in the twisted sector are odd under interchanging the scalar fields ϕ 1 and ϕ 2 of the two N = 2 SU(N) vector multiplets. The SA model is related to the Q 2 model by an additional orientifold projection (see, e.g., [57]). In particular, the scalar fields ϕ 1 and ϕ 2 at the two nodes of the Q 2 quiver are related to the scalar field ϕ belonging to the N = 2 vector multiplet of the SA theory as ϕ 1 = ϕ and ϕ 2 = −ϕ. This additional discrete modding implies that the chiral primary operators O 2n+1 = tr ϕ 2n+1 (x) and O 2n = tr ϕ 2n (x) belong, respectively, to the twisted and untwisted sectors in the SA theory. This explains why the ratio of the correlators (4.17) is different from 1 for the "odd" operators already at the planar level. 15 Explicitly, in the SA model one finds for (4.17) from the localization matrix model representation , (4.18) where the semi-infinite matrix K [n] is obtained from the matrix K in (1.2) with ℓ = 2 by removing its first (n − 1) rows and columns. In [18], the relation in (4.18) was used to derive the weak coupling expansion of R SA 2n+1 . At strong coupling, R SA 2n+1 can be determined in the leading-order (LO) approximation as [58] R SA 2n+1 (λ) = 8π 2 λ n(2n + 1) + . . . .
According to (4.7), the 1/N 2 term here scales at strong coupling as λ 1/2 and is proportional to C SA in (4.7). In an attempt to improve on the LO value of this coefficient (4.9), a refined analysis was performed by either keeping the next to leading order terms in the Bessel operator matrix elements or considering a sequence of its numerically exact finite dimensional truncations. The best estimates obtained by these two approximations were respectively [58] C SA NLO = 0.113 , C SA num = 0.130 . The same methods were also applied to some three-point functions of single-trace chiral primaries in the SA model [59]. The normalized extremal correlators can be expressed in terms of the resolvent 1/(1 − K) nm . For even n 1 and n 2 , the ratio tends, in the planar limit, to 1 (as expected). In other cases the leading term in the strong-coupling expansion is given in the LO approximation by  distribution of eigenvalues in the large N limit. This assumption is justified for correlators of "even" chiral primaries O 2n (x) = tr ϕ 2n (x). However, the correlators of "odd" chiral primaries O 2n+1 (x) involve a deformed eigenvalue distribution [18].
In class II models, i.e. Q L quiver theories, a generalization of (4.17) and (4.18) has been worked out in [57]. The operators in the "twisted" sector of the Z L orbifold are defined as where integer α satisfies 1 ≤ α ≤ L − 1 and O (I) n = tr ϕ n I corresponds to the I-th node of the quiver. The two-point functions T α,n (x)T α,n (0) of the twisted chiral operators (4.24) can be computed using the localization technique. In close analogy with (4.17), one can define the ratio (4.25) In the planar limit, it admits the following representation in terms of the matrix (4.1) for even and odd n (see Eq. (5.21) in [57]) , where the matrices X even and X odd coincide with (4.1) for ℓ = 1 and ℓ = 2, respectively. The matrix X [n] is again obtained from X by removing its first (n − 1) rows and columns. The dependence on α enters (4.26) through the parameter 16 s α = sin 2 πα L . (4.27) At strong coupling, the expressions in (4.26) can be evaluated in the LO approximation to give, for both even and odd n, Recently, three-point functions of similar BPS twisted-sector operators in the orbifold Q L theory were computed in the planar limit at LO approximation in strong-coupling expansion and successfully matched to dual string theory (which in this planar limit of correlators of BPS operators is represented by type IIB supergravity) [60,61].

Comments
Let us draw some conclusions from the above discussion.
1. The Bessel operator enters the N = 2 observables in two distinct ways: the free energy and non-planar correction to circular Wilson loop depend on its Fredholm determinant F loc ℓ = log det(1 − K) (see Eqs.(4.6) and (4.16)). At the same time, the correlators of chiral primary operators are expressed in terms of matrix elements of its resolvent 1/(1 − K) (see Eqs. (4.18) and (4.26)).
2. The results for the strong-coupling expansion of F loc ℓ available so far in the literature were mostly numerical and inconclusive. This applies, in particular, to the value of the coefficient of the leading λ 1/2 term in (4.7) (cf. Eqs. (4.9) and (4.21)), not to mention higher subleading corrections that remained unknown.
3. The matrix elements of 1/(1 − K) are more under control. In particular, the leading O(1/λ) term in the ratios of chiral correlation functions (4.19), (4.23) and (4.28) has been computed numerically, showing a very good agreement with analytic expressions, and, in some cases, was matched to dual supergravity results [60]. Still, subleading corrections at large λ have not been computed yet, as going beyond the LO approximation appears to be non-trivial.
4. Apart from "perturbative" ( 1 √ λ ) n strong-coupling corrections to the observables one expects also exponentially small non-perturbative corrections but no information about them is available except in the case of much simpler FA model [21].

Strong coupling expansion
Let us now apply the results obtained in Sections 2 and 3 to address the open problems mentioned above.

Free energy and circular Wilson loop
As was explained in the previous subsection, these observables can be expressed in terms of the Fredholm determinant of the Bessel operator (2.3) with the symbol (1.6). From the results in Section 2.5 we have .
Here the Widom-Dyson constant B loc ℓ is given by (see (2.24)) where A is Glaisher's constant. Perturbative corrections in 1/g are described by the function f loc ℓ (g + i0), where '+i0' corresponds to a particular prescription for integrating the Borel singularities in (1.11). Its series expansion at large g looks as (see (2.27), (2.11), (2.12) and where g ′ = g − log 2 π . Re-expanding this series in powers of 1/g one would generate terms proportional to powers of log 2 π . Finally, the leading non-perturbative correction to ∆f loc ℓ (g) is given by (3.81) and (3.82).
Substituting (4.29) into the first relation in (4.6) we get the strong-coupling expansion of the free energy in the SA model 17 Let us make a few comments on this result. The perturbative series on the second line of (4.32) has an interesting "homogenous weight" property. Namely, the coefficient in front of 1/λ ′n is proportional to the product of odd Riemann zeta values ζ(2n i + 1) with i n i = n.
Note that both the leading O(λ 1/2 ) term and non-perturbative correction on the last line of (4.32) may also be expressed in terms of λ ′ . 18 Another comment is about the imaginary coefficient of the non-perturbative correction in (4.32). It has the same origin as the coefficient in front of the first term in (3.48). The perturbative part of the strong-coupling expansion in (4.32) is not Borel summable and the Borel singularities should be avoided by slightly moving g in the complex plane. As a consequence, an imaginary part is produced canceling a similar contribution from the non-perturbative correction in (4.32). The latter is dressed by a similar perturbative tail of powers of λ −1/2 . From (4.32) we find, in particular, the exact value of the coefficient C SA defined in (4.7) This demonstrates that (4.9) is a rough estimate and that numerical results (4.21) gives a relatively good approximation.
Similarly, for the orbifold model Q 2 , we obtain from the second relation in (4.6) Notice that the leading non-perturbative correction in (4.34) is suppressed by a factor of λ 1/2 as compared to the one in (4.32), i.e. scales as O(e − √ λ ). This is a consequence of the fact that while both terms in (4.6) receive the leading O( √ λ e − √ λ ) nonperturbative corrections, they cancel against each other in the sum. The normalization coefficient c 1 of the subleading term remains to be determined: to find its value, one has to compute the subleading O(1/g) correction to the functions (3.53).
The leading term of the expansion (4.34) has the expected form (4.7) with the exact value of the coefficient being This may be compared to the previous numerical estimate (4.8). The disagreement is not surprising since the numerical analysis in [19] could only reach λ ≃ 450 where numerical fitting is not yet able to disentangle λ 1/2 from the slowly growing log λ corrections. Interestingly, the Q 2 coefficient (4.35) is twice that of the SA one in (4.33). This may be related to the fact that the SA model is an orientifold projection of the Q 2 model (see Appendix C of [57] for details). In particular, the free-theory content of the SU(N) × SU(N) Q 2 quiver theory is twice that of the SA theory (one free bi-fundamental SU(N) hypermultiplet is the same as the sum of rank-2 symmetric plus antisymmetric hypers). 19 This factor of 2 proportionality is, of course, no longer true for the subleading coefficients in (4.32) and (4.34).
Substituting (4.32) and (4.34) into (4.15) and (4.16) we obtain the strong-coupling expansion of the leading non-planar correction (4.14) to the circular Wilson loop (4.37) 19 In particular, the conformal anomaly coefficients of the two theories thus also differ by factor of 2: Here the first and the second line in each relation defines, respectively, the perturbative and non-perturbative corrections. The equality of the coefficients of the leading λ 3/2 terms in (4.36) and (4.37) follows from (i) the factor of 2 proportionality between (4.33) and (4.35) mentioned above and (ii) the fact that the half-BPS Wilson loop in the Q 2 model is defined in terms of the fields associated with just one of the two SU(N) factors of the gauge group resulting in extra factor of 1/2 in (4.16) as compared to (4.15).

Two-point chiral correlators in SA and Q L models
We can also apply (4.29) to derive the strong coupling expansion of the two-point correlation function of (anti) chiral operators in (4.17).
To this aim, in the SA model, we write (4.18) in the form where we first applied the Cramer's rule to (4.18) and then replaced the determinants by their expressions in terms of the functions F ℓ in (1.1).
At weak coupling, we can use (2.6) and (2.8) to expand R SA 2n+1 (λ) in powers of 't Hooft coupling. We checked that the resulting expressions are in agreement with the results of [18] (see Eq. (6.15) there). At strong coupling, we find from (4.38) and (4.29) where (λ ′ ) 1/2 = λ 1/2 − 4 log 2. Here the expression on the second line yields perturbative corrections in λ −1/2 . As above, using λ ′−1/2 as an expansion parameter is advantageous as it automatically sums up all terms proportional to log 2. The last line in (4.39) represents the leading non-perturbative correction. It comes from the difference ∆f loc 2n+2 − ∆f loc 2n in the exponent of (4.38) and is given by the second relation in (3.82) for ℓ = 2n.
Notice that the strong-coupling expansion (4.39) involves half-integer powers of λ −1 , in agreement with the expectation from the dual string theory side (where λ −1/2 is the inverse string tension). Such terms are not present in the large λ asymptotic expansion of individual matrix elements K ij in (4.38). This illustrates (once again) non-trivial properties of the Fredholm determinant of the Bessel operator as well as the power and efficiency of the techniques for computing it described in the first part of this paper.
The expression on the first line of (4.39) defines the leading behaviour of R SA 2n+1 for λ → ∞. It agrees with (4.19), i.e. confirms the validity of the LO approximation used in [58]. The subleading corrections in (4.39) may be computed to any desired order.
We remark that the relation (4.39) has an interesting behaviour in the double scaling limit λ → ∞ and n → ∞ with the ration = 4πn λ −1/2 held fixed. It is equivalent to the so-called "large bridge" limit previously studied in [10,12] in application to the octagon. From (4.39) we get in this limit For the correlators of even chiral operators in the SA model, we apply (4.20) to find the ratio R SA 2n in terms of the non-planar correction in the free energy ∆F SA (λ) defined in (4.6). The same method can be used to compute the two-point correlation functions (4.25) of twisted-sector operators in the Q L orbifold model. To start with, we notice that semi-infinite matrices in (4.26) coincide with the matrix (1.2) evaluated for special values of ℓ and the symbol replaced with χ(x) = s α χ loc (x), schematically, where χ oct (x) and s α are defined in (1.6) and (4.27), respectively. In close analogy with (4.18) and (4.38), each quantity in (4.26) can be expressed as a ratio of the determinants (1.1) and (2.3) where F loc ℓ (g, α) is a logarithm of the Fredholm determinant of the Bessel operator with the symbol χ = s α χ loc (x). The function F loc ℓ (g, α) vanishes for s α = 0 and coincides with (4.29) for s α = 1.
At strong coupling, we apply (1.7)-(1.10) and (2.25) to obtain the following expansion of the function F loc ℓ (g, α) log (4s α ) + log Γ(ℓ) + B(α) where g ′ ≡ g − 1 2 I 1 (s α ). Here we replaced the expansion coefficients with (2.25) and denoted by I n (s α ) the integrals (2.9) evaluated for χ = s α χ loc (x). Note that I n (s α ) is actually a function of α/L since this is the ratio that appears in s α in (4.27). Also, in the Q 2 model the only twisted sector has α = 1 and thus s α = 1 in (4.27). Hence, in this case, the integrals I n are given by (2.11). For generic values of α/L, the integrals I n (s α ) can be expressed in terms of derivatives of the digamma function The coefficient (2.17)  It does affect, however, the normalization factor b = 4s α in (2.15). The last three terms on the first line of (4.45) come from the Widom-Dyson constant (2.22). The first two terms on the right-hand side of (2.22) do not depend on ℓ and are denoted by B(α). For the sake of simplicity, we did not display non-perturbative corrections to (4.45). They are given by the general expressions (3.68) and (3.69) with x = ±2iπx 1 defined as the smallest root of 1 − s α χ loc (x). For s α = 1, the relation (4.45) coincides with (4.29). Substituting (4.45) into (4.43) we obtain the strong-coupling expansion (4.48) where (λ ′ ) 1/2 = λ 1/2 − 2πI 1 (s α ) and I n (s α ) are given by (4.46). Here the first line defines the leading behaviour at strong coupling and is in agreement with Eq. (5.32) in [57]. Expanding (4.48) in powers of 1/g = 4π/λ 1/2 would produce terms proportional to I 1 (s α ), e.g., As was mentioned above, in the Q 2 theory where s α = 1 the values of I n are given by (2.11) and (2.12). We observe that in this case 1 + ∆ α,k in (4.48) with k = 2n + 1 becomes the same as the 2-point function (4.39) in the SA theory. This equality holds to all orders as follows from the second relation in (4.26), (4.18) and the first relation in (4.42). In general, this is also a consequence of the fact that the SA model is obtained from the Q 2 one by an extra projection that effectively implies relation between these correlators (see discussion below (4.17)).

Conclusions
In 2) after choosing particular values of the non-negative integer ℓ and the "symbol" function χ(x) (see, e.g., (4.5) and (4.6)). 21 The Fredholm determinant representation of the observables is exact in the 't Hooft coupling λ. While their small λ expansion is straightforward, it is quite non-trivial to develop a strong-coupling expansion. The latter is of prime interest as it should be equivalent to the inverse string tension expansion according to the AdS/CFT duality. The advantage of the Fredholm determinant representation is that this difficult problem can be solved by applying the strong Szegő limit theorem. It requires an application of special techniques partially available in mathematical literature.
We found that, for the physically relevant cases, the perturbative expansion of F ℓ (g) in powers of 1/g = 4π/ √ λ reveals a number of remarkable properties that were previously observed in the case of the octagon correlator in planar N = 4 SYM in [9,11,12]. In particular, we demonstrated that F ℓ (g) receives a log g correction. It originates from the Fisher-Hartwig singularity of the symbol χ of the Bessel operator and is given by a simple expression (βℓ + 1 2 β 2 ) log g, which only depends on ℓ and the strength of the singularity β (cf. (1.8)). The structural simplicity of this term calls for its interpretation on the dual string theory side.
Examining the expansion of F ℓ (g) in powers of 1/g, we found that the resulting series can be significantly simplified by changing the expansion parameter to g ′ = g − 1 2 I 1 where a transcendental constant I 1 is given by (2.11) and (2.12). This effectively performs a resummation of all terms containing powers of I 1 . A similar phenomenon was previously noticed in the strong-coupling expansion of the cusp anomalous dimension in planar N = 4 SYM [39] and, thus, may be a universal feature of strong-coupling expansion of observables in superconformal theories.
The resulting perturbative expansion of the determinant in (1.1) takes the following factorized form (see Eqs. (1.7), (2.9), (2.17) and (2.27)) where ℓ β = ℓ + β. Here the two factors depend separately on g and g ′ . Notice that, due to our choice of the g ′−1 = (g − 1 2 I 1 ) −1 as the expansion parameter, the first few powers of g ′−1 are absent inside the second brackets. The power of g in the first factor in (5.1) does not depend on the symbol of the Bessel operator, whereas that of g ′ in the second factor only depends on the strength of the Fisher-Hartwig singularity β of its symbol. It would be interesting to understand the meaning of the g vs g ′ factorization in (5.1) from the dual string theory perspective.
Let us note that in the simplest ℓ = 0 octagon case we have at strong coupling O 0 = exp [F oct 0 (g)] ∼ g 1/8 g ′3/8 e −gA 0 ∼ g 1/2 e −gA 0 , where A 0 = 2I oct 0 = π [10,11]. This implies that the corresponding planar correlator G 0 of four BPS operators scales as 2π A 0 . The constant A 0 was conjectured to have a dual semiclassical string theory interpretation as minimal area of a world sheet surface in AdS 5 ×S 5 [10]. The prefactor √ λ of the exponential may be given the following heuristic string theory explanation. The dual string theory representation for G 0 is in terms of a correlator of 4 BPS vertex operators on a plane or S 2 . Each of them comes essentially from string action and thus carries a normalization factor proportional to the string tension T = √ λ 2π . There is also a string tension dependence in the Möbius volume for S 2 (including which is like dividing by a 3point function); that gives a factor of T −3 . In total, we get T 4 × T −3 ∼ √ λ, in agreement with the above strong-coupling scaling of G 0 .
In this paper we extended the analysis of the Fredholm determinant of the Bessel operator to a wider class of symbols that appear in special N = 2 superconformal models that are planar equivalent to N = 4 SYM. We derived the strong-coupling expansion of the relevant observables and resolved some issues that existed in earlier work [19,20,58]. In particular, we obtained the strong-coupling expansion of the free energy (4.32) in the SA theory (N = 2 model with matter in rank-two symmetric and antisymmetric representations of SU(N)) determining the exact value (4.33) of the coefficient of the leading O( √ λ) term. We also computed systematically the strong coupling corrections to certain two-point functions of chiral primary operators, confirming the conjecture for the coefficient of the leading term [58].
Apart from analytically deriving the coefficients of the perturbative strong-coupling expansion, one main result of this paper is the identification of the non-perturbative (exponentially small at large λ) corrections (1.12) to the Fredholm determinant exp(F ℓ (λ)) and, thus, to the related observables in N = 4 and N = 2 theories. Such corrections are known to be present, in particular, in the cusp anomalous dimension in planar N = 4 SYM where they drive the transition from strong to weak coupling [39,40]. In that case, the leading non-perturbative correction to the cusp anomalous dimension has a clear physical meaning on the dual string theory side [62]: it coincides with the dynamical mass gap of the effective two-dimensional bosonic O(6) sigma model describing massless excitations of the Gubser-Klebanov-Polyakov string. We showed that for the Bessel operator with a general symbol χ(x), the leading non-perturbative term (1.12) depends on the smallest root of 1 − χ(2iπx 1 ) = 0 and its multiplicity n 1 . It would be interesting to understand the dual string theory origin of this correction in the N = 2 superconformal models and whether it also admits some dynamical mass scale interpretation. 22 In general, the wealth of new results for the strong-coupling expansion of various observables obtained in this paper calls for a detailed comparison with dual string theory. The required computations on the string theory side are, however, appear to be very nontrivial. For example, the leading non-planar corrections to free energy in SA (4.32) or Q 2 (4.34) theories should be reproduced by the torus correction in type IIB superstring theory on the corresponding orientifold/orbifold of AdS 5 × S 5 . It is unclear if even the leading √ λ term (4.32) and (4.34) can be reproduced from the one-loop string effective action . . or one needs to compute the torus partition function exactly before expanding in α ′ ∼ 1 √ λ (for some related discussion see [20,21]). 23 National Agency for Research grant ANR-17-CE31-0001-01. AAT was supported by the STFC grant ST/T000791/1.

Note added in v2:
The first few terms of the strong coupling expansion of the observables in the SA theory presented in Section 4.2 were reproduced in a very recent paper [65] by using a high precision numerical calculation.

A Matrix model representation
The localization technique allows us to evaluate the partition function of N = 2 superconformal models on the 4-sphere as a matrix integral [2]. For the models defined in Table 1, this integral (normalized to N = 4 SYM expression) can be expressed in the large N limit as a Fredholm determinant of the Bessel operator with the symbol of the form (1.6). Below we review a derivation of this relation [20] and, then, construct the matrix integral representation of the observable (1.1) and (2.3) for an arbitrary symbol χ(x). The partition function of N = 2 theory with the SU(N) gauge group on S 4 is given by , where integration goes over Hermitian traceless matrices a of dimension N. Here we neglected the instanton contribution because we are interested in the large N limit. The interaction potential is given by (see also [32]) where the function H(x) can be expressed in terms of the Barnes function The two terms in the first relation in (A.2) involve traces over the matter representation R and the adjoint representation of the SU(N) group, respectively. Their difference is denoted as a trace over R ′ . In N = 4 theory, R is the adjoint representation and the interaction term (A.2) vanishes. As a consequence, the partition function Z N =4 is given by a Gaussian integral whose evaluation gives the free energy (4.3), Z N =4 = exp(−F N =4 ). In N = 2 superconformal models of type SA (see Table 1), the trace (A.2) can be evaluated in the large N limit using the identity tr R ′ a 2k = Note that the sum involves traces of odd powers of matrix a. In particular, tr R ′ a 2 = 0 and, as a consequence, the O(x 2 ) term on the right-hand side of (A.3) does not contribute to the partition function (A.1). Neglecting this term and introducing new variables ω k (a) (with k ≥ 1) that satisfy tr(a 2n+1 ) = g 2n+1 the interaction term (A.2) can be rewritten in the planar limit as Here a semi-infinite matrix X nm is given by where g 2 = λ/(4π) 2 = g 2 YM N/(4π) 2 and the function χ loc (x) is defined in (1.6). The equivalence of the two representations of the interaction term, Eqs. (A.2) and (A.6), relies on the relation satisfied by the function χ loc (x) (A. 8) In the absence of the interaction term in (A.1), i.e. S int (a) = 0, the ω n -variables have diagonal two-point expectation values with respect to a Gaussian integration measure: ω n (a)ω m (a) = δ nm in the large N limit. Adding the interaction term (A.6), one obtains the following representation of the partition function (A.1) in the SA model Comparing (A.7) with the analogous matrix defined in (1.2) we observe that they coincide for ℓ = 2 and χ = χ loc Combining together the above relations with (1.1) and (2.3), we conclude that .
where the sum of the three terms in the parenthesis scales as O(t 4 ) at small t.

B Bessel kernel
In this Appendix, we establish the relation between the semi-infinite matrix (1.2) and the integral Bessel operator defined in (2.1) and (2.2). We show below that the matrix K nm represents the Bessel operator B ℓ on a space spanned by basis functions Here the kernel B ℓ (t 1 , t 2 ) is defined in (2.2) and the functions ψ m (t) (with m ≥ 1) have the following form where χ(x) is the symbol of the Bessel operator.
To prove (B.1) it proves convenient to introduce an auxiliary function Here the second relation is known as the Bessel kernel, see e.g. [24]. Defining an integral operator K ℓ with the kernel K ℓ (x, y), it is straightforward to verify using the first relation in (B.3) that where the functions φ n (x) are given by We can use the well-known orthogonality property of the Bessel functions to check that the functions φ n (x) are orthogonal with respect to the scalar product The integral operator K ℓ has previously appeared in the study of the octagon in planar N = 4 SYM [11,12]. Going back to (B.1), we replace the kernel B ℓ (t 1 , t 2 ) with its expression (2.2) to get where the notation was introduced for Here in the first relation ψ m (t 2 ) was replaced with its expression (B.2) and in the second relation the integral over t 2 was evaluated using the last relation in (B.3). Changing variable y → y/(2g) in (B.8) and taking into account (B.4), we arrive at Substituting this expression into (B.7), we observe that the integral over x is proportional to ψ n (t 1 ) with the proportionality coefficient being K nm . In this way, we reproduce (B.1).

C General solution for half-integer β
In this Appendix, we present some details of the derivation of the relations (3.12) and (3.13) for the function Γ(x, y). For the sake of simplicity we concentrate on the case of even ℓ, for odd ℓ analysis goes along the same lines.
For even ℓ, we use the relation Γ(−x, y) = Γ(x, y) to simplify (3.11) as whereΓ(k, y) is an even function of k. Substituting this relation into (3.10) and taking into account the identity where U k (x) is Chebyshev polynomial of the second kind, we arrive at an infinite system of equations for the functionΓ(k, y) Taking advantage of the fact that U m (k) = (−1) k U m (−k) are orthogonal polynomials we can write a general solution to (C.3) as where a n (y) are arbitrary coefficient functions. SinceΓ(k, y) is even in k, the sum involves the Chebyshev polynomials with even indices. The relation (C.5) only holds for k 2 < 1.
To find the functionΓ(k, y) for k 2 > 1, we invert (3.11) and replace Γ(x, y) with its Here the function γ(x, y) is given by a double Neumann series over the Bessel functions (3.3). For each term in this series, e −ikx γ(x, y) vanishes as x → ∞ for k 2 > 1. For k > 1 and k < −1 this allows us to deform the integration contour in (C.6) to the lower and upper half-plane, respectively, and pick up the residues at the poles of the function 1 − χ(x/(2g)) defined in (2.13) and (2.18). The poles are located at x = ±4πigy n and their contribution to (C.6) takes the form where c n (y) e 4πgyn (with n ≥ 1) are given by the residue of (C.6) at the poles.
To find the function Γ(x, y) we split the integral in (C.1) into k 2 < 1 and k 2 > 1 and replaceΓ(k, y) in each of the regions by the corresponding expression, Eqs. (C.5) and (C.7), respectively. Integrating (C.5) over −1 < k < 1, we encounter the integrals where p ℓ (x) is a polynomial in x of degree ℓ defined in (3.14). Then, the two terms on the right-hand side of (C.5) give rise to the first two lines of (3.12). Finally, integrating (C.7) over k 2 > 1 we obtain the last line of (3.12). Repeating the above calculation for odd ℓ, one can derive (3.13).
For ℓ ≥ 2 the relations (3.12) and (3.13) also involve ⌊ℓ/2⌋ functions a n (y). They can be expressed in terms of the functions c j (y) by imposing the condition (3.6). For instance, for β = −1/2 and ℓ = 2, the relation (3.12) combined with Γ(2gx, 2gy) ∼ x 2 as x → 0 leads to a 0 (y) = − sin(y) y − 1 4g j≥1 c j (y) y j . (D.11) Replacing the functions a n (y) in (3.12) and (3.13) with their expressions and imposing the relation (3.8), we obtain an infinite system of equations for the functions c j (y). As in the previous case, its solution takes the form (D.10).

E Borel singularities
The strong-coupling expansion (1.7) involves perturbative Borel non-summable series (1.10) and non-perturbative exponentially small term, Eq. (1.12). To define unambiguously the corresponding functions f (g) and ∆f (g), we have to specify a prescription for integrating the Borel singularities of (1.10).
To illustrate the procedure, here we discuss the strong-coupling expansion (3.48). In close analogy with (3.46), we define the ratio of the modified Bessel functions r(x) = I 0 (x)/I 1 (x) . (E.1) The expression on the left-hand side of (3.48) involves r n = r(4πgx n ).
Recall that for integer index n, I n (x) are entire functions of x with I n (−x) = (−1) n I n (x) , (E.2) that lead to r(−x) = −r(x). At large x, they have the following asymptotic expansion where a k = 1 2 − n k 1 2 + n k /(2 k k! ). Note that for x > 0 the second, exponentially small term on the right-hand side of (E.3) is pure imaginary. It is this term that gives rise to non-perturbative correction in (3.48). Its appearance is closely related to the fact that the series which accompanies the first term in (E.3) is not Borel summable. The problem here is similar to that of separation of the perturbative f (g) and non-perturbative ∆f (g) terms in (1.7) mentioned above.
To define the two terms on the right-hand side of (E.3), we apply the identities between the modified Bessel functions of the first and second kind In contrast to I n (z), the function K n (z) has a cut on the complex z-plane running along the negative real axis. The argument of the Bessel functions on the right-hand sides of (E.4) is shifted slightly away from the real axis to avoid the cut. The function K n (z) admits an integral representation where in the second relation we integrated by parts. At large z, we can expand the hypergeometric function in powers of 1/z and integrate term-by-term to obtain a power series representation of K n (z). In this way, we verify that the second term on the right-hand side of (E.3) arises from K n (x − i0) in (E.4). The first term in (E.3) comes from the large x expansion of K n (−x + i0) in (E.4). The Borel singularities arise in this term because the function K n (−x) has a discontinuity across the cut x > 0. The simplest way to see this is to notice that the integral in (E.5) takes the form of the Borel transform (1.11). The hypergeometric function in (E.5) has a pole at σ = −2z 2 F 1 3 2 − n, 3 2 + n; 2; − σ 2z = 2z (σ + 2z) Γ 3 2 − n Γ 3 2 + n + . . .

(E.6)
For z < 0 this pole is located on the integration contour in (E.5) making K n (z) ill defined. Replacing z → z = −x ± i0 in (E.5) amounts to deforming the integration contour in the vicinity of the pole as σ → σ ± i0, schematically, At large positive x, we can identify the first term inside the brackets on both lines of (E.4) as defining a perturbative correction to I n (x). Notice that this correction is different for the first and the second relation in (E.4). The difference is proportional to the sum of the functions K n (−x + i0) + K n (−x − i0). It is given by the residue at the pole σ = 2x and yields an exponentially small correction ∼ e −x / √ x that matches the second term in (E.3).
Due to the equivalence of the two representations (E.4), it is compensated by the analogous exponentially small correction coming from K n (x) = K n (x − i0) = K n (x + i0).
Taking an average of the two relations in (E.4), we can get another equivalent representation which is valid for x > 0. It corresponds to the principal value prescription for integrating the Borel pole (E.7). We conclude from above analysis that the relations (E.4) and (E.8) correspond to three different prescriptions of regularizing Borel singularities in the asymptotic expansion (E.3).
Let us return to (E.1) and apply the first relation in (E.4) where in the last equality we used the identity for the Bessel function. Applying (E.5) we find that for large positive x the first term in (E.9) defines a perturbative correction to r(x). The second term in (E.9) is exponentially suppressed , (E.11) where in the second line we replaced the Bessel functions by their leading large x expansion. The series inside the brackets suffer from Borel singularities. As explained above, the prescription −x + i0 in the argument of the Bessel functions is equivalent to deforming the integration contour in their Borel transform slightly above the real axis. Had we applied the second relation in (E.4) we would get a similar relation for r(x) with the sign in front of the second term in (E.11) reversed.
Applying (E.11) we can obtain the strong-coupling expansion of r n = r(4πgx n ) and, then, use it to reproduce the relation (3.48).

F Free energy in Q 2 model
The weak-coupling expansion of non-planar correction (4.16) to the circular Wilson loop ∆q Q 2 looks as ∆q Q 2 (λ) = λ To find ∆F Q 2 (λ) at arbitrary coupling, we use the determinant representation of the partition function in the quiver model Q L derived in [57,66]. Specialising the results of these papers to L = 2, we get where the matrices X even and X odd coincide with (4.1) for ℓ = 1 and ℓ = 2, respectively, The product (1 − X even )(1 − X odd ) arises in (F.4) due to the even-odd block factorization of the quadratic form associated with the underlying matrix model at large N, (cf. Eq. (5.9) of [57]). The second relation in (F.4) follows from (1.1) by noticing that X even and X odd coincide with the matrix (1.2) evaluated for ℓ = 1 and ℓ = 2, respectively, and for the symbol given by (1.6). Notice that 1 2 F loc ℓ=2 in (F.4) gives the free energy (4.5) in the SA model. This is a consequence of the fact that the SA model may be viewed as an orientifold projection of the L = 2 quiver [57]. After the projection, the even-even double trace terms in the matrix model interaction potential (A.2) proportional to X even are removed and the matrix X reduces to X odd .