Recently, ground state eigenvectors of the reduced Bardeen-Cooper-Schrieffer (BCS) Hamiltonian, Richardson-Gaudin (RG) states, have been employed as a wavefunction ansatz for strong correlation. This wavefunction physically represents a mean-field of pairs of electrons (geminals) with a constant pairing strength. To move beyond the mean-field, one must develop the wavefunction on the basis of all the RG states. This requires both practical expressions for transition density matrices and an idea of which states are most important in the expansion. In this contribution, we present expressions for the transition density matrix elements and calculate them numerically for half-filled picket-fence models (reduced BCS models with constant energy spacing). There are no Slater-Condon rules for RG states, though an analog of the aufbau principle proves to be useful in choosing which states are important.