Application of the variational principle to a coherent-pair condensate: The HFB case
L. Y. Jia






  liyuan.jia@usst.edu.cn
Department of Physics, University of Shanghai for
Science and Technology, Shanghai 200093, P. R. China










Recently we proposed a scheme that applies the variational principle to a coherent-pair condensate in the BCS case <cit.>. This work extends the scheme to the HFB case by allowing variation of the canonical single-particle basis. The result is equivalent to that of the so-called variation after particle-number projection in the HFB case, but now the particle number is always conserved and the time-consuming projection is avoided. Specifically, we derive the analytical expression for the gradient of the average energy with respect to changes of the canonical basis, which is then used in the family of gradient minimizers to iteratively minimize energy. In practice, we find the so-called ADAM minimizer (adaptive moments of gradient), borrowed from the machine-learning field, is very effective. The new algorithm is demonstrated in a semi-realistic example using the realistic V_low-k interaction and large model spaces (up to 15 harmonic-oscillator major shells). It easily runs on a laptop, and practically the computer time cost (in the HFB case) is several times that of solving HF by the gradient minimizers. We hope the new algorithm could become a common practice because of its simplicity.








Application of the variational principle to a coherent-pair condensate: The HFB case
    L. Y. Jia
    December 30, 2023
====================================================================================







§ INTRODUCTION





The self-consistent mean-field theory is widely used in nuclear structure <cit.>. Its simplicity allows application to the whole nuclear chart, and the self-consistent single-particle levels are the starting point for nuclear shell model (configuration interaction), especially for nuclei far from stability where the evolution of single-particle levels is unknown a priori. The Hartree-Fock (HF) method is the microscopic mean-field theory.

The collective (coherent) pairing effect <cit.> is important across the nuclear chart <cit.>. To incorporate pairing into the mean field, one introduces quasiparticles and the microscopic theory is the Hartree-Fock-Bogoliubov (HFB) method <cit.>. HFB is simply the variational principle using the quasiparticle vacuum as the trial wavefunction <cit.>. Quasiparticles are highly successful and wildly used. As fermions, they have clear physical picture and allow easy computation. However, quasiparticles break the exact particle number <cit.>. Only the average particle number is guaranteed by the chemical potential; effectively, one replaces the target nucleus by an average of neighbouring nuclei. This average may cause errors, especially in phase transition regions of sharp property changes.

To cue the problems, one projects the HFB quasiparticle vacuum onto good particle number; the standard way to project is by numerical integration over the gauge angle <cit.>. The projection can be done after or before the variation <cit.>. For examples of projection after variation (PAV), see Refs. <cit.> and references therein. For examples of variation after projection (VAP), see Refs. <cit.> in the BCS case and Refs. <cit.> in the HFB case. When feasible, VAP is preferred <cit.> over PAV. The practical difficulty of VAP is that numerical projection by integration is time-consuming <cit.> and needed many times when performing VAP. In the literature there are far fewer realistic applications of VAP+HFB than those of the HFB theory without projection.



It is easy to analytically project <cit.> the HFB quasiparticle vacuum onto good particle number — the result is a coherent-pair condensate [see Eq. (<ref>)]. This work applies the variational principle directly to this coherent-pair condensate (VDPC). The particle number is always conserved and the time-consuming projection is avoided. (This feature is emphasized by the word “directly” in the name “VDPC”. We name the new method VDPC because VAP may be misleading: there is no projection at all.) The variational parameters are the canonical single-particle basis and the coherent-pair structure v_α on this basis. Ref. <cit.> has proposed a fast algorithm for varying v_α at fixed canonical basis, which is VDPC+BCS and the result is equivalent to that of VAP+BCS. This work considers VDPC+HFB that varies v_α and the canonical basis together, and the result is equivalent to that of VAP+HFB.



Specifically, we first consider the subproblem of varying the canonical basis at fixed v_α, which we call VDPC+CB (plus canonical basis). This generalizes the HF theory, which uses special values of v_α: v_α = 1 for occupied levels in the HF Slater determinant and v_α = 0 for empty levels. To solve VDPC+CB, we derive the analytical expression for the gradient of the average energy with respect to changes of the canonical basis, which is then used in the family of gradient minimizers <cit.> to iteratively minimize energy. In practice, we find the so-called ADAM minimizer (adaptive moments of gradient) <cit.>, borrowed from the machine-learning field <cit.>, is very effective (energy converges quicker than other minimizers).

The new VDPC+HFB algorithm combines VDPC+BCS and VDPC+CB. Specifically, we insert varying v_α (VDPC+BCS) into several places in the process of varying the canonical basis (VDPC+CB). (It is unnecessary to vary v_α after every iteration of VDPC+CB; after every 20 iterations, for example, is enough.)

We demonstrate the new VDPC+CB and VDPC+HFB algorithms in a semi-realistic example using the realistic V_low-k interaction <cit.> and large model spaces (up to 15 harmonic-oscillator major shells). (The same example/Hamiltonian has been used in Ref. <cit.> for VDPC+BCS.) The energy-convergence pattern and actual computer time cost are given in detail. The new algorithms easily run on a laptop, and practically the computer time cost to solve VDPC+HFB (time to solve VDPC+CB) is about 4 ∼ 6 times (2 ∼ 3 times) that to solve HF by the gradient minimizers. We hope the new VDPC+HFB algorithm could become a common practice because of its simplicity.



This work relates to Refs. <cit.>. The average energy of the coherent-pair condensate (<ref>) has been derived for the pairing Hamiltonian <cit.> and a general Hamiltonian <cit.>. However, the gradient of energy with respect to v_α and to changes of the canonical basis have not been derived. (The gradient with respect to v_α was recently derived in Ref. <cit.>.) For successes and limitations of these works <cit.>, please see the introduction part of Ref. <cit.>. For modern mean-field theories using large model spaces in the pairing channel, VAP has been done only by the numerical gauge-angle integration <cit.>. This work aims to propose the new VDPC algorithm.


We should also mention the Lipkin-Nogami prescription to restore the particle number approximately <cit.>. It is widely used because the exact VAP (by the numerical gauge-angle integration) is computationally expensive <cit.>. Ongoing efforts exist to improve the Lipkin method <cit.>.

This work is organized as follows. Section <ref> briefly reviews
the formalism for the coherent-pair condensate — the trial wavefunction. Section <ref> derives analytically the gradient of the average energy with respect to changes of the canonical single-particle basis. What this gradient expression becomes when reducing the valence space in the pairing channel (make it smaller than the full space) is discussed in Sec. <ref>. We explain how to solve VDPC+CB by the family of iterative gradient minimizers in Sec. <ref>, and describe the VDPC+HFB algorithm in Sec. <ref>. Section <ref> applies the new VDPC algorithms to a semi-realistic example. Section <ref> summarizes the work.









§ COHERENT-PAIR CONDENSATE AS TRIAL WAVEFUNCTION  


This work uses the coherent-pair condensate as the trial wavefunction in the variational principle. This section briefly reviews the relevant definitions and formulas. For clarity we consider one kind of nucleon, the extension to active protons and neutrons is simple: the existence of protons simply provides a correction to the neutron single-particle energy, through the two-body proton-neutron interaction. We assume time-reversal self-consistent symmetry <cit.>, and the single-particle basis state |α⟩ is Kramers degenerate with its time-reversed partner |α̃⟩ (|α̃̃̃⟩ = - |α⟩). No other symmetry is assumed.

For the ground state of the 2N-particle system, the trial wavefunction is an N-pair condensate,

    |ϕ_N⟩ = 1/√(χ_N) (P^†)^N |0⟩
    ,

where

    χ_N = ⟨0| P^N (P^†)^N |0⟩

is the normalization factor, and P^† creates a coherent pair

    P^† = 1/2∑_αβ v_αβ
    a_α^† a_β̃^† .

Requiring |ϕ_N⟩ to be time-even implies that the pair
structure matrix v_αβ = v_βα^* is Hermitian (because T̂
v_αβ a_α^† a_β̃^† |0⟩
= v_αβ^* a_α̃^† a_β̃̃̃^†
|0⟩ = v_αβ^* a_β^† a_α̃^†
|0⟩). The unitary transformation to the canonical single-particle basis diagonalizes v_αβ. Below we always use the canonical basis, thus v_αβ = δ_αβ v_α, v_α is real, and the coherent pair (<ref>) becomes

    P^† = 1/2∑_α v_α
    a_α^† a_α̃^† = ∑_α∈Θ v_α
    P^†_α ,

where

    P^†_α = a_α^† a_α̃^† = P^†_α̃

creates a pair on |α⟩ and |α̃⟩.
In Eq. (<ref>), Θ is the set of pair-indices that picks only one from each degenerate pair |α⟩ and |α̃⟩. For example, with axial symmetry we can choose Θ to be orbits of a positive magnetic quantum number. ∑_α and ∑_α∈Θ means summing over single-particle indices and pair indices.



The trial wavefunction (<ref>) is specified by two sets of variational parameters: the canonical single-particle basis and the pair structure v_α (<ref>) on this basis. Ref. <cit.> has proposed a fast algorithm for varying v_α at fixed canonical basis that is VDPC+BCS. This work first considers varying the canonical basis at fixed v_α in Sec. <ref>, which is VDPC+CB. Then Sec. <ref> considers varying the canonical basis and v_α together that is VDPC+HFB.

If one insists time-reversal self-consistent symmetry <cit.>, and analytically projects the time-even HFB quasiparticle vacuum onto good particle number, one gets the same trial wavefunction (<ref>) (see Eq. (7.18) of Ref. <cit.>). Therefore, the result of VDPC+HFB is equivalent to that of VAP+HFB, and the result of VDPC+BCS <cit.> is equivalent to that of VAP+BCS. The difference is that VDPC always conserves the particle number and avoids the time-consuming numerical projection (gauge-angle integration) of VAP.

The coherent-pair condensate (<ref>) includes the HF Slater determinant as a special case, when v_α is fixed to 1 for the 2N occupied HF orbits and to 0 for other empty HF orbits. From this perspective, VDPC extends the HF theory to incorporate pairing correlations, but not by introducing quasi-particles that break particle number (what HFB does); instead, VDPC always converses the particle number.



For convenience, we introduce {[γ_1 γ_2 …γ_r]} to represent a subspace of the original single-particle space, by removing Kramers pairs of single-particle levels γ_1, γ̃_1, γ_2, γ̃_2, …γ_r, γ̃_r from the latter. Later we will express the gradient of energy in terms of Pauli-blocked normalizations <cit.>,

    χ_N^[γ_1 γ_2 …γ_r]≡⟨0| P^N P_γ_1 P_γ_2… P_γ_r P_γ_1^† P_γ_2^†… P_γ_r^† (P^†)^N |0⟩ ,

which is the normalization in the blocked subspace {[γ_1 γ_2 …γ_r]}. Given v_α, how to compute χ_N^[γ_1 γ_2 …γ_r]? It is explained in detail in Ref. <cit.>. For example, χ_N and χ_N^[α] are computed by recursive relations,

    χ_N = N ∑_α∈Θ (v_α)^2 χ_N-1^[α] ,  
    χ_N - χ_N^[α] = (N v_α)^2 χ_N-1^[α] = χ_N ⟨ϕ_N | n̂_α | ϕ_N ⟩ ,

with initial value χ_N=0^[α] = 1. Knowing χ_N-1^[α]'s, we compute χ_N by Eq. (<ref>), and then χ_N^[α]'s by Eq. (<ref>). ⟨ϕ_N | n̂_α | ϕ_N ⟩ = ⟨ 0 | P^N a_α^† a_α (P^†)^N | 0 ⟩ / χ_N is the average occupation number. Equations (<ref>) and (<ref>) are also valid in the blocked subspaces {[γ_1 γ_2 …γ_r]}, which could be used to compute χ_N^[αβ] and χ_N^[αβγ]; but another method is better (by Eqs. (15) and (16) of Ref. <cit.>). Later we will also need

    | ϕ_N^[β]⟩≡1/√(χ_N^[β]) (P^† - v_β P_β^†)^N |0⟩

that is the pair condensate with β and β̃ blocked.

This section discusses the “kinematics” of the VDPC formalism, next we discuss the “dynamics”.








§ GRADIENT OF ENERGY  


In this section we derive the gradient of the average energy with respect to changes of the canonical single-particle basis. The anti-symmetrized two-body Hamiltonian is

    H = ∑_αβϵ_αβ a_α^† a_β + 1/4∑_αβγμ V_αβγμ a_α^† a_β^† a_γ a_μ .

Note the ordering of αβγμ, thus V_αβγμ = - ⟨αβ| V | γμ⟩.
I assume time-even H (ϵ_αβ = ϵ_β̃α̃, V_αβγμ = V_μ̃γ̃β̃α̃) and real ϵ_αβ and V_αβγμ. No other symmetry of H is assumed. The Hamiltonian matrix elements ϵ_αβ and V_αβγμ have been transformed to the canonical basis (α, β, γ, μ run over the canonical-basis states).

The average energy of the coherent-pair condensate E̅ = ⟨ϕ_N | H | ϕ_N ⟩ has been derived in Eq. (25) of Ref. <cit.>,

    ⟨ϕ_N | H | ϕ_N ⟩ = N^2/χ_N( ∑_α∈Θ (2ϵ_αα + G_αα) (v_α)^2 χ_N-1^[α]
    
    + ∑^αβ_α,β∈Θ G_αβ v_α v_βχ_N-1^[αβ]
    
    + (N-1)^2 ∑^αβ_α,β∈ΘΛ_αβ (v_α v_β)^2 χ_N-2^[αβ]) ,

where we introduce the paring matrix elements G_αβ and the “monopole” matrix elements Λ_αβ as

    G_αβ≡ V_αα̃β̃β ,  
    Λ_αβ≡ V_αββα + V_αβ̃β̃α .

Note G_αβ = G_βα = G_αβ̃, Λ_αβ = Λ_βα = Λ_αβ̃, and G_αα = Λ_αα.

The gradient of the average energy with respect to v_α (at fixed canonical basis), ∂E̅/∂ v_α = ∂[ ⟨ϕ_N | H | ϕ_N ⟩ ]/∂ v_α, has been derived in Eqs. (31) and (32) of Ref. <cit.>. Based on this gradient, Ref. <cit.> proposed a fast algorithm for varying v_α at fixed canonical basis, which is VDPC+BCS.







Now we derive the gradient of the average energy with respect to changes of the canonical basis (at fixed v_α). For this purpose we first parameterize changes of the single-particle basis. We have assumed the Hamiltonian matrix elements ϵ_αβ and V_αβγμ (<ref>) are real, so the eigen wavefunctions are also real. Therefore we restrict the trial wavefunction (<ref>) (approximates the ground state) to be real, that is, v_αβ in Eq. (<ref>) to be real. Thus the unitary transformation to the canonical basis (<ref>) becomes an orthogonal transformation (a real unitary matrix). We parameterize the mixing of two single-particle basis states as (real mixing coefficients)

    |α'⟩ = cosθ |α⟩ + sinθ |β⟩ ,  |β'⟩ = cosθ |β⟩ - sinθ |α⟩ ,

consequently for their time-reversal partners

    |α'⟩ = cosθ |α̃⟩ + sinθ |β̃⟩ ,  |β'⟩ = cosθ |β̃⟩ - sinθ |α̃⟩ .

|α⟩ and |β⟩ belong to different Kramers pairs (P_α P_β). Mixing |α⟩ and |α̃⟩ has no effect on P_α^† = a_α^† a_α̃^† (<ref>), so P^† (<ref>) and |ϕ_N⟩ (<ref>) stay unchanged. If one requires additional self-consistent symmetry (for example, parity and axial symmetry), |α⟩ and |β⟩ have the same values of the corresponding good quantum numbers (parity and angular momentum projection onto the symmetry axis). For infinitesimal mixing (infinitesimal θ≈ 0), keeping the first-order in θ, Eqs. (<ref>) and (<ref>) imply variations of the basis states,

    δ |α⟩ = |α'⟩ - |α⟩≈θ |β⟩ , δ |β⟩ = |β'⟩ - |β⟩≈ - θ |α⟩ , 
    δ |α̃⟩ = |α'⟩ - |α̃⟩≈θ |β̃⟩ , δ |β̃⟩ = |β'⟩ - |β̃⟩≈ - θ |α̃⟩ .




Variation of the average energy E̅ (<ref>) under an infinitesimal θ has been derived in Sec. V of the arXiv manuscript <cit.> (but has not been published in a journal yet). We repeat the result here: 

    δE̅ = δ (⟨ϕ_N | H | ϕ_N ⟩) = 4 θ f_αβ ,




where


    f_αβ = N^2 (v_α - v_β)/χ_N ( [ (v_α + v_β) ϵ_αβ + v_α V_αβ̃α̃α
    + v_β V_αβ̃β̃β ] χ_N-1^[αβ]
    
    + ∑_γ∈Θ^P_γ P_α, P_β v_γ V_αβ̃γ̃γ [ χ_N-1^[αβγ] - (N-1)^2 v_α v_βχ_N-2^[αβγ]]  
    
    + (N-1)^2 (v_α + v_β) ∑_γ∈Θ^P_γ P_α, P_β (v_γ)^2 ( V_αγγβ + V_αγ̃γ̃β ) χ_N-2^[αβγ]) .

f_αβ = - f_βα is skew-symmetric. We pull out the factor 4 when defining f_αβ in Eq. (<ref>), so that f_αβ reduces to the off-diagonal part of the HF mean field when |ϕ_N⟩ (<ref>) reduces to a HF Slater determinant, as shown in Sec. <ref>. Equation (<ref>) means that 4 f_αβ is the partial derivative at θ = 0,

    ∂E̅/∂θ |_θ = 0 = ∂ (⟨ϕ_N | H | ϕ_N ⟩)/∂θ |_θ = 0 = 4 f_αβ ,

where θ is the angle mixes the two canonical basis states |α⟩ and |β⟩ according to Eq. (<ref>). Using (N-1)^2 (v_γ)^2 χ_N-2^[αβγ] = χ_N-1^[αβ]⟨ϕ_N-1^[αβ] | n̂_γ | ϕ_N-1^[αβ]⟩
[Eq. (<ref>) with N → N-1, α→γ, then blocking the α and β Kramers pairs], an equivalent form of f_αβ is

    f_αβ = N^2 (v_α - v_β)/χ_N ( [ (v_α + v_β) (ϵ_αβ + ∑_γ^P_γ P_α, P_β V_αγγβ⟨ϕ_N-1^[αβ] | n̂_γ | ϕ_N-1^[αβ]⟩ ) + v_α V_αβ̃α̃α
    + v_β V_αβ̃β̃β ] χ_N-1^[αβ]
    
    + 1/2∑_γ^P_γ P_α, P_β v_γ V_αβ̃γ̃γ [ χ_N-1^[αβγ] - (N-1)^2 v_α v_βχ_N-2^[αβγ]]
     ) ,


where γ sums over single-particle index. In Eq. (<ref>), the ∑_γ V_αγγβ⟨ϕ_N-1^[αβ] | n̂_γ | ϕ_N-1^[αβ]⟩ term is the monopole correction to ϵ_αβ. 

In this section, we derive the gradient (<ref>) when assuming all the single-particle levels (canonical basis states) participate in the coherent pairing (<ref>). Usually, only the levels near the Fermi surface are important, and they compose the active valence space in the pairing channel (smaller dimension than the full space). In the next section, we consider what the gradient (<ref>) becomes when using such a limited valence space in the pairing channel.








§ REDUCING VALENCE SPACE IN PAIRING CHANNEL  


In the previous section, we derive the energy gradient expression (<ref>) when assuming all the single-particle levels (canonical basis states) participate in the coherent pairing (<ref>). In this section, we consider what this gradient expression becomes when reducing the valence space in the pairing channel.

Usually, only the single-particle levels (canonical basis states) near the Fermi energy E_F  actively participate in the coherent pairing (<ref>) and are partially filled. The levels well below and well above E_F are almost full and empty (n_α≈ 1 and 0). So it is a good approximation to keep these levels completely full and empty (n_α = 1 and 0), and restrict the coherent pairing (<ref>) to those levels near the Fermi surface. The latter compose the active valence space in the pairing channel, which should be large enough (include enough levels) for a desired accuracy of the approximation.



Specifically, we divide the full single-particle space (canonical basis) into three subspaces: L (low), V (valence), and H (high). For α∈ H, we set v_α = 0 in Eq. (<ref>), so n_α = 0. For α∈ L, we set v_α = v_L and will take the limit v_L →∞, so n_α→ 1. For α∈ V, v_α is a finite nonzero number, so 0 < n_α < 1. The criterion for division (into L, V, H subspaces) is not necessarily by single-particle energy; in fact, this work will use the criterion by occupation number n_α. We assume L has dimension 2 N_L, thus the valence particle number in V is 2(N-N_L) ≡ 2 N̅. We introduce the V-subspace coherent pair

    P̅^†≡1/2∑_α∈ V v_α
    P^†_α = ∑_α∈Θ, α∈ V v_α
    P^†_α ,

and the V-subspace coherent-pair condensate

    |ϕ̅_N̅⟩≡1/√(χ̅_N̅) (P̅^†)^N̅ |0⟩ ,

where

    χ̅_N̅≡⟨0| P̅^N̅ (P̅^†)^N̅ |0⟩

is the V-subspace normalization. In summary, we hat the V-subspace symbols by a bar (not to be confused with the average energy E̅ = ⟨ϕ_N | H | ϕ_N ⟩).

It is easy to derive the relations between the full-space quantities and the V-subspace ones. For the coherent pair (<ref>),

    P^† = P̅^† + v_L∑_α∈Θ, α∈ L
    P^†_α .

So for (P^†)^N = (P^†)^N_L + N̅,

    (P^†)^N = A_N^N_L (v_L)^N_L (P̅^†)^N̅ (∏_α∈Θ, α∈ L P_α^†)  
    
    + O[(v_L)^N_L-1] ,

where O[(v_L)^N_L-1] represents terms of the power N_L-1 and lower. A_N^N_L = N! / (N-N_L)! = N! / N̅! is the number of permutations, for selecting the N_L of P_α∈ L^† operators from the N of multiplying P^†. ∏_α∈Θ, α∈ L P_α^† equals ∏_α∈ L a_α^† within a sign that will fully occupy the L-subspace when acting on |0⟩. Using Eq. (<ref>), the normalization (<ref>) is

    χ_N = (A_N^N_L)^2 (v_L)^2N_Lχ̅_N̅ + O[(v_L)^2N_L-1] .

Using Eqs. (<ref>) and (<ref>), the coherent-pair condensate (<ref>) is

    |ϕ_N⟩ = 1/√(χ̅_N̅) (P̅^†)^N̅ (∏_α∈Θ, α∈ L P_α^†) |0⟩ + O[(v_L)^-1]  
    
    = (∏_α∈Θ, α∈ L P_α^†) |ϕ̅_N̅⟩ + O[(v_L)^-1]
    .

When taking the limit v_L →∞, the O[(v_L)^-1] terms vanish, and the trial wavefunction |ϕ_N⟩ becomes

    (∏_α∈Θ, α∈ L P_α^†) |ϕ̅_N̅⟩ ,

which is a V-subspace pair condensate |ϕ̅_N̅⟩ (<ref>) plus a fully occupied L-subspace core ∏_α∈Θ, α∈ L P_α^† |0⟩ and an empty H-subspace. That is, we restrict (reduce) the valence space in the pairing channel to be the V-subspace.



Now we derive the energy gradient expression when using the trial wavefunction (<ref>). This can be done by taking the limit v_L →∞ on Eq. (<ref>) or Eq. (<ref>); here we skip the five-page long derivation and only show the results in Eqs. (<ref>)-(<ref>). For convenience, we introduce new symbols

    K_αβγμ   ≡    V_αβ̃γ̃μ ,  
    
    W_αβγμ   ≡    V_αγμβ + V_αγ̃μ̃β ,

and the HF mean field generated by the fully occupied L-subspace (the core)

    λ_αβ≡ϵ_αβ + ∑_γ∈ L V_αγγβ .

For mixing within the V-subspace,

    f_α∈ V, β∈ V = N̅^2 (v_α - v_β)/χ̅_N̅
    ×( [ (v_α + v_β) λ_αβ + v_α K_αβαα
    + v_β K_αβββ ] χ̅_N̅-1^[αβ]
    
    + ∑_γ∈Θ, γ∈ V^P_γ P_α, P_β v_γ K_αβγγ [ χ̅_N̅-1^[αβγ] - (N̅-1)^2 v_α v_βχ̅_N̅-2^[αβγ]]  
    
    + (N̅-1)^2 (v_α + v_β) ∑_γ∈Θ, γ∈ V^P_γ P_α, P_β (v_γ)^2 W_αβγγχ̅_N̅-2^[αβγ]) .

f_α∈ V, β∈ V means f_αβ with α∈ V, β∈ V. We see that Eq. (<ref>) is just Eq. (<ref>) applied to the V-subspace, while the empty H-subspace contributes nothing, and the full L-subspace changes ϵ_αβ to λ_αβ = ϵ_αβ + ∑_γ∈ L V_αγγβ (monopole correction). For mixing within the L-subspace,

    f_α∈ L, β∈ L = 0 .

This is correct because the L-subspace is full (a determinant), so the trial wavefunction (<ref>) and energy are invariant under mixing. For mixing within the H-subspace,

    f_α∈ H, β∈ H = 0 .

This is correct because the H-subspace is empty, so the trial wavefunction (<ref>) and energy are invariant under mixing. For mixing between the L and H subspaces,

    f_α∈ L,β∈ H = λ_αβ + ∑_γ∈Θ, γ∈ V W_αβγγ⟨ϕ̅_N̅ | n̂_γ | ϕ̅_N̅⟩
    
    = λ_αβ + ∑_γ∈Θ, γ∈ V W_αβγγ (1 - χ̅_N̅^[γ]/χ̅_N̅) .

The two forms are equivalent because of Eq. (<ref>). We see that Eq. (<ref>) is just the HF mean field (<ref>) corrected by the occupation numbers of the V-subspace ⟨ϕ̅_N̅ | n̂_γ | ϕ̅_N̅⟩ (monopole correction). For mixing between the L and V subspaces,

    f_α∈ L,β∈ V = 1/χ̅_N̅( λ_αβχ̅_N̅^[β]
    
    + N̅^2 ∑_γ∈Θ, γ∈ V^P_γ P_β v_γ ( v_γ W_αβγγ - v_β K_αβγγ ) χ̅_N̅-1^[βγ]) .

For mixing between the V and H subspaces,

    f_α∈ V,β∈ H = N̅^2 v_α/χ̅_N̅( v_α ( λ_αβ + K_αβαα ) χ̅_N̅-1^[α]
    
    + ∑_γ∈Θ, γ∈ V^P_γ P_α v_γ K_αβγγχ̅_N̅-1^[αγ]
    
    + (N̅-1)^2 v_α∑_γ∈Θ, γ∈ V^P_γ P_α (v_γ)^2 W_αβγγχ̅_N̅-2^[αγ]) .

Because f_αβ = - f_βα is skew-symmetric, f_β∈ H,α∈ L = - f_α∈ L,β∈ H, f_β∈ V,α∈ L = - f_α∈ L,β∈ V, and f_β∈ H,α∈ V = - f_α∈ V,β∈ H are known by Eqs. (<ref>), (<ref>), and (<ref>). So we have exhausted all possible cases of f_αβ.

In this section, we show that the gradient (<ref>) becomes the expressions (<ref>)-(<ref>) when reducing the valence space in the pairing channel. In the next section, these gradient expressions [Eqs. (<ref>) and (<ref>)-(<ref>)] are used in the iterative gradient minimizers to solve VDPC+CB.








§ ITERATIVE GRADIENT MINIMIZERS FOR VDPC+CB  


The family of iterative gradient minimizers <cit.> is frequently used to solve the variational principle (minimize average energy). It can easily deal with one or several constraints on the solution (subsidiary conditions) <cit.>. This is useful in many cases; for example, when we plot the potential energy surface as a function of the constrained multipole (quadrupole, octupole, …) deformation <cit.>, which is needed by the generator coordinator method to go beyond the mean field <cit.>. In this section, we briefly review three gradient minimizers from the family — the steepest descent <cit.>, the preconditioned gradient <cit.>, and ADAM <cit.> — in the content of VDPC+CB. Applying them to realistic examples will be in Sec. <ref>.








 §.§ Steepest Descent


VDPC+CB varies the canonical basis at fixed v_α. We rename θ in Eq. (<ref>) to be θ_αβ (the angle mixes |α⟩ and |β⟩), and introduce the vector θ⃗ that collects all the independent θ_αβ. So Eq. (<ref>) is rewritten as

    ∂E̅/∂θ_αβ |_θ⃗ = 0 = ∂ (⟨ϕ_N | H | ϕ_N ⟩)/∂θ_αβ |_θ⃗ = 0 = 4 f_αβ ,

and its vector form is

    ∇⃗E̅ = 4 f⃗ .

In one iteration, how to go from the current canonical basis {|α⟩} to the new canonical basis {|α'⟩}? The steepest descent minimizer <cit.> goes in the direction opposite to the gradient ∇⃗E̅. (The direction of the steepest descent for infinitesimal step size.) That is, we use the mixing angle

    θ⃗ = - η∇⃗E̅ = - 4 ηf⃗ ,

or for each component

    θ_αβ = - 4 η f_αβ .

The parameter η sets the step size; η may vary from iteration to iteration <cit.>. Known θ⃗, how to find the unitary transformation U_αα' = ⟨α | α' ⟩ from {|α⟩} to {|α'⟩}? In Eq. (<ref>), because f_αβ = - f_βα is a skew-symmetric matrix, we can also treat θ_αβ = - θ_βα as a skew-symmetric matrix. Then the matrix exponential of θ

    U = exp(θ)

is the desired unitary transformation U_αα'. (The matrix exponential of a skew-symmetric matrix is a unitary matrix.) If all the mixing angles are small, θ≈ 0, the transformation (<ref>) becomes

    U ≈ I + θ ,

which is consistent with Eq. (<ref>). This happens when the step size η≈ 0, or when near the energy minimum so the gradient ∇⃗E̅ = 4 f⃗≈ 0; then θ≈ 0 because of Eq. (<ref>) or Eq. (<ref>).



 §.§ Preconditioned Gradient  


The steepest decent minimizer performs poorly <cit.> when the Hessian matrix has a poor condition number: the partial derivative of E̅ changes slowly in one direction (one θ_αβ) but rapidly in another. In this case it is hard to select the step size η: η should be small enough to avoid overshooting the energy minimum and going uphill in the direction with rapid derivative change (large positive curvature), but then the progress is too small in directions with slow derivative change (small curvature). So many iterations are needed to converge to the energy minimum.

To solve this problem, the preconditioned gradient minimizer <cit.> uses different step sizes in different directions. Specifically, we introduce a precondition factor p_αβ for each direction θ_αβ; the vector p⃗ collects them. Then we replace Eq. (<ref>) by

    θ⃗ = - η∇⃗E̅/p⃗ = - 4 ηf⃗/p⃗ ,

or for each component [replacing Eq. (<ref>) by]

    θ_αβ = - 4 ηf_αβ/p_αβ .

If all p_αβ = 1, we go back to the steepest decent minimizer. Section <ref> will use p_αβ = |e_α - e_β| + 1 MeV, where e_α = d_α / 2 is the single-particle energy defined in Eq. (34) of Ref. <cit.>.





 §.§ ADAM


ADAM means adaptive moments of the gradient <cit.>. It uses the decaying sum of the historical gradient and the historical (element-wise) squared gradient to compute the mixing angle θ⃗. Specifically, ADAM keeps two vector variables m⃗ and s⃗, initialized to zero. In each iteration, we accumulate the gradient ∇⃗E̅ = 4 f⃗ into m⃗ and s⃗ according to

    m⃗   ←   γ_1 m⃗ + (1-γ_1) ∇⃗E̅ ,  
    s⃗   ←   γ_2 s⃗ + (1-γ_2) (∇⃗E̅)^2 .

In Eq. (<ref>), (∇⃗E̅)^2 means element-wise square. The two parameters γ_1 and γ_2 control the decaying rate for weight of the historical gradient and squared gradient. In this work we use γ_1 = 0.8 and γ_2 = 0.995 (slightly smaller than their default values, 0.9 and 0.999, suggested in the original ADAM paper <cit.>). For each component, Eqs. (<ref>) and (<ref>) read

    m_αβ   ←   γ_1 m_αβ + 4 (1-γ_1) f_αβ ,  
    
    s_αβ   ←   γ_2 s_αβ + 16 (1-γ_2) (f_αβ)^2 .


There is one problem about normalization <cit.>. At the n-th iteration, using Eq. (<ref>) recursively gives

    m⃗^(n)   =   γ_1 m⃗^(n-1) + (1-γ_1) (∇⃗E̅)^(n)
       =    (γ_1)^2 m⃗^(n-2) + γ_1 (1-γ_1) (∇⃗E̅)^(n-1)
           + (1-γ_1) (∇⃗E̅)^(n)
       =   …
       =    (γ_1)^n m⃗^(0) + (1-γ_1) [ (γ_1)^n-1 (∇⃗E̅)^(1)
           + … + γ_1 (∇⃗E̅)^(n-1) + (∇⃗E̅)^(n) ]  
       =    (1-γ_1) ∑_k = 0^n-1 (γ_1)^k (∇⃗E̅)^(n-k) ,

where m⃗^(k) and (∇⃗E̅)^(k) mean m⃗ and ∇⃗E̅ at the k-th iteration, and we have used the initial value m⃗^(0) = 0. So m⃗^(n) is a weighted sum of the historical (∇⃗E̅)^(k), k = 1, 2, …, n, and the weights sum to

    (1-γ_1) ∑_k = 0^n-1 (γ_1)^k = 1 - (γ_1)^n .

The normalized version is (n = 1, 2, 3, …)

    m̂⃗̂^(n) = m⃗^(n)/1 - (γ_1)^n .

Thus m̂⃗̂ is the weighted average of the historical gradients, because the weights sum to 1. Similarly, the normalized version of s⃗^(n) is

    ŝ⃗̂^(n) = s⃗^(n)/1 - (γ_2)^n .

ŝ⃗̂ is the weighted average of the historical squared gradients.

In each step, ADAM replaces Eq. (<ref>) by

    θ⃗ = - η m̂⃗̂/√(ŝ⃗̂+δ) ,

or for each component [replacing Eq. (<ref>) by]

    θ_αβ = - η m̂_αβ/√(ŝ_αβ + δ) ,

where δ is a small positive constant for numerical stabilization (to avoid dividing by zero when ŝ_αβ≈ 0). In this work we use δ = 10^-16 (MeV)^2.

Equation (<ref>) shows that ADAM uses different step sizes in different directions. If in some direction θ_αβ, one frequently overshoots the energy minimum because of a large positive curvature, the derivatives in successive iterations will frequently change sign so cancel each other, thus m̂_αβ will be small, which will damp the step size in this direction θ_αβ. The denominator √(ŝ_αβ) in Eq. (<ref>) also damps the step size in directions with large magnitude of gradient.

This section explains how to solve VDPC+CB by the family of iterative gradient minimizers. The next section explains how to solve VDPC+HFB.





§ VDPC+HFB  


VDPC+HFB varies v_α and the canonical basis together to minimize E̅. The new algorithm accomplishes this by combining VDPC+BCS (of Ref. <cit.>, vary v_α) and VDPC+CB (of Sec. <ref>, vary canonical basis). Specifically, we insert varying v_α (VDPC+BCS) into several places in the process of varying the canonical basis (VDPC+CB). It is unnecessary to vary v_α after every iteration of VDPC+CB; after every 20 iterations, for example, is enough. In the next section, we will apply the new VDPC+HFB algorithm to a semi-realistic example.









§ REALISTIC EXAMPLE  


The semi-realistic example of the rare-earth nucleus ^158_ 64Gd_94 has been used in our recent paper <cit.> for VDPC+BCS. In this section, we use the same example to demonstrate the VDPC+CB and VDPC+HFB algorithm. The purpose is to show the effectiveness of the algorithms under realistic interactions, not to accurately reproduce the experimental data. For simplicity, I consider only the neutron degree of freedom, governed by the anti-symmetrized two-body Hamiltonian

    H = ∑_αϵ_α a_α^† a_α + 1/4∑_αβγδ V_αβγδ a_α^† a_β^† a_γ a_δ .

The single-particle levels ϵ_α are eigenstates of the Nilsson model <cit.>. The Nilsson parameters are the same as in Ref. <cit.>; here I only repeat β = 0.349 (the experimental quadrupole deformation <cit.>). The neutron residual interaction V_αβγδ is the low-momentum NN interaction V_low-k <cit.> derived from the free-space N^3LO potential <cit.>.

Specifically, I use the code distributed by Hjorth-Jensen <cit.> to compute the two-body matrix elements of V_low-k in the spherical harmonic oscillator basis up to (including) the 𝒩 = 14 major shell, with the standard momentum cutoff 2.1 fm^-1. (𝒩 = 2n_r+l is the major-shell quantum number.) Then the Nilsson model is diagonalized in this spherical 𝒩≤ 14 basis (its dimension is D = 2 Ω = 1360, and Ω = 680 is the number of vacancies for Kramers pairs), the eigen energies are ϵ_α and the eigen wavefunctions transform the spherical two-body matrix elements into those on the Nilsson basis as used in the Hamiltonian (<ref>). The Hamiltonian (<ref>) of this work is exactly the same as that of Ref. <cit.> discussing VDPC+BCS.

This Hamiltonian (<ref>) has axial symmetry, so parity and angular-momentum projection onto the symmetry axis are good quantum numbers. We assume they are self-consistent symmetry; so when varying the canonical basis, the two single-particle basis states being mixed [see Eq. (<ref>)] must have the same parity and angular-momentum projection. [The trial wavefunction |ϕ_N⟩ (<ref>) has positive parity and zero angular-momentum projection.]

The VDPC+BCS code is written in Mathematica and runs in serial (no parallel computing, this code is from Ref. <cit.>). The VDPC+CB code is written in Matlab and runs in parallel. The VDPC+HFB code combines the two by calling Mathematica from Matlab. All the numerical calculations of this work were done on a laptop that has one quad-core CPU (Intel Core i7-4710MQ @ 2.5 GHz). All time costs plotted in the figures or given in the text are the actual time costs spent on this laptop. This work uses Matlab R2015a and Mathematica 10.2, to give the actual software version.





Now we discuss VDPC+CB, which varies the canonical basis at fixed v_α. How to fix v_α? We perform VDPC+BCS as in Ref. <cit.>. Specifically, we divide the full single-particle space (𝒩≤ 14) of dimension 1360 (Ω = 680) into L, V, H subspaces. This work uses an empty L subspace, a V subspace of dimension 700 that consists of all the levels having occupation number n_α > 1.63 × 10^-7, and an H subspace of dimension 660. Then we fix the 350 variational parameters v_α (350 = 700/2) by VDPC+BCS as in Ref. <cit.>; the computer time cost is about 43 seconds. Because the V subspace is smaller than the full space, there is a cutoff error, which is approximately 2 keV as can be read from Fig. 5 of Ref. <cit.>.

We distinguish between the energy cutoff error and the energy convergence error. It is a cutoff to divide the full single-particle space into L, V, H subspaces. (Because only the V subspace is active.) We define E(final) to be the converged energy (variational minimum) using a given cutoff (given L, V, H subspaces). E(exact) is the converged energy in the full space (no cutoff, the V subspace is the full space). The cutoff error is defined to be E(final) - E(exact). Using a given cutoff (given L, V, H subspaces), the energy iteratively converges to E(final). The convergence error is defined to be E(iter) - E(final), where E(iter) is the energy after iter of iterations. In the following, “error” means “energy convergence error”, unless we explicitly write “energy cutoff error”.

In passing, we list the parameters of VDPC+BCS used in this work. Skipping this paragraph does not affect reading this work; and understanding this paragraph needs Ref. <cit.>. This work performs VDPC+BCS similarly to Fig. 1 of Ref. <cit.>, but slightly changes the parameters of valence-space dimension. In step (iii), the dimension is increased to 240 (still do 10 iterations, more than enough for convergence). In step (v), the dimension is increased to 700 (still do 5 iterations, more than enough for convergence). Specifically, we sort the canonical single-particle basis states by their occupation numbers n_α estimated in step (iv), then from large n_α to small n_α we select 700 basis states, which form the V subspace. Because the V subspace is larger, the energy cutoff error (about 2 keV) is smaller than that in Fig. 1 of Ref. <cit.> (about 2.5 keV). The same VDPC+BCS parameters (240, 10 iterations; 700, 5 iterations) are used later in the VDPC+HFB calculations of this work. Note “iteration” in this paragraph means iterations in VDPC+BCS; everywhere else (outside this paragraph) “iteration” means iterations in VDPC+CB.

We tried to solve VDPC+CB by the three minimizers in Sec. <ref> — the steepest descent, the preconditioned gradient, and ADAM — and found ADAM is the most effective. (The preconditioned gradient minimizer (<ref>) uses p_αβ = |e_α - e_β| + 1 MeV. For α∈ V, e_α = d_α / 2 is the single-particle energy defined in Eq. (34) of Ref. <cit.>. For α∈ L and α∈ H, e_α is the HF single-particle energy [including the correction by the occupation numbers of the V-subspace, see Eq. (<ref>)].) The steepest descent and the preconditioned gradient are conventionally used to solve HF <cit.>; for VDPC+CB they are less effective, because there could be tiny partial derivative f_αβ≈ 0 if v_α≈ v_β, as shown by Eq. (<ref>). (Thus it is harder to select the step size as explained in the beginning of Sec. <ref>.) Figure <ref> shows the results by ADAM. Overall, the energy-error curve is linear on the log-scale plot, so energy converges exponentially with the number of iterations. The energy error drops to less than 1 keV at the 65th iteration, and less than 0.01 keV at the 100th iteration. The accumulated computer time cost increases linearly with the number of iterations, so each iteration costs the same time approximately (2.67 seconds in average). In Fig. <ref>, the ADAM parameters are γ_1 = 0.8, γ_2 = 0.995, δ = 10^-16 (MeV)^2, and a decaying step size — η = 0.05 / 2.5^k/100 at the k-th iteration. [See Eqs. (<ref>), (<ref>), and (<ref>) for definitions of these ADAM parameters.] The performance (speed of convergence) of the ADAM minimizer is not very sensitive to values of these ADAM parameters.





Now we discuss VDPC+HFB, which varies v_α and the canonical basis together. Specifically, we insert varying v_α (VDPC+BCS) into several places in the process of varying the canonical basis (VDPC+CB of Fig. <ref>). We ran VDPC+HFB twice, and the results are shown in Fig. <ref>. The first run inserts varying v_α at (immediately after) the 0th, 20th, 40th, 60th, 80th, 100th, 120th, … iteration (every 20 iterations); the second run inserts varying v_α at the 0th, 40th, 70th, 100th, 120th, … iteration. For each run, the accumulated computer time cost is a linear curve superimposed with sudden jumps where varying v_α is inserted. The slope of the linear curve (neglecting the jumps) is almost the same as that in Fig. <ref> (also 2.67 seconds per iteration in average). Each jump has a similar size, around 52 seconds, which includes the time cost of VDPC+BCS (about 43 seconds) and some overheads. For the first run, overall the energy-error curve is linear on the log-scale plot, so energy converges
exponentially with the number of iterations. The energy error drops to less than 1 keV at the 63th iteration, and less than 0.1 keV at the 82th iteration. For the second run, varying v_α is inserted less frequently, waiting until the energy curve flattens out (VDPC+CB converges). In summary, the first run is more efficient than the second run if look at the two energy curves of Fig. <ref> (energy error versus iterations); but in terms of energy error versus time cost, the two runs have similar efficiency when the energy error is bigger than 1 keV. Comparing Fig. <ref> and Fig. <ref>, the time cost of VDPC+HFB is roughly twice that of VDPC+CB, to achieve the same accuracy (energy error). This extra time cost of VDPC+HFB is mainly spent on the several times of varying v_α, which is not needed in VDPC+CB. Fine-tuning the VDPC+BCS parameters can decrease this extra time.



We benchmark the speed of the VDPC algorithms against that of HF by iterative gradient minimizers. We tried to solve HF by the three minimizers in Sec. <ref> — the steepest descent, the preconditioned gradient, and ADAM — and found the preconditioned gradient is the most efficient (energy converges the fastest). (The preconditioned gradient minimizer (<ref>) uses p_αβ = |e_α - e_β| + 1 MeV, where e_α is the HF single-particle energy.)  For the Hamiltonian (<ref>), HF by the preconditioned gradient typically (we tried different particle numbers) needs 20 ∼ 30 iterations to converge within 1 keV error, and needs 30 ∼ 50 iterations to converge within 0.01 keV error. So HF needs less number of iterations than VDPC+CB or VDPC+HFB. On the other hand, the computer time cost per iteration in HF (2.2 seconds in average) is similar to that in VDPC+CB or VDPC+HFB (2.67 seconds in average). In each iteration, most of time is spent on basis transformation by matrix multiplication (transform the needed Hamiltonian matrix elements into the new canonical basis of the next iteration), this part is the same for HF and for VDPC+CB. The gradient expression of VDPC+CB (<ref>) is more complicated than that of HF, but this part costs insignificant computer time. In summary, based on Fig. <ref> and Fig. <ref>, the total computer time cost of VDPC+CB and of VDPC+HFB are typically 2 ∼ 3 times and 4 ∼ 6 times that of HF, respectively, to achieve the same accuracy (energy error).

Future works can try to decrease the total computer time cost of the VDPC+HFB algorithm in two ways. First, we can decrease the time cost of VDPC+CB as suggested in Ref. <cit.>, for example, enable parallel computing (currently the code runs in serial). Second, we can find a better minimizer for VDPC+CB — needs less iterations to converge. For example, we can fine-tune the ADAM parameters or introduce additional ADAM parameters (third moment, fourth moment, …). However, it is hard to decrease the computer time cost per iteration, which is mainly spent on basis transformation by matrix multiplication (transform the needed Hamiltonian matrix elements into the new canonical basis of the next iteration).







In passing, if we do HF then VDPC+BCS, how good is the energy? (Sometimes, “HF then BCS” is used to approximate HFB.) For the Hamiltonian (<ref>), the (converged) HF energy is 1.225 MeV higher than the (converged) VDPC+HFB energy. VDPC+BCS (after HF) further lowers the energy by 0.326 MeV, so the (converged) energy of “HF then VDPC+BCS” is still 0.899 MeV higher than VDPC+HFB. (The pairing energy 0.326 MeV is small, because the density of the HF single-particle levels is accidentally small near the Fermi surface.) Therefore, “HF then VDPC+BCS” is not good enough, and VDPC+HFB is necessary.











§ CONCLUSIONS  


Recently Ref. <cit.> proposed a scheme that applies the variational principle directly to the coherent pair condensate in the BCS case (VDPC+BCS). This work extends the scheme to the HFB case (VDPC+HFB). The result is equivalent to that of the so-called variation after particle-number projection in the HFB case (VAP+HFB), but now the particle number is always conserved and the time-consuming projection is avoided. The HFB theory is frequently criticized for breaking the exact particle number. Meanwhile, VAP+HFB by the numerical gauge-angle integration may not be very easy, and in the literature there are far fewer realistic applications of VAP+HFB than those of HFB without projection. We hope the new VDPC+HFB algorithm could become a common practice because of its simplicity.

Specifically, this work derives the analytical expression for the gradient of the average energy with respect to changes of the canonical basis. The VDPC+CB algorithm supplies this gradient expression to the family of gradient minimizers to iteratively minimize energy. In practice, we find the so-called ADAM minimizer, borrowed from the machine-learning field, is very effective. The VDPC+HFB algorithm combines VDPC+BCS and VDPC+CB, by inserting VDPC+BCS (vary v_α) into several places in the process of VDPC+CB (vary the canonical basis). The new algorithms are demonstrated in a semi-realistic example using the realistic V_low-k interaction and large model spaces (up to 15 harmonic-oscillator major shells). They easily run on a laptop, and practically the computer time cost to solve VDPC+HFB (time to solve VDPC+CB) is about 4 ∼ 6 times (2 ∼ 3 times) that to solve HF by the iterative gradient minimizers. Future works can further optimize the code for less time cost, as discussed in Sec. <ref>.

The parameters in nuclear mean-field interactions <cit.> are usually fitted (to experimental data) in the theory of HF, HFB, or approximate VAP+HFB (for example, the Lipkin-Nogami prescription followed by particle-number projection <cit.>). These fitting parameters should be fine-tuned to generate optimum interactions for VDPC+HFB. In addition, higher-order correlations beyond VDPC can be included, for example, by the generalized-seniority truncation of the shell model <cit.>.




§ ACKNOWLEDGEMENTS


Support is acknowledged from the National Natural Science
Foundation of China No. 11405109.

















































































































50

Jia_2019  L. Y. Jia, Phys. Rev. C 99, 014302 (2019).

Bender_2003  M. Bender, P.-H. Heenen, and P.-G. Reinhard, Rev. Mod. Phys. 75, 121 (2003).

Bohr_1958 A. Bohr, B. R. Mottelson, and D. Pines, Phys. Rev.
110, 936 (1958).



Bohr_book  A. Bohr and B. Mottelson, Nuclear Structure
(Benjamin, New York, 1975).

Ring_book  P. Ring and P. Schuck, The Nuclear Many-Body
Problem (Springer-Verlag, Berlin, 1980).

Belyaev_1959 S. T. Belyaev, K. Dan. Vidensk. Selsk. Mat. Fys.
Medd. 31(11), 641 (1959).

Dietrich_1964  K. Dietrich, H. J. Mang, and J. H. Pradal, Phys. Rev. 135, B22
(1964).



Dukelsky_2000  J. Dukelsky and G. Sierra, Phys. Rev. B 61, 12302 (2000).

Dussel_2007  G. G. Dussel, S. Pittel, J. Dukelsky, and P. Sarriguren, Phys. Rev. C 76, 011302(R) (2007).

Sandulescu_2008  N. Sandulescu and G. F. Bertsch, Phys. Rev. C 78, 064318 (2008).

Sheikh_2000  J. Sheikh and P. Ring, Nucl. Phys. A665, 71 (2000).

Anguiano_2001  M. Anguiano, J. Egido, and L. M. Robledo, Nucl. Phys. A696, 467 (2001).

Anguiano_2002  M. Anguiano, J. L. Egido, and L. M. Robledo, Phys. Lett. B545, 62 (2002).

Stoitsov_2007  M. V. Stoitsov, J. Dobaczewski, R. Kirchner, W. Nazarewicz, and J. Terasaki, Phys. Rev. C 76, 014308 (2007).

Hupin_2012  G. Hupin and D. Lacroix, Phys. Rev. C 86, 024309 (2012).

Wang_2014  X. B. Wang, J. Dobaczewski, M. Kortelainen, L. F. Yu, and M. V. Stoitsov, Phys. Rev. C 90, 014312 (2014).

Reinhard_1982  P.-G. Reinhard and R. Y. Cusson, Nucl. Phys. A378, 418 (1982).

Robledo_2011  L. M. Robledo and G. F. Bertsch, Phys. Rev. C 84, 014312 (2011).

Adam_2014  D. Kingma and J. Ba, arXiv:1412.6980 [cs.LG] (2014).

Goodfellow_book  I. Goodfellow, Y. Bengio, and A. Courville, Deep Learning (MIT Press, 2016), <http://www.deeplearningbook.org/>.

Bogner_2003  S. Bogner, T. T. S. Kuo, and A. Schwenk, Phys. Rep. 386, 1
(2003).





Hupin_2011_2  G. Hupin and D. Lacroix, Phys. Rev. C 83, 024317 (2011).

Hupin_2011  G. Hupin, D. Lacroix, and M. Bender, Phys. Rev. C 84, 014309 (2011).

Lipkin_1960  H. J. Lipkin, Ann. Phys. 9, 272 (1960).

Nogami_1964  Y. Nogami, Phys. Rev. 134, B313 (1964).

Pradhan_1973  H. C. Pradhan, Y. Nogami, and J. Law, Nucl. Phys. A201, 357
(1973).





Goodman_1979  A. L. Goodman, Adv. Nucl. Phys. 11, 263 (1979).






Jia_2017  L. Y. Jia, Phys. Rev. C 96, 034313 (2017).








Warda_2002  M. Warda, J. L. Egido, L. M. Robledo, and K. Pomorski, Phys. Rev. C 66, 014310 (2002).

Delaroche_2010  J.-P. Delaroche, M. Girod, J. Libert, H. Goutte, S. Hilaire, S. Peru, N. Pillet, and G. F. Bertsch, Phys. Rev. C 81, 014303 (2010).



Nilsson_1955  S. G. Nilsson, Mat. Fys. Medd. Dan. Vid. Selsk. 29, 16 (1955).


NNDC  Data retrieved from the National Nuclear Data Center (Brookhaven National Laboratory), <http://www.nndc.bnl.gov/>.


Entem_2003  D. R. Entem and R. Machleidt, Phys. Rev. C 68, 041001(R)
(2003).


Morten_code  <https://github.com/ManyBodyPhysics/ManybodyCodes/CENS>













Samyn_2004  M. Samyn, S. Goriely, M. Bender, and J. M. Pearson, Phys. Rev. C 70, 044309 (2004).

Talmi_book  Igal Talmi, Simple Models of Complex Nuclei: The Shell Model and Interacting Boson Model (Harwood Academic, Chur, Switzerland, 1993).

Allaart_1988 K. Allaart, E. Boeker, G. Bonsignori, M. Savoia, Y.K. Gambhir,
Phys. Rep. 169, 209 (1988).

Zhao_2014  Y.M. Zhao, A. Arima, Phys. Rep. 545, 1 (2014).

Jia_2015  L. Y. Jia, J. Phys. G: Nucl. Part. Phys. 42, 115105 (2015).

Jia_2016_ph  L. Y. Jia, Phys. Rev. C 93, 064307 (2016).

Jia_2016_Sn  L. Y. Jia and C. Qi, Phys. Rev. C 94, 044312 (2016).


Qi_2016  C. Qi, L. Y. Jia, and G. J. Fu, Phys. Rev. C 94, 014312 (2016).

Jia_arXiv_2018  L. Y. Jia, arXiv:1808.03729 [nucl-th] (2018).





















	
















































































































































































































































































































































































































 













