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]

$$\vert\Phi(\phi)\rangle = \exp\left(i\phi(\hat{N}-N_0)\right)\vert\Phi\rangle \, ,$$ (1)

where $\phi$ is the gauge angle, $\hat{N}$ is the particle-number operator, and $N_0=\langle\Phi\vert\hat{N}\vert\Phi\rangle$ 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. [14], the overlap and energy kernels are defined as

$$I(\phi) = \langle\Phi\vert\Phi(\phi)\rangle \, ,$$ (2)

$$H(\phi) = \langle\Phi\vert\hat{H}\vert\Phi(\phi)\rangle \, ,$$ (3)

and kernels of $(\hat{N}-N_0)^m$ as

$$N_m(\phi) =\langle\Phi\vert(\hat{N}-N_0)^m\vert\Phi(\phi)\rangle \, .$$ (4)

The kernels of Eq. (4) can be calculated as derivatives of the overlap kernel with respect to the gauge angle

$$N_m(\phi) = (-i)^m\frac{d^m}{d\phi^m}I(\phi) \, .$$ (5)

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. [8].

As demonstrated by Lipkin [13], the minimized energy, obtained by the full variation after the particle-number projection (VAPNP), can also be obtained through an auxiliary Routhian,

$$\hat{H}' = \hat{H}-\hat{K}\{\hat{N}-N_0\} \, ,$$ (6)

where Lipkin operator $\hat{K}$, which is a function of the shifted particle-number operator $\hat{N}-N_0$, is chosen so as to ``flatten'' the $N$-dependence of average Routhians calculated for the particle-number projected states [13,14]. Had these projected Routhians been exactly $N$-independent (perfectly flat), the exact projected energy $E_{N_0}$ could have been obtained by minimizing the average value of the Routhian for the unprojected state $\vert\Phi\rangle$, that is,

$$E_{N_0}=\langle\Phi\vert\hat{H}-\hat{K}\{\hat{N}-N_0\}\vert\Phi\rangle \, .$$ (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 $\hat{K}$.

Similarly, after the Lipkin method is executed, the particle-number-projection (PNP) of the final Lipkin state $\vert\Phi\rangle$ 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 for obtaining 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 manageable ansatz for the Lipkin operator $\hat{K}$ has the form of a power expansion,

$$\hat{K}\{\hat{N}-N_0\}= \sum_{m=1}^M k_{m} (\hat{N}-N_0)^m \, ,$$ (8)

where $k_{1}\equiv\lambda$ is the Fermi energy, which is used as a Lagrange multiplier to fix the average particle number. The higher-order Lipkin parameters $k_{m}$ for $m>1$, which cannot be regarded as Lagrange multipliers, are used to best describe the particle-number dependence of the average energies of projected states. Auxiliary equations are needed to determine these higher-order parameters.

Up to now, the LN method was frequently used to estimate values of $k_{2}$ (traditionally denoted by $\lambda_{2}$). However, this method relies on calculating the average values of $\langle\Phi\vert{\hat H} \hat{N}^m\vert\Phi\rangle$ and $\langle\Phi\vert\hat{N}^m\vert\Phi\rangle$, 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 $\phi$-dependence of the reduced Routhian kernel $h'(\phi)$, that is,

$$h'(\phi) = h(\phi) - \sum_{m=1}^M k_{m} n_m(\phi) \, ,$$ (9)

where

$$h' (\phi)=\frac{H' (\phi)}{I(\phi)} \, , ~ h (\phi)=\fra... ...i)}{I(\phi)} \, , ~ n_m(\phi)=\frac{N_m(\phi)}{I(\phi)} \, .$$ (10)

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(\phi)$ and $N_m(\phi)$. Explicit expressions for $k_{m}$ are 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]

$$\vert\Psi_{N_0}\rangle \equiv \hat{P}_{N_0}\vert\Phi\rangle= \frac{1}{2\pi}\int_0^{2\pi} d \phi e^{i \phi (\hat{N}-N_0)}\vert\Phi\rangle \, ,$$ (11)

where $\hat{P}_{N_0}$ is the projection operator for $N_0$ particles and $\vert\Phi\rangle$ is the HFB wave function. For a perfectly flat ($\phi$-independent) reduced Routhian kernel $h'(\phi)\equiv{C}$, we then have the exact average value of the Routhian evaluated for the state projected on particle number $N_0$,

$$E'_{N_0} = \frac{\langle \Phi \vert\hat{ H}' \hat{P}_{N_0} \vert \Phi \rangl... ...} = \frac{\int_0^{2\pi} H'(\phi){\rm d}\phi} {\int_0^{2\pi} I(\phi){\rm d}\phi}$$

$$= \frac{\int_0^{2\pi} h'(\phi)I(\phi){\rm d}\phi } {\int_0^{2\pi} I(\phi){\rm d}\phi} = C \, .$$ (12)

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

$$E_{N_0} = C,$$ (13)

and thus the minimization of the average Routhian (7) is equivalent to the exact VAPNP. Again, any imperfection in the $\phi$-independence of $h'(\phi)$ 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 $k_{m}$ using 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,

$$k_{2}=\frac{h(\phi_2)-k_{1}n_1(\phi_2)-h(0)}{n_2(\phi_2)-n_2(0)} \, ,$$ (14)

where $\phi_2$ is a pre-selected small value of the gauge angle, and the flattened energy reads

$$E_{N_0} =\frac{h(0)n_2(\phi_2)-h(\phi_2)n_2(0)+k_{1}n_1(\phi_2)n_2(0)}{n_2(\phi_2)-n_2(0)} \, .$$ (15)

Had the expansion up to second order been exact, values of $k_{2}$ and $E_{N_0}$ obtained from Eqs. (14) and (15) would have been independent of $\phi_2$. Thus, their eventual dependence on $\phi_2$ indicates the necessity of going beyond second order.

Similarly, at order $M$, we evaluate Lipkin parameters $k_{m}$, $m=1,\ldots,M$, using a set of $M$ small gauge angles $\phi_i$, $i=1,\ldots,M$. In practice, in this work, we use equally spaced values of $\phi_i=i\phi_1$, and at each order we check the eventual dependence of results on the maximum gauge angle used, $\phi_{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 $\pi/2$) [4,8,9,10]. In such a case, none of the standard methods, like VAPNP, PAV, LN, or Kamlah, nor the Lipkin method discussed here, are strictly valid, and a construction of regularized functionals is mandatory [9,10].

In this sense, the Lipkin method that employs appropriately small maximum gauge angles, which do not approach the hypothetically dangerous region of $\pi/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 angels, 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 $\vert\Phi\rangle$ 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 $\vert\Phi\rangle$), 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. [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 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.