  Background: The nuclear many-body system is a strongly correlated quantum
system, posing serious challenges for perturbative approaches starting from
uncorrelated reference states. The last decade has witnessed considerable
progress in the accurate treatment of pairing correlations, one of the major
components in medium-sized nuclei, reaching accuracies below the 1% level of
the correlation energy.
  Purpose: Development of a quantum many-body method for pairing correlations
that is (a) competitive in the 1% error range, and (b) can be systematically
improved with a fast (exponential) convergence rate.
  Method: The present paper capitalizes upon ideas from Richardson-Gaudin
integrability. The proposed method is a two-step approach. The first step
consists of the optimization of a Richardson-Gaudin ground state as variational
trial state. At the second step, the complete set of excited states on top of
this Richardson-Gaudin ground state is used as an optimal basis for a
Configuration Interaction method in an increasingly large effective Hilbert
space.
  Results: The performance of the variational Richardson-Gaudin (varRG) and
Richardson-Gaudin Configuration Interaction (RGCI) method is benchmarked
against exact results using an effective $G$-matrix interaction for the Sn
region. The varRG already reaches accuracies around the 1% level of the
correlation energies, and the RGCI step sees an additional improvement scaling
exponentially with the size of the effective Hilbert space.
  Conclusions: The Richardson-Gaudin models of integrability provide an
optimized complete basis set for pairing correlations.
