Treatment of the Intrinsic Hamiltonian in Particle-Number Nonconserving Theories

We discuss the implications of using an intrinsic Hamiltonian in theories without particle-number conservation, e.g., the Hartree-Fock-Bogoliubov approximation, where the Hamiltonian's particle-number dependence leads to discrepancies if one naively replaces the particle-number operator by its expectation value. We develop a systematic expansion that fixes this problem and leads to an a posteriori justification of the widely-used one- plus two-body form of the intrinsic kinetic energy in nuclear self-consistent field methods. The expansion's convergence properties as well as its practical applications are discussed for several sample nuclei.


Introduction
Symmetry breaking is a powerful concept in quantum manybody theory. In nuclear phyiscs, well-known examples are the breaking of the rotational and translational symmetries of the many-body Hamiltonian by the use of localized single-particle states in the construction of the many-body Hilbert space; the latter, in particular, can cause sizable center-of-mass contaminations of the energies and the many-body wave functions unless the symmetry is restored (see e.g. [1,2] and Refs. therein). In contrast, the breaking of the particle-number symmetry by quasi-particle methods like the Hartree-Fock-Bogoliubov (HFB) [1] approach is a useful tool because it leads to an efficient treatment of the important nuclear pairing correlations, although one will ultimately want to restore this symmetry in a finite system like the nucleus.
Since nuclei are self-bound objects the proper starting point of a nuclear many-body calculation is the translationally invariant intrinsic Hamiltonian The use of the intrinsic kinetic energy T int in a simple Hartree-Fock (HF) calculation has consequences for the validity of Koopmans' theorem, and thereby the interpretation of the HF eigenvalues as single-particle energies. A detailed analysis was given by Khadkikar and Kamble in Ref. [3] and referenced repeatedly over the past few decades (see e.g. [4] or our own work [5]). However, this analysis makes explicit use of the properties of the HF Slater determinant |Ψ , including the assumption that it is an eigenstate of the particle-number operator with eigenvalue A:Â |Ψ = A |Ψ .
Email addresses: Heiko.Hergert@physik.tu-darmstadt.de (H. Hergert), Robert.Roth@physik.tu-darmstadt.de (R. Roth) Naturally, this condition does not hold in a method like HFB, where the particle number is not conserved. We will analyze this case in detail in the following.

The General Case
Since we want to deal with theories that do not conserve particle number, we consider operators in Fock space. In this case, the intrinsic kinetic energy operator can be expressed either as a sum of one-and two-body operators, or a sum of two-body operators alone, i.e., the relative kinetic energies of each nucleon pair: The equality of these two expressions follows from the relation i< j where theÂ resulting from the summation over the second independent particle i or j in the first two terms is again a Fock space operator measuring the total particle number. This distinction is inconsequential as long as one works in a Hilbert space with fixed particle number, becauseÂ can then be replaced by the corresponding eigenvalue. Naturally, one is tempted to use a similar replacementÂ → Â in Fock space as well, but we will demonstrate in the following that this naive treatment of the particle-number dependence leads to discrepancies. Consider a many-body state |Ψ without fixed particle number. Taking the energy expectation value of the intrinsic Hamiltonian with Eqs. (3) or (4) in this state, we obtain the energy functionals and respectively. Since |Ψ does not satisfy the eigenvalue equation (2), we have to consider the operatorÂ −1 directly in all expectation values, which will require a series expansion in practice.
To this end, we note that where we have introduced Eq. (8) defines a formal expansion of the energy functionals in powers of Â −1 . Applying this expansion to E (a) , we obtain Denoting truncations of this series containing all terms up to a given order k by E (a) k , the leading (LO) and next-to-leading order (NLO) functionals are We note that the NLO functional is the one we would obtain by naively replacingÂ with Â in Eq. (6).
Plugging the expansion (8) into Eq. (7), we have In this case, a naive power counting in Â −1 breaks down because terms at a given order of the series are enhanced by factors of Â as a direct consequence of applying Eq. (5) to the expansion, i.e., and The term T that is formally of order Â 0 first appears in E (b) 1 , a linear term in E (b) 2 , and so on. Comparing E (a) 1 with E (b) 1 , we note that i.e., if we simply replaceÂ with the number Â in Eqs. (6) and (7), the functionals are no longer equivalent! To restore a proper power counting to the expansion of E (b) , we first note that the enhanced linear term appearing in E (b) 2 exactly cancels the one in Eq. (17). Likewise, an enhanced quadratic term in E (b) 3 cancels T (∆Â) 2 / Â 2 in Eq. (16), and similar cancellations occur for all higher orders. The cancellation can be enforced explicitly if we define and applying Eq. (5) we find that i.e., the functionals E (a) and E (b) are identical at any given order of the expansion.

Hartree-Fock-Bogoliubov Approximation
We now apply the expansion developed in the previous subsection to the HFB approximation [1], i.e., we assume that |Ψ is a quasi-particle Slater determinant and introduce the density matrix and the pairing tensor We first consider the one-plus two-body form of the intrinsic kinetic energy. At next-to-leading order, the functional reads where t is the single-particle kinetic energy, and we assume that all two-body matrix elements are antisymmetrized in the following. Varying E (a) 1 w.r.t. the density matrix, we obtain the particle-hole field where we have used Varying the energy with respect to the pairing tensor, we find that the pairing field is given by If we start from the pure two-body form (4) of the intrinsic kinetic energy, we follow the analysis in the previous section and apply (18) to obtain the NLO functional The expectation value of the correction term is but due to the properties of the Bogoliubov transformation [1] ρ − ρ 2 = −κκ * , and there is an ambiguity regarding to which field the correction term contributes when E (b) is varied. To make contact with related works in the literature (see [6] and Refs. therein) as well as the HF treatment of Ref. [3], we split it evenly between the particle-hole and particle-particle channels and express the latter contribution in a manifestly real form: We stress that the specific choice (29) only affects unobservable quantities like the particle-hole and pairing fields, while observables like the energy expectation value or the energy differences discussed in Sect. 3 do not depend on it at all. Noting that (cf. Eq.(15)), we find that the particle-hole field is given by and for the pairing field, we obtain Assuming that the single-particle states satisfy we obtain the following relation for the antisymmetrized matrix element of the relative kinetic energy: Plugging this into Eqs. (31) and (32), we find that for the NLO functionals E (a) 1 and E (b) 1 and the specific choice (29) for the correction term.
The equivalence of the fields guarantees that a solution of the HFB equations for the functional E (a) 1 will also solve the equations for E (b) 1 and vice versa: Here, H and R are the usual HFB Hamiltonian and generalized density matrices [1], and the Lagrange multiplier satisfies In this context, some additional remarks are in order. Since the expansion parameter Â −1 is directly linked to the variational degrees of freedom ρ kk ′ , the derivatives (24) generate Â −2 terms in the fields that cause global shifts in the diagonal matrix elements of h, i.e., the underlying single-particle spectrum. In contrast to Â −2 contributions that arise from varying the NNLO functionals E (a/b) 2 , these global shifts are stateindependent, and can be absorbed in the Lagrange multiplier λ in a self-consistent calculation. If they are included implicitly in λ, one cannot directly interpret λ as the Fermi energy of the system, as it is usually done in the literature (see e.g. [7]).
For higher orders of the expansion, the explicit evaluation of the expectation values T ∆Â n and p i · p j ∆Â n occurring in higher-order functionals E (a/b) n becomes increasingly cumbersome and time-consuming in practical calculations. Fortunately, the expansion of the energy functional converges rapidly, and it is sufficient to truncate the expansion of the energy functional at the linear order in practical calculations, as demonstrated in Sect. 3.

Hartree-Fock Limit
Starting from the analysis of the HFB approximation in the previous section, we can easily take the limit of vanishing pairing correlations to obtain the Hartree-Fock limit. Since the HF ground state |Ψ is an eigenstate ofÂ, we immediately see that and therefore all higher-order terms in the expanded functionals E (1) and E (2) automatically vanish. Of course, one could have directly used the eigenvalue equation (2) and avoided the expansion in the first place. Moreover, this implies that the "uncorrected" HF functional defined in analogy to (15) yields the same energies as E (a) and E (b) , while we obtain the relation for the corresponding particle-hole Hamiltonian by moving the "correction" terms appearing in h (b) to the left-hand side of Eq. (35). Equation (42) is exactly the relation Khadkikar and Kamble obtained by plugging (34) in the expression for h b [3]. From our analysis, we now see how this relation for the fields follows directly from the energy functional, and that the correction terms in particular are derivatives of vanishing energy contributions.
The HFB equations (37) for the functionals E (a) and E (b) are reduced to their HF counterparts, i.e., a HF solution ρ obtained with either functional is also a valid solution for the other functional. In addition, it was demonstrated in [3] that ρ also solves the HF equations for the uncorrected functional E (b) , i.e., can write down and solve HF-like equations for such a case, the EFA density matrix represents a statistical mixture rather than a genuine Slater determinant [8]. In this case, we cannot use the eigenvalue equation (2) to construct the energy functional, and we have to resort to the 1/ Â -expansion again, treating all expectation values in the statistical sense.

Particle-Number Projection
Since nuclei are finite systems, one eventually wants to carry out a particle-number projection (PNP) to obtain a state with fixed particle number. The PNP energy functional is constructed from a quasi-particle Slater determinant that is explicitly projected onto the particle number A via the hermitian projector P A |Ψ A = P A |Ψ (44) (see e.g. [1,9] for details). Since |Ψ A satisfies the eigenvalue equation (2), an expansion is not required, just as in the HF case. All three previously defined functionals are equivalent when derived from (44), but E (b) PNP will in general lead to different non-observable quantities like the projected fields h A and ∆ A or their individual contributions to the energy expectation value. In a variation after particle-number projection (VAP), the fields and densities obtained by solving the projected HFB equations [9] are associated with an auxiliary intrinsic state without physical meaning, whereas the expectation values of observables are well defined, even if they are A-dependent. Thus, only these expectation values of observables or derived quantities like separation energies (in the sense of energy differences) should be compared to experiment.

Discussion & Numerical Results
To test the proposed expansion, we have performed spherical HFB calculations using the phenomenological Gogny D1S interaction [10]. We are employing a spherical harmonic oscillator configuration space, and explicitly minimize all energies with respect to the oscillator length a HO ; more details can be found in Ref. [6]. Unless noted otherwise, we include 13 major oscillator shells in our calculations, which leads to a satisfactory convergence in the considered cases. Odd nuclei are treated in a self-consistent blocking method [1], where the odd nucleon is distributed evenly over all magnetic substates of a given j-shell according to the equal-filling approximation (see e.g. [8]).
Let us first consider the convergence of the 1/ Â -expansion. Formally, an A-quasi-particle HFB state can contain states with sharp particle numbers from 0, . . . , 2A. In the extreme cases, this would mean that the operator ∆Â/ Â in Eq. (8) could acquire the operator norm 1 on the space of HFB wavefunctions, causing the breakdown of the series expansion. In the nuclear many-body problem, this breakdown could only occur before the first major shell is fully occupied, and for these very light nuclei the use of mean-field methods is questionable in the first place. At the major shell closure itself, the HFB wavefunction collapses onto the HF solution, and there is no need for an expansion. As we progress along the nuclear chart, we find that particle-number fluctuations only occur within a given major shell, as shown exemplary for the tin isotopes with N = 50, . . . , 82 in Fig. 1. This implies that the operator norm of ∆Â/ Â remains below 1 in practical applications, and guarantees the convergence of the series.
To test the rate of the convergence, we display the quantity for the open-shell nuclei 18 O, 64 Ni, and 120 Sn in Fig. 2. We have picked these specific nuclei for their large values of ∆Â 2 in the respective isotopic chains. We find that is essentially sufficient to include the linear terms in 1/ Â in the energy functional. Beyond the linear order, the largest variations occur for 18 O, where the successive inclusion of terms up to third order causes changes of 100 − 200 keV (see the inset of Fig. 2). As expected, the effect of higher orders rapidly diminishes with increasing masses, and amounts to a few keV for the nickel and tin isotopes. For this reason, we will consider all functionals at linear (i.e. next-to-leading) order, and drop the subscripts in the following.
Information about the effect of the intrinsic kinetic energy on the nuclear pairing correlations can be extracted by considering the three-point binding energy difference formula along an isotopic chain (or equivalently, an isotonic chain with N replaced by Z). By calculating the ground states of odd nuclei self-consistently, we avoid the complications arising in perturbative analyses, as discussed in detail in [12,13]. This is particularly relevant for A-dependent Hamiltonians [6], which add another layer of complexity to the perturbative approach. In Fig. 3, we compare the theoretical ∆ (3) for the properly constructed functional to the uncorrected functional E (b) (15). The latter are typically lower than the values for the proper functional by 200 − 400 keV, with the exceptions occurring around the 1d 5/2 and 1d 3/2 sub-shell closures at N = 56 and N = 70, respectively. The gaps between the relevant sub-shell and the next available one are notably larger for E (b) than for E (a/b) , leading to more pronounced effects when the odd nucleon is added In Sect. 2.3, we have pointed out that a spherical HF treatment of open-shell nuclei will suffer from the same problems as the HFB method due to the use of the equal-filling approximation. This is explicitly demonstrated for the tin isotopic chain in Fig. 4, where we have used the density matrix ρ obtained by minimizing the functional E (a) HF to calculate the energies E (b) HF and E (b) HF . At sub-shell closures, we have a genuine HF problem, and in this case all three functionals are equivalent and ρ is a solution for each set of HF equations [3]. For open-shell nuclei, in contrast, the binding energies obtained with E (b) HF are reduced by several hundred keV, and the difference is given by to numerical accuracy. Finally, we compare the tin ground-state energies obtained with E (b) PNP to those of E (a/b) PNP in calculations where the variation is carried out after particle-number projection (VAP). We are aware that the VAP method is highly problematic for densitydependent interactions like Gogny D1S (see e.g. [14] and references therein). Since we do not aim for a comparison with experiment in the present discussion, we assume that the potential energies of both calculations will be equally affected by any spurious behavior due to the density-dependence, hence any differences in the VAP energy expectation values are due to the different forms of the kinetic energy. Looking at the top half of Fig. 5 and noting the eV scale, we find no such differences: the tin ground-state energies obtained in VAP calculations with either the properly constructed functional E (a/b) PNP or the naive, uncorrected functional E (b) PNP are identical, confirming the discussion of Sect. 2.4. At the same time, the comparison of the pairing energies in the bottom half of Fig. 5 illustrates that the individual particle-hole and particle-particle contributions to the total energy expectation value may well be different in both calculations.

Conclusions
We have presented a detailed analysis of theÂ-dependent intrinsic kinetic energy operator in theories with and without particle-number conservation, in particular the mean-field type Hartree-Fock-Bogoliubov method and its number-conserving extension via particle-number projection. We have shown that a naive treatment where the number operatorÂ is replaced by its expectation value Â causes discrepancies in expectation values obtained with the otherwise equivalent operator forms of T int . We have developed a systematic expansion to fix this problem, but we emphasize that this expansion does not restore, nor is it intended to restore, either the particle number or translational symmetries of the nucleus.
Our discussion provides an a posteriori justification for using the one-plus two-body form of the intrinsic kinetic energy since it is automatically consistent with the power counting of the developed expansion. As a byproduct, we also clarify how differences and ambiguities in non-observable quantities which had been discussed in the context of the HF approximation [3] arise systematically in the presented framework. While we have discussed the specific case of the intrinsic kinetic energy operator in the present article, we point out that the same treatment should be applied to allÂ-dependent observables.