Background
Consider 1-D Ising model of n lattice points with periodic boundary condition,
$\beta H(\sigma_1,\sigma_2,…,\sigma_N) = -\sum_{i=1}^nk(\sigma_i\sigma_{i+1})-\sum_{i=1}^n\sigma_i$
$k=\beta J$ and $h=\beta B$ where $J$ is the coupling constant and $B$ is the applied magnetic field.
Partition function, $Z=\sum_{\sigma_1,\sigma_2,…,\sigma_N}e^{-\beta H(\sigma_1,\sigma_2,…,\sigma_N)}$
where $\sum_k$ is the sum over all the possible microstates and $H_k$ is the hamiltonian of the $k-th$ microstate.
$Z = \sum_{\sigma_1,\sigma_2,…,\sigma_N}exp\Big(k(\sigma_1\sigma_2+\sigma_2\sigma_3+…+\sigma_N\sigma_1)\Big)exp\Big(h(\sigma_1+\sigma_2+…+\sigma_N)\Big)$
$\implies Z = \sum_{\sigma_1,\sigma_2,…,\sigma_N}exp\Big(k\sigma_1\sigma_2+\frac{(\sigma_1+\sigma_2)}{2}\Big)exp\Big(k\sigma_2\sigma_3+\frac{(\sigma_2+\sigma_3)}{2}\Big)…exp\Big(k\sigma_N\sigma_1+\frac{(\sigma_N+\sigma_1)}{2}\Big)$
The $Z$ can be written in terms of transfer matrix as
$Z=\sum_{\sigma_1,\sigma_2,…,\sigma_N}\langle\sigma_1|T|\sigma_2\rangle\langle\sigma_2|T|\sigma_3\rangle…\langle\sigma_N|T|\sigma_1\rangle$
where $T=\begin{pmatrix}e^{k-h}&e^{-k}\\e^{-k}&e^{k+h}\end{pmatrix}$
$\sigma_1=-1$ corresponds to $|\sigma_1\rangle=\begin{pmatrix}1\\0\end{pmatrix}$ and $\sigma_1=+1$ corresponds to $|\sigma_1\rangle=\begin{pmatrix}0\\1\end{pmatrix}$
So, $Z = \sum_{\sigma_1}\langle\sigma_1|T^N|\sigma_1\rangle=\sum_{\sigma_1}\langle O\sigma_1|OT^NO^T|O\sigma_1\rangle\tag{1}$
where $O$ is the orthogonal matrix that diagonalize $T$ and thus $T^N$.
$OTO^T=\begin{pmatrix}\lambda_1&0\\0&\lambda_2\end{pmatrix}$
where $\lambda_{1,2}=e^k\cosh h\pm\sqrt{e^{2k}\sinh^2 h+e^{-2k}}$
So, $Z=Tr(T^N)=Tr(OT^NO^T)=\lambda_1^N+\lambda_2^N\tag{2}$
Work
Now, I am trying to calculate $\langle\sigma_i\rangle$ when $h=0$ using transfer matrix.
We know that when $h=0$, $\langle\sigma_i\rangle=0$
$\langle\sigma_i\rangle=\sum_{\sigma_1,\sigma_2,…,\sigma_N}\langle\sigma_1|T|\sigma_2\rangle\langle\sigma_2|T|\sigma_3\rangle…\langle\sigma_{i-1}|T|\sigma_i\rangle\sigma_i\langle\sigma_i|T|\sigma_{i+1}\rangle…\langle\sigma_N|T|\sigma_1\rangle$
$\langle\sigma_i\rangle=\sum_{\sigma_1,\sigma_i}\langle\sigma_1|T^{i-1}|\sigma_i\rangle\sigma_i\langle\sigma_i|T^{N-i+1}|\sigma_1\rangle$
As quantities inside $\sum$ are just numbers (matrix elements and $\sigma_i$ are numbers),
So, $\langle\sigma_i\rangle=\sum_{\sigma_1,\sigma_i}\sigma_i\langle\sigma_i|T^{N-i+1}|\sigma_1\rangle\langle\sigma_1|T^{i-1}|\sigma_i\rangle$
$\implies\langle\sigma_i\rangle=\sum_{\sigma_i}\sigma_i\langle\sigma_i|T^{N-i+1}\sum_{\sigma_1}\Big(|\sigma_1\rangle\langle\sigma_1|\Big)T^{i-1}|\sigma_i\rangle$
$\implies\langle\sigma_i\rangle=\sum_{\sigma_i}\sigma_i\langle\sigma_i|T^{N-i+1}T^{i-1}|\sigma_i\rangle$
$\implies\langle\sigma_i\rangle=\sum_{\sigma_i}\sigma_i\langle\sigma_i|T^N|\sigma_i\rangle$
$\implies\langle\sigma_i\rangle=\sum_{\sigma_i}\sigma_i\langle O\sigma_i|OT^NO^T|O\sigma_i\rangle$
$\implies\langle\sigma_i\rangle=\sum_{\sigma_i}\sigma_i\langle O\sigma_i|\begin{pmatrix}\lambda_1^N&0\\0&\lambda_2^N\end{pmatrix}|O\sigma_i\rangle\tag{3}$
Doubt
But now how to show that $\langle\sigma_i\rangle=0$?
One way is to calculate $O$ explicitly and plug this in $(3)$ and show that it is $0$. This I have verified.
But I have seen that in books, they introduce Pauli matrix to show that it is $0$. I am not able to understand how that follows from $(3)$.
Best Answer
As I explained in the comments, you don't need to introduce any Pauli matrix to see that $\langle\sigma_i\rangle=0$. In fact, you don't even need to introduce the transfer matrix at all, but just observe that the Gibbs measure gives the same probability to the configurations $(\sigma_1,\dots,\sigma_N)$ and $(-\sigma_1,\dots,-\sigma_N)$.
So, I'll rather answer the underlying technical question that is spelled out in your comments, namely how to rewrite your final expression in terms of one of the Pauli matrices. (This is actually useful if you want to compute, for instance, the 2-point function; see this answer for such an application.)
The idea is very simple, just rewrite \begin{align*} \sum_{\sigma_i=\pm 1} \sigma_i (T^N)_{\sigma_i,\sigma_i} &= \sum_{\sigma_i=\pm 1} \sum_{\sigma'_i=\pm 1} \delta_{\sigma_i\sigma'_i} \sigma_i (T^N)_{\sigma_i,\sigma_i}\\ &= \sum_{\sigma_i=\pm 1} \sum_{\sigma'_i=\pm 1} S_{\sigma_i,\sigma'_i}(T^N)_{\sigma'_i,\sigma_i}\\ &= \sum_{\sigma_i=\pm 1} (ST^N)_{\sigma_i,\sigma_i}, \end{align*} where I have introduced the matrix $S=\begin{pmatrix}1&0\\0&-1\end{pmatrix}$, using the fact that $\delta_{\sigma_i,\sigma'_i}\sigma_i = S_{\sigma_i,\sigma'_i}$. In the second identity, I use the fact that the summand vanishes when $\sigma_i\neq\sigma'_i$ to replace the second instance of $\sigma_i$ by $\sigma'_i$.