A Notebook for Subsystems and Entanglement

Small numerical examples make tensor-product conventions tangible. The downloadable Python notebook computes a complex qubit–qutrit reduction, compares three kinds of correlation, tests a noisy singlet, and truncates a two-mode squeezed state. Each result can also be checked by hand.

Helpful background. Product-basis ordering explains array indices; partial trace and Schmidt decomposition explain the calculations. Basic Python arrays are enough to follow the code.

Composite states as arrays

The file uses # %% cell markers, so a compatible notebook editor can run one cell at a time. It also runs as an ordinary script with Python 3 and NumPy:

python composite-state-calculations.py

The examples were executed with NumPy 2.1.2. They use dense arrays for small systems; a joint density matrix requires $(d_A d_B)^2$ entries, so the method is not a scalable many-body simulation.

The basis index is $a d_B+b$, with subsystem $A$ first. For two qubits this means $00,01,10,11$. Reshaping a density matrix exposes four indices in the order $(a,b,a',b')$:

r = rho.reshape(dA, dB, dA, dB)
rho_A = np.trace(r, axis1=1, axis2=3)
rho_B = np.trace(r, axis1=0, axis2=2)

Tracing over $B$ means setting its ket and bra indices equal and summing them. It does not mean keeping one block of the matrix. NumPy’s trace documentation explains the axis convention.

The downloadable functions check matrix shape, finite entries, Hermiticity, trace one, and positivity before computing entropies. A tolerance of $10^{-10}$ accommodates roundoff in these small examples; it is not an experimental uncertainty. Values outside that tolerance cause an error. Only accepted roundoff asymmetry or tiny negative eigenvalues are removed, rather than silently turning an invalid matrix into a state.

A complex state tests conjugation

Consider

$$|\psi\rangle=\frac{|00\rangle+i|11\rangle+|12\rangle}{\sqrt3}, \qquad C=\frac1{\sqrt3}\begin{pmatrix}1&0&0\\0&i&1\end{pmatrix}.$$

With $C_{ab}$ defined as the amplitude of $|a\rangle|b\rangle$, the two contractions are

rho_A = C @ C.conj().T
rho_B = C.T @ C.conj()

The second formula is easy to mistype. The correct result is

$$\rho_B=\frac13\begin{pmatrix}1&0&0\\0&1&i\\0&-i&1\end{pmatrix}.$$

Its nonzero eigenvalues are $1/3,2/3$, so either subsystem has entropy $0.918295834\ldots$ bits. But eigenvalues alone cannot detect every conjugation error: $C^\dagger C$ has the same eigenvalues and opposite imaginary entries. A local observable distinguishes them:

$$Y_{12}=-i|1\rangle\langle2|+i|2\rangle\langle1|, \qquad \langle Y_{12}\rangle=-\frac23.$$

The notebook also reconstructs the state from the singular-value decomposition:

U, s, Vh = np.linalg.svd(C, full_matrices=False)
restored = sum(s[k] * np.kron(U[:, k], Vh[k, :])
               for k in range(len(s)))

Here the components of the $B$ Schmidt ket are a row of Vh, used without further conjugation. Check the reconstructed vector itself, not just the singular values.

Three states with the same local statistics

The next cell compares

$$\rho_{\mathrm{Bell}}=|\Phi^+\rangle\langle\Phi^+|, \qquad \rho_{\mathrm{cc}}=\frac{|00\rangle\langle00|+|11\rangle\langle11|}{2}, \qquad \rho_{\mathrm{ind}}=I_4/4.$$

All three reduce to $I_2/2$ on either qubit. Their joint entropies are respectively $0,1,2$ bits, while their mutual informations are $2,1,0$ bits. These outputs distinguish entanglement, classical correlation, and independence despite identical single-qubit measurements.

For these finite matrices the entropy routine uses the Hermitian eigenvalues and computes $-\sum_{p_k>0}p_k\log_2p_k$. The zero eigenvalues contribute zero. The local entropy of a mixed joint state is not, by itself, its entanglement.

Partial transpose of a noisy singlet

For $0\le v\le1$, define

$$\rho_v=v|\Psi^-\rangle\langle\Psi^-|+(1-v)I_4/4.$$

The notebook transposes only the two $B$ indices:

rho_PT = r.transpose(0, 3, 2, 1).reshape(dA*dB, dA*dB)

At $v=0.6$, the eigenvalues of $\rho_v^{T_B}$ are $-0.2,0.4,0.4,0.4$. The negative eigenvalue proves entanglement, and its absolute value gives negativity $0.2$. Try $v=0.2$: all eigenvalues become nonnegative. For this two-qubit system PPT is also sufficient for separability; that last conclusion is dimension dependent. See negativity and PPT before applying it to other systems.

Partial transpose need not produce a positive matrix. Consequently the notebook uses a Hermitian eigensolver directly on that result, rather than passing it to the entropy function that expects a physical density matrix.

Occupation cutoffs and discarded probability

A two-mode squeezed state has Schmidt probabilities $p_n=(1-z)z^n$, where $z=\tanh^2r<1$. Keeping $n=0,\ldots,N$ discards exactly

$$\delta=\sum_{n=N+1}^\infty p_n=z^{N+1}.$$

For $z=1/2$ and $N=9$, the discarded probability is $1/1024$. The normalized truncated pure state has squared overlap $1-\delta=1023/1024$ with the exact state. This is a controlled approximation to the state vector; an unbounded observable can require additional tail estimates.

A finite occupation basis also changes operator identities. If $a_N$ is the annihilation operator projected onto $|0\rangle,\ldots,|N\rangle$, then

$$[a_N,a_N^\dagger]=I-(N+1)|N\rangle\langle N|.$$

The top-state correction is a truncation effect, not a floating-point mistake. Increase the cutoff and check the quantity you actually want, rather than treating a small normalization error as proof of convergence. Two-mode entanglement explains the exact squeezed-state distribution.

Exercises

1. Replace $i|11\rangle$ by $-i|11\rangle$ in the complex example. What changes, and what stays fixed?

Solution

The off-diagonal entries of $\rho_B$ are complex conjugated and $\langle Y_{12}\rangle$ becomes $+2/3$. The reduced eigenvalues $1/3,2/3$, purity $5/9$, and entropy stay fixed. This pair of states illustrates why matching entropies does not establish matching measurement statistics.

2. Derive the partial-transpose eigenvalues of $\rho_v$, and compare the output at $v=0.2$ with the separability threshold.

Solution

The singlet’s partial transpose has eigenvalues $-1/2,1/2,1/2,1/2$. Adding the identity contribution gives $(1-3v)/4$ once and $(1+v)/4$ three times. At $v=0.2$ these are $0.1,0.3,0.3,0.3$. The state is PPT, hence separable for two qubits. The exact threshold is $v=1/3$; negativity above it is $(3v-1)/4$.

3. For $z=1/2$, choose the smallest cutoff $N$ with discarded probability below $10^{-6}$.

Solution

We need $2^{-(N+1)}<10^{-6}$. Since $2^{-19}\approx1.907\times10^{-6}$ and $2^{-20}\approx9.537\times10^{-7}$, the smallest cutoff is $N=19$. There are then twenty retained number states per mode.

4. Construct $a_N$ for $N=2$ and verify the boundary correction in its commutator.

Solution

In the basis $|0\rangle,|1\rangle,|2\rangle$,

$$a_2=\begin{pmatrix}0&1&0\\0&0&\sqrt2\\0&0&0\end{pmatrix}.$$

The products are $a_2a_2^\dagger=\operatorname{diag}(1,2,0)$ and $a_2^\dagger a_2=\operatorname{diag}(0,1,2)$. Their difference is $\operatorname{diag}(1,1,-2)=I-3|2\rangle\langle2|$. Every finite-dimensional commutator has trace zero, so it cannot equal the identity.

References

• G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins University Press, 2013.
• M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information, 10th anniversary ed., Cambridge University Press, 2010.
• NumPy Developers, NumPy reference: trace, singular-value decomposition, and Hermitian eigenvalues, accessed 13 October 2026.
• J. Watrous, The Theory of Quantum Information, Cambridge University Press, 2018; author’s text.