Normal-ordered $k$-body approximation in particle-number-breaking theories

TL;DR

This paper introduces a particle-number-conserving normal-ordered k-body approximation for use in symmetry-breaking many-body nuclear theories, enabling more accurate and consistent ab initio calculations involving particle-number symmetry restoration.

Contribution

It develops a novel PNOkB approximation that maintains particle-number symmetry in symmetry-breaking frameworks, improving the consistency of ab initio nuclear calculations.

Findings

PNOkB approximation effectively conserves particle number.

Numerical tests show improved symmetry handling in calculations.

Enables safe use of symmetry-breaking and restored methods in nuclear theory.

Abstract

The reach of ab initio many-body theories is rapidly extending over the nuclear chart. However, dealing fully with three-nucleon, possibly four-nucleon, interactions makes the solving of the A-body Schr\"odinger equation particularly cumbersome, if not impossible beyond a certain nuclear mass. Consequently, ab initio calculations of mid-mass nuclei are typically performed on the basis of the normal-ordered two-body (NO2B) approximation that captures dominant effects of three-nucleon forces while effectively working with two-nucleon operators. A powerful idea currently employed to extend ab initio calculations to open-shell nuclei consists of expanding the exact solution of the A-body Schr\"odinger equation while authorizing the approximate solution to break symmetries of the Hamiltonian. In this context, operators are normal ordered with respect to a symmetry-breaking reference state such that proceeding to a naive truncation may lead to symmetry-breaking approximate operators. The purpose of the present work is to design a normal-ordering approximation of operators that is consistent with the symmetries of the Hamiltonian while working in the context of symmetry broken (and potentially restored) methods. Focusing on many-body formalisms in which U(1) global-gauge symmetry associated with particle number conservation is broken (and potentially restored), a particle-number-conserving normal-ordered k-body (PNOkB) approximation of an arbitrary N-body operator is designed on the basis of Bogoliubov reference states. A numerical test based on particle-number projected Hartree-Fock-Bogoliubov calculations permits to check the particle-number conserving/violating character of a given approximation to a particle-number conserving operator. Using the presently proposed PNOkB approximation, ab initio calculations based on symmetry-breaking and restored formalisms can be safely performed.

Tables

Table 1. Table 1: Contributions to the three-body operator O 𝑂 O in normal-ordered form with respect to the particle vacuum | 0 ⟩ ket 0 |0\rangle and expressed in { c , c † } 𝑐 superscript 𝑐 † \{c,c^{\dagger}\} . The o i ​ j superscript 𝑜 𝑖 𝑗 o^{ij} contributions are sorted horizontally according to i − j 𝑖 𝑗 i-j and vertically according to i + j 𝑖 𝑗 i+j .

$\{c,c^{†}\},|0⟩$	-6	-4	-2	0	+2	+4	+6
$o^{\left[0\right]}$				$o^{00}$
$o^{\left[2\right]}$				$o^{11}$
$o^{\left[4\right]}$				$o^{22}$
$o^{\left[6\right]}$				$o^{33}$

Table 2. Table 2: Contributions to the one-body operator F 𝐹 F in normal-ordered form with respect to the particle vacuum | 0 ⟩ ket 0 |0\rangle and expressed in single-particle basis { c , c † } 𝑐 superscript 𝑐 † \{c,c^{\dagger}\} . The f i ​ j superscript 𝑓 𝑖 𝑗 f^{ij} contributions are sorted horizontally according to i − j 𝑖 𝑗 i-j and vertically according to i + j 𝑖 𝑗 i+j .

$\{c,c^{†}\},|0⟩$	-6	-4	-2	0	+2	+4	+6
$f^{\left[0\right]}$				$f^{00}$
$f^{\left[2\right]}$			$f^{02}$	$f^{11}$	$f^{20}$
$f^{\left[4\right]}$
$f^{\left[6\right]}$

Table 3. Table 3: Contributions to the three-body operator O 𝑂 O in normal-ordered form with respect to the Bogoliubov vacuum | Φ ⟩ ket Φ |\Phi\rangle and expressed in the single-particle basis { c , c † } 𝑐 superscript 𝑐 † \{c,c^{\dagger}\} . The Λ i ​ j superscript Λ 𝑖 𝑗 \Lambda^{ij} contributions are sorted horizontally according to i − j 𝑖 𝑗 i-j and vertically according to i + j 𝑖 𝑗 i+j .

$\{c,c^{†}\},|\Phi ⟩$	-6	-4	-2	0	+2	+4	+6
${\Lambda }^{\left[0\right]}$				${\Lambda }^{00}$
${\Lambda }^{\left[2\right]}$			${\Lambda }^{02}$	${\Lambda }^{11}$	${\Lambda }^{20}$
${\Lambda }^{\left[4\right]}$			${\Lambda }^{13}$	${\Lambda }^{22}$	${\Lambda }^{31}$
${\Lambda }^{\left[6\right]}$				${\Lambda }^{33}$

Table 4. Table 4: Contributions to the three-body operator O 𝑂 O in normal-ordered form with respect to the slater determinant | SD ⟩ ket SD |\text{SD}\rangle and expressed in the single-particle basis { c , c † } 𝑐 superscript 𝑐 † \{c,c^{\dagger}\} . Λ i ​ j superscript Λ 𝑖 𝑗 \Lambda^{ij} contributions are sorted horizontally according to i − j 𝑖 𝑗 i-j and vertically according to i + j 𝑖 𝑗 i+j .

$\{c,c^{†}\},|\text{SD}⟩$	-6	-4	-2	0	+2	+4	+6
${\Lambda }^{\left[0\right]}$				${\Lambda }^{00}$
${\Lambda }^{\left[2\right]}$				${\Lambda }^{11}$
${\Lambda }^{\left[4\right]}$				${\Lambda }^{22}$
${\Lambda }^{\left[6\right]}$				${\Lambda }^{33}$

Table 5. Table 5: Contributions to the three-body operator O 𝑂 O in normal-ordered form with respect to the quasi-particle vacuum | Φ ⟩ ket Φ |\Phi\rangle and expressed { β , β † } 𝛽 superscript 𝛽 † \{\beta,\beta^{\dagger}\} . The O i ​ j superscript 𝑂 𝑖 𝑗 O^{ij} contributions are sorted horizontally according to i − j 𝑖 𝑗 i-j and vertically according to i + j 𝑖 𝑗 i+j .

$\{\beta ,{\beta }^{†}\},|\Phi ⟩$	-6	-4	-2	0	+2	+4	+6
$O^{\left[0\right]}$				$O^{00}$
$O^{\left[2\right]}$			$O^{02}$	$O^{11}$	$O^{20}$
$O^{\left[4\right]}$		$O^{14}$	$O^{13}$	$O^{22}$	$O^{31}$	$O^{40}$
$O^{\left[6\right]}$	$O^{06}$	$O^{15}$	$O^{24}$	$O^{33}$	$O^{42}$	$O^{51}$	$O^{60}$

Table 6. Table 6: Contributions to the three-body operator O 𝑂 O in normal-ordered form with respect to the slater determinant | SD ⟩ ket SD |\text{SD}\rangle and expressed in { β , β † } 𝛽 superscript 𝛽 † \{\beta,\beta^{\dagger}\} . The O i ​ j superscript 𝑂 𝑖 𝑗 O^{ij} contributions are sorted horizontally according to the value of i − j 𝑖 𝑗 i-j and vertically according to the value of i + j 𝑖 𝑗 i+j .

$\{\beta ,{\beta }^{†}\},|\text{SD}⟩$	-6	-4	-2	0	+2	+4	+6
$O^{\left[0\right]}$				$O^{00}$
$O^{\left[2\right]}$			$O^{02}$	$O^{11}$	$O^{20}$
$O^{\left[4\right]}$		$O^{14}$	$O^{13}$	$O^{22}$	$O^{31}$	$O^{40}$
$O^{\left[6\right]}$	$O^{06}$	$O^{15}$	$O^{24}$	$O^{33}$	$O^{42}$	$O^{51}$	$O^{60}$

Equations

$${\beta }_k$$

$${\beta }_k^{†}$$

$${\rho }_{l_1{l}_2}$$

$${\kappa }_{l_1{l}_2}$$

$${\rho }_{l_2{l}_1}$$

$${\kappa }_{l_2{l}_1}$$

$$A$$

$$O$$

$$+\frac{1}{1!1!}\sum _{l_1{l}_2}{o}_{l_1{l}_2}^{11}{c}_{l_1}^{†}{c}_{l_2}$$

$$+\frac{1}{2!2!}\sum _{l_1{l}_2{l}_3{l}_4}{o}_{l_1{l}_2{l}_3{l}_4}^{22}{c}_{l_1}^{†}{c}_{l_2}^{†}{c}_{l_4}{c}_{l_3}$$

$$+\frac{1}{3!3!}\sum _{l_1{l}_2{l}_3{l}_4{l}_5{l}_6}{o}_{l_1{l}_2{l}_3{l}_4{l}_5{l}_6}^{33}{c}_{l_1}^{†}{c}_{l_2}^{†}{c}_{l_3}^{†}{c}_{l_6}{c}_{l_5}{c}_{l_4}+\dots$$

$$\equiv o^{00}+o^{11}+o^{22}+o^{33}+\dots$$

$$=\sum _{n=0}^N{o}^{nn},$$

$$o_{l_1\dots l_n{l}_{n+1}\dots l_{2n}}^{nn}=ϵ(\sigma )o_{\sigma (l_1\dots l_n|l_{n+1}\dots l_{2n})}^{nn},$$

$$o_{l_1\dots l_n{l}_{n+1}\dots l_{2n}}^{nn\ast }$$

$$O$$

$$\times :c_{l_1}^{†}\dots c_{l_i}^{†}{c}_{l_{i+j}}\dots c_{l_{i+1}}:$$

$$\equiv \sum _{n=0}^N{O}^{[2n]}$$

$${\Lambda }_{l_1\dots l_i{l}_{i+1}\dots l_{i+j}}^{ij}$$

$${\Lambda }_{l_1\dots l_i{l}_{i+1}\dots l_{i+j}}^{ij\ast }$$

$$O^{[2n]}$$

$${\Lambda }^{ij}$$

$$\times :c_{l_1}^{†}\dots c_{l_i}^{†}{c}_{l_{i+j}}\dots c_{l_{i+1}}:.$$

$${\Lambda }_{l_1\dots l_{i+j}}^{ij}$$

$${\Lambda }_{l_1\dots l_{i+j}}^{ij(nn)}$$

$$\times \text{Tr}[o^{nn}\rho {\kappa }^{\ast }\kappa {]}_{l_1\dots l_{i+j}}^{(n_{\rho },n_{{\kappa }^{\ast }},n_{\kappa })},$$

$$\text{Tr}[o^{nn}\rho {\kappa }^{\ast }\kappa {]}_{l_1\dots l_{i+j}}^{(n_{\rho },n_{{\kappa }^{\ast }},n_{\kappa })}$$

$$\times {\kappa }_{l_{i+j+1}{l}_{i+j+2}}^{\ast }\dots {\kappa }_{l_{j+n-n_{\rho }-1}{l}_{j+n-n_{\rho }}}^{\ast }{\kappa }_{l_{j+n+1}{l}_{j+n+2}}\dots {\kappa }_{l_{2n-n_{\rho }-1}{l}_{2n-n_{\rho }}}.$$

$$O$$

$$\times {\beta }_{k_1}^{†}\dots {\beta }_{k_i}^{†}{\beta }_{k_{i+j}}\dots {\beta }_{k_{i+1}}$$

$$\equiv \sum _{n=0}^N{O}^{[2n]},$$

$$O_{k_1\dots k_i{k}_{i+1}\dots k_{i+j}}^{ij}$$

$$O_{k_1\dots k_i{k}_{i+1}\dots k_{i+j}}^{ij\ast }$$

$$O^{[2n]}$$

$$O^{ij}$$

$$\times {\beta }_{k_1}^{†}\dots {\beta }_{k_i}^{†}{\beta }_{k_{i+j}}\dots {\beta }_{k_{i+1}}.$$

$$O^{\text{NO2B}}$$

$$=O^{[0]}+O^{[2]}+O^{[4]},$$

$$\left[O^{\text{nNO2B}},A\right]$$

Full text

1 Introduction

The reach of ab initio many-body theories is rapidly extending over the nuclear chart. In particular, the development of many-body methods based on a systematic expansion of the exact solution of the A-body Schroedinger equation around a conveniently chosen reference state allows an efficient description of medium-mass (semi-)magic nuclei. Examples of such methods are many-body perturbation theory (MBPT) Langhammer:2012jx ; Hu:2016txm ; Tichai:2017rqe ; Tichai:2018mll ; Arthuis:2018yoo ; Hu:2018dza , self-consistent Green’s function (SCGF) Dickhoff:2004xx ; Soma:2011aj ; Soma:2013xha ; Carbone:2013eqa ; Lapoux:2016exf ; Duguet:2016wwr ; Raimondi:2018aa ; Raimondi:2018mtv , coupled-cluster (CC) Hagen:2013nca ; Ja14 ; Signoracci:2014dia ; Morris:2017vxi or in-medium similarity renormalization group (IMSRG) Hergert:2015awm ; Hergert:2016iju ; Parzuchowski:2017wcq ; Morris:2017vxi approaches. To extend their reach to open-shell nuclei and capture noticeably challenging static correlations, these single-reference methods employ a symmetry-breaking reference state Soma:2011aj ; Signoracci:2014dia ; Tichai:2018mll and possibly restore the broken symmetry Duguet:2014jja ; Duguet:2015yle ; Qiu:2018edx .

The computational load is particularly acute due to the relevance of three-nucleon interactions, knowing that including even more demanding four-body forces could eventually be mandatory Bazak:2018qnu . As a matter of fact, dealing with three-nucleon interactions in full makes the solving of the A-body Schrödinger equation rapidly too costly as the mass of the system grows. To circumvent the explicit treatment of three-body operators, ab initio calculations of mid-mass nuclei have been performed on the basis of the so-called normal-ordered two-body (NO2B) approximation. This approximation captures dominant effects of three-nucleon forces, while effectively working with two-body operators. In large-scale no-core shell model (NCSM) calculations, the error induced by the NO2B approximation to the Hamiltonian was estimated to be of the order of $1$-$3\%$ RoBi12 ; Geb16 up to the oxygen region.

The NO2B approximation consists of normal ordering the operator with respect to a many-body reference state and discarding the normal-ordered three-body component. While typically formulated with respect to an uncorrelated reference state, i.e. a Slater determinant, normal-ordering techniques and associated NO2B approximations can be also be devised with respect to a correlated reference state KuMu97 , as is employed in the multi-reference IMSRG approach Hergert:2015awm ; Hergert:2016iju , or even with respect to the fully correlated solution of the problem as is done in SCGF theory Carbone:2013eqa . In any case, the NO2B approximation has been employed so far on the basis of symmetry-conserving states. Only recently such an approximation has been employed in Bogoliubov MBPT (BMBPT) Tichai:2018mll ; Tichai:2019ksh in which U(1) symmetry associated with particle-number conservation is spontaneously broken by the (approximate) many-body state. In this context, the normal ordering of operators at play is performed with respect to a particle-number-breaking Bogoliubov reference state such that proceeding to a naive truncation may lead to approximating a particle-number-conserving operator by a particle-number-breaking one. A similar situation shall occur when using a SU(2)-breaking, i.e. deformed, reference state such that proceeding to a naive normal-ordered truncation of a rotationally invariant operator may lead to an angular-momentum non-conserving approximation.

The purpose of the present work is thus to design a general normal-ordering approximation of operators that is consistent with symmetries of the Hamiltonian while working with symmetry broken (and restored) reference state111The symmetry-breaking nature of the many-body methods of present interest concerns the approximate many-body state while making use of a symmetry-conserving Hamiltonian. Approximating the Hamiltonian in a symmetry-violating way is conceptually different and more problematic as one wishes to eventually employ a symmetry-restored approximate many-body state Duguet:2014jja ; Duguet:2015yle ; Qiu:2018edx .. Focusing on many-body formalisms in which U(1) symmetry associated with particle number conservation is broken (and potentially restored), a particle-number-conserving normal-ordered $k$-body (PNOkB) approximation of an arbitrary N-body operator is designed on the basis of Bogoliubov reference states. A numerical test based on particle-number projected Hartree-Fock-Bogoliubov calculations is designed and employed to check the particle-number conserving/violating character of the approximate operator.

2 Formalism

2.1 Bogoliubov reference state

Many-body formalisms of present interest, e.g. BMBPT Duguet:2015yle ; Tichai:2018mll ; Arthuis:2018yoo , Bogoliubov coupled cluster (BCC) Signoracci:2014dia , Gorkov SCGF Soma:2011aj ; Soma:2013xha , Projected BMBPT (PBMBPT) Duguet:2015yle and Projected BCC (PBCC) Duguet:2015yle ; Qiu:2018edx , rely on the use of a particle-number breaking Bogoliubov reference state $|\Phi\rangle$.

The state $|\Phi\rangle$ is a vacuum for a complete set of quasi-particle operators $\{\beta_{k},\beta^{\dagger}_{k}\}$ obtained from an arbitrary single-particle basis $\{c_{l},c^{\dagger}_{l}\}$ of the one-body Hilbert space through a unitary Bogoliubov transformation

[TABLE]

such that $\beta_{k}|\Phi\rangle=0$ for all $k$. Although not presently concerned with these technical details, note that the columns of the transformation matrices $(U,V)$ are typically obtained as eigenvectors of the Hartree-Fock-Bogoliubov (HFB) eigenvalue problem RiSc80 . More importantly for the discussion given below, $|\Phi\rangle$ reduces to a Slater determinant, presently denoted as $|\text{SD}\rangle$, whenever U(1) symmetry is not spontaneously broken, e.g. in closed shell systems where the HFB eigenproblem typically reduces to the Hartree-Fock (HF) one.

Normal and anomalous elementary contractions of single-particle creation and annihilation operators are defined with respect to $|\Phi\rangle$ by

[TABLE]

such that

[TABLE]

2.2 Particle-number conserving operator

We are typically concerned with the treatment of a N-body operator $O$ that commutes with the particle-number operator

[TABLE]

and whose general form in the basis $\{c_{l},c^{\dagger}_{l}\}$ is

[TABLE]

where matrix elements $o^{nn}_{l_{1}\ldots l_{n}l_{n+1}\ldots l_{2n}}$ are fully antisymmetric with respect to the permutation of the $n$ first, resp. $n$ last, indices

[TABLE]

where $\epsilon(\sigma)$ refers to the signature of the permutation $\sigma$. The notation $\sigma(\ldots|\ldots)$ denotes a separation between the $n$ first and the $n$ last indices such that permutations are only considered between members of the same group. If $O$ is hermitian, the n-body matrix elements $o^{nn}_{l_{1}\ldots l_{n}l_{n+1}\ldots l_{2n}}$ satisfy

[TABLE]

A graphical representation of the contributions to a three-body operator $O$ is given in Tab. 1 where the $o^{ij}$ contributions are sorted horizontally according to $i-j$ and vertically according to $i+j$. As $O$ commutes with $A$, there is no contribution with $i\neq j$ in the present case. This graphical representation will become useful below when normal-ordering the operator with respect to the quasi-particle vacuum $|\Phi\rangle$.

2.3 Normal ordering

In order to address (in)appropriate approximations to the operator $O$, its normal-ordered forms with respect to $|\Phi\rangle$ expressed in two different bases are needed.

2.3.1 Single-particle basis

The normal ordering of $O$ with respect to the Bogoliubov vacuum $|\Phi\rangle$ leads to re-expressing the operator under the form

[TABLE]

where $:\ldots:$ denotes the normal-ordered product with respect to $|\Phi\rangle$. Matrix elements $\Lambda^{ij}_{l_{1}\ldots l_{i+j}}$ are fully antisymmetric with respect to the permutation of the $i$ first, resp. $j$ last, indices

[TABLE]

If $O$ is hermitian, field matrix elements $\Lambda^{ij}_{l_{1}\ldots l_{i}l_{i+1}\ldots l_{i+j}}$ satisfy

[TABLE]

The class $O^{[2n]}$ groups all the terms containing a normal-ordered product of $2n$ single-particle operators, i.e. terms possibly containing different numbers of single-particle creation and annihilation operators according to

[TABLE]

where $\Lambda^{ij}$ is a $n$-body normal, resp. anomalous, field if $i=j$, resp. $i\neq j$, containing all the terms with a normal-ordered product of $i$, resp. $j$, single-particle creation, resp. annihilation, operators and reading as

[TABLE]

Applying standard Wick’s theorem Wi50 to $O$, $\Lambda^{ij}$ receives contributions from all $n$-body terms $o^{nn}$ with $n\geq\max(i,j)$, i.e.

[TABLE]

such that $\Lambda^{ij(nn)}_{l_{1}\ldots l_{i+j}}$ accounts for all appropriate contraction patterns and is given by

[TABLE]

where the notation $\text{Tr}[o^{nn}\rho\kappa^{\ast}\kappa]^{(n_{\rho},n_{\kappa^{\ast}},n_{\kappa})}_{l_{1}\ldots l_{i+j}}$ denotes a trace over the indices of the $n_{\rho}$ normal contractions $\rho$, the $n_{\kappa^{\ast}}$ anomalous contractions $\kappa^{\ast}$ and the $n_{\kappa}$ anomalous contractions $\kappa$, i.e. {strip}

[TABLE]

Whenever the Bogoliubov vacuum reduces to a Slater determinant, all anomalous contractions are null, i.e. $\kappa=\kappa^{\ast}=0$, such that $\Lambda^{ij}=0$ for $i\neq j$. The proof leading to Eqs. (2.3.1-15), together with the explicit form of the $\Lambda^{ij}$ matrix elements associated with a three-body operator, are given in App. A.1.

2.3.2 Quasi-particle basis

One can rather choose to express the normal-ordered form of $O$ in the quasi-particle basis $\{\beta_{k},\beta^{\dagger}_{k}\}$. Applying standard Wick’s theorem leads to

[TABLE]

where matrix elements $O^{ij}_{k_{1}\ldots k_{i+j}}$ are fully antisymmetric with respect to the permutation of the $i$ first, resp. $j$ last, indices

[TABLE]

If $O$ is hermitian, matrix elements $O^{ij}_{k_{1}\ldots k_{i+j}}$ satisfy

[TABLE]

The class $O^{[2n]}$ groups all terms containing a normal-ordered product of $2n$ quasiparticle operators that can be further differentiated according to the number of quasiparticle creation and annihilation operators they contained, i.e.

[TABLE]

where $O^{ij}$ gathers all terms with $i$, resp. $j$, quasiparticle creation, resp. annihilation, operators and reads

[TABLE]

Extending to $O^{[6]}$ the results provided in Ref. Signoracci:2014dia , the explicit form of the $O^{ij}$ matrix elements associated with a three-body operator $O$ is given in App. A.2 as an example.

2.4 Approximation

The need to handle (at least) three-nucleon interactions makes the solving of the A-body Schrödinger equation particularly cumbersome beyond the lightest nuclei. At the mean-field level one is able to handle the full Hamiltonian. Beyond the mean-field, however, three-body interactions prevent the implementation of many-body methods in their full glory.

To tackle this difficulty, the so-called normal-ordered two-body (NO2B) approximation of three-body operators has been designed and used with success RoBi12 ; Geb16 . This approximation amounts to employing the full Hamiltonian to generate the mean-field Slater determinant reference state $|\text{SD}\rangle$ and to truncating the normal-ordered form of $H$ with respect to $|\text{SD}\rangle$ at the effective two-body level for the beyond mean-field step. Thus, the NO2B approximation of a three-body operator $O$ commuting with $A$ reads as

[TABLE]

where all terms containing a normal-ordered product of six single-particle operators have been discarded. A graphical representation of this approximation is given in Fig. 1.

2.4.1 Naive extension

Whereas this truncation is transparent and appropriate in the symmetry-conserving context, it has not been discussed in connection with symmetry-breaking reference states. It happens that, even if $O$ commutes with $A$, the approximate operator $O^{\text{NO2B}}$ obtained via a naive generalization of Eq. (2.4) on the basis of a Bogoliubov reference state $|\Phi\rangle$, does not commute with $A$, i.e.

[TABLE]

A graphical representation of this naive extension of the NO2B approximation is given in Fig. 2.

For $O=H$ the naive extension of the NO2B approximation would imply working with a Hamiltonian $H^{\text{nNO2B}}$ whose exact eigenstates are not eigenstates of $A$. While this can be considered as part of a systematic error induced by the approximation, it eventually leads, when restoring the symmetry of the approximate wave function, to energies that problematically depend on the particular way the symmetry restoration is formulated. This key feature is discussed and illustrated later on.

In this context, it is desirable to generalize the NOkB approximation of an arbitrary N-body operator in such a way that the truncated operator is particle-number conserving even if the reference state breaks $U(1)$ symmetry.

2.4.2 Particle-number-conserving NOkB approximation

In the following, $O$ denotes a normal-ordered $N$-body operator such that its naive NOkB (nNOkB) approximation reads as

[TABLE]

The nNOkB approximation suffers from the same pathology as the naive NO2B approximation regarding particle-number violation. One now wishes to design an approximation to $O$ fulfilling the following three requirements

All normal-ordered terms of ranks higher than $k$ must be discarded as the practical goal of normal-ordered approximations is to eventually work with an effective operator characterized by a maximum rank $k<N$. 2. 2.

The approximate operator must commute with $A$. 3. 3.

The error generated by the approximation must be minimal.

While the nNOkB approximation fulfills the first condition, it violates the second by fully retaining the normal-ordered contributions with ranks lower or equal than $k$. One can thus anticipate that fulfilling the second condition in addition to the first one requires to further approximate specific parts of the operator displayed in Eq. 23. In the following, a systematic procedure to achieve this goal is devised for any $N$ and any $k<N$. It must be noted that there is no unique way to do so222In App. B.3, an alternative particle-number conserving approximation based on the quasi-normal ordering proposed in Ref. kong10a is investigated.. Eventually, while the third condition may be anticipated based on reasonable arguments, it can solely be validated through benchmark calculations.

The particle-number-conserving NOkB (PNOkB) approximation to $O$ is given by

[TABLE]

which is manifestly particle-number conserving. The $n$-body matrix elements $\tilde{o}^{nn}_{l_{1}\ldots l_{2n}}$, $n\leq k$, are recursively defined in decreasing order, from $n=k$ down to $n=0$, by

[TABLE]

where $\tilde{\Lambda}^{nn(mm)}$, $m>n$, is the $m$-body contribution to the $n$-body normal field associated with $O^{\text{PNOkB}}$ defined in the same way as the $m$-body contribution to the $n$-body normal field $\Lambda^{nn(mm)}$ associated with the full operator $O$ was introduced in Eq. (2.3.1).

It is clear that $O^{\text{PNOkB}}$ contains, through the fields $\Lambda^{nn}$, $n\leq k$, contributions from $o^{nn}$ up to $n=N$. It is the usual benefit of a normal-ordered approximation to capture the dominant part of all contributions to the original operator while working with effective operators of lower ranks.

To better appreciate the content of $O^{\text{PNOkB}}$, let us consider its normal-ordered form in the single-particle basis. As proven in App. B, $O^{\text{PNOkB}}$ is obtained from $O$ via a two-step process, i.e. by

removing all $\Lambda^{ij}$ fields with $\max(i,j)>k$, 2. 2.

adding $\breve{\Lambda}^{ij}$ defined below for $\max(i,j)\leq k$.

This leads to rewriting the approximate operator as

[TABLE]

with

[TABLE]

where the extra term $\breve{\Lambda}^{ij}$ is given by {strip}

[TABLE]

for $i>j$, by

[TABLE]

for $i<j$ and by

[TABLE]

for $i=j$. In the end, $\breve{\Lambda}^{ij}$ is non-zero only if $i-j\geq 2$ and $n_{\kappa}\geq 1$, or $j-i\geq 2$ and $n_{\kappa^{\ast}}\geq 1$. This is only possible for an initial $N$-body operator $O$ with $N\geq 4$, i.e. up to an initial three-body operator, no extra term is to be considered and the retained fields $\tilde{\Lambda}^{ij}$ are trivially given by the original fields $\Lambda^{ij}$.

The leading criterion behind the definition provided in Eqs. (24-2.4.2) is to ensure that $O^{\text{PNOkB}}$ shares the same normal fields as $O$ for $n\leq k$, i.e.

[TABLE]

as those are believed to be the dominant contributions to the normal-ordered operator. This requirement possibly induces extra contributions $\breve{\Lambda}^{ij}$ to the anomalous fields. This feature relates to the fact that $\tilde{\Lambda}^{ij}$ is obtained from the original operator through two successive normal orderings steps rather than one for $\Lambda^{ij}$. When proceeding to anomalous contractions in the application of Wick’s theorem, it possibly leads to the non-equivalent combinatorial prefactors

[TABLE]

in both procedures.

Eventually, and as best illustrated through the examples worked out in App. B, the effect of the correction terms defined through Eqs. (28-30) is nothing but to modify the numerical prefactors of specific contributions to the normal-ordered anomalous fields expressed in the single-particle basis. Noticing that plainly omitting $\Lambda^{ij}$ with $\max(i,j)>k$ can also be viewed as the mere replacement of its original prefactor by an approximate one (0 in such cases), the function of the extra terms $\breve{\Lambda}^{ij}$ is in fact not different, i.e. it corresponds to modifying the prefactors of a specific set of fields such that the initial objectives are fulfilled. The quantitative performance of the procedure can only be judged a posteriori by comparing many-body results obtained using the full operator and its PNOkB approximation.

2.4.3 PNO1B approximation of a two-body operator

Let us exemplify the PNOkB procedure by building the PNO1B approximation of the two-body operator

[TABLE]

This constitutes the simplest possible case. While it is not of interest for realistic ab initio calculations, it is used throughout the present paper to illustrate the difference between a nNOkB and the newly designed PNOkB.

The normal ordering of $O$ with respect to $|\Phi\rangle$ reads in the single-particle basis as

[TABLE]

and involves the following normal and anomalous fields

[TABLE]

where each contribution, in fact denoting matrix elements of the fields, is written in a compact form without indices such that the key informations, i.e. the original $o^{nn}$ term, the set of elementary contractions involved and the overall prefactor, can be easily extracted.

The nNO1B approximation leads to just dropping $O^{[2]}=\Lambda^{22}$ in Eq. (33). Using properties of elementary contractions

[TABLE]

one can write

[TABLE]

where the matrix elements are given by

[TABLE]

Equations (36-37) explicitly demonstrate that the nNO1B approximation leads to an operator that does not conserve particle number as it contains non-zero two-particle-addition $\tilde{o}^{20}$ and two-particle removal $\tilde{o}^{02}$ contributions.

The application of the PNO1B approximation is more involved and we now proceed to the construction of the corresponding operator

[TABLE]

where the terms are to be obtained recursively on the basis of Eq. (2.4.2). The procedure starts with $\tilde{o}^{kk}=\Lambda^{kk}$, which is normal-ordered in the single-particle basis to generate $\tilde{o}^{(k-1)(k-1)}$ etc. until $\tilde{o}^{00}$ is reached. In the present case, the procedure is trivial since it starts with $k=1$ and leads directly to $\tilde{o}^{00}$ in the first step. Eventually, the result reads as

[TABLE]

where the only needed intermediate quantity is

[TABLE]

One first observes that, even if $O$ has no constant term $o^{00}$ to begin with, $O^{\text{PNO1B}}$ does acquire one. The corresponding contributions originate from the two-body operator whose normal-ordered two-body part is omitted. When normal ordering $O^{\text{PNO1B}}$ in the single-particle basis, this term combines with the fully contracted part obtained from $\tilde{o}^{11}$ to generate the complete fully contracted part of the original two-body operator. One further notes that $\tilde{o}^{00}$ entering $O^{\text{PNO1B}}$ is different from the one appearing in Eq. (37), which underlines the fact that the procedure to obtain the PNO1B does not correspond to performing the nNO1B before dropping particle-number non-conserving terms.

Equations (38-39) fully define the PNO1B approximation of the two-body operator. It is particle-number conserving by construction. It is useful to further characterize the operator by closely inspecting its normal-ordered contributions. Focusing for example on the normal field $\tilde{\Lambda}^{11}$, one obtains

[TABLE]

which indeed satisfies the systematic property $\tilde{\Lambda}^{ii}=\Lambda^{ii}$. Eventually, the complete set of normal-ordered contributions to $O^{\text{PNO1B}}$ expressed in the single-particle basis relates to those of the original two-body operator through

[TABLE]

One thus observes that the non-zero fields are strictly equal to those associated with the original operator, i.e. there is no extra term when approximating a two-body operator.

A graphical representation of the PNO1B approximation of a two-body operator $O$ is given in Fig. 3 .

2.4.4 Other applications

The PNOkB approximation procedure is worked out in details for three other cases in App. B.2, i.e.

the PNO2B approximation of a 3-body operator, 2. 2.

the PNO2B approximation of a 4-body operator, 3. 3.

the PNO3B approximation of a 4-body operator,

The first case applies to state-of-the-art ab initio calculations Tichai:2018mll ; Tichai:2019ksh . The second case illustrates how the extra terms $\breve{\Lambda}^{ij}$ come into play for the first time and may become of practical interest in the future in case four-nucleon interactions have to be accounted for at the two-body effective level. If one becomes capable of working at the effective three-body level, the third case provides the PNO3B approximation of a 4-body operator.

Finally, the main outcomes of the PNO3B and PNO4B approximations of a 5-body operator are briefly compared in App. B.2 to illustrate how the extra terms depend on the rank $k$ of the approximate operator.

2.5 Testing $[O^{\text{(P)NOkB}},A]=0$

A method is now introduced to test (formally and numerically) whether or not a given operator $F$ is commuting with the particle-number operator $A$. Starting from a particle-number-conserving operator $O$, the test is meant to be applied to both333The test constitutes a necessary but not sufficient condition to prove that $F$ commutes with $A$, i.e. the test can thus be used to prove that $O^{\text{nNOkB}}$ does not commute with $A$ but can only indicate that $O^{\text{PNOkB}}$ probably commutes with $A$. This limitation is not problematic in the present case given that $O^{\text{PNOkB}}$ does commute with $A$ by construction. $F\equiv O^{\text{nNOkB}}$ and $O^{\text{PNOkB}}$.

2.5.1 Particle-number projection

The test is based on the fact that, if $F$ commutes with $A$, it must also commute with the particle-number projection operator $P^{\text{A}}$ defined as

[TABLE]

where $R(\varphi)\equiv e^{iA\varphi}$ is the rotation operator in gauge space. Given that $P^{\text{A}}$ is idempotent, i.e. $(P^{\text{A}})^{2}=P^{\text{A}}$, the commutation of $F$ with $P^{\text{A}}$ can be re-expressed as

[TABLE]

Based on the above, the ratio of singly over doubly projected mean-field matrix elements

[TABLE]

is formed such that the particle-number conserving (violating) character of $F$ corresponds to $Q_{F}^{\text{A}}=1(\neq 1)$.

2.5.2 Computation of $Q_{F}^{\text{A}}$

The computation of $Q_{F}^{\text{A}}$ relies on the representation of $P^{\text{A}}$ given in Eq. (43) such that

[TABLE]

where singly- and doubly-rotated mean-field off-diagonal norm and operator kernels are given by

[TABLE]

and by

[TABLE]

where the gauge-rotated Bogoliubov state is defined through $\langle\Phi(\varphi)|\equiv\langle\Phi|R(\varphi)$.

2.5.3 One-body operator

The numerical illustrations provided in the present paper rely on the (P)NO1B of a two-body nuclear Hamiltonian, i.e. $O\equiv H$. Let us thus characterize $Q_{F}^{\text{A}}$ for a generic, i.e. possibly particle-number violating, one-body operator

[TABLE]

A graphical representation of the operator $F$ is given in Fig. 2. Starting from a two-body operator $O$, the explicit expression of the matrix elements $(f^{00},f^{20}_{l_{1}l_{2}},f^{11}_{l_{1}l_{2}},f^{02}_{l_{1}l_{2}})$ associated to $O^{\text{nNO1B}}$ were provided in Sec. 2.4.3, where it was formally proven that the particle-number non-conserving terms $f^{20}$ and $f^{02}$ are indeed non zero in this case. For $O^{\text{PNO1B}}$, the sole non-zero terms are $f^{00}$ and $f^{11}$.

By virtue of the off-diagonal Wick’s theorem BaBr69 , singly- and doubly-rotated connected operator kernels write respectively as

[TABLE]

where singly and doubly gauge-rotated contractions are defined and calculated in App. C.2. Using the change of variables $\phi\equiv\varphi-\varphi^{\prime}$, $f^{(0)}(\varphi,\varphi^{\prime})$ can be expressed

[TABLE]

i.e. it displays a Fourier decomposition in $\varphi^{\prime}$ whose components are nothing but the singly gauge-rotated kernels associated with the various contributions to $F$. The appearance of non-trivial Fourier components, i.e. irreducible representation of $U(1)$, constitute a fingerprint of the particle-number non-conserving character of $F$. In the present case, two such non-trivial modes appear in connection with $f^{20}$ and $f^{02}$. This result can obviously be extended to higher-body operators. Considering a general operator $F$ containing arbitrary combinations of (an even number of) single-particle creation and annihilation operators, its doubly gauge-rotated mean-field connected kernel takes the form

[TABLE]

such that the Fourier component labelled by $k=m-n$ receives contributions from all $k$-particle addition operators contributing to $F$.

With Eqs. (49-50) at hand, and as demonstrated in App. C.3, one obtains

[TABLE]

which is (a priori not) equal to $1$ when the particle non-conserving parts of $F$ are (non) vanishing.

3 Results

3.1 Calculations set up

The nuclear Hamiltonian used in this work includes a chiral two-nucleon (2N) interaction at next-to-next-to-next-to leading order with a cutoff of $\Lambda_{2N}=500\,\text{MeV}$ EnMa03 . The three-nucleon interaction is omitted to be able to compare results obtained with the full Hamiltonian to those obtained via its nNO1B and PNO1B approximations. The Hamiltonian is further softened using a Similarity Renormalization Group (SRG) transformation with a flow parameter $\alpha=0.08\,\text{fm}^{4}$ BoFu07 ; HeRo07 ; RoRe08 ; RoLa11 ; JuMa13 such that only up to transformed two-body operators are retained.

Calculations are performed using the one-body eigenbasis of the spherical harmonic oscillator (HO) Hamiltonian with frequency $\hbar\Omega=20$ MeV. One- and two-body operators are represented using all single-particle states up to $e_{\text{max}}=(2n+l)_{\text{max}}=10$. Spherical Hartree-Fock-Bogoliubov (HFB) calculations are performed in J-coupled scheme while the (single and double) particle-number projection (PHFB) is performed after the variation (PAV) on the basis of $n_{{\text{i}nt}}=500$ equally-spaced integration points over the interval444Working with a Bogoliubov state carrying good, e.g. even, number-parity quantum number, the integration over the gauge angles can indeed be reduced to the interval $[0,\pi]$ in Eq. (2.5.2). $\varphi\in[0,\pi]$.

3.2 Energetics in O isotopes

Before coming to normal-ordered approximations, let us discuss the ground-state energetics obtained with the full Hamiltonian to set the orders of magnitude at play. Since three-body forces are presently discarded for the sake of the demonstration, computed energies are not meant to reproduce experimental data.

In Fig. 4, ground-state energetics along the Oxygen isotopic chain are displayed. Panel (a) provides HFB energies ranging from about $-100$ MeV in 14O to about $-220$ MeV in 26O. Panel (b1) displays static correlations associated with particle number projection that are of the order of555In energy density functional (EDF) calculations, particle-number projection typically lowers the energy, which is not the case with the presently used chiral 2N interaction. Because PHFB calculations with realistic nuclear interactions are rather novel, this feature is briefly analyzed in App. D by looking at the decomposition of the HFB vacuum into its particle-number projected components. $\pm 1$MeV. This order of magnitude is to be compared with dynamical correlations displayed in panel (b2). Calculated via, e.g., BMBPT(2) Duguet:2015yle ; Tichai:2018mll on the basis of the same Hamiltonian and model space, dynamical correlations range from $-26$ to $-38$ MeV in O isotopes. Eventually, both types of correlations can be captured consistently via PBMBPT Duguet:2015yle ; ripoche19b or PBCC Duguet:2015yle ; Qiu:2018edx .

With results from the full Hamiltonian at hands, one can now analyze the effect of approximating it. Panel (c1) of Fig. 4 displays the difference of PHFB energies obtained with the two-body Hamiltonian and with its nNO1B or PNO1B approximations. While results are the same for single or double particle-number projection when using $H^{\text{PNO1B}}$, it is not the case for $H^{\text{nNO1B}}$. This observation is confirmed in Panel(c2) where the ratios $Q_{H^{\text{nNO1B}}}^{\text{A}}$ and $Q_{H^{\text{PNO1B}}}^{\text{A}}$ are displayed. While $Q_{H^{\text{PNO1B}}}^{\text{A}}=1$ for all isotopes, $Q_{H^{\text{nNO1B}}}^{\text{A}}\neq 1$ except in doubly closed-shell isotopes. This key result proves that $H^{\text{nNO1B}}$ is particle-number violating, thus making PHFB energies dependent on the way the particle-number projection is defined. Contrarily, $H^{\text{PNO1B}}$ does behave as a particle-number conserving operator.

One further observes that the effect of the PNO1B approximation is significant, i.e. it is of the same order as the effect of the projection itself displayed in Panel (b1). This is not surprising given that ignoring the residual part of the two-body interaction has never believed to be an appropriate approximation. In order to quantitatively gauge the quality of the PNOkB approximations, one must at least test the PNO2B approximation of a three-nucleon interaction, which is beyond the scope of the present paper that is rather focusing on the symmetry violating/conserving character of a given approximation method.

3.3 Doubly gauge-rotated operator kernel

To further characterize the normal-ordered approximations, the doubly gauge-rotated mean-field connected Hamiltonian kernel $h^{(0)}(\varphi,\varphi^{\prime})$ is analyzed along with those associated with $H^{\text{nNO1B}}$ and $H^{\text{PNO1B}}$. The real and imaginary parts of the kernels are shown in Fig. 5 for 18O. They are displayed as contour plots with respect to the variables $(\phi=\varphi-\varphi^{\prime},\varphi^{\prime})$.

Because $H$ is particle-number conserving, its kernel solely depends on $\phi$ and is independent of $\varphi^{\prime}$, i.e. the Fourier expansion of $h^{(0)}(\varphi,\varphi^{\prime})$ with respect to $\varphi^{\prime}$ in Eq. (51) only contains the trivial component, i.e. the irreducible representation of $U(1)$, characterized by $k=0$. While this feature is already manifest in the upper panels of Fig. 5, it is confirmed in Fig. 6 where the Fourier components are numerically extracted.

Performing the nNO1B approximation, the kernel displayed in the middle panels of Fig. 5 is obtained. It is clear that both the real and imaginary parts now vary with $\varphi^{\prime}$. This indicates that non-trivial components are present in the Fourier decomposition of Eq. (50) due to particle-number non-conserving contributions to $H^{\text{nNO1B}}$. As visible from the middle panel of Fig. 6 non-zero Fourier components confirm the presence of $f^{20(0)}(\phi)$ and $f^{02(0)}(\phi)$, in agreement with the analytical derivation provided in Eqs. (36-37).

Moving from $H^{\text{nNO1B}}$ to $H^{\text{PNO1B}}$, the kernel displayed in the bottom panels of Fig. 5 are obtained. The independence of the kernel with respect to $\varphi^{\prime}$ is recovered, thus testifying of the particle-number conserving nature of the approximate Hamiltonian $H^{\text{PNO1B}}$. This is confirmed in the lower panel of Fig. 6.

3.4 Systematics

The analysis provided in Sec. 3.2 is now extended to Ca isotopes and to the variance operator. The energetics provided in Fig. 7 confirms the conclusions drawn earlier, i.e. while PHFB energies in O and Ca isotopes originating from $H^{\text{PNO1B}}$ are identical when using single or double projection, it is not the case for $H^{\text{nNO1B}}$. It confirms that $H^{\text{nNO1B}}$ ($H^{\text{PNO1B}}$) is particle-number violating (conserving).

As a curiosity, nNO1B and PNO1B approximations are further applied to the particle-number variance operator. Corresponding results are displayed in Fig. 8. When using the full operator, PHFB calculations obviously deliver a null variance, independently of whether the single or double projection is employed. Next, results obtained via the approximate nNO1B (PNO1B) one-body operator do (not) depend on the way the projection is performed, thus confirming that the approximate operator is particle-number violating (conserving).

Furthermore, results do depart significantly from the correct one in the PNO1B approximation, which is again not surprising given the expected crudeness of the one-body approximation to a two-body operator in general. Interestingly, nNO1B delivers essentially identical results to PNO1B when the single projection is employed. More surprisingly, nNO1B does provide essentially exact (i.e. null) results when the double projection is used. The reasons for this unexpected result are analyzed in App. C.4 where they are shown to be specific to the one-body approximation of the particle-number variance operator and thus to be accidental.

4 Conclusion

In view of dealing efficiently with three-, possibly four-, nucleon interactions in ab initio calculations of open-shell nuclei, the present paper addresses approximations based on normal-ordering techniques in the context of many-body methods in which the exact solution of the A-body Schrödinger equation is expanded around a symmetry-breaking reference state. Because a naive extension of the standard normal-ordered approximation, originally designed on the basis of symmetry-conserving states, may lead to symmetry-breaking approximate operators, a systematic approximation techniques delivering symmetry-conserving operators is necessary in this context.

Focusing here on many-body formalisms in which U(1) global-gauge symmetry associated with particle number conservation is broken (and potentially restored), a particle-number-conserving normal-ordered $k$-body (PNOkB) approximation of an arbitrary N-body operator has been designed on the basis of Bogoliubov reference states. After laying down the general formalism, the explicit form of the approximate operator has been provided for various relevant combinations of N and $k$. Furthermore, numerical tests based on particle-number projected Hartree-Fock-Bogoliubov calculations have allowed to check (i) the particle-number violating character of a naive extension of the standard normal-ordered approximation and (ii) the particle-number conserving character of the newly designed PNOkB approximation. Using the PNOkB approximation, ab initio calculations based on formalisms exploiting the breaking and restoration of particle-number can thus be safely performed. The future formulation of an angular-momentum-conserving normal-ordered $k$-body approximation based on deformed Slater determinant or Bogoliubov reference states is envisioned.

5 Acknowledgements

The authors thank H. Hergert for pointing out to them the existence of the quasi-normal ordering of Ref. kong10a .

Appendix A Normal-ordered operator

A.1 Particle basis

A.1.1 Derivation

Applying Wick’s theorem Wi50 to the n-body part of $O$ gives

[TABLE]

The n-body matrix elements $o^{nn}_{l_{1}\ldots l_{n}l_{n+1}\ldots l_{2n}}$ are fully antisymmetric with respect to the permutation of the n first, resp. n last, indices, so that all terms coming from $n_{\rho}$ normal $\rho$ contractions can be recasted into a single term with an appropriate numerical factor, and similarly for all terms coming from $n_{\kappa^{\ast}}$ anomalous $\kappa^{\ast}$ contractions or $n_{\kappa}$ anomalous $\kappa$ contractions. For the term with $n_{\rho}$ normal $\rho$ contractions, there are $n_{\rho}$ creation, resp. annihilitation, operators among $n$ to be considered. Furthermore there are $n_{\rho}!$ ways to pair $n_{\rho}$ creation and $n_{\rho}$ annihilation operators. Thus, the numerical factor resulting from performing $n_{\rho}$ normal $\rho$ contractions and no anomalous contraction of an initial n-body operator is

[TABLE]

Adding $n_{\kappa^{\ast}}$, resp. $n_{\kappa}$, anomalous $\kappa^{\ast}$, resp. $\kappa$ contractions, the numerical factor becomes {strip}

[TABLE]

The term whose prefactor is given by Eq. (A.1.1) contributes to $\Lambda^{ij}$ with $i=n-n_{\rho}-2n_{\kappa^{\ast}}$ and $j=n-n_{\rho}-2n_{\kappa}$. This prefactor accounts for the factor $1/i!j!$ involved in the definition of the operator $\Lambda^{ij}$ in Eq. (12) and for the factor appearing in Eq. (2.3.1).

A.1.2 Example

Taking a three-body operator $O$ commuting with $A$ as an example, the matrix elements $\Lambda^{ij}_{l_{1}\ldots l_{i}l_{i+1}\ldots l_{i+j}}$ are given by

[TABLE]

The graphical representation of $O$ expressed in normal-ordered form with respect to $|\Phi\rangle$ is given in Fig. 3. Contrarily to the case displayed in Fig. 1, non-diagonal contributions $\Lambda^{20}$, $\Lambda^{02}$, $\Lambda^{31}$ and $\Lambda^{13}$ appear.

In the Slater determinant limit, for which $\rho_{ip}=\delta_{ip}$, $\rho_{ap}=0$ and $\kappa_{pq}=0$ 666Labels $i,j,\ldots$ denotes hole indices, i.e. occupied single-particle states in the Slater determinant, whereas $a,b,\ldots$ denote particle indices, i.e. occupied single-particle states in the Slater determinant. Indices $p,q,\ldots$ denotes either of those., the set of matrix elements reduce to

[TABLE]

The graphical representation of $O$ in normal-ordered form with respect to $|\text{SD}\rangle$ is given in Fig. 4. In this case, no non-diagonal contribution appears.

A.2 Quasi-particle basis

The derivation relates to the application of Wick’s theorem with respect to $|\Phi\rangle$, which is particularly straighforward using the set of quasi-particle operators of which $|\Phi\rangle$ is a vacuum. This procedure has already been discussed and exemplified in Ref. Signoracci:2014dia for a three-body operator $O$ commuting with $A$. However, the expressions of the matrix elements $O^{ij}_{k_{1}\ldots k_{i}k_{i+1}\ldots k_{i+j}}$ were only provided for $O^{[0]}$, $O^{[2]}$ and $O^{[4]}$. Extending for completeness the results of Ref. Signoracci:2014dia , matrix elements up to $O^{[6]}$ are given by {strip}

[TABLE]

The graphical representation of the normal-ordered operator $O$ with respect to $|\Phi\rangle$, resp. $|\text{SD}\rangle$, and expressed in the quasi-particle basis, is given in Tab. 5, resp. Tab 6. One notices the presence of non-diagonal contributions in Tab. 5 (as in Tab. 3) and in Tab 6 (even if there is none in Tab. 4). This is due to the use of the quasi-particle algebra that is agnostic regarding the symmetry-breaking or the symmetry-preserving character of the vacuum.

Appendix B Particle-number-conserving NOkB approximation

B.1 Derivation

For each $(i,j)$ such that $\max(i,j)\leq k$ one has to prove that

[TABLE]

The n-body normal fields, for which $i=j=n$, are given by

[TABLE]

where the property $\tilde{\Lambda}^{nn(nn)}=\tilde{o}^{nn}$ and the definition of $\tilde{o}^{nn}$ in Eq. (2.4.2) have been used. Consequently, $\breve{\Lambda}^{nn}=0$ (30) as required.

Let us now focus on the n-body contributions to anomalous fields associated with the full operator $O$ and for which $n\geq\max(i,j)$ and, e.g., $i>j$. One can write {strip}

[TABLE]

where $\Lambda^{ij(nn)(n_{\kappa})}$ denotes the contribution from $o^{nn}$ to $\Lambda^{ij}$ containing $n_{\kappa}$ anomalous contractions $\kappa$. Indeed, several contractions patterns associated with different numbers of normal and anomalous contractions can lead from $o^{nn}$ to $\Lambda^{ij}$, knowing that the minimal number of $\kappa$ contractions is $(i-j)/2$. The goal of the above rewriting is to factorize $(i-j)/2$ anomalous contractions $\kappa$ in order to express each contribution $\Lambda^{ij(nn)(n_{\kappa})}$ in terms of its diagonal (i.e. $i=j$) partner $\Lambda^{ii(nn)(n_{\kappa})}$. Summing (A.34) over n-body contributions with $n=i,\ldots,N$, one obtains

[TABLE]

A similar relation holds for $\tilde{\Lambda}^{ij}$ associated with $O^{\text{PNOkB}}$

[TABLE]

the difference being that the sum over $n$ only extends up to $k$ instead of $N$ for $\Lambda^{ij}$ given that $\tilde{o}^{nn}$ is null for $n>k$ to begin with. Using Eq. (B.1), the combination of the two above identities allows one to relate both sets of fields through

[TABLE]

with

[TABLE]

where $\tilde{\Lambda}^{ii(nn)(n_{\kappa})}$ is in fact zero for $n>k$. The case of anomalous field with $i<j$ is obtained through the same procedure by factorizing $\kappa^{\ast}$ contractions instead of $\kappa$ ones.

Of course, the standard NOkB approximation is recovered from the PNOkB one whenever the reference state is particle-number conserving, i.e. whenever $|\Phi\rangle$ reduces to a Slater determinant $|\text{SD}\rangle$.

B.2 Examples

B.2.1 PNO2B approximation of a three-body operator

Let us exemplify the case of largest interest in current nuclear structure ab initio calculations, i.e. the PNO2B approximation of the three-body operator

[TABLE]

The normal-ordered form of $O$ with respect to $|\Phi\rangle$ reads in the single-particle basis as

[TABLE]

and involves the set of normal and anomalous fields

[TABLE]

The nNO2B approximation leads to just dropping $O^{[3]}=\Lambda^{33}$. The application of the PNO2B approximation is more involved and we now proceed to the construction of the corresponding operator

[TABLE]

where the different terms are to be obtained recursively on the basis of Eq. (2.4.2). In the present case, it leads to

[TABLE]

where the needed intermediate quantities, i.e. the diagonal normal-ordered fields originating from the successive contributions to $O^{\text{PNO2B}}$, are nothing but

[TABLE]

The two above equations fully define the PNO2B approximation of the three-body operator. It is particle-number conserving by construction. Recent ab initio BMBPT calculations Tichai:2018mll have been performed on the basis of this approximation although this was not explicited at the time.

It is useful to further characterize the operator by closely inspecting its normal-ordered contributions. Focusing for example on the normal field $\tilde{\Lambda}^{11}$, one obtains

[TABLE]

which indeed satisfies the systematic property $\tilde{\Lambda}^{ii}=\Lambda^{ii}$. Focusing now on the anomalous field $\tilde{\Lambda}^{20}$, one obtains

[TABLE]

where again no extra term arises.

Eventually, the complete set of normal-ordered contributions to $O^{\text{PNO2B}}$ expressed in the single-particle basis relates to those of the original three-body operator through

[TABLE]

One thus observes that the non-zero fields are strictly equal to those associated with the original operator, i.e. there is no so-called extra term when approximating a three-body operator.

A graphical representation of the PNO2B approximation of a three-body operator $O$ is given in Fig. 9.

B.2.2 PNO3B approximation of a four-body operator

Let us now build the PNO3B approximation of the four-body operator

[TABLE]

The normal ordering of $O$ with respect to $|\Phi\rangle$ reads in the single-particle basis as

[TABLE]

and involves the following normal and anomalous fields

[TABLE]

The nNO3B approximation leads to just dropping $O^{[4]}=\Lambda^{44}$. The application of the PNO3B approximation is more involved and we now proceed to the construction of the corresponding operator

[TABLE]

where the different terms are to be obtained recursively on the basis of Eq. (2.4.2). In the present case, it leads to

[TABLE]

where the needed intermediate quantities, i.e. the diagonal normal-ordered fields originating from the successive contributions to $O^{\text{PNO3B}}$, are

[TABLE]

The two above equations fully define the PNO3B approximation of the four-body operator. It is particle-number conserving by construction. It is useful to further characterize the operator by closely inspecting its normal-ordered contributions. Focusing for example on the normal field $\tilde{\Lambda}^{22}$, one obtains

[TABLE]

which indeed satisfies the systematic property $\tilde{\Lambda}^{ii}=\Lambda^{ii}$. Focusing now on the anomalous field $\tilde{\Lambda}^{20}$, one obtains

[TABLE]

where the extra term is nothing but $\breve{\Lambda}^{20}$. Four-body operators are thus the first for which such an extra term appears.

In order to check the consistency of the derivation, the extra term can also be obtained via the application of Eq. (A.39) to the present case of interest. This gives

[TABLE]

which is indeed consistent with Eq. (A.55). Inspecting Eq. (A.55), one realizes that the difference between $\tilde{\Lambda}^{20}$ and $\Lambda^{20}$ eventually boils down to using a prefactor $1/8$ instead of $1/16$ in the term $\text{Tr}[o^{44}\kappa^{\ast}\kappa\kappa]$. This originates from the non-equivalent combinatorial prefactor at play when generating contributions of the type $\text{Tr}[o^{44}\kappa^{\ast}\kappa\kappa]$ with three anomalous contractions to $\tilde{\Lambda}^{20}$ and $\Lambda^{20}$ via the application of Wick’s theorem to $O^{\text{PNO3B}}$ and $O$, respectively. The difference between the normal-ordered fields of the PNOkB operator and of the original operator, if any, can always be rephrased in terms of modified prefactors of certain contributions involving strings of anomalous contractions.

Eventually, the complete set of normal-ordered contributions to $O^{\text{PNO3B}}$ expressed in the single-particle basis relates to those of the original four-body operator through

[TABLE]

Among the non-zero fields, one (plus its hermitian conjugate) contains an extra contribution, i.e. a term with a modified prefactor compared to the original operator $O$.

B.2.3 PNO2B approximation of a four-body operator

Let us now further approximate a four-body operator by building its PNO2B approximation. Starting from the operator defined in Eq. (A.48), the nNO2B approximation now leads to dropping $O^{[3]}+O^{[4]}=\Lambda^{42}+\Lambda^{33}+\Lambda^{24}+\Lambda^{44}$. The application of the PNO2B approximation is more involved and we now proceed to the construction of the corresponding operator

[TABLE]

where the different terms are to be obtained recursively on the basis of Eq. (2.4.2). In the present case, it leads to

[TABLE]

where the needed intermediate quantities, i.e. the diagonal normal-ordered fields originating from the successive contributions to $O^{\text{PNO2B}}$, are

[TABLE]

The two above equations fully define the PNO2B approximation of the four-body operator. It is particle-number conserving by construction. It is useful to further characterize the operator by closely inspecting its normal-ordered contributions. Focusing for example on the normal field $\tilde{\Lambda}^{11}$, one obtains

[TABLE]

which indeed satisfies the systematic property $\tilde{\Lambda}^{ii}=\Lambda^{ii}$. Focusing now on the anomalous field $\tilde{\Lambda}^{20}$, one obtains

[TABLE]

where the extra term is nothing but $\breve{\Lambda}^{20}$. One notices that it is the same as in the PNO3B approximation of the four-body operator.

Eventually, the complete set of normal-ordered contributions to $O^{\text{PNO2B}}$ expressed in the single-particle basis relates to those of the original four-body operator through

[TABLE]

Among the non-zero fields, one (plus its hermitian conjugate) contains an extra contribution, i.e. a term with a modified prefactor compared to the original operator $O$.

B.2.4 PNO3B and PNO4B approximations of a five-body operator

We have also worked out the PNO3B and the PNO4B approximations to a five-body operator in full details. As the expressions become lengthy, they are not reported here. Still, it is interesting to note that the particular form of the extra terms depends on which PNOkB, e.g. PNO3B or PNO4B, approximation is performed starting from the same original five-body operator. This point can be illustrated by only reporting how the normal-ordered contributions to both $O^{\text{PNO3B}}$ and $O^{\text{PNO4B}}$ relate to those of the original five-body operator. One has

[TABLE]

for the PNO3B approximation and

[TABLE]

for the PNO4B approximation. While in both cases the extra term entering $\tilde{\Lambda}^{31}$ is the same, one observes that $\tilde{\Lambda}^{20}$ does acquire an extra term in the PNO3B approximation but it does not in the PNO4B approximation.

For illustration, a graphical representation of the PNO3B approximation of the five-body operator $O$ is given in Fig. 10.

B.3 Alternative approximation

The procedure to define the PNOkB approximation is not unique. Given the three objectives stated in Sec. 2.4.2, two independent approaches cannot differ drastically in their philosophy and cannot produce very different approximate operators. In the present section, an alternative approximation procedure based on the quasi-normal ordering of Ref. kong10a is briefly investigated to highlight the similarities and differences.

B.3.1 Quasi normal-ordering

The quasi normal-ordering introduced in Ref. kong10a is stipulated through the set of equations

[TABLE]

where $\text{N}[\ldots]$ denotes the quasi normal-ordering defined by the given of the elementary contraction $\gamma$. The antisymmetrization operator $\mathcal{A}(\ldots)$ is the one defined in Ref. shamasundar09a . At this point, the contraction is defined as any two-index tensor such that the above procedure is very general. It happens that such a definition is sufficient to have a generalized Wick’s theorem governing the way the product of two quasi normal-ordered operators can be systematically decomposed into a sum of quasi normal-ordered operators weighted by a set of contractions kong10a .

In the present context, the quasi-normal ordering is only used as a systematic way to reshuffle as a large part as possible of an initial string of single-particle operators into operators of lower ranks in view of neglecting the quasi normal-ordered term having the same rank as the original string. Compared to the normal ordering with respect to the Bogoliubov vacuum associated with standard Wick’s theorem, it is presently possible to only invoke a symmetry-conserving elementary contraction such that each individual term in the quasi normal-ordered form is a scalar under $U(1)$ transformations. For the rest, the expectation value of a quasi-normal-ordered operator defined in such a general way does not vanish a priori, i.e.

[TABLE]

As a result, the quasi normal ordering does not present the practicalities of more traditional Wick’s theorems.

Thus, the approach proceeds in two steps. First, the original operator expressed in quasi normal-ordered form is truncated to produce the particle-number conserving quasi normal-ordered k-body (PQNOkB) approximation. Second, the resulting operator is brought into a normal-ordered form with respect to the Bogoliubov vacuum according to standard Wick’s theorem in view of using it in the many-body formalism of interest.

B.4 PQNO1B of a two-body operator

Applying quasi normal-ordering to a particle-number-conserving two-body operator $O$ leads to

[TABLE]

where the n-body fields $\Lambda^{nn(\text{qno})}$ are given by

[TABLE]

The PQNO1B approximation of $O$ is obtained by neglecting the effective two-body part $\Lambda^{22(\text{qno})}$

[TABLE]

where

[TABLE]

The above quantities are independent of the specific choice made for the elementary contraction $\gamma$. In the present case, $\gamma$ is taken as the normal density $\rho$ of the Bogoliubov vacuum. With this definition at hand, the PQNO1B operator is brought into a normal-ordered form with respect of the Bogoliubov state on the basis of standard Wick’s theorem

[TABLE]

where

[TABLE]

Comparing Eq. (A.73) to Eq. (34), one observes that the contribution to $\Lambda^{00}$ originating from the two-body operator via two anomalous contractions is not present in the PQNO1B operator. The latter thus corresponds to a more drastic appoximation of the original operator $O$ than the PNO1B operator. Performing the PQNO2B approximation of a three-body operator, the PQNO2B operator similarly misses a three-body contribution to $\Lambda^{11}$ associated with two anomalous contractions that is present in the PNO2B operator.

Appendix C Double particle-number projection

C.1 Rotated Bogoliubov vacuum

The rotated Bogoliubov state

[TABLE]

is a vacuum for the set of rotated quasiparticle operators defined through

[TABLE]

where the associated Bogoliubov transformation reads as

[TABLE]

The rotated vacuum can also be expressed via a non-unitary Thouless transformation of $|\Phi\rangle$ according to

[TABLE]

The Thouless operator appearing in Eq. (A.89)

[TABLE]

involves the Thouless matrix

[TABLE]

which is expressed in terms of the Bogoliubov transformation connecting the quasiparticle operators annihilating $|\Phi(\varphi)\rangle$ to those annihilating $|\Phi\rangle$

[TABLE]

where

[TABLE]

C.2 Gauge-rotated contractions

Matrix elements of the doubly gauge-rotated contractions are defined through

[TABLE]

such that singly-rotated ones are nothing but

[TABLE]

Doubly gauge-rotated contractions can be expressed in terms of the Bogoliubov transformation defining $|\Phi\rangle$ and of the gauge angles via RiSc80

[TABLE]

Setting $\varphi^{\prime}=0$ in Eq. (A.105), singly and doubly gauge-rotated contractions appear to be trivially related through

[TABLE]

C.3 Projection

With Eq. (A.106) at hand, the singly- (left-) projected mean-field matrix element on the value A of the one-body operator $F$ reads as

[TABLE]

while the doubly- (left- and right-) projected, on possibly two different values A and A’, is {strip}

[TABLE]

where $\mathcal{N}^{(0)}(\phi)$, $\rho(\phi)$, $\bar{\kappa}^{\ast}(\phi)$ and $\kappa(\phi)$ are $2\pi$-periodic functions.

C.4 Particle-number variance

In Fig. 8, a small (non-zero) value was obtained from the nNO1B approximation of the particle-number variance operation the basis of the double particle-number projection. While a zero variance must be obtained for the exact particle-number variance operator, the fact that its nNO1B approximation systematically delivers a (almost) zero value is not immediately obvious. This feature, for which the double projection is essential, is now briefly analyzed.

C.4.1 Doubly-projected naive NO1B expectation value

The doubly-projected expectation value of the nNO1B approximation of an operator $O$ can be systematically written as

[TABLE]

where

[TABLE]

and where the explicit form of $\tilde{o}^{00}$ and $\tilde{o}^{11}$ are given in Eq. (37).

C.4.2 Double projection and nNO1B approximation

If $A^{2}$ is approximated at the nNO1B level, Eq. (A.109) specifies to

[TABLE]

where the expressions

[TABLE]

were used along with the identity

[TABLE]

originating from the unitarity of the Bogoliubov transformation. Having constrained the average particle number in the Bogoliubov vacuum to equate the targeted particle number A, i.e.

[TABLE]

the doubly-projected nNO1B approximation to the particle-number variance is eventually given by

[TABLE]

The approximate particle-number variance is thus given in this case by the trace over $\tilde{(a^{2})}^{11}$ and the difference between PHFB and HFB normal density matrices $\Delta\rho^{\text{proj}}$. The latter being expected to be small, the doubly-projected nNO1B particle-number variance is expected to be small as well.

Appendix D Effect of particle number projection

As seen in panel (b) of Fig. 4, the particle-number projection provides further/no/lesser binding in 14,22O/16,24O/18,20,26O, knowing that the effect is essentially negligeable in 22,26O.

In order to analyze this trend, the HFB state and energy can be decomposed into their particle-number projected components according to the sum rules777These sum rules are obtained from the Fourier expansion of uncorrelated singly-rotated off-diagonal kernels $\mathcal{N}^{(0)}(\varphi)=\langle\Phi(\varphi)|\Phi\rangle$ and $O^{(0)}(\varphi)\equiv o^{(0)}(\varphi)\mathcal{N}^{(0)}(\varphi)=\langle\Phi(\varphi)|O|\Phi\rangle$ computed at $\varphi=0$.

[TABLE]

where the projected weights and energies are defined as

[TABLE]

The decomposition underlying Eq. (A.115a) is displayed in Fig. 11 along the sequence of oxygen isotopes. One first observes that 16O and 24O display a single component corresponding to their number of particles, i.e. the HFB vacuum reduces to a Slater determinant in these doubly closed-shell nuclei such that the subsequent particle number projection does not provide any static correlations. While in 14O the distribution is slightly skewed towards smaller particle numbers than the average value $14$, it is the opposite in 18,20O. As seen in Fig. 12, the skewness of the distribution eventually impacts the energy associated with each component obtained under the constraint that the energy sum rule (Eq. (A.115b)) must be fulfilled. In particular, the projected energy associated with the largest component, i.e. with the targeted particle number, ends up being more/less negative than the average, i.e. the HFB energy, in 14O/18,20O.