  We present an extension of the pair coupled cluster doubles (p-CCD) method to
quasiparticles and apply it to the attractive pairing Hamiltonian. Near the
transition point where number symmetry gets spontaneously broken, the proposed
BCS-based p-CCD method yields significantly better energies than existing
methods when compared to exact results obtained via solution of the Richardson
equations. The quasiparticle p-CCD method has a low computational cost of
$\mathcal{O}(N^3)$ as a function of system size. This together with the high
quality of results here demonstrated, points to considerable promise for the
accurate description of strongly correlated systems with more realistic pairing
interactions.
