Bounds on 2-body reduced density matrix for a singlet state consisting of N fermions with spin S.

Question

Problem setup: consider a many-body system containing $N$ fermions of half-integer spin $S$ and focus on the states with total spin $S_{tot}=0$. For example, when $N=2$, $S=1/2$, we will get the usual spin singlet state. For a pure $S_{tot}=0$ state \(|\psi\rangle\) (or more generally a mixed state consisting only of $S_{tot}=0$ pure states), calculate the 2-body reduced density matrix \(\rho^2\). The condition $S_{tot}=0$ dictates that the only non-zero elements are

\[\rho^2_{J} \equiv \rho^2_{JM,JM} \equiv\tfrac{1}{2} \langle \psi | \xi^\dagger_{JM}\xi_{JM}| \psi\rangle, \quad \xi^\dagger_{JM} = C^{SSJ}_{m_1m_2M} c^\dagger_{m_1} c^\dagger_{m_2}\]

where \(\xi^\dagger_{JM}\) is the creation operator for a pair of fermions with an even total spin \(J=2k, k\in[0,(2S-1)/2]\) and z-component \(M \in [-J,J]\),  \(c^\dagger_{m}\)  is the creation operator for a fermion of z-component \(m \in [-S,S]\), and $C$ is the CG coefficient.

I'd like to find some \(\mathcal O(1)\) bounds on these \(\rho^2_J\)'s (suggested by numerical experiments). For example for \(J=2S-1\), \(\rho^2_{2S-1} = \langle \psi | c^\dagger_S c^\dagger_{S-1} c_{S-1} c_S | \psi \rangle \leq 1\).  I sense this problem might have something to do with the N-representability of a 2-body density matrix, i.e., whether it can be viewed as reduced from a full N-body density matrix, plus the $S_{tot}=0$ constraint. In that context, it might be possible to characterize the high-dimensional manifold containing all possible combinations of the  \(\rho^2_J\)'s. Here are two constraints that might be helpful:

\[\begin{align} &\sum_J (2J+1)\rho^2_J = N(N-1)/2\\ &\sum_J (2J+1)J(J+1)\rho^2_J = N(N-2)S(S+1) \end{align}\]

The first one is just saying that the trace of the two-body reduced density matrix is equal to the number of pairs. The second one comes from the $S_{tot}=0$ constraint.

I'd appreciate any suggestions you might have.