Because the resulting matrix is sparse, we implement its numerical diagonalization employing the well known ARPACK solver which uses implicitly restarted Arnoldi methods. Convergence of the many-body solutions is guaranteed by using a large enough value for the dimension K of the single-particle basis; we used here K ∼ 50. Naturally, an important ingredient in this respect are the matrix elements of the two-body interaction,