quantum mechanics - To calculate the correlation functions of an XX spin chain, Wick's theorem is used. But is it valid for a chain of any size?

The correlation functions found in Barouch and McCoy's paper (PRA 3, 2137 (1971)) for the XX spin chain use a method which uses Wick's theorem. For the zz correlation function, this gives

$\langle \sigma_l^z \sigma_{l+R}^z \rangle = \langle \sigma_l^z \rangle^2 - G_R^2$

where for $R=1$, $G_1 = -\langle \sigma_l^x \sigma_{l+1}^x+ \sigma_l^y \sigma_{l+1}^y \rangle/2$.

If I calculate $\langle \sigma_l^z \sigma_{l+1}^z \rangle$ both explicitly and using the equation above for 8 qubits, I get different answers.

So is Wick's theorem still valid for 8 qubits, which means I've just made a mistake? Or is it valid only in the thermodynamic limit?

Thanks

Edit:

Thanks for your replies everyone. @lcv However, I haven't used the analytical diagonalisation for this - I have simply used Mathematica to diagonalise the 8 qubit chain numerically after substituting arbitrary values for the coupling strength, magnetic field and temperature. Hence it can't be an error in the diagonalisation. It is the thermal average I have calculated, that is $\langle \sigma^z_l \rangle=tr(\rho \sigma^z_l )$ where $\rho=e^{−H/T}/tr(e^{−H/T})$ and T is temperature. But in doing this, I find that $\langle \sigma^z_l \sigma^z_{l+R} \rangle \neq \langle \sigma^z_l \rangle^2 - G_1^2$ where I've defined $G_1$ above.

Edit2 (@marek @lcv @Fitzsimons @Luboš) I'm going to try to clarify - The open XX Hamiltonian in a magnetic field is

$$\begin{equation} H=-\frac{J}{2}\sum_{l=1}^{N-1} (\sigma^x_l \sigma^x_{l+1} + \sigma^y_l \sigma^y_{l+1})- B \sum_{l=1}^N \sigma^z_l \end{equation}$$

In Mathematica, I have defined the Pauli spin matrices, then the Hamiltonian for 8 qubits. I then put in values for $J$, $B$ and $T$, and calculate the thermal density matrix,

$$\begin{equation} \rho = \frac{e^{-H/T}}{tr(e^{-H/T})} \end{equation}$$

So now I have numerical density matrix. I then calculate $\langle \sigma^z_l \sigma_{l+1}^z \rangle=tr(\rho \sigma^z_l \sigma_{l+1}^z )$ using the definitions of the Pauli spin matrices and $\rho$.

Next I calculate $\langle \sigma_l^z \sigma_{l+R}^z \rangle$ using the result from Wick's theorem which gives $\langle \sigma_l^z \rangle^2 - G_R^2$ where for $R=1$, $G_1 = -\langle \sigma_l^x \sigma_{l+1}^x+ \sigma_l^y \sigma_{l+1}^y \rangle/2$. I again use the Pauli spin matrices I defined and the same numerical $\rho$ to calculate them.

But I get a different (numerical) answer for each of these.