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.