How to find the reduced density matrix of a four-qubit system?

Question

I have the state vector $|p\rangle$ made up of 4 qubits. Say system A is made up of the first and second qubits while system B is made up of qubits 3 and 4. I want to find the reduced density matrix of system A.

I know I could separately extract qubits 1,2 and 3,4 into their own state vectors then find their density matrices and compute the reduced density matrix for system A.

I want to figure out how to do this without having to extract and separate the systems. First I would find the density matrix of $|p\rangle$ and then do a partial trace with respect to system B. I am not sure how to do the partial trace of system B since the system contains 2 qubits.

Can anyone help me figure this out? I am using Python and NumPy for reference.

Answer 1

Wikipedia entry has a nice description of the partial trace and how to compute it.

In your case the partial trace of $M=|p\rangle\langle p|$ over $B$ can be computed as $$ \text{Tr}_B(M) = \text{Tr}_B\bigg(\sum_{k,l} M_{k,l} \otimes |k\rangle\langle l| \bigg) = \sum_{k} M_{k,k} $$
where $k,l \in \{00,01,10,11\}$.

Those $M_{k,l}$ are $4 \times 4$ submatrices of the $16 \times 16$ matrix $M$. What kind of submatrices depends on the actual encoding. In the big-endian encoding – where the matrix row order corresponds to $|0000\rangle, |0001\rangle, ... , |1111\rangle$ – matrix $M_{k,l}$ has element in the position $(i,j)$ that equals to an element in the position $(k,l)$ of the $(i,j)$-th block of the big matrix $M$. So the resulting sum is actually equivalent to a matrix whose elements are traces of $4 \times 4$ blocks of the matrix $M$.

Partial trace over the first subsystem A in the big-endian encoding is easier – it's just the sum of 4 blocks from the diagonal.

Answer 2

Computationally, the easiest way to do this is probably as follows:

Let your state be $$ |\psi\rangle=\sum_{i,j,k,l}c_{ijkl}|ij\rangle_A|kl\rangle_B $$ Rewrite this as a matrix $$ C=\sum_{i,j,k,l}c_{ijkl}|ij\rangle\langle kl| $$ Effectively, you just have to reshape your numpy array.

Then, you can calculate $$ \rho_A=CC^\dagger $$ or $$ \rho_B=C^\dagger C. $$