Lipkin method of particle-number restoration to higher orders

This is an electronic reprint of the original article. This reprint may differ from the original in pagination and typographic detail. Author(s): Title: Year: Version: Please cite the original version: All material supplied via JYX is protected by copyright and other intellectual property rights, and duplication or sale of all or part of any of the repository collections is not permitted, except that material may be duplicated by you for your research use or educational purposes in electronic or print form. You must obtain permission for any other use. Electronic or print copies may not be offered, whether for sale or otherwise to anyone who is not an authorised user. Lipkin method of particle-number restoration to higher orders Wang, Xiaobao; Dobaczewski, Jacek; Kortelainen, Markus; Yu, Lingfei; Stoitsov, M. V. Wang, X., Dobaczewski, J., Kortelainen, M., Yu, L., & Stoitsov, M. V. (2014). Lipkin method of particle-number restoration to higher orders. Physical Review C, 90(1), Article 01431. https://doi.org/10.1103/PhysRevC.90.014312 2014 PHYSICAL REVIEW C 90, 014312 (2014) Lipkin method of particle-number restoration to higher orders X. B. Wang (),1J. Dobaczewski,1,2,*M. Kortelainen,1,2L. F. Yu (),1and M. V. Stoitsov3,4 1Department of Physics, P.O. Box 35 (YFL), University of Jyv¨ askyl¨ a, FI-40014 Jyv¨ askyl¨ a, Finland 2Helsinki Institute of Physics, P.O. Box 64, University of Helsinki, FI-00014 Helsinki, Finland 3Department of Physics & Astronomy, University of Tennessee, Knoxville, Tennessee 37996, USA 4Physics Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA (Received 25 March 2014; revised manuscript received 16 June 2014; published 17 July 2014) Background: On the mean-field level, pairing correlations are incorporated through the Bogoliubov-Valatin transformation, whereby 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 employing 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 particlenumber 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 nonanalytic kink in the ground-state energies as a 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. DOI: 10.1103/PhysRevC.90.014312 PACS number(s): 21.10.−k,21.60.Jz,74.20.−z,71.15.Mb I. INTRODUCTION Ground states of atomic nuclei exhibit nucleonic pairing correlations. These manifest 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], whereby 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. The situation is the same also for other symmetries broken on the mean-field level; e.g., by allowing the nucleus to deform, quadrupole-type correlations are effectively incorporated in the mean-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 the *[email protected] solution. This is the well-known variation-after-projection (VAP) [3–6] method. The drawback of VAP is the computationally expensive integration over the gauge angles of symmetries to be restored, applied during every iteration step. Therefore, the usual practice is to revert to the 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 the 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 the widely used Skyrme-like EDFs, several pathologies exists. With particle-number restoration, poles and nonanalytic behavior preclude obtaining a unique solution [4,8–10]. These difficulties are traced to the noninteger powers of density in the employed EDFs [8,11]. The same also holds for the restoration of angular momentum [12]. To circumvent the prohibitive computational cost of the VAP method, an approximate method is called for. The central rationale in the Lipkin method, which fulfills 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 particular, for the particle-number restoration, a power series expansion as a function of the particle-number fluctuation was 0556-2813/2014/90(1)/014312(11) 014312-1 ©2014 American Physical Society WANG, DOBACZEWSKI, KORTELAINEN, YU, AND STOITSOV PHYSICAL REVIEW C 90, 014312 (2014) 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]. The quantitative effect of the particle-number-restoration largely depends on whether the pairing correlations are strong (midshell) or weak (near closed shell). As pointed out in Refs. [18–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 the 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 (PNP) [23,24]. However, the response term, which has a large influence on the potential energy surface [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 [5,25]. In this work, we propose a different way to derive these expansion parameters, namely, starting from the nondiagonal energy kernels. This paper is organized as follows: In Sec. II, we cover our theoretical framework of the Lipkin method for PNP. In Sec. III, we present numerical results, and in Sec. IV,we give the summary and outlook. Appendix Acontains explicit expressions of the Lipkin method applied in this work and Appendix Bprovides 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,7] |(φ)=exp(iφ(ˆ N−N0))|,(1) where φis the gauge angle, ˆ Nis the particle-number operator, and N0=|ˆ N|is the average particle number. In what follows, for the sake of clarity, we present expressions for a system composed only of one kind of particle. Generalization to two types of particles, that is, to protons and neutrons, is straightforward and is discussed briefly later. Similarly to Ref. [14], the overlap and energy kernels are defined as I(φ)=|(φ),(2) H(φ)=|ˆ H|(φ),(3) and kernels of ( ˆ N−N0)mas Nm(φ)=|(ˆ N−N0)m|(φ).(4) The kernels of Eq. (4) can be calculated as derivatives of the overlap kernel with respect to the gauge angle Nm(φ)=(−i)mdm dφmI(φ).(5) Explicit expressions for these kernels are presented in Appendix A.InEqs.(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. [8]. As demonstrated by Lipkin [13], the minimized energy, obtained by the full variation after PNP (VAPNP), can also be obtained through an auxiliary Routhian, ˆ H=ˆ H−ˆ K{ˆ N−N0},(6) where the Lipkin operator ˆ K, which is a function of the shifted particle-number operator ˆ N−N0, is chosen so as to “flatten” the Ndependence of the average Routhians calculated for the particle-number projected states [13,14]. Had these projected Routhians been exactly Nindependent (perfectly flat), the exact projected energy EN0could have been obtained by minimizing the average value of the Routhian for the unprojected state |; that is, EN0=|ˆ H−ˆ K{ˆ N−N0}|.(7) Otherwise, the Lipkin method gives an approximate VAPNP energy, and its accuracy depends on the quality of the choice made for the Lipkin operator ˆ K. Similarly, after the Lipkin method is executed, the PNP of the final Lipkin state |gives an approximation to the exact VAPNP state. The advantage here is that the time-consuming exact PNP calculation is performed only once; that is, the Lipkin method allows us to obtain the full VAP result by effectively performing only the PAV calculation. Apart from the total energy, other observables must be calculated by using the PNP of the Lipkin state. As suggested by Lipkin [13], the simplest and most manageable ansatz for the Lipkin operator ˆ Khas the form of a power expansion, ˆ K{ˆ N−N0}= M  m=1 km(ˆ N−N0)m,(8) where k1≡λis the Fermi energy, which is used as a Lagrange multiplier to fix the average particle number. The higher-order Lipkin parameters kmfor m>1, which cannot be regarded as Lagrange multipliers, are used to best describe the particlenumber dependence of the average energies of projected states. Auxiliary equations are needed to determine these higher-order parameters. Up to now, the LN method has frequently been used to estimate values of k2(traditionally denoted by λ2). However, this method relies on calculating the average values of |ˆ Hˆ Nm|and |ˆ Nm|, 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 kmthat “flatten” the φdependence of the reduced Routhian kernel h(φ); 014312-2 LIPKIN METHOD OF PARTICLE-NUMBER RESTORATION . . . PHYSICAL REVIEW C 90, 014312 (2014) that is, h(φ)=h(φ)− M  m=1 kmnm(φ),(9) where h(φ)=H(φ) I(φ), h(φ)=H(φ) I(φ),(10) nm(φ)=Nm(φ) I(φ). 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 Nm(φ). Explicit expressions for kmare presented in Appendix A. The equivalency of the energy obtained by minimizing the auxiliary Routhian with that resulting from the exact VAPNP can be demonstrated as follows. In the HFB frame, the PNP state can be obtained in a standard way [2], N0≡ˆ PN0|= 1 2π2π 0 dφeiφ(ˆ N−N0)|,(11) where ˆ PN0is the projection operator for N0particles and | is the HFB wave function. 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 N0: E N0=|ˆ Hˆ PN0| |ˆ PN0|=2π 0H(φ)dφ 2π 0I(φ)dφ =2π 0h(φ)I(φ)dφ 2π 0I(φ)dφ =C. (12) Since for the state projected on N0, the average value of the Lipkin operator, (8), is, by definition, equal to 0, we also have that EN0=C, (13) 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. (12) come from the vicinity of the origin due to the largest weight [17]. Therefore, we can evaluate Lipkin parameters kmusing the gauge-rotated intrinsic states near the origin. This also avoids the singularities caused by vanishing overlaps [8]. As an example, at second order one obtains the Lipkin parameter, k2=h(φ2)−k1n1(φ2)−h(0) n2(φ2)−n2(0) ,(14) where φ2is a preselected small value of the gauge angle, and the flattened energy reads EN0=h(0)n2(φ2)−h(φ2)n2(0) +k1n1(φ2)n2(0) n2(φ2)−n2(0) .(15) Had the expansion up to second order been exact, the values of k2and EN0obtained from Eqs. (14) and (15) would have been independent of φ2. Thus, their eventual dependence on φ2indicates the necessity of going beyond second order. Similarly, at order M, we evaluate Lipkin parameters km, m=1,...,M, using a set of Msmall 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. If at the given order M, the convergence of the expansion of Lipkin operator, (8), is reached, the resulting parameters do not depend on the choice of the maximum gauge angle. We test the convergence based on this philosophy. The above derivations are strictly valid only in the case of energy kernels given by average values of the Hamiltonian. However, in the nuclear EDF approach, most often density-dependent interactions and interactions different in the particle-hole and particle-particle channels are used, and thus poles may occur when the overlaps between gauge rotated intrinsic states vanish (it may happen at gauge angle π/2) [4,8–10]. In this case, none of the standard methods, like VAPNP, PAV, LN, or Kamlah, or the Lipkin method discussed here, are strictly valid, and a construction of regularized functionals is mandatory [9,10]. In this sense, the Lipkin method, which employs appropriately small maximum gauge angles, not approaching the hypothetically dangerous region of π/2, can be regarded as a certain regularization method. By doing so, we regularize the energy kernels in terms of the analytic continuation of the Lipkin energy kernels to the full range of gauge angles. Obviously, at large gauge angles, the calculated and regularized energy kernels can then be different. Thus the tests of convergence of the Lipkin operator are meaningful only in the region of gauge angles where the energy kernels are not ill defined. 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 higherorder terms added (see Appendix A). We also note that the Lipkin parameters kmmust 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. the discussion of the LN and Kamlah methods in Ref. [19]). 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 Bshow that in such a schematic model, the higher-order Lipkin VAPNP method is able to reproduce correctly the exact VAPNP ground-state energies, in both weak and strong pairing regimes, everywhere apart from the immediate vicinity of the closed shell. This 014312-3 WANG, DOBACZEWSKI, KORTELAINEN, YU, AND STOITSOV PHYSICAL REVIEW C 90, 014312 (2014) gives us confidence in the application of this method to 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) [26,27]. This code solves the HFB equations in a three-dimensional Cartesian harmonic oscillator basis. Within this implementation, we tested the Lipkin VAPNP method using the Skyrme SIII parametrization [28] in the particle-hole channel and the volume zero-range pairing interaction in the particleparticle channel. SIII parametrization was selected for this study due to the fact that it contains only integer powers of densities. In the calculation of nondiagonal terms, the densitydependent interaction was treated within the “mixed density prescription,” as discussed in Refs. [5,8,29,30]. The neutron pairing strength, V0=−155.45 MeV fm3, was adjusted within the LN method to reproduce the empirical neutron pairing gap of n=1.245 MeV in 120Sn. 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, the pairing strength was set to 0; 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, we first study the convergence of the Lipkin VAPNP method when terms up to sixth order in expansion (8) are incorporated. At present, we limit our analysis to 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 N0.InFigs.1and 2, we show, for 120Sn and 100Sn, respectively, the dependence of the Lipkin parameters on the maximum gauge angle φM(see the previous section). The figures also show the total Lipkin VAPNP energy EN0and Lipkin correction energy Ecorr Ecorr =|− ˆ K{ˆ N−N0}|=− M  m=1 kmnm(0); (16) 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 EN0and correction energy Ecorr 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 100Sn, FIG. 1. (Color online) Lipkin parameters (a) k2,(b)k4,and (c) k6, (d) Lipkin correction energy Ecorr,(16), and (e) Lipkin VAPNP energy EN0,(7), determined in 120Sn at the second, fourth, and sixth orders, as functions of the maximum gauge angle φM. Note that Lipkin parameters k2,k4,andk6are shown in units of MeV, keV, and eV, respectively, which illustrates the rapid convergence of the Lipkin expansion. the convergence is slightly slower, and the values of Lipkin parameters are significantly higher than those for 120Sn. 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 the Lipkin operator, (8), in 100Sn and 120Sn 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 the fourth and sixth orders already give satisfactory description of the energy kernels. Figure 4shows the same kernels as plotted in Fig. 3for the whole range of gauge angle, up to φ=2π. We see that in 100Sn 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, FIG. 2. (Color online) Same as Fig. 1,butfor100Sn. 014312-4 LIPKIN METHOD OF PARTICLE-NUMBER RESTORATION . . . PHYSICAL REVIEW C 90, 014312 (2014) FIG. 3. (Color online) Reduced energy kernel h(φ)−h(0) (filled squares) and reduced kernels of the Lipkin operator M m=1km(nm(φ)−nm(0)) at orders M=2, 4, and 6 (open symbols), as functions of the gauge angle up to φ=1, calculated in 100Sn and 120Sn. which appears at the magic shell closure, and which cannot be properly described by a polynomial expansion [5,19]. In this case, a good quality of the Lipkin expansion obtained at small gauge angles is not sufficient to guarantee a good convergence at larger gauge angles. The situation is entirely different in the open-shell nucleus 120Sn, where the particlenumber dependence of the projected energies is given by a smooth function, which can be approximated very well 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 3s1/2orbital [8]. In Figs. 5and 6, we show the results of the Lipkin VAPNP method for tin and lead isotopes, respectively. For comparison, the figures also show the results obtained using the LN method, similarly as in Ref. [31], and the projected LN (PLN) method, as in Ref. [5], where the exact PNP energy is obtained via projection from the HFB +LN self-consistent solution. Away from the closed shells, at fourth and sixth orders, the results of the Lipkin VAPNP method are very similar to the PLN results. As pointed out by the VAPNP calculations in Ref. [5], for open-shell nuclei, PLN results are very close to exact VAPNP results. Again, the fourth and sixth orders give similar results, signaling the convergence of the Lipkin FIG. 4. (Color online) Same as Fig. 3, but for gauge angles up to φ=2π. FIG. 5. (Color online) LN, PLN, and Lipkin VAPNP energies of tin isotopes relative to those obtained within the standard HFB method. 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 nonconvergent power series of the Lipkin operator. Once again, this is related to the kink in the particle-number dependence on the projected energies [5,19]. In Fig. 7, we show results obtained by projecting good particle numbers from the states obtained either by the Lipkin VAPNP or the LN method. It is very gratifying to see that, irrespective whether one uses the Lipkin VAPNP or the LN method, the projected energies, shown in Figs. 7(a) and 7(e), 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 the fact that the PNP calculation does not have to be performed at all. Then, as shown in Figs. 7(b)–7(d) FIG. 6. (Color online) Same as Fig. 5, but for lead isotopes. 014312-5 WANG, DOBACZEWSKI, KORTELAINEN, YU, AND STOITSOV PHYSICAL REVIEW C 90, 014312 (2014) FIG. 7. (Color online) (a, e) PLN energies compared with exact PNP energies determined for states obtained by solving the Lipkin equations at second, fourth, and sixth orders. (b–d, f–h) Comparisons at three orders of the Lipkin VAPNP and exact PNP energies. All energies are plotted relative to those obtained within the standard HFB method. and 7(f)–7(h), the obtained energies approximate the PNP energies very well. 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 to the LN method (see Figs. 5and 6) gives results that are very different from those obtained by the PNP. On the contrary, the fourthand sixth-order Lipkin VAPNP method gives results almost identical to the PNP. Finally, at fourth and sixth orders, the nonanalytic behavior of the PNP energies at closed shells causes the largest discrepancies for two or four 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 the 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 this method is insufficient in some cases. This is illustrated in Fig. 8(a), which shows the reduced energy kernel of 124Xe calculated in two dimensions, as a function of the neutron φnand proton φpgauge angles. We clearly see that energy kernel is tilted with respect to the main axes of the neutron and proton gauge angles. Evidently, the Lipkin operator, here being the sum of the neutron and proton contributions separately, leads to a nontilted energy kernel, as shown in Fig. 8(b). Therefore, FIG. 8. (Color online) Reduced energy kernel h(φν,φπ)−h(0,0) (left) compared to the reduced kernel of the Lipkin operator M m=1km(nm(φν,φπ)−nm(0,0)) at sixth order (right), calculated for 124Xe. 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 the proton-neutron pairing scheme, combined with the PLN, these kinds of cross-terms are required [32,33]. Implementations of such more complicated forms of the Lipkin operator will be the subject of future study. IV. SUMMARY In the framework of the nuclear energy DFT, we have derived the Lipkin method of approximate particle-number symmetry restoration up to sixth order. The Lipkin parameters were determined from nondiagonal energy kernels, resulting in a more manageable approach compared to the traditional LN approach. Convergence of the Lipkin VAPNP method was tested by investigating the gauge dependence of expansion parameters. Taking 120Sn 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 the dependence on the gauge angle significantly. With the inclusion of sixth-order terms of the expansion, the overall change is minimal, indicating a converging series. The 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 100Sn and 120Sn show that the second-order Lipkin expansion is typically not sufficiently converged. Within the fourth order, the series already mimics the reduced energy kernels rather well. With inclusion of the sixth-order term, the results stay practically the same, again indicating a well-converged series. In chains of tin and lead isotopes, we have compared the Lipkin VAPNP method to the LN and PLN methods. As pointed out in Ref. [5], for midshell nuclei, the PLN is a very good approximation to the exact VAP method. Our results show that for midshell nuclei, the Lipkin VAPNP method already at 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 014312-6 LIPKIN METHOD OF PARTICLE-NUMBER RESTORATION . . . PHYSICAL REVIEW C 90, 014312 (2014) method is unable to reproduce the exact projected energy. Also, near the closed-shell region, the pairing correlations have a dynamic nature [2]. Within the Lipkin VAPNP method, this kind of feature 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 the necessity to include in the Lipkin operator cross terms, which depend simultaneously on the neutron and proton number operators. For the case of 124Xe, the contour lines of the reduced energy kernel, with respect to the neutron’s and proton’s gauge angles, show tilted shapes. Without the cross terms, this kind of behavior cannot be reproduced. A study of the cross terms will be the 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. We need to point out that this method cannot replace the exact projection for the case of, e.g., calculation of transition matrix elements or evaluation of selection rules, where good quantum numbers are mandatory. The Lipkin method can also be applied to approximate the restoration of other symmetries, broken at the mean-field level, at a small or no extra computational cost. This is important, for example, at the stage of adjusting the functionals, where rapid evaluations of penalty functions are required. This can be particularly significant for adjustments of EDFs tailored for beyond-mean-field multireference studies [34]. 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 Finland and University of Jyv¨ askyl¨ a within the FIDIPRO programme, by the Polish National Science Center under Contract No. 2012/07/B/ST2/03907, by the Academy of Finland under Centre of Excellence Programme 2012–2017 (Nuclear and Accelerator Based Physics Programme at JYFL), by the European Union’s Seventh Framework Programme ENSAR (THEXO) under Grant No. 262010, and by the U.S. Department of Energy under Contract Nos. DE-FC0209ER41583 (UNEDF SciDAC Collaboration) and DE-FG0296ER40963 (University of Tennessee). We acknowledge the CSC–IT Center for Science Ltd., Finland, for the allocation of computational resources. APPENDIX A: KERNELS AND LIPKIN PARAMETERS UP TO SIXTH ORDER To calculate kernels Nm(φ) from derivatives of the overlap kernel I(φ), in Eq. (5), we use the Onishi theorem [2]: |(φ)=det1/2(1 +e2iφCC+) det1/2(1 +CC+)e−iN0φ,(A1) where Cis the Thouless matrix. In the canonical basis of the Bogoliubov transformation, the above overlap is given as |(φ)= μ>0u2 μ+v2 μe2iφe−iN0φ,(A2) where the label μ>0 is associated with one of the pair-conjugated canonical states {μ, ¯μ}. Then the density matrix in the canonical basis can be obtained as ρμν(φ)=|a+ νaμ|(φ) |(φ)=e2iφv2 μ u2 μ+v2 μe2iφ δμν.(A3) We then have the reduced kernels nm(φ)=Nm(φ)/I(φ), given as n1(φ)=Trρ(φ)−N0≡R0(φ),(A4) n2(φ)=R2 0(φ)−id dφR0(φ)≡R2 0(φ)+R1(φ),(A5) n3(φ)=R3 0(φ)+3R0(φ)R1(φ)+R2(φ),(A6) n4(φ)=R4 0(φ)+6R2 0(φ)R1(φ)+4R0(φ)R2(φ)+3R2 1(φ)+R3(φ),(A7) n5(φ)=R5 0(φ)+10R3 0(φ)R1(φ)+10R2 0(φ)R2(φ)+15R0(φ)R2 1(φ)+5R0(φ)R3(φ)+10R1(φ)R2(φ)+R4(φ),(A8) n6(φ)=R6 0(φ)+15R4 0(φ)R1(φ)+20R3 0(φ)R2(φ)+45R2 0(φ)R2 1(φ)+15R2 0(φ)R3(φ)+60R0(φ)R1(φ)R2(φ) +6R0(φ)R4(φ)+15R3 1(φ)+15R1(φ)R3(φ)+10R2 2(φ)+R5(φ),(A9) where we used the definition Rn(φ)=(−i)ndn dφn(Trρ(φ)−N0).(A10) 014312-7 WANG, DOBACZEWSKI, KORTELAINEN, YU, AND STOITSOV PHYSICAL REVIEW C 90, 014312 (2014) Derivatives of the density matrix ρ(φ)(A3) can be calculated as, −id dφρ(φ)=2ρ(φ)(1 −ρ(φ)),(A11) (−i)2d2 dφ2ρ(φ)=4ρ(φ)(1 −ρ(φ))(1 −2ρ(φ)),(A12) (−i)3d3 dφ3ρ(φ)=8ρ(φ)(1 −ρ(φ)) −48ρ2(φ)(1 −ρ(φ))2,(A13) (−i)4d4 dφ4ρ(φ)=16ρ(φ)(1 −ρ(φ))2−16ρ2(φ)(1 −ρ(φ)) −192ρ2(φ)(1 −ρ(φ))3+192ρ3(φ)(1 −ρ(φ))2,(A14) (−i)5d5 dφ5ρ(φ)=32ρ(φ)(1 −ρ(φ))3−128ρ2(φ)(1 −ρ(φ))2+32ρ3(φ)(1 −ρ(φ)) +2304ρ3(φ)(1 −ρ(φ))3 −768ρ2(φ)(1 −ρ(φ))4−768ρ4(φ)(1 −ρ(φ))2,(A15) which gives R0(φ)=Trρ(φ)−N0,(A16) R1(φ)=2Trρ(φ)(1 −ρ(φ)),(A17) R2(φ)=4Trρ(φ)(1 −ρ(φ))(1 −2ρ(φ)),(A18) R3(φ)=8Trρ(φ)(1 −ρ(φ)) −48Trρ2(φ)(1 −ρ(φ))2,(A19) R4(φ)=16Trρ(φ)(1 −ρ(φ))2−16Trρ2(φ)(1 −ρ(φ)) −192Trρ2(φ)(1 −ρ(φ))3+192Trρ3(φ)(1 −ρ(φ))2,(A20) R5(φ)=32Trρ(φ)(1 −ρ(φ))3−128Trρ2(φ)(1 −ρ(φ))2+32Trρ3(φ)(1 −ρ(φ)) +2304Trρ3(φ)(1 −ρ(φ))3 −768Trρ2(φ)(1 −ρ(φ))4−768Trρ4(φ)(1 −ρ(φ))2.(A21) Up to the sixth order, the average Routhian of Eq. (7) to be minimized reads EN0=|ˆ H−k1(ˆ N−N0)−k2(ˆ N−N0)2−k3(ˆ N−N0)3−k4(ˆ N−N0)4−k5(ˆ N−N0)5−k6(ˆ N−N0)6|.(A22) Average values of powers of the particle-number operator are given in Eqs. (A4)–(A9) taken at φ=0. Moreover, in all terms with m⩾2, one can set R0≡0. This gives EN0=|ˆ H|−k1(Trρ−N0)−2k2Trρ(1 −ρ)−4k3Trρ(1 −ρ)(1 −2ρ)−12k4(Trρ(1 −ρ))2 −8k4Trρ(1 −ρ)+48k4Trρ2(1 −ρ)2−80k5(Trρ(1 −ρ))(Trρ(1 −ρ)(1 −2ρ)) −16k5Trρ(1 −ρ)2+16k5Trρ2(1 −ρ)+192k5Trρ2(1 −ρ)3−192k5Trρ3(1 −ρ)2 −120k6(Trρ(1 −ρ))3−240k6(Trρ(1 −ρ))2+1440k6(Trρ(1 −ρ))(Trρ2(1 −ρ)2) −160k6(Trρ(1 −ρ)(1 −2ρ))2−32k6Trρ(1 −ρ)3+128k6Trρ2(1 −ρ)2 −32k6Trρ3(1 −ρ)−2304k6Trρ3(1 −ρ)3+768k6Trρ2(1 −ρ)4+768k6Trρ4(1 −ρ)2.(A23) Hence, the corresponding mean-field Routhian (to be used in the HFB equations) reads1 h=h−k1−(2k2+24k4Trρ(1 −ρ)+8k4)(1 −2ρ)−4k3((1 −2ρ)2−2ρ(1 −ρ)) +96k4ρ(1 −ρ)(1 −2ρ)−80k5(Trρ(1 −ρ)(1 −2ρ))(1 −2ρ)−80k5(Trρ(1 −ρ))(1 −2ρ)2 +160k5(Trρ(1 −ρ))ρ(1 −ρ)−16k5(1 −ρ)2+64k5ρ(1 −ρ)−16k5ρ2+384k5ρ(1 −ρ)3 −1152k5ρ2(1 −ρ)2+384k5ρ3(1 −ρ)−32k6−360k6(Trρ(1 −ρ))2−320k6(Trρ(1 −ρ)(1 −2ρ)) −480k6(Trρ(1 −ρ)) +1440k6(Trρ2(1 −ρ)2)+1984k6ρ+720k6(Trρ(1 −ρ))2ρ+1920k6(Trρ(1 −ρ)(1 −2ρ))ρ +3840k6(Trρ(1 −ρ))ρ−2880k6(Trρ2(1 −ρ)2)ρ−17280k6ρ2−8640k6(Trρ(1 −ρ))ρ2 −1920k6(Trρ(1 −ρ)(1 −2ρ))ρ2+49 920k6ρ3+5760k6(Trρ(1 −ρ))ρ3−57 600k6ρ4+23040k6ρ5.(A24) 1Mean fields, (A24), and reduced energy kernels, (10), are traditionally denoted by the same symbol, h, but should not be confused with one another. Although the expression for the mean-field Routhian looks complicated, it is easy to compute. Additional terms appearing within the Lipkin method are simply composed of powers of density matrices, and the manipulations of these terms can be easily performed numerically. The particle-particle mean fields 014312-8