Lipkin method of particle-number restoration to higher orders

Background: On the mean-field level, pairing correlations are incorporated through the Bogoliubov-Valatin transformation, whereupon the particle degrees of freedom are replaced by quasiparticles. This approach leads to a spontaneous breaking of the particle-number symmetry and mixing of states with different particle numbers. In order to restore the particle number, various methods have been employed, which are based on projection approaches before or after variation. Approximate variation-after-projection (VAP) schemes, utilizing the Lipkin method, have mostly been used within the Lipkin-Nogami prescription. Purpose: Without recurring to the Lipkin-Nogami prescription, and using instead states rotated in the gauge space, we derive the Lipkin method of particle-number restoration up to sixth order and we test the convergence and accuracy of the obtained expansion. Methods: We perform self-consistent calculations using the higher-order Lipkin method to restore the particle-number symmetry in the framework of superfluid nuclear energy-density functional theory. We also apply the Lipkin method to a schematic exactly solvable two-level pairing model. Results: Calculations performed in open-shell tin and lead isotopes show that the Lipkin method converges at fourth order and satisfactorily reproduces the VAP ground-state energies and energy kernels. Near closed shells, the higher-order Lipkin method cannot be applied because of a non-analytic kink in the ground-state energies in function of the particle number. Conclusions: In open-shell nuclei, the higher-order Lipkin method provides a good approximation to the exact VAP energies. The method is computationally inexpensive, making it particularly suitable, for example, for future optimizations of the nuclear energy-density functionals and simultaneous restoration of different symmetries.


I. 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 meanfield 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 symmetryrestoration 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 know variationafter-projection (VAP) [3][4][5] 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-aftervariation (PAV) [2,6] 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 [7][8][9]. These difficulties are traced to the non-integer powers of density in the employed EDFs [7,10]. The same also holds for the restoration of angular momentum [11].
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 [12]. This allows us to approximately evaluate ground-state properties of the corresponding symmetry-restored system without actually performing any projection [12,13]. In particularly for the particle number restoration, a power series expansion as a func-tion of particle number fluctuation was suggested. Based on this idea, Nogami [14,15] 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 [6,16]. 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. [17][18][19][20], the parabolic approximation, which corresponds to a sum up to the second order of the Lipkin or Kamlah [21] 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 [13]. For the LN method, the second-order parameter is calculated through the diagonal matrix elements [14,15,22,23], which requires us to calculate the linear response of the mean field to the particle-number projection [22,23]. However, the response term, which has large influence on potential energy surface (PES) [22], 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 [4,24]. 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. II, we cover our theoretical framework of the Lipkin method for particle number projection. In Sec. III, we present numerical results, and in Sec. IV, 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.

II. LIPKIN METHOD
To start, we first recall some of the standard definitions available in the literature, which are required in the present work. Within the Hartree-Fock-Bogoliubov (HFB) framework, the wave function rotated in the gauge space is defined as [2,6] where φ is the gauge angle,N is the particle-number operator, and N 0 = Φ|N |Φ is the average particle number. In what follows, for a sake of clarity, we present expressions for a system composed only of one kind of particles. Generalizations to two types of particles, that is, to protons and neutrons, is straightforward and is discussed briefly later. Similarly as in Ref. [13], the overlap and energy kernels are defined as and kernels of (N − N 0 ) m as The kernels of Eq. (4) can be calculated as derivatives of the overlap kernel with respect to the gauge angle Explicit expression for these kernels are presented in Appendix A. In Eqs. (3) and (4), kernels are defined in terms of matrix elements. However, within the EDF methods they have to be understood as standard functions of transition density matrices, see, e.g., discussion in Ref. [7]. As demonstrated by Lipkin [12], variation after the particle-number projection (VAPNP) can be performed using an auxiliary Routhian, where Lipkin operatorK, which is a function of the shifted particle-number operatorN − N 0 , is chosen so as to "flatten" the N -dependence of average Routhians calculated for the particle-number projected states [12,13]. Had these projected Routhians been exactly Nindependent (perfectly flat), the exact projected energy E N0 could have been obtained by minimizing the average value of the Routhian for the unprojected state |Φ , that is, Otherwise, the Lipkin method constitutes a suitable approximate VAPNP method, and its accuracy depends on the quality of the choice made for the Lipkin operatorK.
As suggested by Lipkin [12], the simplest and manageable ansatz for the Lipkin operatorK has the form of a power expansion, where k 1 ≡ λ is the Fermi energy, which is used to fix the average particle number. Furthermore, higher-order Lipkin parameters k m are used to best describe the particlenumber dependence of the average energies of projected states.
Up to now, the LN method was frequently used to estimate values of k 2 (traditionally denoted by λ 2 ). However, this method relies on calculating the average values of Φ|ĤN m |Φ and Φ|N m |Φ , and, thus, at higher orders (m > 2) evaluation of these terms becomes cumbersome and impractical.
The essence of the original Lipkin method is different, namely, it relies on deriving expressions for k m that "flatten" the φ-dependence of the reduced Routhian kernel h ′ (φ), that is, where Up to any order, this is a perfectly manageable setup, because for an arbitrary value of the gauge angle, the generalized Wick theorem [2] allows for a simple determination of the energy and overlap kernels H(φ) and N m (φ). Explicit expressions for k m are presented in Appendix A. For a perfectly flat (φ-independent) reduced Routhian kernel h ′ (φ) ≡ C, we then have the exact average value of the Routhian evaluated for the state projected on particle number N 0 , Since for the state projected on N 0 , the average value of the Lipkin operator (8) is, by definition, equal to zero, we also have that and thus the minimization of the average Routhian (7) is equivalent to the exact VAPNP. Again, any imperfection in the φ-independence of h ′ (φ) amounts to a certain approximation of the exact VAPNP. However, since it is now relatively easy to go to higher orders in the power expansion of Eq. (8), we can systematically test the convergence of this expansion. The largest contributions to integrals in Eq. (11) come from the vicinity of the origin due to the largest weight [16]. Therefore, we can evaluate Lipkin parameters k m using the gauge-rotated intrinsic states near the origin. This also avoids the singularities caused by vanishing overlaps [7]. As an example, at second order one obtains the Lipkin parameter, where φ 2 is a pre-selected small value of the gauge angle, and the flattened energy reads Had the expansion up to second order been exact, values of k 2 and E N0 obtained from Eqs. (13) and (14) would have been independent of φ 2 . Thus, their eventual dependence on φ 2 indicates the necessity of going beyond second order.
Similarly, at order M , we evaluate Lipkin parameters k m , m = 1, . . . , M , using a set of M small gauge angles φ i , i = 1, . . . , M . In practice, in this work, we use equally spaced values of φ i = iφ 1 , and at each order we check the eventual dependence of results on the maximum gauge angle used, φ M .
We note here that the minimization of the average Routhian (7) with respect to the HFB state |Φ can be performed by solving the standard HFB equation with additional higher-order terms added, see Appendix A. We also note that Lipkin parameters k m must be determined in each HFB iteration (for each current state |Φ ), in such a way that at the end of the HFB convergence they correspond to the final self-consistent solution, and thus parametrically depend on it. However, this dependence does not give rise to any additional terms in the HFB equation, because the derivation of the Lipkin method is based on treating them as constants, cf. discussion of the LN and Kamlah methods in Ref. [18]. An exactly solvable two-level pairing model offers an ideal environment to test qualitative properties of the Lipkin VAPNP method. The results presented in Appendix B show that in such a schematic model, the higher-order Lipkin VAPNP method is able to reproduce correctly the exact VAPNP ground-state energies, both in weak and strong pairing regimes, everywhere apart from the immediate vicinity of the closed shell. This gives us confidence in applications of this method in more involved cases of actual nuclei, which is discussed in the next section.

III. RESULTS AND DISCUSSION
We have implemented the Lipkin VAPNP method, presented in Sec. II, in computer code HFODD (v2.68c) [25,26]. This code solves the HFB equations in a threedimensional Cartesian harmonic oscillator basis. Within this implementation, we tested the Lipkin VAPNP method using the Skyrme SIII parameterization [27] in the particle-hole channel and volume zero-range pairing interaction in the particle-particle channel. SIII parameterization was selected for this study due to the fact that it contains only integer powers of densities. Used neutron pairing strength, V 0 = −155.45 MeV fm 3 , was adjusted within the LN method to reproduce the empirical neutron pairing gap of ∆ n = 1.245 MeV of 120 Sn.
In principle, at each given order of the Lipkin VAPNP method, this adjustment should be repeated. However, for the sake of meaningful comparison of the results obtained at different orders, we use the same pairing strength throughout all calculations.
For protons, pairing strength was set to zero, that is, the proton subsystem is described by unpaired states. This setup allows us to test the Lipkin VAPNP method  in the neutron paired subsystem, resulting in a clearer interpretation of the obtained results and allowing for a better evaluation of the efficiency of the Lipkin VAPNP method. Because of the used zero-range pairing interaction, we adopted the commonly used equivalent-spectrum cutoff of 60 MeV, applied in the quasiparticle configuration space. All calculations were performed in the spherical basis of 14 major harmonic-oscillator shells.
To begin with, we first study convergence of the Lipkin VAPNP method when terms up to sixth-order in expansion (8) are incorporated. At present, we limit our analysis to the even powers only, that is, we take into account terms with m = 2, 4, and 6. This corresponds to a symmetric approximation around the central value of the particle number N 0 . In Figs. 1 and 2, we show, respectively for 120 Sn and 100 Sn, dependence of the Lipkin parameters on the maximum gauge angle φ M (see the previous section). The figures also show the total Lipkin  VAPNP energy E N0 and Lipkin correction energy E corr cf. Eqs. (7) and (9).
At second order, the obtained results show a clear dependence on φ M , indicating insufficient expansion. On the other hand, at fourth and sixth orders, total energy E N0 and correction energy E corr are already rather insensitive to φ M . Thus, we can conclude that at sixth order, the expansion is well converged, and at least fourth order is required for sufficiently precise results. We also note that for the magic nucleus 100 Sn, the convergence is slightly slower, and the values of Lipkin parameters are significantly higher than those for 120 Sn.
In what follows, we have used the same maximum gauge angle of φ M = 2π 51 ≃ 0.123 in all expansions, regardless of the expansion order. In Fig. 3, convergence of the reduced kernels of Lipkin operator (8) in 100 Sn and 120 Sn is shown. Kernel values at φ = 0 were subtracted, in order to illustrate how well the reduced Routhian kernels (9) stay constant, that is, independent of the gauge angle φ. Again, we clearly see that the second-order expansion is insufficient, whereas fourth and sixth orders already give satisfactory description of the energy ker- nels. Figure 4 shows the same kernels as those plotted in Fig. 3 for the whole range of gauge angle, up to φ = 2π. We see that in 100 Sn the energy kernel near φ = π/2 is poorly described by the Lipkin expansion. This is directly related to the kink in the particle-number dependence of the projected energies, which appears at the magic shell closure, and which cannot be properly described by a polynomial expansion [4,18]. In this kind of case, a good quality of the Lipkin expansion obtained at small gauge angles is not sufficient enough to guarantee a good convergence at larger gauge angles. The situation is entirely different in the open-shell nucleus 120 Sn, where particle-number dependence of the projected energies is given by a smooth function, which can very well be approximated by a polynomial expansion. Here, for all gauge angles, we obtain a perfectly converging Lipkin expansion of the exact energy kernel, even in the vicinity of the pole related to the nearly half-filled 3s 1/2 orbital [7].
In Figs. 5 and 6, we show results of the Lipkin VAPNP method for tin and lead isotopes, respectively. For com- parison, figures also show results obtained using the LN method, similarly as in Ref. [28], and PLN method, as in Ref. [4], where the exact PNP energy is obtained via projection from HFB+LN self-consistent solution.
Away from the closed shells, at fourth and sixth orders, results of the Lipkin VAPNP method are very similar to those of PLN results. As pointed out by the VAPNP calculations in Ref. [4], for open shell nuclei, PLN results are very close to exact VAPNP results. Again, fourth and sixth orders give similar results, signaling the convergence of the Lipkin expansion. We can thus conclude that the fourth-order Lipkin VAPNP method is a good approximation of the exact VAPNP method. Near shell closure, differences between various orders of the Lipkin VAPNP method are large, indicating a non-convergent power series of the Lipkin operator. Once again, this is related to the kink in the particle-number dependence on the projected energies [4,18].
In Fig. 7, we show results obtained by projecting good particle numbers from the states obtained either by the Lipkin VAPNP or LN methods. It is very gratifying to see that irrespective on whether one uses the Lipkin VAPNP or LN methods, the projected energies, shown in the two top panels, are very similar. This fact means that all approximate methods analyzed in this study lead to similar pair condensates, whereas they differ in the determination of corrective mean-field energies. The main advantage of using the Lipkin VAPNP method is in the fact that the PNP calculation does not have to be performed at all. Then, as shown in lower panels of Fig. 7, the obtained energies very well approximate the PNP energies. This is particularly true near the middle of the shell, where the influence of closed-shell kinks in the projected energies is weaker.
We also see that at closed shells, the second-order Lipkin VAPNP method, similarly as the LN method -see Figs. 5 and 6 -gives results that are very different from those obtained by the PNP. On the contrary, the fourthand sixth-order Lipkin VAPNP methods gives results almost identical to the PNP. Finally, at fourth and sixth orders, the non-analytic behavior of the PNP energies at closed shells causes the largest discrepancies for 2 or 4 particles away from the closed shell.
Our current implementation of the Lipkin method in the computer code HFODD allows us to treat pairing correlations simultaneously for neutrons and protons. However, this has been implemented such that the Lipkin operator is simply a sum of neutron and proton contributions of Eq. (8), with Lipkin parameters determined by independent gauge-angle rotations for neutrons and protons. Although this implementation works perfectly well, we have realized that such a method is insufficient in some cases. This is illustrated in Fig. 8(a), which shows the reduced energy kernel of 124 Xe calculated in two dimensions, as a function of the neutron φ n and proton φ p gauge angles. We clearly see that energy kernel is tilted with respect to main axes of neutron and proton gauge angle. Evidently, the Lipkin operator, here being a sum of neutron and proton contributions separately, leads to a non-tilted energy kernel, as shown in Fig. 8(b). Therefore, to fully reproduce the true energy kernel, one has to use the Lipkin operator that contains cross terms, which depend on products of neutron and proton particle numbers. Within proton-neutron pairing scheme, combined with PLN, these kind of cross-terms are required [29,30]. Implementations of such more complicated forms of the Lipkin operator will be subject of future study.

IV. SUMMARY
In the framework of the nuclear energy-densityfunctional theory, we derived the Lipkin method of approximate particle-number-symmetry restoration up to sixth order. The Lipkin parameters were determined from non-diagonal energy kernels, resulting in a more manageable approach as compared to the traditional Lipkin-Nogami approach.
Convergence of the Lipkin VAPNP method was tested by investigating gauge dependence of expansion parameters. Taking 120 Sn as an example, the Lipkin expansion up to the second order was found to have explicit gauge angle dependence. Inclusion of fourth order terms subsequently diminished dependence on the gauge angle significantly. With the inclusion of sixth-order terms of the expansion, overall change is minimal, indicating a converging series. Accuracy of the Lipkin VAPNP method was tested by comparing the reduced energy kernel and Lipkin operator approximated by a power series. It was found that the chosen Lipkin operator describes well the small gauge-angle rotation of the intrinsic wave function. The results obtained for 100 Sn and 120 Sn showed that the second-order Lipkin expansion is typically not sufficiently converged. Within fourth order, the series already mimics the reduced energy kernels rather well. With inclusion of sixth order term, the results stay practically the same, indicating again a well converged series.
In chains of tin and lead isotopes, we compared the Lipkin VAPNP method to LN and PLN methods. As pointed in Ref. [4], for mid-shell nuclei, PLN is a very good approximation to the exact VAP method. Our results show that for the mid-shell nuclei, the Lipkin VAPNP method already at the second order gives rather well converged results. When advancing to higher orders, the results are improved. Near closed shells, because of the kink in the particle-number dependence on the projected energy, the Lipkin VAPNP method is unable to reproduce the exact projected energy. Also, near the closed-shell region, the pairing correlations have dynamic nature [2]. Within the Lipkin VAPNP method, these kind of features cannot be reproduced with a well converging series expansion. The main features of the results obtained for nuclei were corroborated within the exactly-solvable two-level pairing model.
When neutrons and protons were treated simultaneously within the Lipkin VAPNP method, we observed a necessity to include in the Lipkin operator cross terms, which depend simultaneously on neutron and proton number operators. For the case of 124 Xe, the contour lines of the reduced energy kernel, with respect to neutron's and proton's gauge angles, show tilted shapes. Without the cross terms, this kind on behavior cannot be reproduced. A study of the cross terms will be subject of future work.
The Lipkin VAPNP method presented in this work allows for a computationally inexpensive way to approximate the exact VAPNP energy of the ground state. The Lipkin method can be also applied to approximate the restoration of other symmetries, broken on the meanfield level, with a small or no extra computational cost. This is important, for example, in future adjustments of novel EDFs, tailored for beyond mean-field multireference studies [31]. Work towards restoring isospin and rotational symmetries within the framework of the Lipkin method is currently in progress.

Acknowledgments
This work was supported in part by the Academy of To calculate kernels N m (φ) from derivatives of the overlap kernel I(φ), in Eq. (5), we use the Onishi theorem [2]: where C is the Thouless matrix. Next, we use the identity which is valid for any matrix A, provided that dA dφ commutes with A, which is our case here. This gives, where ρ(φ) = e 2iφ CC + (1 + e 2iφ CC + ) −1 .
We then have the reduced kernels n m (φ) = N m (φ)/I(φ) given as To lighten the notation, we omitted the arguments of ρ(φ) and R n (φ) from here on. Derivatives of the density matrix ρ(φ) (A4) can be calculated as, Up to the sixth order, the average Routhian of Eq. (7) to be minimized reads Average values of powers of the particle-number operator are given in Eqs. (A5)-(A10) taken at φ = 0. Moreover, in all terms with m ≥ 2, one can set R 0 ≡ 0. This gives Hence, the corresponding mean-field Routhians (to be used in the HFB equations) reads [32] h Lipkin parameters k m for m = 1, . . . , M can be determined from Eq. (9) by requiring that it is fulfilled at gauge angle φ = φ 0 = 0 and also at all M other nonzero values of the gauge angle φ i . This gives where C is the flattened Routhian. Then, at sixth order, Lipkin parameters k m can be easily obtained by inverting the matrix built of coefficients n m (φ i ) as 1 n 1 (0) · · · n 6 (0) 1 n 1 (φ 1 ) · · · n 6 (φ 1 ) 1 n 1 (φ 2 ) · · · n 6 (φ 2 ) 1 n 1 (φ 3 ) · · · n 6 (φ 3 ) 1 n 1 (φ 4 ) · · · n 6 (φ 4 ) 1 n 1 (φ 5 ) · · · n 6 (φ 5 ) 1 n 1 (φ 6 ) · · · n 6 (φ 6 ) At lower orders, or when neglecting odd orders, a smaller number of the gauge angle points can be used. In fact, the value of C obtained from the first row of Eq. (A27) does not appear in the mean-field Routhian (A25) and can be ignored. Anyhow, at convergence it is calculated from Eq. (A24). Moreover, the value of k 1 obtained from the second row of Eq. (A27) can be ignored too, because, when iterating the HFB equation, it is anyhow used to properly adjust the average particle number, and, at convergence, one has Trρ = N 0 and thus it has no influence on the value of the right-hand side of Eq. (A24).
Appendix B: Lipkin method applied to the two-level pairing model In this appendix, we apply the Lipkin method to the standard two-level pairing model, which is characterized by two Ω-fold degenerate levels with the single-particle energy difference 2ǫ and pairing strength G. Below we closely follow the notations and definitions presented in Refs. [17,18], where the results obtained within the LN method have been studied. In Fig. 9, we show particle-number dependence of the   [17,18] as differences between the total and Hartree-Fock energies. Note that the results are exactly symmetric with respect to the mid shell, that is, those for particle numbers of N and 2Ω− N are exactly identical.
The ratio of R = 1, that is, perfect agreement, is for all particle numbers reached in the strong-pairing regime. For weak pairing, the largest discrepancies appear at mid shell, N = 20, and they gradually decrease towards smaller (or larger) particle numbers. This is related to the kink in particle-number dependence of ground-state energies [18], cf. Fig. 9, which disappears with increasing pairing correlations.
For N = 20, with increasing order of the Lipkin ex-  pansion, the agreement with exact results gradually increases, and the Lipkin VAPNP method, even at second order, is here visibly superior to the LN method. Note that at N = 20, the odd orders of expansion (third and fifth) do not bring any improvement -this is owing to the symmetry of the model with respect to the mid shell. For N = 18, the Lipkin expansion cannot reproduce the kink appearing at the adjacent particle number of N = 20 (see Fig. 9), and it does not seem to converge to the exact result, whereas the LN results are clearly superior. For smaller particle numbers, this pattern gradually changes, and for N ≤ 12 the Lipkin expansion does converge to the exact result and at orders higher than four becomes better than the LN method.
We stress here that in the realistic cases discussed in Sec. III, the pattern of comparison between the LN and Lipkin VAPNP methods pertains to moderately high pairing strengths, certainly beyond the pairing phase transition, which in the two-level model appears at x c = 1/(Ω − 1) ≃ 0.053.
Finally, in Fig. 11, we show dependence of results on the maximum gauge angle φ M used in the Lipkin VAPNP method, see Secs. II and III. We see that for all particle numbers, the second-order results do depend on φ M , indicating an insufficient order of expansion. For N ≤ 12 we see that with increasing order of expansion, the results become perfectly independent of φ M , which characterizes a converging expansion. On the other hand, closer to the mid shell, even at sixth order a visible dependence on φ M still remains.