How can you (computationally) calculate the halfchain entanglement entropy of a spinchain?

Question

I am simulating a (small) spinchain with exact diagonalization and dynamics. I would like to track the entanglement entropy of half the chain with the other part of the chain.

I have the vectors of my state $|\Psi(t)>$ in the basis $|s_1, s_2, .., s_i, .. , s_n>$ which are the spins at each sites at a certain time. I know you can get the density matrix of this by taking the outer product $\rho = |\Psi(t)>< \Psi(t)|$ which is a $n \times n$ matrix, but this is a pure state so entanglement entropy will be zero, so now I want the density matrix of the subsystem which includes the first (or last) half of the spin chain. This probably includes tracing out certain parts of the density matrix $\rho$ to $\rho_A$, but how do I do this specifically, what parts?

This is specifically a computational question because I am simulating in Python, but it's not the coding part that's troublesome, I really have no idea how to get the "first half of the chain" out of this full density matrix.

Answer 1

It all boils down to a two-partite system, $$ \vert\psi\rangle = \sum c_{ij} \vert i,j\rangle\ , $$ where $i$ and $j$ are all indices in the left and right part, respectively. Then, the reduced density matrix is given by $$ \rho = \sum_{ii'} \left(\sum_j c_{ij}c^*_{i'j}\right)\vert i\rangle\langle i'\vert\ . $$ That is, in the standard basis $\rho = CC^\dagger$, where $C=(c_{ij})$ is the coefficients of the pure state $\vert\psi\rangle$, arranged in a matrix.

So computationally, what you want to do is to take your coefficient vector $c_{ij}$, reshape it into a matrix $C$ (this comes at zero computational cost as it just changes the way the object is indexed), and then compute its singular values. The square of those singular values are then the eigenvalues of $CC^\dagger$, and from those, you can easily compute the entropy.