The paper Osborne and Nielsen, Phys. Rev. A 66, 032110 gives an exact solution for the two-site reduced density matrix for the transverse-field Ising model with periodic boundary conditions (Eq. 26): $$ \rho_{0r} =\frac{I_{0r}+\left<\sigma^z\right>\left(\sigma_0^z+\sigma_r^z\right)+\sum_{k=1}^3\left<\sigma_0^k\sigma_r^k\right>\sigma_0^k\sigma_r^k}{4} $$ where $0r$ indicates that site 0 and $r$ are in subsystem A and the rest of the infinite chain in subsystem B.
Question: When I evaluate this at criticality and compute the entanglement entropy $S_{0r}=Tr(\rho_{0r}\log{\rho_{0r}})$, I obtain a complex number... How does one set up this computation correctly? (see below)
The correlation functions at criticality are known (Eqs. 49-52): $$ \begin{align} \left<\sigma_0^x\sigma_r^x\right>&=\left(\frac{2}{\pi}\right)^r 2^{2r(r-1)}\frac{H(r)^4}{H(2r)}\\ \left<\sigma_0^y\sigma_r^y\right>&=-\frac{\left<\sigma_0^x\sigma_r^x\right>}{4r^2-1}\\ \left<\sigma_0^z\sigma_r^z\right>&=\frac{4}{\pi^2}\frac{1}{4r^2-1}\\ \left<\sigma^z\right>&=\frac{2}{\pi} \end{align} $$ where $H(r)=1^{r-1}2^{r-2}\dots (r-1)$.
Evaluating at $r=1$:
\begin{align}
\left<\sigma_0^x\sigma_1^x\right>&=\frac{2}{\pi}\approx0.63662\\
\left<\sigma_0^y\sigma_1^y\right>&=-\frac{2}{3\pi}\approx -0.212207\\
\left<\sigma_0^z\sigma_1^z\right>&=\frac{4}{3\pi^2}\approx 0.135095\\
\left<\sigma^z\right>&=\frac{2}{\pi}
\end{align}
These values match Pfeuty (1970) (where $n=r$):
So far so good. Using the usual $2\times 2$ identity matrix and Pauli matrices $\sigma^k$, I use the Kronecker product to construct the full matrices: $$ \begin{align} \sigma_0^k &= \sigma^k \otimes I\\ \sigma_r^k &= I \otimes \sigma^k\\ I_{0r} &= I\otimes I \end{align} $$
Finally, the reduced density matrix is: $$ \begin{align} \rho_{01} &=\frac{I_{01}+\left<\sigma^z\right>\left(\sigma_0^z+\sigma_1^z\right)+\sum_{k=1}^3\left<\sigma_0^k\sigma_1^k\right>\sigma_0^k\sigma_1^k}{4}\\ &=\begin{pmatrix} \frac{1}{4}+\frac{1}{3\pi^2}+\frac{1}{\pi} & 0 & 0 & \frac{2}{3\pi}\\ 0 & \frac{1}{4}-\frac{1}{3\pi^2} & \frac{1}{3\pi} & 0\\ 0 & \frac{1}{3\pi} & \frac{1}{4}-\frac{1}{3\pi^2} & 0\\ \frac{2}{3\pi} & 0 & 0 & \frac{4+3(\pi-4)\pi}{12\pi^2} \end{pmatrix} \end{align} $$ Note that $\rho^*=\rho$ and $Tr(\rho)=1$, so that is good. However, the two-site entanglement entropy is now $$ S_{01}=Tr(\rho_{01}\log{\rho_{01}})\approx0.468674 + 0.223869 i. $$
What is going wrong here?
I noticed that the references Pfeuty and Osborne/Nielsen differ by a factor of $1/\pi$ for $\left<\sigma_0^z\sigma_r^z\right>$, which appears to be a typo, but both these values lead to a complex number for $S_{01}$, so the problem is not caused by that.
This is a follow-up question on Entanglement entropy of the infinite transverse field Ising chain at the critical point, but it is a separate technical problem, hence the new post here.
