Multi-reference many-body perturbation theory for nuclei

Abstract

Perturbative and non-perturbative expansion methods already constitute a tool of choice to perform ab initio calculations over a significant part of the nuclear chart. In this context, the categories of accessible nuclei directly reflect the class of unperturbed state employed in the formulation of the expansion. The present work generalizes to the nuclear many-body context the versatile method of Burton and Thom (J Chem Theory Comput 16(4):5586, 2020) by formulating a perturbative expansion on top of a multi-reference unperturbed state mixing deformed non-orthogonal Bogoliubov vacua, i.e. a state obtained from the projected generator coordinate method (PGCM). Particular attention is paid to the part of the mixing taking care of the symmetry restoration, showing that it can be exactly contracted throughout the expansion, thus reducing significantly the dimensionality of the linear problem to be solved to extract perturbative corrections. While the novel expansion method, coined as PGCM-PT, reduces to the PGCM at lowest order, it reduces to single-reference perturbation theories in appropriate limits. Based on a PGCM unperturbed state capturing (strong) static correlations in a versatile and efficient fashion, PGCM-PT is indistinctly applicable to doubly closed-shell, singly open-shell and doubly open-shell nuclei. The remaining (weak) dynamical correlations are brought consistently through perturbative corrections. This symmetry-conserving multi-reference perturbation theory is state-specific and applies to both ground and excited PGCM unperturbed states, thus correcting each state belonging to the low-lying spectrum of the system under study. The present paper is the first in a series of three and discusses the PGCM-PT formalism in detail. The second paper displays numerical zeroth-order results, i.e. the outcome of PGCM calculations. Second-order, i.e. PGCM-PT(2), calculations performed in both closed- and open-shell nuclei are the object of the third paper.

Appendices

Appendix A: Permutation operators

Many of the algebraic expressions derived below can be economically written via the use of so-called permutation operators that perform appropriate anti-symmetrizations of the matrix element they act on. A permutation operator $P(s_1/s_2/\ldots /s_n)$, where $s_i$ ($i=\{1,\ldots ,n\}$) denotes a given set of indices, permutes the indices belonging to the different sets in all possible ways, without permuting the indices within each set. Furthermore, the sign given by the signature of each permutation multiplies the corresponding term. In the present work, the needed permutation operators read as

$$\begin{aligned} P(k_1/k_2)&\equiv {} 1 - P_{k_1k_2} ,\end{aligned}$$

(62a)

$$\begin{aligned} P(k_1/k_2k_3)&\equiv {} 1 - P_{k_1k_2} - P_{k_1 k_3},\end{aligned}$$

(62b)

$$\begin{aligned} P(k_1/k_2k_3k_4)&\equiv {} 1 - P_{k_1k_2} - P_{k_1k_3} - P_{k_1k_4} ,\end{aligned}$$

(62c)

$$\begin{aligned} P(k_1k_2/k_3k_4)&\equiv {} 1 - P_{k_1k_3} - P_{k_1k_4} - P_{k_2k_3} - P_{k_2k_4} + P_{k_1k_3}P_{k_2k_4} ,\end{aligned}$$

(62d)

$$\begin{aligned} P(k_1/k_2/k_3k_4)&\equiv {} P(k_1k_2/k_3k_4) P(k_1/k_2) \nonumber \\&= 1 - P_{k_1k_3} - P_{k_1k_4} - P_{k_2k_3} - P_{k_2k_4} + P_{k_1k_3}P_{k_2k_4}+ P_{k_3k_4} + P_{k_1k_3} P_{k_3k_4}\nonumber \\&\quad + P_{k_1k_4}P_{k_3k_4} + P_{k_2k_3}P_{k_3k_4} +P_{k_2k_4}P_{k_3k_4} - P_{k_1k_4}P_{k_2k_3} ,\end{aligned}$$

(62e)

$$\begin{aligned} P(k_1/k_2/k_3/k_4)&\equiv P(k_1k_2/k_3k_4) P(k_1/k_2)P(k_3/k_4) , \end{aligned}$$

(62f)

where the exchange operator $P_{k_ik_j}$ commutes indices $k_i$ and $k_j$.

Appendix B: Symmetry group

The symmetry group of H underlines the symmetry quantum numbers carried by its many-body eigenstates. In the present context, the group²⁸

$$\begin{aligned} \text {G}_{H} \equiv \text {SU(2)} \times \text {I} \times \text {U(1)}_N \times \text {U(1)}_Z \end{aligned}$$

(63)

associated with the conservation of total angular momentum, parity and neutron/proton numbers is explicitly considered. The group is a compact Lie group but is non Abelian as a result of SU(2).

1.1 B.1 Unitary representation

Each subgroup is represented on Fock space ${\mathcal {F}}$ via the set of unitary rotation operators

$$\begin{aligned} R_{\mathbf {J}}(\varOmega )&\equiv e^{-\imath \alpha J_z} e^{-\imath \beta J_y} e^{-\imath \gamma J_z} , \end{aligned}$$

(64a)

$$\begin{aligned} R_{N}(\varphi _{n})&\equiv e^{-i\varphi _{n} N} , \end{aligned}$$

(64b)

$$\begin{aligned} R_{Z}(\varphi _{p})&\equiv e^{-i\varphi _{p} Z} , \end{aligned}$$

(64c)

$$\begin{aligned} \varPi (\varphi _\pi )&\equiv e^{-i\varphi _{\pi } F} , \end{aligned}$$

(64d)

where $\varOmega \equiv (\alpha ,\beta ,\gamma )$, $\varphi _\pi$ and $\varphi _{n}$ ($\varphi _{p}$) denote Euler, parity and neutron- (proton-)gauge angles, respectively. The one-body operators entering the unitary representations of interest denote the generators of the group made out of the three components of the total angular momentum $\mathbf {J} = (J_x,J_y,J_z)$, neutron- (proton-)number N (Z) operators as well as of the one-body operator

$$\begin{aligned} F&\equiv \sum _{ab} f_{ab} c^{\dagger }_a c_b \end{aligned}$$

(65)

defined through its matrix elements [52]

$$\begin{aligned} f_{ab} \equiv \frac{1}{2}\left( 1- \pi _a \right) \delta _{ab} , \end{aligned}$$

(66)

where $\pi _a$ denotes the parity of one-body basis states that are presently assumed to carry a good parity. The eigenstates of H are characterized by

$$\begin{aligned}&J^2 |\varPsi ^{\sigma }_{\mu } \rangle \equiv \hbar ^2 J(J+1) |\varPsi ^{\sigma }_{\mu } \rangle , \end{aligned}$$

(67a)

$$\begin{aligned}&J_z |\varPsi ^{\sigma }_{\mu } \rangle \equiv \hbar M |\varPsi ^{\sigma }_{\mu } \rangle , \end{aligned}$$

(67b)

$$\begin{aligned}&N |\varPsi ^{\sigma }_{\mu } \rangle \equiv \text {N} |\varPsi ^{\sigma }_{\mu } \rangle , \end{aligned}$$

(67c)

$$\begin{aligned}&Z |\varPsi ^{\sigma }_{\mu } \rangle \equiv \text {Z} |\varPsi ^{\sigma }_{\mu } \rangle , \end{aligned}$$

(67d)

$$\begin{aligned}&\varPi (\pi ) |\varPsi ^{\sigma }_{\mu } \rangle \equiv \varPi |\varPsi ^{\sigma }_{\mu } \rangle , \end{aligned}$$

(67e)

where $\sigma \equiv (\text {J}\text {M}\varPi \text {NZ})$ and where the operator $\varPi (\pi )$ is nothing but the parity operator and $J^2 \equiv \mathbf {J} \cdot \mathbf {J}$ is the Casimir of SU(2).

The irreducible representations (IRREPs) of the group are given by [53]

$$\begin{aligned} \langle \varPsi ^{\sigma }_{\mu } | R(\theta ) |\varPsi ^{\sigma '}_{\mu '} \rangle \equiv D_{\text {MM'}}^{\tilde{\sigma }}(\theta ) \delta _{\tilde{\sigma }\tilde{\sigma }'} \delta _{\mu \mu '} , \end{aligned}$$

(68)

with $\tilde{\sigma } \equiv (\text {J}\varPi \text {NZ})$ and

$$\begin{aligned} D_{\text {MM'}}^{\tilde{\sigma }}(\theta )\equiv D_{\text {MM'}}^{\text {J}}(\varOmega ) e^{-i\varphi _{n} \text {N}}e^{-i\varphi _{p} \text {Z}} e^{-\frac{i}{2}(1-\varPi )\varphi _\pi } , \end{aligned}$$

(69)

and where the rotation operators have been gathered into

$$\begin{aligned} R(\theta ) \equiv R_{\mathbf {J}}(\varOmega )R_{N}(\varphi _{n})R_{Z}(\varphi _{p})\varPi (\varphi _\pi ) , \end{aligned}$$

(70)

with

$$\begin{aligned} \theta \equiv (\varOmega , \varphi _n,\varphi _p, \varphi _\pi ) \end{aligned}$$

(71)

encompassing all rotation angles. The domain of definition of the group is thus

$$\begin{aligned} D_{\text {G}_H}&\! \equiv \!D_{\alpha } \times D_{\beta } \times D_{\gamma } \times D_{\varphi _{n}}\times D_{\varphi _{p}} \times D_{\varphi _{\pi }} \nonumber \\&\!=\![0,4\pi ]\!\times \![0,\pi ]\!\times \![0,2\pi ]\!\times \![0,2\pi ]\!\times \![0,2\pi ]\!\times \!\{0,\pi \}. \end{aligned}$$

(72)

In Eq. (69), $D^{\text {J}}_{\text {MM'}}(\varOmega )$ denotes Wigner D-matrices that can be expressed in terms of (real) reduced Wigner d-functions through $D^{\text {J}}_{\text {MM'}}(\varOmega )\equiv e^{-i\text {M}\alpha } \, d^{\text {J}}_{\text {MM'}}(\beta ) \, e^{-i\text {M'}\gamma }$.

Given that the degeneracy of the IRREPs is $d_{\tilde{\sigma }} = 2J+1$ and the volume of the group is

(73)

the orthogonality of the IRREPs read as

$$\begin{aligned} \int _{D_{\text {G}}} \text {d}\theta D_{\text {MK}}^{\tilde{\sigma }*}(\theta ) D_{\text {M'K'}}^{\tilde{\sigma }'}(\theta )= \frac{v_{\text {G}_H}}{d_{\tilde{\sigma }}} \delta _{\tilde{\sigma }\tilde{\sigma }'} \delta _{MM'} \delta _{KK'} . \end{aligned}$$

(74)

Furthermore, the unitarity of the symmetry transformations, i.e. $R^{\dagger }(\theta )R(\theta )=R(\theta )R^{\dagger }(\theta )=1$, induces

$$\begin{aligned} \sum _{M} D_{\text {MK}}^{\tilde{\sigma }*}(\theta ) D_{\text {MK'}}^{\tilde{\sigma }}(\theta )&= \delta _{KK'} , \end{aligned}$$

(75a)

$$\begin{aligned} \sum _{K} D_{\text {MK}}^{\tilde{\sigma }}(\theta )D_{\text {M'K}}^{\tilde{\sigma }*}(\theta )&= \delta _{MM'} . \end{aligned}$$

(75b)

An irreducible tensor operator $T^{\tilde{\sigma }}_{K}$ of rank J and a state $| \varPsi ^{\tilde{\sigma }K}_{\mu } \rangle$ transform under rotation according to

$$\begin{aligned} R(\theta ) \, T^{\tilde{\sigma }}_{K} \, R(\theta )^{-1}= & {} \sum _{M} T^{\tilde{\sigma }}_{M} \, D^{\tilde{\sigma }}_{MK}(\theta ), \end{aligned}$$

(76a)

$$\begin{aligned} R(\theta ) \, | \varPsi ^{\tilde{\sigma } K}_{\mu } \rangle= & {} \sum _{M} | \varPsi ^{\tilde{\sigma } M}_{\mu } \rangle \, D^{\tilde{\sigma }}_{MK}(\theta ). \end{aligned}$$

(76b)

Peter–Weyl’s theorem ensures that any function $f(\theta ) \in L^2(\text {G}_H)$ can be expanded according to

$$\begin{aligned} f(\theta ) \equiv \sum _{\tilde{\sigma }} \sum _{MK} \, f^{\tilde{\sigma }}_{MK} \, \, D^{\tilde{\sigma }*}_{MK}(\theta ) , \end{aligned}$$

(77)

such that the set of complex expansion coefficients $\{f^{\tilde{\sigma }}_{MK}\}$ can be extracted thanks to the orthogonality of the IRREPs through

$$\begin{aligned} f^{\tilde{\sigma }}_{MK} = \frac{d_{\tilde{\sigma }}}{v_{\text {G}_H}} \int _{D_{\text {G}_H}} \text {d}\theta D_{\text {MK}}^{\tilde{\sigma }}(\theta ) f(\theta ) . \end{aligned}$$

(78)

1.2 B.2 Projection operators

The operator

$$\begin{aligned} P^{\sigma } \equiv P^{\text {J}}_{\text {M}} P^{\text {N}} P^{\text {Z}} P^{\varPi } \end{aligned}$$

(79)

collects the projection operators on good symmetry quantum numbers

$$\begin{aligned} P^{\text {J}}_{\text {M}}&\equiv \sum _K g_K P^{\text {J}}_{\text {MK}} \nonumber \\&\equiv \sum _K g_K \frac{2J+1}{16\pi ^2}\int _{[0,4\pi ]\times [0,\pi ]\times [0,2\pi ]} \mathrm {d}\varOmega \, D^{\text {J}*}_{\text {MK}}(\varOmega ) R_{\mathbf {J}}(\varOmega ) , \end{aligned}$$

(80a)

$$\begin{aligned} P^{\text {N}}&\equiv \frac{1}{2\pi }\int _{0}^{2\pi } \mathrm {d}\varphi _{n} e^{i\varphi _{n} \text {N}} R_{N}(\varphi _{n}) , \end{aligned}$$

(80b)

$$\begin{aligned} P^{\text {Z}}&\equiv \frac{1}{2\pi }\int _{0}^{2\pi } \mathrm {d}\varphi _{p} e^{i\varphi _{p} \text {Z}} R_{Z}(\varphi _{p}) , \end{aligned}$$

(80c)

$$\begin{aligned} P^{\varPi }&\equiv \frac{1}{2} \sum _{\varphi _\pi =0,\pi } e^{\frac{i}{2}(1-\varPi )\varphi _\pi } \varPi (\varphi _\pi ), \end{aligned}$$

(80d)

such that one can write in a compact way

$$\begin{aligned} P^{\sigma }&=\, \frac{d_{\tilde{\sigma }}}{v_{\text {G}}} \sum _K g_K \int _{D_{\text {G}}} \text {d}\theta D_{\text {MK}}^{\tilde{\sigma }*}(\theta ) R(\theta ) \nonumber \\&\equiv \sum _K g_K P^{\tilde{\sigma }}_{\text {MK}} . \end{aligned}$$

(81)

The so-called transfer operator $P^{\tilde{\sigma }}_{\text {MK}}$ fulfills

$$\begin{aligned}&P^{\tilde{\sigma }}_{\text {MK}} = \sum _{\mu } | \varPsi ^{\tilde{\sigma }M}_{\mu } \rangle \langle \varPsi ^{\tilde{\sigma }K}_{\mu } | , \end{aligned}$$

(82a)

$$\begin{aligned}&P^{\tilde{\sigma }\dagger }_{\text {MK}} = P^{\tilde{\sigma }}_{\text {KM}} , \end{aligned}$$

(82b)

$$\begin{aligned}&P^{\tilde{\sigma }}_{\text {MK}} P^{\tilde{\sigma }'}_{\text {M'K'}} = \delta _{\tilde{\sigma }\tilde{\sigma }'}\delta _{KM'} P^{\tilde{\sigma }}_{\text {MK'}} , \end{aligned}$$

(82c)

along with the identity

$$\begin{aligned} P^{\sigma } R(\theta ) = \sum _K g_K \sum _{M'} D_{\text {KM'}}^{\tilde{\sigma }}(\theta ) P^{\tilde{\sigma }}_{\text {MM'}} . \end{aligned}$$

(83)

The present paper is eventually interested in the particular case where $g_K = \delta _{K0}$.

Appendix C: Bogoliubov algebra

1.1 C.1 Operators definition

Given a basis ${\mathcal {B}}_1 \equiv \{| l \rangle \}$ of the one-body Hilbert space ${\mathcal {H}}_1$ whose associated set of particle creation and annihilation operators is denoted as $\{c_l^\dagger , c_l\}$, an arbitrary particle-number-conserving operator $O$ is represented as

$$\begin{aligned} O\equiv \sum _{n=0}^{r} {O^{nn}} , \end{aligned}$$

(84)

where each $n$-body component²⁹

$$\begin{aligned} O^{nn}\equiv \frac{1}{n!} \frac{1}{n!} \sum _{\begin{array}{c} a_1\cdots a_n\\ b_1\cdots b_n \end{array}} o^{a_1\cdots a_n}_{b_1\cdots b_n} \, C^{a_1\cdots a_n}_{b_1\cdots b_n} . \end{aligned}$$

(85)

Here

$$\begin{aligned} C^{a_1\cdots a_n}_{b_1\cdots b_n}\equiv c^\dag _{a_1}\cdots c^\dag _{a_n}c_{b_n}\cdots c_{b_1} \end{aligned}$$

(86)

defines a string of n particle creation and n particle annihilation operators such that

$$\begin{aligned} \left( C^{a_1\cdots a_n}_{b_1\cdots b_n}\right) ^{\dagger }=C_{a_1\cdots a_n}^{b_1\cdots b_n} . \end{aligned}$$

(87)

The string is in normal order with respect to the particle vacuum $|0\rangle$

$$\begin{aligned} N(C^{a_1\cdots a_n}_{b_1\cdots b_n})=C^{a_1\cdots a_n}_{b_1\cdots b_n} , \end{aligned}$$

(88)

where $N(\ldots )$ denotes the normal ordering with respect to $|0\rangle$, and it is anti-symmetric under the exchange of any pair of upper or lower indices, i.e.

$$\begin{aligned} C^{a_1\cdots a_n}_{b_1\cdots b_n} = \epsilon (\sigma _u) \epsilon (\sigma _l) \, C^{\sigma _u(a_1\cdots a_n)}_{\sigma _l(b_1\cdots b_n)} , \end{aligned}$$

(89)

where $\epsilon (\sigma _u)$ ($\epsilon (\sigma _l)$) refers to the signature of the permutation $\sigma _u(\ldots )$ ($\sigma _l(\ldots )$) of the n upper (lower) indices.

In Eq. (85), the n-body matrix elements $\{o^{a_1\cdots a_n}_{b_1\cdots b_n}\}$ constitute a mode-2n tensor, i.e. a data array carrying 2n indices, associated with the string they multiply. The n-body matrix elements are also fully anti-symmetric under the exchange of any pair of upper or lower indices, i.e.

$$\begin{aligned} o^{a_1\cdots a_n}_{b_1\cdots b_n} = \epsilon (\sigma _u) \epsilon (\sigma _l) \, o^{\sigma _u(a_1\cdots a_n)}_{\sigma _l(b_1\cdots b_n)} . \end{aligned}$$

(90)

1.2 C.2 Bogoliubov state

The linear Bogoliubov transformation ${\mathcal {W}}(q)$ connects the set of particle creation and annihilation operators to a set of quasi-particle creation and annihilation operators obeying fermionic anti-commutation rules

$$\begin{aligned} \{\beta _k(q),\beta _l(q)\}&= 0 , \end{aligned}$$

(91a)

$$\begin{aligned} \{\beta ^\dagger _k(q), \beta ^\dagger _l(q)\}&= 0 , \end{aligned}$$

(91b)

$$\begin{aligned} \{\beta _k(q), \beta ^\dagger _l(q)\}&= \delta _{kl} . \end{aligned}$$

(91c)

Formally the transformation reads [22]

$$\begin{aligned} \begin{pmatrix} \beta (q)\\ \beta ^\dagger (q) \end{pmatrix} \equiv {\mathcal {W}}^\dagger (q) \begin{pmatrix} c\\ c^\dagger \end{pmatrix} , \end{aligned}$$

(92)

with

$$\begin{aligned} {\mathcal {W}}(q)\equiv \begin{pmatrix} U(q)&{}\quad V^*(q)\\ V(q)&{}\quad U^*(q) \end{pmatrix} , \end{aligned}$$

(93)

so that the Bogoliubov transformation in expanded form reads as

$$\begin{aligned} \beta _k(q)\equiv & {} \sum _l {U^*}_{lk}(q) \, c_l + {V^*}_{k}^{l}(q) \, c^\dagger _l , \end{aligned}$$

(94a)

$$\begin{aligned} \beta _k^\dagger (q)\equiv & {} \sum _l U^{lk}(q) \, c^\dagger _l + V^{k}_l(q) \ c_l . \end{aligned}$$

(94b)

The anti-commutation rules (Eq. (91)) constrain ${\mathcal {W}}(q)$ to be unitary

$$\begin{aligned} {\mathcal {W}}^\dagger (q) {\mathcal {W}}(q) = {\mathcal {W}}(q) {\mathcal {W}}^\dagger (q) = 1 , \end{aligned}$$

(95)

which translates into

$$\begin{aligned} U^\dagger (q) U(q) + V^\dagger (q) V(q) =&1 , \end{aligned}$$

(96a)

$$\begin{aligned} V^T(q) U(q) + U^T(q) V(q) =&0 , \end{aligned}$$

(96b)

$$\begin{aligned} U(q) U^\dagger (q) + V^*(q) V^T(q) =&1 , \end{aligned}$$

(96c)

$$\begin{aligned} V(q) U^\dagger (q) + U^*(q) V^T(q) =&0 . \end{aligned}$$

(96d)

The normalized Bogoliubov product state $| \varPhi (q) \rangle$ is defined, up to a phase, as the vacuum of the quasi-particle operators, i.e.

$$\begin{aligned} \beta _k(q) | \varPhi (q) \rangle \equiv 0,\quad \forall k . \end{aligned}$$

(97)

Contrary to Slater determinants, which constitute a subset of Bogoliubov states, the latter are not eigenstates of neutron and proton number operators in general.

1.3 C.3 One-body density matrices

The Bogoliubov vacuum is fully characterized by its normal $\rho (q)$ and anomalous $\kappa (q)$ one-body density matrices whose matrix elements are defined through

$$\begin{aligned} \rho ^{l_1}_{l_2}(q)&\equiv \langle \varPhi (q) | c^\dagger _{l_2} c_{l_1} | \varPhi (q) \rangle = V^*(q) V^T(q) , \ \end{aligned}$$

(98a)

$$\begin{aligned} \kappa ^{l_1l_2}(q)&\equiv \langle \varPhi (q) | c_{l_2} c_{l_1} | \varPhi (q) \rangle = V^*(q) U^T(q) , \end{aligned}$$

(98b)

such that

$$\begin{aligned} \rho ^\dagger (q)&= \rho (q) , \end{aligned}$$

(99a)

$$\begin{aligned} \kappa ^T(q)&= - \kappa (q) . \end{aligned}$$

(99b)

1.4 C.4 Normal-ordered operators

It is eventually useful to normal order operators with respect to the Bogoliubov vacuum $| \varPhi (q) \rangle$ and express them in terms of the associated quasi-particle operators. Applying standard Wick’s theorem leads thus to the rewriting of O according to

$$\begin{aligned} O\equiv \sum _{n=0}^{r} \mathbf{O}^{[2n]}(q)\equiv \sum _{n=0}^{r} \sum _{\begin{array}{c} i,j=0/\\ i+j=2n \end{array}}^{2r} \mathbf{O}^{ij}(q) , \end{aligned}$$

(100)

with³⁰

$$\begin{aligned} \mathbf{O}^{ij}(q)\equiv \frac{1}{i!} \frac{1}{j!}\sum _{\begin{array}{c} k_1\cdots k_i\\ l_1\cdots l_j \end{array}} \mathbf{o}^{k_1\cdots k_i}_{l_1\cdots l_j}(q) \, B^{k_1\cdots k_i}_{l_1\cdots l_j}(q) , \end{aligned}$$

(101)

where

$$\begin{aligned} B^{k_1\cdots k_i}_{l_1\cdots l_j}(q)\equiv \beta ^\dagger _{k_1}(q)\cdots \beta ^\dagger _{k_i}(q) \beta _{l_j}(q)\cdots \beta _{l_1}(q) \end{aligned}$$

(102)

denotes a string of i quasi-particle creation and j quasi-particle annihilation operators such that

$$\begin{aligned} \left( B^{k_1\cdots k_i}_{l_1\cdots l_j}(q)\right) ^{\dagger }=B_{k_1\cdots k_i}^{l_1\cdots l_j}(q) . \end{aligned}$$

(103)

The string is in normal order with respect to the Bogoliubov state $| \varPhi (q) \rangle$

$$\begin{aligned} :B^{k_1\cdots k_i}_{l_1\cdots l_j}(q):=B^{k_1\cdots k_i}_{l_1\cdots l_j}(q) , \end{aligned}$$

(104)

where $:\ldots :$ denotes the normal ordering with respect to $| \varPhi (q) \rangle$, and it is anti-symmetric under the exchange of any pair of upper or lower indices, i.e.

$$\begin{aligned}&B^{k_1\cdots k_i}_{l_1\cdots l_j}(q) = \epsilon (\sigma _u) \epsilon (\sigma _l) \, B^{\sigma _u(k_1\cdots k_i)}_{\sigma _l(l_1\cdots l_j)}(q) . \end{aligned}$$

(105)

In Eq. (101), the matrix elements $\{\mathbf{o}^{k_1\cdots k_i}_{l_1\cdots l_j}(q)\}$ are fully anti-symmetric under the exchange of any pair of upper or lower indices, i.e.

$$\begin{aligned} \mathbf{o}^{k_1\cdots k_i}_{l_1\cdots l_j}(q) = \epsilon (\sigma _u) \epsilon (\sigma _l) \, \mathbf{o}^{\sigma _u(k_1\cdots k_i)}_{\sigma _l(l_1\cdots l_j)}(q) , \end{aligned}$$

(106)

and are functionals of the Bogoliubov matrices (U(q), V(q)) and of the matrix elements $\{o^{a_1\cdots a_n}_{b_1\cdots b_n}\}$ initially defining the operator O. As such, the content of each operator $\mathbf{O}^{ij}(q)$ depends on the rank r of O. For more details about the normal ordering procedure and for explicit expressions of the matrix elements up to $r=3$, see Refs. [18, 19, 24, 46, 47, 54].

1.5 C.5 Constrained Hartree–Fock–Bogoliubov theory

The state $|\varPhi (q)\rangle$ is obtained by minimizing its total energy under the constraints that it satisfies³¹

$$\begin{aligned} \langle \varPhi (q) | A | \varPhi (q) \rangle&= \text {A} , \end{aligned}$$

(107a)

$$\begin{aligned} \langle \varPhi (q) | Q | \varPhi (q) \rangle&= q , \end{aligned}$$

(107b)

where Q is a generic operator of interest. To do so, one considers the Routhian

$$\begin{aligned} R&\equiv H - \lambda _{\text {A}} (A-\text {A}) - \lambda _{q} (Q-q) \end{aligned}$$

(108a)

$$\begin{aligned}&\equiv \varOmega - \lambda _{q} (Q-q) , \end{aligned}$$

(108b)

where $\lambda _{\text {A}}$ and $\lambda _{q}$ denote two Lagrange parameters³². The Routhian reduces to the so-called grand potential $\varOmega$ whenever $\lambda _{q}=0$, i.e. when performing unconstrained calculations with respect to the order parameter q. Minimizing³³

$$\begin{aligned} R(q)&\equiv \langle \varPhi (q) | R | \varPhi (q) \rangle = \mathbf{R}^{00}(q) \end{aligned}$$

(109)

according to Ritz’ variational principle, the Bogoliubov matrices (U(q), V(q)) are found as solutions of the constrained Hartree–Fock–Bogoliubov (HFB) eigenequation [22]

$$\begin{aligned} \begin{pmatrix} {\overline{h}}(q) &{} \quad {\overline{\varDelta }}(q) \\ -{\overline{\varDelta }}^*(q) &{}\quad -{\overline{h}}^*(q) \end{pmatrix} \begin{pmatrix} U(q) \\ V(q) \end{pmatrix}_k&= E_k(q) \begin{pmatrix} U(q) \\ V(q) \end{pmatrix}_k , \end{aligned}$$

(110)

where the eigenvalues $\{E_k(q)\}$ are referred to as quasi-particle energies. The constrained HFB Hamiltonian matrix

$$\begin{aligned} {{\mathcal {H}}}(q)&\equiv \begin{pmatrix} {\overline{h}}(q) &{} {\overline{\varDelta }}(q) \\ -{\overline{\varDelta }}^*(q) &{} -{\overline{h}}^*(q) \end{pmatrix} , \end{aligned}$$

(111)

is built out of the constrained one-body Hartree-Fock and Bogoliubov fields

$$\begin{aligned} {\overline{h}}^{l}_{l'}(q)&\equiv \frac{\partial \mathbf{R}^{00}(q)}{\partial {\rho ^*}^{l}_{l'}(q)} \nonumber \\&= \langle \varPhi (q) |\{[c_l,R],c^{\dagger }_{l'}\}| \varPhi (q) \rangle \nonumber \\&= \frac{\partial \mathbf{H}^{00}(q)}{\partial {\rho ^*}^{l}_{l'}(q)}- \lambda _{q} \frac{\partial \mathbf{Q}^{00}(q)}{\partial {\rho ^*}^{l}_{l'}(q)} - \lambda _{\text {A}} \delta _{ll'} , \end{aligned}$$

(112a)

$$\begin{aligned} {\overline{\varDelta }}_{ll'}(q)&\equiv \frac{\partial \mathbf{R}^{00}(q)}{\partial \kappa ^*_{ll'}(q)} \nonumber \\&= \langle \varPhi (q) |\{[c_l,R],c_{l'}\}| \varPhi (q) \rangle \nonumber \\&= \frac{\partial \mathbf{H}^{00}(q)}{\partial \kappa ^*_{ll'}(q)} - \lambda _{q} \frac{\partial \mathbf{Q}^{00}(q)}{\partial \kappa ^*_{ll'}(q)} , \end{aligned}$$

(112b)

where

$$\begin{aligned} \frac{\partial \mathbf{H}^{00}(q)}{\partial {\rho ^*}^{l}_{l'}(q)}&= f^{l}_{l'}[| \varPhi (q) \rangle ] \end{aligned}$$

(113)

delivers nothing but the matrix elements of the one-body operator defined in Eq. (36) computed from the normal one-body density matrix of $| \varPhi (q) \rangle$.

At convergence, where the constraints are satisfied, the HFB energy is

$$\begin{aligned} \langle \varPhi (q) | H | \varPhi (q) \rangle&= \mathbf{H}^{00}(q) = \mathbf{R}^{00}(q) . \end{aligned}$$

(114)

Furthermore, Eq. (110) implies that

$$\begin{aligned} {\mathcal {W}}^\dagger (q) {{\mathcal {H}}}(q) {\mathcal {W}}(q)&= \begin{pmatrix} \mathbf{R}^{11}(q) &{} \mathbf{R}^{20}(q) \\ -\mathbf{R}^{20*}(q) &{} -\mathbf{R}^{11*}(q) \end{pmatrix} \nonumber \\&= \begin{pmatrix} E(q) &{} 0 \\ 0 &{} -E(q) \end{pmatrix}, \end{aligned}$$

(115)

such that the properties

$$\begin{aligned} \mathbf{R}^{20}(q)&= \mathbf{R}^{02}(q)=0 , \end{aligned}$$

(116a)

$$\begin{aligned} \mathbf{R}^{11}(q)&= \sum _{k} E_{k}(q) \beta ^{\dagger }_k(q) \beta _k(q) , \end{aligned}$$

(116b)

are fulfilled at convergence.

1.6 C.6 Elementary excitations

Given the Bogoliubov state $| \varPhi (q) \rangle$, a complete basis of Fock space ${{\mathcal {F}}}$ is obtained by generating all its elementary excitations

$$\begin{aligned} | \varPhi ^{k_1\cdots k_i}(q) \rangle \equiv B^{k_1\cdots k_i}(q) | \varPhi (q) \rangle , \end{aligned}$$

(117)

where $B^{k_1\cdots k_i}(q)$ defines the subclass of strings defined in Eq. (102) containing quasi-particle creation operators only.

It is interesting to note that each state defined through Eq. (117) is itself a Bogoliubov vacuum whose associated Bogoliubov transformation can be deduced from the one defining $| \varPhi (q) \rangle$ (Eq. (94)). Writing as $K_{n}\equiv \{k_1\cdots k_n\}$ the n-tuple defining a given elementary excitation $| \varPhi ^{K_n}(q) \rangle$, the associated Bogoliubov transformation $(U(q,K_n),V(q,K_n))$ is given by the matrices

$$\begin{aligned} U^{lk}(q,K_n)&\equiv U^{lk}(q)\quad \text {if}\ k \notin K_n , \end{aligned}$$

(118a)

$$\begin{aligned} U^{lk}(q,K_n)&\equiv V^{k*}_{l}(q)\quad \text {if}\ k \in K_n , \end{aligned}$$

(118b)

$$\begin{aligned} V^{l}_{k}(q,K_n)&\equiv V^{l}_{k}(q) \quad \text {if}\ k \notin K_n , \end{aligned}$$

(118c)

$$\begin{aligned} V^{l}_{k}(q,K_n)&\equiv U^*_{lk}(q) \quad \text {if}\ k \in K_n . \end{aligned}$$

(118d)

Such a consideration can be exploited to eventually compute matrix elements of operators between two Bogoliubov states that may differ not only by the value of the collective coordinate q but also by the elementary excitation character. The idea of evaluating matrix elements by redefining each elementary excitation of an original Bogoliubov vacuum as a novel vacuum is a generalization of the so-called generalized Slater–Condon rules [55]. The present work follows a numerically more efficient route where quasi-particle excitation operators are explicitly processed in order to limit the number of reference vacua to those entering the PGCM unperturbed state $|\mathrm \varTheta ^{\sigma }_{\mu } \rangle$ introduced in Eq. (28) and avoid the combinatorics associated with the redefinition of many Bogoliubov transformations through Eq. (118).

1.7 C.7 Rotated Bogoliubov state

Given the Bogoliubov state $| \varPhi (q) \rangle$, its rotated partner

$$\begin{aligned} |\varPhi (q;\theta )\rangle \equiv R(\theta ) |\varPhi (q) \rangle , \end{aligned}$$

(119)

is also a Bogoliubov state whose associated quasi-particle operators $\{\beta (q;\theta ), \beta ^\dagger (q;\theta )\}$ are characterized by the Bogoliubov transformation

$$\begin{aligned} {\mathcal {W}}(q;\theta )&= \begin{pmatrix} r(\theta )&{}\quad 0\\ 0&{}\quad r(\theta )^\dagger \end{pmatrix} {\mathcal {W}}(q) \nonumber \\&\equiv \begin{pmatrix} U(q;\theta )&{}\quad V^*(q;\theta )\\ V(q;\theta )&{}\quad U^*(q;\theta ) \end{pmatrix} , \end{aligned}$$

(120)

where $r(\theta )$ defines the matrix representation of $R(\theta )$ in the one-body Hilbert-space. Its elements are

$$\begin{aligned} r^{l_1}_{l_2}(\theta ) \equiv \langle l_1 | R(\theta ) | l_2 \rangle . \end{aligned}$$

(121)

Given that ${\mathcal {W}}(q;\theta )$ is a unitary Bogoliubov transformation, Eqs. (96) is also satisfied when substituting (U(q), V(q)) for $(U(q;\theta ),V(q;\theta ))$.

Because $R(\theta ) \in \text {G}_{H}$, the energy of the rotated HFB state

$$\begin{aligned} \langle \varPhi (q;\theta ) | H | \varPhi (q;\theta ) \rangle&= \mathbf{H}^{00}(q;\theta ) \, \end{aligned}$$

(122)

is in fact independent of the rotation angle; i.e. $\mathbf{H}^{00}(q;\theta )=\mathbf{H}^{00}(q)$ for all $\theta$.

Elementary excitations of the rotated Bogoliubov state are given by

$$\begin{aligned} | \varPhi ^{k_1\cdots k_i}(q;\theta ) \rangle \equiv B^{k_1\cdots k_i}(q;\theta ) | \varPhi (q; \theta ) \rangle \, \end{aligned}$$

(123)

where the rotated string reads as

$$\begin{aligned} B^{k_1\cdots k_i}(q;\theta )&\equiv R(\theta ) B^{k_1\cdots k_i}(q) R^{\dagger }(\theta ) \nonumber \\&= \beta ^\dagger _{k_1}(q;\theta )\cdots \beta ^\dagger _{k_i}(q;\theta ) , \end{aligned}$$

(124)

such that they are nothing but the rotated elementary excitations

$$\begin{aligned} | \varPhi ^{k_1\cdots k_i}(q;\theta ) \rangle \equiv R(\theta ) | \varPhi ^{k_1\cdots k_i}(q) \rangle \, . \end{aligned}$$

(125)

1.8 C.8 Single-reference partitioning

The presently developed perturbation theory is of MR character due to the fact that the PGCM unperturbed state (Eq. (28)) is a linear combination of several Bogoliubov vacua. However, in the limit where the PGCM state reduces to a single Bogoliubov state, which itself reduces to a Slater determinant whenever U(1) symmetry is conserved, PGCM-PT becomes of single-reference character and must thus entertain some connection with single-reference (B)MBPT [8]. To clarify this connection, the partitioning at play in the latter approaches are now briefly recalled.

1.8.1 C.8.1 BMBPT

Because of the inherent necessity to control the average particle number, the operator driving the perturbation in BMBPT is the grand potential $\varOmega$ [17]. Whenever $| \varPhi (q) \rangle$ results from a constrained HFB calculation (see Sect. C.5), a natural partitioning is given by

$$\begin{aligned} \varOmega = \varOmega _{0}(q) + \varOmega _{1}(q) \ , \end{aligned}$$

(126)

such that

$$\begin{aligned} \varOmega _{0}(q)&\equiv \varvec{\varOmega }^{00}(q) + \overline{\varvec{\varOmega }}^{11}(q) \ , \nonumber \\ \varOmega _{1}(q)&\equiv \varvec{\varOmega }^{20}(q) + \breve{\varvec{\varOmega }}^{11}(q) + \varvec{\varOmega }^{02}(q) \end{aligned}$$

(127a)

$$\begin{aligned}&\quad + \varvec{\varOmega }^{40}(q) \!+\! \varvec{\varOmega }^{31}(q) \!+\! \varvec{\varOmega }^{22}(q) \! +\! \varvec{\varOmega }^{13}(q) \!+ \!\varvec{\varOmega }^{04}(q), \end{aligned}$$

(127b)

with $\breve{\varvec{\varOmega }}^{11}(q)\equiv \varvec{\varOmega }^{11}(q) - \overline{\varvec{\varOmega }}^{11}(q)$ and where the diagonal one-body part of $\varOmega _{0}(q)$

$$\begin{aligned} \overline{\varvec{\varOmega }}^{11}(q) \equiv \sum _{k} E_k(q) \beta ^{\dagger }_k(q) \beta _k(q) , \end{aligned}$$

(128)

is built out of the positive eigenvalues generated through Eq. (110). In general, the partitioning defined in Eqs. (126)–(128) is not canonical. Indeed, while Eq. (116) is fulfilled for the Routhian R, it is not for $\varOmega$ except for $\lambda _q = 0$, i.e. whenever $| \varPhi (q) \rangle$ is the solution of an unconstrained HFB calculation.

The BMBPT expansion is formulated using the eigenbasis of $\varOmega _{0}(q)$ that is given by

$$\begin{aligned} \varOmega _{0}(q)\, | \varPhi (q) \rangle&= \varvec{\varOmega }^{00}(q) \, | \varPhi (q) \rangle , \end{aligned}$$

(129a)

$$\begin{aligned} \varOmega _{0}(q)\, | \varPhi ^{k_1 \ldots }(q) \rangle&= \left[ \varvec{\varOmega }^{00}(q) + E_{k_1}(q)+\ldots \right] | \varPhi ^{k_1 \ldots } \rangle . \end{aligned}$$

(129b)

1.8.2 C.8.2 MBPT

Whenever $q_{\text {U(1)}}=0$, $| \varPhi (q) \rangle$ is a Slater determinant and BMBPT reduces to MBPT. In this situation, the Bogoliubov field ${\overline{\varDelta }}(q)$ is zero and the Lagrange term associated with the particle number constraint entering the Routhian becomes superfluous and can be omitted. As a result, Eq. (110) reduces to

$$\begin{aligned} h(q) \, \left( U(q)\right) _k&= e_{k}(q) \, \left( U(q)\right) _k , \end{aligned}$$

(130)

i.e. to the constrained HF equation where the one-body HF field reads as

$$\begin{aligned} h^{l}_{l'}(q) \equiv f^{l}_{l'}[| \varPhi (q) \rangle ]- \lambda _{q} \frac{\partial \mathbf{Q}^{00}(q)}{\partial {\rho ^*}^{l}_{l'}(q)} . \end{aligned}$$

(131)

Solving Eq. (130) delivers constrained HF single-particle states $\{a^{\dagger }_k(q)\}$ through the unitary one-body basis transformation U(q) along with the associated HF single-particle energies

$$\begin{aligned} e_{k}(q) = f^{k}_{k}[| \varPhi (q) \rangle ]- \lambda _{q} \frac{\partial \mathbf{Q}^{00}(q)}{\partial {\rho ^*}^{k}_{k}(q)} . \end{aligned}$$

(132)

The HF Slater determinant is built by occupying the A lowest HF single-particle states

$$\begin{aligned} | \varPhi (q) \rangle&\equiv A^{i_1\cdots i_{\text {A}}}(q) | 0 \rangle , \end{aligned}$$

(133)

where a string $A^{p_1\cdots }_{h_1\cdots }(q)$ is defined in terms of constrained HF single-particle creation and annihilation operators.

Because the Lagrange term associated with the particle-number constraint is superfluous, the operator driving the perturbative expansion in MBPT is nothing but the Hamiltonian. Given the above, the unperturbed Hamiltonian deriving from Eq. (127) becomes

$$\begin{aligned} H_{0}(q)&\equiv \mathbf{H}^{00}(q) + :h(q): \nonumber \\&= E^{(0)}(q) + \sum _{k} e_{k}(q) :A^{k}_{k}(q):\ , \end{aligned}$$

(134)

where the latter form is given in the eigenbasis of h(q). Equation (132) makes clear that $h(q)=F_{[| \varPhi (q) \rangle ]}$ whenever $\lambda _{q}=0$ such that $\{e_{k}(q)\}$ denotes nothing but the eigenvalues of $F_{[| \varPhi (q) \rangle ]}$ in that particular case.

The A-body eigenbasis of $H_0$ is given by

$$\begin{aligned} H_0(q) | \varPhi (q) \rangle&= E^{(0)}(q) | \varPhi (q) \rangle , \end{aligned}$$

(135a)

$$\begin{aligned} H_0(q) | \varPhi ^{p_1\cdots }_{h_1\cdots }(q) \rangle&= E^{(0)}_{p_1\cdots h_1\cdots }(q) | \varPhi ^{p_1\cdots }_{h_1\cdots }(q) \rangle , \end{aligned}$$

(135b)

with

$$\begin{aligned} E^{(0)}(q)&\equiv \mathbf{H}^{00}(q) = \langle \varPhi (q) | H | \varPhi (q) \rangle , \end{aligned}$$

(136a)

$$\begin{aligned} E^{(0)}_{p_1\cdots h_1\cdots }(q)&\equiv E^{(0)}(q) + e_{p_1}(q) + \ldots - e_{h_1}(q) - \ldots , \end{aligned}$$

(136b)

where elementary particle-hole excitations of the unperturbed Slater determinant are defined through

$$\begin{aligned} | \varPhi ^{p_1\cdots }_{h_1\cdots }(q) \rangle&\equiv A^{p_1\cdots }_{h_1\cdots }(q) | \varPhi (q) \rangle . \end{aligned}$$

(137)

Whenever applied at the minimum of $\mathbf{H}^{00}(q)$, the spectrum of $H_0$ is typically non-degenerate with respect to elementary excitations³⁴, i.e. it displays a gap-full spectrum

$$\begin{aligned} E^{(0)}_{p_1\cdots h_1\cdots }(q) -E^{(0)}(q) > 0 . \end{aligned}$$

(138)

This is schematically illustrated in Fig. 3.

1.9 C.9 Overlap between Bogoliubov vacua

Given two Bogoliubov vacua $|\varPhi (q,\theta )\rangle$ and $|\varPhi (p)\rangle$, their overlap is a key ingredient to the calculation of the needed many-body matrix elements. To express the result, Bloch–Messiah–Zumino decompositions [22] of the Bogoliubov transformations ${{\mathcal {W}}}(p)$ and ${{\mathcal {W}}}(q)$ are invoked, e.g. the matrices defining ${{\mathcal {W}}}(p)$ are expressed as the product of unitary matrices D(p) and C(p) and special block-diagonal matrices ${{\overline{U}}}(p)$ and ${{\overline{V}}}(p)$ according to

$$\begin{aligned} U(p)&\equiv D(p) {\overline{U}}(p) C(p) , \end{aligned}$$

(139a)

$$\begin{aligned} V(p)&= D^*(p) {\overline{V}}(p) C(p) . \end{aligned}$$

(139b)

Further denoting by $v_{k}(p)$ the BCS-like coefficients making up ${\overline{V}}(p)$, the overlap eventually reads as [56]³⁵

$$\begin{aligned}&\langle \varPhi (p)|\varPhi (q;\theta )\rangle \nonumber \\&\quad ={} (-1)^n \frac{\det (C^*(p))\det (C(q))}{\prod _k^n v_{k}(p)v_{k}(q)} \nonumber \\&\qquad \times \mathrm {pf} \left[ \begin{pmatrix} V(p)^T U(p) &{}\quad V^T(p) \mathbf{r}^T(\theta ) V^*(q)\\ -V(q)^\dag \mathbf{r}(\theta ) V(p) &{}\quad U^\dag (q) V^*(q) \end{pmatrix} \right] , \end{aligned}$$

(140)

where 2n denotes the dimension of ${{\mathcal {H}}}_1$ and where the Pfaffian of a symplectic matrix has been considered.

1.10 C.10 Transition Bogoliubov transformation

Given the Bogoliubov vacua $|\varPhi (q,\theta )\rangle$ and $|\varPhi (p)\rangle$, the two sets of quasi-particle operators are related via the Bogoliubov transformation

$$\begin{aligned} \begin{pmatrix} \beta (q;\theta )\\ \beta ^\dagger (q;\theta ) \end{pmatrix}&= {\mathcal {W}}^\dag (q;\theta ) {\mathcal {W}}(p) \begin{pmatrix} \beta (p)\\ \beta ^\dag (p) \end{pmatrix} \nonumber \\&\equiv \begin{pmatrix} D^\dag (p,q;\theta ) &{}\quad E^\dag (p,q;\theta )\\ E^T(p,q;\theta ) &{}\quad D^T(p,q;\theta ) \end{pmatrix} \begin{pmatrix} \beta (p)\\ \beta ^\dag (p) \end{pmatrix} \nonumber \\&\equiv {\mathcal {W}}^\dag (p,q;\theta ) \begin{pmatrix} \beta (p)\\ \beta ^\dag (p) \end{pmatrix} , \end{aligned}$$

(141)

where

$$\begin{aligned} E(p,q;\theta )&\equiv V^T(p) U(q;\theta ) + U^T(p)V(q;\theta ) , \end{aligned}$$

(142a)

$$\begin{aligned} D(p,q;\theta )&\equiv U^\dagger (p) U(q;\theta ) + V^\dagger (p) V(q;\theta ) . \end{aligned}$$

(142b)

Given that ${\mathcal {W}}(p,q;\theta )$ is a unitary Bogoliubov transformation, Eq. (96) is also satisfied when substituting (U(q), V(q)) for $(D(p,q;\theta ),E(p,q;\theta ))$.

1.11 C.11 Similarity transformation

1.11.1 C.11.1 Thouless transformation

The two Bogoliubov vacua $|\varPhi (q,\theta )\rangle$ and $|\varPhi (p)\rangle$ can be connected via a non-unitary Thouless transformation

$$\begin{aligned} |\varPhi (q;\theta )\rangle = \langle \varPhi (p)|\varPhi (q;\theta )\rangle \exp \left[ \mathbf{Z}^{20}(p,q;\theta )\right] |\varPhi (p)\rangle ,\nonumber \\ \end{aligned}$$

(143)

where matrix elements of the Thouless operator

$$\begin{aligned} \mathbf{Z}^{20}(p,q;\theta ) \equiv \frac{1}{2}\sum _{k_1k_2} \mathbf{z}^{k_1k_2}(p,q;\theta ) B^{k_1k_2}(p) \end{aligned}$$

(144)

are expressed in terms of the transition Bogoliubov transformation between both vacua (Eqs. (141) and (142)) according to

$$\begin{aligned} \mathbf{z}(p,q;\theta ) = E^*(p,q;\theta ) D^{*-1}(p,q;\theta ) . \end{aligned}$$

(145)

1.11.2 C.11.2 Similarity-transformed operators

Given $|\varPhi (p)\rangle$, $|\varPhi (q,\theta )\rangle$ and an operator O, the similarity-transformed operator is introduced as

$$\begin{aligned} ^{Z}O \equiv e^{-\mathbf{Z}^{20}(p,q;\theta )} O e^{\mathbf{Z}^{20}(p,q;\theta )} , \end{aligned}$$

(146)

which obviously depends on $(p,q;\theta )$ via $\mathbf{Z}^{20}(p,q;\theta )$. Because the similarity transformation is not unitary, $^{Z}O$ is not hermitian. Such similarity-transformed operators appear repeatedly in the PGCM-PT formalism developed in the present work.

Normal ordering O with respect to $|\varPhi (p)\rangle$ according to Eqs. (100) and (101), $^{Z}O$ is obtained by simply replacing the quasi-particle operators $\{\beta _k^\dagger (p), \beta _k(p)\}$ by the similarity-transformed ones

$$\begin{aligned} \begin{pmatrix} ^Z\beta (p)\\ ^Z\beta ^{\dag }(p) \end{pmatrix}&\equiv e^{-\mathbf{Z}^{20}(p,q;\theta )}\begin{pmatrix} \beta (p)\\ \beta ^\dag (p) \end{pmatrix} e^{\mathbf{Z}^{20}(p,q;\theta )} \nonumber \\&= \begin{pmatrix} 1 &{}\quad \mathbf{z}(p,q;\theta )\\ 0 &{}\quad 1 \end{pmatrix} \begin{pmatrix} \beta (p)\\ \beta ^\dag (p) \end{pmatrix} \nonumber \\&\equiv {^{Z}{\mathcal {X}}}^\dagger (p,q;\theta ) \begin{pmatrix} \beta (p)\\ \beta ^\dag (p) \end{pmatrix} . \end{aligned}$$

(147)

Expressing the result in terms of the initial set $\{\beta _k^\dagger (p), \beta _k(p)\}$ and applying Wick’s theorem allows one to eventually express $^{Z}O$ in normal-ordered form with respect to $|\varPhi (p)\rangle$, i.e. according to Eqs. (100)–(102), where the set of $(p,q;\theta )$-dependent matrix elements are functions of the original set of matrix elements and of the matrix $\mathbf{z}(p,q;\theta )$. The explicit expressions of these matrix elements are provided in Appendix D for a two-body operator O, i.e. an operator with $r=2$ in Eqs. (84)–(86) and/or Eqs. (100)–(102).

As made clear in Appendix E, one also needs the similarity transformation of a de-excitation operator $B_{l_1 \ldots l_i}(p)$ acting on the corresponding vacuum bra $\langle \varPhi (p)|$, i.e.

$$\begin{aligned} \langle \varPhi (p)| \, {^{Z}B}_{l_1 \ldots l_i}(p)&= \langle \varPhi (p)| \prod _{n=i}^1 {^{Z}\beta }_{l_n}(p) \nonumber \\&= \langle \varPhi (p)| \prod _{n=i}^1 \!\left( \!\beta _{l_n}(p)+\sum _m \mathbf{z}^{l_il_m}\beta ^\dag _{l_m}(p)\!\right) , \end{aligned}$$

(148)

where the transformation (147) is used repeatedly and where the dependence of $\mathbf{z}(p,q;\theta )$ on $(p,q;\theta )$ is omitted for simplicity. This gives for a single de-excitation

$$\begin{aligned} \langle \varPhi (p)| \, {^{Z}B}_{l_1l_2}(p)&= \langle \varPhi ^{l_1l_2}(p)| + \mathbf{z}^{l_1l_2} \langle \varPhi (p)| , \end{aligned}$$

(149)

and for a double de-excitation

$$\begin{aligned} \langle \varPhi (p)| \, {^{Z}B}_{l_1l_2l_3l_4}(p) =&\langle \varPhi ^{l_1l_2l_3l_4}(p)| \nonumber \\&+ P(l_1l_2/l_3l_4) \, \mathbf{z}^{l_3l_4} \, \langle \varPhi ^{l_1l_2}(p)| \nonumber \\&+ P(l_1/l_3l_4) \, \mathbf{z}^{l_1l_2}{} \mathbf{z}^{l_3l_4} \, \langle \varPhi (p)|, \end{aligned}$$

(150)

where the final expressions are obtained by expanding the product of transformed quasi-particle operators, by applying Wick’s theorem and by acting on the bra to eliminate many null terms. The definition of the needed permutation operators can be found in Appendix A. Interestingly, one observes that the excitation rank is not increased through the similarity transformation in Eqs. (149) and (150).

1.11.3 C.11.3 Rotated/similarity-transformed operators

A rotated/similarity-transformed operator associated to an operator O and given $|\varPhi (p)\rangle$, $|\varPhi (q,\theta )\rangle$ is introduced as

$$\begin{aligned} ^{Z}O(\theta ) \equiv e^{-\mathbf{Z}^{20}(p,q;\theta )} R(\theta ) O R^{\dagger }(\theta ) e^{\mathbf{Z}^{20}(p,q;\theta )} , \end{aligned}$$

(151)

where the extra dependence in $\theta$ due to the additional rotation compared to $^{Z}O$ defined in Eq. (146) is made apparent in the newly introduced notation $^{Z}O(\theta )$ such that $^{Z}O(0) = {}^{Z}O$.

Of course, the particular form used to express the initial operator O does not impact the actual content of $^{Z}O$ or $^{Z}O(\theta )$. In the PGCM-PT formalism of present interest, it happens that $^{Z}O$ and $^{Z}O(\theta )$ arise for operators O that are initially normal ordered with respect to $| \varPhi (p) \rangle$ and $| \varPhi (q) \rangle$, respectively, and thus expressed in terms of quasi-particle operators $\{\beta _k^\dagger (p), \beta _k(p)\}$ and $\{\beta _k^\dagger (q), \beta _k(q)\}$, respectively. With this in mind, $^{Z}O(\theta )$ is obtained by simply replacing the quasi-particle operators $\{\beta _k^\dagger (q), \beta _k(q)\}$ by rotated/similarity-transformed ones

$$\begin{aligned} \begin{pmatrix} ^Z\beta (q;\theta )\\ ^Z\beta ^{\dag }(q;\theta ) \end{pmatrix}&\equiv e^{-\mathbf{Z}^{20}(p,q;\theta )} R(\theta ) \begin{pmatrix} \beta (q)\\ \beta ^\dag (q) \end{pmatrix} R^{\dagger }(\theta ) e^{\mathbf{Z}^{20}(p,q;\theta )} \nonumber \\&\equiv e^{-\mathbf{Z}^{20}(p,q;\theta )} \begin{pmatrix} \beta (q;\theta )\\ \beta ^\dag (q;\theta ) \end{pmatrix} e^{\mathbf{Z}^{20}(p,q;\theta )} . \end{aligned}$$

(152)

Eventually, the operator $^{Z}O(\theta )$ needs to be re-expressed in terms of the set $\{\beta _k^\dagger (p), \beta _k(p)\}$, the goal being to express all quantities involved in a many-body matrix element of interest in terms of a single set of quasi-particle operators. To do so, the rotated/similarity-transformed quasi-particle operators are written as

$$\begin{aligned} \begin{pmatrix} ^Z\beta (q;\theta )\\ ^Z\beta ^{\dag }(q;\theta ) \end{pmatrix}&\equiv e^{-\mathbf{Z}^{20}(p,q;\theta )} \begin{pmatrix} \beta (q;\theta )\\ \beta ^\dag (q;\theta ) \end{pmatrix} e^{\mathbf{Z}^{20}(p,q;\theta )} \nonumber \\&= {\mathcal {W}}^\dag (p,q;\theta ) {^{Z}{\mathcal {X}}}^\dag (p,q;\theta ) \begin{pmatrix} \beta (p)\\ \beta ^\dag (p) \end{pmatrix} \nonumber \\&\equiv {^{Z}{\mathcal {Y}}}^\dagger (p,q;\theta ) \begin{pmatrix} \beta (p)\\ \beta ^\dag (p) \end{pmatrix} , \end{aligned}$$

(153)

with

$$\begin{aligned}&{^{Z}{\mathcal {Y}}}^\dagger (p,q;\theta )\nonumber \\&\quad = \begin{pmatrix} D^\dag (p,q;\theta )&{}\quad D^\dag (p,q;\theta ) \mathbf{z}(p,q;\theta ) + E^\dag (p,q;\theta ) \\ E^T(p,q;\theta ) &{}\quad E^T(p,q;\theta ) \mathbf{z}(p,q;\theta ) + D^T(p,q;\theta ) \end{pmatrix} \nonumber \\&\quad = \begin{pmatrix} D^\dag (p,q;\theta )&{}\quad 0 \\ E^T(p,q;\theta ) &{}\quad D^{*-1}(p,q;\theta ) \end{pmatrix} \, \end{aligned}$$

(154)

where the second line is obtained by inserting Eq. (145) into the first one and utilizing Eqs. (96a) and (96b).

Here, as made clear in Appendix E, one only needs to perform the rotation/similarity-transformation of an excitation operator $B^{k_1 \ldots k_i}(q)$ acting on the vacuum $|\varPhi (p)\rangle$

$$\begin{aligned} {^{Z}B}^{k_1 \ldots k_j}(q;\theta ) \, |\varPhi (p)\rangle =\prod _{n=1}^j {^{Z}\beta }^{\dagger }_{k_n}(q;\theta ) |\varPhi (p)\rangle = \prod _{n=1}^j \sum _{k_m}\Big (E_{k_m}^{k_n}\beta _{k_m}(p) + {D^{-1\dag }}^{k_mk_n} \beta _{k_m}^\dag (p)\Big )|\varPhi (p)\rangle , \end{aligned}$$

(155)

where the transformation (153) and (154) has been used repeatedly and where the dependence of $E(p,q;\theta )$ and $D(p,q;\theta )$ on $(p,q;\theta )$ has been omitted for simplicity. This gives for a single excitation

$$\begin{aligned} {^{Z}B}^{k_1k_2}(q;\theta ) \, |\varPhi (p)\rangle = \sum _{j_1j_2} {D^{-1\dag }}^{j_1k_1} {D^{-1\dag }}^{j_2k_2}|\varPhi ^{j_1j_2}(p)\rangle + \sum _{j_1} E_{j_1}^{k_1}{D^{-1\dag }}^{j_1k_2}|\varPhi (p)\rangle , \end{aligned}$$

(156)

and for a double excitation

$$\begin{aligned} {^{Z}B}^{k_1k_2k_3k_4}(q;\theta ) |\varPhi (p)\rangle&= \sum _{j_1j_2j_3j_4} {D^{-1\dag }}^{j_1k_1} {D^{-1\dag }}^{j_2k_2} {D^{-1\dag }}^{j_3k_3} {D^{-1\dag }}^{j_4k_4} |\varPhi ^{j_1j_2j_3j_4}(p)\rangle \nonumber \\&\quad + P(k_1k_2/k_3k_4)\sum _{j_1j_2j_4} {D^{-1\dag }}^{j_1k_1} {D^{-1\dag }}^{j_2k_2} E_{j_4}^{k_3} {D^{-1\dag }}^{j_4k_4} |\varPhi ^{j_1j_2}(p)\rangle \nonumber \\&\quad + P(k_1/k_3k_4)\sum _{j_2j_4} E_{j_2}^{k_1} {D^{-1\dag }}^{j_2k_2} E_{j_4}^{k_3} {D^{-1\dag }}^{j_4k_4} |\varPhi (p)\rangle , \end{aligned}$$

(157)

where the final expressions are obtained by expanding the product of transformed quasi-particle operators, applying Wick theorem and acting on the ket to eliminate many vanishing terms. Interestingly, one observes that the excitation rank is not increased through the rotation and similarity transformation in Eqs. (156) and (157).

Appendix D: Similarity-transformed matrix elements

Let us consider a two-body operator O (see Eqs. (100)–(102) with $r=2$), in normal-ordered form with respect to $| \varPhi (p) \rangle$

$$\begin{aligned} O&\equiv \mathbf{O}^{00}(p) + \Big [\mathbf{O}^{11}(p) + \{\mathbf{O}^{20}(p) + \mathbf{O}^{02}(p)\}\Big ] + \Big [\mathbf{O}^{22}(p) + \{\mathbf{O}^{31}(p) + \mathbf{O}^{13}(p)\} + \{\mathbf{O}^{40}(p) + \mathbf{O}^{04}(p)\}\Big ] \nonumber \\&= \mathbf{O}^{00}(p)+ \frac{1}{(1!)^2}\sum _{k_1 k_2} \mathbf{o}^{k_1}_{k_2}(p) B^{k_1}_{k_2}(p) + \frac{1}{2!}\sum _{k_1 k_2} \mathbf{o}^{k_1 k_2}(p) B^{k_1k_2}(p) + \frac{1}{2!}\sum _{k_1 k_2} \mathbf{o}_{k_1 k_2}(p) B_{k_1k_2}(p)\\&\quad + \frac{1}{(2!)^{2}} \sum _{k_1 k_2 k_3 k_4} \mathbf{o}^{k_1 k_2}_{k_3 k_4}(p) B^{k_1 k_2}_{k_3 k_4}(p) + \frac{1}{3!1!}\sum _{k_1 k_2 k_3 k_4} \mathbf{o}^{k_1k_2k_3}_{k_4}(p) B^{k_1k_2k_3}_{k_4}(p) + \frac{1}{1!3!}\sum _{k_1 k_2 k_3 k_4} \mathbf{o}^{k_1}_{k_2 k_3 k_4}(p) B^{k_1}_{k_2 k_3 k_4}(p) \nonumber \\&\quad + \frac{1}{4!} \sum _{k_1 k_2 k_3 k_4} \mathbf{o}^{k_1 k_2 k_3 k_4}(p) B^{k_1 k_2 k_3 k_4}(p) + \frac{1}{4!} \sum _{k_1 k_2 k_3 k_4} \mathbf{o}_{k_1 k_2 k_3 k_4}(p) B_{k_1 k_2 k_3 k_4}(p) .\nonumber \end{aligned}$$

(158)

Expressing the similarity-transformed partner $^{Z}O$ (Eq. (146)) under the same form, its $(p,q;\theta )$-dependent matrix elements read as

$$\begin{aligned} ^{Z}{} \mathbf{O}^{00}&\equiv {^{Z}{} \mathbf{S}^{00}},\end{aligned}$$

(159a)

$$\begin{aligned} {^{Z}{} \mathbf{o}}_{k_1k_2}&\equiv {^{Z}s}_{k_1k_2},\end{aligned}$$

(159b)

$$\begin{aligned} {^{Z}{} \mathbf{o}}^{k_1}_{k_2}&\equiv {^{Z}{} \mathbf{s}}^{k_1}_{k_2} + \sum _{l_1} {^{Z}{} \mathbf{s}}_{l_1k_2} \mathbf{z}^{l_1k_1},\end{aligned}$$

(159c)

$$\begin{aligned} {^{Z}{} \mathbf{o}}^{k_1 k_2}&\equiv {^{Z}{} \mathbf{s}}^{k_1 k_2} + P(k_1/k_2) \sum _{l_1} {^{Z}{} \mathbf{s}}^{k_1}_{l_1} \mathbf{z}^{l_1k_2} - \sum _{l_1l_2} {^{Z}{} \mathbf{s}}_{l_1l_2} \mathbf{z}^{l_1k_1} \mathbf{z}^{l_2k_2},\end{aligned}$$

(159d)

$$\begin{aligned} {^{Z}{} \mathbf{o}}_{k_1k_2k_3k_4}&\equiv {^{Z}\mathbf{s}}_{k_1k_2k_3k_4},\end{aligned}$$

(159e)

$$\begin{aligned} {^{Z}{} \mathbf{o}}^{k_1}_{k_2k_3k_4}&\equiv {^{Z}\mathbf{s}}^{k_1}_{k_2k_3k_4} + \sum _{l_1} {^{Z}{} \mathbf{s}}_{l_1k_2k_3k_4} \mathbf{z}^{l_1k_1},\end{aligned}$$

(159f)

$$\begin{aligned} {^{Z}{} \mathbf{o}}^{k_1k_2}_{k_3k_4}&\equiv {^{Z}\mathbf{s}}^{k_1k_2}_{k_3k_4} +P(k_1/k_2) \sum _{l_2} {^{Z}\mathbf{s}}^{k_1}_{l_2k_3k_4} \mathbf{z}^{l_2k_2} - \sum _{l_1l_2}{^{Z}\mathbf{s}}_{l_1l_2k_3k_4} \mathbf{z}^{l_1k_1} \mathbf{z}^{l_2k_2},\end{aligned}$$

(159g)

$$\begin{aligned} {^{Z}{} \mathbf{o}}^{k_1k_2k_3}_{k_4}&\equiv {^{Z}\mathbf{s}}^{k_1k_2k_3}_{k_4} + P(k_3/k_1k_2)\sum _{l_3} {^{Z}\mathbf{s}}^{k_1k_2}_{l_3k_4} \mathbf{z}^{l_3k_3}\nonumber \\&\quad - P(k_1/k_2k_3) \sum _{l_2l_3}{^{Z}{} \mathbf{s}}^{k_1}_{l_2l_3k_4} \mathbf{z}^{l_2k_2}{} \mathbf{z}^{l_3k_3} - \sum _{l_1l_2l_3}{^{Z}\mathbf{s}}_{l_1l_2l_3k_4} \mathbf{z}^{l_1k_1} \mathbf{z}^{l_2k_2} \mathbf{z}^{l_3k_3} , \end{aligned}$$

(159h)

$$\begin{aligned} {^{Z}{} \mathbf{o}}^{k_1k_2k_3k_4}&\equiv {^{Z}\mathbf{s}}^{k_1k_2k_3k_4} +P(k_4/k_1k_2k_3) \sum _{l_4}{^{Z}\mathbf{s}}^{k_1k_2k_3}_{l_4} \mathbf{z}^{l_4k_4} - P(k_1k_2/k_3k_4) \sum _{l_3l_4}{^{Z}{} \mathbf{s}}^{k_1k_2}_{l_3l_4} \mathbf{z}^{l_3k_3} \mathbf{z}^{l_4k_4}\nonumber \\&\quad + P(k_1/k_2k_3k_4) \sum _{l_2l_3l_4}{^{Z}{} \mathbf{s}}^{k_1}_{l_2l_3l_4} \mathbf{z}^{l_2k_2}{} \mathbf{z}^{l_3k_3} \mathbf{z}^{l_4k_4} + \sum _{l_1l_2l_3l_4}{^{Z}{} \mathbf{s}}_{l_1l_2l_3l_4} \mathbf{z}^{l_1k_1} \mathbf{z}^{l_2k_2} \mathbf{z}^{l_3k_3} \mathbf{z}^{l_4k_4} . \end{aligned}$$

(159i)

The transformed matrix elements depend on $(p,q;\theta )$ through their dependence on Z and further depend on p through the matrix elements defining the original normal-ordered operator in Eq. (159a). All these dependencies have been dropped in Eqs. (159) and (160) for the sake of readability.

The intermediate matrix elements entering Eqs. (159) and (160) are defined through

$$\begin{aligned} ^{Z}{} \mathbf{S}^{00}&\equiv \mathbf{O}^{00} + \frac{1}{2} \sum _{l_1l_2} \mathbf{o}_{l_1l_2}{} \mathbf{z}^{l_1l_2} + \frac{1}{8} \sum _{l_1l_2l_3l_4} \mathbf{o}_{l_1l_2l_3l_4} \mathbf{z}^{l_1l_2}{} \mathbf{z}^{l_3l_4},\end{aligned}$$

(160a)

$$\begin{aligned} {^{Z}{} \mathbf{s}}_{l_1l_2}&\equiv \mathbf{o}_{l_1l_2} + \frac{1}{2} \sum _{l_3l_4} \mathbf{o}_{l_1l_2l_3l_4} \mathbf{z}^{l_3l_4},\end{aligned}$$

(160b)

$$\begin{aligned} {^{Z}{} \mathbf{s}}^{l_1}_{l_2}&\equiv \mathbf{o}^{l_1}_{l_2} + \frac{1}{2} \sum _{l_3l_4} \mathbf{o}^{l_1}_{l_2l_3l_4} \mathbf{z}^{l_3l_4}, \end{aligned}$$

(160c)

$$\begin{aligned} {^{Z}{} \mathbf{s}}^{l_1l_2}&\equiv \mathbf{o}^{l_1l_2} + \frac{1}{2} \sum _{l_3l_4} \mathbf{o}^{l_1l_2}_{l_3l_4} \mathbf{z}^{l_3l_4}, \end{aligned}$$

(160d)

$$\begin{aligned} {^{Z}{} \mathbf{s}}_{l_1l_2l_3l_4}&\equiv \mathbf{o}_{l_1l_2l_3l_4}, \end{aligned}$$

(160e)

$$\begin{aligned} {^{Z}{} \mathbf{s}}^{l_1}_{l_2l_3l_4}&\equiv \mathbf{o}^{l_1}_{l_2l_3l_4}, \end{aligned}$$

(160f)

$$\begin{aligned} {^{Z}{} \mathbf{s}}^{l_1l_2}_{l_3l_4}&\equiv \mathbf{o}^{l_1l_2}_{l_3l_4}, \end{aligned}$$

(160g)

$$\begin{aligned} {^{Z}{} \mathbf{s}}^{l_1l_2l_3}_{l_4}&\equiv \mathbf{o}^{l_1l_2l_3}_{l_4}, \end{aligned}$$

(160h)

$$\begin{aligned} {^{Z}{} \mathbf{s}}^{l_1l_2l_3l_4}&\equiv \mathbf{o}^{l_1l_2l_3l_4} , \end{aligned}$$

(160i)

and incorporate all terms where both indices of a given matrix $\mathbf{z}$ are contracted with indices of a matrix element of O.

Appendix E: PGCM-PT(2) matrix elements

From a technical viewpoint, the building blocks of the PGCM-PT formalism presented in Sect. 3 are many-body matrix elements of the following kind

$$\begin{aligned} O^{\tilde{\sigma }}_{pIqJ}&\equiv \langle \varPhi ^I(p)| O P^{\tilde{\sigma }}_{00} |\varPhi ^J(q)\rangle \nonumber \\&=\frac{d_{\tilde{\sigma }}}{v_{\text {G}}} \sum _\theta D_{\text {00}}^{\tilde{\sigma }*}(\theta ) \, \langle \varPhi (p)| B_{l_1\ldots l_i}(p) \,O \,R(\theta ) \,B^{k_1\ldots k_j}(q)| \varPhi (q) \rangle \nonumber \\&=\frac{d_{\tilde{\sigma }}}{v_{\text {G}}} \sum _\theta D_{\text {00}}^{\tilde{\sigma }*}(\theta ) \, \langle \varPhi (p)| \, {^{Z}B}_{l_1 \ldots l_i}(p) \,\, ^{Z}O {^{Z}B}^{k_1 \ldots k_j}(q;\theta ) \, | \varPhi (p) \rangle \, \langle \varPhi (p)|\varPhi (q;\theta )\rangle \nonumber \\&\equiv \frac{d_{\tilde{\sigma }}}{v_{\text {G}}} \sum _\theta D_{\text {00}}^{\tilde{\sigma }*}(\theta ) \, O_{pIqJ}(\theta ) \,\langle \varPhi (p)|\varPhi (q;\theta )\rangle , \end{aligned}$$

(161)

where I and J denote arbitrary i-tuple and j-tuple excitations of the corresponding vacua. Whenever the bra and/or ket is not excited, i.e. whenever the associated excitation operator is the identity, the index is conventionally put to 0. All quantities appearing in Eq. 161 have been introduced and/or worked out in the previous appendices.

While the matrix elements introduced in Eq. (161) are defined (and could be evaluated) for an operator and excitations of arbitrary ranks, the implementation of PGCM-PT(2) on the basis of a two-body Hamiltonian only requires a subset of them that are now worked out explicitly.

1.1 E.1 Type-1 matrix elements

The first category of matrix elements $O^{\tilde{\sigma }}_{pIqJ}$ is obtained whenever

• 1. O is a two-body operator,

• 2. $I=0,S,D$, where 0/S/D stands for no/single/double excitation,

• 3. $J=0$, i.e. the ket state is fixed to be the vacuum.

Starting thus from a two-body operator O in the form given by Eq. (159a), the many-body matrix elements of interest are worked out by exploiting Eqs. (149), (150), (159) and (160) and by applying Wick’s theorem with respect to $| \varPhi (p) \rangle$ such that

$$\begin{aligned} O_{p0q0}(\theta )&= {^Z\mathbf{O}}^{00} , \end{aligned}$$

(162a)

$$\begin{aligned} O_{pSq0}(\theta )&= \mathbf{z}^{k_1k_2} \, {^Z\mathbf{O}}^{00} + {^Z\mathbf{o}}^{k_1k_2} , \end{aligned}$$

(162b)

$$\begin{aligned} O_{pDq0}(\theta )&= P(k_1/k_3k_4)\mathbf{z}^{k_1k_2} \mathbf{z}^{k_3k_4} {^Z\mathbf{O}}^{00}\nonumber \\&\quad + P(k_1k_2/k_3k_4) \mathbf{z}^{k_3k_4} {^Z\mathbf{o}}^{k_1k_2} + {^Z\mathbf{o}}^{k_1k_2k_3k_4} , \end{aligned}$$

(162c)

where the $(p,q;\theta )$ dependencies of the various quantities appearing on the right-hand side have been omitted.

1.2 E.2 Type-2 matrix elements

The second category of matrix elements $O^{\tilde{\sigma }}_{pIqJ}$ is obtained whenever

• 1. O is a one-body operator,

• 2. $I=0,S,D$,

• 3. $J=0,S,D$.

Starting from a one-body operator O, i.e. a sub-part of the operator given by Eq. (159a), the evaluation of this second category of many-body matrix elements further requires the use of Eqs. (156) and (157) given that excitations of the ket are now in order.

Vacuum-to-vacuum and $(I=0,J=0)$ and excitation-to-vacuum $(I\ne 0,J=0)$ matrix elements can be deduced from Eq. (162) and are thus not repeated here. Vacuum-to-excitation $(I=0,J\ne 0)$ matrix elements are given by

$$\begin{aligned} O_{p0qS}(\theta )&= {^Z\mathbf{O}}^{00} \sum _{j_1} E_{j_1}^{k_1}{D^{-1\dag }}^{j_1k_2} + \sum _{j_1j_2} {^Z\mathbf{o}}_{j_1j_2} {D^{-1\dag }}^{j_1k_1} {D^{-1\dag }}^{j_2k_2} \end{aligned}$$

(163a)

$$\begin{aligned} O_{p0qD}(\theta )&= P(k_1/k_3k_4) {^Z\mathbf{O}}^{00} \sum _{j_2j_4} E_{j_2}^{k_1} {D^{-1\dag }}^{j_2k_2} E_{j_4}^{k_3} {D^{-1\dag }}^{j_4k_4} \nonumber \\&\quad + P(k_1k_2/k_3k_4) \sum _{j_1j_2j_4} {^Z\mathbf{o}}_{j_1j_2} {D^{-1\dag }}^{j_1k_1} {D^{-1\dag }}^{j_2k_2} E_{j_4}^{k_3} {D^{-1\dag }}^{j_4k_4} . \end{aligned}$$

(163b)

Excitation-to-excitation $(I\ne 0,J\ne 0)$ matrix elements are of course the most involved ones. Single-to-single $(I= S'=\{i_1 i_2\},J= S=\{k_1 k_2\})$ ones read as

$$\begin{aligned} O_{pS'qS}(\theta )&= P(i_1/i_2) {^Z\mathbf{O}}^{00} {D^{-1\dag }}^{i_1k_1}{D^{-1\dag }}^{i_2k_2} + P(i_1/i_2)P(k_1/k_2) \sum _{j_1} {^Z\mathbf{o}}^{i_1}_{j_1} {D^{-1\dag }}^{j_1k_1} {D^{-1\dag }}^{i_2k_2} \nonumber \\&\quad + P(k_1/k_2) {\mathbf{z}}^{i_1i_2} \sum _{j_1j_2} {D^{-1\dag }}^{j_1k_1} {^Z\mathbf{o}}_{j_1j_2} {D^{-1\dag }}^{j_2k_2} + P(i_1/i_2) {^Z\mathbf{o}}^{i_1i_2} \sum _{j_2} E_{j_2}^{k_1} {D^{-1\dag }}^{j_2k_2}\nonumber \\&\quad + {^Z\mathbf{O}}^{00} {\mathbf{z}}^{i_1i_2} \sum _{j_2} E_{j_2}^{k_1} {D^{-1\dag }}^{j_2k_2} . \end{aligned}$$

(163c)

Double-to-single $(I= D=\{i_1 i_2i_3 i_4\},J= S=\{k_1 k_2\})$ matrix elements read as

$$\begin{aligned} O_{pDqS}(\theta )= & {} P(i_1i_2 / i_3i_4) P(k_1/k_2) {D^{-1\dag }}^{i_1k_1} {D^{-1\dag }}^{i_2k_2} {^Z\mathbf{o}}^{i_3i_4} + P(i_1i_2 / i_3i_4) P(k_1 / k_2){\mathbf{z}}^{i_3i_4} {D^{-1\dag }}^{i_1k_1} {D^{-1\dag }}^{i_2k_2} {^Z\mathbf{O}}^{00} \nonumber \\&+ P(i_1/i_2 / i_3i_4) P( k_1 / k_2 ) {\mathbf{z}}^{i_3i_4} {D^{-1\dag }}^{i_1k_1}\sum _{j_2} {^Z\mathbf{o}}_{i_2j_2}{D^{-1\dag }}^{j_2k_2} + P(i_1i_2 / i_3i_4) {\mathbf{z}}^{i_1i_2} {^Z\mathbf{o}}^{i_3i_4} \sum _{j_2}E_{j_2}^{k_2} {D^{-1\dag }}^{j_2k_1} \nonumber \\&+ P(i_1i_2 / i_3i_4) {\mathbf{z}}^{i_1i_2}{\mathbf{z}}^{i_3i_4} \sum _{j_1j_2} {^Z\mathbf{o}}_{j_1j_2} {D^{-1\dag }}^{j_1k_1} {D^{-1\dag }}^{j_2k_2} + P(i_1/ i_3i_4) {\mathbf{z}}^{i_1i_2}{\mathbf{z}}^{i_3i_4} {^Z\mathbf{O}}^{00}\sum _{j_2}E_{j_2}^{k_2} {D^{-1\dag }}^{j_2k_1} . \end{aligned}$$

(164a)

Single-to-double $(I= S=\{i_1 i_2\},J= D=\{k_1 k_2k_3 k_4\})$ matrix elements read as

$$\begin{aligned} O_{pSqD}(\theta )&= P(k_1k_2 / k_3k_4) P(i_1 / i_2) {^Z\mathbf{o}}_{j_3j_4} {D^{-1\dag }}^{j_3k_3} {D^{-1\dag }}^{j_4k_4} {D^{-1\dag }}^{i_1k_1} {D^{-1\dag }}^{i_2k_2} \nonumber \\&\quad + P(k_1k_2 / k_3k_4)P(i_1 / i_2) {^Z\mathbf{O}}^{00} {D^{-1\dag }}^{i_1k_1} {D^{-1\dag }}^{i_2k_2} \sum _{j_3} E_{j_3}^{k_4}{D^{-1\dag }}^{j_3k_3}\nonumber \\&\quad + P(k_1 / k_2 / k_3k_4) P(i_1 / i_2) {D^{-1\dag }}^{i_1k_1} \left( \sum _{j_3} E_{j_3}^{k_3} {D^{-1\dag }}^{j_3k_4}\right) \left( \sum _{j_2} {D^{-1\dag }}^{j_2k_2} {^Z\mathbf{o}}^{i_2}_{j_2}\right) \nonumber \\&\quad + P(k_1k_2 / k_3k_4) {\mathbf{z}}^{i_2i_1} \left( \sum _{j_3} E_{j_3}^{k_4} {D^{-1\dag }}^{j_3k_3}\right) \left( \sum _{j_1j_2} {^Z\mathbf{o}}_{j_1j_2} {D^{-1\dag }}^{j_1k_1} {D^{-1\dag }}^{j_2k_2} \right) \nonumber \\&\quad + P(k1/ k_3k_4) {^Z\mathbf{o}}^{i_1i_2} \left( \sum _{j_1} E_{j_1}^{k_1} {D^{-1\dag }}^{j_1k_2}\right) \left( \sum _{j_3} E_{j_3}^{k_4} {D^{-1\dag }}^{j_3k_3}\right) \nonumber \\&\quad + P(k_1/ k_3k_4) {^Z\mathbf{O}}^{00} {\mathbf{z}}^{i_2i_1} \left( \sum _{j_3} E_{j_3}^{k_4} {D^{-1\dag }}^{j_3k_3}\right) \left( \sum _{j_1} E_{j_1}^{k_2} {D^{-1\dag }}^{j_1k_1}\right) . \end{aligned}$$

(164b)

Double-to-double $(I= D=\{i_1 i_2i_3 i_4\},J= D'=\{k_1 k_2k_3 k_4\})$ matrix elements read as

$$\begin{aligned} O_{pDqD'}(\theta )&= P(i_4/i_3/i_2/i_1) {^Z\mathbf{O}}^{00} {D^{-1\dag }}_{k_4i_4} {D^{-1\dag }}^{i_3k_3} {D^{-1\dag }}^{i_2k_2} {D^{-1\dag }}^{i_1k_1} \nonumber \\&\quad + P(k_4/k_3k_2k_1) P(i_4/i_3/i_2/i_1) \left( \sum _{j_4} {^Z\mathbf{o}}^{i_4}_{j_4} {D^{-1\dag }}^{j_4k_4} \right) {D^{-1\dag }}^{i_3k_3} {D^{-1\dag }}^{i_2k_2} {D^{-1\dag }}^{i_1k_1}\nonumber \\&\quad + P(i_1i_2/i_3i_4) {\mathbf{z}}^{i_1i_2} P(k_1k_2/k_3/k_4) {D^{-1\dag }}^{i_4k_4} {D^{-1\dag }}^{i_3k_3}\sum _{j_2j_1} {D^{-1\dag }}^{j_1k_1} {^Z\mathbf{o}}_{j_1j_2} {D^{-1\dag }}^{j_2k_2} \nonumber \\&\quad + P(k_1k_2/k_3k_4) P(i_1i_2/i_3/i_4) {^Z\mathbf{o}}^{i_1i_2} {D^{-1\dag }}^{i_4k_4} {D^{-1\dag }}^{i_3k_3} \left( \sum _{j_1} E_{j_1}^{k_2} {D^{-1\dag }}^{j_1k_1} \right) \nonumber \\&\quad + P(i_1i_2/i_3/i_4) P(k_1k_2/k_3k_4) {\mathbf{z}}^{i_1i_2} {^Z\mathbf{O}}^{00} {D^{-1\dag }}^{i_4k_4} {D^{-1\dag }}^{i_3k_3}\left( \sum _{j_1} E_{j_1}^{k_2} {D^{-1\dag }}^{j_1k_1} \right) \nonumber \\&\quad + P(i_1i_2/i_3/i_4) P(k_1k_2/k_3/k_4) {\mathbf{z}}^{i_1i_2} \left( \sum _{j_1} E_{j_1}^{k_2} {D^{-1\dag }}^{j_1k_1} \right) \left( \sum _{j_4} {^Z\mathbf{o}}^{i_4}_{j_4} {D^{-1\dag }}^{j_4k_4} {D^{-1\dag }}^{i_3k_3} \right) \nonumber \\&\quad + P(i_1/i_3i_4) P(k_1k_2/k_3k_4) {\mathbf{z}}^{i_1i_2} {\mathbf{z}}^{i_3i_4} \left( \sum _{j_1} E_{j_1}^{k_2} {D^{-1\dag }}^{j_1k_1} \right) \left( \sum _{j_4j_3} {D^{-1\dag }}^{j_3k_3} {^Z\mathbf{o}}_{j_3j_4} {D^{-1\dag }}^{j_4k_4} \right) \nonumber \\&\quad + P(i_1i_2/i_3i_4) P(k_1/k_3k_4) {\mathbf{z}}^{i_1i_2}{^Z\mathbf{o}}^{i_3i_4} \left( \sum _{j_3} E_{j_3}^{k_4} {D^{-1\dag }}^{j_3k_3} \right) \left( \sum _{j_1} E_{j_1}^{k_2} {D^{-1\dag }}^{j_1k_1} \right) \nonumber \\&\quad + P(i_1/i_3i_4) P(k_1/k_3k_4) {\mathbf{z}}^{i_1i_2} {\mathbf{z}}^{i_3i_4} {^Z\mathbf{O}}^{00} \left( \sum _{j_3} E_{j_3}^{k_4} {D^{-1\dag }}^{j_3k_3} \right) \left( \sum _{j_1} E_{j_1}^{k_2} {D^{-1\dag }}^{j_1k_1} \right) . \end{aligned}$$

(164c)

1.3 E.3 Type-3 matrix elements

The third category of matrix elements is obtained from $O^{\tilde{\sigma }}_{pIqJ}$ whenever

• 1. O is a zero-body operator, i.e. the identity operator multiplied by the number $O^{00}$,

• 2. $I=0,S,D$,

• 3. $J=0,S,D$.

All these matrix elements can be deduced from the previous cases by solely keeping the terms proportional to ${^Z\mathbf{O}}^{00} = O^{00}$ in the appropriate expressions.