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.
 A: XX (and more generally XY) spin models became very popular lately, especially in connection to quantum-information theoretic problems. The XY model (without magnetic field though) has been introduced and properly solved in 1961 by Lieb, Schultz and Mattis (Ann. Phys. 16, 407, (1961)). By the way this is a great reference to start with, if you didn't read this paper already. Nonetheless, in many recent articles, this model is solved incorrectly, and this might be the origin of the discrepancy that you observed. 
Let me address your point.
First I imagine that the average that you mean is the ground state average: i.e. you digonalize the Hamiltonian, take the eigenvector related to the lowest eigenvalue, and compute statistical averages with this vector. 
This said, to specify completely the problem you must assign some boundary conditions (BC). For many reasons (which I won't enter here), the most appropriate boundary condition is the periodic one (PBC) for the spins. With this I mean $\sigma^\alpha_{L+1}=\sigma^\alpha_1$, $L$ being the system size. When you map the spin problem to a Fermionic one through the Jordan-Wigner (JW) mapping, the Hamiltonian becomes that of free Fermions plus a peculiar border term. This is because the JW string accumulates in the bond which couples spins 1 and L. The result is that an XY (or XX) spin model with PBC is equivalent to a free-fermion model with a parity dependent boundary conditions (the parity operator is $P\sim \prod \sigma^z_i \sim (-1)^N$, and $N$ is the total number operator in the Fermi picture, and the $\sim$ symbol ignores potential factors of $i$). In the parity 1 sector (which comprises the ground state) BC are antiperiodic (i.e. $c_{L+1}=-c_1$ for the Fermi operators $c_i$). This boils down to antiperiodic quasimomenta when you do the Fourier transform: $k_n = \pi (2n+1)/L$ with $n=0,1,\ldots,L-1$. 
If you instead choose open boundary conditions, there is no funky border term in the fermionic problem either, but the diagonalization of an open chain is slightly more cumbersome. 
This could be one source of errors.
Another potential source of error is degeneracy (but this shouldn't be your case). When you diagonalize the Hamiltonian you should check that the ground state is unique. In a quasi-free systems you have degenerate ground states when the one-particle dispersion becomes exactly zero for some quasimomentum. This for instance, is known to happen in a pure free fermionic chain with nearest neighbour interaction and PBC,  when the size $L$ is a multiple of 4. 
As I said, a great reference is Lieb, Schultz and Mattis article, where the precise  diagonalization procedure is presented. Boundary conditions are also discussed in a recent paper of mine (PRA 81, 060101 (2010)) or in the book of M. Henkel, Conformal invariance and critical phenomena.
Hope this helps.
Best,
Lorenzo
A: Wick's theorem works for an arbitrary number of states equipped with their creation and annihilation oscillators. In the same way, Wick's theorem for spins applies to a collection of an arbitrary number of spins - or qubits, which is the same thing if their spin is $1/2$. Wick's theorem is just a mathematical identity that works even for - and especially for - a finite number of qubits and/or number of creation and annihilation operators.
So it can't be the case that the source of the discrepancy hides in Wick's theorem. However, to uncover the reason of the discrepancy, you would probably have to be more specific about what you calculated (the expectation value in what state? Why isn't the original form of the expectation value already good as a "result"? In which form do you want the result?) for others to isolate the mistake or subtlety. Or others could even calculate the right result, avoiding mistakes - but it's not clear from the text what exactly you want them to calculate.
Best wishes
Luboš
A: As Lubos mentions Wick's Theorem is true for any number of particles. This is simply because it is giving a different way to express a two qubit operator, and hence has nothing to do with what system it is applied to.
If the magnetic field is 0 and the couplings are constant, then an XX spinchain can be diagonalized by a Fourier transform, so you may want to try the calculation via this approach. This will only work for a ring and not a linear chain.
I believe that the problem you are having is indeed due to the finite length of the chain, but is hiding in a different place. Let me explain why: For a 2 qubit chain, the correlation operator commutes with the Hamiltonian, however this is not true for longer chains, and so we would expect to obtain different results for the $R=1$ correlation in the case of a two qubit chain and the case of a 3 qubit chain. Therefore, it seems clear that you are burying the chain length somewhere. If I had to guess, I would say it is probably an edge effect that has been skipped in doing the Jordan-Wigner transformation or in your explicit calculation, but obviously I can't be sure with out seeing both calculations explicitly.
Hope this is useful.
