Anisotropic scaling of the two-dimensional Ising model I: The torus

We present detailed calculations for the partition function and the free energy of the finite two-dimensional Ising model with periodic and antiperiodic boundary conditions, variable aspect ratio, and anisotropic couplings on the square lattice, as well as for the corresponding free energy scaling functions. Therefore we discuss the dimer mapping, the interplay between its topology and the different types of boundary conditions. As a central result we show how both the finite system as well as the scaling form decay into contributions for the bulk, a characteristic finite-size part, and - if present - the surface tension and its emergence due to at least one antiperiodic boundary in the system. For the scaling limit we expand the proper finite-size scaling theory to the anisotropic case and show how this anisotropy can be absorbed into suitable scaling variables.


Introduction
The two-dimensional Ising model on the square lattice is by far the most examined and best understood non-trivial system in statistical physics. Albeit there are still open questions, there is a whole plethora of properties known exactly, starting with the exact partition function on the torus in the thermodynamic limit calculated by Onsager [1,2], over the universal finite-size scaling at its continuous phase transition from a ferromagnetic low-temperature to a paramagnetic high-temperature phase [3,4], to the exact solutions at criticality due to conformal field theory for arbitrary geometries and boundary conditions (BCs) [5][6][7][8]. There are as many ways to calculate all these properties as there are people working on the topic, but a few are truly worth mentioning here. Beneath Onsager's expansion of the transfer matrix calculation to two dimensions, there are at least two major ideas, which are repeatedly used over last decades; one is the mapping onto spinors, as done lately by Baxter for the case of a rectangular geometry with open boundaries [9] and originally introduced by Kaufman [2]. The other one is the dimer mapping, introduced by Kasteleyn [10][11][12], refined by Fisher [13], and rigorously and exhaustively examined by McCoy & Wu [14][15][16], see [17] for a comprehensive presentation of the topic. Only quite recently, a connection between these two methods was established [18,19], which reduces the Pfaffians emerging in the dimer method to matrices corresponding to the spinor picture even for arbitrary couplings, therefore preserving the possibility to apply arbitrary boundary conditions in both directions [20]. Note that this correspondence goes beyond the simpler case of translational invariant couplings, where both methods are known to lead to the same 2 × 2 matrices. The dimer mapping will be our starting point, and we will be putting a lot of effort into its analysis relating the topology of the underlying graph and the boundary conditions of the spin system.
In this work we will focus on the finite-size contributions to the free energy in systems with periodic BCs in both directions and arbitrary couplings in one direction. This is in contrast to the works on the infinitely extended "layered Ising model" [16,21,22], were the configuration of couplings varies within a given number of layers and repeats periodically, without any finite-size contributions. In our analysis we will separate the leading contribution of the thermodynamic limit and analyse its first finite-size correction, which is responsible for the critical Casimir effect, and calculate the corresponding finite-size scaling functions. Boundaries that destroy the translational invariance in one direction as well as additional surface fields were discussed in the literature for the case of the half plane [14,15,23] and the slab geometry [24][25][26][27]; its generalisation to the cylinder with finite aspect ratios will be the topic of a subsequent paper [28].
The original inspiration for this work lies within the aforementioned universal finite-size scaling; in the vicinity of criticality the behaviour of many systems can be sorted into associated universality classes, categorised only by some rough properties like its spatial and spin dimensions, and split up into subcategories treating the BCs confining the system. Accompanying the universality, Fisher & de Gennes [29,30] predicted that the diverging correlation length ξ (∞) at criticality gives, in finite systems, rise to the thermodynamic analogue of the quantum-electrodynamic (QED) Casimir effect [31,32]. The critical fluctuations present in such a geometrically confined system lead to effective forces between the systems surfaces, namely the critical Casimir forces. In contrast to the always attractive QED Casimir effect, these thermal forces may be attractive or repulsive depending on the BCs of the surfaces, and can even change their behaviour with the temperature. This fact together with its steering temperature dependence makes them especially interesting as experimental model systems, e. g., they may be used to control the interaction strength in colloidal suspensions in order to investigate the aggregation processes [33][34][35][36]. The effect itself was first measured by Garcia & Chan as a critical thinning of a 4 He film near the λ-transition [37]. The universality of this measurement was first proven by Monte Carlo simulations of the X Y -model in three dimensional thin films with Dirichlet boundary conditions [38]. Other experiments were made on thin films of 3 He- 4 He mixtures near the tricritical point [39], again in the X Y -universality class, and binary liquids, whose demixing transition is in the Ising-universality class [40]. The latter experimental system was expanded to the direct measurement of interactions between spherical particles and a potentially chemically striped surface as well as the observation of aggregation processes in the above-mentioned colloidal suspensions [41].
For the present work, in Section 2 we first recapitulate the mapping between the twodimensional Ising model and closest-packed dimers on a carefully chosen lattice as proposed by Kasteleyn, and expand the argumentation to arbitrary nearest-neighbour couplings between each pair of sites. Therefore we will refine the argumentation by McCoy & Wu and, in contrast to the calculations for the rectangular system with open boundary conditions, we will implement translational invariance in one direction to simplify the problem further as done in [16], which leaves us with the toroidal and the cylindrical topology, where the latter one will be the topic of a subsequent paper [28]. Allowing anisotropic couplings, we will show how the introduction of dual couplings leads to a much more natural calculation compared to one of our previous works [73]. Later on, the general calculation of the emerging determinant gives us the opportunity to implement several different boundary conditions, but also makes a recapitulation of the scaling theory necessary, especially for the anisotropic case, which will be done in Section 3. Here we will also recapitulate the relation between the different points of view concerning the preferred direction within the system, e. g., the direction in which the critical Casimir force is measured, as well as the relations between the scaling functions for the free energy and the critical Casimir force.
Afterwards we will start with the calculation of the partition function for the anisotropic toroidal case in Section 4, which will be a crucial point for the calculation for the symmetric (++) and the antisymmetric (+−) boundary fields in [28], as it enlightens the interdependency of the dimer mapping and the distinction of periodic and anti-periodic boundary conditions. Additionally it gives us a good opportunity to identify the bulk contribution, which is present in all cases examined later on. We apply the anisotropic scaling theory to the Onsager dispersion relation in order to obtain the scaling forms of the relevant terms and introduce its hyperbolic parametrisation, which allows us to regularise the infinite sums by rewriting them as contour integrals over singly periodic functions in the complex plane, giving a modified Abel-Plana formula, and reproduce the free energy scaling function for the toroidal case [68,73,74]. With at least one antiperiodic boundary the system forms a domain wall, which introduces a surface tension contribution to the residual free energy, giving us a first glimpse on its general form. The cylindric geometry and boundary field effects will be discussed in [28].

The dimer representation
We will start with a brief summary of the dimer representation of the two-dimensional Ising model as it was introduced by Kasteleyn [10,11] and refined by Fisher [13] and McCoy & Wu [16], as well as an explanation of how the corresponding matrices are constructed. Then we will assume translational invariance in one direction and reduce the calculation of the determinant successively by two Schur reductions. This makes it convenient to introduce the dual couplings, which lead to a more readable and natural form. Eventually we obtain the determinant of a (quasi-)cyclic tridiagonal matrix, which can be calculated in terms of a simple 2 × 2 transfer matrix. This step will be our starting point for all further calculations, as the various BCs can be simply introduced in the reduced matrix and thus form such 2 × 2 matrices for the boundary terms. Finally we assume translational invariance in both direction, reproducing the result of McCoy & Wu.
The two-dimensional Ising model on a finite L × M square lattice with periodic boundary conditions (p) in both directions as depicted in Fig. 1, i. e., on the torus, is described by the reduced Hamiltonian (in units of k B T , with Boltzmann constant k B ) where K ⊥ ,m and K ,m are the reduced couplings between the nearest neighbours in perpendicular and parallel direction, respectively, and σ ,m ∈ {−1, +1} are spin variables with periodic indices, i. e., σ +L,m ≡ σ ,m ≡ σ ,m+M . The partition function Z (p,p) = tr e −H (p,p) can be rewritten into a high-temperature expansion [75] with z δ ,m = tanh K δ ,m (δ =⊥, ) and with the non-singular part Here we used the identity e Kσ i σ j = cosh K + σ i σ j sinh K, (2.4) as σ i σ j = ±1 depending on whether the spins are aligned or unaligned. Expanding the sum and both products leaves only even powers of the σ ,m , which are always equal to one, and a factor of 2 LM , which we already included in Z (p,p) 0 . This can be interpreted as follows: Each term in the sum gives a possible configuration of closed and commonly intersecting polygons on the original lattice, where every z ⊥/ ,m gives the position of a perpendicular/parallel link that has to be drawn, see Fig. 2 for an example.

Dimers
If we replace each site of the lattice with a properly chosen cluster, this problem is equivalent to finding the generating function of the closest-packed dimer configuration on the expanded lattice [10,11,14]. In this context a dimer is a two-atomic construct, which always occupies two sites and the connecting bond of the lattice. Fisher introduced a six-site cluster for this replacement [13], see Fig. 3(a), which can be further reduced to a four-site cluster [11], see Fig. 3(b). The dimers can be arranged in such a way that they reflect the polygon structure in a biunique way, occupying either an original lattice edge or lying within the cluster, thus leaving the edge unoccupied, see again Fig. 2.
The problem of finding the generating function of the closest-packed dimer configurations on an arbitrary planar graph was solved by Kasteleyn in terms of Pfaffians, as it gives the number of perfect matchings of a given directed planar graph with an even number of sites. This is especially powerful because of the connection between the Pfaffian and the determinant, namely Figure 3: The six-site (a) and the four-site (b) cluster that may be used to expand the lattice for the dimer mapping.
To apply this method, we have to construct a suitable skew-symmetric matrix, which corresponds to the directed graph of the lattice; the construction of this matrix will be the topic of this section, were we choose a different way than the original works of Kasteleyn and Mc-Coy & Wu (c. f. chapter IV and section 2 of chapter V in [17]) and start with matrices with only positive entries and later antisymmetrise them. The directed graph we need has to fulfil some restrictions, the most important and fundamental being that each elementary polygon is clockwise odd, that is, the number of edges pointing in clockwise direction has to be odd. This is naturally satisfied for the cluster expanded lattice, and guarantees that every term in the Pfaffian will have the same sign on any given planar graph. Consequently, for the toroidal topology of the original lattice, one Pfaffian is not sufficient, since the argument aforementioned originates in the need that every transitions cycle, which describes the transition from one dimer configuration to another, needs to be odd, see [17, Chapter 4.3-4.5] for details. On a torus this cannot be accomplished by a single but by a superposition of four Pfaffians, corresponding to the four possible combinations of periodic and anti-periodic boundary conditions in the two directions, that cancels out the redundant terms and besides gives the correct signs for the transition cycles that wind around the torus in either or both directions [76]. These four combinations of periodic and anti-periodic BCs are implemented by reversing the direction of the directed graphs in the lines responsible for the periodicity. This last point will be discussed in detail in Section 4, where the interplay of the dimer mapping and periodic and anti-periodic boundaries will be discussed in more detail, as well as later on in the upcoming second part of this paper [28], where we will map the (++) and the (+−) boundary conditions onto periodic and anti-periodic BCs on a suitably extended lattice.
To construct the matrix we will first take a look at the adjacency matrix of the graph, which decomposes into the clusters replacing the original sites and the two direction of the lattice. Using the common labels for the cluster as shown in Fig. 3(a), the graph can be represented as To simplify the problem further, the graph can be reduced into a non-planar one with only four sites, eliminating the two interior sites labeled 1 and 2, e. g., as the Schur complement (2.25) of the antisymmetric form of (2.6) with respect to the central 2 × 2 block. The related cluster is depicted in Fig. 3(b), and its adjacency matrix is Nevertheless, this reduction introduces two additional terms for each lattice site of the original system that is not part of any polygon, as there are now three combinations on the Kasteleyn cluster for those cases. But if each bond within the clusters has a weight equal to 1, those additional terms cancel each other and we ensure that no cluster has an influence on (2.2). On this resolution, the coupling in parallel and perpendicular direction are represented by the two matrices respectively. The nearest-neighbour structure in a row and a column are both represented by the n × n matrix , 0, −1} accounting for different BCs, i. e., b = 0 for open, b = +1 for periodic, and b = −1 for anti-periodic boundaries of the Ising system, concerning the necessity of transition cycles on the graph to be odd. Note, that b = −1 accounts for periodic boundaries on the directed graph, i. e., in the dimer system, in the sense that all edges are likewise aligned, while it accounts for antiperiodic BCs in the Ising model. However, the topology of the underlying directed graph is not representative for the Ising model, which is emphasised by the fact that the Ising partition function is a combination of four Pfaffians. The identification with periodicity and antiperiodicity for the values of b stems solely from the identification within the open cylinder, see [28,73], and thus we will use "+" to mark periodic and "−" to mark antiperiodic BCs. If we give each bond a different weight, see Fig. 1, we can represent the parallel and perpendicular couplings by with Z δ = diag(z δ 1,1 , z δ 1,2 , . . . , z δ L,M ) and the n × n identity matrix 1 n . The final 4L M × 4L M adjacency matrix of the graph then reads and, since two vertices are connected at most by a single edge, we can expand it to a skewsymmetric block form 1 to calculate the partition function of the torus as

Translational invariance
Since we are interested only in the boundary conditions along the parallel direction, we now assume translational invariance in this direction and thus present an alternative solution to the problem of the layered Ising model presented in [16]. In contrast to the original calculations we do not restrict ourselves to the cylindrical case but include it as the special case of the torus, with β = 0. The Ising model on the rectangle with open boundaries in both directions is discussed in [9,18,19] for all temperature and in [77,78] for the critical case. To account for the translational invariance in parallel direction, all M couplings in each column are the same, z δ ,m ≡ z δ ∀ m, and the matrices Z δ simplify to therefore it commutes with its transposed, and thus the unitary matrix U ± that diagonalises H ± also does so for H ± . For convenience we decompose the matrix A αβ into its three contributions (2.10) as where In the expanded space, U β=± becomes and has no influence on the first two contributions as they commute, For the parallel contribution, we obtain the diagonal matrices for m ∈ {0, 1, 2, . . . , M −1}, and thus we will call β = −1 even and β = +1 odd. Note that the eigenvalues lie equidistantly on the unit circle and thus we have a free shifting parameter for the spectrum. We have chosen it in such a way, that the eigenvalue ϕ (−) 0 = 0 appears in the even spectrum; a shift by −π on the other hand would have given rise to a dependency on whether M is even or odd. For completeness we present the characteristic polynomials here, as we will need them later, too. Note that we have chosen to use the normalised version of P ± β to omit the aforementioned dependency on whether M is even or odd.
Since det(A⊗ B) = det(B ⊗ A), we can rearrange the matrices such that the product is block diagonal, simplifying its determinant to the form where we defined J α = 1 + z ⊥ H α for the perpendicular direction and the diagonal matrices At this point we introduce the dual couplings t ≡ (z ) * , which can be defined within the low-temperature expansion of the partition function analogous to (2.2), but on the dual lattice. They are connected to the couplings z via the self-inverse duality transform x * of an arbitrary quantity x, Thus we will omit the superscripts in the following and write z ≡ z ⊥ and rewrite the parallel couplings z through their duals t * . In the following we will, if necessary, mark a dual coupling by an asterisk as shown in (2.24). Consequently we may also rewrite the two matrices z ⊥ → z and z → t * . As shown in [73]  where we will make the reduction with respect to upper left 2 × 2 block matrix, as for the first step this part solely represents the couplings in the translationally invariant, parallel direction. The first reduction gives For the second determinant in (2.26) we apply a second Schur reduction with respect to the upper left block matrix, which yields Again, the first determinant is rather simple, and can be combined with (2.27) to We finally find The L × L matrixC α is a symmetric tridiagonal (quasi-)cyclic matrix of the form Here we have used the convenient abbreviation [18] x for an arbitrary quantity x, as well as which solely covers the ϕ (β) m -dependency. The matrixC α in (2.33) will be our starting point for the different BCs, but as we are hardly interested in the case of every line having a different coupling, we will apply the changes due to the corresponding boundary conditions and then frankly assume homogenous but still anisotropic couplings. The determinant of a matrix of the form (2.33) can be easily calculated with a 2 × 2 transfer matrix ansatz, namely analogous to [14] for a system with open boundaries in one direction. Note that this formula is formally equivalent to the more general form in [18], where no translational invariance in the direction perpendicular to the open boundaries is assumed. Luckily the b are all negative, giving an additional factor (−1) L to (2.36) and thus eliminating the dependency on whether L is even or odd. Thus our results are correct for arbitrary integer L and M . If we assume translational invariance in both directions, we may reduce the coupling matrices further to z ⊥ = z ⊥ 1 and z = z 1 and diagonalise the perpendicular direction, too. Then the determinant becomes

Scaling theory
The two-dimensional Ising model is one of the simplest systems showing a temperature-driven continuous phase transitions in the absence of a bulk magnetic field. Albeit this work focuses on the two-dimensional case, the following statements can be generalised to d > 2. We assume an L × M rectangular system, where L is its extent in the perpendicular (⊥) direction and M the extent in parallel ( ) direction, for reasons that will be clear later on. In the vicinity of the phase transition the thermal fluctuations of the medium become long-ranged, but by imposing some sort of BCs, even (anti-)periodicity, they are confined by the geometry. The reduced free energy 2 for arbitrary BCs α in perpendicular direction and translationally invariant BCs where the first term with bulk and surface free energy densities f b and f (α) s , describes the leading behaviour in the thermodynamic limit L, M → ∞, while the latter one covers the finite-size effects.
Another way to decompose the free energy [18,19] is into contributions F cf. [19]. Note that it is possible to impose boundaries in more than one direction, introducing edges and corners and suitable contributions depending on the dimension of the hypercuboid; for the two dimensional case of the Ising model on the open rectangle the necessary calculations were done lately with the dimer approach presented above [18,19], while the corresponding leading term F ∞ (L, M ) was analysed in detail within the spinor representation by Baxter [9].

General scaling behaviour
Near criticality the bulk correlation length in direction δ diverges as where τ = T /T c − 1 is the reduced temperature andξ δ is the correlation length amplitude in the unordered phase, while ν is the associated scaling exponent, with ν = 1 for the twodimensional case. In the region around T c where the correlation length is the dominant length scale, the residual free energy F (α,β) ∞,res only depends on the two length ratios L/ξ where we assume the bulk correlation length amplitudes differ depending on the direction because of the anisotropic couplings. Following Fisher & de Gennes [29], the residual free energy thus fulfils the scaling ansatz where Θ (α,β) is the total residual free energy scaling function, depending on the two temperature scaling variables which are related by the reduced aspect ratio via the relation Additionally, it is sometimes [68] advantageous to consider scaling functions depending on a volume-like scaling variable Consequently, we may change our focus onto either one of the two directions or the volume and rewrite the scaling functions accordingly into a perpendicular (⊥), a parallel ( ), or a volume (•) form The critical Casimir force is defined as derivative of the residual free energy with respect to the corresponding system length, e. g. in perpendicular direction, ∞,res (L, M ) (3.14) and analogously to (3.13) scales as The residual free energy and the Casimir force scaling functions are connected by the relations [68,79] In the following we will focus on the scaling function Θ (α,β) (x , ρ) for the parallel direction, as it arises naturally within our calculations from the product in (2.22). Just like the residual free energy can be decomposed into its bulk, surface, and strip contribution in (3.6), the scaling function consists of such parts [19], which can be written as For the torus and the other systems with translational invariance in both directions there is no surface and thus Θ (α,β) s (x ) ≡ 0 for α, β ∈ {p, a}, all other cases will be discussed in the second part of this paper [28].

Anisotropic scaling
Now we need to discuss the influence of the weakly anisotropic couplings on the scaling behaviour. Here the critical point expands to a critical line as we may rewrite the condition of criticality [80] as (3.19) which again becomes a point by fixing the couplings by a factor κ as In terms of t and z this becomes and we can identify the critical point with the equation The correlation lengths along the two directions read [17] ξ (∞) with the dual couplings z * and t * from (2.24). This form emphasises the choice of the reduced temperatures in the two directions as The amplitude ratio r ξ ≡ξ ⊥ /ξ ξ [81] can be expanded around criticality t = z using log(1 + x) x to give r ξ τ /τ ⊥ , such that the scaling variables fulfil

The torus
After we have calculated the partition function on the finite lattice, we will carve out the bulk properties and calculate its thermodynamic limit. Afterwards we will see how the dimer approach and the choice of signs distinguishes between periodicity and antiperiodic BCs as a consequence of its interplay with the directed graph being subject to the Pfaffian. Let us return to (2.33), for which we now assume independent homogeneity in both directions, i. e., z ≡ z ∀ and t ≡ t ∀ , as well as µ ± (ϕ) ≡ µ ± (ϕ) ∀ . Then it is comfortable to factorise the homogeneous factor −z/t − , leaving every non-zero off-diagonal element but b L equal to −1. Due to this procedure, the diagonal entries simplify to which can be parameterised into the anisotropic Onsager dispersion relation Thus the matrix for the torus has the form which is connected to (2.33) by To calculate this determinant, we use the transfer matrix approach (2.36) to find with the transfer matrix The eigenvalues of T are e ±γ (β) m and, depending on α = ±1, the determinant thus reads det C (4.7) In the following we will omit the superscript of the γ, since the parity is encoded in the product over either even or odd numbers. Having said that, we denote the parity in the other direction with ± to account for the dependency on α, and thus we use where we have incorporated the square root of (2.5). Then the non-regular part of the partition function reads where the sgn(t − z) stems from the root in (2.5), too, and assures that for arbitrary anisotropy κ the contribution is a monotonic function in temperature. This can be easily understood by looking at the critical value of Z − e : There γ 0 = 0 and thus one of the factors becomes zero, but as we have taken the square root in (4.9), we have to correct the respective sign. Of course this is just exactly the result by Kaufman [2].

Bulk and finite-size contribution
In section 3 we discussed how the free energy of a finite system may be decomposed into summands describing the different contributions from volume, surfaces, and finiteness. As a matter of fact, the toroidal geometry has no surfaces and thus the aforementioned decomposition includes only the bulk contribution and the finite-size part. Therefore it is perfectly suitable to identify the former one and calculate its thermodynamic limit.
We follow the procedure of Ferdinand & Fisher [74] and start by splitting the Z ± e/o into their exponentially growing and decaying parts as Now we may factorise the odd product over the leading exponential in (4.9) from all four contributions, as it is slightly larger than the one over the even numbers. The corresponding bulk part of the free energy then simply reads thus leaving the residual part as To calculate the bulk contribution to the free energy in the thermodynamic limit, 14) Figure 4: Bulk free energy density f b (z, t). The dashed line marks criticality for arbitrary anisotropy κ = K ⊥ /K . For fixed anisotropy the two couplings are connected by t = (z * ) κ , see (3.21), and the black lines mark the run of the corresponding curve, where the isotropic case κ = 1 is shown in red.
we use the Euler-Maclaurin sum formula to obtain the integral representation and omit all corrections of O(M −1 ) or higher. This is in perfect agreement with former results by McCoy & Wu [14]. Fig. 4 shows the bulk free energy density as function of the two coupling variables z and t. For the isotropic critical case we find where G is Catalan's constant. Together with the contribution from the regular part of the partition function, this coincides perfectly with the result by Izmailian [60].

Scaling form
As discussed in Sec. 3, in a finite system there is always not only a residual contribution present, but this contribution can be split up into different parts as well. We start with the scaling form of γ as it is the central quantity for the toroidal geometry, afterwards we will first calculate the residual bulk free energy scaling function and then the one for the residual strip free energy. Therefore we introduce the hyperbolic parametrisation based on the scaling form of γ, which allows us to regularise the arising terms properly. At criticality the bulk residual scaling function coincides with the change of the free energy of the associated conformal field theory due to the projection from the plane onto the cylinder.
For the scaling limit of γ we start with the definition for the anisotropic temperature scaling variables in (3.24) and (3.25) and solve them for z and t to get the somehow unhandy formulas (4.17b) Now we use them in (4.2) together with ϕ = Φ/M , substitute M = ε −1 and expand cos(ε Φ) up to the second order in ε around ε = 0 to obtain where we have introduced With a Puiseux series expansion around ε = 0, we finally find, with L = ρ r ξ M , As we made only the most general assumptions, we can use this result throughout the whole rest of this work. We finally note that (4.19), rewritten as is the finite-size scaling form of the Onsager dispersion relation (4.2) just as in the isotropic case [19]. Now we will turn to the characteristic polynomials (2.21). Their scaling form is quite simple as we only have to replace ϕ = Φ/M to obtain P ± e (Φ) = e ±iΦ − 1, (4.23a) where the ± accounts for the two possible choices of the eigenvalue e ±iϕ (β) m in (2.19). This freedom is essential to the regularisation of the residual free energies, because if we assume Φ to be complex to calculate the sum in terms of contour integrals, P + e/o (Φ) diverges in the lower and P − e/o (Φ) diverges in the upper half-plane. To avoid these divergences we can switch between the two possible realisations if we cross the real axis. We obtain suitable counting polynomials for the even and odd sums as integration kernel from the logarithmic derivative of the P ± e/o (Φ) as 24b) Figure 5: Complex structure of the integrand of (4.27). The two contours C + in the upper and C − in the lower half-plane are connected by non-contributing paths at ℜ{ω} → ±∞ to form the closed contour C. The complex phase is colour coded from −π to π, while the lines of constant absolute value c are shown as black dotted (c < 1), dashed (c = 1) or solid (c > 1) lines ranging from 2 −3 to 2 +3 . The zeros at ω = ±iπ/2 are marked as green dots, while the poles on the lines with ℑ{ω} = nπ, n ∈ , are marked as red dots. Additionally, the phase jump due to the different choices for the counting polynomial for the upper and lower half-plane are marked as red lines. Note that due to the transformation into the hyperbolic ω-plane the structure is 2π-periodic along the imaginary axis. The contributions along the paths C + and C − are equal, as the phase in the upper half-period is reversed with respect to the lower half-period. and consequently for the alternating sum solving the equation for arbitrary ω ∈ . Note that this is the critical limiting case of the elliptic parametrisation of (4.2), which is used off criticality for finite systems with more complex BCs, e. g., for the open rectangle [9]. Combining (4.21), (4.24b) and (4.25) we now have the tools to calculate the scaling func- Figure 6: Complex structure of the alternating sum over the Γ m in the hyperbolic ω-plane (integrand of (4.31a)), together with the contour C = C + + C − . Here no phase jump is present and the 2π-periodicity along the imaginary axis is continuous, but nevertheless the upper and lower contour are equal due to the reversed phase in the upper and lower half-period.
tion of the residual bulk free energy By rewriting the sum as contour integral over Γ with the counting polynomial K ± o and substituting the hyperbolic parametrisation we get Note that the additional term ±i in K ± o can be interpreted as the integral in (4.26) for the bulk free energy, as and thus a shift of the integration path along the imaginary axis can be interpreted as smooth interpolation between the sum on the real axis and the integral at infinity, which is mapped onto the lines at ±iπ/2 within the hyperbolic parametrisation. The integrand is shown in Fig. 5 together with an appropriately chosen contour. For |Φ| → ∞ the integrand vanishes, thus making the integral over the paths ℑ{ω} = ±π/2 the only relevant parts of the contour, which are equal as they only differ by a phase of π. This leaves us with where we first resubstituted to the Φ-plane and then integrated by parts, reproducing the result from [73] calculated for the isotropic case.
A similar procedure will be our way to go in every upcoming calculation, that is, use the associated counting polynomial as integral kernel in the contour integral representation of the sum, transform it into the hyperbolic parametrisation, shift the integration path to ℑ{ω} = ±π/2, and perform the integration. Applying this scheme to the alternating sum in the exponent of (4.13), then gives us again reproducing the result of [73]. Note that the alternating sum regularises itself and thus there is only the well known kernel for alternating sums present, which is shown in Fig. 6 in the same manner as Fig. 5.
Lastly we have to calculate the scaling limit of the p ± e/o (L, M ), which is really straight forward, as we only have to expand the product symmetrically to infinity and replace Lγ by its scaling form. Thus we get with where we have used the discrete form of Γ , Due to the periodic BCs the system has the symmetry Θ (p,p) (x , ρ) = Θ (p,p) ⊥ (x ⊥ , 1/ρ). The limit (4.39) for x δ → −∞ due to the degeneracy is shown for corresponding values of ρ as dotted lines.
with Φ m = mπ, and used the symmetries of γ m with respect to m = 0 and m = M to expand the product symmetric around zero. Now we combine them into an even and an odd scaling function and calculate the analogue of an alternating strip contribution Finally we find the scaling function for the torus to be see Fig. 7, which has some remarkable properties. First and foremost this form shows the symmetry as well as δΘ b and δΨ (p,p) exchange their roles under this transformations, which may be seen best in the two limiting cases where for the latter one P + o becomes dominant near x ⊥ = 0 and equal to P − o for |x ⊥ | 1, thus Ψ Due to the antiperiodic BCs the system has the symmetry Θ (a,a) (x , ρ) = Θ (a,a) ⊥ (x ⊥ , 1/ρ). For negative values of x δ the scaling function diverges linearly, marked as red (dark) dashed line. The scaling function converges very fast to its limiting form at ρ = ∞, while its biggest contribution is at ρ = 1 as there the domain wall lies diagonal in the system and is thus longest compared to the two length scales L and M .
which stems from the different order of the limits in the scaling regime and the thermodynamic limit. While in an infinitely large system as supposed by the latter one, the system is frozen in either of the two possible magnetised states because a transition would require infinitely much energy, the finiteness of the system in the scaling limit allows such a transition, resulting in a factor 2 in the partition function and thus this topological contribution in the scaling function. However, whenever this symmetry is broken, e. g., by a surface field, this contribution is not present and the scaling functions decay to zero for |x δ | → ∞. Indeed our results for the torus coincide with the results by Ferdinand & Fisher [74], as well as with the more recent results in [68,73], although we did not assume isotropy. Thus we showed explicitly for finitesize scaling functions of the Ising model that a coupling anisotropy can be absorbed into a generalised aspect ratio ρ as proposed, e. g., in [81], and that therefore the universal finitesize scaling functions are independent of the coupling anisotropy.

Antiperiodicity and surface tension
Let us now turn back to the finite systems solution using the dimer approach and consider antiperiodic boundary conditions in the parallel direction, i. e, σ ,m ≡ −σ ,m+M . They can be implemented as one line of anti-ferromagnetic parallel couplings, which corresponds to a change in the orientation of the oriented lattice concerning this particular line within the dimer representation. Here we need to emphasise that the oriented graph from the dimer mapping is absolutely independent of the choice of sign of any coupling. Nevertheless, changing the sign of all couplings in one row (or column) can be interpreted as reversing the direction of all corresponding edges and leaving the couplings unchanged. We now may choose the line, which imposes these BCs to be the line concerning the handling of the transition circles An important consequence of an antiperiodic boundary is the formation of a domain wall and thus a surface tension contribution σ to the free energy. This contribution can be calculated as difference between the free energy of the system with antiperiodic and the one of the  b (x δ ) is the bulk contribution to the surface tension, cf. [28], and H(x) is the Heaviside step function with H(x) ≡ d dx max{0, x}.

Conclusion
We presented a systematic calculation of the universal free energy finite-size scaling functions for anisotropic Ising systems with translationally invariant boundary conditions. Therefore we started with the commonly known dimer representation of the two-dimensional square lattice Ising model, which we recapitulated in some detail to show the connection between the Pfaffian representation and the boundary conditions of the underlying system. Using the coupling constant z and its dual counterpart t turned out to be beneficial for the calculation, as it compresses the formulas to a clearly arranged form. The calculation of the necessary Pfaffians of the original 4L M × 4L M matrices was condensed down to a straightforward formalism for periodic and antiperiodic boundaries assuming translational invariance in one direction. Within this approach we implemented an anisotropic scaling theory and showed analytically that a coupling anisotropy can be absorbed into a generalised aspect ratio ρ. For the transition from the finite system to the scaling form we introduced the hyperbolic parametrisation of the scaling form of the Onsager dispersion relation, which is the scaling limit of the elliptic parametrisation of finite systems. The results are in perfect agreement with former calculations, but the method is far more general and we will use it to apply boundary fields to the open cylinder in a subsequent paper [28].