Fishnet four-point integrals: integrable representations and thermodynamic limits

We consider four-point integrals arising in the planar limit of the conformal"fishnet"theory in four dimensions. They define a two-parameter family of higher-loop Feynman integrals, which extend the series of ladder integrals and were argued, based on integrability and analyticity, to admit matrix-model-like integral and determinantal representations. In this paper, we prove the equivalence of all these representations using exact summation and integration techniques. We then analyze the large-order behaviour, corresponding to the thermodynamic limit of a large fishnet graph. The saddle-point equations are found to match known two-cut singular equations arising in matrix models, enabling us to obtain a concise parametric expression for the free-energy density in terms of complete elliptic integrals. Interestingly, the latter depends non-trivially on the fishnet aspect ratio and differs from a scaling formula due to Zamolodchikov for large periodic fishnets, suggesting a strong sensitivity to the boundary conditions. We also find an intriguing connection between the saddle-point equation and the equation describing the Frolov-Tseytlin spinning string in $AdS_{3}\times S^{1}$, in a generalized scaling combining the thermodynamic and short-distance limits.


Introduction
Over the past few years, we have witnessed tremendous progress in our understanding of mathematical and physical properties of Feynman integrals, notably in the massless limit relevant for the description of high-energy processes. A great deal of this progress has centred around the most symmetrical and solvable conformal field theory in four dimensions, planar N = 4 super-Yang-Mills (SYM) theory, where dualities, conformal symmetry and integrability paved the way to the development of a variety of new methods for efficient higher-loop calculations. It has also been realized that these methods, including integrability, are not restricted to supersymmetric planar field theories. They extend to a large class of interacting non-supersymmetric planar theories, the fishnet conformal field theories [1][2][3][4][5][6][7][8][9], which admit formulations in any spacetime dimensions and may be realized in four dimensions by deforming N = 4 SYM, in ways that preserve conformal symmetry and integrability. This web of theories allows us to shed light on properties of a broad family of higher-loop planar integrals, with regular iterative structures, the so-called fishnet graphs. These diagrams are endowed with integrable structures, such as Yangian symmetry [10][11][12][13], which put severe constraints on their analytic expressions and facilitate their determination.
In this paper, we will examine a simple class of planar fishnet integrals in four dimensions, associated with regular square lattices attached to four generic spacetime points, as shown in the left panel of figure 1. They generalize the well-known ladder conformal integrals, calculated many years ago [14]. They are subject to stringent analytic properties and were argued to admit integrability-based integral and determinantal representations [15], some of which have been proven recently [16][17][18]. The first goal of this paper will be to place the proposals in ref. [15] on a firmer footing by proving their mathematical equivalence. To be precise, we will show the equivalence of the so-called Berenstein-Maldacena-Nastase (BMN) [19] and flux-tube (FT) integral representations with the determinant of a Hankel matrix of ladder integrals, using exact summation and integration techniques.
The integrability of fishnet graphs was noticed many years ago by Zamolodchikov [20], who considered them as an example of exactly solvable lattice models and calculated their leading behaviour in the thermodynamic limit for periodic boundary conditions. Taking this "continuum" limit was partly motivated by viewing fishnet graphs as approximating the worldsheet of a string, see e.g. ref. [21]. Recent studies have given support to this picture, in a more holographic setting, connecting the conformal fishnet integrals to quantum mechanical systems in Anti-de Sitter space (AdS) [22][23][24][25][26]. In particular, in the continuum limit, evidence was found that tubular fishnet diagrams, like the one shown in the right panel of figure 1, admit an effective description in terms of a two-dimensional (2d) non-linear sigma model in AdS.
Motivated by this consideration, in the second half of this paper, we study the thermodynamic limit of the four-point fishnet integrals. We do so by taking the dimensions of the four operators, or alternatively the lengths of each side of a rectangular fishnet, to be large, holding the aspect ratio (the ratio of side lengths) fixed. In this limit, the integrability-based integrals can be evaluated in a saddlepoint approximation, and the saddle-point equations coincide with solvable problems encountered in the study of matrix models and integrable systems. This connection will allow us to obtain a closed parametric expression for the bulk free energy in terms of complete elliptic integrals. Interestingly, the free energy is independent of the cross ratios parametrizing the positions of the boundary operators. However, we will find that it depends in a complicated manner on the rectangular aspect ratio. For any value of the aspect ratio, it differs from Zamolodchikov's result for periodic graphs, indicating a strong dependence of the thermodynamic limit on the boundary conditions. We shall also study a more general scaling regime, which combines the short-distance and thermodynamic limits; this limit 1 2 3 4 Figure 1. Specimens of planar fishnet graphs in position space with different boundary conditions. Bulk points stand for quartic integration points, d 4 x/π 2 , while edges are massless propagators, 1/(x − y) 2 . For the fishnet four-point functions considered in this paper, the m vertical and n horizontal external legs end on four fixed spacetime points x1,2,3,4, as shown in the left panel (m lines end on x3,4; n lines end on x1,2). In the right panel, we show an m-by-n fishnet graph with periodic boundary conditions in the angular direction. Such graphs are expected to scale as ∼ g −2mn c when m, n → ∞, with gc a numerical constant [20]. The four-point function considered in this paper (left panel) scales differently, with a free energy per vertex that depends on the aspect ratio k = n/m. reveals an intriguing connection with the equation describing a classical spinning string in AdS 3 × S 1 , studied by Frolov and Tseytlin [27].
This paper is structured as follows. In section 2, we analyze the integrability-based representations for the fishnet integrals and prove their equivalence with the determinant of ladder integrals. Our analysis relies on exact methods which may be applicable to more general fishnet graphs. We proceed with the analysis of the thermodynamic limit in section 3. Applying the saddle-point method to the matrix-model integrals gives us first order differential equations for the free-energy density, which can be solved simultaneously in terms of complete elliptic functions. We compare the obtained expression with a direct numerical evaluation of the determinant and study particular limits analytically. We conclude with a brief discussion of the dependence of the free energy on the boundary conditions. In section 4 we consider a more general scaling which combines the thermodynamic and short-distance limits and unveils a connection with the singular equation describing a classical spinning string in AdS 3 × S 1 . Section 5 contains concluding remarks.

Integral representations
The fabric of the fishnet graphs -the conformal fishnet theory -is a theory of two N c × N c matrices of complex scalar fields φ 1,2 and a quartic interaction, with (Euclidean) Lagrangian density [1,2,5] with the trace running over the "colour" indices and with g 2 a marginal coupling constant, which is kept fixed in the large N c limit [5,28]. 1 x 3 x 4 x 1 x 2 u m … Figure 2. A transformation maps the space-time integral into a momentum space amplitude, shown in red lines, with external massive momenta p1 = x23, p2 = x31, p3 = x14 and p4 = x42, and with massless internal lines. The latter amplitude is subject to important analyticity requirements such as the Steinmann relations.
The fishnet four-point function of interest is defined as the leading colour-ordered contribution to the correlator with the trace embracing all the fields. It receives a single contribution at large N c , which is given by the mn-loop Feynman diagram shown in figure 2. The graph is both ultraviolet (UV) and infrared (IR) finite, for generic x 1,2,3,4 . Dropping a colour factor and stripping off overall weights, one has where Φ m,n is a coupling-independent function of the two conformal cross ratios parametrized here in terms of two numbers z,z. The latter are complex conjugates in Euclidean spacetime signature, and are treated as independent real numbers in the Minkowskian case.
A representation for Φ m,n was proposed in ref. [15] using integrability and analyticity. The proposal is that for n m where I m,n (z,z) is a pure function of weight m × n given by the determinant of an m × m Hankel matrix of ladder functions. It reads and normalization factor The particular case m = 1 coincides with the well-known ladder integral L p (z,z), defined for a number of rungs p 0 by [14,29] L p (z,z) = with Li j (z) the polylogarithm of weight j, where the Taylor series converges for |z| < 1.
The determinant formula (2.6), when interpreted as the result for a scattering amplitude rather than a correlation function, nicely satisfies the Steinmann relations [30,31], which forbid sequential discontinuities in overlapping channels. In figure 2, the red lines indicate the scattering interpretation, and the two channels are (p 1 + p 2 ) 2 = x 2 12 and (p 2 + p 3 ) 2 = x 2 34 . The Steinmann relations imply that Disc x 2 12 Disc x 2 34 I m,n = 0, (2.11) in the region where x 2 12 and x 2 34 both vanish, which is the neighborhood of z,z → 1. Conformal invariance then translates eq. (2.11) into the vanishing of the double discontinuity in (1 − z)(1 −z): (2.12) (For example, one can wrap z around 1, and keepz fixed, or wrap both in the same direction.) The single discontinuity of the ladder integral (2.9) is [15,32] Since the single discontinuity contains only log z and logz, the double discontinuity of any ladder integral vanishes, [Disc 1 ] 2 L p = 0. Furthermore, the discontinuities of any two columns (or two rows) of the matrix M p defined in (2.7) are proportional, because which guarantees the vanishing of the double discontinuity of the determinant (2.6) [15]. (A more general solution to the Steinmann relations is given by the minors of a semi-infinite matrix of ladders, as found in the study of large-charge correlators in N = 4 SYM [33,34].) As alluded to before, along with the determinant (2.6), two matrix-model-like integrals, the BMN integral and FT integral, were given in ref. [15]. They originate from two different conjectural N = 4 SYM integrability-based tiling constructions, the "hexagonalization" [35][36][37][38][39] and the "pentagon OPE" [40][41][42]. The former has been receiving a lot of attention recently, as it is closely related to the method of separation of variables for conformal integrable spin chains, which enabled the calculation of fishnet correlators in the 2d fishnet conformal field theories [7]. This method was recently extended to four dimensions in refs. [16,17,26,43], allowing to put on firm ground the results obtained from the N = 4 SYM integrability-based tiling constructions. According to these analyses, the BMN representation discussed below is no longer a conjecture -it is proved to be equivalent to the Feynman integral. In this section we will close the remaining gaps to the determinant formula by proving its equivalence with both the BMN and FT integral representations, as schematized in figure 3.

BMN Integral
FT Integral Feynman Integral Figure 3. Representations of the fishnet four-point integral. The Feynman integral was recently proved to be equivalent to the BMN integral (2.26) in refs. [16,17]. In this paper we concentrate on the pyramid at the top and derive the determinant formula from the BMN and FT integrability-based representations.

Dual integral
To begin, let us recast the determinant formula (2.5) into an integral form which will prove very useful later on. Introducing the parametrization we will show in this subsection that one may write the determinant in the form where I Dual m,n stands for the m-fold integral 17) with N given in eq. (2.8), and with the kinematical prefactor Note that this prefactor differs from the one in eq. (2.5), such that The representation (2.16) may be seen as a reduction of the more general functional determinant formulae obtained in refs. [44][45][46][47][48][49] for large-charge correlators in N = 4 SYM theory, which resum infinite series of fishnet diagrams and generalizations thereof [33,34]. In the following, we shall refer to eq. (2.17) as the "dual integral", as it relates in our set-up to Fourier-like transforms of the integrability-based integrals to be studied later on.
The dual representation given in eqs. (2.16) and (2.17) holds for σ ∈ R and Im ϕ ∈ (−π, π), which includes both the Minkowskian and the Euclidean kinematics, corresponding respectively to purely real and purely imaginary values of ϕ. It can be derived from the determinant (2.6) by first replacing the ladder functions by their integral representation [29,44,45] . (2.20) This identity follows directly from the definition (2.10) of the polylogarithms. Namely, plugging eq. (2.10) into the sum in eq. (2.9), with 2σ = log (zz), and relabelling j → 2p − j, one finds , (2.21) making use of the combinatorial identity (for p > 0) and changing variables, x = r − σ. One recovers eq. (2.20) by manipulating the contour of integration, noting that the integrand is antisymmetric under x → −x.
One then inserts eq. (2.20) inside the determinant in eq. (2.6) and applies the so-called Cauchy-Binet-Andréief formula (Lemma 1 in appendix A) to bring the determinant below the integral sign. Using eq. (2.19) to remove the factors of (z −z)/ √ zz, one finds (2.23) One encounters two copies of the determinant which is actually independent of σ, and is the Vandermonde determinant

BMN integral
The first integrability-based integral to be discussed is the BMN integral in ref. [15]. This representation yields the pure function I m,n in the form of a sum of integrals, The translation to the variables σ and ϕ follows from eq. (2.15), after fixing the square-root ambiguity appropriately. The representation (2.26) is initially defined in Euclidean kinematics, that is, for σ ∈ R and −iϕ ∈ (−π, π), and then analytically continued. Its direct evaluation is straightforward for low values of m and the result is seen to reproduce the ladder determinant. In particular, one easily verifies the agreement with the ladder series [35], when m = 1. In this subsection, we shall prove the general identity with I Dual m,n the integral form of the determinant (2.17), for all m and n, with n m. To begin, we note that the bulk of the integrand in eq. (2.26) can be written concisely as a Vandermonde determinant. Namely, defining conjugate variables one immediately finds that the integrand can be cast into the compact form where ∆ 2m (ξ) is the Vandermonde determinant for the set {ξ 1 , ξ 2 , . . . , ξ 2m }, (2.32) One can make immediate use of this remarkable simplification to disentangle pairs of ξ variables and cast the integral in the form of the Pfaffian of a 2m × 2m skew-symmetric matrix B, with the sum running over the permutations of the set {1, . . . , 2m}, by plugging into eq. (2.32), and defining A drawback is that the elements B ij (z,z) are not ladders, for j = i ± 1, but derivatives thereof. 3 Further non-trivial algebraic identities are needed to "purify" the expression and reduce the integral to a determinant of an m × m matrix, as discussed in refs. [44][45][46][47][48][49] in the related context of large-charge correlators in N = 4 SYM.
Interestingly, it turns out to be possible to bypass this difficulty by decoupling all ξ variables from the onset. (This promotes the full permutation symmetry of the integrand to the integral level, if not for the contours, see below.) The preliminary step is to relax the constraints i(ξ 2 − ξ 2 −1 ) = a ∈ N among conjugate variables (2.30), by turning the sums over positive a's into integrals. This task may be done efficiently with the help of the Mellin-Barnes summation technique. The key formula is for any suitably smooth function f and with C parallel to the imaginary axis, C = ε+iR with ε ∈ (0, 1). This representation stems from the combination of the Sommerfeld-Watson transform, with an integral transform for the inverse sine factor, see eq. (A.4) in appendix A. This combination is well-suited to our problem as it puts in place the dual variable x entering the integral form of the ladder determinant. It also removes the constraints on ϕ, which may now be chosen anywhere inside the extended kinematics Im ϕ ∈ (−π, π). Applying eq. (2.37) to each sum in eq. (2.31), we get with all dependence on σ factorizing neatly into the integral (2.40) 3 One finds, for r 1, The key simplification is that ξ 2 −1 and ξ 2 can now be regarded as independent variables. They would be real for u ∈ R and a ∈ iR. Instead, they run here parallel to the real axis, along the contours because a has a small positive real part ε. Hence, after taking into account the Jacobian (= −i) for the change of variables (2.30). We may now write I 1 (σ, {x}) as a bunch of decoupled ξ-integrals using eq. (2.34) and recast it as the determinant of a 2m × 2m matrix A, with elements, and Closing the contour of integration in the upper/lower half ξ plane, depending on whether ±x − σ ≷ 0, one finds with θ(t) the Heaviside step function, θ(t) = 1 for t 0 and θ(t) = 0 otherwise. One proceeds with straightforward linear algebra. Plugging (2.46) into eq. (2.43), we first factor out ∓[i(±x − σ)] n−m θ(−x ± σ), column by column, and then 1/(m + n − k)!, row by row. The leftover factor is a Vandermonde determinant for the set {i( As such, it does not depend on σ and can be written concisely as (2.47) Assembling all factors together, we obtain and, therefore, It is worth mentioning that the above analysis would also apply to the anisotropic fishnet diagrams [6,20] associated with a theory in which the two scalar fields are given non-canonical dimensions 1 ∓ 2ω, with ω ∈ [0, 1/2) playing the role of an anisotropy parameter. The equivalent of the "BMN representation" for the m × n deformed fishnet diagram in four dimensions was worked out in refs. [16,17], along with other generalizations of the correlator, using the method of separation of variables. It was found to take the same form as for isotropic fishnets, that is (2.26), up to the replacement with f (y) = Γ(y + ω)/Γ(y + 1 − ω) and with Γ the Euler Gamma function. 5 Importantly, since the ξ-integrals stay factorized under this deformation, one can repeat the above analysis and derive a 2m × 2m determinant representation for the correlator by slightly modifying the integrals in eq. (2.45), using It is not clear, however, if the result may be further simplified and cast into the form of a determinant of an m × m matrix, as in the undeformed case, since the factorization of the determinant observed in eq. (2.48) depends very much on the explicit form of the ξ-integrals. In contrast, 2d fishnet correlators were found [7] to admit a uniform m × m determinant representation for all values of the anisotropy parameter.

FT integral
The second integrability-based representation is the so-called flux-tube representation. It arises in the Minkowskian kinematics when one considers the correlator as sitting inside a null conformal frame, as depicted in figure 4. The end-points of the correlator are then interpreted as producing beams of flux-tube particles at the positions with σ 1,2 ∈ R, along two orthogonal lightcone directions, n 2 =n 2 = 0, n ·n = 1, as explained in the caption of figure 4. A nice feature of the flux-tube representation is that it treats symmetrically the m and n lines of the fishnet graph. Each set of lines is assigned to a set of rapidities, {u i=1,...,m } and {v j=1,...,n }, representing the momenta of the flux-tube excitations, moving across the conformal square, and conjugate to the coordinates (2.52). The price to pay is that it gives rise to more complicated integrals involving both rational and hyperbolic functions. As explained in ref. [15], in this representation, the correlator takes the form Incidentally, similar generalizations appear in the context of the exact solutions to the Kardar-Parisi-Zhang equation in 1+1 dimensions, where the initial geometry provides the anisotropy, see e.g. refs. [50][51][52] for instances where f is a ratio of Gamma functions. 6 We simplified the original expression in ref. [15] by using sinh (π(x − y)) = cosh (πx) cosh (πy) (tanh (πx)−tanh (πy)) and stripping off powers of π.
x 3 x 4 x 1 x 2 n n 3 2 1 4 Figure 4. The flux-tube factorization relies on lightcone coordinates to parametrize the positions of the four points x1,2,3,4 along a null conformal frame (square). In the case of interest, i.e. for spacelike separations, x 2 ij < 0, the null square is the Poincaré patch of (conformally compactified) R 1,1 ⊂ R 1,3 , as shown in the right panel. Two points are placed at the boundary of the patch, with n · x1 =n · x4 = 1 and n,n two orthogonal lightlike vectors, n 2 =n 2 = 0, n ·n = 1. The remaining points are sitting at n · x2 = z andn · x3 =z, with z,z both real and negative.
, and with the ∆'s referring as before to Vandermonde determinants, and similarly for the v's. Integrals of this type were discussed in refs. [17,53,54] in connection with the separation of variables for SL(2, R) spin chains. In particular, in the cases where m = 0 or n = 0, the integral matches with the norm of the separated wave function, and it can be evaluated straightforwardly using eq. (54) of ref. [53]. Below we explain how to calculate the integral in the general case, m, n = 0. The mixing among hyperbolic and rational interactions makes the flux-tube integral harder to treat than the BMN one. The problem can be circumvented, however, by collecting these two types of interactions into suitable determinants.
The rational part, for instance, can be encoded in a generalized Cauchy determinant [55], where C is a "centaur" n × n matrix, mixing Cauchy and Vandermonde entries, defined for n m and two sets of variables {u 1 , . . . , u m } and {v 1 , . . . , v n } by The hyperbolic factors, on the other hand, can be seen as coming from the action of the Vandermonde differential operator on the factor i,j sech(πu i ) sech (πv j ), as shown in appendix A (Lemma 4). Accordingly, the integral can be written more concisely as (2.59) We may then perform an integration by parts, using with the arrows indicating on which side of the integrand in eq. (2.59) the derivative is acting. Note that there are no boundary terms, since the integrand is exponentially small at infinity. Using the Leibniz formula for multi-linear forms, which applies to any function F({u i , v j }), we get The next step is to disentangle the integrals in {u i } and {v j }. We may do so by exponentiating the first m lines of the matrix C using Schwinger's trick, with i0's introduced to ensure convergence. We shall omit the latter prescription in the following, since the integrand is regular when u i = v j , ∀i, j. Factoring out e iriui , row by row, we get , and where the matrixC has the same structure as C up to the (2.65) The integral I 2 is fully factorized, so it can be evaluated directly using a Fourier transform, We then treat the action of the Vandermonde differential operator on detC. We find , and defining the polynomial with σ = σ 1 − σ 2 , see eq. (2.52). Acting then on the columns ofC with P (∂ vj ) and collecting factors row by row, we find We proceed with the computation of the v-integrals, Applying the Cauchy-Binet-Andréief formula (Lemma 1) with the measure dν = dv e 2ivσ2 sech(πv)/(2π), and the function g , we obtain The n × n matrix N so defined has the special structure and F is an m × m Vandermonde matrix with elements . (2.74) Calculating I 4 using the block determinant formula and collecting the numerous constant factors, we finally get (2.76) We may now conclude. Plugging P (−ir i ) = i m m j=1 (2σ +r j +r i ) into the above expression and noting that we obtain the sought-after relation after changing variables, r i = x i − σ, and using eq. (2.52). We conclude that the flux-tube representation (2.54) is equivalent to the dual representation (2.23).

Thermodynamic limit
We now turn to the thermodynamic limit of the rectangular fishnet correlator, which we define as the limit n, m → ∞, holding fixed the aspect ratio of the rectangle, k ≡ n/m ∈ (1, ∞). We will see that the correlator scales properly in this limit, with a finite "free energy" per site, 7 7 Note that we do not include a minus sign in the definition of the free energy, for convenience.
The free energy is independent of the cross ratios (z,z or ϕ, σ), as long as they are held fixed at generic values, but it is a nontrivial function of the aspect ratio k. This analysis can be done by applying the large m saddle-point method to either of the integrals introduced earlier. However, an analytic form for the free energy F is more easily found if one combines several solutions together. We analyze here two of the previously discussed integrals, the dual and BMN integrals, leaving aside the flux-tube integral which is substantially harder to study.

Dual integral
As seen earlier, the determinant formula (2.5) can be recast as a matrix-model-like integral (2.16), or with the prefactor defined in eq. (2.18) and with the effective potential The latter comprises the Vandermonde interactions and the potential and k = n/m. This integral is similar to known matrix models, such as the O(−2) model, see [56][57][58] and references therein, and is readily amenable to the large m saddle-point method.
The key simplification at large m is that the effective potential scales large, Γ = O(m 2 ), and thus the integral is sharply peaked at its minimum, 0 = δΓ/δx i . The latter equations can be interpreted as the conditions for the static equilibrium of a symmetric collection of charges at positions {±x i } subject to Coulomb-like interactions in the presence of the external potential V (x).
In this regime, only the first and last terms on the right-hand side of eq. (3.4) remain. The latter is superficially small, but should nonetheless be kept, since it ensures the convergence of the integral at large x. Without it, the charges would run away, x i → ∞. Instead, their course stops at large distance ∼ m, owing to the linear scaling at large x, V (x) ≈ |x|/m, which in turn determines the large m scaling, The first term in eq. (3.4) is of order O(m 0 ), for β ∈ (0, ∞), and it acts at the lower bound of the domain. It generates a logarithmic repulsion, which may be interpreted as coming from n − m nondynamical charges which push the dynamical charges away from x 2 = σ 2 . Altogether, we expect that the charges are, after rescaling, densely distributed on a compact domain, when β is finite. 8 Furthermore, it readily follows from eq. (3.6) that the spacetime parameters ϕ, σ drop out in the large m limit, as long as ϕ, σ stay of order O(1) at large m. Hence, one may approximate the potential by in the thermodynamic limit discussed here. A more general scaling where the spacetime parameters scale large with m will be discussed later on, in section 4.

Saddle-point equation
We may now develop the standard analysis and introduce a density which we normalized for convenience such that (3.10) We expect it to admit a smooth description at large m over the support of the distribution C. Since the scaled potential (3.8) is strictly convex downward and diverges at the origin, the contour is a single interval C = (a, b) with 0 < a b < ∞, for finite β. In particular, for a small interval b ∼ a, the charges fill the bottom of the potential described by the Gaussian matrix model. It corresponds to the low density regime, β ∼ 0, that is, k → ∞. In the opposite regime, when k → 1, the logarithmic repulsion at x = 0 disappears and the left boundary a → 0.
The large m scaling of the integral becomes manifest after rewriting the effective potential (3.3) in terms of the density, and approximated using the asymptotic behavior, for large z 1. Defining then the scaling function 14) Figure 5. Plot of the density ρ and shifted potentialṼ (x) = 2πx − log x + λ at low density β ≡ 0.015.
with the limit taken at fixed k, we get where with λ a Lagrange multiplier. A derivative in x eliminates λ and yields the singular equation with − denoting Cauchy principal value. The equation can also be cast in a more classic form, after unfolding the interval, (x, y) ∈ C = (−b, −a) ∪ (a, b), and imposing ρ(−x) = ρ(x).

Spin-chain mapping
Equation (3.19) describes a well-known two-cut matrix-model problem, which appears in various contexts and can be solved exactly (see e.g. refs. [59][60][61][62][63][64] for reviews of general methods). In particular, it was found to control the classical limit of the XXX −1/2 Heisenberg spin chain [65,66]. In this case, one considers a collection of S magnons propagating on a closed spin chain of length J. The magnons' momenta are parametrized by a set of roots {u j , j = 1, . . . , S}, quantized by the Bethe ansatz equations (3.20) The classical limit corresponds to J, S → ∞ while keeping the density S/J fixed. In this limit, the roots scale as u i = O(J) and the logarithms of the Bethe ansatz equations simplify (see [ One also verifies that the normalizations agree if β = S/J. Thanks to this identification, one may read off the solution to our problem from the solution given in refs. [65,66]. It is expressed in terms of elliptic integrals. In particular, the parameters a, b and β are given parametrically as with K(q) and E(q) the complete elliptic integrals of the 1st and 2nd kind, respectively, Here, q ∈ (0, 1) with the lower/upper bound corresponding respectively to β → 0 and β → ∞. One may also express the density ρ(x) in closed form, in terms of the incomplete elliptic integral of the third kind. (See eqs. (C.7) and (C.8) in appendix C of ref. [65].) We plot ρ(x) in figure 5 for illustration. We will not need its explicit analytic expression here, but it will be useful to know its derivative with respect to β, δρ = ∂ β ρ.
The derivative of the density is rather simple to construct. Differentiating both sides of eq. (3.17) with respect to β, the potential drops out. Using also the vanishing of ρ at the boundaries, ρ(a) = ρ(b) = 0, one is left with a simpler (homogeneous) problem In the Coulomb gas picture, this equation describes a collection of charges with no external potential, but subject to hard-wall boundary conditions at a and b. So, despite the absence of a potential, the charges do not run away but accumulate close to the boundaries, c = a, b, where δρ(x) ∼ 1/ |x − c|.
Taking a derivative of eq. (3.25) with respect to x removes the left-hand side of the equation, and the solution is immediately obtained using a standard inversion formula, The only freedom here is the overall normalization, which is fixed using taking into account again that the boundary terms vanish.

Differential equation I
The next step is to calculate the scaling function (3.15). This may be done in principle using the exact solution for ρ given in refs. [65,66]. Its direct integration proves difficult however. So, here we shall consider the much simpler problem of determining df /dβ using the density derivative δρ = ∂ β ρ.
Observe first that eq. (3.15) simplifies when evaluated on the saddle point. Acting with dxρ(x) on both sides of eq. (3.17) and using the normalisation condition (3.10), one finds the relation which may be used to eliminate the double integral in eq. (3.15), with f 0 given in eq. (3.16). A derivative with respect to β then yields using that δ = δ since ρ vanishes at the boundary. One may further simplify the right-hand side of eq. (3.30), by acting on eq. (3.17) and eq. (3.25) with 2 dxδρ(x) and 2 dxρ(x), respectively, and taking their difference, The integral on the right-hand side of the equation for λ, as well as the one defining ∂ β λ, see eq. (3.25), can be taken immediately using eq. (3.26), Finally, eliminating λ yields with a, b, β related to each other through eq. (3.23).
The differential equation (3.33) can be integrated directly around particular points, such as β = 0 or ∞, after evaluating its right-hand side using eqs. (3.23). It is less straightforward to integrate it for general β, given the complicated dependence of the coefficients a, b on β. Fortunately, we will not need to deal with this problem here, as the solution will come for free after combining eq. (3.33) with a similar equation controlling the thermodynamic limit of the BMN integral.

BMN integral
Consider now the BMN integral (2.26). In this case, instead of a single matrix-model integral, we have an ensemble of integrals, labelled by integers {a }. However, at large m, we expect it to be dominated by the lowest modes, with a = 1, ∀ = 1, . . . , m. This is because the external potential V a (u) = (m + n) log (u 2 + a 2 /4) (3.34) has its minimum at a = 1 and generates power suppressions for the higher a's in the thermodynamic limit. Owing to this truncation, there is no dependence on ϕ, and similarly σ drops out, for σ, ϕ = O(1).
In summary, focusing on the leading saddle, one may write Furthermore, one notes that the potential (3.34) scales linearly with m, for all u's. Hence, there is no need to rescale the rapidities here and we can proceed directly to the thermodynamic limit with the density which we expect to be smooth for u = O(m 0 ). Varying the integral (3.35), we get assuming that the density is even and supported on a compact interval (−B, B), as suggested by the symmetry and convexity of the potential. We seek a solution obeying the usual boundary conditions and normalized such that  (3.41) The problem gets harder in the intermediate regime k ∈ (1, ∞). It can nonetheless be solved exactly by following a method introduced in ref. [67] to tackle a similar looking matrix-model equation.
To begin, we define the function It is analytic in C, except along the intervals (−B ± i/2, +B ± i/2), where it has square-root branch cuts, and at u = 0, where it has a simple pole, r(u) ∼ (k + 1)/u. We also note that the function has a pole at u = ∞, with residue The function r relates to the standard (single-cut) resolvent implying that one may recover the density ρ(u) from the discontinuities of r(u) across either cut, In particular, it follows from this relation that r(u) must be continuous at the branch points u = ±i/2 ± B, since the density vanishes at u = ±B, see eq. (3.38).  The main motivation for introducing a resolvent in this way is that the saddle-point equation (3.37) now takes the simple form It implies that r(u) extends to an anti-periodic function of period i on the second Riemann sheet, that is, after continuing below the cuts, see figure 6.
In addition, the resolvent r possesses important reality and reflection properties, with * denoting complex conjugation, which stem from the fact that ρ is a real symmetric function. It then suffices to find r on the first quadrant, in order to determine it everywhere else. Taking the cuts into account, we should in fact consider the open domain D shown in the left panel of figure 7, which is obtained by removing the segment between u = i/2 and u = B + i/2 from the first quadrant.

Schwarz-Christoffel mapping
The function r obeys simple functional relations, but is defined over a complicated domain. The next step is to make the domain more regular with the help of a Schwarz-Christoffel (SC) transformation. Namely, D can be seen as the interior of a polygon (with a vertex at infinity) and as such can be mapped conformally to the upper half-plane H = {w ∈ C : Im w > 0} with the boundary ∂D mapping to the real w axis. The map is constructed canonically, using the interior angles at the vertices of ∂D to determine the exponents in the SC transform w → u. Fixing u = w = 0 and u = w = ∞ for simplicity, we get with w 4 < w 3 < w 2 < 0 the images of u 4 = i/2 + i0 , u 3 = i/2 + B and u 2 = i/2 − i0, respectively, and with A > 0 an arbitrary normalization factor, introduced for later convenience.
We may now easily express the resolvent r as a function of w. The key observation is that the map F : w → F(w) = 1/r(u(w)) 2 (3.49) is a holomorphic bijective function of H to itself, that is, a real Möbius transformation. This is manifest in the neighbourhood of w = 0 or w = ∞. At these points, r(u) ∼ (k ± 1)/u and √ w/u is constant.
Hence, one can always write for a suitable choice of the parameters A, w 2,3,4 , see eq. (3.52) below. The claim is that eq. (3.50) holds globally, for all w ∈ H.
To justify the claim, it is enough to show that F extends to a continuous function of the closed upper half-plane taking real values on the real axis. The Schwarz reflection principle can then be used to conclude that F extends to an entire function on C 9 and thus F(w) = w given the asymptotic behaviours at 0 and ∞. Now, it is easy to show that F preserves the boundary ∂H. This follows directly from the reality and periodicity properties of r. Namely, r(u) is respectively real for u ∈ (0, ∞) and purely imaginary for u ∈ (0, i∞), because of eq. (3.47). It is also purely imaginary just below the cut, i.e. u ∈ (u 2 , u 3 ), owing to eq. (3.46), and similarly for u ∈ (u 3 , u 4 ), thanks to the discontinuity formula (3.46) and the reality of ρ. 10 Furthermore, F is continuous on ∂H, including at the fold-point w = w 3 , since ρ(u) vanishes at u = u 3 , and at the origin w = 0, where F → 0. Hence, as desired, F maps ∂H into itself. Finally, one may determine the unfixed parameters (A, w 2,3,4 ) in the SC transformation, by matching both sides of eq. (3.50) at the marked points on the boundary. The analysis is deferred to appendix B. The results can be written concisely in terms of the elliptic data a, b, q introduced in eqs. (3.23), with β = 1/(k − 1) and with the modulus q = 1 − w 2 /w 4 ∈ (0, 1).

Duality map
Prior to moving on to the calculation of the free energy, let us open a parenthesis here to comment on the relation between the above solution and the one constructed earlier. The two solutions should be images of each other, as they are based on two equivalent representations, for any m and n. In practice, however, it is not easy to follow the Mellin-Barnes transformation used in section 2.2 all the way to the large m limit. One may nonetheless find the duality map directly in this limit from a comparison of the two saddle-point solutions.
Recall first of all that the interval (a, b) determines the support of the dual density, denoted ρ dual (x) here, to avoid confusion. Looking at (3.52), we see that this interval can be identified with the segment (w 2 , w 3 ) of the w plane. More precisely, the above relations suggest to set thus identifying the image of the resolvent r(u) with the dual x plane. The reciprocal relation identifies the u plane with the image of a suitable resolvent r dual (x) for the dual problem, This relation can be worked out from the elliptic parametrization (3.48) together with (3.53). One may also proceed as follows. Using the map (3.53) and properties of r(u), one easily sees that the inverse map r dual (x) = −r dual (−x) is an analytic function of x with square-root cuts along (−b, −a) ∪ (a, b) and poles at x = 0 and x = ∞, with residues respectively. It must also obey the equation for any x ∈ (a, b), since on this interval Im u = i/2 according to eq. (3.53), see table in figure 7. These requirements define a Riemann-Hilbert problem whose solution is given by , with ρ dual (x) the density constructed in the previous section. In particular, equations (3.55) and (3.56) are easily checked using (3.10) and the dual saddle-point equation (3.18). At last, r dual (x) is seen to be a resolvent for the dual problem, since for any x ∈ (a, b). Hence, nicely, the large m duality map merely interchanges rapidities and resolvents of the two problems.

Differential equation II
We may now calculate df /dk using the density of the BMN integral (3.35). The starting point is the saddle-point expression obtained from eq. (3.35), 59) which is extremal on the just-constructed solution ρ, with µ a Lagrange multiplier for ρ(u)du = 1. These relations imply that with δρ = dρ/dk, after eliminating the double integration with eq. (3.60) and using δ ρ = δρ = 0. We then introduce an auxiliary function, which returns −df /dk when evaluated at u = 0. This function is nothing but an integrated version of the resolvent r(u), with given large u behavior, This equation can be integrated with the help of eqs. (3.48) and (3.50). It yields (3.64) up to a constant of integration, fixed by the large u asymptotics, using u ∼ (k − 1) √ w 1 and eq. (3.52). Lastly, the limit u → 0 yields using u ∼ (1 + k) √ w 1 (see appendix B) and switching to variables a, b with eq. (3.52).

Free energy
We shall now derive a parametric representation for the scaling function F (k) = f (k)/k. From eq. (3.14) we see that F (k) is normalized by the rectangle area mn rather than m 2 .

Parametric representation
A key simplification arises when combining our two differential equations, eqs. (3.33) and (3.66). Namely, one finds that the function f can be extracted directly, without any integration. This is possible because the equations are linearly independent and have only one solution in common. Taking the sum of the two equations eliminates the derivatives of f , after using eqs. (3.33), (3.66), (3.16) and β = 1/(k − 1). 11 We assume here that the two differential equations are set in the same parametrization, that is, that the moduli q are the same. This is the case in the range β ∈ (0, ∞) ↔ q ∈ (0, 1), where the elliptic parametrization is bijective, with the lower bound being observed in the limit q → 0 and with β ranging between 0 and ∞ for q ∈ (0, 1). Hence, all values of k = 1 + 1/β ∈ (1, ∞) are covered once and only once. As a double check, one may verify that this expression solves each differential equation separately. For this purpose, it is useful to have the β derivatives of a and b, which follow from eq. (3.23) and the chain rule. Introducing F (k) = f (k)/k and using eq. (3.23), one arrives at where we recall that with K(q) and E(q) the complete elliptic integrals, see eq. (3.24), and q ∈ (0, 1).
The system of equations (3.71) and (3.72) is our final expression for the free-energy density. It determines F (k) parametrically over the whole domain k ∈ (1, ∞). Figure 8 shows the comparison between this scaling function and a numerical estimate of the ladder determinant (2.6) for large m and n. As one can see, the agreement is excellent for all accessible values of k.
Note that one may first simplify the determinant to facilitate the comparison. Since the cross ratios drop out in this regime, one can conveniently set ϕ = σ = 0, or equivalently send z,z → −1, in eq. (2.5), giving Below, we study in detail the behaviour of the free energy at k = ∞ and k = 1, using known asymptotic expressions for the elliptic integrals at q = 0 and q = 1. These points are of particular interest since they correspond to the ladder limit n m and to the symmetric point n = m, respectively.

Ladder limit
The limit k → ∞ corresponds to q → 0 + . The elliptic integrals are regular at this point and admit convergent expansions in integer powers of q, from which it follows that using eq. (3.72). Expanding further and plugging the series inside eq. (3.71), one finds (3.76) Notice that the half-integer powers of 1/k in the expansion of q cancel out in F . Hence, F admits a regular expansion in 1/k, up to a logarithm ∼ log k/k. The series is seen to alternate and it approximates the exact curve extremely well throughout the entire domain, including the point k = 1, where (3.76) is found to converge rather quickly towards the exact value F (1) = log (π 2 /4). It suggests that the series converges for all k ∈ [1, ∞). The leading term in this expansion has a simple interpretation. As seen earlier, the limit k → ∞ corresponds to a dilute gas approximation, which means that the free energy is a multiple of the free energy of an individual ladder L p at large p. This is easily verified using the definition (2.9). At large weight, the polylogarithm series (2.9) truncates, assuming |z| < 1 for convenience, and the sum 2p j=p in eq. (2.9) is dominated by the terms with j ∼ 2p. One can thus take the large p limit of the summand at fixed l = 2p − j, using Stirling's approximation for the factorials, and perform the sum over l = 0, 1, . . . , It yields lim p→∞ 1 p log [L p (z,z)] = log 4 = 1.386... , in agreement with F (∞).

Square limit
The limit k → 1 corresponds to q → 1 − and in this limit Subleading corrections have a complicated pattern, involving not only powers of (k − 1) but also of its logarithms. This feature goes back to logarithms in the expansion of the elliptic integrals around q = 1, with the dots standing for regular series in (1 − q). A double expansion is easily generated as follows. First, define the variable which is small when k → 1 and changes sign under k → 1/k, and then introduce the ansatz with the expansion coefficients, α and α n , being functions of log ζ. The latter coefficients are determined self-consistently, by plugging (3.82) inside eq. (3.81), expanding in powers of ζ and matching the coefficients on both sides of the equation. For instance, at the leading order in ζ, one finds 83) or equivalently The full solution for α = α(log ζ) can also be given in closed form, with W −1 (x) the second branch of the Lambert W function, defined for x ∈ (−1/e, 0) and such that when x → 0 − . Cancelling the terms at higher orders in eq. (3.83) determines the higher coefficients α n recursively. One finds that α n = 0 when n is odd, in line with the k ↔ 1/k (dihedral) symmetry of the correlator, and one observes that the non-vanishing coefficients can be expressed in terms of α. The first few terms read which is thus an even regular series, up to the tails of logarithms encoded in eq. (3.84). We remark that the expansion (3.87) has a very small range of validity if α is given by eq. (3.84).

Comparison with periodic fishnets
It is instructive to compare our formula with the large-order behaviour of fishnet graphs with doubly periodic boundary conditions. In his pioneering paper, A. Zamolodchikov [20] obtained the expression for the m × n fishnet torus partition function Z m,n , in any spacetime dimensions D, using integrability and bootstrap conditions. Assuming the thermodynamic scaling, he found that where the so-called critical coupling g c is a constant which depends on D but not on the fishnet lengths m, n. The explicit expression for D = 4 is given by [20] − log g 2 c = log Clearly, this result bears no similarity at all with our general expression for F (k). It is also numerically different, over the whole physical range k ∈ (1, ∞); it is just over half of the lowest value achieved by F (k), namely F (1) = log (π 2 /4) = 0.903.... From the 2d lattice perspective, this discrepancy and, more generally, the dependence of the scaling function F on the aspect ratio k appear surprising at first sight. They go against the common wisdom that the thermodynamic free energy is extensive and independent of the boundary conditions. 12 Put differently, they suggest that the boundary conditions set by the four-point fishnet correlator are atypical and "severe" enough to penetrate deeply into the bulk of the lattice, in such a way that the uniform periodic distribution may hold only locally, far away from the boundary, if at all. This sensitivity to the boundary conditions is likely to hold for other correlators. In fact a similar phenomenon may be observed for half-periodic fishnets, mixing open and closed boundary conditions, such as the one shown in the right panel of figure 1. An interesting aspect of these mixed boundary conditions is that they provide a regularization of the otherwise ill-defined doubly periodic fishnets, which suffer from UV/IR divergences. They yield finite correlation functions, which depend on several cross ratios and which may be expanded over a complete basis of states in principal series representations, as discussed in detail in ref. [5] in a particular case. Doing so, one projects over states carrying arbitrary scaling dimension ∆ along the radial direction. 13 . Fishnets of this type were considered in the thermodynamic limit in ref. [22], for certain states, and their free energy was calculated for a large range of scaling dimensions ∆ using the Thermodynamic Bethe Ansatz (TBA) equations of the fishnet theory. Importantly, they were found to scale as in (3.88) only for ∆'s that are much smaller than the fishnet size. The latter requirement is understood as a condition for the validity of the continuum description of large fishnet graphs, which is governed in the case at hand by the 2d non-linear AdS 5 sigma model. If on the contrary the dimension ∆ scales large with the system size, then the corresponding state in the sigma model is very excited, with a finite energy density, and the free energy density departs significantly from the ground-state behavior (3.88).
Another way of phrasing the situation is in terms of the so-called graph building operator T m [4,20], which is the integral kernel which adds a rung to periodic fishnet graphs with m radial lines. It may be interpreted as a finite-time evolution operator for an integrable spin chain, with spins in a suitable representation of the four-dimensional conformal group [4,17]. Its logarithm defines a Hamiltonian for a system of length m, which is well approximated by the Hamiltonian of the AdS 5 sigma model, in the large volume limit m → ∞, up to higher-derivative interactions associated with 1/m corrections. The point is that this expansion can be trusted, and the latter corrections neglected, only when the quantum numbers of the state, and in particular the scaling dimension, are much smaller than the length m. One may expect this condition to be observed for m × n half-periodic fishnets, at generic values of the cross ratios, as long as both m and n are large -with the limit of large number of rungs (or large discrete time), n → ∞, being used here to project on the 2d low-energy states. However, it is not obeyed, even for large m and large n, for the special configurations of the boundary points which do not have a smooth 2d limit and extract contributions from the high-energy part of the spectrum of log T m . With that in mind, it appears less surprising to find that the boundary conditions set by our four-point function are strong enough to drive the system away from the scaling (3.88). Indeed, the fishnet lengths m, n coincide in this case with the scaling dimensions of the operators sitting at the boundary, which therefore scale large in the thermodynamic limit. In other words, the interpretation would be that it is because of the "roughness" of the boundary conditions that we observe a loss of universality in the thermodynamic scaling.
We should stress that strong sensitivity to the boundary conditions has also been observed in other solvable statistical models. The most famous example is given by the six-vertex model and related exactly solvable models of statistical mechanics, see [70,71] and references therein, where the partition function with domain-wall boundary conditions was found to scale very differently from the torus partition function in the large volume limit. The interpretation here is that the former partition function is subject to the formation of a so-called arctic curve [72,73], which separates two different phases of the model. Namely, the configurations appear nearly frozen across a macroscopic domain near the boundary and only start fluctuating deep inside the lattice, where they approach the disordered phase. The free energy is not distributed evenly throughout the lattice but instead undergoes a (more or less) sharp transition at the place where the two phases co-exist. It has also been suggested that this phenomenon could result from the roughness of the boundary conditions [70,71] at least in some situations.
It would be interesting to see if this analogy can be made more precise for the fishnet lattices. It would also be interesting to see if the phenomenon discussed here admits a dual AdS interpretation, despite the tension with the 2d low-energy description (3.90). In the next section, we will see some indirect evidence of a connection with classical string theory in AdS 3 , as we consider a generalized thermodynamic scaling limit which combines large fishnet and short spacetime limits.

Short-distance limit and spinning string
Let us now discuss a more general scaling limit, in which the spacetime parameters scale large with the fishnet lengths m, n. We shall focus on the situation where |σ| ∼ m, which corresponds to the Euclidean short-distance regime, x 1 → x 3 or x 2 → x 3 , depending on the sign of σ, and use the dual integral (3.2). 14

Mapping with spinning string
The dual equations stay essentially the same in the general scaling limit, if not for the potential (3.8), which is replaced by Introducing the density (3.9) and the parameter which is held fixed, with β = 1/(k − 1), in the limit |σ|, m, n → ∞, the saddle-point equation (3.18) becomes with ξ a b. Our previous discussion relates then to the limit ξ → 0. Now, equation (4.3) turns out to be the same as the finite-gap equation [74][75][76] for a classical string in AdS 3 × S 1 , constructed explicitly by Frolov and Tseytlin [27]. 15 It describes a folded string moving in AdS 3 with spin S and global time energy E and boosted with a momentum J along a big circle in S 1 . In this setup the parameter ξ plays the role of the so-called BMN coupling, with √ λ/2π the string tension and with λ the 't Hooft coupling of the dual gauge theory. At weak BMN coupling, the string theory description merges smoothly with the spin-chain analysis reviewed in section 3.1.2, as discussed in detail in refs. [65,66,74]. However, as the coupling gets bigger, the classical string analysis takes over.
The solution to eq. (4.3) is known for any ξ and may be found explicitly in refs. [74][75][76]. Here, we simply quote the expressions for the parameters of the distribution, and q = 1 − a 2 /b 2 ∈ (0, 1), which generalize to ξ = 0 the formulae in eqs. (3.23). Recall that β controls the normalization of the density in our convention. In the string theory, it relates to a linear combination of the string quantum numbers [74][75][76] It reduces to the spin-chain expression, see comment below (3.22), in the weak BMN coupling limit ξ → 0, that is, when the "anomalous dimension" E − S − J ∼ 0.

Short-string regime
One may then derive a differential equation for the scaling function, by following the same lines of analysis as before. One only needs to replace the potential used in eq. (3.30) by its deformed version, with δρ = ∂ β ρ given in eq. (3.26). It yields the differential equation with the derivative taken at fixed ξ and with f 0 given explicitly in eq. (3.16). The extra piece on the left-hand side stems from the prefactor of the dual integral (2.16), d(z,z) m ≈ e m|σ| = e 4πm 2 ξ/β , (4.10) in the limit m, |σ| → ∞.
Integrating eq. (4.9) for general β and ξ is beyond the scope of this paper. Here, we will focus on taking the strong coupling limit ξ 1 at fixed β. This limit maps to the short-string domain on the string theory side, with the spin, energy and momentum obeying the flat-space dispersion relation S ∼ E 2 − J 2 [27]. 16 The small and large β limit coincides then with J 2 S and J 2 S, respectively. At the level of eq. (4.3), the short-string regime corresponds to an expansion around the bottom of the potential V (c) = 0 in the limit ξ → ∞. Expanding eqs. (4.5) around this point yields where c = (1+ 1 + 16π 2 ξ 2 )/4π ∼ ξ. Inserting the expressions in the right-hand side of the differential equation (4.9), using 12) one finds that the terms linear in ξ cancel out, such that This equation is significantly simpler than the one encountered earlier, for small ξ, and it can be integrated directly, using eq. (3.16) for f 0 and the boundary condition lim β→0 (β 2 f ) = 0. Re-expressing it in terms of the original variables, we find log Φ m,n = mn log (2|σ|) + 1 2 (m 2 log m + n 2 log n − (m + n) 2 log (m + n)) + 3mn which holds when |σ| m, n 1. It is manifestly symmetric under m ↔ n. One reads from it that Φ m,n ∼ |σ| mn in the short-distance limit, |σ| → ∞. Note that it sends the exponent infinitely far away from the predicted one (3.88) for periodic fishnets. This scaling is nonetheless in perfect agreement with the UV power counting of the fishnet integral, with each loop giving rise to a logarithmic divergence ∼ |σ|. It is quite intriguing to see this simple field theory behaviour arising here from a short string analysis.
Let us add finally that subleading terms can be easily produced by keeping higher orders in (4.11). For illustration, one easily finds for the first few corrections to (4.15). They constitute corrections to the flat-space regime, coming from the curvature of AdS, in the string theory interpretation.

Hankel determinant of factorials
We can compare the "stringy" prediction (4.15) with a direct evaluation of the determinant of ladders. As we will now explain, the latter simplifies drastically at large σ and can be given in closed form for all m and n. We first note the large σ behaviour of the ladders, choosing σ → −∞ for convenience, The highest power of |σ| = 1 2 [− log (zz)] dominates in the sum (2.9), using the fact that Li j (z) ≈ z in the limit |σ| = −σ → ∞, or z,z 1. Inserting the above estimate inside the determinant (2.6) then yields the simple expression for the leading large |σ| behaviour. Now, determinants of Hankel matrices of factorials form nice sequences, which can be written concisely in terms of Barnes' G-function, see eq. (3.12). In particular, a classic mathematical result yields for n = m, (4.19) The general formula, valid for any n m, is given by [78] det((i + j + n − m − 2)!) 1 i,j m = G(m + 1)G(n + 1) G(n − m + 1) , (4.20) and it can be derived by evaluating the determinant with the method of orthogonal polynomials, as explained in appendix C. It combines neatly with the prefactor (3.12) such as to give Φ m,n ≈ (2|σ|) mn G(m + 1)G(n + 1) which is symmetric under m ↔ n, as it should be. The large m, n limit follows immediately from the asymptotic behaviour of Barnes' G-function at large argument, eq. (3.13), reproducing perfectly the stringy result (4.15). One may also find concise expressions for the subleading logarithmic corrections at large |σ|. They originate from the lower terms in the polynomials in |σ| in the ladders, with K the modified Bessel function of the second kind. Plugging this expression inside the determinant and expanding at large σ, for several values of m, n, we found that the corrections can be parametrized in terms of symmetric polynomials in m and n of increasing degrees, with Φ 0 m,n denoting the leading order expression (4.21). These polynomials are readily seen to match with the stringy curvature corrections (4.16), when m, n 1.

Conclusion
In this paper, we examined simple four-point fishnet integrals in four dimensions, for operators carrying finite or large scaling dimensions, m, n. The equivalence between the various integral representations allowed us to obtain closed-form expressions for the free-energy density controlling the large-order behaviour of these diagrams. It was found to be independent of the cross ratios (away from singular points), in line with common wisdom on the scaling of boundary data in the thermodynamic limit. However, the "universality property" of the free energy appears incomplete, owing to non-trivial dependence on the fishnet aspect ratio k = n/m. The dependence on k of the free energy raises a number of questions. In particular, it suggests that the fishnet partition function studied in this paper may be subject to the formation of an arctic curve, as seen in the six-vertex model and related exactly solvable models of statistical mechanics [70][71][72][73] with domain wall boundary conditions [79]. The analogy is quite neat at the mathematical level, with the latter set of boundary conditions giving rise to a determinant, while the periodic ones are found to be governed by a more involved (linearized) TBA equation.
It would be interesting to see if this analogy holds up under a more detailed analysis. Our measure of the free energy density is rather crude, since it averages over the entire graph. It would be interesting to dive more deeply into the diagram and extract information about the local distribution. This may perhaps be done by deforming the boundary conditions, adding excitations on the local operators at the boundary, or by considering higher-point functions, with additional operators. (For instance, a five-point function with a Lagrangian insertion might help exploring local properties of the fishnet graphs.) One may also draw inspiration from the methods used to construct the arctic curve in the six-vertex model, see e.g. refs. [80,81].
Recently, a lot of progress has been made in obtaining alternative descriptions of the fishnet diagrams in terms of dual systems in AdS. In particular, we mentioned that cylindrical fishnet graphs are expected to relate to the non-linear AdS sigma model in the continuum limit. In light of the aforementioned discrepancy, it is not immediately clear if this description applies to the thermodynamic limit of the open-string-like correlators studied in this paper. The question remains as to whether a dual AdS description exists in these cases as well. This description may not be smooth throughout the entire system or given in terms of a local 2d quantum field theory, but it may nonetheless be formulated entirely in the AdS space. The AdS string-bit formulation of the fishnet graphs developed recently in refs. [23][24][25]82] -dubbed the fish-chain model -may shed light on it, as it is does not assume a low-energy approximation and is naturally tied to the regime of large scaling dimensions. Still, some work appears needed to cast it in a form that is suitable for studying the thermodynamic limit.
The generalized scaling limit discussed in section 4 is an interesting corner for exploring a potential new connection with a dual AdS description. In this case we saw that the saddle-point equation agrees perfectly with the finite-gap equation for a rigid string spinning in AdS 1 × S 1 . Some aspects of the mapping, such as the emergence of the internal circle S 1 or the relation between the spacetime parameter of the fishnet correlator and the BMN coupling of the string, are quite intriguing and it is not yet clear if this mapping results from a mathematical accident or stands as an indication of a genuine correspondence.

A Technical lemmas
In this appendix, we provide the key lemmas used for the exact evaluation of the BMN and FT integrals.
The following Lemma yields an integral relation for the product of two determinants defined by two sets of functions.
The next Lemma permits one to convert sums over discrete integers into contour integrals.

Lemma 2 (Mellin-Barnes summation formula)
Let C = ε + iR for some ε ∈ (0, 1) be a contour in the complex plane, parallel to the imaginary axis and y ∈ C such that Re y > 0. The Mellin-Barnes summation formula allows one to replace a summation over integers by an integral over C, using Resw=a[π/ sin (πw)] = (−1) a . Furthermore, we note the second identity for Re y > 0 and w ∈ C, π sin (πw) where the right-hand side can be viewed as a Mellin transform of a Fermi factor.
We introduce now a lemma about derivatives of the inverse of hyperbolic cosines which will be essential to prove Lemma 4.
Proof. We prove it by induction. Assume it is valid for some j 1, we get
The next Lemma allows us to transform a Vandermonde determinant over hyperbolic tangents into a Vandermonde determinant over partial derivatives.

B Elliptic parameters
In this appendix, we determine the parameters A, w 2 , w 3 , w 4 in the elliptic parametrization (3.48), (3.50) by fixing the mapping at the boundary points.
Matching the asymptotics. The asymptotic behaviours at small u and at large u are easily matched on both sides of eq. (3.50). For u small, w is small and in this limit eq. Hence, We can proceed similarly for the large u regime, which maps to large w. Equation Matching boundary points. The remaining parameters are fixed by evaluating the solution at special points on the boundary. Namely, enforcing that the pre-images of {0, w 2 , w 3 , w 4 } are {0, i 2 − i0, i 2 + B, i 2 + i0}, we obtain the relations √ −w 4 = E(p)K(q) + K(p)E(q) − K(p)K(q) K(q) = π 2K(q) , (B.10) using the Legendre relation for elliptic functions of conjugate moduli, p + q = 1, to simplify the numerator. Summary. We found where E(γ, q), F (γ, q) are the incomplete elliptic integrals obtained by changing the upper limit in eq. (3.24) from π/2 to γ. We arrive at eq. (3.52) by replacing the elliptic integrals by the parameters a, b, β introduced in eq. (3.23).

C Orthogonal polynomials
In this appendix, we calculate the short-distance limit of the determinant of ladders, eq. (4.21), using the method of orthogonal polynomials, see ref. [59] for a review. First, recall the integral representation (2.21) of the ladder L p (z,z). In the limit σ → −∞, or equivalently z,z → 0, it yields recalling that = n − m and using G(z + 1) = z−1 i=0 i!.