From Holonomy of the Ising Model Form Factors to n-Fold Integrals and the Theory of Elliptic Curve

We recall the form factors $ f^{(j)}_{N,N}$ corresponding to the $\lambda$-extension $C(N,N; \lambda)$ of the two-point diagonal correlation function of the Ising model on the square lattice and their associated linear differential equations which exhibit both a ``Russian-doll'' nesting, and a decomposition of the linear differential operators as a direct sum of operators (equivalent to symmetric powers of the differential operator of the complete elliptic integral $E$). The scaling limit of these differential operators breaks the direct sum structure but not the ``Russian doll'' structure, the ``scaled'' linear differential operators being no longer Fuchsian. We then introduce some multiple integrals of the Ising class expected to have the same singularities as the singularities of the $n$-particle contributions $\chi^{(n)}$ to the susceptibility of the square lattice Ising model. We find the Fuchsian linear differential equations satisfied by these multiple integrals for $n=1,2,3,4$ and, only modulo a prime, for $n=5$ and 6, thus providing a large set of (possible) new singularities of the $\chi^{(n)}$. ...


Introduction
This paper displays a selection of works and results that have been obtained by the authors in collaboration with B.M. McCoy, W. Orrick and J.-A. Weil. It also provides new ideas and viewpoints at the end of Subsection 3.4, in Section 5 and in the Conclusion 10. We also give new results of linear differential operators modulo prime that had not been published before in Appendix C.3 and Appendix C. 4.
The two dimensional Ising model in zero magnetic field is, historically, the most important solvable model in all of theoretical physics. The free energy [91], the partition function on the finite lattice [67] and the spontaneous magnetization [92,125] were computed long ago by Onsager, Kaufman and Yang. These computations, and subsequent studies of the correlation functions [68,123], form the basis of scaling theory and of the renormalization group approach to critical phenomena.
Let us first recall the form factors [25] of the lattice Ising model. Our starting point will be the expansions of the diagonal correlations in an exponential form [123], both for T < T c where E h and E v are the horizontal and vertical interaction energies of the Ising model. We will restrict in the following to the isotropic Ising model. For diagonal correlation functions, there is no difference between the isotropic and anisotropic models: the diagonal correlations are functions of the modulus k = sinh(2E v /k B T ) sinh(2E h /k B T ). The difference comes with off-diagonal modes and is sketched in [26].
When the exponentials in (1.1) and (1.2) are expanded, the correlations can also be written in what is called a "form factor" expansion The form factor f (j) N,N is interpreted as the "j-particle" contribution to the two-point correlation function. It is natural to consider λ-extensions [83,123]  which weight each f (j) N,N by some power of λ, and to interpret λ as being analogous to a coupling constant in a quantum field theory expansion. Such λ-extensions naturally emerge from the Fredholm determinant framework in [123]. We will present new integral representations for F N,N in Section 2. We will see that they are much simpler, and more transparent, than the forms obtained from C(M, N ) of [123] by specializing to M = N .
On another hand, Jimbo and Miwa introduced in [60] an isomonodromic λ-extension of C(N, N ). Remarkably this more general function C(N, N ; λ) also satisfies [27,93] the Painlevé VI equation (1.7). The motivation of introducing an isomonodromic parameter λ, in the framework of isomonodromy deformations, is, at first sight, quite different from the "coupling constant Fredholm-expansion" motivation at the origin of the form factor λ-extensions (1. 5) and (1.6). In [27] we have shown that these two λ-extensions are actually the same by demonstrating that the recursive solutions of (1.7), analytic 2 in t 1/2 , agree with (1.5) and (1.6) where the f (j) N,N 's are obtained from C ± (N, N ; λ), the λ-extension of C ± (N, N ). The normalization condition (1.8) fixes one integration constant in the solution to (1.7). We find that the second integration constant is a free parameter, and, denoting that parameter by λ, we find that our one-parameter family of solutions for C − (N, N ) can be written in a form structurally similar to the right hand side of (1.5). Furthermore, we have confirmed, by comparison with series expansions of the multiple integral formulas for f (j) N,N derived in Section 2, that this family of solutions is, in fact, identical to C − (N, N ; λ) as defined in (1.5). Similarly, the condition (1.9) gives rise to a one-parameter family of solutions for C + (N, N ) that is identical to (1.6). The form factor expressions for the two-point correlation functions C(M, N ) of [87,88,94,93,123,124] are obtained by expanding the exponentials in (1.1), and (1.2), in the form given in [123] as multiple integrals, and integrating over half the variables. The form of the result depends on whether the even, or odd, variables of [123] are integrated out. For the general anisotropic lattice, one form of this result is given, for arbitrary M and N , in [93]. When specialized to the isotropic case the result is where s denotes sinh(2K), and where [93] C j (M, For T < T c , let us first recall equation (3.15) of Wu's paper [122], which reduces, for the diagonal correlations C(N, N ), to where α 2 is t 1/2 . Comparing with (1.3) we see that the second term in (2.2) is f (2) N,N = F (2) N,N . Performing the change of variables ξ = z 1 and ξ ′ = 1/z 2 , deforming the contour of integration for both z 1 and z 2 (one has to consider only the discontinuity across the branch cut 3 running from 0 to α 2 ), and rescaling z 1 and z 2 , in, respectively, x 1 = z 1 /α 2 and x 2 = z 2 /α 2 , we obtain: Similarly, when T > T c , the leading term for G (1) N,N is given by equation (2.29) of [122] f (1) N,N = G (1 − t 1/2 z)(1 − t 1/2 z −1 ) 1/2 which, after deforming the contour of integration to the branch cut, and scaling z = t 1/2 x, becomes f (1) N,N (t) = G (1) N,N (t) = where 2 F 1 (a, b; c; z) is the hypergeometric function [37]. When the low temperature expansion of Section 3 of Wu [122] is performed to all orders, we find that (1.1) holds with F (2n) from which, after deformation of integration contours and rescaling, one obtains, for T < T c , the following new integral representation of F Similarly for T > T c the expansion of Section 2 of Wu [122] is performed to all orders and we find that (1.2) holds with F (2n) N,N given by (2.4) and Changing variables and deforming contours, we obtain: N,N (t) = (−1) n t N (2n+1)/2+2n π 2n+1 The form factor expressions are then obtained by expanding the exponentials. Thus we find, for T < T c , that the form factors in (1.5) read and, for T > T c , the odd form factors in (1.6) read where the last product in (2.6) has to be taken to be equal to unity for n = 0, 1. We note that the factors 1/(n!) 2 and 1/(n!(n + 1)!) in (2.5) and (2.6), arise because the integrands are symmetric functions of the variables x 2j and x 2j−1 , separately. This is to be contrasted with (2.1), where there is no separation in the odd and even integrals φ j . In the simplest case the previous integral representation (2.6) gives f (1) N,N (t) defined by (2.3) where one recognizes the Euler representation of an hypergeometric function.
Do note that the (G is not unique. In contrast, the form factor expressions (2.5), (2.6) are unique and well-defined.

Fuchsian linear differential equations for
We use formal computer algebra to study the functions f (j) N,N . We obtain the Fuchsian linear differential equations satisfied by the f (j) N,N for fixed j ≤ 9 and arbitrary N . We also find the truly remarkable result that the f N,N are each solutions of linear differential operators which have a nested "Russian-doll" structure. Beyond this "Russian doll" structure, each linear differential operator is the direct sum of linear differential operators equivalent 4 to symmetric powers of the second order differential operator corresponding to f (1) N,N , (or equivalently to the second order differential operator L E , corresponding to the complete elliptic integral E). A direct consequence is that the form factors f A simple example is f (2) 0,0 = K(K − E)/2. In previous studies on the Ising susceptibility [126,127,128,129], efficient programs were developed which, starting from long series expansions of a holonomic function, produce the linear ordinary differential equation (in this case Fuchsian) satisfied by the function. In order for these programs to be used to study the f  N,N in terms of theta functions of the nome of elliptic functions, presented in [93].
We have obtained the Fuchsian linear differential equations satisfied by the (diagonal) form factors f (j) N,N up to j = 9. The analysis of these linear differential operators shows a remarkable Russian-doll structure similar to the nesting of (the differential operators of) theχ (j) 's found in [126,127,128,129]. Specifically we find that the expressions f (1) N,N , f N,N , f (5) N,N , f (7) N,N are actually solutions of the linear ODE for f (9) N,N , and that f N,N , f (4) N,N , f (6) N,N are actually solutions of the linear ODE for f (8) N,N . In addition, we find that all the linear differential operators for the f (j) N,N 's have a direct sum decomposition in operators equivalent to symmetric powers of the linear differential operator corresponding to f (1) N,N . Consequently, all the f (j) N,N 's can also be written as polynomials in terms of the complete elliptic integrals E and K. The remainder of this section is devoted to the presentation of these results.

Fuchsian linear differential equations for
The linear differential operator F 9 (N ) which annihilates f (9) N,N has the following factorized form where the linear differential operators L r (N ) are of order r. The first two operators read with Dt = d/dt and: The expressions (or forms) of L 6 (N ), L 8 (N ) and L 10 (N ) are given in [25]. The linear differential operators F 2n+1 (N ), which annihilate f (2n+1) N,N for n = 0, . . . , 3, are such that: Thus we see that the linear differential operator for f (2n−1) N,N right divides the linear differential operator for f (2n+1) N,N for n ≤ 3. We conjecture that this property holds for all values of n. We thus have a "Russian-doll" (telescopic) structure of these successive linear differential operators.

Fuchsian linear differential equations for
The linear differential operator F 8 (N ) (corresponding to f (8) N,N ) has the following factorized form where the linear differential operators L r (N ) are of order r. The first two read: The expressions (or forms) of the linear differential operators L 5 (N ), L 7 (N ) and L 9 (N ) are given in [25].
Similarly to (3.3) there is also a Russian-doll (telescopic) structure of these successive linear differential operators: Again, we see that the linear differential operator for f (2n−2) N,N right divides the linear differential operator for f (2n) N,N for n ≤ 4. We conjecture that this property holds for all values of n.

Direct sum structure
Not only do the linear differential operators L j (N ) have a factorized Russian-doll structure, but we have found that they also have a direct sum decomposition when the integer N is fixed.
To illustrate this direct sum decomposition, let us write the corresponding linear differential operator for f where L 2 (N ) is the linear differential operator for f (1) N,N and where the fourth order operator M 4 (N ) is displayed in [25] for successive values of N . One remarks on these successive expressions that the degree of each polynomial occurring in these linear differential operators M 4 (N ) grows linearly with N .
As a further example consider f (5) (N, N ), where we find that the corresponding linear differential operator decomposes as where L 2 (N ) is the linear differential operator for f (1) N,N , M 4 (N ) is the previous fourth order differential operator, and the sixth order operator M 6 (N ) has again coefficients whose degrees grow with N for successive values of N . There is nothing specific to f   N,N 's, n being even or odd. In contrast with the Russian-doll way of writing the linear differential operators for f (n) N,N , the direct sum structure, as a consequence of this growing degree, cannot, for generic N , be written in a closed form as operators with polynomials coefficients in front of the derivatives. This "non-closure" of the direct sum structure will have some consequences when performing the scaling limit of these linear differential operators (see Section 4 below).

Equivalence of various L j (N )'s and M j (N )'s linear differential operators
We find that the symmetric square 5 of L 2 (N ) and the linear differential operator L 3 (N ) are equivalent 6 L 3 (N )U (N ) = V (N )Sym 2 (L 2 (N )) 5 The symmetric j-th power of a second order linear differential operator having two solutions f1 and f2 is the linear differential operator of order j + 1, which has f j 1 , . . . , f j−k 1 f k 2 , . . . , f j 2 as solutions. 6 For the equivalence of linear differential operators, see [108,114,115]. with the following intertwiners: Similarly, with the symmetric cube of L 2 (N ), we have the equivalence with: More generally, all the L m (N )'s are (m − 1)-symmetric-power of L 2 (N ). As a consequence their solutions are (m − 1)-homogeneous polynomials of the two hypergeometric solutions of L 2 (N ).
Similarly, for the linear differential operators occurring in the direct sum, one easily verifies, for every integer N , that, for instance, the M 4 (N )'s are equivalent to the cubic-symmetric-power of L 2 (N ) where, for N = 0, 1, 2: As a further example, one can verify, for every value of the integer N , that the sixth order operator M 6 (N ) is equivalent to the fifth symmetric power of L 2 (N ). The solutions of the linear differential operators M m (N ) are also (m − 1)-homogeneous polynomials of the two hypergeometric solutions of L 2 (N ). As a consequence of this direct sum decomposition, the solutions f (n) (N, N ) are (non-homogeneous) polynomials of the two hypergeometric solutions of L 2 (N ) or, equivalently, f (1) N,N (or the hypergeometric solution of (3.2)) and its first derivative.
The second order linear differential operator L 2 (N ) is equivalent [27] to the second order linear differential operator L E corresponding to the complete elliptic integral of the second kind E. As a consequence of the previously described direct sum decomposition, the f (n) N,N 's can also be written as polynomial expressions of the complete elliptic integral of the second kind E and its first derivative E ′ , or alternatively, E and the complete elliptic integral 7 of the first kind K.
Let us just give here a set of miscellaneous examples of polynomial expressions of various form factors. For f (2) N,N , one has 2f (2) where E and K are given by (3.1). Other examples are given in [25].
Miscellaneous remarks. All these remarkable structures are not restricted to diagonal twopoint correlation functions. We keep on restricting to the isotropic Ising model: for the anisotropic Ising model one has (for the correlations and may have for the form factors) similar but more complicated results involving the complete elliptic integral of the third kind Π (see for instance equation (3.35) in H. Au-Yang and J.H.H. Perk [6], or pp. 23-48 in [2], more recently [121] and for a sketch of how the algebro-differential structures generalize in that anisotropic case [26]).
• Further, one can calculate various j-particle contributions f (j) M,N of the off-diagonal two point correlation functions, and verify, again, that they are, in the isotropic case, also polynomial expressions of the complete elliptic integrals E and K. For instance: where s = sinh(2K).
Other miscellaneous examples of such off-diagonal j-particle contributions are displayed in [25]. In the anisotropic case polynomial expressions of E and K and complete elliptic integral of the third kind Π could take place for j-particle contributions f (j) M,N . The occurrence of elliptic integral of the third kind and not more involved hyperelliptic integrals is still not clear (see after equation (3.20) in [93] the remark on Glasser's unpublished work). This work is still in progress.
• The products of the two-point correlation functions C(N, N ) are also solutions of Fuchsian linear ODE's. As a consequence the equal-time xx-correlations [98] of the free-fermion zero-field XY quantum chain, which are, alternatingly, C(N, N ) 2 and C(N, N )C(N +1, N +1), also satisfy a Fuchsian linear ODE.
• Far beyond, recalling Boel, Kasteleyn and Groeneveld papers [20,21,49] one can see that all the two-point correlation functions of Ising models (not necessarily free-fermion Ising models!) can be expressed as sums, weighted with ± signs, of products of two-point correlation functions. Consequently all the n-point correlation functions of the square Ising model are (simple) polynomial expressions of the complete elliptic integrals E and K and, of course, the npoint correlation functions of the square Ising model are solutions of Fuchsian linear ODE's. For the anisotropic Ising model the n-point correlation functions are solutions of PDE's associated with complete elliptic integrals of the third kind (see [26] for a sketch).
• Recalling the relations (1.3), (1.4) between the two-point correlation functions and the form factors we see that, since the isotropic two-point correlation functions and the form factors are both polynomial expressions of the complete elliptic integral E and K, relations (1.3), (1.4) can be interpreted as an infinite number of quite non trivial identities on the complete elliptic integral E and K, for instance: We have similar identities for the (isotropic) off-diagonal two-point correlations C(M, N ). These linear relations on an infinite number of polynomial expressions of the complete elliptic integrals 8 E and K have to be compared with the infinite number of (non-linear) relations on a finite number of polynomial identities on the complete elliptic integral E and K which correspond to (A.1), the quadratic finite difference relations [82,84,96,97] on the two-point correlation functions displayed in Appendix A.
• At criticality, k = 1, many remarkable and much simpler identities can be obtained, for instance the formula 9 (2.34) in [3] on the next to the diagonal (anisotropic) two-point correlations (see also [121]): where F is the hypergeometric function.

The elliptic representation of Painlevé VI
The results we have underlined in this section, namely the unexpectedly simple and remarkable polynomial expressions for the form factors f (j) N,N , correspond to the fact that the associated linear differential operators are direct sums of operators equivalent to symmetric powers of the second order differential operator L E . We already encountered this central key role played by the linear differential operator L E , or the hypergeometric second order linear differential operator (46) given in [27], in our previous holonomic analysis of the two-point correlation functions of the Ising model [27]. In order to understand the key role played by L E , or equivalently operator L 2 (N ), it is worth recalling (see [78], or for a review [50]) the so-called "elliptic representation" of Painlevé VI. This elliptic representation of Painlevé VI amounts to seeing Painlevé VI as a "deformation" (see equation (33) in [50]) of the hypergeometric linear differential equation associated with the linear differential operator One easily verifies that this linear differential operator is actually equivalent (in the sense of the equivalence of differential operators) with L E , or equivalently L 2 (N ). This deep relation between elliptic curves and Painlevé VI explains the occurrence of Painlevé VI on the Ising model, and on other lattice Yang-Baxter integrable models which are canonically parametrized in term of elliptic functions (like the eight-vertex Baxter model, the RSOS models, see for instance [15]). One can see in Section 6 of [25], other examples of this deep connection between the transcendent solutions of Painlevé VI and the theory of elliptic functions, modular curves and quasi-modular functions. Along this line one should note that other linear differential operators, not straightforwardly linked to L E but more generally to the theory of elliptic functions and modular forms (quasimodular forms . . . ), also emerge in the analysis of the λ-extensions of the two-point correlation functions of the Ising model, for selected 10 values of λ: λ = cos(πm/n). This is detailed in Appendix B. 4 The scaling limit of f The closed (exact) formulae (3.3), (3.4) we obtain for the linear differential operators in these nested "Russian doll" structures, enable us to take the scaling limit of these linear operators. We study this scaling limit in this section and show that the "Russian-doll" structure remains valid. The linear differential operators in that "scaled" nested Russian-doll structure remain equivalent to the symmetric power of a singled-out second order linear differential operator (corresponding to the modified Bessel function). In contrast, in the scaling limit, the direct sum of operators decomposition structure is lost, and we explain why.
The scaling of the f (n) N,N 's amounts, on the functions, and on the corresponding differential operators, to taking the limit N → ∞ and t → 1, keeping the limit x = N (1 − t) finite, or in other words, to performing the change of variables t = 1 − x/N , keeping only the leading term in N . Performing these straightforward calculations, the linear differential operators in t for the f (n) N,N 's where N was a parameter, become linear differential operators in the only scaling variable x.
Calling F scal j the scaling limit of the operator F j (N ) we find for j even that and L scal 9 , L scal 7 are given in [25]. Similarly, for j odd, we have and L scal 10 , L scal 8 , L scal 6 are given in [25]. 10 For generic values of λ, the λ-extension C(M, N ; λ) are not holonomic.
Thus, we see that the scaled operators F scal j have a "Russian-doll" structure straightforwardly inherited from the one for the lattice operators F j (N ).
Consider the linear second-order differential operator corresponding to the modified Bessel function K n (x/2) for n = 0, namely: We recognize, in this linear differential operator, the exact identification with the scaled differential operator F scal 1 = L scal 2 . We find that the symmetric square of the linear differential operator B, and the scaled operator L scal 3 are equivalent: Similarly, the symmetric third power of the linear differential operator B, and the scaled operator L scal 4 are equivalent, and, more generally, the symmetric j-th power of (4.3) and the scaled operator L scal j+1 are equivalent: Recall that the linear differential operators F j (N ), corresponding to the form factors f (j) N,N , can be written as direct sums only when the integer N is fixed. At the scaling limit, this feature disappears for the scaled linear differential operators F scal j which have no direct sums. Therefore while the scaling limit preserves the Russian-doll (telescopic) structure (see (3.3), (4.2)) and also preserves the fact that the various operators in this Russian-doll (telescopic) structure are equivalent to symmetric powers of an operator (4.3) which replaces the operator L E , the direct sum structure is lost. As a consequence the scaling of the f There is one exception that concerns f (2) N,N . Its scaled linear differential operator F scal 2 , has the non shared property of being equivalent to the direct sum of Dx with the symmetric square of (4.3), namely: From this equivalence, one immediately deduces the expression of the scaling of the f (2) N,N as a quadratic expression of the modified Bessel functions of x/2 which actually identifies with formula (2.31b)-(3.151) in [123].
The occurrence of modified Bessel functions, emerging from a confluence of two regular singularities of the complete elliptic integrals E and K, or from the hypergeometric function 2 F 1 , should not be considered as a surprise if one recalls the following limit of the hypergeometric function 2 F 1 yielding confluent hypergeometric functions 1 F 1 . These confluent hypergeometric functions, 1 F 1 , are nothing but modified Bessel functions [37] 2 F 1 a, p, b; Remark. It was shown, in Section 3, as a consequence of the decomposition of their linear differential operators in direct sums of operators equivalent to symmetric powers of L E , that the functions f (n) N,N are polynomial expressions of E and K functions. Therefore their singularities are only the three regular points t = 0, t = 1 and t = ∞. The scaling limit (t = 1 − x/N , t → 1, N → ∞) corresponds to the confluence of the two regular singularities t = 0 and t = ∞, yielding the, now, irregular singularity at x = ∞. The occurrence of irregular singularities with their Stokes phenomenon, and, especially, the loss of a remarkable direct sum structure, shows that the scaling limit is a quite non-trivial limit.
Contrary to the common wisdom, the scaling limit does not correspond to more "fundamental" symmetries and structures (more universal . . . ): this limit actually destroys most of the remarkable structures and symmetries of the lattice models 11 .
5 Bridging with other formula of form factors in the scaling limit: work in progress The Ising Form Factors in the scaling limit as they can be found in Wu, McCoy, Tracy, Barouch [123] read and: The Ising Form Factors in the scaling limit are also given, in many field theory papers, as follows (y j = cosh θ j , see (9) in [7], see also Mussardo [86]): It remains to show that those expressions are actually solutions of the scaled linear differential operators displayed in the previous Section 4, namely (4.1) and (4.2). A straight check yields too large formal calculations. Our strategy should rather be to obtain the series expansions of the n-fold integrals (5.1) or (5.2), in the t variable for (5.1), or the r variable for (5.2), and check that these series expansions are actually solutions of the non-Fuchsian linear differential operators (4.1) and (4.2) of Section 4.
These checks will show that the expressions (5.1), (5.2) of the n-fold integrals of the scaled form factors are actually solutions of linear differential operators with an irregular singularity at infinity and with a remarkable Russian-doll structure but no direct sum structure. This will indicate that such expressions (5.1), (5.2) generalise modified Bessel functions but cannot be simply expressed in terms of polynomial expressions of modified Bessel functions. The interpretation of these expressions in terms of τ -functions (Hirota equations, hierarchies, . . . ) and the link between these Russian-doll structures and Bäcklund transformations or Hirota transformations remains to be done in detail.
6 Other n-fold integrals: from diagonal correlation functions to the susceptibility of the Ising model The study of two-point correlation functions (even n-points correlations . . . ) can be seen as a "warm-up" for the truly challenging problem of the study of the susceptibility of the Ising model and its associated n-fold integrals, the χ (n) (see next Section 7 below). Staying close to the diagonal correlation functions we have introduced a simplification of the susceptibility of the Ising model by considering a magnetic field restricted to one diagonal of the square lattice [28]. For this "diagonal susceptibility" model [28], we benefited from the form factor decomposition of the diagonal two-point correlations C(N, N ), that has been recently presented [25], and subsequently proved by Lyberg and McCoy [74]. The corresponding n-fold integrals χ (n) d were found to exhibit remarkable direct sum structures inherited from the direct sum structures of the form factor [25,28]. The linear differential operators of the form factor [25] being closely linked to the second order differential operator L E (resp. L K ) of the complete elliptic integrals E (resp. K), this "diagonal susceptibility" [28] is also closely linked to the elliptic curves of the two-dimensional Ising model. By way of contrast, we note that the singularities of the linear ODE's for these n-fold integrals [28] are quite elementary (consisting of only n-th roots of unity) in comparison with the singularities we will encounter below with the quite simple integrals (7.4).
Using the form factor expansions (1.5) and (1.6), the λ-extension of this diagonal susceptibility can be written as where the sum is over j even (resp. odd) for T below (resp. above) T c and wherẽ By use of the explicit expressions (2.5), and (2.6), for f (j) N,N we find explicitly, for T < T c , that and, for T > T c , that: We have also found [28], for j = 1, . . . , 4, that theχ d± 's satisfy Fuchsian linear differential equations which have a Russian-doll nesting just as was found for theχ n 's in [126,127,128,129]. In the case of these j-particle components of the "diagonal" susceptibility, we can see that this Russian-doll nesting of the corresponding linear differential operators is straightforwardly inherited [25], not from the Russian-doll nesting of the diagonal form factors f (j) (N, N )'s (this is not sufficient), but from their direct sum (of operators equivalent to symmetric powers) decomposition.

Direct sum decompositions of theχ
are straightforwardly inherited [28] from direct sums of the f (j) N,N thus yielding a scenario for the direct sum decompositions of the "true" χ (n) . However, recalling the non-unicity of the (G Section 2), the direct sum decomposition of the χ (2n) can be seen as a canonical one, when the direct sum decomposition of the χ (2n+1) is not.
7 Other n-fold integrals linked to the susceptibility of the Ising model The susceptibility χ of the square lattice Ising model has been shown by Wu, McCoy, Tracy and Barouch [123] to be expressible as an infinite sum of holomorphic functions, given as multiple integrals, denoted χ (n) , that is kT χ = χ (n) . B. Nickel found [87,88] that each of these χ (n) 's is actually singular on a set of points located on the unit circle |s| = | sinh(2K)| = 1, where K = J/kT is the usual Ising model temperature variable. These singularities are located at solution points of the following equations From now on, we will call these singularities of the "Nickelian type", or simply "Nickelian singularities". The accumulation of this infinite set of singularities of the higher-particle components of χ(s) on the unit circle |s| = 1, leads, in the absence of mutual cancelation, to some consequences regarding the non holonomic (non D-finite) character of the susceptibility, possibly building a natural boundary for the total χ(s). However, it should be noted that new singularities, that are not of the "Nickelian type", were discovered as singularities of the Fuchsian linear differential equation associated [126,128,129] with χ (3) and as singularities of χ (3) itself [29] but seen as a function of s. They correspond to the quadratic polynomial 1 + 3w + 4w 2 where 2w = s/(1 + s 2 ). In contrast with this situation, the Fuchsian linear differential equation, associated [127] with χ (4) , does not provide any new singularities. Some remarkable "Russian-doll" structure, as well as direct sum decompositions, were found for the corresponding linear differential operators for χ (3) and χ (4) . In order to understand the "true nature" of the susceptibility of the square lattice Ising model, it is of fundamental importance to have a better understanding of the singularity structure of the n-particle contributions χ (n) , and also of the mathematical structures associated with these χ (n) , namely the infinite set of (probably Fuchsian) linear differential equations associated with this infinite set of holonomic functions. Finding more Fuchsian linear differential equations having the χ (n) 's as solutions, beyond those already found [126,127] for χ (3) and χ (4) , probably requires the performance of a large set of analytical, mathematical and computer programming "toursde-force".
As an alternative, and in order to bypass this "temporary" obstruction, we have developed, in parallel, a new strategy. We have introduced [29] some single (or multiple) "model" integrals as an "ersatz" for the χ (n) 's as far as the locus of the singularities is concerned. The χ (n) 's are defined by (n − 1)-dimensional integrals [88,94,124] (omitting the prefactor 12 ) 3) The two families of integrals we have considered in [29] are very rough approximations of the integrals (7.2). For the first family 13 , we considered the n-fold integrals corresponding to the product of (the square 14 of the) y i 's, integrated over the whole domain of integration of the φ i (thus getting rid of the factors G (n) and R (n) ). Here, we found a subset of singularities occurring in the χ (n) as well as the quadratic polynomial condition 1 + 3w + 4w 2 = 0.
For the second family, we discarded the factor G (n) and the product of y i 's, and we restricted the domain of integration to the principal diagonal of the angles φ i (φ 1 = φ 2 = · · · = φ n−1 ). These simple integrals (over a single variable), were denoted [29] where x(φ) is given by (7.3).
Remarkably these very simple integrals both reproduce all the singularities, discussed by Nickel [87,88], as well as the quadratic roots of 1 + 3w + 4w 2 = 0 found [126,129] for the linear ODE of χ (3) . One should however note that, in contrast with the χ (n) , no Russian-doll, or direct sum decomposition structure, is found for the linear differential operators corresponding to these simpler integrals Φ We return to the integrals (7.2) where, this time, the natural next step is to consider the following family of n-fold integrals x i (7.5) 12 The prefactor reads (1 − s 4 ) 1/4 /s for T > Tc and (1 − s −4 ) 1/4 for T < Tc and in terms of the w variable. 13 Denoted Y (n) (w) in [29].
14 Surprisingly the integrand with ( Q n j=1 yj) 2 yields second order linear differential equations [29], and consequently, we have been able to totally decipher the corresponding singularity structure. By way of contrast the integrand with the simple product ( Q n j=1 yj ) yields linear differential equations of higher order, but with identical singularities [29].
which amounts to getting rid of the (fermionic) factor (G (n) ) 2 in the n-fold integral (7.2). This family is as close as possible to (7.2), for which we know that finding the corresponding linear differential ODE's is a huge task. The idea here is that the methods and techniques we have developed [126,129] for series expansions calculations of χ (3) and χ (4) , seem to indicate that the quite involved fermionic term (G (n) ) 2 in the integrand of (7.2) should not impact "too much" on the location of singularities of these n-fold integrals (7.2). This is the best simplification of the integrand of (7.2) for which we can expect to retain much exact information about the location of the singularities of the original Ising problem. However, we certainly do not expect to recover from the n-fold integrals (7.5) the local singular behavior (exponents, amplitudes of singularities, etc . . . ). Getting rid of the (fermionic) factor (G (n) ) 2 are we moving away from the elliptic curves of the two-dimensional Ising model? Could it be possible that we lose the strong (Russian-doll, direct sum decomposition) algebro-differential structures of the corresponding linear differential operators inherited from the second order differential operator L E (resp. L K ) of the complete elliptic integrals E (resp. K), but keep some characterization of elliptic curves through more "primitive" (universal) features of these n-fold integral like the location of their singularities?
In the sequel, we give the expressions of Φ H , Φ H and the Fuchsian linear differential equations for Φ (n) H for n = 3 and n = 4. For n = 5, 6, the computation (linear ODE search of a series) becomes much harder. Consequently we use a modulo prime method to obtain the form of the corresponding linear ODE with totally explicit singularity structure. These results provide a large set of "candidate singularities" for the χ (n) . From the resolution of the Landau conditions [29,35] for (7.5), we have shown that the singularities of (the linear ODEs of) these multiple integrals actually reduce to the concatenation of the singularities of (the linear ODEs of) a set of one-dimensional integrals. We discuss the mathematical, as well as physical, interpretation of these new singularities. In particular we can see that they correspond to pinched Landau-like singularities as previously noticed by Nickel [89]. Among all these polynomial singularities, the quadratic numbers 1 + 3w + 4w 2 = 0 are highly selected. We will show that these selected quadratic numbers are related to complex multiplication for the elliptic curves parameterizing the square Ising model.
We present the multidimensional integrals Φ

(n)
H and the singularities of the corresponding linear ODE for n = 3, . . . , 6, that we compare with the singularities obtained from the Landau conditions. We have shown [30] that the set of singularities associated with the ODEs of the multiple integrals Φ (n) H reduce to the singularities of the ODEs associated with a finite number of one-dimensional integrals. Section 9 deals with the complex multiplication for the elliptic curves related to the singularities given by the zeros of the quadratic polynomial 1 + 3w + 4w 2 .  (2 − δ k,0 )(2 − δ p,0 )w n(k+p) a n (k, p), where a(k, p) is a 4 F 3 hypergeometric series dependent on w.
The advantage of using these simplified integrals (7.5) instead of the original ones (7.2) is twofold. Using (8.1) the series generation is straightforward compared to the complexity related to the χ (n) . As an illustration note that on a desk computer, Φ (n) H are generated up to w 200 in less than 10 seconds CPU time for all values of n, while the simplest case of the χ (n) , namely χ (3) , took three minutes to generate the series up to w 200 . This difference between the Φ (n) H and χ (n) increases rapidly with increasing n and increasing number of generated terms. We note that for the Φ (n) H quantities and for a fixed order, the CPU time is decreasing 15 with increasing n. For χ (n) the opposite is the case. The second point is that, for a given n, the linear ODE can be found with less terms in the series compared to the linear ODE for the χ (n) . Indeed for χ (3) , 360 terms were needed while 150 terms were enough for Φ (3) H . The same feature holds for χ (4) and Φ (4) H (185 terms for χ (4) and 56 terms 16 for Φ (4) H ). With the fully integrated sum (8.1), a sufficient number of terms is generated to obtain the linear differential equations. We succeeded in obtaining the linear differential equations, respectively of minimal order five and six, corresponding to Φ H (n ≥ 5), the calculations, in order to get the linear ODEs become really huge 17 . For this reason, we introduce a modular strategy which amounts to generating long series modulo a prime and then deducing the ODE modulo that prime. Note that the ODE of minimal order is not necessarily the simplest one as far as the required number of terms in the series expansion to find the linear ODE is concerned. We have already encountered such a situation [127,28]. For Φ (5) H (resp. Φ (6) H ), the linear ODE of minimal order is of order 17 (resp. 27) and needs 8471 (resp. 9272) terms in the series expansion to be found.
Actually, for Φ H (resp. Φ H ), we have found the corresponding linear ODEs of order 28 (resp. 42) with only 2208 (resp. 1838) terms from which we have deduced the minimal ones. The form of these two minimal order linear ODEs obtained modulo primes is sketched in Appendix C. In particular, the singularities (given by the roots of the head polynomial in front of the highest order derivative), are given with the corresponding multiplicity in Appendix C. Some details about the linear ODE search are also given in Appendix C.
We have also obtained very long series (40000 coefficients) modulo primes for Φ H , but, unfortunately, this has not been sufficient to identify the linear ODE (mod. prime) up to order 100.
The singularities of the linear ODE for the first Φ (n) H are respectively zeros of the following polynomials (besides w = ∞): For n = 7 and n = 8, besides modulo primes series calculations mentioned above, we also generated very long series from which we obtained in floating point form, the polynomials given in Appendix D (using generalised differential Padé methods).
If we compare the singularities for the ODEs for the Φ (n) H to those obtained with the "Diagonal model" 18 presented in [29], i.e. for the ODEs for the Φ (n) D , one sees that the singularities of the linear ODE for the "Diagonal model" are identical to those of the linear ODE of the Φ (n) H for n = 3, 4 (and are a proper subset to those of Φ (n) H for n = 5, 6). The additional singularities for n = 5, 6 are zeros of the polynomials: For n = 7, the zeros of the following polynomials (among others) are singularities which are not of Nickel's type (7.1) and do not occur for Φ (n) The linear ODEs of the multiple integrals Φ We found it remarkable that the linear ODEs for the integrals Φ (n) D display all the "Nickelian singularities" (7.1) , as well as the new quadratic numbers 1 + 3w + 4w 2 = 0 found for χ (3) . It is thus interesting to see how the singularities for Φ   H (resp. Φ (6) H ). In [30] it was shown how this comes about and how it generalizes. For this, we had to solve the Landau conditions [30] for the n-fold integrals (7.5).

Bridging physics and mathematics
In a set of papers [26,27] and in the previous sections, we have underlined the central role played by the elliptic parametrization of the Ising model, in particular the role played by the second order linear differential operator L E (or L K ) corresponding to the complete elliptic integral E (or K), and the occurrence of an infinite number of modular curves [25], canonically associated with elliptic curves. We are getting close to identify the lattice Ising model, (or more generally Baxter model), with the theory of elliptic curves. In such an identification framework one may seek for "special values" of the modulus k that could have a "physical meaning", as well as a "mathematical interpretation" (beyond just being singularities), as singularities of the χ (j) .

Revisiting the theory of elliptic curves with a physics viewpoint
The deep link between the theory of elliptic curves and the theory of modular forms is now well established [105]. More simply the crucial role of the modular group in analysing elliptic curves is well known. For that reason seeking "special values [71]" of the modulus k, that might have a "physical meaning" as well as a mathematical meaning, as singularities of the χ (j) , it may be interesting to, alternatively, introduce the modular function called the j-function which corresponds to Klein's absolute invariant multiplied by (12)  and, alternatively, seek for "special values" of the j-function (9.1), since it automatically takes into account the modular symmetry group of the problem. The modular group requires one to introduce the period ratio and the nome of the elliptic functions. The elliptic nome, defined in terms of the periods of the elliptic functions, reads where τ is the half period ratio 19 .
The SL(2, Z) transformations of the modular group which preserve the j-function (9.1), should not be confused with isogenies of elliptic curves like the Gauss or Landen transformations and, more generally (n integer) which actually modify the j-function (9.1), but are "compatible" with the lattice of periods (the inclusion of one lattice into the other one). Roughly speaking, and as far as the elliptic curves of the Ising model (resp. Baxter model) are concerned, the SL(2, Z) transformations of the modular group are invariance symmetries (reparametrizations), while the transformations (9.4) are highly non-trivial covariants that we will see as exact representations of the renormalization group.

Landen and Gauss transformations as generators of the exact renormalization group
Let us consider the complete elliptic integral K(k) defined as: K(k) = 2 F 1 1/2, 1/2; 1; k .
Two relations between K(k), evaluated at two different modulus, can be found in, e.g. [36] and read 20 : The arguments in K, on the right-hand-side of (9.5), (9.6), are the square of the modulus k transformed by the so-called (descending) Landen or (ascending) Landen (or Gauss) transformations: A sequence of such transformations can be used to evaluate (numerically), in a rapidly convergent way, the elliptic integrals from iterations of (9.7) or of (9.8). Changing k to the complementary modulus k ′ = √ 1 − k 2 , and likewise for the transformed k, the half period ratio transforms through (9.7), (9.8), like (9.3).
The real fixed points of the transformations (9.7) and (9.8) are k = 0 (the trivial infinite or zero temperature points) and k = 1 (the ferromagnetic and antiferromagnetic critical point of the square Ising model). Iterating (9.7) or (9.8), one converges, respectively to k = 0 or k = 1. In terms of the half period ratio, this reads, respectively, τ = ∞ and τ = 0 which correspond to a degeneration of the elliptic parametrization into a rational parametrization. In view of these fixed points, it is natural to identify the transformations (9.7) or (9.8), and, more generally, any transformation 21 τ → nτ or τ → τ /n (n integer), as exact generators of the renormalization group. It is a straightforward exercise, using the identities (9.5), (9.6), to write a "renormalization recursion" on the internal energy U of the Ising model where c and s denote cosh(2K) and sinh(2K) respectively.

Complex multiplication of elliptic curves and fixed points of Landen transformations
Since we are interested in singularities in the complex plane of some "well-suited" variable (s, k, w), one should not restrict (9.7) and (9.8) as transformations on real variables, restricting to real fixed points of these transformations, but actually consider the fixed points of these transformations seen as transformations on complex variables. For instance, if one considers (9.8) as an algebraic transformation of the complex variable k and solve k 2 1 − k 2 = 0, one obtains: The roots of k 2 + 3k + 4 = 0, (9.9) are (up to a sign) fixed points of (9.8). We thus see the occurrence of additional non-trivial complex selected values of the modulus k, beyond the well-known values k = 1, 0, ∞, corresponding to degeneration of the elliptic curve into a rational curve, and physically, to the critical Ising model and to (high-low temperature) trivializations of the model. Of course, when extending (9.8) to complex values, one can be concerned about keeping track of the sign of k 1 in (9.8) in front of the square root √ k. Reference [30] provides a similar fixed point calculation for (9.8) extended to complex values, but for a representation of (9.8) in term of the modular j-function. Such calculations single out the remarkable integer value j = (−15) 3 , which is known to be one of the nine Heegner numbers (see [30]). It is important to note that this representation of (9.8) in term of the modular j-function is the well-known fundamental modular curve symmetric in j and j 1 (see [30,51,75]) j 2 j 2 1 − (j + j 1 )(j 2 + 1487jj 1 + j 2 1 ) − 12 · 30 6 (j + j 1 ) + 3 · 15 3 (16j 2 − 4027jj 1 + 16j 2 1 ) + 8 · 30 9 = 0.
which represents, at the same time, the Landen and Gauss transformations (9.3) as a consequence of the modular invariance (τ ↔ 1/τ ). A straightforward calculation of the elliptic nome (9.2) gives, for the polynomial (9.9) and the polynomial deduced from the Kramers-Wannier duality k → 1/k respectively, exact values for τ , the half period ratio, as very simple quadratic numbers: These quadratic numbers correspond to complex multiplication and to j = (−15) 3 . These two quadratic numbers are such that 2τ 1 ∓ 1 = τ 2 . Let us focus on τ 2 for which we can we write Taking into account the two modular group involutions τ → 1 − τ and τ → 1/τ , we find that 1 − 2/τ is, up to the modular group, equivalent to τ /2. The quadratic relation τ 2 − τ + 2 = 0 thus amounts to looking at the fixed points of the Landen transformation τ → 2τ up to the modular group. This is, in fact a quite general statement: the complex multiplication values can all be seen as fixed points, up to the modular group, of the generalizations of Landen transformation, namely τ → nτ for n integer (here ≃ denotes the equivalence up to the modular group): Complex multiplication corresponds to integer values of the modular j-function (as in the case of the Heegner numbers see [30]). For elliptic curves in field of characteristic zero, the only well-known selected set of values for k corresponds (besides k = 0, 1, ∞) to the values for which the elliptic curve has complex multiplication [77], and we see these selected values, here, as fixed points, in the complex plane, of transformations (isogenies) that are exact representations of generators of the renormalization group.
It is now totally natural to see if the singularities we have obtained for the n-fold integrals (7.5), can be interpreted in the framework of elliptic curve theory, in terms of this physically, and mathematically, highly selected set of values for elliptic curves, namely complex multiplication values.

Complex multiplication for 1 + 3w + 4w 2 = 0
Let us consider the first unexpected singularities 1 + 3w + 4w 2 = 0 we found [126,129] for the Fuchsian linear differential equation of χ (3) , and also found in other n-fold integrals of the Ising class [29]. This polynomial condition reads in the s variable, 2s 2 + s + 1 s 2 + s + 2 = 0. We have shown [29] that χ (3) itself is not singular at the roots of the first polynomial whose roots are such that |s| < 1, but is actually singular at the roots of the second polynomial. In the variable k = s 2 , these singularities read: (4k 2 + 3k + 1)(k 2 + 3k + 4) = 0. (9.10) The second polynomial has actually been seen to correspond to fixed points of the Landen transformation (see (9.9)). Note that the two polynomials in (9.10) are related by the Kramers-Wannier duality k → 1/k (and therefore both correspond to the same value of the modular j-function: j = (−15) 3 ).
In other words we see that the selected values 1 + 3w + 4w 2 = 0, occurring in the (hightemperature) susceptibility of the Ising model as singularities of the three-particle term χ (3) , actually correspond to the occurrence of complex multiplication on the elliptic curves of the Ising model, and can also be seen as fixed points of the renormalization group when extended to complex values of the modulus k.
Let us note that the occurrence of Heegner numbers and complex multiplication has already occurred in other contexts, even if the statement was not explicit. In the framework of the construction of Liouville field theory, Gervais and Neveu suggested [46] new classes of critical statistical models (see Appendix E), where, besides the well-known N -th root of unity situation, they found the following selected values of the multiplicative crossing [103] t: If one wants to see this multiplicative crossing [40,52,73,80] as a modular nome (see [30]), the two previous situations actually correspond to selected values of the modular j-function namely j((1 + i √ 3)/2) = (0) 3 for (9.11), and j(1 + i) = (12) 3 for (9.12), which actually correspond to Heegner numbers and complex multiplication [77]. It is however important not to feed the confusion already too prevalent in the literature, between a "temperature-like" nome like (9.2) and a multiplicative crossing modular nome (see Appendix E). In the Baxter model [11], the first is denoted by q and the second one by x. In fact one probably has, not one, but two modular groups taking place, one acting on the "temperature-like" nome q and the other one acting on the multiplicative crossing x. We will not go further along this quite speculative line which amounts to introducing elliptic quantum groups [79] and (see Appendix E) elliptic gamma functions (generalization of theta functions 22 ) which can be seen [39] as "automorphic forms of degree 1", when the Jacobi modular forms are "automorphic forms of degree 0" and are associated (up to simple semi-direct products) to SL(3, Z) instead of SL(2, Z). 9.5 Beyond 1 + 3w + 4w 2 = 0 As a consequence of the fact that the modular j-function is a function of w 2 , the quadratic polynomial condition 1 − 3w + 4w 2 = 0, corresponds to the same selected values of the modular j-function as 1 + 3w + 4w 2 = 0, namely j = (−15) 3 . The quadratic polynomial 1 − 3w + 4w 2 = 0 actually occurs in the singularities of the linear ODE for Φ (6) H (and all the higher Φ (2n) H , if one believes formulas (28) and (29) in [30]).
In view of the remarkable mathematical (and physical) interpretation of the quadratic values 1+3w+4w 2 = 0, (and also 1−3w+4w 2 = 0) in terms of complex multiplication, or fixed points of the renormalization group, it is natural to see if such a "complex multiplication" interpretation also exists for other singularities of χ (n) , and as a first step, for the singularities of the linear differential equations of our n-fold integrals (7.5), that we expect to identify, or at least, have some overlap with the singularities of the χ (n) .
Among the singularities of the linear ODE for Φ (n) H given in (8.2), (8.3) or obtained from the formula (29) given in [30] up to n = 15, we have found no other singularity identified with the remarkable Heegner numbers [100] or, more generally, with other selected values of the modular j-function, associated to complex multiplication.
Could it be that the (non-Nickelian) singularities (8.2), (8.3), which do not match with complex multiplication of the elliptic curves, are actually selected for mathematical structures more complex or more general than elliptic curves (possibly linked [39] to SL(3, Z) instead of SL(2, Z) modular group)? This could amount to moving away from the isotropic Ising model towards the Baxter model. At first sight the analysis of the anisotropic Ising model [26] could be considered as a first step in that "Baxter-model" direction. The selected situations for elliptic functions and complete elliptic integrals, would thus, be generalized to the search of "selected situations" of their multidimensional generalizations (Lauricella, Kampé de Ferié, Appell, . . . ) that we have actually seen to occur in the anisotropic Ising model [26] and even in our series expansions of χ (3) and χ (4) .
Along similar lines, one may recall the n-fold integrals introduced by Beukers, Vasilyev [116,117] and Sorokin [109,110] and other well-poised hypergeometric functions or the Goncharov-Manin integrals [43] which occur in the moduli space of curves [19,48]. These integrals [31,41,42,55,56,69,101] look almost the same as the ones we have introduced and analyzed in the study of the diagonal susceptibility of the Ising model [28].
It is worthy to recall that ζ(3) appeared in some of our "connection matrix method" results for the differential Galois group [128] of the Fuchsian linear ODE for χ (3) and χ (4) , and the occurrence of zeta functions in many n-fold integrals. Also recall that Feynman amplitudes can be seen as periods in the "motivic sense" [19], and are often linked to multiple zeta numbers. Along this line, the following integral [41,43] deals with ζ(3): From the series expansion of this holonomic n-fold integral, we have obtained [30] an order four Fuchsian linear differential equation (see Appendix F). On such linear differential operators the "logarithmic" nature of these integrals becomes clear. The occurrence of linear differential operators is not a complete surprise if one recalls that in Apéry's proof of the irrationality of ζ(3) a crucial role is played by the linear differential operator [18] this operator being linked to the modularity of the algebraic variety: x These n-fold integrals have to be compared with the (more involved) n-fold Ising integrals corresponding to the χ (n) , and to the theory of elliptic curves (rather than rational curves CP 1 in the previously cited examples [31,34,41,42,55,56,69,101]), we try to underline in this paper.
With these new singularities, are we exploring some remarkable "selected situations" of some moduli space of curves [1,32] corresponding to pointed [9,38,54] (marked) curves [10], instead of simple elliptic curves [53]? In practice this will probably just correspond to considering a product of n times a rational or elliptic curve minus some sets of remarkable algebraic varieties [28],

Conclusion
We have displayed several examples of n-fold holonomic integrals associated with the twodimensional Ising model on a square lattice [8]. The corresponding linear differential operators with polynomial coefficients are shown to be very closely linked to the theory of elliptic curves (and modular forms) and display many remarkable structures (Russian-doll structure, direct sum structure, complex multiplication as selected singular values for these operators, . . . ). These linear differential operators are not only Fuchsian operators, they are Fuchsian operators with rational exponents: the various indicial polynomials corresponding to all the regular singularities of these linear differential operators have only rational (or integer) roots. It is tempting to try to understand these deep algebraico-differential structures as a consequence of the underlying elliptic curve in the Ising model, or more generally, of some algebraic varieties built from this elliptic curve (product of curves, . . . ), or corresponding to the integrands of these n-fold integrals. Could it be possible that these large number of remarkable properties have a geometrical interpretation (generalisation of hypergeometric functions and Picard-Fuchs systems, Griffiths cohomology of hypersurface of CP n , rigid local systems [17,45,62,63,64,65,66,99], . . . ) with a strong background of algebraic geometry? One could, for instance, imagine that these various n-fold holonomic integrals might be interpreted as periods of some algebraic varieties, all the strong and deep algebraico-differential structures we have displayed in this paper, being a consequence of this very rigid geometrical framework. The central role played by the theory of elliptic curves and their isomonodromic deformations (Painlevé equations) for the Ising model on a lattice is also underlined in the fundamental finite-difference (non-linear overdetermined) system of quadratic functional relations [82,84,96,97] (see (A.1) in Appendix A) for the two-point correlation functions of the Ising model on the square lattice. As Painlevé and (discrete) integrability specialists call it, these lattice equations are finite-difference generalisation of Painlevé equations and they have a lot of very deep consequences : they are, for instance, the very reason why the susceptibility series can be calculated from a program with polynomial growth [93]. Such an overdetermined system (A.1) can be seen as generating an infinite number of non-trivial identities on the complete elliptic integrals of the first and second kind.
It is important to note that all these remarkable structures and deep symmetries (remarkable functional identities, algebraico-differential structures, modular forms, continuous [81,90] and discrete Painlevé structures, . . . ), underline the central role played by the theory of elliptic curves for the two-dimensional Ising model on a lattice. Note that a large part of these remarkable structures and deep (lattice) symmetries is lost in the scaling limit. In the scaling limit some of these remarkable structures remain (the Russian-doll telescopic embedding of the linear differential operators), but, for instance, the direct-sum structure is lost. The scaling limit yields the occurrence of an irregular singularity at infinity: the Fuchsian character of the linear operators is lost, as well as most of the remarkable structures associated with the underlying elliptic curve theory. For instance for two-point correlation functions, the complete elliptic integrals of first and second kind K and E are replaced by modified Bessel functions (with their irregular singularity at infinity), but the fact that form factors are simple polynomial expressions of E and K is lost: the form factors, in the scaling limit, are not simple polynomial expressions of modified Bessel functions. In the scaling limit, a large part of the strong background of algebraic geometry that exists on the lattice model, and yields so many remarkable deep and strong structures and symmetries, seems to disappear. If the geometrical interpretation we suggested for the lattice model exist, could it be possible that it is essentially lost in the scaling limit, the underlying algebraic varieties necessary for this geometrical interpretation being lost, or becoming some complicated analytical manifolds? Recalling the emergence of an irregular singularity (at infinity), an irregular singularity can, in principle, be understood [44] as a confluence of two regular singularities (for complete elliptic integrals of first and second kind we have the confluence of the two regular singularities 0, ∞ among the three regular singularities 0, 1, ∞). To our knowledge we have not often seen 23 in the litterature the structures associated to irregular singularities (Stokes multipliers, singular behaviours, . . . ) be obtained as a "confluent limit" of the structures associated with the two regular singularities. From a general viewpoint, in a desire to see analytical manifolds as a confluent limit of algebraic varieties, one can imagine that the structures of the Ising model in the scaling limit could, in principle, be obtained from a (very involved) "confluent limit" of the remarkable structures deeply linked to the theory of elliptic curves that exist for the Ising model specifically on the lattice. This remains to be done. In practice we see that, paradoxically from a criticality-universality mainstream viewpoint, the (off-critical, non-universal) Ising model on a lattice has much deeper, and fundamental, structures than the Ising model in the scaling limit. Note that the results of holonomic and algebraic-geometry nature we have displayed in this paper, are not specific of the two-dimensional Ising model, or, even, of free-fermion models. We have not used the free-fermion character of the Ising model. We have heavily used the elliptic parametrisation of the two-dimensional lattice Ising model. One can imagine that many of these results, and structures, exist for Yang-Baxter integrable models with an elliptic parametrisation (the Baxter model [22,23,24], . . . ), and, more generally, for any Yang-Baxter integrable model 24 , the central role of the elliptic curve being replaced by the relevance of the algebraic variety emerging in the Yang-Baxter equations (higher genus curves [5], Abelian varieties, . . . ). Integrable models on a lattice are probably deeper, and dressed with much more symmetries and remarkable structures 25 , than their scaling limits. Such an apparently paradoxical (for the field theory mainstream) conclusion is certainly not a surprise for Painlevé and (discrete) integrability specialists who are used to see, and understand, lattice equations as deeper, and more fundamental [85,106], than the differential equations.

A Quadratic partial difference Painlevé generalisations
Quadratic partial difference equations were shown [82,84,96,97] to be satisfied by two-point correlation functions of the two-dimensional Ising model on the square lattice. These quadratic partial difference equations (valid in the anisotropic case), are actually valid for the λ-extension of the two-point correlation functions C(M, N ; λ) for any value of λ : N,N of Sections 2-4, and the corresponding remarkable differential structures, may be used to obtain many further results. We displayed some of these results in Section 6 of [25]. Recalling that, when λ = 1, the Ising correlation functions C(N, N ; 1) satisfy Fuchsian linear differential equations [27] with an order that grows with N , it is quite natural to inquire whether there are any other values of λ for which C(N, N ; λ) will satisfy a Fuchsian linear differential equation. One such family of λ is motivated by the work of Cecotti and Vafa [33] on N = 2 supersymmetric field theories where they encountered λ extensions of the Ising correlations in the scaling limit [83] with (m and n are integers): Indeed, we have found that for n = 3, . . . , 20, the functions C(N, N ; λ) satisfy Fuchsian linear differential equations whose orders, in contrast with those of the λ = 1 equations [27], do not depend on N .
The function C(N, N ; λ) is such that its log-derivative is actually a solution of the sigma form of Painlevé VI : it is a transcendent function "par excellence". However, the unexpectedly simple expressions for these form factors f Actually these special values (B.2) of λ already occurred in a study of N = 2 supersymmetric field theories [33] in a similar series construction of solutions of the Painlevé V (or Painlevé III for a ratio of functions) equation for the scaling limit of the Ising model [123].
Recalling the quadratic finite difference equations [3,84] (A.1) we can deduce that the offdiagonal terms C(M, N ; λ) are, in the isotropic case, algebraic expressions of sum of ratios of theta functions and their derivatives. For the singled-out values λ = cos(πm/n) (n and m integers), the off-diagonal terms C(M, N ; λ) are, in the isotropic case, algebraic expressions of the variable t: do these algebraic expressions also correspond to modular curves? Actually they clearly single out t = 0, 1, ∞, . . . .

Remark.
For this set of selected values of λ the λ-extension C(N, N ; λ) are seen to be algebraic expressions of the variable t and, more remarkably, associated with a modular curve P (C, t) = 0 (P denotes a polynomial with integer coefficients, C denotes C(N, N ; λ) for λ = cos(πm/n), for certain integer values of n and m, and the only branching points are for t = 0, 1, ∞). The fact that the only singular points are t = 0, 1, ∞ can be seen to be inherited from the fact that the λ-extension C(N, N ; λ) is actually solution of (1.7) for any λ: the sigma form of Painlevé VI, namely (1.7), naturally singles out t = 0, 1, ∞ (and only t = 0, 1, ∞).
C Linear differential equations of some Φ The minimal order linear differential equation satisfied by Φ and where the apparent singularities polynomial P 5 (w) in the head polynomial reads the other polynomials P n (w) are given in [30].

C.2 Linear ODE for Φ (4) H
The minimal order linear differential equation satisfied by Φ where the apparent singularities polynomial P 6 (x) in the head polynomial reads: The other polynomials P n (w) are given in [30].
C.3 Linear ODE modulo a prime for Φ (5) H The linear differential equation of minimal order seventeen satisfied by Φ (5) H is of the form 17 n=0 a n (w) d n dw n F (w) = 0, with a 17 (w) = (1 − 4w) 12 (1 + 4w) 9 (1 − w) 2 (w + 1) (1 + 2w) 1 + 3w + 4w 2 2 a 14 (w) = w 9 (1 − 4w) 9 (1 + 4w) 6 P 14 (w), · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · where the 430 roots of P 17 (w) are apparent singularities. The degrees of these polynomials P n (w) are such that the degrees of a i (w) are decreasing as: deg(a i+1 (w)) = deg(a i (w)) + 1. In fact, with 2208 terms we have found the ODE of Φ ((5) H at order q = 28 using the following ansatz for the linear ODE search (Dw denotes d/dw) where α(n) = min(0, n) and: the p(i) being the unknown polynomials. The minimal order ODE is deduced from the set of linearly independent ODEs found at order 28. Instead of these linear ODEs with quite large apparent singularities polynomials, we can provide an alternative linear ODE of higher order with no apparent singularities. This is the so-called "desingularization" procedure of a linear ODE. The price to pay to get rid of the large apparent polynomial can be that the higher order ODE with no apparent polynomial may not be Fuchsian anymore (because of an irregular singularity at infinity). One can also consider desingularizations preserving Fuchsianity.
Remark. We sketch such a quite tedious result (if we give explicitly the undefined polynomials the result would be really huge . . . ) to give the reader some hint of how such an exact result modulo a prime looks like: the exact expressions of the various polynomials, which are the coefficients in front of the derivatives, can actually be factorized modulo prime without any ambiguity. For instance the factor (w + 27448) 2 in the head polynomial (coefficient of Dw 25 ) is nothing but (w − 1) 2 modulo the prime 27449. The interest of such an exact calculation is that we can exactly compare the various factors in the head polynomial with a set of polynomials we have conjectured to be singularities of the linear ODE. We can totally confirm the existence of some of (or all) these conjectured polynomials, and discriminate between apparent singularities and "true" singularities. The prime is large enough to avoid any ambiguity corresponding to accidental factorisations (because the prime would be too small): this is confirmed by the same calculations performed for other similar large enough primes.

C.4 Linear ODE modulo a prime for Φ (6) H
The linear differential equation of minimal order (namely twenty-seven), satisfied by Φ a 24 (x) = (1 − 16x) 13 x 18 P 24 (x), · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · where the 307 roots of P 27 (x) are apparent singularities. The degrees of the P n (w) polynomials are such that the degrees of a i (w) are decreasing as: deg(a i+1 (w)) = deg(a i (w)) + 1 In fact, with 1838 terms we have found the linear ODE of Φ H at order q = 42 using the following ansatz for the linear ODE search (Dx denotes d/dx) where α(n) = min(0, n) and the p(i) being the unknown polynomials. The minimal order ODE is deduced from the set of linearly independent ODEs found at order 42.
Here also, instead of this linear ODE with quite a large apparent singularities polynomial, we can provide an alternative linear ODE of higher order with no apparent singularities (but it may not be Fuchsian anymore). We give in the following the linear differential operator, modulo the prime 32749, of order 30. At this order, the linear differential operator has no apparent singularities (x + 8187) 3 x 2 + 10234x + 22515 (x + 10234) 16  where the P n and Q n a "short" notation for polynomials of degree n (that may be different from one order D m x to another). The factor (x + 32748) in the head polynomial (coefficient of Dx 30 ) is nothing but the factor (x − 1) modulo the prime 32749. D Singularities in the linear ODE for Φ (7) H and Φ H , we generated long series, unfortunately, insufficient to obtain the corresponding linear ODE. Actually, we have also generated very long series modulo a prime (40000 coefficients) and we have not been able to find a linear ODE when the order of the ODE is less than 100. However, by steadily increasing the order q of the ODE and the degrees n of the polynomials in front of the derivatives, one may recognize, in floating point form, the singularities of the linear ODE as the roots of the polynomial in front of the higher derivative. A root is considered as singularity of the still unknown linear ODE, when as q and n increase, it persists with more stabilized digits.
Note that the stabilized digits in these singularities can be as low as two digits.

E Selected values for Liouville theory and Potts models
New classes of critical statistical models where suggested [46] by Gervais and Neveu from the construction of Liouville field theory. With the Q-state standard scalar Potts model notations (see (1.3) in [46]), they introduced y, such that Q 1/2 = 2 cos(πy/2). Rational values of y correspond to selected values of Q (Tutte-Beraha numbers see Section 4 of [104]) for which the standard scalar Potts model has rational critical exponents. At this step, and in order to make explicit the selected role of these particular values, we can recall the expression (see (3.3) in [103]) of the partition function per site of the Q-state standard scalar Potts model on the checkerboard lattice in terms of Eulerian products (see (3.5) in [103]) like (with the notations of [103]): This Eulerian product form made very clear the fact that the partition function can be seen as some automorphic function with respect to an infinite discrete group generated by the inverse relation and the symmetries of square [57,58]. Such Eulerian product over an infinite discrete group also made very clear the fact that these singled-out values of Q actually correspond 26 to N -th root of unity situation that occur in some many domains of theoretical physics [111,113] (dilogarithms, Kac determinant, . . . ). Do note that such situation generalizes, mutatis mutandis, to the Baxter model: the partition function per site can actually be written as an infinite discrete product [11,13,14] over a group generated by the inverse relation and geometrical symmetries of lattice [76], expressions like (E.1) being replaced by (with Baxter's notations [11,13,14]) where q = exp(−πI ′ /I), x = exp(−πλ/2I), z = exp(−πv/2I).
Such an expression of the partition function per site of the Baxter model as infinite product can also be found in [39] in terms of product and ratio of theta and elliptic gamma functions.
In [46] Gervais and Neveu underlined that they had built Liouville field theory for other singled-out values of Q than N -th root of unity situations like (E.2), namely (see (2.3)  that is t 2 = −e −π √ 3 and t 2 = −e −2π respectively. Actually the variable t in (E.1) or in [103] is exactly what is called the multiplicative crossing in conformal theory [40,52,73,80]. Conformal field theoreticians are keen on introducing modular group structure for which the multiplicative crossing is seen as a modular nome q. If we follow this line recalling the relation q = exp(iπτ ) between the nome and the half period ratio τ (see [100]), we find that the two previous situations actually correspond to singled-out values of the modular j-function namely j((1+i √ 3)/2) = (0) 3 26 Note that this Q, corresponding to the number of state of the Potts model, should not be confused with a nome q. It was unfortunately denoted q in [46]. 27 One has to be careful with the various notations in the literature where, as far as nomes are concerned, one moves from q to q 2 . In [100] the nomeq corresponds to t 2 . Relation (9.12) readsq = q 2 = e 2iπτ = −e −π √ 3 .
for (E.3), and j(1 + i) = (12) 3 for (E.4), which actually correspond to Heegner numbers and complex multiplication [77,100]. Considering λ-extensions of two-point diagonal correlation function of the Ising model, we found [25] modular curves corresponding to polynomial relations between a (modular) function and its first derivative, this (modular) function being a very simple ratio of Jacobi theta functions (see Section 6.1 in [25]). Along this line it is worth recalling the "special value" −e −π √ 3 of the nome of Jacobi theta functions (at zero argument) for which a ratio of Jacobi theta functions becomes a simple algebraic expression [120] At this step it is fundamental to raise an important confusion that overwhelms the theoretical physics literature. In many domains of theoretical physics the existence of a modular group and/or N -th root of unity situations in some "nome" always denoted q, is underlined and analyzed. In Liouville theory this nome q is the exponential 28 of , in conformal field theory 29 two q's, and two modular group structures, can be introduced, the second one corresponding to finite size analysis with the introduction of a modular parameter for partition function on a (finite size l × l ′ ) torus (see for instance (3.33) in [95]).
Sticking with Baxter's notations the complex multiplication situation, we see in this paper with selected values like 1 + 3w + 4w 2 = 0, corresponds to selected values of the modulus of the elliptic curves, or of the nome q which measures the distance to criticality (temperature-like variable) of the off-critical lattice model. In contrast the selected values (B.1) of λ (for which modular curves are seen to occur for the λ-extensions of the correlation functions) correspond to N -th root of unity situations for the multiplicative crossing x. Most of the field theory papers (QFT, CFT, . . . ) where selected values (N -th root of unity situations) occur correspond to models at criticality: for these models there is no (temperature-like off-critical) variable like our previous nome q (the elliptic curve is gone, being replaced by a rational curve). All the selected situations encountered are in the multiplicative crossing variable x within a rational parametrization of the model.

F Factorisations of multiple integrals linked to ζ(3)
From the series expansion of the triple integral (9.13) we have obtained the corresponding order four Fuchsian linear differential equation (Dx denotes d/dx) L n = Dx 4 + 2(3x − 1) (x − 1)x Dx 3 + 7x 2 + (n 2 + n − 5)x − 2n(n + 1) (x − 1) 2 x 2 Dx 2 + x 2 + 2n(n + 1) (x − 1) 2 x 3 Dx + n(n + 1) (n 2 + n + 1)x + (n − 1)(n + 2) (x − 1) 2 x 4 which has the following factorization in order-one differential operator: 28 Not to be confused with the q of the q-state Potts model in the paper that cope with Liouville theory and Potts model in the same time! 29 They are, of course, many other occurrences of modular groups and/or occurrences of a nome q (quantum dilogarithms, q-deformation theories, q-difference equations, q-Painlevé, q-analogues of hypergeometric functions, . . . ). The confusion is increased with the dilute AL models and their relations with the Ising Model in a Field for which the corresponding nome q could be associated with the magnetic field of the Ising model [16,107,118,119].
Such factorization in order-one differential operator having rational solutions is characteristic of the strong geometrical interpretation we are seeking for (interpretation of n-fold integrals as periods of some algebraic variety) for the Fuchsian linear differential operators we have obtained for many n-fold integrals (of the "Ising class" [8]). Such a factorization in order-one linear differential operator having rational solutions does not seem to take place in general for our Fuchsian linear differential operators, but seems actually to occur modulo many primes for the Fuchsian linear differential operators of the χ (n) . Such calculations, mixing geometrical interpretation and "modular" calculations on our n-fold integrals, remain to be done.