Why are all dualities conformal? Theory and practical consequences

We relate duality mappings to the"Babbage equation"F(F(z)) = z, with F a map linking weak- to strong-coupling theories. Under fairly general conditions F may only be a specific conformal transformation of the fractional linear type. This deep general result has enormous practical consequences. For example, one can establish that weak- and strong- coupling series expansions of arbitrarily large finite size systems are trivially related, i.e., after generating one of those series the other is automatically determined through a set of linear constraints between the series coefficients. This latter relation partially solve or, equivalently, localize the computational complexity of evaluating the series expansion to a simple fraction of those coefficients. As a bonus, those relations also encode non-trivial equalities between different geometric constructions in general dimensions, and connect derived coefficients to polytope volumes. We illustrate our findings by examining various models including, but not limited to, ferromagnetic and spin-glass Ising, and Ising gauge type theories on hypercubic lattices in 1


I. INTRODUCTION
The utility of weak and strong coupling expansions and of dualities in nearly all branches of physics can hardly be overestimated. This article is devoted to two interrelated fundamental questions: (1) What information does the existence of finite order complementary weak and strong coupling series expansion of given physical quantities (e.g., partition functions, matrix elements, etc.) provide? (2) To what extent can dualities be employed to partially solve those various problems?
Consider a theory of (dimensionless) coupling strength g for which weak and strong coupling expansions may, respectively, be performed in powers of g and 1/g or in other monotonically increasing/decreasing functions f +/− (g). Common wisdom asserts that as ordinary expansion parameters (e.g., g and 1/g) behave very differently, weak and strong coupling series cannot, generally, be simply compared. On a deeper level, if these expansions describe different phases (as they generally do) then the series must become non-analytic (in the thermodynamic limit) at finite values of g (where transitions occur) and thus render any equality between them void. A duality may offer insightful information on a strong coupling theory by relating it to a system at weak coupling that may be perturbatively examined. As is well known, when they are present, selfdualities are manifest as an equivalence of the coefficients in the two different series; this leads to an invariance under an inversion that is qualitatively (and in standard field theories, e.g., QED/Electroweak/QCD is exactly) of the canonical form "g ↔ 1/g" (or, more generally, f + (g) ↔ f − (g)). In vacuum QED, with Lagrangian density L = [ 0 E 2 /2 − B 2 /(2µ 0 )], the ratio g = 0 µ 0 of the couplings in front of the E 2 and B 2 terms relates to a g ↔ 1/g reciprocity. This reciprocity is evident from the invariance of Maxwell's equations in vacuum under the exchange of magnetic and electric fields 1 , E → B; B → − E and the Lagrangian density that results. In Yang-Mills (YM) theories, such an exchange between dual fields has led to profound insights from analogies between the Meissner effect and the behavior of vortices in superconductors to confinement and flux tubes -a hallmark of QCD [2][3][4][5] . Abstractions of dualities in electromagnetism and in YM theories produced powerful tools such as those in Hodge and Donaldson theories 6 .
In both classical and quantum models, dualities (and the f + (g) ↔ f − (g) inversion) are generated by linear transformations (appearing, e.g., as unitary transformations or more general isometries relating one local theory to another in fundamental "bondalgebraic" [7][8][9][10][11][12][13] incarnations or, in the standard case, Fourier transformations [14][15][16][17][18] ). Such linear transformations lead to an effective inversion of the coupling constant g. Dual models share, for instance, their partition functions (and thus the same series expansion). As realized by Kramers and Wannier (KW) [19][20][21][22][23][24][25] , self-dualities provide structure that enables additional information allowing, for instance, the exact computation of transition points. This does not imply that the full partition function is determined with complexity polynomial in the size of the system, it is solvable, via self-dualities alone (and indeed as we illustrate in this work, self-dualities do not suffice).
The principal message of this work is that these linear duality transformations allow a broad investigation of questions (1) and (2) above. Specifically, we will examine arbitrarily large yet finite size systems for which no phase transitions appear and consequently (1') Equate the weak and strong coupling expansions to find linear constraints on the expansion coefficients, and (2') When possible, invoke self-duality to obtain yet further linear equations that those coefficients need to satisfy. This analysis will lead to the concept of partial solvability, that is the ability to compute a specific physical quantity with complexity polynomial in the size of the system, given partial information that is determined by other means. To illustrate, we consider nearest neighbor Ising (H = − ab J ab s a s b , s a = ±1) and various other models on hypercubic lattices Λ of N = L D sites in D dimensions (with even length L), endowed with periodic boundary conditions. Unless stated otherwise, we will focus on uniform ferromagnetic systems (J ab = J > 0 for all lattice links ab ). In 26 we consider other boundary conditions, system sizes and lattice aspect ratios, and show that our results are essentially unchanged for large systems with random J ab = ±J.

II. SERIES EXPANSIONS OF ISING MODELS
To demonstrate our concept, we use standard expansions [22][23][24][27][28][29][30] , where the coupling constant is (g ≡)K ≡ βJ with β the inverse temperature. Defin- leads to a canonical loop structure for the expansion, where C l is the number of (not necessarily connected) loops of total perimeter l = 2l (l = 1, 2, · · · ) that can be drawn on the lattice and C 0 = 1. For each such loop, i.e., Γ = (ab) · · · (mn) formed by the bonds (nearest neighbor pair products {(s a s b )}) appearing in Eq. (2), there is a complementary loop Γ = Λ − Γ for which the sum of Eq. (2) remains unchanged. Consequently, the H-T series coefficients are trivially symmetric, C DN −l = C l . We next briefly review the low temperature (L-T), or strong coupling, expansion. There are two degenerate ground states (with s a = +1 for all sites a or s a = −1) of energy E 0 = −JDN . All excited states can be obtained by drawing closed surfaces marking domain wall boundaries. The domain walls have a total (D − 1) dimensional surface area s , the energy of which is E = E 0 + 2s J. Taking into account the two-fold degeneracy, the L-T expansion of the partition function in powers of (f − (K) ≡)e −2K is with C s the number of (not necessarily connected) closed surfaces of total area s = 2l (C 0 = 1). That is, the L-T expansion is in terms of (D − 1)-dimensional "surface areas" enclosing D-dimensional droplets. Geometrically, there are no closed surfaces of too low areas s . Thus, in the L-T expansion of Ising ferromagnets, The L-T coefficients exhibit a trivial complementarity symmetry akin to that in the H-T series. Given any spin configuration {s a }, there is a unique correspondence with a staggered spin configuration s a = (−1) D α=1 aα s a where a α are the (integer) Cartesian components of the hypercubic lattice site a (i.e., a = (a 1 , a 2 , · · · , a D )). Domain walls associated with such staggered configuration are inverted relative to those in the original spin configuration s a . That is, if a particular domain wall appears in s a then it will not appear in s a and vice versa. As a result, C DN −s = C s (for the even L hypercubic lattices that we consider).

III. EQUATING WEAK AND STRONG COUPLING SERIES
Our approach is to compare H-T and L-T series expansions of arbitrary models by means of a duality transformation. In the Ising model, the mapping relates expansions inT to those in e −2K . In the complex planes of either of the expansion parameters f ± (K) (i.e.,T or e −2K ), Eqs. (5) constitute fractional linear transformations.T is the magnetization of a single Ising spin immersed in an external magnetic field of strength h = K/β when there is a minimal coupling (a Zeeman coupling) between the dual fields: the Ising spin and the external field. This transformation may be applied to Ising models in all dimensions D -not only to the D = 2 model for which the KW correspondence holds. These transformations emulate, yet are importantly different from, a g ↔ 1/g correspondence (the latter never enables an equality of two finite order polynomials in the respective expansion parameters). Employing the second of Eqs. (5), By virtue of Eq. (2), this can be cast as a finite order series inT multiplying (cosh K) DN . Indeed, by invoking 1 −T 2 = 1 (cosh K) 2 and the binomial theorem, where Analogously, Equating Eqs. (2) and (7) and invoking Eq. (4) leads to a linear relation among expansion coefficients, where V and P are, respectively, DN −component and (DN + D − 1)−component vectors defined by In Eq. (10), the rectangular matrix where the DN × DN matrix M D is equal to with a square matrix A D DN 2 × DN 2 whose elements A D k,l (1 ≤ k, l ≤ DN/2) are given by Eq. (8). Constraints (4) are captured by T D in Eq. (12), where the matrix elements B D k,l = 1, if k = l, and B D k,l = 0 otherwise; O is a (D − 1) × DN 2 null matrix. Apart from the direct relations captured by Eq. (10) that relate the H-T and L-T series coefficients to each other, there are additional constraints including those (i) originating from equating coefficients of odd powers ofT and e −2K to zero and (ii) of trivial symmetry related to complimentary loops/surfaces in the H-T and L-T expansion that we discussed earlier, C i = C DN −i and C i = C DN −i . It may be verified that these restrictions are already implicit in Eq. (10). Notably, as the substitutions i ↔ (2k − i), (2l) ↔ (DN − 2l) in Eq. (8) show, Eqs. (7) and (9) are, respectively, invariant under the two independent symmetries C 2l ↔ C DN −2l and C 2l ↔ C DN −2l and thus the linear relations of Eq. (10) adhere to these symmetries. Thus, the equalities between the lowest (small 2l) and highest (i.e., (DN − 2l)) order coefficients are a consequence of the duality given by Eq. (5) that relates expansions in the weak and strong coupling parameters.
The total number of unknowns (series coefficients) in Eq. (10) is U = DN with 1/2 of these unknowns being the H-T expansion coefficients and the other 1/2 being the L-T coefficients (the components V i ). In 26 (in particular, Table 1 therein), we list the rank (R) of the matrix W D appearing in Eq. (10) for different dimensions D and number of sites N . As seen therein, for the largest systems examined the ratio R/U tends to 3/4 suggesting that in all D only ∼ 1/4 of the combined L-T and H-T coupling series coefficients need to be computed by combinatorial means. The remaining ∼ 3/4 are determined by Eq. (10). This fraction might seem trivial at first sight. If, for instance, the first 1/2 of the H-T coefficients C 2l are known (i.e., those with l ≤ DN/4) then the remaining H-T coefficients C 2l (with l > DN/4) can be determined by the symmetry relation C DN −2l = C 2l and once all of the H-T series coefficients are known (and thus the partition function fully determined), the partition function may be written in the form of Eq. (3) and the L-T coefficients {C 2l } extracted. Thus by the symmetry relations alone knowing a 1/4 of the coefficients alone suffices. The symmetry relations are a rigorous consequence of the duality relations for any value of N . As the duality relations may include additional information apart from symmetries, it is clear that R/U ≥ 3/4 for finite N (i.e., knowing more than a 1/4 of the coefficients is not necessary in order to evaluate all of the remaining H-T and L-T coefficients with the use of duality). For a given aspect ratio, the smaller N is (and the smaller the number of unknowns U ), the additional relations of Eqs. (4) carry larger relative weight and the ratio R/U may become larger. Thus, 3/4 is its lower bound. Indeed, this is what we found numerically for all (non self-dual) systems that we examined 26 . As D increases, the lowest non-vanishing orders in the L-T expansion become more separated and Eqs. (10) become more restrictive for small N systems 31 .
The H-T and L-T series are of the form of Eqs. (2) and (3) for all geometries that share the same minimal D dimensional hypercube (i.e., of minimal size L = 2) of 2 D sites. Thus, equating the series gives rise to linear relations of the same form for both a hypercube of size N = L D (with general even L) as well as a tube of N/2 D−1 hypercubes stacked along one Cartesian direction. However, although the derived linear equations are the same, the partition functions for systems of different global lattice geometries are generally dissimilar (indicating that the linear equations can never fully specify the series). Thus, the set of coefficients not fixed by the linear relations depends on the global geometry.
Parity and boundary effects may influence the rank R of the matrix W D in Eq. (10). As demonstrated in 26 for D = 2 lattices in which (at least) one of the Cartesian dimensions L is odd, as well as systems with non-periodic boundary conditions, R/U ∼ 2/3. That is, in such cases ∼ 1/3 of the coefficients need to be known before Eq. (10) can be used to compute the rest. As explained in 26 , symmetry and duality arguments can be enacted to show that in these cases, R/U ≥ 2/3 for finite N , i.e., its lower bound is 2/3. A further restriction is that of discreteness -the coefficients C 2l , C 2l (counting the number of loops/surfaces of given perimeter/surface area) must be non-negative integers for the ferromagnetic Ising model.
Let us illustrate the concepts above with a minimal periodic 2 × 2 ferromagnetic (J > 0) system with Hamil- Invoking Eqs. (11,12), There are U = 8 unknown coefficients in Eq. (16); the rank (R) of the matrix W is eight. Thus, in this minimal finite system, the Eqs. (10) are linearly independent (R/U = 1) and all coefficients may be determined ( Generally, not all coefficients may be determined by duality alone. As we discussed, in the large system limit, R/U → 3/4. A 4 × 4 example appears in 26 .

IV. PARTIAL SOLVABILITY AND BINARY SPIN GLASSES
If J ab = ±J independently on each lattice link ab , then Eqs. (4) need not hold. Instead of Eq. (10), we have 26 This less restrictive equation (by comparison to Eq. (10)), valid for all J ab = ±J, is of course still satisfied by the ferromagnetic system. For the matrix S D , a large system value of R/U ∼ 3/4 is still obtained 26 (see Table  2 therein). The partition functions for different J ab = ±J realizations will be obviously different. Nevertheless, all of these systems will share these linear relations 32 . Unlike the ferromagnetic system, the integers C l , C s may be negative. Computing the partition function of general binary spin glass D = 2 Ising models is a problem of polynomial complexity in the system size. When D ≥ 3, the complexity becomes that of an NP complete problem 33,34 . Therefore, our equations partially solve and "localize" NP-hardness to only a fraction of these coefficients; the remaining coefficients are determined by linear equations. The complexity of computing n m , required for each element of S D , is O(n 2 ). Our equations enable a polynomial (in N ) consistency checks of partition functions. In performing the expansions of Eqs. (2) and (3), the complexity of determining the number of loops (or surfaces) of given size l (or s ) (i.e., the coefficients C l or C s ) increases rapidly with l (or s ).
Our relations may be applied to systematically simplify the calculation of these coefficients. As we now explain, the situation becomes exceedingly transparent in the Ising models discussed thus far. For these theories, the coefficients are symmetric:

V. GEOMETRICAL CONSEQUENCES OF DUALITIES
Wegner's duality 35 relates interactions between {s a } on the boundaries of "d dimensional cells" to generalized Ising gauge type models with interactions between {s a } on the boundaries of "(D − d) dimensional cells". Here, a "d = 1 dimensional cell" corresponds to a (one-dimensional) nearest neighbor edge (i.e., one whose boundary is ab ) associated with standard s a s b interactions that we largely focused on thus far, d = 2 corresponds to a product of four s a 's at the centers of the four edges which form the boundary of a two-dimensional plaquette (as in standard hypercubic lattice gauge theories), d = 3 corresponds to the product of six s a 's at the center of the six two-dimensional faces which form the boundary of a three-dimensional cube, etc. The Hamiltonian is the sum of products of (2d) s a 's on the boundaries of all of the d dimensional hypercubes in the lattice (in a lattice ofÑ sites there are N c =Ñ D d such hypercubes and N s =Ñ D d−1 Ising variables s a at the centers of their faces). If the dimensionless interaction strength for a d dimensional cell is K d then the couplings in the two dual models will be related by Eqs. (5) or, equivalently, sinh 2K d sinh 2K D−d = 1. The D = 3, d = 1 duality corresponds to the duality between the D = 3 Ising model and the D = 3 Ising gauge theory. The D = 2, d = 1 case is that of the KW self-duality. For general d, Wenger derived his duality from an equivalence between the H-T and L-T coefficients.
We now turn to new, and rather universal, geometrical results obtained by our approach that hold in general dimensions. If the ground state degeneracy is 2 Ng (e.g., N g = 1 for the standard (d = 1) Ising models, N g =Ñ +2 in D = d = 2 Ising gauge theories with periodic boundary conditions), then we find 36 that the H-T and L-T series for these models are given by Eqs. (2) and (3) with the following substitution N = Nc 2l . Thus, Eq. (17) obtained for standard (d = 1) Ising models also holds for general d following this substitution. In systems with d dimensional cells, C denote, respectively, the number of closed surfaces of total d and (D − d) dimensional surface areas equal to 2l. By building on our earlier results, we observe that, when the hypercubic lattice length L is even, Eq. (17) universally relates, in any dimension D (and for any d), these numbers to each other leaving only ∼ 1/4 of these undetermined. By comparison to Eq. (17), additional geometrical conditions that hold for d = 1 (Eqs. (4)) produce the slightly more restrictive Eq. (10). Similar additional constraints appear for d > 1. A KW type self-duality present for D = 2d leads to linear equations 2l } (the number of surfaces of total D/2dimensional surface area (2l)) to themselves. We next explicitly discuss the D = 2, d = 1 case (i.e, the standard D = 2 Ising model). Similar considerations hold for any D = 2d system.

VI. DUALITIES VERSUS SELF-DUALITIES
More information can be gleaned for self-dual systems, e.g., the KW self-duality of the D = 2 Ising model. In this model, C 2l ∼ C 2l (as C 2l and C 2l are both the number of closed d = 1 dimensional loops of length 2l) when Eqs. (10) are applied to large systems (L l), see 26 . Consequently, the number of coefficients that need to be explicitly evaluated is nearly 1/2 of those obtained by matching the H-T and L-T expansions without invoking self-duality 26 : R/U ∼ 7/8 of the coefficients are determined by self-duality once ∼ 1/8 of the coefficients are provided. We caution that the relation C 2l ∼ C 2l is only asymptotically correct in the limit of large system sizes. Consequently, we find 26 that R/U asymptotically approaches 7/8 from below (and not from above as it would have if this relation were exact for finite size systems 37 ) as N becomes larger.

VII. CONSTRAINTS ON DUALITY TRANSFORMATIONS
We briefly comment on general transformations {F q } (i.e., not only those of standard weak to strong coupling dualities) and couch our guiding principle more generally. For Ising, Ising gauge, and several other theories, the pertinent mapping between the strong and weak coupling expansion parameters f +/− (K), is the Möbius transformation of Eq. (5), F (z) = (1 − z)/(1 + z), with real z (0 ≤ z ≤ 1). The duality (Möbius) transformation trivially satisfies Babbage's equation F (F (z)) = z for all values of z. For self-dual models, such as the D = 2 Ising model or D = 4 Ising gauge theories, we can easily find the critical (self-dual) point, z * , by solving the equation In general, for arbitrary models, one may find such transformations, represented by a function F (z), in terms of some parameter z (a coupling constant which can be complex-valued). Richer transformations appear in diverse arenas including Renormalization Group (RG) calculations. Based on these considerations we may have Self-dual fixed point F (F (z)) = z, Self-duality/duality F (· · · F (F (z * )) · · · ) = z * , RG fixed points.
More general transformations F 1 (F 2 (· · · F n (z) · · · )) may yield linear equations in a manner identical to those appearing for the Ising theories studied in the current work. Expansion parameters z in self-dual theories satisfy F (F (z)) = z; this yields a constraint on all possible self-dualities. Self-dual solutions are afforded by fractional linear maps F (z) = (b + az)/(−a + cz) with a 2 + bc = 0. As we will expand on elsewhere, aside from the Möbius transformation, an equivalent duality appearing in Ising and all Potts models has F 1 (z) = (az + b)/(cz + d) and F 2 (z) = (−dz + b)/(cz − a) such that F 1 (F 2 (z)) = z.
In 26 , we prove results concerning the uniqueness of the solutions. Specifically, all self-duality mappings that can be made meromorphic by a change of variables, can only be of the fractional linear type. This illustrates that all such self-dual points solve a quadratic equation. This uniqueness may rationalize the appearance of fractional linear (dual) maps in disparate arenas ranging from statistical mechanics models, such as the ones that we study here, to S-dualities in electrodynamics and YM theories.

VIII. SUMMARY
The mere existence of two or more finite order series expansions of a theory, related by dualities, may "partially solve" that theory. By "partial solvability" we allude to the ability to compute, with complexity polynomial in the system size, the full partition function given partial information (e.g., a finite fraction (1 − R/U ) of all series coefficients in the examples discussed in this work). Stated equivalently, we saw how to systematically exhaust all of the information that duality relations between disparate systems provides. This yields restrictive linear equations on the combined set of series coefficients of the dual systems. These equations allow for more than the computation of one set of (e.g., L-T or strong coupling) coefficients in terms of the other half (e.g., H-T or weak coupling). In Ising models and generalized Ising gauge (i.e., Wegner type) theories on even length hypercubic lattices in general dimensions D, only ∼ 1/4 of the coefficients were needed as an input to fully determine the partition functions; in the self-dual planar Ising model only ∼ 1/8 of the coefficients were needed as an input -the self-duality determined all of the remaining coefficients by linear relations. For an Ising chain, the H-T series expansion contains only one (two) term(s) for open (periodic) boundary conditions, i.e., Z = 2(2 cosh K) L−1 (Z = [(2 cosh K) L + (2 sinh K) L ]), thus trivially all coefficients are determined. As Ising models on varied D > 1 lattices and random Ising spin glass systems all solve a common set of linear equations, our analysis demonstrates that properties such as critical exponents cannot be determined by dualities alone.
For the even size hypercubic lattices with periodic boundary conditions studied in this work there are no closed loops (surfaces) of an odd length. Consequently, C l = C s = 0 for odd l or odd s as we have invoked. If we were to formally allow for additional odd l or s coefficients then the ratio R/U = 1/2 instead of the values of R/U that we derived (see Table I). However, when the conditions C l = 0 for odd l are imposed for the H-T coefficients these lead (via duality) to non-trivial constraints on the L-T series coefficients C l = f l ({C s }) = 0 with f l linear functions. (Similarly, a vanishing of the L-T series coefficients leads to non-trivial relations amongst the H-T coefficients.) These constraints lead to R/U > 1/2 and to universal geometric equalities discussed earlier. Not all of these equations are linearly independent; only a subset of them are. We earlier obtained lower bounds on R/U using a complementarity symmetry; the linear constraints may relate to the complementarity of the coefficients.   Table I summarizes our main findings for numerous models on even size lattices in D dimensions when these lattices are endowed with periodic boundary conditions. As noted earlier, in the supplemental material 26 we discuss other lattice sizes and boundary conditions. With the aid of our linear equations, the NP hardness of models such as the Ising spin glass in finite dimensions D > 2 is confined to a fraction (1 − R/U ) of determining all O(N ) coefficients in these models. As underscored in this work, once these are computed, the remaining fraction R/U of the coefficients are given by rather trivial linear equations. A similar matching of series, as performed in this work for the partition function, may be done for any physical quantity, such as matrix elements of operators, admitting a finite series expansion. We demonstrated that in all self-dualities (and multiple dualities), the form of the mapping can only be of a fractional linear type in appropriate coupling variables 26 . This Möbius map provides another interpretation and strategy to explore the linear constraints relating finite order series expansions of two dual models. These conclusions apply to both finite and infinite size systems. We speculate that in models with a high number of isometries (similar to those in N = 4 supersymmetric YM theories 38 ), a dominant part of the theory might become encoded in relations analogous to the linear equations studied here. Employing Cramer's rule and noting that the determinants of the matrices appearing therein are volumes of polytopes spanned by vectors comprising the columns of these matrices, relates series amplitudes to polyhedral volumes 26 . In N = 4 supersymmetric YM theories, polyhedral volume correspondences for scattering amplitudes recently led to much activity 39 .

IX. RANK OF FERROMAGNETIC ISING MODELS IN GENERAL D DIMENSIONS
Below, in Table II, we display the rank of the matrix W D of Eq. (12) in the main text. As is seen, for the larger systems examined, in all spatial dimensions D, the ratio (R/U ) tend to ∼ 3/4. For any finite N , the rank while the number of unknowns, The additive contribution of D − 1 to the rank originates from the number of conditions of Eqs. (4) of the main text. We note that although we tabulate here results for symmetric hypercubes, similar values appear in L 1 × L 2 × · · · × L D lattices for which all sides L i are even (and consequently all possible loop lengths l and surface areas s are even and lead to Eqs. (2) and (3) in the main text). As discussed in the main text, for small system size N , the conditions of Eqs. (4) therein are more restrictive. This further leads to the monotonic decrease of R/U as N is increased. For the lowest terms, when D is large, these additional conditions play a more prominent role and the problem becomes more solvable.  (4) in the main text apply), R/U tends to ∼ 3/4 for large systems.

X. RANK OF BINARY SPIN-GLASS ISING MODELS IN GENERAL D DIMENSIONS
In Table II, we examined the uniform (ferromagnetic) Ising model. We now briefly turn to Ising models with general (random) signs, J ab = ±J. In such a case, the restrictions of Eq. (4) of the main text no longer apply. Furthermore, while C 0 = 1, in the L-T expansion C 0 need not be unity (as, apart from global Z 2 (i.e., s a → (−s a ) at all lattice sites a), the ground state may have additional degeneracies). Thus, the equation to be solved is Here, V is a (DN + 1)−component vector, and Q is a (DN + 2)−component vector defined by In Eq. (21), the matrix where the (DN + 1) × (DN + 1) matrix H D is equal to , with matricesÃ D ,B D , whose elements are given bỹ The vector G is given by The last row of the matrix G and the vector Q capture the constraint As is seen in Table III, the effect of the restrictions of Eqs. (4) of the main text diminishes for large system sizes; the ratio (R/U ) ∼ 3/4 for N 1 also in this case when Eqs. (4) of the main text can no longer be invoked. For any finite N , the rank (in this case Eq. (4) of the main text no longer holds) while (as stated above), the number of unknowns or components of V is We remark that although we tabulate above the results for symmetric hypercubes, i.e., those in which along all Cartesian directions i, L i = L, similar values appear for general lattices of size L 1 × L 2 × · · · × L D with even L i (for all 1 ≤ i ≤ D). As expected, for finite N , the ratio (R/U ) in the more restricted ferromagnetic case (Table II) is always larger than that in the corresponding random ±J system (Table  III).

XI. RANK OF SELF-DUAL RELATIONS FOR THE FERROMAGNETIC D = 2 ISING MODEL
We now turn to finite size (vector) self-dualities of the D = 2 Ising model 1 with both periodic (p)/anti-periodic (a) boundary conditions, (28) whereK is the coupling dual to K (as determined by f + (K) = f − (K), i.e., sinh 2K sinh 2K = 1). Inserting Eqs. (2) and (3) of the main text into Eq. (28), allowing C γ;δ 0 , C γ;δ 0 to be different from unity ( γ = p, a; δ = p, a ), and repeating our earlier steps we obtain (with where 1 ≤ k, l ≤ N + 1. Furthermore, Setting Z p;a = Z a;p in Eq. (28), we have M self V = 0, with V a (6N + 6) dimensional vector whose components are and where M self is Here, 2 N ±1 are multiples of the identity (1) matrix, O is the null matrix, the elements of A are given by A k,l = A D=2 k−1,l−1 (1 ≤ k, l ≤ N + 1) and 1, O and A are all (N + 1) × (N + 1) matrices.  The rank R of the matrix M self is provided in Table IV. In the D = 2 self-dual Ising model the ratio R/U ∼ 7/8. In the large system (N ) limit, C 2l ∼ C 2l . Thus, by comparison to the non self-dual Ising models, there is indeed a reduction by a factor of two in the number of coefficients needed before the rest are trivially determined by the linear relations M self V = 0. It is noteworthy that, contrary to the behavior for non-self dual systems, the values of (R/U ) are monotonically increasing in N , approaching their asymptotic value of 7/8 from below.

XII. DIFFERENT SYSTEM SIZES, ASPECT RATIOS, AND BOUNDARY CONDITIONS
We now examine the analog of Eq. (12) of the main text for planar systems of size L 1 ×L 2 in which the parity of L 1,2 can be either even or odd. Systems with only even L i were examined in the main text (and Table II) to find a ratio of R/U ∼ 3/4. As seen below, when either (or both) L 1 , L 2 are odd, R/U ∼ 2/3 (i.e., 1/3 of the coefficients are needed as an input to determine the full series by linear relations).  In a system with periodic boundary conditions, when at least one of the Cartesian dimensions L i of the lattice is odd, closed loops of odd length may appear: the H-T series is no longer constrained to be of the form of Eq. (2) of the main text with even l . By contrast, as the total number of lattice links is even, the difference between low energy or "good"' (i.e., s a = s b ) and high energy "bad" (s a = −s b ) bonds s a s b is even and Eq. (3) of the main text still holds.

L1 L2 N U R
In Tables V and VI we provide the results for D = 2 periodic lattices. As seen, when at least one of the lattice dimensions is odd, R/U ∼ 2/3. The reduction in the value of R/U relative to the even size lattice with periodic boundary conditions originates from the fact that the H − T expansion admits both even and odd powers ofT ; there are more undetermined H-T coefficients. More potently, a lower bound of R/U ≥ 2/3 can be proven by invoking the symmetries captured by the duality relations. The H-T expansion is symmetric (C DN −l = C l ). In the H-T expansion both odd and even powers l appear. By contrast, in the L-T expansion only even powers s arise. The number of non-vanishing H-T series coefficients {C l } is DN ; the number of L-T expansion coefficients C s in the low temperature form of the partition function, Z L−T is DN 2 . Given {C s }, by invoking duality (i.e.,by equating Z H−T = Z L−T ) we may compute all {C l }. Thus if we compute these DN 2 L-T coefficients (a 1/3 of the combined H-T and L-T coefficients), the partition function and all remaining H-T coefficients can be fully determined.
We next turn to systems with open boundary conditions. In these systems, no odd power shows up in the H-T series expansion. By contrast, odd powers may appear in the L-T expansion (as, e.g., when there is a single s a = −1 at the boundary on a D = 2 lattice in which all other sites b have s b = +1 for which there are three "bad" bonds). In Table VII we list the matrix rank for open boundary conditions. Once again, we find that for these R/U ∼ 2/3. The reduction in the value of R/U by comparison to the even size lattice with periodic boundary conditions in this case has its origins in the fact that here the L − T expansion admits both even and odd powers of (e −2K ); there are, once again, more undetermined coefficients. We remark that the symmetries l ↔ DN − l and s ↔ DN − s of the H-T and L-T expansions that were present for periodic boundary conditions no longer appear when the system has open boundaries. In the H-T expansion, only even powers appear. In the L-T expansion both odd and even orders s are present. Thus, the number of H-T coefficients is essentially 1/2 of that of the L-T coefficients. Thus, if all of the H-T coefficients C l were known (i.e., a 1/3 of all of the combined coefficients), it is clear (even without doing any calculations) that the remaining L-T coefficients C s can be determined by duality by setting Z H−T = Z L−T and computing {C s }. Thus, for any finite size system, the fraction of coefficients required to compute the rest via the duality relations of Eqs. (10) is bounded from above, i.e., (1 − R/U ) ≤ 1/3. Asymptotically, for large N , the fraction R/U approaches 2/3 from below (that the approach is from below and not from above follows from this inequality).
Defining the parity P ≡ 0 if both L 1 and L 2 are even 1 otherwise , we find the ranks listed in Table VII.

XIII. EXPLICIT TEST CASES
In what follows, we examine specific small lattice systems via our general method.
A. Random ±J Ising systems on a periodic 2 × 2 plaquette In the main text we examined a periodic 2 × 2 ferromagnetic (J > 0) system with Hamiltonian H = −2J[s 1 s 2 +s 1 s 3 +s 2 s 4 +s 3 s 4 ]. If instead of the ferromagnetic Hamiltonian, we have an Ising model with general (random) J ab = ±J on each of the links ab then Eq. (4) of the main text will no longer hold. Instead of Eq. (10) of the main text, the equation satisfied is given by Eq. (21) where the DN dimensional (or, in this case, eightdimensional) matrix M differs from W by the omission of the matrix T (the single last column of W in this 2 × 2 example). That is, rather explicitly, S is given by As noted in the main text, the DN dimensional vector Q in the ±J Ising system differs from P in the ferromagnetic case in that it does not have an additional D − 1 (which equals one in this D = 2 example) last entries being equal to zero. It is clear that the ferromagnetic system investigated above fulfills Eq. (21) (the ferromagnetic system satisfies one further constraint related to the last row of W ). We wish to underscore that any 2 × 2 Ising system with general couplings J ab = ±J on each of the links ab will comply with Eq. (21).

B. A periodic 4 × 4 ferromagnetic system
For a periodic 4 × 4 ferromagnetic system, there are 32 coefficients and the rank of the matrix W is 26. Thus, the coefficients are linear functions of a subset of six coefficients. Choosing these to be, e.g., C 2 , C 4 , C 6 , C 4 , C 6 , C 8 the remaining coefficients are given by and C 20 = C 12 , In both of these D = 2 dimensional examples (and far more generally) the trivial reciprocity relations must be (and indeed are) obeyed.

XIV. CRAMER'S RULE AND AMPLITUDES AS POLYTOPE VOLUME RATIOS -A 2 × 2 TEST CASE ILLUSTRATION
In this section, we explicitly illustrate how the computation of the remaining series coefficients from the smaller number of requisite ones using our linear equations is, trivially, related to a volume ratio of polytopes (high dimensional polyhedra). In the main text, we remarked on the two key ingredients of this correspondence: (1) given known coefficients, we may solve for the remaining ones via our linear equations by applying Cramer's rule wherein the coefficients are equal to the ratio of two determinants. (2) the determinants (appearing in Cramer's rule) as well as those of any other matrices are equal to volumes of polytopes spanned by the vectors comprising the columns or rows of these matrices.
To make this lucid, we consider the periodic 2×2 ferromagnetic system of the main text. As it was mentioned earlier, the first six rows and last two rows of matrix W are linearly independent. Thus, in this example. we may choose these rows to construct an 8 × 8 matrixW . The corresponding vectorP that needs to be solved for satisfies the equationW V +P = 0. Here,  We may invoke Cramer's rule and find all of the undetermined coefficients, As the denominator in Eq. (38) is common to all V i , we see that V i is essentially given by the determinant det B i .
The matrix B i is obtained by replacing the i-th column ofW with (−P ). We summarize the results in Table  VIII. Putting all of the pieces together, we obtain (as we must) the exact partition function, It is hardly surprising that Cramer's rule can be applied -that is obvious given the linear equations. What we wish to highlight is that each of the determinants appearing in Cramer's rule (Eq. (38)) can be interpreted as the volume of a high-dimensional parallelepiped spanned by the vectors comprising the matrix columns. In the case above, the dimension d of each of the matrices B i andW is equal to their rank R = 8. The volume of the corresponding high dimensional tetrahedron (or polytope) spanned by the vectors forming B i is given by (det B i )/d !. In systems in which Cramer's rule may be applied, the dimensionality d is given by the rank of R of the system of linear equations, d = R.
In Figure 1, we schematically depict such a high dimensional volume.

XV. A PROOF THAT ALL MEROMORPHIC DUALITY TRANSFORMATIONS MUST BE OF THE FRACTIONAL LINEAR FORM
Charles Babbage, "the father of the computer", 2 and others since, e.g, 3,4 , have shown that the functional equation F (F (x)) = 0 has an infinite number of solutions. This observation can be summed as follows. Given any solution f , a very general class of solutions can be written as where f is a particular solution of this equation and φ is an arbitrary (or in a physics type parlance,"gauge like") function with a well defined inverse φ −1 . If we have the particular solution we can find other solutions using above equation just using a function φ with its inverse defined for a specific region. To make Babbage's observation clear, we note that if, as an example, we examine the Möbius transformation ( Figure 2) of Eq. (5), f (x) = (1 − x)/(1 + x) and consider φ(x) = x 2 and a particular branch φ −1 (x) = √ x for complex x (or the standard √ x function for real x ≥ 0) then it is clearly seen F = (1 − x 2 )/(1 + x 2 ) is also a solution to the equation F (F (x)) = x. Similarly, if we choose φ(x) = e −2x then φ −1 (f (φ(x))) = − 1 2 ln((1 − e −2x )/(1 + e −2x )) which the astute reader will recognize as the transformation of Eq. (5). That Eq. (40) constitutes a solution is seen by explicit longhand substitution. That is, where we have used the fact that f (f (x)) = x. with |z| ≤ 1, as a conformal mapping in the complex plane that maps circles onto new shifted circles with a different radius. Let us consider a circle of radius r with its center at origin. Using the transformation above, it would be mapped to a new circle of radius 2r/(1−r 2 ) with its center shifted to the point (1 + r 2 )/(1 − r 2 ) (on the real axis). Three of such circles with different colors are shown in the figure above on the lefthand side. On the righthand side we see these three circles (with the same color as on the lefthand side) after transformation. The green dot represents the self-dual point (z * = √ 2 − 1).
We now turn to a rather trivial yet as far as we are aware new result concerning this old equation that we establish here. We assert that if there exists a transformation φ that maps a complex numbers z on the Riemann sphere, z → φ(z) such that the resulting function F is meromorphic then any such function F must be of a fractional linear form, with a 2 + bc = 0. Of course, a broad class of functions of the form of Eq. (40) may be generated by choosing arbitrary φ that have an inverse yet all possible rational functions will be of the fractional linear form. For instance, the function F = (1 − x 2 )/(1 + x 2 ) discussed in an example above is, obviously, not of a fractional linear form.
Proof: The proof below is done by contradiction. A general meromorphic function on the Riemann sphere is a rational function, i.e., with P (z) and Q(z) relatively prime polynomials. (If the polynomials P and Q are not relatively prime then we can obviously divide both by any common factors that they share to make them relatively prime in the ratio appearing in Eq. (43)). As a first step, we may find the solution(s) w to the equation Unless both P and Q are linear in z, there generally will be (by the fundamental theorem of algebra) more than one solution to this equation (or, alternatively, a single solution may be multiply degenerate). That is, unless P and Q are both linear, the polynomial will be of order higher than one (m > 1) in w and will, for general z, have more than one different (non-degenerate) zero. When varying z over all possible complex values, it is impossible that the polynomial W (w) will always have only degenerate zero(s) for the relatively prime P and Q (we prove this in the rather simple Multiplicity Lemma below). We denote the general zeros of the polynomial W by w 1 , w 2 , ...., w m . That is, Now if F (F (z)) = z, then all solutions {z ji } to the equations F (z ji ) = w i (for which the polynomial (in z), W i (z) ≡ F (z) − w i Q(z) vanishes) will, for all i, solve the equation In the last equation above, on the righthand side there is a single (arbitrary) complex number z whereas on the lefthand side there are multiple (see, again, the multiplicity Lemma) viable different solutions z ji . Thus, at least one of the solutions in this set z ji = z. We denote one such solution by Z. Putting all of the pieces together, the equation F (F (t)) = t cannot be satisfied for all complex t (in particular, it is not satisfied for t = Z). Thus, both P and Q must be linear and the fractional linear form of Eq. (42) follows once it is restricted to this class. A simple multiplicity Lemma Below, we prove (by contradiction) that it is impossible that W (w) (Eq. (45)) will have an m − th order (m > 1) degenerate root for all z. We assume, on the contrary, that W z (w) = A(z)(w − B(z)) m (48) with A(z) and B(z) functions of z. At z + δz (with infinitesimal δz), the root is given by That is, by definition, 0 = W z+δz (B(z) + δB).
Therefore, w = B(z) is a root of Q(w). As the root of Q is independent of z, this implies that the assumed multiply degenerate root (i.e., B(z)) of W is independent of z. Recall (Eq. (45)) that W z = P (w) − zQ(w). As w = B(z) is (for all z) a root of both W z and Q, it follows that w = B(z) is also a root of P (w). It follows that both P (w) and Q(w) share a root (and a factor of (w − B(Z)) when factorized to their zeros), e.g., when written as P = C a (w − p a ), with C and D constants and with {p a } and {q b } the roots of P and Q respectively, at least one of the zeros ({p a }) of P must be equal to one of the zeros ({q b }) of Q. Thus, P and Q are not relatively prime if m > 1. This, however, is a contradiction and therefore establishes our assertion and proves this Lemma.