Abstract Based on the results obtained in (Hucht 2017 J. Phys. A: Math. Theor. 50 065201), we show that the partition function of the anisotropic square lattice Ising model on the L × M rectangle, with open boundary conditions in both directions, is given by the determinant of an M/2 × M/2 Hankel matrix, that equivalently can be written as the Pfaffian of a skew-symmetric M × M Toeplitz matrix. The M − 1 independent matrix elements of both matrices are Fourier coefficients of a certain symbol function, which is given by the ratio of two characteristic polynomials. These polynomials are associated to the different directions of the system, encode the respective boundary conditions, and are directly related through the symmetry of the considered Ising model under exchange of the two directions. The results can be generalized to other boundary conditions and are well suited for the analysis of finite-size scaling functions in the critical scaling limit using Szegő’s theorem.