# Revisiting the variational two-particle reduced density matrix for nuclear
systems

J. G. Li  School of Physics, and State Key Laboratory of Nuclear Physics and
Technology, Peking University, Beijing 100871, China  N. Michel  Institute of
Modern Physics, Chinese Academy of Sciences, Lanzhou 730000, China  School of
Nuclear Science and Technology, University of Chinese Academy of Sciences,
Beijing 100049, China  W. Zuo  Institute of Modern Physics, Chinese Academy of
Sciences, Lanzhou 730000, China  School of Nuclear Science and Technology,
University of Chinese Academy of Sciences, Beijing 100049, China  F. R. Xu
<EMAIL_ADDRESS>School of Physics, and State Key Laboratory of Nuclear Physics
and Technology, Peking University, Beijing 100871, China

###### Abstract

In most nuclear many-body methods, observables are calculated using many-body
wave functions explicitly. The variational two-particle reduced density matrix
method is one of the few exceptions to the rule. Ground-state energies of both
closed-shell and open-shell nuclear systems can indeed be evaluated by
minimizing a constrained linear functional of the two-particle reduced density
matrix. However, it has virtually never been used in nuclear theory, because
nuclear ground states were found to be well overbound, contrary to those of
atoms and molecules. Consequently, we introduced new constraints in the
nuclear variational two-particle reduced density matrix method, developed
recently for atomic and molecular systems. Our calculations then show that
this approach can provide a proper description of nuclear systems where only
valence neutrons are included. For the nuclear systems where both neutrons and
protons are active, however, the energies obtained with the variational two-
particle reduced density matrix method are still overbound. The possible
reasons for the noticed discrepancies and solutions to this problem will be
discussed.

###### pacs:

##  I Introduction

In all areas of modern theoretical quantum mechanics, which comprises nuclear
physics, quantum chemistry, and condensed matter physics, many-body physics
consist in the study of correlated systems. They are induced by the two-body
interactions present between the fundamental particles building the many-body
system. The main problem arising when one aims at solving the many-body
problem in quantum mechanics is practical: the number of the occupied orbitals
at basis level increases linearly with system size, which causes the number of
configurations necessary to precisely represent the correlated many-body wave
function to increase exponentially.

In nuclear physics, most of the many-body methods are based on the use of an
explicit many-body wave function, which is a solution of the many-body
Schrödinger equation. Most popular methods are based on shell-model Caurier
_et al._ (2005); Otsuka  _et al._ (2020), variational quantum Monte Carlo
Carlson  _et al._ (2015) and valence space in-medium similarity
renormalization group methods Hergert  _et al._ (2016); Stroberg  _et al._
(2019). Other methods have been developed to recapture the features of the
many-body wave function with objects of reduced size, such as with the
coupled-cluster theory Hagen  _et al._ (2016), where one uses a similarity
transformed Hamiltonian acting on a simple reference state, for example a
Hartree-Fock state, using an exponentiated operator. However, the coupled-
cluster theory is typically limited to the use of singles-and-doubles
approximations Hagen  _et al._ (2016). Many-body wave functions are not
strictly necessary in order to calculate observables of many-body quantum
systems. Indeed, all one-body and two-body observables can be exactly
calculated using the two-particle reduced density matrix (2RDM) Mayer (1955);
Nakatsuji and Yasuda (1996); Coleman (1963); Mazziotti (2004); Hammond and
Mazziotti (2005); Mazziotti (2005, 2006a); Verstichel  _et al._ (2009); Nakata
and Yasuda (2009); Mazziotti (2012); Verstichel  _et al._ (2012); Mazziotti
(2016). In fact, the number of variables of the 2RDM functional scales only
polynomially with the number of basis orbitals, which is far less than the
exponential increase of basis configurations in configuration interaction
methods.

The 2RDM method is equivalent to solving the Schr$\ddot{\text{o}}$dinger
equation in order to obtain the many-body wave functions. The ground-state
energy is indeed a well-defined linear functional of the Hamiltonian matrix
and 2RDM elements in the variational 2RDM method. For a $N$-body system, the
Hamiltonian can be written as :

$\displaystyle\hat{H}=\sum_{ij}t_{ij}a_{i}^{\dagger}a_{j}+\frac{1}{2}\sum_{ijkl}V_{ij;kl}a_{i}^{\dagger}a_{j}^{\dagger}a_{l}a_{k}.$
(1)

And, the energy of the system can be expressed as Verstichel  _et al._ (2009):

$\displaystyle E(\Gamma)=\text{Tr}~{}\Gamma
H^{(2)}=\frac{1}{4}\sum_{ijkl}\Gamma_{ij;kl}H^{(2)}_{ij;kl},$ (2)

where the 2RDM elements appear :

$\displaystyle\Gamma_{ij;kl}=\langle\Psi^{N}|a_{i}^{\dagger}a_{j}^{\dagger}a_{l}a_{k}|\Psi^{N}\rangle,$
(3)

along with the reduced Hamiltonian two-body matrix elements :

$\displaystyle H_{ij;kl}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{1}{N-1}(\delta_{ik}t_{jl}-\delta_{il}t_{jk}-\delta_{jk}t_{il}+\delta_{jl}t_{ik})$
(4) $\displaystyle+$ $\displaystyle V_{ij;kl},$

where $t$ is the one-body kinetic energy operator and $V$ is the two-body
interaction, containing the Coulomb force for atomic systems and the
Coulomb+nuclear force for nuclear systems, and $N$ is the number of electrons
or nucleons.

2RDM was introduced by L$\ddot{\text{o}}$wdin, Mayer, and Coleman in the early
days of many-body quantum physics Löwdin (1955); Mayer (1955); Coleman (1963).
One hereby exploits the fact that a trial 2RDM can be considered as a variable
to obtain the ground-state energy without using many-body wave functions. The
2RDM can then be obtained from a minimization procedure, arising from the
variational principle, and hence being equivalent to solving the full many-
body Schrödinger equation Löwdin (1955). However, contrary to initial
expectations, the trial 2RDM, obtained by minimizing energy and subject only
to wave function normalization, is very different from the exact 2RDM arising
from many-body wave function eigenstates. In fact, even in weakly correlated
atoms and molecules, correlation energy is too large by orders of magnitude
Coleman (2000). It then became clear that the 2RDM obtained by energy
minimization must be constrained by additional rules (or conditions) in order
to reproduce the physical 2RDM. Coleman described these necessary and
sufficient conditions as $N\text{-}representability$ conditions Coleman
(1963), i.e., that the 2RDM functional is associated to a physical many-body
wave function. However, it was demonstrated that, in practice, demanding for
an 2RDM functional to fulfill exactly $N\text{-}representability$ conditions
is as difficult as solving the many-body Schrödinger equation Garrod and
Percus (1964). Approximate methods thus had to be devised for 2RDM-based
methods to be useful in practice.

Research on 2RDM methods has been systematically done for 50 years on atomic
and molecular electronic systems Coleman (1963); Hammond and Mazziotti (2005);
Zhao  _et al._ (2004), which allowed to identify the most important nontrivial
conditions needed to constrain the 2RDM functional. These conditions consist
in positive semidefinite operators, whose expectation value thus has to be
non-negative when evaluated with a 2RDM. Moreover, they can be systematically
arranged in a hierarchy where each level yields an increasingly tighter lower
bound on the exact ground-state energy Nakata and Yasuda (2009); Mazziotti
(2011, 2016, 2012); Anderson  _et al._ (2013); Nakata  _et al._ (2001).
However, limitations in optimization software and computer resources prevented
the practical development of the 2RDM method for several years. Significant
progress has been made over the last decade in 2RDM methods, as the possible
use of powerful computers allowed to devise numerical schemes otherwise
impossible to apply in practice Nakata  _et al._ (2001); Nakatsuji and Yasuda
(1996); Mazziotti (2004); Verstichel  _et al._ (2009); Nakata and Yasuda
(2009); Mazziotti (2012); Verstichel  _et al._ (2012); Mazziotti (2016, 2006b,
2020); Li and DePrince (2019); Maradzike  _et al._ (2017). Two approaches have
been developed to calculate the ground state 2RDM : (i) variational 2RDM with
energy minimization Mazziotti (2004); Verstichel  _et al._ (2012); Mazziotti
(2016) and (ii) non-variational approaches based on solutions of the
contracted Schr$\ddot{\text{o}}$dinger equation Nakatsuji and Yasuda (1996);
Mazziotti (2006b, 2020). In the present work, we focus on the variational 2RDM
method. The variational 2RDM minimization problem with positivity constraints
can be treated by a special class of optimization techniques known as
semidefinite programming Mazziotti (2004); Verstichel  _et al._ (2009); Nakata
and Yasuda (2009); Kamarchik and Mazziotti (2007); Mazziotti (2011). Note that
these optimization schemes also have applications in control theory Antipin
and Vasilieva (2015), combinatorial optimization Mahdavi and Shiri (2015) and
finance Zhang  _et al._ (2006).

The first 2RDM calculations of nuclear systems were devised in 1970s
Mihailović and Rosina (1975). However, the obtained ground-state energies were
typically overbound by several MeV, and even sometimes tens of MeV. It was
then clear that the variational 2RDM method could not compete with
configuration interaction methods effected in small model spaces, for example,
with which tremendous success was obtained in the description of $p$ Cohen and
Kurath (1965), $sd$ Halbert  _et al._ (1968) and $psd$ crossed-shell nuclei
Zuker  _et al._ (1968); Zuker (1969). Many works Assié and Lacroix (2009);
Tohyama and Schuck (2014, 2017, 2019); Barton  _et al._ (2021) have been tried
to implement the 2RDM technique using a truncated time-dependent Bogoliubov-
Born-Green-Kirkwood-Yvon approach together with an adiabatic method for the
ground state of nuclear systems. However, there has been no variational 2RDM
calculation of nuclear system for a long time. Nevertheless, shell-model
approaches have not solved the many-body problem in the nuclear chart. Indeed,
even with the most powerful machines, shell-model calculations in large model
spaces cannot be performed at present without model space truncations Caurier
_et al._ (2005). The situation is even more drastic for heavy nuclei, where
only mean-field based methods can be applied in systematic calculations Erler
_et al._ (2012); Kortelainen  _et al._ (2010); Stoitsov  _et al._ (2003);
Bender  _et al._ (2003); Kortelainen  _et al._ (2010); Bulgac  _et al._
(2018). Consequently, one can still find interest in the use of the
variational 2RDM method even in our era of intense numerical calculations.

The 2RDM method is a powerful many-body method for the description of atomic
and molecular systems Nakatsuji and Yasuda (1996); Mazziotti (2004);
Verstichel  _et al._ (2009); Nakata and Yasuda (2009); Mazziotti (2012);
Verstichel  _et al._ (2012); Mazziotti (2016, 2006b, 2020). The precision of
2RDM calculations is, in fact, comparable with coupled-cluster with singles-
and-doubles approximations for these weakly correlated systems. Furthermore,
it has been shown that variational 2RDM calculations can capture many-electron
excitations in the calculations of atomic systems Mazziotti (2004). The
variational 2RDM has also been successfully applied in the context of the
pairing model Alcoba1  _et al._ (2018). In the present paper, we will employ
the variational 2RDM for nuclear systems. The partial 3-particle $T_{1}$ and
$T_{2}^{\prime}$ conditions, are first considered in the variational 2RDM
calculations for the nuclear systems. Systematic calculations are performed
using variational 2RDM with different conditions for nuclear systems.

This paper is constructed as follows: we will introduce the variational 2RDM
method in Sec.(II) at theoretical level. The used N-representability
conditions, fundamental to determine the 2RDM by energy minimization, will be
described. The standard 2-positivity and 3-positivity N-representability
conditions developed for the study of molecular systems and which we will use
afterwards in nuclei will be presented. The computational method used for
constrained energy minimization, based on the semidefinite optimization
boundary point method, will be shortly dealt with. Results will be shown in
Sec.(III), where ground states of nuclear Hamiltonians based on cluster-
orbital shell-model (COSM) and standard shell-model (SM) frameworks are
considered. The energies obtained with both 2RDM and exact full configuration
interaction methods will be compared. One will see that the quality of results
largely differ according to the nuclear system considered. The influences of
the 2-positivity and 3-positivity N-representability conditions on nuclear
ground states will be emphasized for that matter. Plans for future use of the
variational 2RDM will then be depicted. Conclusion will be made afterwards.

##  II Method

In this section, we introduce the theoretical details underlying the
variational 2RDM framework. Standard N-representability conditions, as well as
the energy optimization procedure, belonging to the class of semidefinite
optimization boundary point methods, will be presented and extended to the
case of nuclear systems.

###  II.1 Variational 2RDM method

The energy of a many-body system is a linear functional of the reduced density
matrix. For a two-body Hamiltonian, the ground-state energy $E_{g}$ can be
written as

$\displaystyle
E_{g}=\text{Tr}[H~{}{{}^{2}\Gamma}]=\text{Tr}[H^{(2)}~{}{{}^{2}\Gamma}],$ (5)

where $H^{(2)}$ is the two-body reduced Hamiltonian and ${}^{2}\Gamma$ is two-
particle reduced density matrix. The definition of the $H^{(2)}$ and
${}^{2}\Gamma$ operators can be found in Eqs. (3) and (4). The only difference
from atomic and molecular systems is that one has both valence proton and
neutron occupied states Mihailović and Rosina (1975). In nuclear physics, we
typically treat the many-body problem by including the nuclear interaction as
a two-body operator (the three-body force is present by way of a normal-order
approach in realistic calculations Hagen  _et al._ (2007)). The 2RDM is then
sufficient to calculate all the observables related up to a one-body or two-
body operator.

An arbitrary 2RDM is, in general, not associated to a many-body wave function,
so that the use of an underconstrained 2RDM typically leads to overbound
ground states Coleman (2000); Verstichel (2012); Mihailović and Rosina (1975).
Conversely, except for the trivial case of non-interacting particles, the
minimization of the energy in Eq. (5) with respect to the 2RDM provides with
the exact ground-state energy if the energy minimization is effected by
constraining Hamiltonian expectation values over all possible
$N$-representable 2RDM:

$\displaystyle
E_{g}=\mathop{\text{min}}\limits_{{}^{2}\Gamma\in^{2}\mathcal{D_{N}}}\text{Tr}[H^{(2)}~{}{{}^{2}\Gamma}],$
(6)

where ${}^{2}\mathcal{D_{N}}$ contains all the 2RDM functionals generated by a
many-body wave function. In order to show that a given 2RDM belongs to
${}^{2}\mathcal{D_{N}}$, one would have to check that the expectation value of
all combinations of creation and annihilation operators is positive or equal
to zero. Clearly, this is as difficult as solving the many-body Schrödinger
equation, so that Eq. (6) is only of theoretical interest. Therefore, in
practice, one replaces ${}^{2}\mathcal{D_{N}}$ by the set of 2RDM functionals
verifying a subset of N-representability conditions, which consist in the
N-representability conditions, developed in molecular 2RDM variational
approaches Garrod and Percus (1964); Erdahl (1978); Zhao  _et al._ (2004).

###  II.2 N-representability conditions and many-body quantum numbers

In 1932, von Neumann showed that a general matrix is a density matrix only if
it fulfills the following conditions von Neumann (1995): the matrix must be
(i) Hermitian, (ii) normalized (fixed trace), (iii) antisymmetric (fermions)
or symmetric (bosons) under particle exchange, and (iv) positive semidefinite
to yield non-negative probabilities to find particles. Variational
minimization of the energy in Eq. (6) with respect to the 2RDM constrained by
von Neumann conditions thus provides with the exact ground-state energy. The
positive semidefinite character of $\Gamma$ is called the $P$ condition Garrod
and Percus (1964); Coleman (1963).

The 2RDM can be mapped onto the 1RDM through operator contractions:

$\displaystyle\rho^{(\nu)}_{ij}=\frac{1}{A-1}\sum_{k}\Gamma_{ik;jk},$ (7)

where $\nu$ stands for proton or neutron and $A$ is the number of nucleons.
Consequently, one obtains that the 1RDM is positive semidefinite, i.e.,
$\Gamma\geq 0\Rightarrow\rho^{(\nu)}\geq 0$. Moreover, the density
normalization implies that the trace of the proton (neutron) 1RDM must be
equal to the proton (neutron) number:

$\displaystyle\text{Tr}\rho^{(\nu)}=N_{\nu}.$ (8)

For simplicity, one will no longer state in the following whether a
creation/annihilation operator is related to proton or neutron states, even
though this must be effected in practice.

Additional ensembles of N-representability conditions can be obtained by
considering the positivity of operators related to the 2RDM. Different
representations of $p$RDMs can be easily defined using second quantization
Mazziotti (2012). For example, let us consider the metric (or overlap)
matrices $M$ defined by:

$\displaystyle
M_{ij}=\langle\Phi_{i}|\Phi_{j}\rangle=\langle\Psi|C_{i}C_{j}^{\dagger}|\Psi\rangle,$
(9)

from the set of basis functions

$\displaystyle\langle\Phi_{i}|=\langle\Psi|C_{i},$ (10)

where each $C_{i}$ is a product of $p$ creation and/or annihilation operator.
When $C_{i}$ are products of $p$ creation operators, the metric matrix in Eq.
(9) is the $p$RDM. Added to that, one can generate $p$ additional metric
matrices if the $C_{i}$ become products of $p$ second-quantized operators
having different numbers of creation and annihilation operators. These $p+1$
metric matrices must be positive semidefinite Mazziotti (2012); Coleman and
Yukalov (2000); Mazziotii (2007). The demanded positivity of all the
considered metric matrices of order $p$ is donated as the $p$-positivity
condition.

Two independent classes of matrices can be constructed to impose the
1-positivity condition, defined as

$\displaystyle\rho_{ij}=\langle\Phi_{i}^{\rho}|\Phi_{j}^{\rho}\rangle=\langle\psi|a_{i}^{\dagger}a_{j}|\psi\rangle,~{}~{}\rho\geq
0,$ (11) $\displaystyle
q_{ij}=\langle\Phi_{i}^{q}|\Phi_{j}^{q}\rangle=\langle\psi|a_{i}a_{j}^{\dagger}|\psi\rangle,~{}~{}q\geq
0,$ (12)

where $|\Phi_{j}^{\rho}\rangle$ and $|\Phi_{j}^{q}\rangle$ are $(N-1)$\- and
$(N+1)$\- particle basis function, respectively. $\rho_{ij}$ corresponds to
the 1RDM, whereas $q_{ij}$ can be expressed as a linear function of
$\rho_{ij}$ by using the anticommutation relations:

$\displaystyle q_{ij}=\delta_{ij}-\rho_{ij},$ (13)

The 1-positivity condition, embodied by the two metric matrices ($\rho_{ij}$
and $q_{ij}$), implies that the eigenvalues of a fermionic 1RDM are between 0
and 1.

Three independent metric matrices can also be generated from the operators
$C_{i}$ out of the 2-positivity condition. They consist of the $P$-,$Q$-, and
$G$\- matrices:

$\displaystyle
P_{ij;kl}=\langle\Phi_{ij}^{P}|\Phi_{kl}^{P}\rangle=\langle\psi|a_{i}^{\dagger}a_{j}^{\dagger}a_{l}a_{k}|\psi\rangle,~{}~{}P\geq
0,$ (14) $\displaystyle
Q_{ij;kl}=\langle\Phi_{ij}^{Q}|\Phi_{kl}^{Q}\rangle=\langle\psi|a_{i}a_{j}a_{l}^{\dagger}a_{k}^{\dagger}|\psi\rangle,~{}~{}Q\geq
0,$ (15) $\displaystyle
G_{ik;lj}=\langle\Phi_{ij}^{G}|\Phi_{kl}^{G}\rangle=\langle\psi|a_{i}^{\dagger}a_{k}a_{j}^{\dagger}a_{l}|\psi\rangle,~{}~{}G\geq
0,$ (16)

where $|\Phi_{kl}^{P}\rangle,|\Phi_{kl}^{Q}\rangle$, and
$|\Phi_{kl}^{G}\rangle$, are $(N-2)$-, $(N+2)$\- , and $N$\- particle basis
functions, respectively. The $P$ matrix is the 2RDM, and the $Q$ and $G$
matrices can be rewritten as a function of the 1RDM and 2RDM by reordering the
creation and annihilation operators using anticommutation relations. One then
obtains the following linear maps:

$\displaystyle Q_{ij;kl}$ $\displaystyle=$
$\displaystyle\delta_{ik}\delta_{jl}-\delta_{jk}\delta_{il}+P_{ij;kl}$ (17)
$\displaystyle-$
$\displaystyle(\delta_{ik}{\rho_{jl}}-\delta_{il}{\rho_{jk}}-\delta_{jk}{\rho_{il}}+\delta_{jl}{\rho_{ik}}),$
$\displaystyle G_{ij;kl}$ $\displaystyle=$
$\displaystyle\delta_{jl}{\rho_{ik}}-P_{il;kj}.$ (18)

The metric matrices associated to the 3-positivity condition read similarly to
Eqs. (14), (15), and (16):

${}^{3}D_{ijk;lmn}=\langle\Psi|a_{i}^{\dagger}a_{j}^{\dagger}a_{k}^{\dagger}a_{n}a_{m}a_{l}|\Psi\rangle,~{}~{}^{3}D\geq
0,$ (19)
${}^{3}E_{ijk;lmn}=\langle\Psi|a_{i}^{\dagger}a_{j}^{\dagger}a_{k}a_{n}^{\dagger}a_{m}a_{l}|\Psi\rangle,~{}~{}^{3}E\geq
0,$ (20)
${}^{3}F_{ijk;lmn}=\langle\Psi|a_{i}^{\dagger}a_{j}a_{k}a_{n}^{\dagger}a_{m}^{\dagger}a_{l}|\Psi\rangle,~{}~{}^{3}F\geq
0,$ (21)
${}^{3}Q_{ijk;lmn}=\langle\Psi|a_{i}a_{j}a_{k}a_{n}^{\dagger}a_{m}^{\dagger}a_{l}^{\dagger}|\Psi\rangle,~{}~{}^{3}Q\geq
0,$ (22)

where ${}^{3}D$ and ${}^{3}Q$ are the 3p-RDM and 3h-RDM, respectively, and
${}^{3}E$ and ${}^{3}F$ are generalizations of the $G$ matrix. In practice,
however, the full 3-positivity condition is too cumbersome to apply. As a
consequence, the two less stringent conditions $T_{1}$ and $T_{2}$, originally
proposed by Erdahl Erdahl (1978), and implemented for small atoms and
molecules by Zhao $et~{}al.$ Zhao  _et al._ (2004) are used in practical
calculations. They read:

$\displaystyle(T_{1})_{ijk;lmn}={{}^{3}}D_{ijk;lmn}+{{}^{3}}Q_{ijk;lmn},~{}~{}~{}T_{1}\geq
0,$ (23)
$\displaystyle(T_{2})_{ijk;lmn}={{}^{3}}E_{ijk;lmn}+{{}^{3}}F_{ijk;lmn},~{}~{}~{}T_{2}\geq
0,$ (24)

which are both positive semidefinite. The fundamental interest of Eqs. (23)
and (24) is that the effected summation of the considered 3-positive metric
matrices cause their connected (or cumulant) part to cancel exactly, so that
both the $T_{1}$ and the $T_{2}$ metric matrices can be rewritten as a
function of the 1RDM and 2RDM only:

$\displaystyle(T_{1})_{ijk;lmn}$ $\displaystyle=$
$\displaystyle\delta_{kn}\delta_{jm}\delta_{il}-\delta_{km}\delta_{il}\delta_{jn}+\delta_{in}\delta_{km}\delta_{jl}$
(25) $\displaystyle-$
$\displaystyle\delta_{kn}\delta_{im}\delta_{jl}+\delta_{jn}\delta_{im}\delta_{kl}-\delta_{in}\delta_{jm}\delta_{kl}$
$\displaystyle-$
$\displaystyle(\delta_{kn}\delta_{jm}-\delta_{jn}\delta_{kn}){\rho_{il}}+(\delta_{kn}\delta_{im}-\delta_{in}\delta_{km}){\rho_{jl}}$
$\displaystyle-$
$\displaystyle(\delta_{jn}\delta_{im}-\delta_{in}\delta_{jm}){\rho_{kl}}+(\delta_{km}\delta_{jl}-\delta_{jn}\delta_{kl}){\rho_{im}}$
$\displaystyle-$
$\displaystyle(\delta_{kn}\delta_{il}-\delta_{in}\delta_{kl}){\rho_{jm}}+(\delta_{jn}\delta_{il}-\delta_{in}\delta_{jl}){\rho_{km}}$
$\displaystyle-$
$\displaystyle(\delta_{jl}\delta_{km}-\delta_{jm}\delta_{kl}){\rho_{in}}+(\delta_{km}\delta_{il}-\delta_{im}\delta_{kl}){\rho_{jn}}$
$\displaystyle-$
$\displaystyle(\delta_{jm}\delta_{il}-\delta_{im}\delta_{jl}){\rho_{kn}}+\delta_{kn}{P_{ij;lm}}$
$\displaystyle-$
$\displaystyle\delta_{jn}{P_{ik;lm}}+\delta_{in}{P_{jk;lm}}-\delta_{km}{P_{ij;lm}}$
$\displaystyle-$
$\displaystyle\delta_{jm}{P_{ik;ln}}-\delta_{im}{P_{jk;ln}}+\delta_{kl}{P_{ij;mn}}$
$\displaystyle-$ $\displaystyle\delta_{jl}{P_{ik;mn}}+\delta_{il}{P_{jk;mn}},$
$\displaystyle(T_{2})_{ijk;lmn}$ $\displaystyle=$
$\displaystyle(\delta_{il}\delta_{jm}-\delta_{im}\delta_{jl}){\rho_{kn}}$ (26)
$\displaystyle-$ $\displaystyle\delta_{kn}{P_{ij;lm}}-\delta_{il}{P_{km;nj}}$
$\displaystyle-$
$\displaystyle\delta_{jl}{P_{km;ni}}+\delta_{im}{P_{kl;nj}}-\delta_{jm}{P_{kl;ni}}.$

The $T_{1}$ and $T_{2}$ conditions are, however, inferior to the full
3-positivity condition, since they imply the positive semidefinite character
of the expectation value of the sum
$C_{i}C_{j}^{\dagger}+C_{j}C_{i}^{\dagger}$ with a trial 2RDM, and not that of
its individual terms, equal to $C_{i}C_{j}^{\dagger}$ and
$C_{j}C_{i}^{\dagger}$. Usually, an enhancement of the $T_{2}$ condition,
denoted as the the $T_{2}^{\prime}$ condition, is used instead of the $T_{2}$
condition. The $T_{2}^{\prime}$ condition arises from the addition of a one-
particle operator in the metric defining $T_{2}$. In the $T_{2}^{\prime}$
condition, the positive semidefinite form is equal to:

$\displaystyle(C_{i}+B_{i})^{\dagger}(C_{j}+B_{j})+C_{j}C_{i}^{\dagger},$ (27)

where $B_{i}$ is a one-body operator. The positive semidefinite character
arising from Eq. (27) translates into a matrix positivity condition:

$\displaystyle
T_{2}^{\prime}=\left(\begin{array}[]{cc}(T_{2})_{ijk;lmn}&P_{ij;rk}\\\
P_{sn;lm}^{\dagger}&\rho_{sr}\end{array}\right)\geq 0.$ (30)

One can check that the use of $T_{2}$ or $T_{2}^{\prime}$ virtually has the
same numerical cost. For more theoretical and historical details, we refer the
reader to Ref. Mazziotii (2007).

Nuclear wave functions are subject to particle number conservation, which can
be expressed with the particle and pair number operators $\hat{N}_{\nu}$ and
$\hat{N}_{\nu\nu}$, respectively, reading

$\displaystyle\hat{N}_{\nu}$ $\displaystyle=$
$\displaystyle\sum_{\alpha_{\nu}}a^{\dagger}_{\alpha_{\nu}}a_{\alpha_{\nu}},$
(31) $\displaystyle\hat{N}_{\nu\nu}$ $\displaystyle=$
$\displaystyle\sum_{\alpha_{\nu}<\beta_{\nu}}a^{\dagger}_{\alpha_{\nu}}a^{\dagger}_{\beta_{\nu}}a_{\beta_{\nu}}a_{\alpha_{\nu}},$
(32)

where $\nu$ represents proton or neutron. Many-body systems obey rotational
symmetries as well, i.e. the total angular momentum operator $J$ and total
isospin operator $T$ (in the absence of the Coulomb force for the latter) are
conserved quantities. As the expectation values of the operators associated to
conservation laws are linear with respect to the 2RDM, they can be directly
included as constraints:

$\displaystyle\langle\Phi|\hat{N}_{\nu\nu}|\Phi\rangle$ $\displaystyle=$
$\displaystyle\text{Tr}~{}\Gamma_{\nu\nu}=\frac{N_{\nu}(N_{\nu}-1)}{2},$ (33)
$\displaystyle\langle\Phi|\hat{N}\hat{Z}|\Phi\rangle$ $\displaystyle=$
$\displaystyle\text{Tr}~{}\Gamma_{pn}=NZ,$ (34)
$\displaystyle\langle\Phi|\hat{J}^{2}|\Phi\rangle$ $\displaystyle=$
$\displaystyle\text{Tr}~{}\hat{J}^{2}\Gamma=J(J+1),$ (35)
$\displaystyle\langle\Phi|\hat{T}^{2}|\Phi\rangle$ $\displaystyle=$
$\displaystyle\text{Tr}~{}\hat{T}^{2}\Gamma=T(T+1),$ (36)

where $N_{\nu}=Z$ or $N_{\nu}=N$, and the isospin conservation is used only in
the absence of the Coulomb force. Note that proton and neutron number
conservation, embodied in Eq. (8), is a direct consequence of Eqs. (7), (33)
and (34). In particular, this implies that
$\langle\Phi|\hat{N}_{\nu}^{2}|\Phi\rangle=\langle\Phi|\hat{N}_{\nu}|\Phi\rangle^{2}$,
as $2\hat{N}_{\nu\nu}=\hat{N}_{\nu}^{2}-\hat{N}_{\nu}$ (see Eqs. (31) and
(32)).

Variational 2RDM calculations can be constrained using different conditions.
Calculations using the $P$, $Q$ and $G$ 2-positive conditions, combined with
the quantum number conservation conditions of Eqs. (33), (34), (35), and (36)
will be denoted by $PQG$. The inclusion of $T_{1}$ or $T_{2}^{\prime}$, or
$T_{1}$ and $T_{2}^{\prime}$ will be denoted by $PQGT_{1}$,
$PQGT_{2}^{\prime}$, and $PQGT_{1}T_{2}^{\prime}$, respectively.

Even though Eqs. (14), (15), (16), (23), and (30) make only use of 1RDM and
2RDM, they cannot be used directly in practice. This arises because their
computational cost increases quickly with the number of one-body basis states.
Indeed, for a number of one-body basis states equal to $s$, the memory needed
to store matrices using $PQG$ ($PQGT_{1}$, $PQGT_{2}$,
$PQGT_{1}T_{2}^{\prime}$) matrices is $\mathcal{O}(s^{4})$
($\mathcal{O}(s^{6})$) while the time of calculation scales as
$\mathcal{O}(s^{6})$ ($\mathcal{O}(s^{9})$), and is independent of the number
of valence nucleons. These estimates arise because up to four (six)
creation/annihilation operators are present in $PQG$ ($PQGT_{1}$, $PQGT_{2}$,
$PQGT_{1}T_{2}^{\prime}$) operators, on the one hand, and because one has to
diagonalize the $PQG$ ($PQGT_{1}$, $PQGT_{2}$, $PQGT_{1}T_{2}^{\prime}$)
matrices in the constrained energy minimization procedure (see Sec.(II.3)), on
the other hand, along with a linear system to solve of similar computational
cost Verstichel (2012).

Consequently, it is necessary in practice to couple angular momenta in Eqs.
(14), (15), (16), (23), and (30). Atomic and molecular systems are considered
in the $LS$ scheme as spin $S$ and orbital $L$ angular momenta can be
considered as good quantum numbers Verstichel (2012). Conversely, in nuclear
systems, only the total angular momentum $J$ is conserved, so that one uses
$jj$ coupling scheme, i.e., one works in $J$-scheme. Coupling algebra for that
matter is standard, as it makes use of the Wigner-Eckart theorem only, so that
the obtained equations are defined from the $J$-coupled 1RDM and 2RDM, Wigner
signs and three-body coefficients of fractional parentage Verstichel (2012).
The obtained reduction in the demanded memory and time of calculation is huge,
as the memory and time of calculation for the $PQG$ ($PQGT_{1}$, $PQGT_{2}$,
$PQGT_{1}T_{2}^{\prime}$) case is now $\mathcal{O}(r^{4})$
($\mathcal{O}(r^{6})$) and $\mathcal{O}(r^{6})$ ($\mathcal{O}(r^{9})$),
respectively, with $r$ the number of orbitals. Indeed, the ratio $s/r$, which
is the average number of states in one basis orbital, is around 5, so that the
gain in computational cost can easily reach 10,000-100,000 Verstichel (2012).
However, the use of $PQGT_{1}$, $PQGT_{2}$ or $PQGT_{1}T_{2}^{\prime}$
conditions still remains more expensive that of $PQG$ by a large factor, which
can reach 100-1000, as one can have tens of orbitals in practical calculations
Verstichel (2012). Consequently, it is customary to impose only $PQG$
conditions in atomic and molecular conditions Mazziotti (2011); Verstichel
(2012). However, we will see in the following that $T_{1}$ and
$T_{2}^{\prime}$ conditions can barely be ignored in nuclear systems (see
Sec.(III)).

###  II.3 Semidefinite optimization

In the variational reduced density matrix theory, the energy of Eq. (5) is
minimized with respect to the 2RDM subject to the constraints outlined above.
The solution of the $E_{g}$ problem is then a semidefinite linear optimization
problem with constraints. The primal formulation of this problem, i.e., where
the 2RDM is considered as a variable to optimize, can be expressed as:

$\displaystyle\text{minimize}~{}~{}E_{g}=\text{Tr}~{}H\Gamma,$
$\displaystyle\text{such~{} that}~{}~{}A\Gamma=X,$
$\displaystyle\text{with}~{}~{}X\geq 0,$ (37)

where $\Gamma$ is the primal solution vector and $H$ is the Hamiltonian which
contains the one-body and two-body matrix elements. The constraint matrix $A$
and constraint vector $X$ encode the $N$-$representability$ conditions and the
equalities of Eqs. (33), (34), (35), and (36) associated to conserved quantum
numbers that the $\Gamma$ matrix must satisfy. The operator $X$ maps the
primal solutions onto the set of positive semidefinite RDMs:

$\displaystyle X=\begin{pmatrix}~{}P&0&0&0&0&0\\\ 0&Q&0&0&0&0\\\
0&0&G&0&0&0\\\ 0&0&0&T_{1}&0&0\\\ 0&0&0&0&T_{2}^{\prime}&0\\\
0&0&0&0&0&C~{}\end{pmatrix},$ (38)

where $C$ represents the action of the operators associated to conserved
quantum numbers, as it can be written as $C\Gamma=0$ (see Eqs. (33), (34),
(35), and (36)). In order to find the $\Gamma$ matrix, one will apply the
boundary point positive semidefinite algorithm. The boundary point method is
actually an instance of a more general class of augmented Lagrangian
approaches for solving the positive semidefinite problem Mazziotti (2004,
2011); Verstichel  _et al._ (2011); Verstichel (2012). The standard Lagrangian
for the positive semidefinite problem can be written as,

$\displaystyle\mathcal{L}=\text{Tr}~{}H\Gamma+\Lambda\cdot(X-A\Gamma),$ (39)

to minimize the $H\Gamma$ functional and where a Lagrangian multiplier matrix
$\Lambda$ is introduced to account for constraints. While Eq. (39) would
provide the $\Gamma$ matrix from a theoretical point of view, it cannot be
used in practice because the convergence would be very poor, or even
nonexistent numerically Staszczak  _et al._ (2010). This arises because the
functional of Eq. (39) is not convex, i.e., it does not resemble a paraboloid
in the vicinity of the sought minimum Staszczak  _et al._ (2010).

In order to account for this problem, the augmented Lagrangian method for
primal problems has been developed. For this, one adds a quadratic penalty for
the so-called infeasibility, i.e., the fact that $\Gamma$ does not belong to
the set of acceptable solutions during optimization because $A\Gamma\neq X$
Zhao  _et al._ (2004):

$\displaystyle\mathcal{L}=\text{Tr}~{}H\Gamma+\Lambda\cdot(X-A\Gamma)+\frac{\sigma}{2}||X-A\Gamma||^{2},$
(40)

where a penalty parameter, $\sigma>0$, is introduced and determines how
stringently the $N$-$representability$ constraint is enforced. As Eq. (40) is
quadratic in $\Gamma$, its minimum is straightforward to find with standard
minimization methods if $X$ is kept fixed Verstichel (2012). $X$ is then
obtained as a second step also from the minimization of Eq. (40), but by
keeping $\Gamma$ fixed Verstichel (2012). In particular, $X$ positivity is
enforced by way of a diagonalization procedure and matrix manipulations, where
the negative eigenvalues of the $P$, $Q$, $G$, $T_{1}$, $T_{2}^{\prime}$
metric matrices and the non-zero eigenvalues of the $C$ matrix in Eq. (38) are
replaced by zero Verstichel (2012); Mazziotti (2011). $\sigma$ is modified at
each iteration according to the values of $\nabla\mathcal{L}$ and
$||X-A\Gamma||$ to accelerate convergence Burer and Monteiro (2003);
Verstichel (2012). One then obtains an iterative procedure, where, at each
iteration, $\Gamma$ is the minimum of a slightly different problem as that of
Eq. (II.3) with $||A\Gamma-X||\rightarrow 0$ at convergence (see
Refs.Verstichel (2012); Mazziotti (2011) for details).

The boundary point method has been successfully used in the variational 2RDM
calculations for atomic systems, as the calculated ground-state energies agree
well with those obtained from those implemented from full configuration
interaction approaches Mazziotti (2011); Verstichel  _et al._ (2012);
Kamarchik and Mazziotti (2007).

The number of iterations to reach convergence when using the augmented
Lagrangian method is about one to a few thousands for atomic and molecular
systems Verstichel (2012); Mazziotti (2011). The situation is similar when
applying this method to nuclei according to our numerical tests. However, the
number of iterations can easily reach 10,000 due to the rather large strength
of the nucleon-nucleon interaction. Added to that, it can happen that the
augmented Lagrangian method does not converge when $N$-$representability$
constraints are not sufficiently strong (see Sec.(III)).

Because the variational 2RDM conditions are independent of a reference wave
function, they can treat both moderate and strong coupling systems, such as
atomic and nuclear systems. Furthermore, closed-shell and open-shell systems,
as well as nuclei with an even or odd number of protons and neutrons, are
treated equally in the variational 2RDM calculations, so that it is as general
as configuration interaction methods from that point of view. This is in
contrast with coupled-cluster theory, for example, where a reference state is
used and only closed-shell systems plus or minus one to two nucleons can be
used Hagen  _et al._ (2016).

##  III Results

###  III.1 Variational 2RDM calculations with cluster-orbital shell-model
Hamiltonian

The first example will be considered with a phenomenological Hamiltonian for
simplicity, where one uses the core + valence particle picture. We will then
apply the variational 2RDM method using the cluster-orbital shell-model (COSM)
Hamiltonian, which is the natural framework when using core and valence
nucleons with an effective Hamiltonian Suzuki and Ikeda (1988). The COSM
coordinates of the valence nucleons are defined from the center of the core
nucleus, and the kinetic energy of the center of mass motion is subtracted
from the total $A$-body Hamiltonian. As a consequence, no center of mass
motion can arise in COSM. In fact, the COSM coordinates form a set of Jacobi
coordinates. The Hamiltonian in the COSM frame reads:

$\displaystyle\hat{H}=\sum_{i=1}^{N_{\text{val}}}\left(\frac{\hat{\mathbf{p}}_{i}^{2}}{2\mu_{i}}+\hat{U}^{(c)}_{i}\right)+\sum_{i<j}^{N_{\text{val}}}\left(\hat{V}_{ij}+\frac{\hat{\mathbf{p}}_{i}\cdot\hat{\mathbf{p}}_{j}}{M_{c}}\right),$
(41)

where $N_{\text{val}}$ is the number of valence nucleons, $\mu_{i}$ is the
reduced mass of the $i$-th nucleon with respect to the core, and $M_{c}$ is
the mass of the core. The single-particle potential $\hat{U}^{(c)}_{i}$ mimics
the field generated by the core on each nucleon. It is represented by a Woods-
Saxon (WS) potential with a spin-orbit term. $\hat{V}_{ij}$ is the two-body
residual interaction, which is composed of the nucleon-nucleon interaction and
the two-body Coulomb interaction. The additional two-body kinetic term in Eq.
(41), depending on $M_{c}$, is the recoil term of the COSM Hamiltonian and
arises from the translationally invariant character of the COSM framework.

As we aim at comparing the numerical results arising from the variational 2RDM
method in nuclear systems to exact eigenenergies, it is sufficient to consider
a simple nucleon-nucleon interaction for $V_{ij}$. Thus, the surface Gaussian
interaction (SGI) Jaganathen  _et al._ (2014) is used :

$\displaystyle V_{\text{SGI}}(\mathbf{r}_{i}-\mathbf{r}_{j})$
$\displaystyle=V_{0}~{}\text{exp}\left[-\left(\frac{\mathbf{r}_{i}-\mathbf{r}_{j}}{\mu_{I}}\right)^{2}\right]\delta(r_{i}-r_{j}-2R_{0}),$
(42)

where $\mu_{I}$ is the interaction range; $V_{0}$ is the strength of the
interaction and $R_{0}$ is the nuclear radius. In the present work, we fix the
interaction range and nucleus radius to $\mu_{I}$ = 1 fm and $R_{0}=2$ fm,
respectively. $V_{0}$ is adjusted to reproduce the selected experimental data.

As a first application of the variational 2RDM method, we will compare the
variational eigenenergies of 6,8He ground states obtained in the variational
2RDM method with shell-model results, which can be deemed as exact for that
matter. 4He is selected as the inner core, and the parameters of the one-body
WS potential mimicking the core are taken from Ref.Jaganathen  _et al._
(2017). Calculations are performed in a model space spanned by the $sp$
partial waves. For each partial wave, six harmonic oscillator (HO) shells are
used. The model space is consists of 17 orbitals, excluding the core space.
The $V_{0}$ coupling constant of the SDI interaction is adjusted to reproduce
the binding energy of 6He ens . Its strength is equal to $-435$ MeV.

Table 1:  Ground-state energies of 6,8He (in MeV), calculated with SM and variational 2RDM methods using the COSM Hamiltonian. Variational 2RDM calculations are constrained with $PQG$ conditions. Experimental data ens  are also provided. Nuclei | $PQG$ | SM | Exp
---|---|---|---
6He | $-1.090$ | $-0.974$ | $-0.975$
8He | $-5.845$ | $-5.795$ | $-3.289$
Figure 1: Energy and normalized occupations of 8He ground state using SM and
variational 2RDM calculations. COSM calculations (brown diamond) are performed
using different truncation schemes from 0p0h to 4p4h (exact). Variational 2RDM
calculations (red diamond) are constrained with $PQG$ conditions.

Results obtained with variational 2RDM and shell-model are shown in Tab. 1,
along with experimental data ens . Due to the small model space used, the
obtained 8He ground-state energy is overbound in our calculations when
compared with experimental data ens . However, this is not important for our
purpose as we only aim at benchmarking the variational 2RDM results with exact
energies. We can see that the variational 2RDM calculations of 6,8He are close
to the exact COSM results, where the largest difference is about 100 keV.

Let us now compare partial wave occupation in variational 2RDM and COSM
frameworks for 8He. For this, one performs several COSM calculations in
truncated spaces, where one allows 0, 1, $\dots$, 4 neutrons in orbitals
outside the $0p$-shell. These truncated spaces are labeled by 0p0h, $\dots$,
4p4h, respectively. Normalized occupation of the $sp$ partial waves are shown
in Fig. 1(b). Normalized occupations are obtained by being divided by the
number of valence particles. As can be expected, the occupation of $s_{1/2}$
and $p_{1/2}$ partial waves increases with the size of the model space, as
they are generated by the configuration mixing induced by the nucleon-nucleon
interaction. Occupations of $sp$ partial waves in the variational 2RDM and
COSM methods become close in 3p3h and 4p4h. At the same time, the ground-state
energy of 8He is almost the same in both methods. This means that the 2RDM
calculation of the ground state obtained with $PQG$ conditions is close to the
exact 2RDM. Hence, in this case, the variational 2RDM allows to calculate
almost exactly the ground-state energy of 8He, even when the simplest
variational 2RDM framework is used, where only the N-representability $PQG$
conditions are constrained.

Figure 2: Similar to Fig. 1(a), but for 8Be. Figure 3: Similar to Fig. 1(b),
but for proton and neutron occupations in 8Be.

In the calculations of 6,8He ground states, only neutrons are considered in
the valence space. Consequently, it is interesting to consider a nucleus in
which both valence protons and neutrons are present. Thus, we performed
variational 2RDM calculations for the 8Be ground state, where two valence
protons and two valence neutrons are present outside the 4He core. The same
COSM Hamiltonian as for 6,8He is used for 8Be, except that the Coulomb force
is added in proton space. The model space used for protons is also the same as
that of neutrons. The ground state of 8Be obtained in the variational 2RDM
with $PQG$ conditions are presented in Fig. 2, along with COSM calculations
with different truncations. We can see that the ground-state energy of 8Be
calculated using variational 2RDM with $PQG$ condition is close to the exact
results. However, the calculated normalized occupations (see Fig. 3) show that
large differences occur between COSM and variational 2RDM calculations. A
large percentage of $s_{1/2}$ partial wave occupation occurs using variational
2RDM, of about 12%, while it is only about 2% in the exact wave function
provided by diagonalization within COSM. Consequently, the occupation of
$p_{3/2}$ partial waves is reduced due to the large $s_{1/2}$ component in the
2RDM calculation, while that of the $p_{1/2}$ partial is almost the same in
both variational RDM and COSM calculations. This discrepancy translates into
an energy overbinding of the 8Be ground state using the variational 2RDM
method (see Fig. 2). Indeed, the 8Be ground state obtained with the
variational 2RDM method is too bound by 1.7 MeV compared with the COSM
calculation. Thus, this energy overbinding is generated by spurious couplings
mainly involving the $s_{1/2}$ and $p_{3/2}$ partial waves.

Hence, in order to obtain precise partial wave occupations of the 8Be ground
state using the variational 2RDM, it is necessary to include higher-order
conditions Erdahl (1978); Zhao  _et al._ (2004) or Hamiltonian-based
constraints, i.e., inequalities involving the expectation value of two-body
operators deduced from the variational principle Shenvi and Izmaylov (2010).
However, the computational cost of the variational 2RDM induced by the $T_{1}$
and $T_{2}^{\prime}$ conditions is expensive when more than a few basis
orbitals are present (see Sec.(II.2)), so that we could not include these
constraints in the variational 2RDM calculations using the COSM Hamiltonian.

The variational 2RDM calculations with $PQG$ condition could provide good
descriptions of 6,8He ground-state energies. The partial wave occupations
obtained with the variational 2RDM method are close to those arising from COSM
calculations using many-body wave functions. For 8Be, the obtained energy with
the variational 2RDM method is close to the exact value. However, the partial
wave occupations obtained therein show discrepancies when compared to exact
values.

###  III.2 Variational 2RDM calculations with USDB shell-model interactions

The variational 2RDM calculations with $PQG$ conditions can provide the
ground-state energies of 6,8He and 8Be using a COSM Hamiltonian within a model
space built from $sp$ partial waves. However, the partial wave occupations
obtained for 8Be are different from the exact values by about 10 %. This
discrepancy might arise from the fact that one has both valence protons and
neutrons or because we only constrained the variational 2RDM calculation using
$PQG$ conditions. Hence, one would like to be able to apply $T_{1}$ and
$T_{2}^{\prime}$ conditions as well. For this, one can only have a few basis
orbitals, as otherwise calculations become too expensive numerically.
Consequently, we will study the variational 2RDM method using the USDB
interaction Brown and Richter (2006) as a testing ground. The USDB interaction
is a standard shell-model interaction for $sd$-shell nuclei, which can provide
good descriptions of the nuclei, such as their the binding energies and
spectra Brown and Richter (2006). The model space for valence protons and
neutrons consists of the $0d_{5/2}$, $1s_{1/2}$ and $0d_{3/2}$ shells in the
USDB interaction, so that the small number of basis orbitals makes this
framework convenient to study the effect of $T_{1}$ and $T_{2}^{\prime}$
conditions.

For the USDB interaction, the Hamiltonian is written as

$\displaystyle\hat{H}=\hat{H}_{0}+\hat{H}_{1},$ (43)

where $\hat{H}_{0}$ represents the single-particle energies of orbitals in the
model space, $\hat{H}_{1}$ is the residual two-body interaction. Isospin-
breaking effects are not considered in the USDB interaction Brown and Richter
(2006). In the present work, we will first consider the ground-state energies
of even-even oxygen isotopes, in which the nuclei only have valence neutrons;
after that, we will focus on the $N=Z$ even-even nuclei, thus bearing both
valence protons and neutrons, which consist of 20Ne, 24Mg, 28Si, and 32S. Due
to the isospin conserving character of the USDB interaction, the
configurations of valence neutrons and protons are the same in $N=Z$ nuclei.
In the latter results, we then just show the normalized occupations of valence
neutrons.

Table 2:  Ground-state energies of oxygen isotopes, calculated with the variational 2RDM method using the USDB interaction. Variational 2RDM calculations are constrained with $PQG$, $PQGT_{1}$, $PQGT_{2}^{\prime}$, or $PQGT_{1}T_{2}^{\prime}$ conditions. The star indicates that the variational 2RDM does not converge using the considered constraints. Results are compared with SM calculations. Nuclei |  20O |  22O |  24O |  26O
---|---|---|---|---
$PQG$ | $-30.344$ | $-58.055$ | $*$ | *
$PQGT_{1}$ | $-29.701$ | $-36.404$ | $-42.317$ | *
$PQGT_{2}^{\prime}$ | $-23.712$ | $-34.554$ | $-41.257$ | $-40.869$
$PQGT_{1}T_{2}^{\prime}$ | $-23.712$ | $-34.554$ | $-41.256$ | $-40.866$
SM | $-23.632$ | $-34.498$ | $-41.225$ | $-40.868$
Figure 4: Variational 2RDM calculations of oxygen isotopes with different
N-representability conditions, i.e., $PQG$, $PQGT_{1}$, $PQGT_{2}^{\prime}$,
and $PQGT_{1}T_{2}^{\prime}$ constraints. Results are compared with SM
calculations. Figure 5: Calculated normalized occupations of the ground-state
energies of oxygen isotopes using variational 2RDM, with different
N-representability conditions, and SM.

The calculated ground-state energies of oxygen isotopes using variational 2RDM
with different N-representability conditions are presented in Tab. 2 and Fig.
4. The results are also compared with the exact values computed from a SM
calculation. We focus on the ${}^{20\text{-}26}$O even-even isotopes in the
oxygen chain. For 24,26O, the variational 2RDM calculations using $PQG$ and/or
$PQGT_{1}$ conditions are not shown because one could not obtain convergence
of the variational 2RDM method. From Fig. 4 and Tab. 2, we can see that the
results using variational 2RDM calculations with $PQG$ conditions are
overbound when compared with the exact results. For example, the ground-state
energy of 22O is overbound by about 25 MeV, which is excessively large. The
situation can be slightly improved by including the $T_{1}$ condition in
$PQG$, but the results of variational 2RDM calculations with $PQGT_{1}$ are
still overbound by a few MeV. The energy overbinding is fully overcome in the
variational 2RDM by including the $T_{2}^{\prime}$ condition. The results of
the variational 2RDM with $PQGT_{2}^{\prime}$ and $PQGT_{1}T_{2}^{\prime}$ are
indeed very close to the exact results, with differences within 100 keV. One
can also note that the results obtained using the $PQGT_{1}T_{2}^{\prime}$ and
$PQGT_{2}^{\prime}$ conditions are nearly the same. This is consistent with
the conclusions obtained in the study of molecular systems Mazziotti (2005,
2006a); Verstichel (2012).

Let us now consider partial-wave occupations (see Fig. 5). The variational
2RDM calculations obtained with $PQGT_{2}^{\prime}$ or
$PQGT_{1}T_{2}^{\prime}$ are close to the exact results. However, the
calculations with $PQG$ or $PQGT_{1}$ do not give accurate calculations. In
Tab. 2, the detailed values of the obtained ground-state energies are
presented. We can see that the largest difference is about 80 keV (which
happens in 20O) when the variational 2RDM calculations using
$PQGT_{1}T_{2}^{\prime}$ condition are compared with the exact shell-model
results. For 26O, the variational 2RDM calculation with
$PQGT_{1}T_{2}^{\prime}$ conditions is slightly higher than the exact results,
which may be caused by the computational error in the variational
calculations.

In order to ponder out the influence of the isospin $T=0$ part of the nucleon-
nucleon interaction, we now consider nuclei in which protons and neutrons are
both active. As testing grounds for our benchmarking purpose, we calculate the
ground states of the $N=Z$ even-even $sd$-nuclei, 20Ne, 24Mg, 28Si, and 32S.
The results are shown in Tab. 3 and Fig. 6. Similar to oxygen isotopes, we see
that the variational 2RDM calculations using $PQG$ conditions are too
overbound. The situation is even worse when one has both valence protons and
neutrons, because overbinding can reach 50-100 MeV. This is clearly due to the
presence of $T=0$ matrix elements in the nucleon-nucleon interaction, which
are about twice larger than $T=1$ matrix elements. However, the situation is
not improved when including the $T_{1}$ condition along with $PQG$ conditions.
In fact, only the consideration of the $T_{2}^{\prime}$ constraint could
significantly ameliorate the quality of calculations. Note that the
variational 2RDM calculations with $PQGT_{1}T_{2}^{\prime}$ constraints are
similar to those obtained with $PQGT_{2}^{\prime}$ conditions. However, the
discrepancies between variational 2RDM using $PQGT_{1}T_{2}^{\prime}$ and
exact results using SM are still large. The largest energy difference occurs
in 28Si, where it is about 6.5 MeV. Conversely, the energy discrepancy is
small in 20Ne, as it is only about 150 keV. Normalized partial-wave
occupations are presented in Fig. 7. The variational 2RDM calculations with
$PQG$ and $PQGT_{1}$ conditions give similar results, as both of them show
large discrepancies when compared with exact results. Partial-wave occupations
are largely improved by including the $T_{2}^{\prime}$ condition. The
variational 2RDM calculations using $PQGT_{2}^{\prime}$ and
$PQGT_{1}T_{2}^{\prime}$ provide nearly the same occupations. The results are
close to SM calculations, but the variational 2RDM calculations, including the
$T_{2}^{\prime}$ condition, still show small discrepancies. This explains
ground-state energy overbinding (See Fig. 6).

Table 3:  Similar to Tab. 2, but for the $N=Z$ even-even nuclei 20Ne, 24Mg, 28Si, and 32S. Nuclei |  20Ne |  24Mg |  28Si |  32S
---|---|---|---|---
$PQG$ | $-44.160$ | $-129.809$ | $-243.363$ | $-219.061$
$PQGT_{1}$ | $-44.167$ | $-129.808$ | $*$ | $-219.058$
$PQGT_{2}^{\prime}$ | $-41.536$ | $-91.656$ | $-142.332$ | $-186.186$
$PQGT_{1}T_{2}^{\prime}$ | $-41.537$ | $-91.672$ | $-142.314$ | $-186.314$
SM | $-41.397$ | $-87.104$ | $-135.861$ | $-182.452$
Figure 6: Similar to Fig. 4, but for the $N=Z$ even-even nuclei, 20Ne, 24Mg,
28Si, and 32S. Figure 7: Similar to Fig. 5, but for the $N=Z$ even-even nuclei
20Ne, 24Mg, 28Si, and 32S.

The variational 2RDM calculations using USDB Hamiltonian have shown that the
$PQGT_{1}T_{2}^{\prime}$ conditions tremendously improve the precision of
calculations. While $PQG$ conditions can lead to an energy overbinding which
may even reach 100 MeV, it is decreased by more than an order of magnitude
with $PQGT_{1}T_{2}^{\prime}$ conditions. However, a precision of 5-10 MeV is
still not sufficient for nuclear binding energy. Stronger or additional
physical conditions Mazziotti (2012) should be included to improve the
calculations of those nuclei.

Nevertheless, $PQG$ conditions can provide good descriptions for atomic or
molecular systems Hammond and Mazziotti (2005), where numerical precision is
comparable to that of coupled-cluster calculations with singles and doubles
excitations Mazziotti (2006a); Hammond and Mazziotti (2005). The ground states
of 6,8He can also be well reproduced using only $PQG$ constraints with a COSM
Hamiltonian within an $sp$ model space.

The crux of the problem seems to be how soft the inter-particle interaction
is. In atomic or molecular systems, one only has an interaction of the Coulomb
type, whose induced inter-particle correlations are weak. Indeed, the Hartree-
Fock approximation provides with more than 99% of the exact wave function,
while electron-electron correlation energy is typically 0.1% of the total
binding energy Mazziotti (2016). This is in sharp contrast with the nucleon-
nucleon interaction, which provides correlations so strong that they can
generate nuclear deformation and clustering Freer  _et al._ (2018). In fact,
mean-field and residual interaction have comparable effects on nuclear many-
body wave functions, as SM calculations have shown Caurier  _et al._ (2005).

Table 4:  Similar to Tab. 3, but using the USDB(1/2) interaction in the Hamiltonian. Nuclei |  20Ne |  24Mg |  28Si |  32S
---|---|---|---|---
$PQG$ | $-28.098$ | $-77.242$ | $-132.453$ | $-130.003$
$PQGT_{1}$ | $-28.089$ | $-72.090$ | $-96.748$ | $-124.090$
$PQGT_{2}^{\prime}$ | $-27.150$ | $-58.852$ | $-90.886$ | $-119.759$
$PQGT_{1}T_{2}^{\prime}$ | $-27.011$ | $-58.848$ | $-90.885$ | $-119.756$
SM | $-26.627$ | $-57.007$ | $-89.294$ | $-119.298$
Figure 8: Similar to Fig. 6, but where the USDB(1/2) interaction is utilized
as Hamiltonian. Figure 9: Similar to Fig. 7, but where the USDB(1/2)
interaction is utilized as Hamiltonian.

Consequently, we also investigate the dependence on interaction strength by
considering a softer nucleon-nucleon interaction. For this, we multiply the
matrix elements of the USDB effective interaction by 1/2, while keeping the
same single-particle energies. We denote this interaction by USDB(1/2). As the
contribution of the USDB(1/2) residual two-body interaction is much smaller
than that of the initial USDB interaction, we are able to quantify the error
induced by the approximate N-representability character of the 2RDM issued
from energy minimization. Clearly, the USDB(1/2) interaction can only be used
for benchmarking purposes and is not suitable to describe $sd$-shell nuclei.
We then reconsider the same nuclei with the variational 2RDM calculations with
the USDB(1/2) interaction as those with the initial USDB interaction, i.e.,
the $N=Z$ $sd$-shell nuclei 20Ne, 24Mg, 28Si, and 32S.

The calculated energies of these nuclei using the variational 2RDM and shell-
model method are shown in Fig. 8 and Tab. 4. Even when using a softer nucleon-
nucleon interaction, the variational 2RDM calculations with $PQG$ condition
still cannot provide an adequate description of those nuclei. Indeed, the
obtained energies are overbound when compared with the exact values by up to
about 45 MeV. Similar to the previous studies, the $T_{1}$ condition gives
only but a small contribution to binding energies with the USDB(1/2)
interaction, while the inclusion of the $T_{2}^{\prime}$ condition largely
improves the calculations. Moreover, the results obtained with
$PQGT_{2}^{\prime}$ and $PQGT_{1}T_{2}^{\prime}$ are nearly the same based on
the USDB(1/2) interaction. The largest energy discrepancy found in the binding
energies of the considered nuclei using the USDB(1/2) 2RDM with
$PQGT_{1}T_{2}^{\prime}$ conditions is about 1.8 MeV. This value is smaller
than that obtained with the USDB interaction in the same conditions, where it
is about 6.5 MeV. The calculated occupations in the variational 2RDM using
USDB(1/2) interaction follow a similar trend as that with the USDB
interaction. (see Fig. 9). Due to the soft character of the USDB(1/2)
interaction, the occupations obtained in the variational 2RDM calculations
with $PQGT_{1}T_{2}^{\prime}$ or $PQGT_{2}^{\prime}$ conditions are nearly the
same as those provided by SM calculations.

##  IV Towards a practical variational 2RDM method for atomic nuclei

The error induced by the approximate N-representability of the variational
2RDM strongly depends on the correlations induced by the two-body interaction
of the considered Hamiltonian. As mentioned in Ref.Verstichel (2012), the 1RDM
provided by the variational RDM calculation is always exactly N-representable,
because this property is equivalent to have the 1RDM eigenvalues to be between
0 and 1. Consequently, the ground-state energy arising from a one-body
Hamiltonian has to be exactly reproduced in a variational RDM calculation.
Conversely, overbinding arises when including two-body forces, which increases
along with the strength of the interaction. While this discrepancy is mild in
atomic and molecular calculations in the presence of $PQG$ conditions, errors
can still reach several MeVs in nuclear binding energies by using
$PQGT_{1}T_{2}^{\prime}$ constraints, even after making the nuclear
interaction softer by a factor of two. Indeed, the maximal error is reduced by
about a factor of 3.5 after dividing the matrix elements of the USDB
interaction by 2, i.e., they are of the same order of magnitude.

Therefore, a direct implementation of the variational 2RDM method cannot be
used in practice for atomic nuclei, contrary to atomic and molecular systems,
where it can directly provide almost exact binding energies. While the
situation has been much improved since the original work of Mihailović and
Rosina of Ref.Mihailović and Rosina (1975), it is far from being satisfactory,
as the typical error made on nuclear binding energies in systematic
calculations such as density functional theory (DFT) Bender  _et al._ (2003);
Kortelainen  _et al._ (2010); Bulgac  _et al._ (2018) and microscopic-
macroscopic models Bhagwat  _et al._ (2012); Möller  _et al._ (2012) is
typically 0.5 to 1 MeV. Clearly, one has to devise alternative strategies in
order to develop a variational 2RDM method which can be of practical interest
in nuclear physics.

As mentioned in Sec.(III), one of the most interesting paths for that matter
would be to use the variational principle in small model spaces Shenvi and
Izmaylov (2010). For example, one can diagonalize exactly the Hamiltonian
matrix if one takes into account only a few basis orbitals. The obtained
energy provides with a useful constraint for the exact ground state. This
method has been implemented, and one has seen that it accelerates convergence
significantly. However, it has been noticed during the first tests that the
constraint works well only if the correlation energy provided by the rest of
the model space is of a few hundreds of keV at most. As consequence,
constraints based on the variational principle are interesting mainly when
eigenstates are built with a very large part from a few configurations of low
energy. This is especially the case in the Gamow shell model Id Betan  _et
al._ (2002); Michel  _et al._ (2002, 2009); Li  _et al._ (2020), where
continuum coupling is small, so that the $0p$ or $1s0d$ orbitals are mainly
occupied Michel  _et al._ (2009). Moreover, according to the results obtained
in Sec.(III.1) for 6,8He and 8Be, it might also be possible that $PQG$
constraints are sufficient to impose N-representability in the presence of
constraints using the variational principle in small model spaces.

The situation is reversed when one enters the zones of the nuclear chart
associated to medium and heavy nuclei. From a numerical point of view, the
variational 2RDM method is well suited for ground states therein, because one
typically has few basis orbitals but many valence nucleons in model space.
Consequently, it is possible to apply in practice the $PQGT_{1}T_{2}^{\prime}$
constraints. Moreover, the variational 2RDM method can be applied easily to
nuclei having an even or odd number of protons or neutrons, contrary to the
Hartree-Fock-Bogoliubov method, for example, which becomes cumbersome to apply
with odd nuclei Bender  _et al._ (2003). Contrary to light nuclei, however,
the main problem here is that the nuclear interaction acts similarly on all
basis orbitals, and shell-model wave functions spread in a large configuration
space Caurier  _et al._ (2005). Hence, the variational 2RDM method may lead to
overbound ground states. In fact, the calculations depicted in Sec.(III.2)
have shown that this is already the case in the $sd$ shell. Consequently, it
might be of interest here to minimize the ground-state energies of Hamiltonian
whose two-body part is reduced. The overbinding caused by the approximate
N-representability of the 2RDM issued from energy minimization would then be
compensated by the use of smaller two-body matrix elements. The 2RDM obtained
by energy minimization could then become closer to that of the physical ground
state, even though it would arise from the use of a different Hamiltonian.
This procedure can be related to the Hohenberg-Kohn and Kohn-Sham theories
Hohenberg and Kohn (1964); Kohn and Sham (1965), at the origin of DFT Bender
_et al._ (2003), where the nuclear density, as a function of the correlated
ground state, is determined from a mean-field functional. The validity of this
method could be checked in shell-model spaces where diagonalization of
Hamiltonian can be effected, or by using solvable Hamiltonians for that
matter, of SU(3) type for example Elliott (1958a, b).

##  V Summary

We have developed the variational two-particle reduced density matrix method
for atomic nuclei, and have applied it to the ground states of nuclear
systems. The $N\text{-}representability$ conditions, represented by the two-
particle $P$, $Q$, and $G$ metric matrices, as well as the partial three-
particle $T_{1}$ and $T_{2}^{\prime}$ metric matrices, have been derived and
applied for the first time within the $J$-scheme formalism. In addition, the
neutron and proton numbers $N$ and $Z$, total spin $J$ and the isospin $T$ are
considered as constraints in our calculations.

We have firstly applied the variational two-particle reduced density matrix
method within the cluster-orbital shell-model framework. This has allowed to
use a phenomenological interaction, which is convenient for benchmarking
purposes. The 6,8He and 8Be ground-state energies have been calculated using
$PQG$ constraints, and compared with shell-model calculations, which can be
considered as the exact result (since shell-model dimensions are tractable).
The results show that the two-particle reduced density matrix calculations
could provide good descriptions of ground-state energies and partial-wave
occupations in 6,8He. However, in the case of 8Be, the partial-wave
occupations exhibits a clear difference from the shell-model calculation and a
small but noticeable overbinding of ground-state energy. Due to the large
number of basis orbitals in the cluster-orbital shell-model model space, only
the $PQG$ conditions could be imposed, because computational cost becomes too
large if one includes the partial 3-positive $T_{1}$ and $T_{2}^{\prime}$
conditions. Consequently, one had to consider a different Hamiltonian and
shell-model space to study the impact of the $T_{1}$ and $T_{2}^{\prime}$
constraints.

We have then employed the variational two-particle reduced density matrix
method using the standard USDB interaction in the $sd$ shell. Indeed, the
number of basis orbitals is small, on the one hand, and Hamiltonian matrices
can always be diagonalized exactly using the shell-model framework. Hence, we
could perform variational two-particle reduced density matrix calculations
with all possible sets of conditions, i.e., $PQG$, $PQGT_{1}$,
$PQGT_{2}^{\prime}$, and $PQGT_{1}T_{2}^{\prime}$. We firstly studied the
oxygen even-mass isotopes. The results showed that the calculations with $PQG$
conditions lead to overbound ground-state energies. The situation is largely
improved when the $T_{2}^{\prime}$ condition is imposed. The variational two-
particle reduced density matrix calculations with $PQGT_{2}^{\prime}$ and
$PQGT_{1}T_{2}^{\prime}$ can then provide an accurate description of oxygen
isotopes. The largest difference of ground-state energies compared to exact
values is about 100 keV, so that calculations using the variational two-
particle reduced density matrix method can be considered as optimal.
Furthermore, the obtained occupations with $PQGT_{2}^{\prime}$ or
$PQGT_{1}T_{2}^{\prime}$ conditions are nearly the same as those arising from
shell-model calculations.

However, the calculations of $N=Z$ even-even mass isotopes in that same
framework exhibit a large discrepancy with exact results, even in the presence
of $PQGT_{2}^{\prime}$ or $PQGT_{1}T_{2}^{\prime}$ constraints. The largest
error is indeed as large as 6.5 MeV, which occurs in 28Si. In order to
investigate the effect of Hamiltonian correlations in ground-state energy
errors, we made the USDB Hamiltonian softer by halving its two-body matrix
elements. The error in the ground-state energy of 28Si is then decreased to
1.8 MeV with $PQGT_{2}^{\prime}$ or $PQGT_{1}T_{2}^{\prime}$ conditions, which
is smaller than with the initial USDB interaction, but remains too large
nevertheless compared to those obtained with other methods to solve the
nuclear many-body problem.

These results then suggest that the present variational two-particle reduced
density matrix method cannot be applied satisfactorily to nuclear systems as
in atomic or molecular physics. Indeed, it is impossible to reach the same
precision for binding energy due to the strength of inter-nucleon
correlations, even if one uses $PQGT_{1}T_{2}^{\prime}$ constraints.
Theoretical frameworks have then been proposed, which will be the object of
future works. The use of the variational principle in reduced model spaces
would then provide with an additional stringent condition in light nuclei,
because the most important configurations of their ground states are built
with the orbitals of lowest energy. Therefore, it would be more interesting to
optimize the 2RDM of a less correlated Hamiltonian, so that the energy
overbinding arising from the approximate N-representability of the 2RDM would
be mitigated by the use of a softer interaction.

Acknowledgments – This work has been supported by the National Key R&D Program
of China under Grant No. 2018YFA0404401; the National Natural Science
Foundation of China under Grants No. 11835001, 11921006, 12035001, and
11975282; the State Key Laboratory of Nuclear Physics and Technology, Peking
University under Grant No. NPT2020KFY13; the Strategic Priority Research
Program of Chinese Academy of Sciences under Grant No. XDB34000000; and the
CUSTIPEN (China-U.S. Theory Institute for Physics with Exotic Nuclei) funded
by the U.S. Department of Energy, Office of Science under Grant No. de-
sc0009971. We acknowledge the High-Performance Computing Platform of Peking
University for providing computational resources.

## References

  * Caurier  _et al._ (2005) E. Caurier, G. Martínez-Pinedo, F. Nowacki, A. Poves, and A. P. Zuker, Rev. Mod. Phys. 77, 427 (2005).
  * Otsuka  _et al._ (2020) T. Otsuka, A. Gade, O. Sorlin, T. Suzuki, and Y. Utsuno, Rev. Mod. Phys. 92, 015002 (2020).
  * Carlson  _et al._ (2015) J. Carlson, S. Gandolfi, F. Pederiva, S. C. Pieper, R. Schiavilla, K. E. Schmidt, and R. B. Wiringa, Rev. Mod. Phys. 87, 1067 (2015).
  * Hergert  _et al._ (2016) H. Hergert, S. Bogner, T. Morris, A. Schwenk, and K. Tsukiyama, Phys. Rep 621, 165 (2016), memorial Volume in Honor of Gerald E. Brown.
  * Stroberg  _et al._ (2019) S. R. Stroberg, H. Hergert, S. K. Bogner, and J. D. Holt, Annu. Rev. Nucl. Part. Sci 69, 307 (2019).
  * Hagen  _et al._ (2016) G. Hagen, M. Hjorth-Jensen, G. Jansen, and T. Papenbrock, Phys. Scr 91, 063006 (2016).
  * Mayer (1955) J. E. Mayer, Phys. Rev. 100, 1579 (1955).
  * Nakatsuji and Yasuda (1996) H. Nakatsuji and K. Yasuda, Phys. Rev. Lett. 76, 1039 (1996).
  * Coleman (1963) A. J. Coleman, Rev. Mod. Phys. 35, 668 (1963).
  * Mazziotti (2004) D. A. Mazziotti, Phys. Rev. Lett. 93, 213001 (2004).
  * Hammond and Mazziotti (2005) J. R. Hammond and D. A. Mazziotti, Phys. Rev. A 71, 062503 (2005).
  * Mazziotti (2005) D. A. Mazziotti, Phys. Rev. A 72, 032510 (2005).
  * Mazziotti (2006a) D. A. Mazziotti, Phys. Rev. A 74, 032501 (2006a).
  * Verstichel  _et al._ (2009) B. Verstichel, H. van Aggelen, D. Van Neck, P. W. Ayers, and P. Bultinck, Phys. Rev. A 80, 032508 (2009).
  * Nakata and Yasuda (2009) M. Nakata and K. Yasuda, Phys. Rev. A 80, 042109 (2009).
  * Mazziotti (2012) D. A. Mazziotti, Phys. Rev. Lett. 108, 263002 (2012).
  * Verstichel  _et al._ (2012) B. Verstichel, H. van Aggelen, W. Poelmans, and D. Van Neck, Phys. Rev. Lett. 108, 213001 (2012).
  * Mazziotti (2016) D. A. Mazziotti, Phys. Rev. Lett. 117, 153001 (2016).
  * Löwdin (1955) P.-O. Löwdin, Phys. Rev. 97, 1474 (1955).
  * Coleman (2000) A. J. Coleman, in  _Many-electron densities and reduced density matrices_ (Kluwer Academic/Plenum Publishers, Springer, Boston, MA, 2000) pp. 1–17.
  * Garrod and Percus (1964) C. Garrod and J. K. Percus, J. Math. Phys. 5, 1756 (1964).
  * Zhao  _et al._ (2004) Z. Zhao, B. J. Braams, M. Fukuda, M. L. Overton, and J. K. Percus, Jour. Chem. Phys. 120, 2095 (2004).
  * Mazziotti (2011) D. A. Mazziotti, Phys. Rev. Lett. 106, 083001 (2011).
  * Anderson  _et al._ (2013) J. S. Anderson, M. Nakata, R. Igarashi, K. Fujisawa, and M. Yamashita, Comp. Theo. Chem 1003, 22 (2013).
  * Nakata  _et al._ (2001) M. Nakata, H. Nakatsuji, M. Ehara, M. Fukuda, K. Nakata, and K. Fujisawa, Jour. Chem. Phys. 114, 8282 (2001).
  * Mazziotti (2006b) D. A. Mazziotti, Phys. Rev. Lett. 97, 143002 (2006b).
  * Mazziotti (2020) D. A. Mazziotti, Phys. Rev. A 102, 030802(R) (2020).
  * Li and DePrince (2019) R. R. Li and A. E. DePrince, Phys. Rev. A 100, 032509 (2019).
  * Maradzike  _et al._ (2017) E. Maradzike, G. Gidofalvi, J. M. Turney, H. F. Schaefer, and A. E. DePrince, Jour. Chem. Theo. Comp. 13, 4113 (2017).
  * Kamarchik and Mazziotti (2007) E. Kamarchik and D. A. Mazziotti, Phys. Rev. Lett. 99, 243002 (2007).
  * Antipin and Vasilieva (2015) A. Antipin and O. Vasilieva, in  _Analysis, Modelling, Optimization, and Numerical Techniques_ , edited by G. O. Tost and O. Vasilieva (Springer International Publishing, Cham, 2015) pp. 1–36.
  * Mahdavi and Shiri (2015) A. Mahdavi and M. E. Shiri, Comput. Opt. Appl. 60, 263 (2015).
  * Zhang  _et al._ (2006) K. Zhang, X. Yang, and K. Teo, Automatica 42, 1407 (2006).
  * Mihailović and Rosina (1975) M. Mihailović and M. Rosina, Nucl. Phys. A 237, 221 (1975).
  * Cohen and Kurath (1965) S. Cohen and D. Kurath, Nucl. Phys. A 37, 24 (1965).
  * Halbert  _et al._ (1968) E. C. Halbert, J. B. McGrory, and B. H. Wildenthal, Phys. Rev. Lett. 20, 1112 (1968).
  * Zuker  _et al._ (1968) A. P. Zuker, B. Buck, and J. B. McGrory, Phys. Rev. Lett. 21, 39 (1968).
  * Zuker (1969) A. P. Zuker, Phys. Rev. Lett. 23, 983 (1969).
  * Assié and Lacroix (2009) M. Assié and D. Lacroix, Phys. Rev. Lett. 102, 202501 (2009).
  * Tohyama and Schuck (2014) M. Tohyama and P. Schuck, Eur. Phys. Jour. A 50 (2014).
  * Tohyama and Schuck (2017) M. Tohyama and P. Schuck, Eur. Phys. Jour. A 53 (2017).
  * Tohyama and Schuck (2019) M. Tohyama and P. Schuck, Eur. Phys. Jour. A 55 (2019).
  * Barton  _et al._ (2021) M. C. Barton, P. D. Stevenson, and A. Rios, arXiv:2101.07584 [nucl-th] (2021).
  * Erler  _et al._ (2012) J. Erler, N. Birge, M. Kortelainen, W. Nazarewicz, E. Olsen, A. M. Perhac, and M. Stoitsov, Nature 486, 509 (2012).
  * Kortelainen  _et al._ (2010) M. Kortelainen, T. Lesinski, J. Moré, W. Nazarewicz, J. Sarich, N. Schunck, M. V. Stoitsov, and S. Wild, Phys. Rev. C 82, 024313 (2010).
  * Stoitsov  _et al._ (2003) M. V. Stoitsov, J. Dobaczewski, W. Nazarewicz, S. Pittel, and D. J. Dean, Phys. Rev. C 68, 054312 (2003).
  * Bender  _et al._ (2003) M. Bender, P.-H. Heenen, and P.-G. Reinhard, Rev. Mod. Phys. 75, 121 (2003).
  * Bulgac  _et al._ (2018) A. Bulgac, M. M. Forbes, S. Jin, R. N. Perez, and N. Schunck, Phys. Rev. C 97, 044313 (2018).
  * Alcoba1  _et al._ (2018) D. R. Alcoba1, P. Capuzzi, A. Rubio-García, J. Dukelsky, G. E. Massaccesi, O. B. Oña, A. Torre, and L. Lain, J. Chem. Phys. 149, 194105 (2018).
  * Hagen  _et al._ (2007) G. Hagen, T. Papenbrock, D. J. Dean, A. Schwenk, A. Nogga, M. Włoch, and P. Piecuch, Phys. Rev. C 76, 034302 (2007).
  * Verstichel (2012) B. Verstichel,  _Variational determination of the two-particle density matrix as a quantum many-body technique_ (PhD thesis, Gent University, 2012).
  * Erdahl (1978) R. M. Erdahl, Int. Jour. Quantum Chem 13, 697 (1978).
  * von Neumann (1995) J. von Neumann,  _Mathematical Foundations of Quantum Mechanics_ (Princeton University Press, Priceton, 1995).
  * Coleman and Yukalov (2000) A. J. Coleman and V. I. Yukalov,  _Reduced Density Matrix : Coulson’s Challenges_ (Springer, New York, 2000).
  * Mazziotii (2007) D. A. Mazziotii,  _Two-Electron Reduced-Density-Matrix Mechanics_ (Wiley, New York, 2007).
  * Verstichel  _et al._ (2011) B. Verstichel, H. van Aggelen, D. Van Neck, P. Bultinck, and S. De Baerdemacker, Comput. Phys. Comm. 182, 1235 (2011).
  * Staszczak  _et al._ (2010) A. Staszczak, M. Stoitsov, A. Baran, and W. Nazarewicz, Eur. Phys. J. A , 85 (2010).
  * Burer and Monteiro (2003) S. Burer and R. D. C. Monteiro, Math. Prog. 95, 329 (2003).
  * Suzuki and Ikeda (1988) Y. Suzuki and K. Ikeda, Phys. Rev. C 38, 410 (1988).
  * Jaganathen  _et al._ (2014) Y. Jaganathen, N. Michel, and M. Płoszajczak, Phys. Rev. C 89, 034624 (2014).
  * Jaganathen  _et al._ (2017) Y. Jaganathen, R. M. Id Betan, N. Michel, W. Nazarewicz, and M. Płoszajczak, Phys. Rev. C 96, 054316 (2017).
  * (62) https://www.nndc.bnl.gov/ensdf/.
  * Shenvi and Izmaylov (2010) N. Shenvi and A. F. Izmaylov, Phys. Rev. Lett. 105, 213003 (2010).
  * Brown and Richter (2006) B. A. Brown and W. A. Richter, Phys. Rev. C 74, 034315 (2006).
  * Freer  _et al._ (2018) M. Freer, H. Horiuchi, Y. Kanada-En’yo, D. Lee, and U.-G. Meißner, Rev. Mod. Phys. 90, 035004 (2018).
  * Bhagwat  _et al._ (2012) A. Bhagwat, X. Viñas, M. Centelles, P. Schuck, and R. Wyss, Phys. Rev. C 86, 044316 (2012).
  * Möller  _et al._ (2012) P. Möller, W. D. Myers, H. Sagawa, and S. Yoshida, Phys. Rev. Lett. 108, 052501 (2012).
  * Id Betan  _et al._ (2002) R. Id Betan, R. J. Liotta, N. Sandulescu, and T. Vertse, Phys. Rev. Lett. 89, 042501 (2002).
  * Michel  _et al._ (2002) N. Michel, W. Nazarewicz, M. Płoszajczak, and K. Bennaceur, Phys. Rev. Lett. 89, 042502 (2002).
  * Michel  _et al._ (2009) N. Michel, W. Nazarewicz, M. Płoszajczak, and T. Vertse, Jour. Phys. G: Nucl. Part. Phys. 36, 013101 (2009).
  * Li  _et al._ (2020) J. G. Li, B. S. Hu, Q. Wu, Y. Gao, S. J. Dai, and F. R. Xu, Phys. Rev. C 102, 034302 (2020).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Elliott (1958a) J. P. Elliott, Proc. R. Soc. Lond. 245, 128 (1958a).
  * Elliott (1958b) J. P. Elliott, Proc. R. Soc. Lond. 245, 562 (1958b).