Introduction

Ground states of the atomic nuclei exhibit nucleonic pairing correlations. This manifests itself in odd-even mass staggering, properties of low-lying excited states, moments of inertia, etc., to name but a few examples [1,2]. To successfully describe these phenomena, nucleonic pairing is usually introduced within mean-field models and handled through the Bogoliubov-Valatin transformation [2], whereupon the particle degrees of freedom are replaced by quasiparticles. This effectively incorporates the pairing correlations but, as a consequence, leads to particle-number-mixed wave functions. Situation is the same also for other symmetries broken on the mean-field level, e.g., by allowing nucleus to deform, quadrupole-type correlations are effectively incorporated in the mead-field picture, at the expense of breaking the rotational invariance. Nevertheless, true ground states conserve all symmetries of the underlying Hamiltonian, including the particle number.

To link the spontaneous breaking of symmetries to symmetry-conserving states, various symmetry-restoration schemes have been utilized. In principle, in self-consistent approaches solved within iterative methods, broken symmetries should be restored during every step towards solution. This is the well known variation-after-projection (VAP) [3,4,5,6] method. The drawback of VAP is computationally expensive integration over the gauge angles of symmetries to be restored, applied during every iteration step. Therefore, usual practice is to revert to computationally less intensive projection-after-variation (PAV) [2,4,7] scheme, where symmetries are restored at the end, from the converged self-consistent symmetry-broken mean-field solution.

When applying VAP method to superfluid nuclear density functional theory (DFT), some of the energy density functionals (EDFs) seem to be ill-suited for the task. In particularly, with widely used Skyrme-like EDFs, several pathologies exists. With particle number restoration, poles and non-analytic behavior preclude obtaining a unique solution [4,8,9,10]. These difficulties are traced to the non-integer powers of density in the employed EDFs [8,11]. The same also holds for the restoration of angular momentum [12].

To circumvent prohibitive computational cost of the VAP method, an approximate method is called for. The central rationale in the Lipkin method, which fulfill this goal, is to replace the original Hamiltonian by an auxiliary Routhian, making the symmetry-projected states degenerate in energy [13]. This allows us to approximately evaluate ground-state properties of the corresponding symmetry-restored system without actually performing any projection [13,14]. In particularly for the particle number restoration, a power series expansion as a function of particle number fluctuation was suggested. Based on this idea, Nogami [15,16] introduced a prescription to calculate coefficients of the power expansion at second order, which is called the Lipkin-Nogami (LN) method, and which has been widely used in nuclear DFT calculations [7,17].

Quantitative effect of the particle-number-restoration largely depends on whether the pairing correlations are strong (mid-shell) or weak (near closed shell). As pointed out in Refs. [18,19,20,21], the parabolic approximation, which corresponds to a sum up to the second order of the Lipkin or Kamlah [22] approximation, may fail at the limit of weak pairing. In this work we extend the Lipkin method beyond second order used so far, so as to make the first tests of its convergence and accuracy.

The central issue of the Lipkin method is to search for a suitable set of the Lipkin power-expansion parameters [14]. For the LN method, the second-order parameter is calculated through the diagonal matrix elements [15,16,23,24], which requires us to calculate the linear response of the mean field to the particle-number projection [23,24]. However, the response term, which has large influence on potential energy surface (PES) [23], is cumbersome to evaluate. Therefore, usually in calculations involving the LN method, an approximate prescription of seniority pairing is used to obtain the effective pairing strength for the second-order term [25,5]. In this work, we propose a different way to derive these expansion parameters, namely, starting from the non-diagonal energy kernels.

This paper is organized as follows: In Sec. 2, we cover our theoretical framework of the Lipkin method for particle number projection. In Sec. 3, we present numerical results, and in Sec. 4, we give the summary and outlook. Appendix A contains explicit expressions of the Lipkin method applied in this work and Appendix B provides an illustration of the method within the exactly solvable two-level pairing model.