Entanglement entropy of the infinite transverse field Ising chain at the critical point Consider the $1$D transverse-field Ising model,
$$H = -\sum_i \sigma_i^z \sigma_{i+1}^z-h\sum_i \sigma_{i}^x$$
at the critical point $h=1$ with periodic boundary conditions.
Question: What is the (von Neumann) entanglement entropy of the ground state, taking a finite number of sites in subsystem $A$ as the chain length $L\rightarrow \infty$?
In particular, I would like to know the single spin and two neighbouring spin cases, i.e. $S_1$ and $S_{12}$. In these cases, one or two sites are in subsystem $A$ and the rest is in the (infinitely long) subsystem $B$, respectively. In the case of two sites in $A$, I assume they are consecutive (potentially simplifies the problem).
Using exact diagonalization and the symmetries of the Hamiltonian, it is relatively easy to plot $S_1/\log(2)$ and $S_{12}/(2\log(2))$ for $L$ up to $26$:


These entanglement entropies seem to converge to some value, but the question is to what value exactly?
I know conformal field theory gives a result for half-chain entropy (e.g. Understanding entanglement entropy in the transverse field Ising model), or when subsystem $A$ is a fraction of the total length $L$, but this result does not give quantitative answers in the finite case of just one or two spins, right?
 A: These few-site entropies tend to converge fairly fast with finite size, because they ultimately depend on very local physics. If you switch to DMRG you should easily see the values converge fully. You can also analytically obtain the single- and two-site reduced density matrices of transverse-field Ising and XY chains in the thermodynamic limit, see Osborne and Nielsen, Phys. Rev. A 66, 032110 (2002). In their notation the Ising chain reads [Eq. (1) with $\gamma=1$]
$$
H = -\sum_j \left( \lambda \sigma^x_j \sigma^x_{j+1} + \sigma_j^z  \right),
$$
so the case $\lambda=1$ corresponds to your $h=1$.
The one-site reduced density matrix is given by [Eq. (25)]
$$
\rho_1 = \frac{I+\langle \sigma^z\rangle \sigma^z}{2},
$$
where in the ground state [Eq. (35)]
$$
\langle \sigma_z \rangle = \frac{1}{\pi} \int_0^\pi \mathrm{d}\phi \frac{1+\lambda\cos \phi }{\sqrt{1+\lambda^2+2\lambda \cos\phi}}.
$$
Evaluating the entropy at $\lambda=1$, where $\langle \sigma^z\rangle = 2/\pi$, one finds
$$
\begin{align}
S_1 &= -Tr \left[ \rho_1 \log \rho_1 \right] / \log{2}\\
&= \frac{\pi  \left(\log (4)+2 \log (\pi )-\log \left(\pi ^2-4\right)\right)-4 \tanh ^{-1}\left(\frac{2}{\pi }\right)}{\pi  \log (4)} \approx 0.68376.
\end{align}
$$
Your largest system size gets quite close!
The two-site reduced density matrix for two sites separated by a distance $r$ is given by [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}.
$$
The expression involves spin-spin correlations $\langle \sigma_0^k \sigma_r^k\rangle$, which are known at criticalty ($\lambda=1$), see Pfeuty, Ann. Phys. (N.Y.) 57, 79-90 (1970) [Eqs. (3.3)-(3.5)].
$$
\begin{align}
\langle \sigma_0^x\sigma_r^x\rangle_c &=\left(\frac{2}{\pi}\right)^r 2^{2r(r-1)}\frac{H(r)^4}{H(2r)}\\
\langle\sigma_0^y\sigma_r^y\rangle_c &=-\frac{\langle \sigma_0^x\sigma_r^x\rangle }{4r^2-1}\\
\langle\sigma_0^z\sigma_r^z\rangle_c&=\frac{4}{\pi^2}\frac{1}{4r^2-1}
\end{align}
$$
where $H(r)=1^{r-1}2^{r-2}\dots (r-1)$. (Osborne and Nielsen reproduce these as Eqs. (49)-(51), but lack a factor of $1/\pi$ for $\langle \sigma^z_0 \sigma^z_r\rangle_c$.) The subscript ${}_c$ denotes that these are connected correlation functions, i.e.
$$
\langle \sigma_0^k \sigma_r^k \rangle_c = \langle \sigma_0^k \sigma_r^k \rangle - \langle \sigma_0^k\rangle \langle \sigma_r^k\rangle,
$$
but note that the two-site reduced density matrix uses the full expression $\langle \sigma_0^k \sigma_r^k \rangle$, see the answer to this related question. At criticality only $\langle \sigma^z\rangle$ is non-zero, so the disconnected part only matters for the $zz$-correlation.
Approximate values for the nearest-neighbor correlations are $\langle \sigma_0^x\sigma_1^x\rangle_c\approx 0.63662$, $\langle \sigma_0^y\sigma_1^y\rangle_c\approx -0.212207$, and $\langle \sigma_0^z\sigma_1^z\rangle_c\approx 0.135095 \quad \Rightarrow \quad \langle \sigma_0^z\sigma_1^z\rangle\approx 0.54038$.
Evaluation of the nearest-neighbors two-site entropy with normalization as in the question then yields
$$
S_{12} = -Tr \left[ \rho_{01} \log \rho_{01} \right]/(2\log 2) \approx 0.428054,
$$
which is also consistent with the numerical result.
