This work is concerned with zeta functions of two-dimensional shifts of finite type. A two-dimensional zeta function \(\zeta^{0}(s)\), which generalizes the Artin-Mazur zeta function, was given by Lind for \(\mathbb{Z}^{2}\)-action \(\phi\). In this paper, the \(n\)th-order zeta function \(\zeta_{n}\) of \(\phi\) on \(\mathbb{Z}_{n\times \infty}\), \(n\geq 1\), is studied first. The trace operator \(\mathbf{T}_{n}\), which is the transition matrix for \(x\)-periodic patterns with period \(n\) and height \(2\), is rotationally symmetric. The rotational symmetry of \(\mathbf{T}_{n}\) induces the reduced trace operator \(\tau_{n}\) and \(\zeta_{n}=\left(\det\left(I-s^{n}\tau_{n}\right)\right)^{-1}\).

The zeta function \(\zeta=\prod_{n=1}^{\infty} \left(\det\left(I-s^{n}\tau_{n}\right)\right)^{-1}\) in the \(x\)-direction is now a reciprocal of an infinite product of polynomials. The zeta function can be presented in the \(y\)-direction and in the coordinates of any unimodular transformation in \(GL_{2}(\mathbb{Z})\). Therefore, there exists a family of zeta functions that are meromorphic extensions of the same analytic function \(\zeta^{0}(s)\). The natural boundary of zeta functions is studied. The Taylor series for these zeta functions at the origin are equal with integer coefficients, yielding a family of identities, which are of interest in number theory. The method applies to thermodynamic zeta functions for the Ising model with finite range interactions.