{'fraction_non_alphanumeric': 0.03993304638928742, 'fraction_numerical': 0.01250256200040992, 'mean_word_length': 4.2785791561485755, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Pairing correlation involving the continuum states

The Hartree-Fock-Bogoliubov equation for the ground states of even-even atomic nuclei is solved using the canonical representation in the coordinate space for zero range interactions like the Skyrme force. The gradient method is improved for faster convergence to the solutions under constraint of orthogonality between canonical orbitals. Necessity of the cut-off of the pairing interaction is shown even when the number of the canonical orbitals are restricted. A repulsive dependence of the interaction on the pairing density is introduced as an implementation of the cut-off which leaves the HFB super matrix state-independent.


Introduction
In neutron-rich nuclei, the pairing correlation significantly involves the continuum single-particle states. This makes the HF+BCS approximation inadequate due to the unlocalization of the neutron density distribution and demands one to solve the Hartree-Fock-Bogoliubov (HFB) equation without approximations. The solution in the coordinate space was first formulated in Ref. 1) using the quasi-particle states and performed for spherically symmetric states. However, its application to deformed states is difficult because there are quite a large number of quasiparticle states even for a moderate size of the normalization box (i.e., the cavity to confine the nucleons to discretize the single-particle states). Every HFB solution has an equivalent expression of BCS variational form. The single-particle states to construct the BCS type wavefunctions are called the HFB canonical basis or sometimes the natural orbitals. This expression was used to solve the HFB for spherical states originally in Ref. 2) .
Although spherical solutions can be obtained easily with present computers (for zero-range forces), deformed solutions are still difficult to obtain. The two-basis method 3,4,5) is the only one implemented so far for neutron-rich nuclei. Some recent developments like Ref. 6) are also in progress.
We have applied the canonical-basis method to deformed states in a three dimensional cubic mesh representation with density dependent delta interactions. 7) It has turned out to be a very efficient alternative approach to obtain the solutions. The origin of its effectiveness is that the number of necessary single-particle basis states to describe the ground state of a nucleus is proportional to the volume of the nucleus in the canonical-basis method while it commensurates with the volume of the normalization box in conventional methods. The difference of the number of the basis states amounts to a factor of 10 1 -10 3 .
In this paper we discuss the canonical-basis formulation of the HFB, the method to obtain the canonical-basis solutions, faster gradient-method paths than a naive imaginary-time evolution, the necessity of the cut-off of the pairing interaction, and an implementation of the cutoff in terms of an interaction dependent on the pairing density.

HFB in the canonical representation
To begin with, let us formulate the HF and the HFB in the coordinate-space representation in order to elucidate a difficulty of the HFB and to suggest its possible solution in terms of the canonical-basis representation. For the sake of simplicity, we consider only one kind of nucleons and designate the number of nucleons by N in Eqs. (1)- (14), which are in this section. The z-component of the spin of a nucleon is represented by s (= ± 1 2 ).
In the HF, one should minimize Ψ|H|Ψ for single Slaterdeterminant states, by varying {ψ i } i=1,···,N under orthonormality conditions ψ i |ψ j = δ ij . The operator a † i creates a nucleon with a wavefunction ψ i ( r, s). The distribution function of the density of the nucleons is related to the wavefunctions as In the HFB, the solution takes the following form, where b i is the annihilation operator of a negative-energy Bogoliubov quasi-particle with amplitudes φ i ( r, s) for presence and ϕ i ( r, s) for absence. I is the number of the basis states of the employed representation. For a threedimensional Cartesian mesh (3D-mesh) representation, 8) it is the number of the mesh points (times four when spin-orbit interactions are included) and typically 10 4 - and a constraint on the expectation value of the number of nucleons, N = s d ra( r, s)a † ( r, s).
The nucleon density for state (4) is expressed as The essential difference between the HF and the HFB is that one has to consider only N ∼ 10 2 wavefunctions in the former while one has to treat explicitly as many single-particle wavefunctions as the number of the basis in the latter.
Owing to the Bloch-Messiah theorem, 9) the state (4) can be expressed (for the ground states of even-even nuclei) as, where a † i and a † ı create a nucleon with wavefunction ψ i ( r, s) and ψī( r, s), respectively, which are called as the canonical basis 9) or the natural orbitals 2) of the HFB vacuum |Ψ . One must use K = 1 2 I for the exact equivalence between Eqs. (4) and (10) in general cases. When |Ψ is a time-reversal invariant state, which we assume in this paper, ψ i and ψī are the time-reversal state of each other. In this case, only one wavefunction of each time reversal pair should be counted as independent variables of the variational procedure.
To obtain solutions in the canonical-basis framework, one should vary {ψ i , u i , v i } i=1,···,K under three kinds of constraints, i.e, the orthonormality conditions, fixed expectation value of the number of nucleons, and the normalization of the u-v factors The nucleon density is expressed as Reinhard et al. regarded that the advantage of the representation (10) over (4) is that one has to consider only a single set of wavefunctions {ψ i } i=1,···,I unlike a double set {φ i , ϕ i } i=±1,···,±I/2 . 2) However, we expect much greater benefit from the canonical-basis representation. Namely, i may be truncated as i ≤ K = O(N ) ≪ 1 2 I to a very good approximation. It is because ψ i appearing on the right-hand side of Eq. (14) must be a localized function as ρ( r) on the left-hand side, while the orthogonality (12) does not allow many low-energy wavefunctions to exist in the vicinity of the nucleus. For 3D-mesh representations, I is proportional to the volume of the cavity while K is proportional to the volume of the nucleus. The latter is 10 1 -10 3 times as small as the former.
Incidentally, the situation is quite different in the quasiparticle formalism. On the one hand, the localization of the density demands only the localization of ϕ i through Eq. (9) while φ i does not have to be localized in general. On the other hand, the orthogonality condition (6) involves both ϕ i and φ i . This discrepancy enables many quasiparticle states having similar ϕ i to be orthogonal to each other by differing their φ i .

Mean fields for zero-range interactions
Let us present the effective Hamiltonian employed in this paper. We adopt a density-dependent zero-range interaction. Zero-rangeness makes the mean-field potentials local, which is an essential advantage for coordinate-space solutions. On the other hand, the omission of momentum dependences are merely for the sake of simplicity and there will not be essential differences in the formulation if we use the full-form Skyrme force. 9) Our force is expressed using the parameterization of the Skyrme force as,v ( r 1 , s 1 ; r 2 , s 2 ) where r i and s i are the position vector and the spin state of the two interacting nucleons, i = 1, 2. Dependence on the isospin state is redundant because of the zerorangeness. We adopt t 0 = −983.4 MeV fm 3 , t 3 = 13106 MeV fm 3+3α , and α = 0.98 when the force is used to construct the mean-field (HF) potential. 10) We use different strengths to make the pairing potential. We express the pairing force in the parameterization of Ref. 11) aŝ P σ is the exchange operator of spin variables. 1 2 (1 − P σ ) is a projector into spin-singlet states so that the interaction acts only between like nucleons. We use ρ c = 0.32 fm −3 (to roughly vanish the volume-changing effect 11) ) and v p = −440 MeV fm 3 .
When one considers both protons and neutrons, the state of the nucleus is expressed as a product of two BCS type states (10) for the protons and the neutrons: where q distinguishes between protons (p) and neutrons (n). a † i,p creates a proton having a wavefunction ψ i,p ( r, s) while a † i,n creates a neutron with a wavefunction ψ i,n ( r, s). The product form is due to the pairing interaction (16) acting only between like nucleons.
For the sake of simplicity, we treat N =Z nuclei without Coulomb interaction in this paper. In this case, the wavefunctions are the same between protons and neutrons, i.e., ψ i,p ( r, s) = ψ i,n ( r, s), ψī ,p ( r, s) = ψī ,n ( r, s). Moreover, because the potentials are independent of the spin, the wavefunctions ψ i ( r, s) can be factorized into a product of a spin wavefunction and a real function of the position, which we write ψ i ( r) in the following. It holds ψ i,p ( r) = ψī ,p ( r) = ψ i,n ( r) = ψī ,n ( r). (15) and (16), the total energy for state (17) is written as,

With the interactions
where m is the average of the proton and the neutron masses divided by 1−1/A for the correction of the centerof-mass motion. H( r) is called the Hamiltonian density while function of position r in the right-hand side are where g = 4, which is a factor to account for the situation that a wavefunction ψ i takes care of four nucleons for the spin and isospin degeneracy.
The mean-field potential V and the pairing potentialṼ are defined as in which V ,Ṽ , H, ρ, andρ are local functions of r while t 0 , t 3 , α, v p , and ρ c are constants.
The mean-field and the pairing Hamiltonians are The quasiparticle states are the eigenvectors of the socalled HFB super matrix composed of h andh: This is just an eigenvalue problem of a hermitian (because of the time-reversal symmetry) matrix. The canonical orbitals are also determined by h andh but in a more complex way as described in the next section.

Gradient method for canonical-basis HFB
In this section we describe a procedure to obtain the canonical-basis solution of the HFB equation directory, not by way of quasi-particle states. Instead of minimizing E = Ψ|H|Ψ with |Ψ given by Eq. (17) under constraints of Eqs. (12) and (13), one may introduce a Routhian R, and minimize it without constraints. ǫ F is probably the most familiar Lagrange multiplier, whose physical meaning is the Fermi level. In the definition (25), K 2 Lagrange multipliers λ ij obeying hermiticity, are introduced instead of 1 2 K(K + 1) independent multipliers. This hermitization of λ is adopted in order to make R real so that two conditions, δR/δψ i = 0 and δR/δψ * i = 0, become equivalent and thus one has to consider only one of them. Note that δ ij is subtracted from ψ i |ψ j , in contrast to Ref. 2) , which is in order to treat λ ij not as constants like ǫ F but as functionals of the wavefunctions.
Using notations, the stationary conditions of R result in two kinds of equations. One is ∂R/∂v i = 0, which concerns the occupation amplitudes v i and is fulfilled by The minimum of R for the variations of u i and v i is obtained (when ∆ i ≥ 0) by taking the minus sign for the double sign in the right-hand side of Eq. (28) and choosing the same sign as v i for u i = ± 1 − v 2 i . The stationary condition for a wavefunction ψ i ( r, s) gives the following equation: where One can regard H i as a state-dependent single-particle Hamiltonian. This dependence on states makes the orthogonalization conditions essential to the method. For HF, the orthogonalization conditions are easily fulfilled because ψ i are eigenstates of the same hermite operator h and thus are orthogonal between themselves at the solution : ψ j |ψ i · (ǫ j − ǫ i ) = 0. The orthogonalization procedure is needed only because states satisfying the orthogonality are unstable for decaying into Pauliforbidden configurations. On the other hand, for the canonical-basis HFB method, the orthogonalization is essential because the single-particle Hamiltonians H i differs from state to state. Therefore, the determination of the explicit functional form of λ ij is the most important part of the method. Reinhard et al. have proposed Let us reason on which grounds the above definition can be deduced. Understanding these grounds is indispensable in order to modify the definition later for faster convergences. From the requirement that Eq. (29) must hold at the solution (where ψ i |ψ j = δ ij ), one can deduce, Eqs. (31) and (32) are equivalent at the solution because λ ij is defined to be hermite by Eq. (26). Since this hermiticity must hold at any points to ensure the equality between the number of constraints and the number of independent multipliers, one should not adopt Eq. (32) but Eq. (31).
One can utilize the gradient method to obtain the HFB solutions in the canonical-basis formalism. The most naive implementation of the gradient method agrees with the imaginary-time evolution method in its first order approximation of the size of the imaginary-time step ∆τ : We have developed a 3D-mesh canonical-basis HFB program from scratch according to the above formulation.
We take an example of our calculations using the program for the ground state of 40 Ca. The wavefunctions are expressed with 39 × 39 × 39 mesh points with mesh spacing of a=0.8 fm. We employed the 17-point finitedifference approximation to the Laplacian. Note that the requirement of precision is higher for HFB than for HF because one has to treat larger momentum components than the Fermi momentum in HFB. The vanishing boundary conditions are imposed on the boundary (the 0th and the 40th mesh points) and the wavefunctions are anti-symmetrically reflected in the boundary to apply the finite-difference formula. We considered K = 20 canonical basis, which can contain Kg = 80 (=2 × A, A = 40) nucleons.
In Fig. 1, the error of the second equality of Eqs. (29) neglecting the error of orthogonality, i.e., max i=1,···,K |H i ψ i − j λ ij ψ j | is plotted with a solid line versus the number of evolution steps. The corresponding quantity for HF, max i=1,···,A/4 |hψ i − ψ i |h|ψ i ψ i | , is also plotted with a dash line. The figure demonstrates that one can indeed obtain HFB solutions with the natural-orbital HFB method in the 3D-mesh representation. We obtained similar convergence curves for the error of the orthogonality and for the inconsistency between the potential and the densities.
The speed of the convergence is, however, about ten times as slow as the HF case. This is the subject of the next section.  We adopted an additional procedure which is not indispensable to obtain the solutions but effective to make the convergence more robust and somewhat quicker: After every 25th gradient-method steps, we orthogonalize {ψ i } with the Gram-Schmidt algorithm in the ascending order of ǫ i , defined in Eqs. (27), and then diagonalize the super matrix of the HFB Hamiltonian (24) by expanding the quasi-particle wavefunctions (φ i , ϕ i ) (1 ≤ i ≤ K and their negative-energy partners) in a 2K-dimensional basis {ψ i } ⊕ {ψ i }, and finally transform the resulting quasi-particle wavefunctions to canonical orbitals and occupation amplitudes {ψ i , v i } to renew them.
Incidentally, the period of 25 steps may be too frequent because the effect seems to saturate at periods around 100. We adopt the period of 25 in most calculations, however, because the increase in the computation time is only a several percent of the total time with this period.

Acceleration of the gradient method
We show the origin of the slow convergence and present a solution of the difficulty in this section.
Steepest-descent paths, which the gradient method draws, depend on the choice of the independent variables. For example, Eq. (33) is obtained when one uses (Re ψ i , Im ψ i ) as independent variables to define the gradient vector and then express it in coordinates (ψ i , ψ * i ). If one uses scaletransformed wavefunctions χ i ≡ α −1/2 i ψ i , where α i is a scaling factor, a gradient step becomes, which is equivalent to The change from Eq. (33) to Eq. (36) is equivalent to multiplying α i to the single-particle Hamiltonian H i in Eqs. (29).
When one parameterizes the scaling factor as α i = v −2ν i , the modified single-particle Hamiltonian becomes When ν = 0 (i.e., α i = 1), to which the imaginary-time evolution (33) corresponds, the single-particle Hamiltonian α i H i (=H i ) can be very small for canonical orbitals whose ǫ i is much higher than the Fermi level (i.e., ǫ i − ǫ F ≫ ∆ i ). This smallness makes the changes of such orbitals very slow. On the other hand, for ν = 1, the potential can be very deep for such high-lying orbitals due to the factor u i /v i in front ofh. In this case, however, the gradient step may be numerically dangerous. We usually use 1 2 ≤ ν < 1, which provides a fast and numerically stable method of solution.
When one introduces the acceleration factors (i.e., α i > 1), the multipliers λ ij should be modified from Eq. (31) by the following reason. In the computation of a gradient vector given by the first of Eqs. (29), the last term takes much more computing time than the first two terms due to δλ jk /δψ * i . One can forget the last term if the orthogonality relations (12) are fulfilled along the path of the steepest descent. Let's suppose that the relations are satisfied before a gradient-method step is taken and require that they are conserved to the first order in ∆τ after the step, i.e., with Substituting Eq. (39) into Eq. (38) and requiring the hermiticity (26) result in This form of λ ij fulfills the requirement that it should agree with the expression (32) at the solution as the naive form of Eq. (31) does. Two forms differ, however, before reaching the solution if one assumes α i = α j in general. Therefore, in order to conserve the orthogonality to vanish the last term in Eq. (29) one must not use Eq. (31) but Eq. (40). We have indeed suffered from large errors of orthogonality by using the naive form (31). On the other hand, by using the correct form (40), we have observed not only that the error does not grow during the evolution but also that the error decreases without performing explicit orthogonalization procedures periodically during the evolution. This decrease should originate in the second order terms in ∆τ in Eq. (38), whose effects we did not consider.  We compare the results of calculations between ν = 0 and ν = 1 2 in Fig. 2. The wavefunctions are expressed with 19×19×19 mesh points with mesh spacing of a=1.0 fm. We considered K = 16 canonical basis, which can contain Kg = 64 (=2×A) nucleons. The figure shows the convergence history to a HFB solution. Four quantities are plotted as functions of the number of gradient steps. They are, from the top to the bottom, the total energy E, the pairing gap∆ (averaged with weight u i v i ), the size of quadrupole deformation β, and the size of triaxiality of deformation γ. The last two quantities are determined from the mass quadrupole moments. 12) The dot curves were obtained without accelerations, i.e., with α i = 1 or ν = 0 in Eq. (37), while the solid ones were obtained with the acceleration method with α i = 1/v i or ν = 1 2 . One can see that the convergences of these quantities become by far faster by using the acceleration method. This result demonstrates that canonical-basis HFB can be solved without very heavy numerical computations.

Cut-off of the pairing interaction for canonical-basis HFB
Finally let us discuss on the cut-off schemes of the paring interaction in relation to the canonical-basis HFB method.
Delta function forces without cut-off leads to a divergence of the strength of the pairing correlation. 13) In order to circumvent the divergence, in the conventional method of solution, one usually takes only quasiparticles whose excitation energy ǫ qp is lower than some cut-off energy parameter ǫ cut to construct the ground state. 1) Namely Eq. (4) is modified to where θ is the step function.
In the canonical-basis method, the restriction on the number of canonical orbitals may prevent the divergence without introducing explicitly a cut-off energy. We have examined this idea by performing numerical calculations for various situations. Then, we noticed that observables sometimes jumps suddenly in the course of long-time evolution.  An example is shown in Fig. 3. The calculation was done for a nucleus 32 S on a 19×19×19 cubic mesh with a mesh spacing a=1 fm. We considered K = 20 canonical orbitals, to which we gave harmonic oscillator wavefunctions at the beginning. An 11-point formula was employed for the Laplacian. The acceleration parameter in Eq. (37) is taken as ν = 0.7. At first we suspected that these jumps were due to the acceleration method. However, with smaller ν, we still observed jumps; only they come later.
We investigated the origin of the jumps and found that each sudden change was due to a shrinkage of a highlying (i.e., having large ǫ i ) canonical orbital to a mesh point. This shrinkage can decrease the total energy of the nucleus for the following reason: The contribution of a canonical orbital to the kinetic energy density is proportional to its BCS occupation probability v 2 i , while pairing density is proportional to u i v i . Because v i ≪ 1 and u i . = 1 for high-lying orbitals, it holds that u i v i ≫ v 2 i . Therefore the increase in the kinetic energy due to the shrinkage is easily compensated by the gain of the pairing correlation energy at the shrunken point.
The observed jumps indicate that, in most cases (or maybe all the cases), there are no potential barriers between such physically meaningless solutions and the physically reasonable one.
Incidentally, we confirmed that the lack of potential barriers was not due to the discrete approximation of the Laplacian in the kinetic energy term: The approximation was based on the Lagrange polynomial interpolation, which tend to underestimate the expectation value for high momentum components. One might suspect that more accurate treatments of the kinetic energy could restore a barrier between the physical and unphysical solutions. However, we observed similar jumps by using the Fourier transformation with periodic boundary conditions, 14) which gives the exact result up to π/a.
We have decided that it is necessary to introduce cut-off for the reliability of the method.
As the cut-off scheme, we first employed cut-off factors which are dependent on orbitals. In the BCS approximation, a smooth cut-off method 8) is often utilized, in which the interaction is modified aŝ where f (ǫ) is a function of single-particle energy ǫ. The function takes on ∼ 1 well below a chosen cut-off energy and smoothly becomes zero above it. In analogy to the smooth cut-off method, we modified the pairing density asρ with We made the dependence to be on the kinetic energy (∝ k 2 i ), not on the mean-field energy ǫ i as in Eq. (43), to avoid a highly complicated expression of the gradient for the latter case. (In BCS, such complications are just neglected.) The result was successful to prevent the shrinkages to points. However, this cut-off scheme has an disadvantage that the HFB super matrix in Eq. (24) cannot be defined. It follows that we do not have well-defined quasiparticle states and cannot utilize them to express the HFB ground state. This drawback is rather serious because quasiparticles are useful to improve the precision of HFB solutions obtained by the canonical-basis method and to accelerate the convergence further, as described in the fourth section.
As an alternative method, we have introduced a repulsive pairing-density dependence to the pairing force, Eq. (16), in addition to the usual density dependence: With forces (15) and (46), the expectation value of the energy for N = Z systems is expressed as a space integral of a Hamiltonian density: The mean-field potential V remains the same as Eq. (20) while the pairing potentialṼ has an additional term: With this new type of force, the shrinkage problem is completely removed.  In Fig. 4, we show the dependence of the pairing gap (averaged with weight u i v i ) on the parameterρ c , which controls the pairing density dependence. The values are taken after 5,000 gradient steps. The set up of the calculations are the same as in Fig. 3 except forρ c . One can see that the pairing gap has reasonable values withρ c ≤ 1 fm −3 . This is a good news because the pairing-densitydependent term, −(ρ/ρ c ) 2 , can be small for the values ofρ which one finds in physical solutions, in contrast to the usual density-dependent term, −ρ/ρ c , which cancels roughly 50% of the density-independent term inside nucleus. This situation is illustrated in Fig. 5.   We show an example of the time evolution of ǫ i in Fig. 6. The set up of the calculation is the same as in Figs. 3 and 4 except that a 7-point approximation is used for the Laplacian to favor the emergence of unphysical solutions and make this calculation a very severe test of the cutoff scheme. The abscissa is the number of the gradient steps while the ordinate is the expectation value ǫ i of the mean-field Hamiltonian for all the K canonical orbitals. The pairing density dependence parameterρ c is set at the standard value of 0.3 fm −3 before step 3,000 (interval I) and after step 6,000 (interval III). Between steps 3,000 and 6,000 (interval II)ρ c is increased temporally to 3 fm −3 , which results in a too weak dependence to prevent the divergence. One can see that an unphysical solution emerges without the pairing-density dependence term (in interval II) and that the unphysical solution is suppressed (in interval I) and is quickly restored to a physical one (in interval III) with the presence of the term.
Before ending the section, let us mention that our next subject concerning the cut-off is whether repulsive momentum dependences of the Skyrme force help to make the cut-off scheme more natural because it has been known that such momentum dependences make well-developed plateau before the divergence sets in. 13) Figure 6 Fig. 6. Evolution history of ǫ i , i.e., the expectation values of the mean-field Hamiltonian for canonical orbitals. The pairing-density dependent term is switched on in intervals I and III (ρ c = 0.3 fm −3 ). It is suppressed to be very weak in interval II (ρ c = 3 fm −3 ).

Conclusions
We have developed a method to obtain canonical-basis HFB solutions in a coordinate-space three-dimensional (3D) Cartesian mesh representation. The features of our method are summarized as follows.
i) It is not for spherical but for deformed nuclei and thus it can treat both deformation and continuum pairing simultaneously.
ii) It is not based on the oscillator-basis expansion but described in the 3D Cartesian mesh representation and thus can treat e.g. deformed halo-like orbitals.
iii) There is a strong reason to believe that the necessary number of canonical orbitals is much smaller than the number of single-particle basis.
iv) In order to perform variations under constraint of orthogonality between the wavefunctions of the canonical orbitals, Lagrange multipliers were introduced as functionals of the wavefunctions in Ref. 2) . We have clarified the necessary conditions for the form of these Lagrange multiplier functionals. We have modified the functionals appropriately so that they conserve the orthogonality during the course of the accelerated evolutions explained in the next item.
v) We have found that the convergence to HFB solutions is very slow when one employs a naive gradient method. On the other hand, the convergence is quite rapid for Hartree-Fock solutions, which neglects the pairing correlation. We have investigated the origin of the slow convergence and found that the time scale of the gradient evolution is different from one canonical orbital to another depending on their BCS occupation amplitudes v i . The difference ranges over many orders of magnitude.
We have introduced an orbital-dependent acceleration method of the gradient evolutions and could overcome the difficulty of the slow convergence.
vi) We have examined the effects of the cut-off of the pairing interaction. The 3D mesh mean-field methods have a practical use only for zero-range interactions like the Skyrme force presently. One needs a cut-off for zerorange pairing forces, without which the pairing correlation energy diverges to −∞. 13) We have found that zerorange forces need cut-off even when the number of canonical orbitals are finite. The divergence occurs through shrinking of high-energy canonical orbital(s) to (a) mesh point(s).
vii) To suppress the divergence, we have introduced a pairing-density dependent interaction as a better choice than orbital-dependent cut-off factors.
We believe that, by choosing a faster gradient path with the acceleration method and adopting the cut-off scheme in terms of the pairing-density dependence, the canonicalbasis HFB method is now fully understood and has the potential to become the standard method to treat neutronrich nuclei in HFB. The details of this work will be published soon. 15)