I just recently completed a problem in which I had a Hilbert space of the form
$$
H = H_1 \otimes H_2 \otimes H_3
$$
and was tasked with finding the reduced density matrix for the system in the subspace
$$
H^\prime = H_1 \otimes H_2
$$
for some explicit states on an orthonormal basis. For our purposes, we can assume the orthonormal basis to be that of a qubit. I was able to determine that, for a state $|\psi\rangle$, the reduced density matrix follows from computing
$$
\hat{\rho}_{\psi^\prime} = \mathrm{Tr}_3(\hat{\rho}_\psi).$$
For example, if we have
$$
|\psi\rangle = |0\rangle_1\otimes |0\rangle_2 \otimes |0\rangle_3,$$
our density matrix is
$$
\hat{\rho}_\psi = |0\rangle_1\otimes |0\rangle_2 \otimes |0\rangle_3 \;\; \langle 0|_1 \otimes \langle 0|_2 \otimes \langle 0|_3
$$
such that, when we compute the trace in the $H_3$ basis with
$$
\mathrm{Tr}_3(\hat{\rho}_\psi) = \langle 0|_3 \hat{\rho}_\psi|0\rangle_3 + \langle 1|_3 \hat{\rho}_\psi|1\rangle_3,
$$
we get
$$
\hat{\rho}_{\psi^\prime} = |0\rangle_1 \otimes |0\rangle_2 \;\; \langle 0|_1 \otimes \langle 0|_2.
$$
I am quite comfortable with this and wrote things explicitly using tensor products to keep everything nice and tidy. From here, I know how to compute $\mathrm{Tr}_1(\hat{\rho}_\psi)$ and $\mathrm{Tr}_2(\hat{\rho}_\psi)$.
Moreover, I can extrapolate the above procedure to an $n-$dimensional Hilbert space whereby I can compute the reduced density matrix for every $(n-1)-$dimensional subspace. All I need to do is compute the trace in the $H_i$ basis for whatever value of $i$ I choose. However, here is where I am puzzling over how to generalize my procedure.
Suppose we now have
$$H^{\prime \prime} = H_2.$$
Would we compute the reduced density matrix for a system in this subspace via
$$\hat{\rho}_{\psi^{\prime \prime}} = \mathrm{Tr}_1(\mathrm{Tr}_3(\hat{\rho}_\psi)) = \mathrm{Tr}_3(\mathrm{Tr}_1(\hat{\rho}_\psi))?$$
Generalizing further, suppose our initial Hilbert space is now
$$
H = H_1 \otimes H_2 \dots \otimes H_n
$$
for a finite natural number $n$ and we'd like to find the reduced density matrix for a finite subspace
$$
H^\prime = H_i \otimes H_{j}\otimes \dots \otimes H_k
$$
where $i,j,k$, and all numbers in between constitute a finite subset of the set $\{1,\dots,n\}.$ Would it be correct to infer from my above considerations that
$$
\hat{\rho}_{\psi^\prime} = \mathrm{Tr}_i(\mathrm{Tr}_j\dots(\mathrm{Tr}_k(\hat{\rho}_\psi))) = \mathrm{Tr}_{i,j,\dots,k}(\hat{\rho}_\psi)?
$$
Thanks in advance!
Best Answer
I would change your notation for product operators. It's much easier if you write the product-state (i.e., separable) density matrix as $\rho = \rho^{\,}_1 \otimes \rho^{\,}_2 \otimes \rho^{\,}_3$, so for the simplest case above, it's $\rho = \left| 0 \middle\rangle \hspace{-0.4mm} \middle\langle 0 \right|^{\,}_1 \otimes \left| 0 \middle\rangle \hspace{-0.4mm} \middle\langle 0 \right|^{\,}_2 \otimes \left| 0 \middle\rangle \hspace{-0.4mm} \middle\langle 0 \right|^{\,}_3$, and this will make it easier to think about taking the trace over one of the subsystems (e.g., the third one).
The order of the traces is unimportant (the trace operation commutes). In general, suppose you have a Hilbert space $\mathcal{H}$, and suppose that you can divide this into two regions, $A$ and $A^c$, with Hilbert spaces $\mathcal{H}^{\,}_A$ and $\mathcal{H}^{\,}_{A^c}$ such that $\mathcal{H} = \mathcal{H}^{\,}_A \cup \mathcal{H}^{\,}_{A^c}$. Then if you want the reduced density matrix for region $A$, for any subsystem $A$, you simply take the trace over all degrees of freedom in $A^c$ (the complement of $A$). If $A^c$ contains multiple distinguishable regions, trace over all of them in any order.
Note that the cases above correspond to separable density matrices, whose partial trace is trivial. We can extend this to generic separable density matrices of the form $$\rho = \sum_{k} p^{\,}_k \, \rho^{\,}_{A,k} \otimes \rho^{\,}_{A^c,k}$$ so that the reduced density matrix for region $A$ is $$\rho^{\,}_A = \underset{A^c}{\rm tr} \left[ \, \rho \, \right] = \sum_{k} p^{\,}_k \, \underset{A^c}{\rm tr} \left[ \, \rho^{\,}_{A^c,k} \, \right] \, \rho^{\,}_{A,k} = \sum_k \, p^{\prime}_k \, \rho^{\,}_{A,k} ,$$ where $\sum_k p_k' = 1$ for normalization.
But the notation I suggest is helpful in going to the generic case. In general, we can always choose an operator basis that is separable, and decompose the density matrix in that basis. For a system of $N$ qubits, a simple operator basis is the set all $4^N$ "Pauli strings", which are the operators that act on every site as either the identity or one of the three Paulis. Denoting these strings by $\sigma^{\vec{\mu}} = \bigotimes_{j=1}^N \, \sigma^{\mu_j}_j$, where $\mu=0$ corresponds to the identity, we have $$ \rho \, = \, \sum\limits_{\vec{\mu}} \, p^{\,}_{\vec{\mu}} \, \bigotimes\limits_{j=1}^{N} \, \sigma^{\mu^{\,}_j}_j \, ,$$ where the coefficients are given by $$p^{\,}_{\vec{\mu}} \, = \, \frac{1}{2^N} \, {\rm tr} \left[ \, \rho \, \bigotimes\limits_{j=1}^{N} \, \sigma^{\mu^{\,}_j}_j \, \right] \, ,~$$ and then taking the trace over $j \in A^c$ to get $\rho^{\,}_A$ is straightforward: $$ \rho^{\vphantom{\dagger}}_A \, = \, \underset{A^c}{\rm tr} \left[ \, \rho \, \right] \, = \, \sum\limits_{\{ \mu^{\,}_j \}} \, p^{\,}_{\vec{\mu}} \, \bigotimes\limits_{j \in A} \, \sigma^{\mu^{\,}_j}_j \times \bigotimes\limits_{k \in A^c} {\rm tr} \left[ \, \sigma^{\mu^{\,}_k}_k \, \right] \, = \, \sum\limits_{\{ \mu^{\,}_j \}} \, p^{\,}_{\vec{\mu}} \, \bigotimes\limits_{j \in A} \, \sigma^{\mu^{\,}_j}_j \, \left[ \prod\limits_{k \in A^c} \, \delta^{\vphantom{\prime}}_{\mu^{\,}_k,0} \right] \, , ~~$$ so that the trace over $A^c$ merely restricts the sum over Pauli strings to those that act as the identity $\mathbb{1}$ on every site in $A^c$. Note that similar operator bases exist for qudits with $q>2$ states per site; in general, we choose a basis in which all nonidentity operators are traceless, so that $\rho^{\,}_A$ follows from restricting the sum over basis operators that defines the full density matrix $\rho$ to the basis operators that act trivially (as $\mathbb{1}$) on all degrees of freedom in $A^c$.