A new approach for solving the Bethe ansatz (Gaudin-Richardson) equations of the standard pairing problem is established based on the Heine-Stieltjes correspondence. For k pairs of valence nucleons on n different single-particle levels, it is found that solutions of the Bethe ansatz equations can be obtained from one (k + 1) × (k + 1) and one (n − 1) × (k + 1) matrices, which are associated with the extended Heine-Stieltjes and Van Vleck polynomials, respectively. Since the coefficients in these polynomials are free from divergence with variations in contrast to the original Bethe ansatz equations, the approach thus provides with a new efficient and systematic way to solve the problem, which, by extension, can also be used to solve a large class of Gaudin-type quantum many-body problems and to establish a new efficient angular momentum projection method for multi-particle systems.PACS numbers: 21.60. Cs, 21.60.Fw, 03.65.Fd, 71.10.Li, 74.20.Fg, 02.60.Cb It is well known that the pairing force, similar to that in the Bardeen-Cooper-Schrieffer (BCS) theory of superconductors [1], as one of main residual interactions introduced to the nuclear shell model, is key to manifest ground state properties and low energy spectroscopy of nuclei, such as binding energies, odd-even effects, singleparticle occupancies, excitation spectra, electromagnetic transition rates, beta-decay probabilities, transfer reaction amplitudes, low-lying collective modes, level densities, and moments of inertia, and so on [2]- [4]. Unlike electrons in solids, the drawbacks of the application of the BCS theory and its extensions to nuclei are noticeable due to the fact that the number of valence nucleons under the influence of the pairing force is too few to be treated by such particle-number nonconservation (quasiparticle) approximations
