As the Hamiltonian eigenstates are pure states, we can evaluate their entanglement by bypassing the construction of the density operator and use the Schmidt decomposition instead. Specifically, for a given two-body wave function \( \vert \Psi\rangle \) expressed in terms of the Hartree product states, we can write
$$ \begin{align*} \vert\Psi\rangle &= \sum_{k = 0}^{N^L} \sum_{l = 0}^{N^R} C_{kl}\vert\phi^L_k \phi^R_l\rangle=\sum_{p = 0}^{\tilde{N}}\sigma_{p}\vert\psi^L_p \psi^R_p\rangle, \end{align*} $$where \( C_{kl} = \sum_{p = 0}^{\tilde{N}} U_{kp}\sigma_{p} V^{*}_{lp} \) is the singular value decomposition of the two-body coefficients. We have
$$ \begin{gather*} \vert\psi^L_p\rangle\equiv \sum_{k = 0}^{N^L} U_{kp} \vert\phi^L_k\rangle, \qquad \vert\psi^R_p\rangle\equiv \sum_{l = 0}^{N^R} V^{*}_{lp} \vert\phi^R_l\rangle, \end{gather*} $$are the Schmidt states, \( \tilde{N} \) is either \( N^L \) or \( N^R \) depending on the definition of the singular value decomposition, and \( \sigma_p \) are the singular values with \( \sigma_p^2 \) representing the occupation of the pair \( \vert\psi^L_p \psi^R_p\rangle \).