The Exact Susceptibility of the Spin-S Transverse Ising Chain with Next-nearest-neighbor Interactions.

Let us consider the Hamiltonian

$${\cal H}=({\cal H}_{0}+{\cal H}_{1})-H{\cal Q},$$

(1)

where

$\displaystyle{\cal H}_{0}=-J_{0}\sum_{j=1}^{N}s_{j}^{z}s_{j+1}^{z},\hskip 22.76228pt{\cal H}_{1}=-J_{1}\sum_{j=1}^{N}s_{j}^{z}s_{j+2}^{z},\hskip 22.76228pt{\cal Q}=g\mu_{\rm B}\sum_{j=1}^{N}s_{j}^{x},$

(2)

and the periodic boundary condition is assumed. Here ${\cal H}_{\rm Is}={\cal H}_{0}+{\cal H}_{1}$ is the Ising interactions, and ${\cal Q}$ is the transverse term which does not commute with ${\cal H}_{\rm Is}$.

The transverse susceptibility $\chi$ at zero-field is defined as

	$\displaystyle\chi$	$\displaystyle=$	$\displaystyle\frac{\partial}{\partial H}\langle{\cal Q}\rangle|_{H=0}$		(3)
		$\displaystyle=$	$\displaystyle\frac{\partial}{\partial H}\frac{{\rm Tr\>}{\cal Q}\exp(-\beta{\cal H})}{{\rm Tr}\exp(-\beta{\cal H})}\Big{|}_{H=0}.$

The expectation $\langle\>\rangle$ is taken by ${\cal H}$, and $\beta=1/k_{\rm B}T$ where $k_{\rm B}$ is the Boltzmann constant and $T$ is the temperature. For the purpose to calculate (3), let us introduce the formula

$$\frac{d}{dH}\exp(-\beta{\cal H})=\exp(-\beta{\cal H})\int_{0}^{\beta}\exp(\lambda{\cal H}){\cal Q}\exp(-\lambda{\cal H})d\lambda,$$

(4)

which can be derived by multiplying $\exp(\beta{\cal H})$ from the left-hand side and differentiating in terms of $\beta$. Using (4), $\langle{\cal Q}\rangle$ can be expanded in terms of $H$, and we obtain

	$\displaystyle\chi$	$\displaystyle=$	$\displaystyle\frac{\partial}{\partial H}[\frac{\langle{\cal Q}\rangle_{0}+H\int_{0}^{\beta}\langle\exp(\lambda{\cal H}_{\rm Is}){\cal Q}\exp(-\lambda{\cal H}_{\rm Is}){\cal Q}\rangle_{0}d\lambda}{1+H\int_{0}^{\beta}\langle{\cal Q}\rangle_{0}d\lambda}+O(H^{2})]\Big{|}_{H=0}$		(5)
		$\displaystyle=$	$\displaystyle\int_{0}^{\beta}\langle\exp(\lambda{\cal H}_{\rm Is}){\cal Q}\exp(-\lambda{\cal H}_{\rm Is}){\cal Q}\rangle_{0}d\lambda,$

where $\langle\>\rangle_{0}$ is the expectation taken at zero external field $H=0$.

Next let us introduce inner derivatives $\delta_{0}$ and $\delta_{1}$ through the relations

	$\displaystyle\delta_{0}^{0}{\cal Q}={\cal Q},\hskip 22.76228pt\delta_{0}^{n+1}{\cal Q}=[{\cal H}_{0},\delta_{0}^{n}{\cal Q}],$
	$\displaystyle\delta_{1}^{0}{\cal Q}={\cal Q},\hskip 22.76228pt\delta_{1}^{n+1}{\cal Q}=[{\cal H}_{1},\delta_{1}^{n}{\cal Q}],$		(6)

and introduce $\delta=\delta_{0}+\delta_{1}$. The expectation in (5) can be expanded using

$$\exp(\lambda{\cal H}_{\rm Is}){\cal Q}\exp(-\lambda{\cal H}_{\rm Is})=\sum_{p=0}^{\infty}\frac{\lambda^{p}}{p!}\delta^{p}{\cal Q},$$

(7)

and (5) is expressed as

$$\chi=\int_{0}^{\beta}\langle\sum_{p=0}^{\infty}\frac{\lambda^{p}}{p!}\delta^{p}{\cal Q}\>{\cal Q}\rangle_{0}d\lambda.$$

(8)

It can be seen that (8) is equal to

	$\displaystyle\chi$	$\displaystyle=$	$\displaystyle\sum_{p=0}^{\infty}\frac{1}{p!}\int_{0}^{\beta}\lambda^{p}d\lambda\langle\delta^{p}{\cal Q}\>{\cal Q}\rangle_{0}$		(9)
		$\displaystyle=$	$\displaystyle\sum_{p=0}^{\infty}\frac{\beta^{p+1}}{(p+1)!}\langle\delta^{p}{\cal Q}\>{\cal Q}\rangle_{0}.$

This formula was already derived in Ref.\citen96Minami. What we have to do next is to evaluate the correlation functions $\langle\delta^{p}{\cal Q}\>{\cal Q}\rangle_{0}$ with non-zero next-nearest-neighbor interaction, and to perform the infinite sum in (9).

Then let us consider $\langle\delta^{p}{\cal Q}\>{\cal Q}\rangle_{0}$. Because of the commutatibity ${\cal H}_{1}{\cal H}_{0}={\cal H}_{0}{\cal H}_{1}$, we generally find $\delta_{1}\delta_{0}=\delta_{0}\delta_{1}$, and then $\delta^{p}{\cal Q}$ is written as

$\displaystyle\delta^{p}{\cal Q}$

$\displaystyle=$

$\displaystyle(\delta_{0}+\delta_{1})^{p}{\cal Q}=g\mu_{\rm B}\sum_{l=0}^{p}\left(\begin{array}[]{c}p\\ l\end{array}\right)\delta_{0}^{p-l}\delta_{1}^{l}\sum_{j=1}^{N}s^{x}_{j},$

(12)

where

$\displaystyle\left(\begin{array}[]{c}p\\ l\end{array}\right)=\frac{p!}{(p-l)!l!}.$

(15)

It is straightforward to show inductively that

	$\displaystyle\delta_{0}^{p_{0}}\delta_{1}^{p_{1}}\sum_{j=1}^{N}s^{x}_{j}$
	$\displaystyle=\sum_{j=1}^{N}\Big{(}(-J_{0})^{p_{0}}\sum_{l_{0}=0}^{p_{0}}\left(\begin{array}[]{c}p_{0}\\ l_{0}\end{array}\right)(s^{z}_{j-1})^{p_{0}-l_{0}}(s^{z}_{j+1})^{l_{0}}\Big{)}\Big{(}(-J_{1})^{p_{1}}\sum_{l_{1}=0}^{p_{1}}\left(\begin{array}[]{c}p_{1}\\ l_{1}\end{array}\right)(s^{z}_{j-2})^{p_{1}-l_{1}}(s^{z}_{j+2})^{l_{1}}\Big{)}\tau_{p_{0}+p_{1}}s^{k(p_{0}+p_{1})}_{j},$		(20)
			(21)

where

$\displaystyle\tau_{p}=\left\{\begin{array}[]{cl}1&\hskip 8.5359ptp={\rm even}\\ i&\hskip 8.5359ptp={\rm odd}\end{array}\right.$

(24)

and $k(p)=x$ for even $p$, and $k(p)=y$ for odd $p$. Then from (12) and (21), we find

	$\displaystyle\langle\delta^{p}{\cal Q}\>{\cal Q}\rangle_{0}$	$\displaystyle=$	$\displaystyle(g\mu_{\rm B})^{2}\sum_{l=0}^{p}\left(\begin{array}[]{c}p\\ l\end{array}\right)\langle\delta_{0}^{p-l}\delta_{1}^{l}\sum_{j=1}^{N}s^{x}_{j},\sum_{j^{\prime}=1}^{N}s^{x}_{j^{\prime}}\rangle_{0}$		(27)
		$\displaystyle=$	$\displaystyle(g\mu_{\rm B})^{2}\sum_{j=1}^{N}\sum_{l=0}^{p}\left(\begin{array}[]{c}p\\ l\end{array}\right)(-J_{0})^{p-l}(-J_{1})^{l}$		(30)
		$\displaystyle\times$	$\displaystyle\sum_{l_{0}=0}^{p_{0}}\left(\begin{array}[]{c}p_{0}\\ l_{0}\end{array}\right)\sum_{l_{1}=0}^{p_{1}}\left(\begin{array}[]{c}p_{1}\\ l_{1}\end{array}\right)\langle(s^{z}_{j-2})^{p_{1}-l_{1}}(s^{z}_{j-1})^{p_{0}-l_{0}}(\tau_{p}s^{k(p)}_{j}s^{x}_{j})(s^{z}_{j+1})^{l_{0}}(s^{z}_{j+2})^{l_{1}}\rangle_{0},$		(35)
					(36)

where $p_{0}=p-l$, $p_{1}=l$ and thus $p_{0}+p_{1}=p$. Here we used the fact that eigenstates of ${\cal H}_{\rm Is}$ are direct products of eigenstates of each $s_{j}^{z}$, and hence $\langle(s^{z}_{j-2})^{p_{1}-l_{1}}(s^{z}_{j-1})^{p_{0}-l_{0}}\tau_{p}s^{k(p)}_{j}(s^{z}_{j+1})^{l_{0}}(s^{z}_{j+2})^{l_{1}}\cdot s^{x}_{j^{\prime}}\rangle_{0}=0$ when $j^{\prime}\neq j$. The correlation function in (36) can be obtained through the transfer matrix method.

The transfer matrix is defined as the matrix composed of the Boltzmann factors between two sites as its elements:

$\displaystyle V=\sum_{s_{j},s_{j+1}}|s_{j+1}\rangle e^{Ks_{j}s_{j+1}}\langle s_{j}|,$

(37)

and thus

$\displaystyle\langle s_{j+1}|V|s_{j}\rangle=e^{Ks_{j}s_{j+1}}.$

(38)

Here we consider the state $|s_{j}\rangle_{j}$ at site $j$ with the eigenvalue $s_{j}$ $(s_{j}=S,S-1,\ldots,-S)$: $s^{z}_{j}|s_{j}\rangle_{j}=s_{j}|s_{j}\rangle_{j}$ , and $|s_{j}\rangle_{j}$ is written $|s_{j}\rangle$ for abbreviation. If one operates $V$ from the left-hand side, then one additional bond is introduced to the chain; for example $V|s_{j}=s\rangle$ provides all the $(j+1)$-th states $|s_{j+1}\rangle$ multiplied by the Boltzmann factor obtained under the condition that the site $j$ has fixed eigenvalue $s_{j}=s$. In our case, we consider the configurations of adjacent two sites to introduce the next-nearest-neighbor interactions. The transfer matrix, in our notation, is defined as

$\displaystyle V=\sum_{s_{j-1},s_{j},s^{\prime}_{j},s^{\prime}_{j+1}}|s^{\prime}_{j}\rangle|s^{\prime}_{j+1}\rangle e^{K_{0}s_{j-1}s^{\prime}_{j}}e^{K_{1}s_{j-1}s^{\prime}_{j+1}}\delta_{s_{j}s^{\prime}_{j}}\langle s_{j-1}|\langle s_{j}|,$

(39)

and thus

$\displaystyle\langle s^{\prime}_{j}|\langle s^{\prime}_{j+1}|\>V\>|s_{j-1}\rangle|s_{j}\rangle=e^{K_{0}s_{j-1}s^{\prime}_{j}}e^{K_{1}s_{j-1}s^{\prime}_{j+1}}\delta_{s_{j}s^{\prime}_{j}},$

(40)

where $K_{0}=\beta J_{0}$, $K_{1}=\beta J_{1}$, and $s_{j-1},s_{j},s^{\prime}_{j},s^{\prime}_{j+1}=S,S-1,\ldots,-S$.

For example in the case of $S=1/2$, we find

$\displaystyle V=\left(\begin{array}[]{cccc}xy&0&1/xy&0\\ x/y&0&y/x&0\\ 0&y/x&0&x/y\\ 0&1/xy&0&xy\end{array}\right)$

(45)

where $x=e^{\frac{1}{4}K_{0}}$ and $y=e^{\frac{1}{4}K_{1}}$. The configurations $\{|s_{j-1}\rangle|s_{j}\rangle\}$ and $\{|s^{\prime}_{j}\rangle|s^{\prime}_{j+1}\rangle\}$ are both arranged by binary rules as $|\frac{1}{2}\rangle|\frac{1}{2}\rangle,|\frac{1}{2}\rangle|\frac{-1}{2}\rangle,|\frac{-1}{2}\rangle|\frac{1}{2}\rangle,|\frac{-1}{2}\rangle|\frac{-1}{2}\rangle$, in which each configuration corresponds to the $\rho=1,2,3,4$-th element of the matrix, respectively. When one operate the transfer matrix from the left-hand side, then the Boltzmann weights coming from an additional site are introduced. Configurations with the inconsistency $s_{j}\neq s^{{}^{\prime}}_{j}$ are eliminated multiplying $0$ element. Then we find the matrix elements $V_{\rho\rho^{\prime}}$ for example $(V)_{11}=\langle\frac{1}{2}|\langle\frac{1}{2}|\>V\>|\frac{1}{2}\rangle|\frac{1}{2}\rangle=xy$, $(V)_{21}=\langle\frac{1}{2}|\langle\frac{-1}{2}|\>V\>|\frac{1}{2}\rangle|\frac{1}{2}\rangle=x/y$, and $(V)_{12}=\langle\frac{1}{2}|\langle\frac{1}{2}|\>V\>|\frac{1}{2}\rangle|\frac{-1}{2}\rangle=0$. The eigenequation is

	$\displaystyle 0$	$\displaystyle=$	$\displaystyle{\rm det\>}(\lambda E-V)$		(46)
		$\displaystyle=$	$\displaystyle\lambda^{4}-2xy\lambda^{3}+y^{2}(x^{2}-\frac{1}{x^{2}})\lambda^{2}+2\frac{y}{x}(y^{2}-\frac{1}{y^{2}})\lambda-(y^{2}-\frac{1}{y^{2}})^{2}$
		$\displaystyle=$	$\displaystyle\Big{[}\lambda^{2}-y(x+\frac{1}{x})\lambda+(y^{2}-\frac{1}{y^{2}})\Big{]}\Big{[}\lambda^{2}-y(x-\frac{1}{x})\lambda-(y^{2}-\frac{1}{y^{2}})\Big{]}.$

The eigenvalues are obtained as

$\displaystyle\lambda_{1}^{\pm}$

$\displaystyle=$

$\displaystyle\frac{1}{2}\Big{[}y(x+\frac{1}{x})\pm\sqrt{y^{2}(x+\frac{1}{x})^{2}-4(y^{2}-\frac{1}{y^{2}})}\Big{]},$

(47)

and

$\displaystyle\lambda_{2}^{\pm}$

$\displaystyle=$

$\displaystyle\frac{1}{2}\Big{[}y(x-\frac{1}{x})\pm\sqrt{y^{2}(x-\frac{1}{x})^{2}+4(y^{2}-\frac{1}{y^{2}})}\Big{]}.$

(48)

When $J_{1}=0$, then $y=1$, and we find $\displaystyle\lambda_{1}^{+}=x+\frac{1}{x}$, $\lambda_{1}^{-}=0$, $\displaystyle\lambda_{2}^{+}=x-\frac{1}{x}$, and $\lambda_{2}^{-}=0$. When $J_{0}=0$, then $x=1$, and we find $\displaystyle\lambda_{1}^{+}=y+\frac{1}{y}$, $\displaystyle\lambda_{1}^{-}=y-\frac{1}{y}$, and $\lambda_{2}^{\pm}=\pm\sqrt{\displaystyle y^{2}-\frac{1}{y^{2}}}$. The eigenvector corresponding to $\lambda_{1}^{+}$ is obtained as

$\displaystyle{\bf u}_{1}=\left(\begin{array}[]{c}1/xy\\ u_{-}\\ u_{-}\\ 1/xy\end{array}\right),$

(53)

where

	$\displaystyle u_{\pm}$	$\displaystyle=$	$\displaystyle\frac{1}{2}\Big{(}R_{-}\pm y\>(x-\frac{1}{x})\Big{)},$
	$\displaystyle R_{-}$	$\displaystyle=$	$\displaystyle\sqrt{y^{2}(x+\frac{1}{x})^{2}-4(y^{2}-\frac{1}{y^{2}})}.$		(54)

The elements of $V$ are all non-negative, and $V$ is irreducible. Hence from the Perron-Frobenius theorem, the maximum eigenvalue of $V$ is unique. Together with the fact that $\lambda_{1}^{+}$ is the maximum eigenvalue for $J_{1}=0$ and for $J_{0}=0$, $\lambda_{1}^{+}$ is also the maximum for all $J_{0}$ and $J_{1}$. Let $U$ be a matrix whose column vectors are the eigenvectors of $V$, the first column being the eigenvector ${\bf u}_{1}$. Then the first row vector ${\bar{\bf u}}_{1}$ of the matrix $U^{-1}$ is obtained as

$\displaystyle{\bar{\bf u}}_{1}=\frac{xy}{2R_{-}}\Big{(}u_{+},\frac{1}{xy},\frac{1}{xy},u_{+}\Big{)}.$

(55)

It is easy to confirm that ${\bar{\bf u}}_{1}\cdot{\bf u}_{1}=1$.

Let us then consider the expectation value in (36) for general spin $S$. The quantum expectation of this type can be calculated using the transfer matrix method.[20] Let $\{|s_{1},\ldots,s_{N}\rangle\}$ be a complete set of the eigenstates of ${\cal H}_{\rm Is}$, in which each $|s_{1},\ldots,s_{N}\rangle=|s_{1}\rangle|s_{2}\rangle\cdots|s_{N}\rangle$ is a direct product of eigenstates of $s_{j}^{z}$ $(j=1,\ldots,N)$. Then, the expectation value $\langle s_{j}^{k}s_{j}^{x}\rangle_{0}$, for example, can be calculated using the transfer matrix $V$ as

	$\displaystyle\langle s_{j}^{k}s_{j}^{x}\rangle_{0}$	$\displaystyle=$	$\displaystyle\frac{{\rm Tr\>}s_{j}^{k}s_{j}^{x}\exp(-\beta{\cal H}_{\rm Is})}{{\rm Tr\>}\exp(-\beta{\cal H}_{\rm Is})}=\frac{\displaystyle\sum_{\{s_{j}\}}\langle s_{1},\ldots,s_{N}|s_{j}^{k}s_{j}^{x}\exp(-\beta{\cal H}_{\rm Is})|s_{1},\ldots,s_{N}\rangle}{\displaystyle\sum_{\{s_{j}\}}\langle s_{1},\ldots,s_{N}|\exp(-\beta{\cal H}_{\rm Is})|s_{1},\ldots,s_{N}\rangle}$		(56)
		$\displaystyle=$	$\displaystyle\frac{{\rm Tr\>}D_{j}^{(p)}V^{N}}{{\rm Tr\>}V^{N}},$

where $D_{j}^{(p)}$ is an operator defined by

$\displaystyle D_{j}^{(p)}=\sum_{s^{\prime}_{j},s^{\prime}_{j+1}}|s^{\prime}_{j}\rangle|s^{\prime}_{j+1}\rangle\Delta_{\rho}^{(p)}\langle s^{\prime}_{j}|\langle s^{\prime}_{j+1}|,$

(57)

and represented by a diagonal matrix $D^{(p)}$ with the elements

	$\displaystyle(D^{(p)})_{\rho\rho^{\prime}}$	$\displaystyle=$	$\displaystyle\Delta_{\rho}^{(p)}\delta_{\rho\rho^{\prime}},$
	$\displaystyle\Delta_{\rho}^{(p)}$	$\displaystyle=$	$\displaystyle\left\{\begin{array}[]{ll}(s^{x}s^{x})_{ll}=(S(S+1)-(S-l+1)^{2})/2&p={\rm even}\\ (s^{y}s^{x})_{ll}=(S-l+1)/2i&p={\rm odd},\end{array}\right.$		(60)

and $\rho$ is the index which corresponds to the spin configuration $(s^{\prime}_{j},s^{\prime}_{j+1})$ of the adjacent two sites $j$ and $j+1$ (see for example below (45)), and $l$ is defined by $s^{\prime}_{j}=S-l+1$, i.e. $\rho$ denotes a configuration $(s^{\prime}_{j},s^{\prime}_{j+1})$, and $l$ is determined from $s^{\prime}_{j}$. In the case of $S=1/2$, for example, $D^{(p)}$ is a diagonal matrix $D^{(p)}={\rm diag\>}(1/4,1/4,1/4,1/4)$ when $p$ is even, and $iD^{(p)}={\rm diag\>}(1/4,1/4,-1/4,-1/4)$ when $p$ is odd. Let us introduce a parameter $m_{\rho}=s^{\prime}_{j}=S-l+1$ which is the eigenvalue satisfying $s_{j}^{z}|s^{\prime}_{j}\rangle=m_{\rho}|s^{\prime}_{j}\rangle$. The parameter $m_{\rho}$ will be used below.

Similarly the correlation function in (36) for general spin $S$ is then expressed, using the transfer matrix (39) and the operator

$\displaystyle M_{j}=\sum_{s^{\prime}_{j},s^{\prime}_{j+1}}|s^{\prime}_{j}\rangle|s^{\prime}_{j+1}\rangle s^{\prime}_{j}\langle s^{\prime}_{j}|\langle s^{\prime}_{j+1}|,$

(61)

which is represented by a diagonal matrix $M={\rm diag}(m_{1},\ldots,m_{n})$, as

			$\displaystyle\langle(s_{j-2}^{z})^{p_{1}-l_{1}}(s_{j-1}^{z})^{p_{0}-l_{0}}(\tau_{p}s_{j}^{k(p)}s_{j}^{x})(s_{j+1}^{z})^{l_{0}}(s_{j+2}^{z})^{l_{1}}\rangle_{0}=\frac{{\rm Tr\>}M^{p_{1}-l_{1}}VM^{p_{0}-l_{0}}V\tau_{p}D^{(p)}VM^{l_{0}}VM^{l_{1}}V^{N-4}}{{\rm Tr\>}V^{N}}$		(70)
		$\displaystyle=$	$\displaystyle\frac{{\rm Tr\>}M^{p_{1}-l_{1}}VM^{p_{0}-l_{0}}V\tau_{p}D^{(p)}VM^{l_{0}}VM^{l_{1}}U\left(\begin{array}[]{cccc}1&&&\\ &(\lambda_{2}/\lambda_{1})^{N-4}&&\\ &&\ddots&\\ &&&(\lambda_{n}/\lambda_{1})^{N-4}\end{array}\right)U^{-1}}{{\rm Tr\>}V^{4}U\left(\begin{array}[]{cccc}1&&&\\ &(\lambda_{2}/\lambda_{1})^{N-4}&&\\ &&\ddots&\\ &&&(\lambda_{n}/\lambda_{1})^{N-4}\end{array}\right)U^{-1}},$

where $U$ is a matrix whose column vectors are the eigenvectors of $V$, the first column being the eigenvector that corresponds to the unique maximum eigenvalue $\lambda_{1}$, and $n=(2S+1)^{2}$. Matrix $U^{-1}VU=\Lambda$ is a diagonal matrix whose diagonal elements are the eigenvalues of $V$. Here let us consider the thermodynamic limit $N\rightarrow\infty$, then we find $(\lambda_{\rho}/\lambda_{1})^{N-4}\to 0$ for $\rho\neq 1$, and (70) is expressed as

$$\rightarrow\frac{1}{c}\sum_{\iota,\>\kappa=1}^{n}\sum_{\mu,\>\rho,\>\nu=1}^{n}{\bar{u}}_{1\iota}m_{\iota}^{p_{1}-l_{1}}V_{\iota\mu}m_{\mu}^{p_{0}-l_{0}}V_{\mu\rho}\tau_{p}\Delta_{\rho}^{(p)}V_{\rho\nu}m_{\nu}^{l_{0}}V_{\nu\kappa}m_{\kappa}^{l_{1}}u_{\kappa 1}\hskip 17.07182pt(N\to\infty),$$

(71)

where

$\displaystyle c=\sum_{\iota,\>\kappa=1}^{n}\sum_{\mu,\>\rho,\>\nu=1}^{n}{\bar{u}}_{1\iota}V_{\iota\mu}V_{\mu\rho}V_{\rho\nu}V_{\nu\kappa}u_{\kappa 1}=\lambda_{1}^{4},$

(72)

and ${\bar{u}}_{1\iota}=(U^{-1})_{1\iota}$, $u_{\kappa 1}=(U)_{\kappa 1}$. The factor ${\bar{u}}_{1\iota}u_{\kappa 1}$ is the exactly renormalized Boltzmann weight with fixed configurations $\iota$ and $\kappa$. When the transfer matrix $V$ cannot be diagonalized, which occur for example when $y^{2}(x+1/x)/2=1$ with $x\neq 1$ in the case of spin $1/2$, we can consider the Jordan normal form, and obtain an identical result. The susceptibility in the thermodynamic limit is obtained, from (9), (36) and (71), as

	$\displaystyle\chi_{0}$	$\displaystyle=$	$\displaystyle\lim_{N\rightarrow\infty}\frac{\chi}{N}$		(80)
		$\displaystyle=$	$\displaystyle\sum_{p=0}^{\infty}\frac{\beta^{p+1}}{(p+1)!}\lim_{N\rightarrow\infty}\frac{1}{N}\langle\delta^{p}{\cal Q}{\cal Q}\rangle_{0}$
		$\displaystyle=$	$\displaystyle\sum_{p=0}^{\infty}\frac{\beta^{p+1}}{(p+1)!}(g\mu_{\rm B})^{2}\sum_{l=0}^{p}\left(\begin{array}[]{c}p\\ l\end{array}\right)(-J_{0})^{p-l}(-J_{1})^{l}$
			$\displaystyle\hskip 22.76228pt\times\sum_{l_{0}=0}^{p_{0}}\left(\begin{array}[]{c}p_{0}\\ l_{0}\end{array}\right)\sum_{l_{1}=0}^{p_{1}}\left(\begin{array}[]{c}p_{1}\\ l_{1}\end{array}\right)\frac{1}{c}\sum_{\iota,\>\kappa=1}^{n}\sum_{\mu,\>\rho,\>\nu=1}^{n}{\bar{u}}_{1\iota}m_{\iota}^{p_{1}-l_{1}}V_{\iota\mu}m_{\mu}^{p_{0}-l_{0}}V_{\mu\rho}\tau_{p}\Delta_{\rho}^{(p)}V_{\rho\nu}m_{\nu}^{l_{0}}V_{\nu\kappa}m_{\kappa}^{l_{1}}u_{\kappa 1}$
		$\displaystyle=$	$\displaystyle\frac{(g\mu_{\rm B})^{2}}{c}\sum_{p=0}^{\infty}\frac{\beta^{p+1}}{(p+1)!}\sum_{l=0}^{p}\left(\begin{array}[]{c}p\\ l\end{array}\right)(-J_{0})^{p-l}(-J_{1})^{l}$
			$\displaystyle\hskip 22.76228pt\times\sum_{\iota,\>\kappa=1}^{n}\sum_{\mu,\>\rho,\>\nu=1}^{n}{\bar{u}}_{1\iota}V_{\iota\mu}V_{\mu\rho}\tau_{p}\Delta_{\rho}^{(p)}V_{\rho\nu}V_{\nu\kappa}u_{\kappa 1}(m_{\mu}+m_{\nu})^{p-l}(m_{\iota}+m_{\kappa})^{l}$
		$\displaystyle=$	$\displaystyle\frac{(g\mu_{\rm B})^{2}}{c}\sum_{\iota,\>\kappa=1}^{n}\sum_{\mu,\>\rho,\>\nu=1}^{n}\sum_{p=0}^{\infty}\frac{\beta^{p+1}}{(p+1)!}\Big{(}(-J_{0})(m_{\mu}+m_{\nu})+(-J_{1})(m_{\iota}+m_{\kappa})\Big{)}^{p}{\bar{u}}_{1\iota}V_{\iota\mu}V_{\mu\rho}\tau_{p}\Delta_{\rho}^{(p)}V_{\rho\nu}V_{\nu\kappa}u_{\kappa 1}.$

The sum should be classified according to $(-J_{0})(m_{\mu}+m_{\nu})+(-J_{1})(m_{\iota}+m_{\kappa})=J(\iota,\mu,\nu,\kappa)$ and $p$. Note that $J(\iota,\mu,\nu,\kappa)^{0}=1$ even if $J(\iota,\mu,\nu,\kappa)=0$. Then we find

	$\displaystyle\chi_{0}$	$\displaystyle=$	$\displaystyle\frac{(g\mu_{\rm B})^{2}}{c}\sum_{\begin{subarray}{c}\iota,\>\kappa,\>\mu,\>\nu,\>\rho=1\\ (J(\iota,\mu,\nu,\kappa)=0)\end{subarray}}^{n}\frac{\beta^{1}}{1!}{\bar{u}}_{1\iota}V_{\iota\mu}V_{\mu\rho}\tau_{0}\Delta_{\rho}^{(0)}V_{\rho\nu}V_{\nu\kappa}u_{\kappa 1}$		(85)
		$\displaystyle+$	$\displaystyle\frac{(g\mu_{\rm B})^{2}}{c}\sum_{\begin{subarray}{c}\iota,\>\kappa,\>\mu,\>\nu,\>\rho=1\\ (J(\iota,\mu,\nu,\kappa)\neq 0)\end{subarray}}^{n}\sum_{\begin{subarray}{c}p=0\\ (p={\rm even})\end{subarray}}^{\infty}\frac{\beta^{p+1}}{(p+1)!}J(\iota,\mu,\nu,\kappa)^{p}{\bar{u}}_{1\iota}V_{\iota\mu}V_{\mu\rho}\Delta_{\rho}^{(\rm even)}V_{\rho\nu}V_{\nu\kappa}u_{\kappa 1}$
		$\displaystyle+$	$\displaystyle\frac{(g\mu_{\rm B})^{2}}{c}\sum_{\begin{subarray}{c}\iota,\>\kappa,\>\mu,\>\nu,\>\rho=1\\ (J(\iota,\mu,\nu,\kappa)\neq 0)\end{subarray}}^{n}\sum_{\begin{subarray}{c}p=1\\ (p={\rm odd})\end{subarray}}^{\infty}\frac{\beta^{p+1}}{(p+1)!}J(\iota,\mu,\nu,\kappa)^{p}{\bar{u}}_{1\iota}V_{\iota\mu}V_{\mu\rho}i\Delta_{\rho}^{(\rm odd)}V_{\rho\nu}V_{\nu\kappa}u_{\kappa 1},$

where $\Delta_{\rho}^{(\rm even)}=\Delta_{\rho}^{(0)}$ and $\Delta_{\rho}^{(\rm odd)}=\Delta_{\rho}^{(1)}$. Finally we obtain

	$\displaystyle\frac{\chi_{0}}{(g\mu_{\rm B})^{2}}$	$\displaystyle=$	$\displaystyle\frac{\beta}{c}\sum_{\begin{subarray}{c}\iota,\>\kappa,\>\mu,\>\nu,\>\rho=1\\ (J(\iota,\mu,\nu,\kappa)=0)\end{subarray}}^{n}{\bar{u}}_{1\iota}V_{\iota\mu}V_{\mu\rho}\Delta_{\rho}^{(0)}V_{\rho\nu}V_{\nu\kappa}u_{\kappa 1}$		(86)
		$\displaystyle+$	$\displaystyle\frac{1}{c}\sum_{\begin{subarray}{c}\iota,\>\kappa,\>\mu,\>\nu,\>\rho=1\\ (J(\iota,\mu,\nu,\kappa)\neq 0)\end{subarray}}^{n}\frac{\sinh\beta J(\iota,\mu,\nu,\kappa)}{J(\iota,\mu,\nu,\kappa)}{\bar{u}}_{1\iota}V_{\iota\mu}V_{\mu\rho}\Delta_{\rho}^{(\rm even)}V_{\rho\nu}V_{\nu\kappa}u_{\kappa 1}$
		$\displaystyle+$	$\displaystyle\frac{1}{c}\sum_{\begin{subarray}{c}\iota,\>\kappa,\>\mu,\>\nu,\>\rho=1\\ (J(\iota,\mu,\nu,\kappa)\neq 0)\end{subarray}}^{n}\frac{\cosh\beta J(\iota,\mu,\nu,\kappa)-1}{J(\iota,\mu,\nu,\kappa)}{\bar{u}}_{1\iota}V_{\iota\mu}V_{\mu\rho}i\Delta_{\rho}^{(\rm odd)}V_{\rho\nu}V_{\nu\kappa}u_{\kappa 1}.$

Equation (86) provides the exact zero-field transverse susceptibility of the spin-$S$ transverse Ising model (1). The result is written in terms of finite summations.

In the case of $S=1/2$, long and straightforward calculations yield that the susceptibility is

	$\displaystyle\frac{\chi_{0}}{(g\mu_{\rm B})^{2}}$	$\displaystyle=$	$\displaystyle\frac{\beta}{4c}\Big{[}(x^{2}+\frac{1}{x^{2}})+\frac{y}{R_{-}}(x+\frac{1}{x})\Big{(}(x-\frac{1}{x})^{2}+\frac{2}{y^{4}}\Big{)}\Big{]}$		(87)
		$\displaystyle+$	$\displaystyle\frac{1}{J_{0}}\frac{1}{2c}\frac{xy}{R_{-}}(x^{2}-\frac{1}{x^{2}})\Big{[}\frac{R_{-}}{xy}+\frac{1}{x}(x+\frac{1}{x})\Big{]}$
		$\displaystyle+$	$\displaystyle\frac{1}{J_{1}}\frac{1}{2c}\frac{xy}{R_{-}}\frac{1}{y^{2}}(y^{2}-\frac{1}{y^{2}})\Big{[}\frac{R_{-}}{xy}+\frac{1}{x}(x+\frac{1}{x})\Big{]}$
		$\displaystyle+$	$\displaystyle\frac{1}{J_{0}-J_{1}}\frac{1}{2c}\frac{xy}{R_{-}}\frac{1}{x^{2}}(x^{2}\frac{1}{y^{2}}-\frac{1}{x^{2}}y^{2})\Big{[}\frac{y}{x}u_{-}+\frac{1}{y^{2}}\Big{]}$
		$\displaystyle+$	$\displaystyle\frac{1}{J_{0}+J_{1}}\frac{1}{2c}\frac{xy}{R_{-}}(x^{2}y^{2}-\frac{1}{x^{2}y^{2}})\Big{[}xyu_{+}+\frac{1}{y^{2}}\Big{]},$

where

	$\displaystyle c$	$\displaystyle=$	$\displaystyle\frac{1}{2}\Big{(}y^{4}(x^{4}+\frac{1}{x^{4}})+4(x^{2}+\frac{1}{x^{2}})+4+\frac{2}{y^{4}}\Big{)}$		(88)
		$\displaystyle+$	$\displaystyle\frac{1}{2}\frac{y}{R_{-}}(x+\frac{1}{x})\Big{(}(x^{2}+\frac{1}{x^{2}})+\frac{2}{y^{4}}\Big{)}\Big{(}(x-\frac{1}{x})^{2}y^{4}+4\Big{)}$
		$\displaystyle=$	$\displaystyle(\lambda_{1}^{+})^{4},$

and $x=e^{K_{0}/4}$, $y=e^{K_{1}/4}$, and $u_{\pm}$ and $R_{-}$ are defined in (54). Generally in the case of spin $S$, we have to find the eigenvectors of $n\times n=(2S+1)^{2}\times(2S+1)^{2}$ matrix $V$.