  Eigenvectors of the reduced Bardeen-Cooper-Schrieffer Hamiltonian have
recently been employed as a variational wavefunction ansatz in quantum
chemistry. This wavefunction is a mean-field of pairs of electrons (geminals).
In this contribution we report optimal expressions for their reduced density
matrices in both the original physical basis and the basis of the
Richardson-Gaudin pairs. Physical basis expressions were originally reported by
Gorohovsky and Bettelheim. In each case, the expressions scale like
$\mathcal{O}(N^4)$, with the most expensive step the solution of linear
equations. Analytic gradients are also reported in the physical basis. These
expressions are an important step towards practical mean-field methods to treat
strongly-correlated electrons.
