New algorithm in the variation after projection calculations for non-yrast nuclear states
Jia-Qi Wang Zao-Chun Gao Ying-Jun Ma Y. S. Chen







China Institute of Atomic Energy, Beijing 102413,
People's Republic of China
College of Physics, Jilin University, Changchun 130012, China


zcgao@ciae.ac.cn
China Institute of Atomic Energy, Beijing 102413,
People's Republic of China


College of Physics, Jilin University, Changchun 130012, China


China Institute of Atomic Energy, Beijing 102413,
People's Republic of China

             


We present a novel and simple algorithm in the variation after projection (VAP) approach for the non-yrast nuclear states. It is for the first time that the yrast state and non-yrast states can be varied on the same footing. The orthogonality among the calculated states is automatically fulfilled by solving the Hill-Wheeler equation. This avoids the complexity of the frequently used Gram-Schmidt orthogonalization, as adopted by the excited VAMPIR method. Thanks to Cauchy's interlacing theorem in matrix theory, the sum of the calculated lowest projected energies with the same quantum numbers can be safely minimized. Once such minimization is converged, all the calculated energies and corresponding states can be obtained, simultaneously. The present VAP calculations are performed with time-odd Hartree-Fock Slater determinants. It is shown that the calculated VAP energies (both yrast and non-yrast) are very close to the corresponding ones from the full shell model calculations. It appears the present algorithm is not limited to the VAP, but should be universal, i.e., one can do the variation with different forms of the many-body wavefunctions to calculate the excited states in different quantum many-body systems.


21.60.Jz
                             

                              
New algorithm in the variation after projection calculations for non-yrast nuclear states 
    Y. S. Chen
    December 30, 2023
==========================================================================================




For a given many-body quantum system, its spectrum and the corresponding wave functions should be obtained by solving Schrödinger's equation. In nuclear structure studies, this is done in the full shell model(SM) calculations with a given Hamiltonian. However, due to the combinatorial computational cost, full SM calculations up to now have been restricted to rather small model spaces. To tackle the problem of the eigensystems in even larger model spaces, it looks that the only way is to compress the huge configuration space into a relatively small one so that the code can be run on a present-day computer. Unfortunately, the obtained energies and states are approximated ones. To make such approximated solutions as good as possible, various methods, such as shell model truncation<cit.>, stochastic quantum Monte Carlo approaches <cit.>, and the VAMPIR method <cit.> have been developed.

Recently, we implemented a variation after projection (VAP) calculation <cit.>. Instead of using the Hartree-Fock-Bogoliubov (HFB) vacuum state as adopted in the VAMPIR approach, we take a time-odd Hartree-Fock(HF) Slater determinant(SD). This means the particle number projection can be omitted, and only the spin projection is required. Although the VAP wave function is simply a spin projected SD, our calculations show that it can still be very close to the corresponding shell model wave function. This is because the involved spin projection plays a key role in obtaining a good shell model approximation<cit.>. However, that VAP can only be applied to the yrast states. In the present work, we construct a new and easy algorithm to extend our VAP calculations to the non-yrast excited states and show that the calculated non-yrast states are also very close to the exact shell model ones.

We should mention that the EXCITED VAMPIR can be used to calculate the excited states, too. In that method, one should first vary the HFB vacuum state to find the energy minimum for the yrast state. Then one varies the second HFB vacuum for the first excited state, but the Gram-Schmidt orthogonalization must be applied at each iteration to ensure the orthogonality between the first excited state and the yrast one. Similarly, the HFB vacua are added and varied one by one for the higher excited states, and the Gram-Schmidt orthogonalization should be treated through out all the VAMPIR iterations. To incorporate the important correlations, the more general EXCITED FED VAMPIRE (EFV) <cit.> uses several instead of only one symmetry-projected HFB vacuum for the description of each state. Those symmetry-projected HFB vacua for the same nuclear state are still added and varied successively. All the VAMPIR calculations then are ended with a final diagonalization of the given Hamiltonian in the space spanned by all the obtained projected HFB vacua.

The present VAP looks similar to the EFV in the sense that both of them take the projected basis. But here, we developed a quite different strategy in optimizing the low-lying nuclear states. The projected SDs are varied simultaneously and all the calculated states with the same spin (and parity if necessary) can be obtained on the same footing. The orthogonality among the calculated states is automatically fulfilled by solving the Hill-Wheeler equation. Thus, we do not need the Gram-Schmidt orthogonalization here.

To address the present algorithm, let us first introduce Cauchy's interlacing theorem relating to the eigenvalues of Hermitian matrices (Theorem 4.3.17 in <cit.>). Let A_n and B_m be Hermitian matrices of orders n and m, respectively. The eigenvalues of A_n are denoted by μ_1≤μ_2≤⋯≤μ_n. Those of B_m are λ_1≤λ_2≤⋯≤λ_m.
 if B_n-1 is of order m=n-1 and A_n is
    
    A_n
    =([ B_n-1    y⃗;  y⃗^†     a ]),

where, a is a real number and y⃗ is a complex vector of order n-1.
Then Cauchy's interlacing theorem tells us that

    μ_1≤λ_1≤⋯≤μ_j≤λ_j≤μ_j+1≤⋯≤λ_n-1≤μ_n.

One can immediately have the following conclusion  by repeated application of the above inequalities. For any principle submatrix B_m of A_n, we have (Theorem 4.3.28.in <cit.>)

    μ_j≤λ_j≤μ_j+n-m  (1≤ j ≤ m).


Since a unitary transformation does not change the eigenvalues of a given Hermitian matrix, the above conclusion can be naturally generalized to the following Poincaré separation theorem. Supposing that u⃗_1,u⃗_2,...u⃗_m are orthonormal complex vectors in the n-dimensional space where A_n is defined, the matrix elements (B_m)_ij=u⃗^†_i(A_n)u⃗_j form a Hermitian matrix B_m, and one can have the same inequalities (<ref>) from Poincaré separation theorem(see Theorem 4.3.37 in <cit.>)

However, if u⃗_1,u⃗_2,...u⃗_m are not orthonormal but independent (If not, we remove the redundant ones), for instance, the projected states used here, one may solve the following generalized eigenvalue equation of order m,

    ∑_j=1^m[(B_m)_ij-λ_k (N_m)_ij]f^k_j=0,

where (B_m)_ij=u⃗^†_i(A_n)u⃗_j and (N_m)_ij=u⃗^†_iu⃗_j. To solve Eq. (<ref>), the first step is the diagonalization of N_m and one has

    ∑_j=1^m(N_m)_ijR^k_j=n_kR^k_i,

 where n_k>0 and R^k with k=1,2,...m are eigenvalues and the corresponding eigenvectors, respectively. Then one can establish a new set of orthonormal vectors  v⃗_1,v⃗_2,...v⃗_m in the space spanned by u⃗_1,u⃗_2,...u⃗_m
 
    v⃗_k=∑_i=1^m R^k_iu⃗_i/√(n_k).

It is easy to prove that the eigenvalues, λ_k, in Eq.(<ref>) are identical to those of the Hermitian matrix C_m with the elements (C_m)_ij=v^†_i(A_n)v_j. According to the Poincaré separation theorem, we still have inequalities (<ref>).

In the practical shell model calculations, the Lanczos algorithm is adopted. The Lanczos matrix is enlarged with the iterations. With Cauchy's interlacing theorem, one can easily understand why all the calculated energies decrease monotonically until they converge to the full shell model energies.

Let us denote the full shell model energies for a given nuclear Hamiltonian by e_1≤ e_2⋯. Suppose that we are only interested in the lowest m nuclear states. In all approximated shell model methods, efforts have been made to try to find a proper configuration subspace, so that the calculated lowest m states are as close to the exact shell model ones as possible. Denoting the lowest m approximated energies in a configuration subspace by E_1≤ E_2≤⋯≤ E_m, then it is always true that

    E_j≥ e_j   (1≤ j ≤ m),

according to Cauchy's interlacing theorem or Poincaré separation theorem.

Therefore, one can define the non-negative energy differences δ E_j=E_j-e_j and the total energy difference

    Δ E=∑_j=1^m δ E_j=∑_j=1^m E_j-∑_j=1^m e_j.

It is obvious that if Δ E=0, then δ E_j=0 for all calculated states and we have the exact shell model results. However, in the approximated methods, Δ E>0, and one may expect that it should be as small as possible by adjusting the configuration subspace, in which the Hamiltonian is diagonalized. If Δ E reaches a minimum, then all the corresponding m lowest energies are determined, and their approximation can be tested by comparing them with the shell model results.

Since e_j energies in Eq.(<ref>) are fixed for a given Hamiltonian, instead of minimizing the Δ E, one can equivalently minimize the sum of the m lowest E_j energies, S_m≡∑_j=1^m E_j. This makes the calculations more practical in the case of large model space, in which e_j energies are impossible to be obtained.

Generally, it is unnecessary to restrict the form of the approximated nuclear wavefunctions. One may take the HF SDs, or HFB vacua with particle number projection, or their spin projected states, to form the configuration subspace. As a testing example of the present algorithm, here we take n different HF SDs differed by k, |Φ^k⟩, in the model space that the Hamiltonian is defined, and project them onto good spin. These projected states then form a configuration subspace. The form of the nuclear states with good spin J in that subspace can be written as

    |Ψ_JM⟩=∑_k=1^n∑_K=-J^Jf_KkP^J_MK|Φ^k⟩,

where, P^J_MK is the spin (or angular momentum) projection operator, and n is the number of adopted |Φ^k⟩ Slater determinants. |Φ^k⟩ can be established through the Thouless formula <cit.>,

    |Φ^k⟩=𝒩_ke^1/2∑_μνd^k_μνβ^k†_0,μβ^k†_0,ν|Φ^k_0⟩,

where 𝒩_k is the normalization parameter, and β^k†_0,μ is the quasiparticle creation operator for the |Φ^k_0⟩ vacuum. The initial |Φ^k_0⟩ SDs can be obtained in the same way as in Ref. <cit.>, and they are assumed to be different from each other. Once |Φ^k_0⟩ SDs are randomly chosen, they no longer change. Thus |Φ^k⟩ can be conveniently varied only by changing the d^k_μν parameters. Like Ref. <cit.>, we also assume the matrix elements, d^k_μν, to be complex numbers, i.e.,

    d^k_μν=x^k_μν+iy^k_μν,

where, x^k_μν and y^k_μν are real numbers. For simplicity, we use d⃗ to denote the vector of all the independent x^k_μν and y^k_μν parameters.
With a given d⃗, all the |Φ^k⟩ SDs are determined and one can establish the Hill-Wheeler equation

    ∑_Kk(H_K'k', Kk-EN_K'k', Kk) f_Kk=0

and the normalization condition

    ∑_K'k', Kkf^*_K'k'f_KkN_K'k', Kk=1,

where,

    H_K'k', Kk   =   ⟨Φ^k'|ĤP^J_K' K|Φ^k⟩,
    
    N_K'k', Kk   =   ⟨Φ^k'|P^J_K' K|Φ^k⟩.

Solving Eq.(<ref>), one can get a set of energies E_1≤ E_2≤⋯ and the corresponding states |Ψ^1_JM⟩,|Ψ^2_JM⟩,⋯.
Clearly, those energies are the functions of d⃗. Thus the sum of the lowest m energies, S_m, is also determined by d⃗, and we can write S_m=S_m(d⃗).

Therefore, one can minimize S_m by changing d⃗. For simplicity, we still call such energy minimization as variation after projection (VAP). To make the minimization efficient, the gradient of S_m(d⃗) should be calculated,

    ∇ S_m(d⃗)=∑_j=1^m∇ E_j(d⃗),

where, the components of ∇ E_j are actually ∂ E_j/∂ x^k_μν and ∂ E_j/∂ y^k_μν, whose expressions can be easily deduced based on the formulas in Refs. <cit.>.

Once the quantities of S_m(d⃗) and its gradient ∇ S_m(d⃗) are available, one might perform the minimization of S_m(d⃗) using the quasi-Newton method. This may work if the Hessian at the minimum is positive definite. Unfortunately, we have learned that the Hessian of the projected energy at the the minimum is actually semi-positive definite <cit.>. This might be the reason that the traditional quasi-Newton method does not work well in our VAP calculations. To make the VAP calculation converges more reliably, we calculated the exact Hessian matrix, denoted by 𝐇[S_m(d⃗)]. Similar to Eq.(<ref>), we also have

    𝐇[S_m(d⃗)]=∑_j=1^m 𝐇[E_j(d⃗)],

where the matrix elements of 𝐇[E_j(d⃗)] are ∂^2 E_j/∂ x^k_μν∂ x^k'_μ'ν', ∂^2 E_j/∂ x^k_μν∂ y^k'_μ'ν' and ∂^2 E_j/∂ y^k_μν∂ y^k'_μ'ν'. The expressions for such second derivatives can also be derived without much difficulty based on the deductions in Refs. <cit.>.

To check the correctness of the above calculated quantities, one can use the following equations by definition,

    lim_δd⃗→0S_m(d⃗+δd⃗)-S_m(d⃗)/|δd⃗|=∇ S_m(d⃗)·e⃗,
          lim_δd⃗→0∇ S_m(d⃗+δd⃗)-∇ S_m(d⃗)/|δd⃗|·e⃗
       =   e⃗^t·𝐇[S_m(d⃗)]·e⃗,

where, e⃗=δd⃗/|δd⃗|, which shows the direction of δd⃗. In the present numerical calculations, the calculated S_m, its gradient and Hessian matrix indeed fulfill the above equations at arbitrary point d⃗ and e⃗. With a small |δd⃗|=10^-4, the calculated differences between both sides of Eqs. (<ref>) and (<ref>) are usually less than 10^-6.



Once the above quantities are available, we start from d⃗=0 and take the exact trust region method <cit.> to search for the minimum of S_m(d⃗).  At the minimum, the gradient ∇ S_m(d⃗) should be zero. Here, the VAP iteration terminates if |∇ S_m(d⃗)|≤ 0.01keV. This is a very strict condition so that the obtained minimum is precise.

To test the validity of the present algorithm, we performed the VAP calculations in the sd shell and take the USDB interaction <cit.>. The simplest calculation for the non-yrast states is the one with n=2 and m=2, i.e., at a given spin, we take two SDs to calculate the lowest two nuclear states, simultaneously. The calculated nuclei are ^24Mg, ^25Mg, ^26Mg, and ^26Al, whose calculated energies are shown in Fig.<ref>. It looks that all the calculated VAP energies are very close to the shell model ones calculated with the NuSHellX<cit.>. To show the differences between the VAP and the shell model more clearly, we calculated the δ E_j (j=1,2) from Fig.<ref>, and show them in Fig.<ref>. One can see that all the δ E_j are indeed nonnegative, as predicted by Cauchy's interlacing theorem. It is interesting to see that the energy differences of the yrare states, δ E_2, are roughly close to the yrast ones, δ E_1 after the minimization of S_m=2.



 The δ E_i values are expected to be as small as possible so that one can obtain satisfactory
  shell model approximations. This can be realized by adding more SDs to the VAP wave functions.
   As a testing example, we only calculated the yrast states (i.e., m=1) with the cases n=1,2 and 3.
 The calculated δ E_1 energy differences are shown in Fig.<ref>. The case with n=1 is exactly the Fig.2 in Ref.<cit.>. For the n=2 case, all the δ E_1 values (except for the ones at the two highest spins) are smaller than the corresponding ones with n=1. Notice that for the two highest spin states, we have δ E_1=0 in all cases due to the unphysical tiny configuration space. With more SDs added in, the VAP results become more closer to the shell model ones.
 One can see that, again, the δ E_1 with n=3 becomes even smaller.




However, in the practical application, it may not be enough that only two lowest states are calculated. Here, to show the validity of the present algorithm, we calculated the lowest 10 states (m=10) with 10 SDs (n=10) for the J=0 states in even-even sd shell nuclei. The results are shown in Fig.<ref>. Like Fig.<ref>, the VAP energies are also very close to the corresponding SM ones. Again, all the δ E_i are non-negative. Especially for the ^20Ne, ^28Ne and ^36Ar nuclei, we have δ E_i=0 for all the calculated states. This means we got the exact shell model results for these 3 nuclei and confirms the correctness of the present calculations. Most of  δ E_i values are in the range from 0 to 300keV. Some of them are a little larger but the largest one is still within 600keV. Of course, such approximation can be further improved by increasing the number, n, of the included SDs.


It may be interesting to compare the present algorithm with the one taken by the VAMPIR. Let us first calculate a simplest case of m=1 and n=2 for the ground 0^+ state in ^24Mg. The USDB interaction is still adopted and the exact ground state energy in ^24Mg is -87.104MeV. We vary these two SDs simultaneously in the present VAP and obtained the lowest minimum -87.039MeV. Now, let us follow the method of FED VAMPIR <cit.>.  Varying the first SD,  one can get the energy  minimum of -86.936MeV. Then a second SD is added to improve the approximation. Fixing the obtained first SD and varying the second one, this makes the ground state energy lower to -87.008MeV which still lies above our -87.039MeV. This clearly shows the SDs obtained by the FED VAMPIR algorithm are not fully optimized. As can be seen that if one fixes the second SD and comes back to vary the first SD again, then the ground state energy drops further from -87.008MeV to -87.015MeV.

Another comparison is the case of m=2 and n=2 for the ground 0^+ state and the first excited 0^+ state in ^24Mg. Our algorithm gives the lowest S_2 to be -166.474MeV which corresponds to the lowest two 0^+ energies E_1=-86.944MeV and E_2=-79.530MeV. While the exact SM results are e_1=-87.104MeV and e_2=-79.766MeV. If we follow the algorithm in the EXCITED VAMPIR <cit.>, then the ground state energy obtained by varying the first SD is E_1=-86.936MeV as mentioned in the last paragraph. For the minimization of E_2, the first SD will be fixed and the second SD will be varied under the Gram-Schmidt orthogonalization. This is essentially equivalent to the minimization of S_2=E_1+E_2 because E_1 is fixed. After the minimization of E_2, the VAMPIR performs the final dialgonalization in the space spanned by the two optimized projected SD states and the final lowest two energies can be obtained. In the present case of J^π=0^+, such dialgonalization does not change S_2. Therefore, it is clear that the whole procedure of the EXCITED VAMPIR is exactly equivalent to the minimization of S_2 in our VAP but with the first SD fixed. In this way, the obtained S_2=166.341MeV, corresponding to E_1=-86.940MeV and E_2=-79.400MeV. Both energies are above the present VAP results, especially E_2, which is 130keV above our -79.530MeV. This again shows that the projected states in the VAMPIR method are not fully optimized.

Let us study the Hessian matrix by calculating its eigenvalues. In Fig.<ref>, we show the eigenvalues of the Hessian matrix corresponding to the VAP calculations of ^24Mg in Fig.<ref>. Clearly, the Hessian matrix is indefinite at the starting point. But after convergence, it becomes semipositive definite. This is very similar to the case of m=1 and n=1 as discussed in Ref.<cit.>.


Due to the simplicity of the present algorithm, all the above VAP calculations are performed with a common code, except for taking different m,n (n≥ m) and the spin J as input parameters. We hope the present VAP will be an alternative way of extending the shell model calculations in larger model spaces. Calculations for heavier nuclei will be performed after our VAP code is parallelized.

 Equation (<ref>) is essentially a mathematical conclusion. The present method of the S_m minimization may be universal and applicable to other quantum many-body systems; for instance, the excited electronic states in a chemical system <cit.>. One may also change the form of the approximated wavefunction in a convenient way. For example, if one simply takes the superposition of the deformed SDs without spin projection, then the minimization of S_m should be a natural extension of the Hartree-Fock method for the excited states. Such calculation can be implemented by simply removing the spin projection from the present VAP code.

Acknowledgements The work is supported by the
National Natural Science Foundation of China under Contract Nos.
11575290, 11475072 and 11490560.






















0

Horoi94 M. Horoi, B.A. Brown, and V. Zelevinsky, Phys. Rev. C50, 2274(1994).

Koonon97 S.E. Koonin, D.J. Dean, K. Langanke, Phys. Rep. 278,1(1997).

Otsuka01 T. Otsuka, M. Honma, T. Mizusaki, N. Shimizu, and Y. Utsuno, Prog. Part. Nucl. Phys. 47, 319 (2001).

Schmid04 K.W. Schmid, Prog. Part. Nucl. Phys. 52, 565 (2004).

tuya17 T. Ya, Y. He, Z.-C. Gao, J.-Q. Wang, and Y. S. Chen, Physical Review C95, 064307(2017).

gao15 Z.-C. Gao, Mihai Horoi, Y.S.Chen Physical Review C92, 064310(2015).

Schmid89 K.W. Schmid, R.-R. Zheng, F. Gr ümmer and A. Faessler, Nucl. Phys. A499, 63(1989).

mt  R.A. Horn and C.R. Johnson, Matrices Analysis (Cambridge Univeristy Press, 2013)

Ring80 P. Ring and P. Schuck, The Nuclear Many-Body Problem
(Springer Verlag, New York, Heidelberg, Berlin, 1980).

hu14 Q.-L. Hu, Z.-C. Gao, Y. S. Chen, Phys. Lett.  B734, 162 (2014).

numopt J. Norcedal and S. J. Wright, Numerical Optimization
(Springer Verlag, New York, Berlin, Heidelberg, 2006).

usdb B. A. Brown and W. A. Richter, Phys. Rev. C74, 034315 (2006).

nushellx B. A. Brown, W. D. M. Rae, E. McDonald and M. Horoi, NuShellX@MSU,
http://www.nscl.msu.edu/brown/resources/resources.html;
W. D. M. Rae, NuShellX, http://www.garsington.eclipse.co.uk/.

carlos13 C. A. Jimnez-Hoyos, R. Rodrguez-Guzmn, and G. E. Scuseria, J. Chem. Phys.139, 224110 (2013)



