# Mixing of one-particle-one-hole projected states with the variation after
projection wave functions

Xiao Lu  China Institute of Atomic Energy, Beijing 102413, China  Zhan-Jiang
Lian  China Institute of Atomic Energy, Beijing 102413, China  Zao-Chun Gao
China Institute of Atomic Energy, Beijing 102413, China

###### Abstract

In this paper, we study the mixing of one-particle-one-hole projected states
with the variation after projection (VAP) wave functions in attempt to improve
the approximation of this method. It turns out that when minimizing only the
lowest (yrast) energy with given spin and parity, the one-particle-one-hole
projected states can not be mixed into the converged VAP wave function, which
is very similar to the situation of the Hartree-Fock method. However, if one
minimizes the sum of several lowest energies with the same spin and parity,
the one-particle-one-hole mixing can make some improvements to the VAP wave
functions. We expect such one-particle-one-hole mixing may be useful in the
calculations of the low-lying excited states in heavy nuclei with a large
model space.

##  I Introduction

The nuclear shell model (SM) Caurier05  has been very successful in describing
various properties of nuclei, especially in the neighborhood of the closed
shells, where the configuration space is usually small but good enough for the
construction of the nuclear wave function. However, in heavy deformed nuclear
region, the configuration space is huge, in which the full SM calculation is
almost impossible. Such huge configuration space must be compressed so that SM
calculation can be performed on a present-day computer. Unfortunately, the
energies and wave functions obtained in a compressed configuration space are
approximated ones. To make such approximated solutions as close as possible to
the exact shell model ones, various methods, such as shell model truncation
ST94 , stochastic quantum Monte Carlo approaches MC97 ; MC01 , projected
configuration interaction PCI09 , and the class of variation after projection
(VAP) methods VAMPIR04 ; Shimizu21 ; VAP15 ; VAP17 ; VAP18 ; VAP21 , have been
developed.

Among those approximated SM methods, the VAP method is an important one, in
which the nuclear wave functions with good quantum numbers are sufficiently
optimized VAMPIR04 ; Shimizu21 ; VAP15 ; VAP17 ; VAP18 ; VAP21 . In Refs.
VAMPIR04 ; Shimizu21 , the Hartree-Fock-Bogoliubov (HFB) vacuum states are
adopted, and the neutron and proton number projections are performed in
addition to the angular momentum projection. Very good approximations of their
results have been achieved. However, the full projection of the HFB state onto
good quantum numbers (i.e., neutron number, proton number, angular momentum
and parity) requires five-fold integration, which is too much time consuming.
On the other hand, in Ref. VAP15 , we have shown that the angular momentum
projection is crucial in obtaining good shell model approximation. So, a
simpler VAP can be the one in which the Hartree-Fock (HF) Slater determinants
(SDs) are taken and only the angular momentum projection and parity projection
are performed VAP17 ; VAP18 ; VAP21 . Apparently, with the same number of
projected states, the approximation of the VAP with projected SDs is not as
good as that with projected HFB vacua. This is because HFB vacuum incorporates
more correlations especially pairing. However, the approximation of the former
can be continuously improved by adding more and more projected SDs. In this
sense, we are more interested in the VAP with projected SDs due to its
relatively less computational cost. Thus we only discuss the VAP with
projected SDs in this paper.

In Ref. VAP18 , we have proposed a new algorithm of the VAP calculation, in
which the yrast state and non-yrast states can be varied on the same footing.
In this algorithm, the orthogonality among the calculated states is
automatically fulfilled by solving the Hill-Wheeler (HW) equation. This avoids
the complexity of the frequently used Gram-Schmidt orthogonalization, as
adopted in Ref. VAMPIR04 . After that, we further found that the complicated
$K$-mixing can be safely removed from VAP with the adopted SDs fully symmetry-
unrestricted. This considerably simplifies the VAP calculation, especially in
the application to the high-spin states with arbitrary deformation VAP21 .

However, if the VAP method is applied to heavy deformed nuclei, the
computational burden becomes quite heavy. Under this situation, the VAP
calculation with large number of projected states becomes impractical. A
further improvement of such VAP wave function may be the inclusion of the
projected particle-hole states built on top of the converged reference SDs.
This is very similar to the cases of the Projected Shell Model SY95  and the
MONSTER method VAMPIR04 , in which the projeced quasi-particle states on top
of the selected projected vacuum are included to form the nuclear wave
function.

The mixing of particle-hole projected states with the VAP wave functions
requires more complicated projected matrix elements to form a wider Hill-
Wheeler (HW) equation. Fortunately, in the present VAP calculation VAP17 ;
VAP18 ; VAP21 , all the matrix elements needed in the one-particle-one-hole
(1p-1h) mixing are already available. These matrix elements are originally
used to build the Hessian matrix of the VAP energy, so that the VAP
calculation may converge after a very few iterations (see Fig. 1 below). Thus
the improved energies and the nuclear wave functions with 1p-1h mixing can be
calculated without much computational cost. However, if one further considers
the 2p-2h mixing, there are two troubles. The first one is, new projected
matrix elements among the 2p-2h configurations must be calculated
specifically. The other one is the number of the 2p-2h SDs is much larger than
the 1p-1h one. Therefore the calculation including 2p-2h mixing is much more
complicated. Here, for simplicity, we only discuss the 1p-1h mixing in this
paper and see how it affects the VAP wave function.

The paper is organized as follows. Section II provides a general introduction
of the adopted method. Section III discusses the 1p-1h mixing after the VAP
for the yrast states. Section IV discusses the 1p-1h mixing after the VAP for
the non-yrast states. Summary and outlook are presented in Section V.

##  II The VAP Method and the 1p-1h mixing

Let’s start the introduction of the present method with the simplest case that
only a single SD, $|\Phi\rangle$, is considered. At given quantum numbers of
spin, $J$, parity, $\pi$, and magnetic quantum number, $M$, one can generate
$2J+1$ projected states, $P_{MK}^{J\pi}|\Phi\rangle$, from $|\Phi\rangle$.
Here, $P_{MK}^{J\pi}$ stands for the product of the angular momentum
projection operator, $P_{MK}^{J}$, and the parity projection operator,
$P^{\pi}$. Generally, all these $2J+1$ projected states are expected to be
taken to form the trial nuclear wave function,

$|\Psi_{J\pi
M\alpha}\rangle=\sum_{K=-J}^{J}f_{K}^{J\pi\alpha}P_{MK}^{J\pi}|\Phi\rangle,$
(1)

where, $\alpha$ is used to differ the states with the same $J$, $\pi$ and $M$.
The coefficients, $f_{K}^{J\pi\alpha}$, and the corresponding energy,
$E^{J\pi}_{\alpha}$, are determined by solving the following Hill-Wheeler (HW)
equation,

$\displaystyle\sum_{K^{\prime}=-J}^{J}(H^{J\pi}_{KK^{\prime}}-E^{J\pi}_{\alpha}N^{J\pi}_{KK^{\prime}})f^{J\pi\alpha}_{K^{\prime}}=0,$
(2)

where
$H^{J\pi}_{KK^{\prime}}=\langle\Phi|\hat{H}P^{J\pi}_{KK^{\prime}}|\Phi\rangle$
and $N^{J\pi}_{KK^{\prime}}=\langle\Phi|P^{J\pi}_{KK^{\prime}}|\Phi\rangle$.
For convenience, we assume $E^{J\pi}_{1}\leq E^{J\pi}_{2}\leq\cdots\leq
E^{J\pi}_{2J+1}$. The coefficients $f^{J\pi\alpha}_{K}$ should satisfy the
normalization condition,

$\displaystyle\sum_{K,K^{\prime}=-J}^{J}f^{J\pi\alpha*}_{K}N^{J\pi}_{KK^{\prime}}f^{J\pi\alpha}_{K^{\prime}}=1.$
(3)

$|\Psi_{J\pi M\alpha}\rangle$ needs to be optimized so that it can be as close
as possible to the exact shell model one. According to Eq. (2), $|\Psi_{J\pi
M\alpha}\rangle$ is completely determined by $|\Phi\rangle$ provided the
Hamiltonian is given. Thus the optimization of Eq. (1) is actually realized by
varying the $|\Phi\rangle$ state.

To vary the $|\Phi\rangle$ state, one may first need to parameterize it.
Suppose there is a normalized HFB vacuum state, $|\Phi_{0}\rangle$. The
corresponding quasiparticle operators are denoted by $\beta^{\dagger}_{0,\mu}$
and $\beta_{0,\mu}$. Using the Thouless theorem TH60 ; PR80 ,
$|\Phi_{0}\rangle$ can be changed to a new HFB vacuum state $|\Phi\rangle$,
namely,

$|\Phi\rangle=\mathcal{N}e^{\frac{1}{2}\sum_{\mu\nu}d_{\mu\nu}\beta_{0,\mu}^{\dagger}\beta_{0,\nu}^{\dagger}}|\Phi_{0}\rangle=\mathcal{N}e^{\sum_{\kappa}d_{\kappa}A^{\dagger}_{\kappa}}|\Phi_{0}\rangle,$
(4)

where $d$ is a complex skew matrix. The matrix elements $d_{\mu\nu}$ will be
considered as the variational parameters. $\mathcal{N}$ is the normalization
parameter, so that $\langle\Phi|\Phi\rangle=1$. For convenience, the subscript
$\kappa$ is used to stand for the $(\mu,\nu)$ numbers with $\mu<\nu$, and the
particle pair operators are defined as,

$\displaystyle\hat{A}_{\kappa}^{\dagger}$ $\displaystyle=$
$\displaystyle\beta_{0,\mu}^{\dagger}\beta_{0,\nu}^{\dagger},$ (5)
$\displaystyle\hat{A}_{\kappa}$ $\displaystyle=$
$\displaystyle(\beta_{0,\mu}^{\dagger}\beta_{0,\nu}^{\dagger})^{\dagger}=\beta_{0,\nu}\beta_{0,\mu}.$
(6)

Here, the vacuum states $|\Phi_{0}\rangle$ and $|\Phi\rangle$ are reduced to
be SDs, so that the particle number projections can be omitted. Therefore, the
$\hat{A}_{\kappa}^{\dagger}$ operators become 1p-1h operators corresponding to
the $|\Phi_{0}\rangle$ state.

Now, let’s present a brief introduction on how to establish the VAP iteration.
Suppose that we have a $|\Phi\rangle$ state at certain VAP iteration. Then we
set $|\Phi_{0}\rangle=|\Phi\rangle$ and $d=0$. At this $d=0$ point, the
energy, $E^{J\pi}_{\alpha}$, and the corresponding $f_{K}^{J\pi\alpha}$
coefficients are calculated first by solving Eq. (2). This process includes
two successive diagonalizations. The first one is the diagonalization of the
norm matrix $N^{J\pi}$ in Eq. (2) and we have

$\displaystyle\sum_{K^{\prime}=-J}^{J}N^{J\pi}_{KK^{\prime}}R^{k}_{K^{\prime}}=\sigma_{k}R^{k}_{K},$
(7)

where $\sigma_{k}\geq 0$ and $R^{k}$ with $k=1,2,...2J+1$ are eigenvalues and
the corresponding eigenvectors, respectively. Then one can establish a new set
of orthonormal basis states, $|\psi_{k}^{J\pi}\rangle\,(k=1,2,\cdots,2J+1)$,

$\displaystyle|\psi_{k}^{J\pi}\rangle=\frac{1}{\sqrt{\sigma_{k}}}\sum_{K=-J}^{J}R^{k}_{K}P^{J\pi}_{MK}|\Phi\rangle,$
(8)

and Eq. (2) can be transformed into a normal eigenvalue equation

$\displaystyle\sum_{k^{\prime}=1}^{2J+1}\left[\langle\psi_{k}^{J\pi}|\hat{H}|\psi_{k^{\prime}}^{J\pi}\rangle-E^{J\pi}_{\alpha}\delta_{kk^{\prime}}\right]u^{J\pi\alpha}_{k^{\prime}}=0.$
(9)

The second diagonalization is then performed on Eq. (9), so that the energies,
$E^{J\pi}_{\alpha}$, in Eq. (2) are obtained and the coefficients of the wave
functions, $f^{J\pi\alpha}_{K}$, are transformed from $u^{J\pi\alpha}_{k}$,
i.e.,

$\displaystyle
f^{J\pi\alpha}_{K}=\sum_{k=1}^{2J+1}\frac{R^{k}_{K}u^{J\alpha}_{k}}{\sqrt{\sigma_{k}}}.$
(10)

Notice that, in Eq. (8), $\sigma_{k}$ should not be zero, or the corresponding
$|\psi_{k}^{J\pi}\rangle$ should be abandoned. However, if $\sigma_{k}$ is too
small, the $|\psi_{k}^{J\pi}\rangle$ basis states may not be calculated
precisely enough to guarantee the stability of the calculated energies and
wave functions. This is one of the main troubles when all $2J+1$ projected
basis states are included in the VAP calculation, because some of the smallest
$\sigma_{k}$ values could become tiny and tiny and finally the VAP iteration
collapses.

After we have the projected energies, $E^{J\pi}_{\alpha}$, and the
corresponding wave functions, we chose the lowest one, $E^{J\pi}_{1}$, and
minimize it. Thus the gradient and the Hessian of $E^{J\pi}_{1}$ in the space
spanned by the $d_{\kappa}$ parameters are calculated according to the
formulas in Ref. VAP17 , then the next improved SD $|\Phi\rangle$ can be
determined by adopting the trust region Newton’s algorithm numopt . The VAP
iteration terminates when the gradient of $E^{J\pi}_{1}$ becomes less than
0.01keV, which is precise enough to obtain the exact minimum of
$E^{J\pi}_{1}$.

Figure 1: (Color online) Quantities as functions of the VAP iteration for the
$J^{\pi}$= 8+ state in 26Mg. Results with Eq. (1) are shown as open circles.
Those with Eq. (11) are shown as filled dots in different colors. The
calculated quantities are (a) the norms, $\sigma_{k}$; (b)the absolute values
of the energy gradients; and (c) the VAP energies. The iteration terminates
once the absolute value of the energy gradient is less than 0.01keV. The USDB
interaction is adopted.

In Fig. 1, we demonstrate such VAP iteration using the example of the yrast
$8^{+}$ state in 26Mg, which will be further discussed in the next Section.
The initial SD $|\Phi\rangle$ is randomly chosen. The adopted Hamiltonian is
USDB usdb . The results are shown as black circles. From this initial
$|\Phi\rangle$, the VAP converges normally after 43 iterations.

However, according to Ref. VAP21 , the form of the trial VAP wave function
with spin, $J$, can be simplified by adopting just one projected state rather
than adopting all $2J+1$ angular momentum projected states for each selected
reference state. To confirm this point, we chose only one of the $2J+1$
projected states in Eq. (1), and the trial VAP wave function can be
considerably simplified as,

$|\Psi_{J\pi
M}(K)\rangle=\frac{P_{MK}^{J\pi}|\Phi\rangle}{\sqrt{\langle\Phi|P_{KK}^{J\pi}|\Phi\rangle}},$
(11)

where, the $K$ number can be randomly chosen. Here, we still take the above
example and choose $K=0,2$ and $5$ to optimize Eq. (11) independently by using
exactly the same VAP code. The results are shown as filled dots with different
colors in Fig. 1. It is shown that all the converged energies for $|\Psi_{J\pi
M}(K)\rangle$ with different $K$ are almost the same and also very close to
the converged $E^{J\pi}_{1}$.

Now, we have one optimized wave function of Eq. (1) and three optimized ones
of Eq. (11) with $K=0,2$ and $5$, which must be compared with one another.
Overlaps among these wave functions are calculated. The values of
$|\langle\Psi_{J\pi M1}|\Psi_{J\pi M}(K)\rangle|$ corresponding to the
converged results in Fig. 1 are 0.9867,0.9853,0.9898 for $K=0,2$ and $5$,
respectively. Actually, results with other $K$ numbers are almost the same as
in Fig. 1. This clearly shows the equivalence of Eq. (1) and Eq. (11) and
implies the equivalence among the $|\Psi_{J\pi M}(K)\rangle$ states
themselves. The latter is confirmed by calculating the quantities of
$|\langle\Psi_{J\pi M}(K)|\Psi_{J\pi M}(K^{\prime})\rangle|$, which turn out
to be 0.9930, 0.9911 and 0.9875 for $(K,K^{\prime})=$(0,2), (0,5) and (2,5),
respectively.

We should stress that the above equivalences are valid only for the optimized
VAP wave functions. With the same $|\Phi\rangle$, Eq. (1) is expected to be
better than Eq. (11). This can be seen in Fig. 1 that the energy of Eq. (1) is
indeed much lower than the ones of Eq. (11) at the first iteration where the
same initial $|\Phi\rangle$ is used to start all the VAP calculations.
However, in the process of VAP iteration, the $|\Phi\rangle$ state in
$|\Psi_{J\pi M1}\rangle$ and those ones in the $|\Psi_{J\pi M}(K)\rangle$
states change independently, and they may become different from one another.
This clearly tells us the converged VAP wave functions with different
reference SD states can be almost equivalent.

The advantages of introducing Eq. (11) into VAP are obvious. First of all, the
number of all required matrix elements with Eq. (11) is only one
$(2J+1)^{2}$th of that with Eq. (1). Consequently, the computational time with
Eq. (11) is remarkably reduced. In the example of Fig. 1, the elapsed time for
each iteration with Eq. (11) is only about 3s on an Intel Xeon CPU with 20
cores, while such time with Eq. (1) is about 360s on the same CPU. Secondly,
the VAP iteration with Eq. (11) converges more reliably and faster than that
with Eq. (1). As one can see from Fig. 1(a), the smallest $\sigma_{k}$ values
corresponding to Eq. (1) tend to be tiny and tiny as the iteration proceeds.
Fortunately, they come back after the 30th iteration in this calculation.
However, this is not always lucky. In many cases, the $\sigma_{k}$ values
could keep on moving to zero and one can not get a converged VAP wave function
with Eq. (1) (See the examples in the supplemental material for Ref. VAP21 ).
Finally, the VAP calculation with Eq. (11) can be conveniently extended to
arbitrary high-spin states, while the VAP with Eq. (1) is very difficult to do
that.

Recognizing the equivalence of Eq. (1) and Eq. (11) in VAP, we prefer to take
the simplified Eq. (11) and extend it by including more reference SDs,

$|\Psi_{J\pi
M\alpha}(K)\rangle=\sum_{i=1}^{n}f_{i}^{J\pi\alpha}P_{MK}^{J\pi}|\Phi_{i}\rangle.$
(12)

The coefficients $f_{i}^{J\pi\alpha}$ and the corresponding energy,
$E^{J\pi}_{\alpha}$, can be obtained by solving the following HW equation,

$\sum_{i^{\prime}=1}^{n}\langle\Phi_{i}|(\hat{H}-E^{J\pi}_{\alpha})P_{KK}^{J\pi}|\Phi_{i^{\prime}}\rangle
f_{i^{\prime}}^{J\pi\alpha}=0.$ (13)

These $f_{i}^{J\pi\alpha}$ coefficients also should satisfy the normalization
condition,

$\sum_{ii^{\prime}=1}^{n}f_{i}^{J\pi\alpha*}\langle\Phi_{i}|P_{KK}^{J\pi}|\Phi_{i^{\prime}}\rangle
f_{i^{\prime}}^{J\pi\alpha}=1.$ (14)

Here and below, all $E^{J\pi}_{\alpha}$ energies are calculated from Eq. (13)
rather than Eq. (2) since we no longer use Eq. (1). In this paper, we assume
$E^{J\pi}_{1}\leq E^{J\pi}_{2}\leq\cdots\leq E^{J\pi}_{n}$. But in practical
calculations we are only interested in the lowest $m(\leq n)$ energies.

From Eq. (13), one can see that changes of the $|\Phi_{i}\rangle$ states may
directly change the $E^{J\pi}_{\alpha}$ energies. Our target is to try to find
a set of these $n$ $|\Phi_{i}\rangle$ states, so that the calculated lowest
$m$ $E^{J\pi}_{\alpha}$ energies can be as close as possible to the
corresponding exact shell model ones denoted by $e_{\alpha}^{J\pi}$. Here, we
also assume $e^{J\pi}_{1}\leq e^{J\pi}_{2}\leq\cdots\leq e^{J\pi}_{m}$.

Thanks to the Cauchy’s interlacing theorem, which makes it clear that
$E_{\alpha}^{J\pi}\geq e_{\alpha}^{J\pi}$ for any excited energies. This has
been strictly proved in Ref. VAP18 . Recently, we just recognized there is a
famous Hylleraas-Undheim-MacDonald (HUM) theorem widely known in the field of
quantum chemistry hum1 ; hum2 . This theorem exactly tells us the relation of
$E_{\alpha}^{J\pi}\geq e_{\alpha}^{J\pi}$.

Therefore, one can define the non-negative energy difference $\delta
E_{\alpha}=E_{\alpha}^{J\pi}-e_{\alpha}^{J\pi}$. Then the summation of $\delta
E_{\alpha}$ for the lowest $m$ states must be non-negative, too. Namely,

$\displaystyle\Delta E_{m}=\sum_{\alpha=1}^{m}\delta
E_{\alpha}=\sum_{\alpha=1}^{m}E_{\alpha}^{J\pi}-\sum_{\alpha=1}^{m}e_{\alpha}^{J\pi}\geq
0.$ (15)

Clearly, if $\Delta E_{m}=0$, then $\delta E_{\alpha}=0$ for all included
states. This means $E_{\alpha}^{J\pi}=e_{\alpha}^{J\pi}$ and we obtain the
exact eigenenergies of the given $\hat{H}$. But in general, $\Delta E_{m}>0$.
Then our target is to try to find a number ($n$) of basis states through
variation, so that $\Delta E_{m}$ becomes a minimum. Then the corresponding
$E_{\alpha}^{J\pi}$ energies at the minimum of $\Delta E_{m}$ can be
simultaneously obtained. These $E_{\alpha}^{J\pi}$ energies can be compared
with the exact $e_{\alpha}^{J\pi}$ ones to test the quality of the present VAP
method VAP18 . When $m=1$, we come back to the Ritz variational principle.

However, $\Delta E_{m}$ is usually unknown because we don’t know the exact
$e_{i}^{J\pi}$ energies in most cases. Actually, with a given $\hat{H}$, the
$e_{i}^{J\pi}$ energies are determined, and can be considered as constants in
Eq. (15). So minimizing $\Delta E_{m}$ is equivalent to minimizing

$\displaystyle S_{m}=\sum_{\alpha=1}^{m}E_{\alpha}^{J\pi}.$ (16)

Since the involved $|\Phi_{i}\rangle$ states are varied independently, each of
them has its own $d$ matrix, which is denoted by $d_{i}$. The matrix elements
of $d_{i}$ are then denoted by $d_{i\kappa}$, namely

$|\Phi_{i}\rangle=\mathcal{N}_{i}e^{\sum_{\kappa}d_{i\kappa}A^{\dagger}_{i\kappa}}|\Phi_{0,i}\rangle.$
(17)

For simplicity, we use $\vec{d}$ to denote the vector of all the independent
$d_{i\kappa}$ parameters. Given a certain $\vec{d}$, all $|\Phi_{i}\rangle$
states can be uniquely determined. Consequently, the quantity $S_{m}$,
determined by the set of $|\Phi_{i}\rangle$ states, can be considered as a
function of $\vec{d}$. Therefore, the VAP calculation is actually the
minimization of $S_{m}$ in the space spanned by all the $d_{i\kappa}$
parameters.

Let us give a brief introduction to the present VAP iteration. Starting from a
set of randomly chosen $|\Phi_{0,i}\rangle$ states, the gradient and the
Hessian matrix of $S_{m}$ at $\vec{d}=0$ are calculated. Then a new point
$\vec{d}_{\text{m}in}\neq 0$ can be searched so that $S_{m}$ is as low as
possible. This $\vec{d}_{\text{m}in}$ determines a new set of
$|\Phi_{i}\rangle$ states. We then update the $|\Phi_{0,i}\rangle$ states with
these $|\Phi_{i}\rangle$ ones and perform the next VAP iteration. Notice that
we only calculate the gradient and the Hessian matrix at $\vec{d}=0$. This is
because at this special point, $|\Phi_{i}\rangle=|\Phi_{0,i}\rangle$, and the
particle-hole operators $\hat{A}_{i\kappa}^{\dagger}$ become the ones
corresponding to $|\Phi_{i}\rangle$. This may considerably simplify the VAP
formulation and the corresponding calculation.

Now, let us consider the 1p-1h mixing. Suppose that we have a set of reference
SD states, $|\Phi_{i}\rangle$, then the expanded wave function including 1p-1h
components can be written as,

$\displaystyle|\Psi_{J\pi M\alpha}^{\prime}(K)\rangle$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}{f^{\prime}}_{i}^{J\pi\alpha}P_{MK}^{J\pi}|\Phi_{i}\rangle$
(18)
$\displaystyle+\sum_{i=1}^{n}\sum_{\kappa}{f^{\prime}}_{i\kappa}^{J\pi\alpha}P_{MK}^{J\pi}A_{i\kappa}^{\dagger}|\Phi_{i}\rangle.$

Here, the $A_{i\kappa}^{\dagger}$ operators are the 1p-1h operators
corresponding to $|\Phi_{i}\rangle$. The coefficients,
${f^{\prime}}_{i^{\prime}(\kappa^{\prime})}^{J\pi\alpha}$, and the
corresponding energy, ${E^{\prime}}^{J\pi}_{\alpha}$, can be obtained by
solving the following expanded HW equations,

$\left\\{\begin{aligned}
&\sum_{i^{\prime}=1}^{n}\langle\Phi_{i}|(\hat{H}-{E^{\prime}}^{J\pi}_{\alpha})P_{KK}^{J\pi}|\Phi_{i^{\prime}}\rangle{f^{\prime}}_{i^{\prime}}^{J\pi\alpha}\\\
&+\sum_{i^{\prime}=1}^{n}\sum_{\kappa^{\prime}}\langle\Phi_{i}|(\hat{H}-{E^{\prime}}^{J\pi}_{\alpha})P_{KK}^{J\pi}A_{i^{\prime}\kappa^{\prime}}^{\dagger}|\Phi_{i^{\prime}}\rangle{f^{\prime}}_{i^{\prime}\kappa^{\prime}}^{J\pi\alpha}=0,\\\
&\sum_{i^{\prime}=1}^{n}\langle\Phi_{i}|A_{i\kappa}P_{KK}^{J\pi}(\hat{H}-{E^{\prime}}^{J\pi}_{\alpha})|\Phi_{i^{\prime}}\rangle{f^{\prime}}_{i^{\prime}}^{J\pi\alpha}\\\
&+\sum_{i^{\prime}=1}^{n}\sum_{\kappa^{\prime}}\langle\Phi_{i}|A_{i\kappa}(\hat{H}-{E^{\prime}}^{J\pi}_{\alpha})P_{KK}^{J\pi}A_{i^{\prime}\kappa^{\prime}}^{\dagger}|\Phi_{i^{\prime}}\rangle{f^{\prime}}_{i^{\prime}\kappa^{\prime}}^{J\pi\alpha}=0.\\\
\end{aligned}\right.$ (19)

The coefficients, ${f^{\prime}}_{i^{\prime}(\kappa^{\prime})}^{J\pi\alpha}$,
should also satisfy the normalization condition, so that $\langle\Psi_{J\pi
M\alpha}^{\prime}(K)|\Psi_{J\pi M\alpha}^{\prime}(K)\rangle=1$. Again, we
assume $E^{\prime J\pi}_{1}\leq E^{\prime J\pi}_{2}\leq\cdots\leq E^{\prime
J\pi}_{m}$. Notice that all the matrix elements appearing in Eq. (19) are
already available in our VAP calculations, because they are needed in the
evaluation of the gradient and the Hessian matrix of $S_{m}$.

##  III 1p-1h mixing with the yrast states

Let us first study the simplest case of Eq. (11). The energy corresponding to
Eq. (11) can be expressed as

$E=\frac{\langle\Phi|\hat{H}P^{J\pi}_{KK}|\Phi\rangle}{\langle\Phi|P^{J\pi}_{KK}|\Phi\rangle}.$
(20)

Clearly, $E$ is a functional of $|\Phi\rangle$ governed by the complex
$d_{\kappa}$ parameters through Eq. (4). Here, $d_{\kappa}$ can be explicitly
written as

$d_{\kappa}=x_{\kappa}+iy_{\kappa},$ (21)

where, $x_{\kappa}$ and $y_{\kappa}$ are real. We assume $E$ reaches a minimum
at $d=0$, then the gradient of $E$ at this point should be zero. Namely,

$\displaystyle\left.\frac{\partial E}{\partial x_{\kappa}}\right|_{\vec{d}=0}$
$\displaystyle=$
$\displaystyle\left\\{\langle\Phi|A_{\kappa}(\hat{H}-E)P^{J\pi}_{KK}|\Phi\rangle\right.$
$\displaystyle\left.+\langle\Phi|(\hat{H}-E)P^{J\pi}_{KK}A^{\dagger}_{\kappa}|\Phi\rangle\right\\}\frac{1}{\langle\Phi|P^{J\pi}_{KK}|\Phi\rangle}$
$\displaystyle=$ $\displaystyle 0;$ $\displaystyle\left.\frac{\partial
E}{\partial y_{\kappa}}\right|_{\vec{d}=0}$ $\displaystyle=$ $\displaystyle
i\left\\{-\langle\Phi|A_{\kappa}(\hat{H}-E)P^{J\pi}_{KK}|\Phi\rangle\right.$
$\displaystyle\left.+\langle\Phi|(\hat{H}-E)P^{J\pi}_{KK}A^{\dagger}_{\kappa}|\Phi\rangle\right\\}\frac{1}{\langle\Phi|P^{J\pi}_{KK}|\Phi\rangle}$
$\displaystyle=$ $\displaystyle 0.$ (23)

This immediately leads to the following equation

$\displaystyle\langle\Phi|(\hat{H}-E)P^{J\pi}_{KK}A^{\dagger}_{\kappa}|\Phi\rangle=0.$
(24)

On the other side, it is known that in the Hartree-Fock theory, if the HF
energy

$E_{HF}=\frac{\langle\Phi|\hat{H}|\Phi\rangle}{\langle\Phi|\Phi\rangle}$ (25)

reaches a minimum, then one can easily see that no 1p-1h state on top of
$|\Phi\rangle$ can be mixed into $|\Phi\rangle$, which means

$\langle\Phi|\hat{H}A^{\dagger}_{\kappa}|\Phi\rangle=0.$ (26)

This inspired us that similar conclusion in the present VAP could be true. But
here the projection operator is involved. The projected 1p-1h state

$\displaystyle|\Psi_{\kappa}\rangle\equiv
P^{J\pi}_{MK}A^{\dagger}_{\kappa}|\Phi\rangle$ (27)

is no longer orthogonal to $|\Psi_{J\pi M}(K)\rangle$. One can do the Gram-
Schmid orthogonalization to get a new state

$\displaystyle|\Psi^{\prime}_{\kappa}\rangle=|\Psi_{\kappa}\rangle-\langle\Psi_{J\pi
M}(K)|\Psi_{\kappa}\rangle|\Psi_{J\pi M}(K)\rangle,$ (28)

so that

$\displaystyle\langle\Psi_{J\pi M}(K)|\Psi^{\prime}_{\kappa}\rangle=0.$ (29)

If there is no mixing of the 1p-1h projected state $|\Psi_{\kappa}\rangle$
with $|\Psi_{J\pi M}(K)\rangle$, then one should expect that

$\displaystyle\langle\Psi_{J\pi
M}(K)|\hat{H}|\Psi^{\prime}_{\kappa}\rangle=0.$ (30)

It is easy to prove that Eq. (30) is exactly equivalent to Eq. (24). This
clearly tells us that the property of no 1p-1h mixing at an energy minimum can
be generalized from HF theory to VAP.

Figure 2: (Color online) The converged $E^{J+}_{1}$ energies (black circles
for $K=0$ and red diamond for $K=5$) and the corresponding $E^{\prime J+}_{1}$
energies (blak bars for $K=0$ and red bars for $K=5$ ) as functions of the
number of reference SDs, $n$, for the yrast states with (a) $J^{\pi}$=8+ and
(b) $J^{\pi}$=0+ in 26Mg. The horizonal lines show the full SM eneriges. The
blue bar shows the converged energy of Eq. (1) taken from Fig. 1.

When more reference SDs are involved, as shown in Eq. (12), the situation
becomes more complicated. Since the included SDs $|\Phi_{i}\rangle$ and
$|\Phi_{i^{\prime}}\rangle$ are independent, and the associated matrix
elements
$\langle\Phi_{i}|(\hat{H}-{E^{\prime}}^{J\pi}_{\alpha})P_{KK}^{J\pi}A_{i^{\prime}\kappa^{\prime}}^{\dagger}|\Phi_{i^{\prime}}\rangle$
in Eq. (19) could be arbitrary. One may expect the 1p-1h mixing may take some
effect. To check such effect, the $E_{1}^{J\pi}$ energy and its corresponding
${E^{\prime}}_{1}^{J\pi}$ have been calculated with $n>1$. The converged
$E_{1}^{J\pi}$ energy and its corresponding ${E^{\prime}}_{1}^{J\pi}$ as
functions of $n$ are shown in Fig. 2. For the $8^{+}$ state in Fig. 2(a), the
calculations are performed with $K=0$ and $K=5$, respectively, to show the $K$
independence. It is clearly seen that the $E_{1}^{J\pi}$ energies with $K=0$
and $K=5$ are almost coincide and they become closer and closer to the
corresponding shell model energy, $e_{1}^{J\pi}$. However, it is striking that
all the energy differences $\delta E^{J\pi}_{1}$ are still numerically zero
for all calculated $n$ numbers. This implies that there may exist a theorem
for this interesting phenomenon. According to the above calculations, it seems
that, the situation with $E_{1}^{J\pi}={E^{\prime}}_{1}^{J\pi}$ only appears
after the VAP iteration converges, while the gradient of the energy becomes
zero.

Let us try to dig out this theorem. For convenience, we rewrite the VAP wave
function of Eq. (12) for the yrast state in a simpler form,

$|\Psi_{J\pi M1}(K)\rangle=\sum_{i=1}^{n}f_{i}P_{MK}^{J\pi}|\Phi_{i}\rangle.$
(31)

The corresponding energy, $E$, can be written as,

$\displaystyle
E=\frac{\sum_{ii^{\prime}=1}^{n}f_{i}^{*}{f}_{i^{\prime}}\langle\Phi_{i}|\hat{H}P^{J\pi}_{KK}|\Phi_{i^{\prime}}\rangle}{\sum_{ii^{\prime}=1}^{n}f_{i}^{*}{f}_{i^{\prime}}\langle\Phi_{i}|P^{J\pi}_{KK}|\Phi_{i^{\prime}}\rangle}.$
(32)

Clearly, $E$ is a function of all involved variational parameters
$d_{j\kappa}$, which are complex numbers,

$d_{j\kappa}=x_{j\kappa}+iy_{j\kappa},$ (33)

where $x_{j\kappa}$ and $y_{j\kappa}$ are real numbers. Notice that the
$d_{j}$ matrix only determines the $j$-th reference SD state
$|\Phi_{j}\rangle$. Similar to the deductions in Ref. VAP17 , the partial
derivatives $\frac{\partial E}{\partial x_{j\kappa}}$ and $\frac{\partial
E}{\partial y_{j\kappa}}$ at $\vec{d}=0$ can be expressed as

$\displaystyle\left.\frac{\partial E}{\partial
x_{j\kappa}}\right|_{\vec{d}=0}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}f_{j}^{*}\langle\Phi_{j}|A_{j\kappa}(\hat{H}-E)P^{J\pi}_{KK}|\Phi_{i}\rangle
f_{i}$ (34) $\displaystyle+$
$\displaystyle\sum_{i=1}^{n}f_{i}^{*}\langle\Phi_{i}|(\hat{H}-E)P^{J\pi}_{KK}A^{\dagger}_{j\kappa}|\Phi_{j}\rangle
f_{j},$ $\displaystyle\left.\frac{\partial E}{\partial
y_{j\kappa}}\right|_{\vec{d}=0}$ $\displaystyle=$
$\displaystyle-i\sum_{i=1}^{n}f_{j}^{*}\langle\Phi_{j}|A_{j\kappa}(\hat{H}-E)P^{J\pi}_{KK}|\Phi_{i}\rangle
f_{i}$ (35) $\displaystyle+$ $\displaystyle
i\sum_{i=1}^{n}f_{i}^{*}\langle\Phi_{i}|(\hat{H}-E)P^{J\pi}_{KK}A^{\dagger}_{j\kappa}|\Phi_{j}\rangle
f_{j}.$

Assuming the energy $E$ has reached a minimum at $\vec{d}=0$, the partial
derivatives in Eqs. (34) and (35) should be zero. Then one can get

$\displaystyle\sum_{i=1}^{n}f_{i}^{*}\langle\Phi_{i}|(\hat{H}-E)P^{J\pi}_{KK}A^{\dagger}_{j\kappa}|\Phi_{j}\rangle=0,$
(36)

and its Hermitian conjugate,

$\displaystyle\sum_{i=1}^{n}\langle\Phi_{j}|A_{j\kappa}(\hat{H}-E)P^{J\pi}_{KK}|\Phi_{i}\rangle
f_{i}=0.$ (37)

Now, let us come back to the 1p-1h mixing. For a 1p-1h projected state,

$|\Psi_{j\kappa}\rangle=P^{J\pi}_{KK}A_{j\kappa}^{\dagger}|\Phi_{j}\rangle,$
(38)

it may not be orthogonal to the converged VAP wave function $|\Psi_{J\pi
M}(K)\rangle$ in Eq. (31). So one may need to get a new state by performing
the Gram-Schmid orthogonalization,

$\displaystyle|\Psi^{\prime}_{j\kappa}\rangle=|\Psi_{j\kappa}\rangle-\langle\Psi_{J\pi
M1}(K)|\Psi_{j\kappa}\rangle|\Psi_{J\pi M1}(K)\rangle,$ (39)

so that

$\displaystyle\langle\Psi_{J\pi M1}(K)|\Psi^{\prime}_{j\kappa}\rangle=0.$ (40)

If there is no mixing between $|\Psi^{\prime}_{j\kappa}\rangle$ and
$|\Psi_{J\pi M1}(K)\rangle$, then one should have

$\displaystyle\langle\Psi_{J\pi
M1}(K)|\hat{H}|\Psi^{\prime}_{j\kappa}\rangle=0.$ (41)

Therefore, we have

$\displaystyle\langle\Psi_{J\pi M1}(K)|(\hat{H}-E)|\Psi_{j\kappa}\rangle=0,$
(42)

which is exactly the same as the Eq. (36). Thus, we have analytically proved
that 1p-1h mixing for the yrast state does not improve the converged VAP wave
function even with more reference SDs.

##  IV 1p-1h mixing with non-yrast states

As a further exploration, let us study the 1p-1h mixing with the non-yrast
states. Here, the sum of the lowest energies, $S_{m}(m>1)$, with the same
quantum numbers is minimized, as proposed in Ref. VAP18 . When the gradient of
$S_{m}$ becomes zero, the VAP calculation converges. However, the gradient of
$S_{m}$ being zero does not guarantee that the gradients of its members,
$E^{J\pi}_{\alpha}$, are zero. Thus, it is expected the 1p-1h mixing might
play some role in lowering the VAP energies and improving the VAP wave
functions in this case.

To check if such 1p-1h mixing is valid or not, some $fp$ shell nuclei are
calculated using the GXPF1A interaction gxpf1a . First of all, the sum of the
lowest five $J^{\pi}=9^{+}$ energies in 48Cr and the sum of the lowest five
$J^{\pi}=\frac{15}{2}^{-}$ energies in 49Cr are minimized, respectively. The
trial wave function in Eq. (11) is taken. The $K$ values are chosen to be $0$
for 48Cr and $1/2$ for 49Cr. For the number of selected reference SDs, we
chose $n=5,10,15$ and $20$, respectively, so that one can see the changes of
$E_{\alpha}^{J\pi}$ and ${E^{\prime}}_{\alpha}^{J\pi}$ as functions of $n$.
The calculated results are shown in Fig. 3. As expected, all the
$E_{\alpha}^{J\pi}$ energies gradually drop toward the corresponding shell
model ones as $n$ increases. However, different from the results in Fig. 2,
all the calculated ${E^{\prime}}_{\alpha}^{J\pi}$ energies are lower than the
corresponding $E_{\alpha}^{J\pi}$ values. This means the 1p-1h mixing takes
effect and indeed improves the VAP wave functions due to the nonzero gradients
of the calculated VAP energies, ${E}_{\alpha}^{J\pi}$.

The energy differences, $E_{\alpha}^{J\pi}-e_{\alpha}^{J\pi}$ and
${E^{\prime}}_{\alpha}^{J\pi}-e_{\alpha}^{J\pi}$ as functions of $n$ are
plotted in Fig. 4. Obviously, the values of
${E^{\prime}}_{\alpha}^{J\pi}-e_{\alpha}^{J\pi}$ are considerably smaller than
the ones of $E_{\alpha}^{J\pi}-e_{\alpha}^{J\pi}$. All the ratios of
${{E^{\prime}}_{\alpha}^{J\pi}-e_{\alpha}^{J\pi}}$ to
${E_{\alpha}^{J\pi}-e_{\alpha}^{J\pi}}$ are blow 0.8. Actually, most of them
are less than 0.6. This clearly shows the 1p-1h mixing can play an important
role in improving a group of nuclear wave functions with the same $J^{\pi}$
who are varied simultaneously.

Moreover, for all calculated states in Fig. 3, our results show that all the
energies with 1p-1h mixing, ${E^{\prime}}_{\alpha}^{J\pi}$, at $n=10$ are even
lower than the VAP energies, $E_{\alpha}^{J\pi}$, at $n=20$. We should mention
that the computational time for the VAP calculation with $n=20$ is about 3
times longer than the one with $n=10$. Since all the matrix elements in Eq.
(19) have already been evaluated during the VAP calculation, the only time
spent on the 1p-1h mixing is solving the expanded HW equations of Eq. (19).
Thus the 1p-1h mixing may not only improve the approximation of the VAP wave
functions but also substantially save the computational time. However, the
price is that the number of the included projected basis states is much
larger.

Figure 3: (Color online) The converged VAP energies, $E^{J\pi}_{\alpha}$,
(black circles) and their corresponding $E^{\prime J\pi}_{\alpha}$ ones (red
bars) calculated with different $n$ numbers for the lowest five states with
(a)$J^{\pi}=9^{+}$ in 48Cr and (b) $J^{\pi}=\frac{15}{2}^{-}$ in 49Cr. The SM
energies $e^{J\pi}_{\alpha}$ are shown as black bars. The GXPF1A interaction
is adopted.  Figure 4: (Color online) The energy differences,
$E_{\alpha}^{J\pi}-e_{\alpha}^{J\pi}$ (black squares) and
${E^{\prime}}_{\alpha}^{J\pi}-e_{\alpha}^{J\pi}$ (red dots), with different
$n$ numbers. The values of $E_{\alpha}^{J\pi}$, ${E^{\prime}}_{\alpha}^{J\pi}$
and $e_{\alpha}^{J\pi}$ are shown in Fig. 3.  Figure 5: (Color online) The
lowest five VAP energies, $E_{\alpha}^{J\pi}$ (blue circles) and the
corresponding ${E^{\prime}}_{\alpha}^{J\pi}$ (red bars) in 56Ni with $J^{\pi}$
running from $0^{+}$ to $12^{+}$. The lowest five VAP energies are obtained by
minimizing the sum of them at each $J^{\pi}$. The SM energies (black dots) are
shown for comparison.

Finally, the VAP plus 1p-1h mixing calculations for the states with spin
ranging from $J=0$ to $12$ in 56Ni are performed. For each spin, we still
minimize the sum of the lowest five energies with $n=10$. The GXPF1A
interaction is again adopted. The calculated results are shown in Fig. 5. It
is seen that with only 10 projected SDs for each spin, the VAP method is still
able to reproduce the schematic pattern of the low-lying energy levels of the
full shell model. The 1p-1h mixing again makes more or less improvements for
the calculated states. One may see that for some states far from the yrast
line, such improvements seem more apparent.

##  V Summary and outlook

Particle-hole mixing is a natural way of improving the quality of the nuclear
wave functions, as has been adopted in the random-phase approximation (RPA)
PR80 . However, the 2p-2h mixing involves too large matrices that is not our
present interest. In the present work, we only consider the possible 1p-1h
mixing with the VAP wave functions. It is well known that in the Hartree-Fock
method, the 1p-1h excited states can not be mixed into the converged HF vacuum
state. Very similarly, in our VAP calculation, when only the lowest energy
with a given spin and parity is minimized, we recognized that the 1p-1h
projected states on top of the projected reference SDs can not be mixed into
the converged VAP wave function. This interesting phenomenon has been
analytically proved in the present work and can be considered as a natural
generalization from the Hartree-Fock method. The common point is that for both
HF and VAP, the gradients of their converged energies are zero. However, after
one minimizes the sum of a group of lowest energies with the same quantum
numbers in VAP, the later 1p-1h mixing may indeed give the energies lower than
the corresponding VAP ones. This is because the final gradient of each energy
is not guaranteed to be zero. With these results, one can easily imagine such
1p-1h mixing is indeed somewhat limited, especially for those states which
energy gradients are relatively small. Of course, one may expect that the
approximation of the present work could be further improved by including a
number of the most important 2p-2h projected basis states. Such work will be
done in the future.

On the other side, it is well known that the pairing correlation is important
in the nuclear system. Such correlation can be incorporated by introducing the
HFB transformation. However, in the present method, we adopt the projected HF
SDs rather than the projected HFB vacuum states in order to save the
computational cost. Nevertheless, it is interesting that the obtained VAP
energies in the present examples are still very close to the exact shell model
ones even without the explicit HFB transformation. This might because the
adopted model space is rather small. However, in large model spaces, the BCS
correlation may be more prominent. This could make the approximation of our
method not as good as the present examples. One may think the mixing of the
1p-1h projected states not enough since the BCS correlation is non-
perturbative. In this situation, we consider three possible ways to improve
the approximation of the VAP wave functions. The first way is further
considering the mixing of the most important $n$p-$n$h projected SDs on top of
the VAP wave functions, which is very similar to our previous work PCI10 . The
second way is to add more reference SDs (i.e. $n$ in Eq. (12) is large) into
the VAP wave function, this certainly requires more computational cost. The
third way is directly taking the projected HFB vacuum states such as the
VAMPIR VAMPIR04  or the QVSM Shimizu21 , but this definitely increases the
computational time at least by two orders of magnitude due to the extra
projections for the numbers of both proton and neutron. It would be
interesting to compare these different methods in heavy deformed nuclei if all
such calculations can be implemented in large model spaces in the future.

###### Acknowledgements.

This work is supported by the National Natural Science Foundation of China
under Grant Nos. $11975314$,11575290, by the Key Laboratory of Nuclear Data
foundation (JCKY2022201C158) and by the Continuous Basic Scientific Research
Project Nos. WDJC-$2019$-$13$, BJ$20002501$.

##  Appendix A Basic matrix elements for the 1p-1h mixing

Matrix elements necessary for the 1p-1h mixing are presented as a complement
of those provided in Ref. lian22 , so that all basic matrix elements required
in VAP calculations are available. Following the notations in Ref. lian22 ,
let’s assume $|\Phi_{a}\rangle$ and $|\Phi_{b}\rangle$ as Slater determinants
of an $N-$particle system, i.e.,

$\displaystyle|\Phi_{a}\rangle$ $\displaystyle=$
$\displaystyle\hat{a}^{\dagger}_{1}\cdots\hat{a}^{\dagger}_{N}|\rangle,$ (43)
$\displaystyle|\Phi_{b}\rangle$ $\displaystyle=$
$\displaystyle\hat{b}^{\dagger}_{1}\cdots\hat{b}^{\dagger}_{N}|\rangle.$ (44)

According to the generalized Wick theorem, the overlap between
$|\Phi_{a}\rangle$ and $|\Phi_{b}\rangle$ can be written as,

$\displaystyle\langle\Phi_{a}|\Phi_{b}\rangle=\langle|\hat{a}_{N}\cdots\hat{a}_{1}\hat{b}^{\dagger}_{1}\cdots\hat{b}^{\dagger}_{N}|\rangle=\det(R),$
(45)

where $R$ is an $N\times N$ matrix with entries

$\displaystyle R_{ij}={\langle|\hat{a}_{i}\hat{b}^{\dagger}_{j}|\rangle}.$
(46)

In VAP calculations, we need to calculate the following explicit matrix
elements for the energy, the corresponding gradient and Hessian matrix,

$\displaystyle\langle\Phi_{a}|\hat{O}|\Phi_{b}\rangle,\,\langle\Phi_{a}|\hat{O}\hat{b}_{\mu}\hat{b}^{\dagger}_{\nu}|\Phi_{b}\rangle,\,\langle\Phi_{a}|\hat{a}_{\nu}\hat{a}^{\dagger}_{\mu}\hat{O}|\Phi_{b}\rangle,$
$\displaystyle\langle\Phi_{a}|\hat{a}_{\nu}\hat{a}^{\dagger}_{\mu}\hat{O}\hat{b}_{\mu^{\prime}}\hat{b}^{\dagger}_{\nu^{\prime}}|\Phi_{b}\rangle,\,\langle\Phi_{a}|\hat{O}\hat{b}_{\mu}\hat{b}^{\dagger}_{\nu}\hat{b}_{\mu^{\prime}}\hat{b}^{\dagger}_{\nu^{\prime}}|\Phi_{b}\rangle,$
$\displaystyle\langle\Phi_{a}|\hat{a}_{\nu^{\prime}}\hat{a}^{\dagger}_{\mu^{\prime}}\hat{a}_{\nu}\hat{a}^{\dagger}_{\mu}\hat{O}|\Phi_{b}\rangle,$
(47)

where,$1\leq\mu(\mu^{\prime})\leq N$ and $N+1\leq\nu(\nu^{\prime})\leq M$.
Here, $M$ is the dimension of the model space. $\hat{O}$ refers to a one-body
$\hat{T}$ or two-body operator $\hat{V}$, which can be written as

$\displaystyle\hat{T}$ $\displaystyle=$
$\displaystyle\sum_{\alpha\gamma}T_{\alpha\gamma}\hat{c}^{\dagger}_{\alpha}\hat{c}_{\gamma},$
(48) $\displaystyle\hat{V}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\sum_{\alpha\beta\delta\gamma}V_{\alpha\beta\gamma\delta}\hat{c}^{\dagger}_{\alpha}\hat{c}^{\dagger}_{\beta}\hat{c}_{\delta}\hat{c}_{\gamma}.$
(49)

In Ref. lian22 , we have presented the formulas for the following matrix
elements:

$\displaystyle\langle\Phi_{a}|\hat{O}|\Phi_{b}\rangle,\,\langle\Phi_{a}|\hat{a}_{\nu}\hat{a}^{\dagger}_{\mu}\hat{O}|\Phi_{b}\rangle,$
$\displaystyle\langle\Phi_{a}|\hat{a}_{\nu^{\prime}}\hat{a}^{\dagger}_{\mu^{\prime}}\hat{a}_{\nu}\hat{a}^{\dagger}_{\mu}\hat{O}|\Phi_{b}\rangle,$

which are necessary for the construction of the energy variance. Then
$\langle\Phi_{a}|\hat{O}\hat{b}_{\mu}\hat{b}^{\dagger}_{\nu}|\Phi_{b}\rangle$
and
$\langle\Phi_{a}|\hat{O}\hat{b}_{\mu}\hat{b}^{\dagger}_{\nu}\hat{b}_{\mu^{\prime}}\hat{b}^{\dagger}_{\nu^{\prime}}\hat{|}\Phi_{b}\rangle$
can be obtained through the following relations,

$\displaystyle\langle\Phi_{a}|\hat{O}\hat{b}_{\mu}\hat{b}^{\dagger}_{\nu}|\Phi_{b}\rangle$
$\displaystyle=$
$\displaystyle\langle\Phi_{b}|\hat{b}_{\nu}\hat{b}^{\dagger}_{\mu}\hat{O}|\Phi_{a}\rangle^{*},$
(50)
$\displaystyle\langle\Phi_{a}|\hat{O}\hat{b}_{\mu}\hat{b}^{\dagger}_{\nu}\hat{b}_{\mu^{\prime}}\hat{b}^{\dagger}_{\nu^{\prime}}\hat{|}\Phi_{b}\rangle$
$\displaystyle=$
$\displaystyle\langle\Phi_{b}|\hat{b}_{\nu^{\prime}}\hat{b}^{\dagger}_{\mu^{\prime}}\hat{b}_{\nu}\hat{b}^{\dagger}_{\mu}\hat{O}|\Phi_{a}\rangle^{*}.$
(51)

Therefore the only left matrix elements in Eq. (47) are the type of
$\langle\Phi_{a}|\hat{a}_{\nu}\hat{a}^{\dagger}_{\mu}\hat{O}\hat{b}_{\mu^{\prime}}\hat{b}^{\dagger}_{\nu^{\prime}}|\Phi_{b}\rangle$
which will be reused in the 1p-1h mixing in addition to those
$\langle\Phi_{a}|\hat{O}|\Phi_{b}\rangle$,
$\langle\Phi_{a}|\hat{O}\hat{b}_{\mu}\hat{b}^{\dagger}_{\nu}|\Phi_{b}\rangle$,
and
$\langle\Phi_{a}|\hat{a}_{\nu}\hat{a}^{\dagger}_{\mu}\hat{O}|\Phi_{b}\rangle$
ones. Thus, we only present the formulas for the
$\langle\Phi_{a}|\hat{a}_{\nu}\hat{a}^{\dagger}_{\mu}\hat{O}\hat{b}_{\mu^{\prime}}\hat{b}^{\dagger}_{\nu^{\prime}}|\Phi_{b}\rangle$
matrix elements. Based on the formulation in Ref. lian22 , one can get

$\displaystyle\langle\Phi_{a}|\hat{a}_{\nu}\hat{a}_{\mu}^{\dagger}\hat{T}\hat{b}_{\mu^{\prime}}\hat{b}_{\nu^{\prime}}^{\dagger}|\Phi_{b}\rangle$
(52) $\displaystyle=$
$\displaystyle\mathbb{T}_{\nu\nu^{\prime}}\overline{R}\\{\mu|\mu^{\prime}\\}+R_{\nu\nu^{\prime}}\overline{T}\\{\mu|\mu^{\prime}\\}$
$\displaystyle-\sum_{i,j}(\mathbb{T}_{\nu j}R_{i\nu^{\prime}}+R_{\nu
j}\mathbb{T}_{i\nu^{\prime}})\overline{R}\\{i\mu|j\mu^{\prime}\\}$
$\displaystyle-\sum_{i,j}R_{\nu
j}R_{i\nu^{\prime}}\overline{T}\\{i\mu|j\mu^{\prime}\\},$

where $\overline{R}\\{i|j\\}$, $\overline{R}\\{ij|kl\\}$ and
$\overline{R}\\{ijk|lmn\\}$ can be efficiently calculated using the formulas
in Ref. lian22 . $\mathbb{T}_{ij}$ is defined by

$\mathbb{T}_{ij}=\langle|\hat{a}_{i}\hat{T}\hat{b}^{\dagger}_{j}|\rangle=\sum_{\alpha\gamma}T_{\alpha\gamma}S^{+}_{i\alpha}S^{-}_{\gamma
j},$ (53)

where,
$S^{+}_{i\alpha}={\langle|\hat{a}_{i}\hat{c}^{\dagger}_{\alpha}|\rangle}$,
$S^{-}_{\gamma j}={\langle|\hat{c}_{\gamma}\hat{b}^{\dagger}_{j}|\rangle}$.
$\overline{T}\\{i|j\\}$ and $\overline{T}\\{ij|kl\\}$ can be written as,

$\displaystyle\overline{T}\\{i|j\\}$ $\displaystyle=$
$\displaystyle\sum_{i^{\prime}j^{\prime}}\mathbb{T}_{i^{\prime}j^{\prime}}\overline{R}\\{ii^{\prime}|jj^{\prime}\\},$
(54) $\displaystyle\overline{T}\\{ij|kl\\}$ $\displaystyle=$
$\displaystyle\sum_{i^{\prime}j^{\prime}}\mathbb{T}_{i^{\prime}j^{\prime}}\overline{R}\\{iji^{\prime}|klj^{\prime}\\}.$
(55)

For the two-body operator, we have

$\displaystyle\langle\Phi_{a}|\hat{a}_{\nu}\hat{a}_{\mu}^{\dagger}\hat{V}\hat{b}_{\mu^{\prime}}\hat{b}_{\nu^{\prime}}^{\dagger}|\Phi_{b}\rangle$
$\displaystyle=$ $\displaystyle
R_{\nu\nu^{\prime}}\overline{V}\\{\mu|\mu^{\prime}\\}-\sum_{ik}R_{\nu
k}R_{i\nu^{\prime}}\overline{V}\\{i\mu|k\mu^{\prime}\\}$
$\displaystyle+\sum_{ik}\mathbb{V}_{i\nu
k\nu^{\prime}}\overline{R}\\{i\mu|k\mu^{\prime}\\}$
$\displaystyle-\sum_{ijkl}[R_{\nu
l}\mathbb{V}_{ijk\nu^{\prime}}+R_{j\nu^{\prime}}\mathbb{V}_{i\nu
kl}]\overline{R}\\{ij\mu|kl\mu^{\prime}\\},$

with

$\displaystyle\mathbb{V}_{ijkl}$ $\displaystyle=$
$\displaystyle\sum_{\alpha\beta\gamma\delta}V_{\alpha\beta\gamma\delta}S_{i\alpha}^{+}S_{j\beta}^{+}S_{\gamma
k}^{-}S_{\delta l}^{-},$ (57) $\displaystyle\overline{V}\\{i|l\\}$
$\displaystyle=$
$\displaystyle\sum_{j<k,m<n}\mathbb{V}_{jkmn}\overline{R}\\{ijk|lmn\\},$ (58)
$\displaystyle\overline{V}\\{ij|lm\\}$ $\displaystyle=$
$\displaystyle\sum_{i^{\prime}<j^{\prime},l^{\prime}<m^{\prime}}\mathbb{V}_{i^{\prime}j^{\prime}l^{\prime}m^{\prime}}\overline{R}\\{iji^{\prime}j^{\prime}|lml^{\prime}m^{\prime}\\}.$

## References

  * (1) E. Caurier, G. Martínez-Pinedo, F. Nowacki, A. Poves, A. P. Zuker, Rev. Mod. Phys. 77 427 (2005).
  * (2) M. Horoi, B. A. Brown, and V. Zelevinsky, Phys. Rev. C 50, R2274 (1994).
  * (3) S. E. Koonin, D. J. Dean, and K. Langanke, Phys. Rep. 278, 1 (1997).
  * (4) T. Otsuka, M. Honma, T. Mizusaki, N. Shimizu, and Y. Utsuno, Prog. Part. Nucl. Phys. 47, 319 (2001).
  * (5) Z. C. Gao and M. Horoi, Phys. Rev. C 79, 014311 (2009).
  * (6) Z. C. Gao, Mihai Horoi and Y. S. Chen, Phys. Rev. C 80, 034325 (2009).
  * (7) K. W. Schmid, Prog. Part. Nucl. Phys. 52, 565 (2004).
  * (8) N. Shimizu, Y. Tsunoda, Y. Utsuno, and T. Otsuka, Phys. Rev. C 103, 014312 (2021).
  * (9) Z. C. Gao, M. Horoi, and Y. S. Chen, Phys. Rev. C 92, 064310 (2015).
  * (10) T. Ya, Y. He, Z. C. Gao, J. Q. Wang, and Y. S. Chen, Phys. Rev. C 95, 064307 (2017).
  * (11) J. Q. Wang, Z. C. Gao, Y. J. Ma, and Y. S. Chen, Phys. Rev. C 98, 021301(R) (2018).
  * (12) Z. C. Gao, Phys. Lett. B 824, 136795 (2022).
  * (13) K. Hara and Y. Sun, Int. J. Mod. Phys. E 4, 637 (1995).
  * (14) D. J. Thouless, Nucl. Phys. 21, 225 (1960).
  * (15) P. Ring and P. Schuck, The Nuclear Many-Body Problem, (Springer-Verlag, New York, 1980).
  * (16) J. Norcedal and S. J. Wright, Numerical Optimization (Springer Verlag, New York, Berlin, Heidelberg, 2006).
  * (17) E. A. Hylleraas, B. Undheim, Z. Phys. 65, 759 (1930).
  * (18) J. K. L. MacDonald, Phys. Rev. 43, 830 (1933).
  * (19) B. A. Brown and W. A. Richter, Phys. Rev. C 74, 034315 (2006).
  * (20) M. Honma, T. Otsuka, B.A. Brown, T. Mizusaki, Eur. Phys. J. A 25(1), 499 (2005).
  * (21) Z. J. Lian, X. Lu, and Z. C. Gao, Phys. Rev. C 106, 044308 (2022).