Efficient computation of matrix elements of generic Slater determinants
Javier Rodríguez-Laguna, Luis Miguel Robledo, Jorge Dukelsky
EEfficient computation of matrix elements of generic Slater determinants
Javier Rodriguez-Laguna, Luis Miguel Robledo,
2, 3 and Jorge Dukelsky Dpto. F´ısica Fundamental, Universidad Nacional de Educaci´on a Distancia (UNED), Madrid, Spain. Dpto. F´ısica Te´orica, Universidad Aut´onoma de Madrid (UAM), Madrid, Spain. Center for Computational Simulation, Universidad Polit´ecnica de Madrid,Campus de Montegancedo, Boadilla del Monte, 28660-Madrid, Spain Instituto de Estructura de la Materia (IEM), CSIC, Madrid, Spain. (Dated: September 9, 2019)We present an extension of the L¨owdin strategy to find arbitrary matrix elements of generic Slaterdeterminants. The new method applies to arbitrary number of fermionic operators, even in the caseof a singular overlap matrix.
I. INTRODUCTION
Many developments in quantum many body physicsrequire the efficient computation of matrix elements offermionic operators between Slater determinant states.In many relevant cases, the matrix elements should becomputed between Slater determinant states which arenot based on the same set of single-body fermionic states.An early example of this kind of calculations can be foundin the seminal paper by the Swedish physicist Per-OlovL¨owdin [1] who developed in 1955 a very smart strategy,based on a careful application of the properties of thedeterminants. The strategy was developed in full for thecase of two and four fermionic operators, and requiredthe overlap matrix between the two Slater determinantsto be non-singular [1]. Further developments, aimed tosimplify the original complexity of the formulas and tofacilitate their use in the framework of the valence bondtheory, can be found in the literature [2, 3]. There aremany applications in quantum chemistry requiring over-laps of operators between Slater determinant includingConfiguration Interaction (CI) and symmetry restorationmethods (see Refs [4–9] as an example). More recentlyBrouder [10] proposed a method to reduce the combina-torial complexity of Wick’s theorem to a more manage-able algebraic complexity. As applications of the method,formulas for the overlap of a general product of creationand annihilation operators between arbitrary Slater de-terminants were proposed and used to compute, for in-stance, the generating function of the Green function or k -density correlation operators. The method uses gen-eral ideas coming from the world of quantum groups andHopf algebras but leads to rather involved expressions.Slater determinants are also used to expand the wavefunctions of the fractional quantum Hall effect (FQHE)as discussed in Ref [11]. In lattice QCD the study of phys-ical systems involving several hadrons [12] as it is the casein the description of collisions [13] or when the hadronsaggregate to form atomic nuclei [14–16], require the eval-uation of matrix elements involving the product of (3 N ) ( N is the number of hadrons) creation or annihilationoperators coming from the quark fields. If the matrix el-ements are computed in terms of Wick contractions (seebelow), the required number of terms grows exponen- tially fast with N [12]. General overlaps of Slater deter-minants are also required in nuclear physics, in the frame-work of the MonteCarlo Shell Model (MCSM) [17, 18], orin the field of symmetry restored quasiparticle excitations[19, 20].In all the above situations the matrix elements can becomputed with the help of Wick’s theorem or its gen-eralizations, but the number of contractions to considergrows with the factorial of the number of operators in-volved and therefore it becomes unmanageable very soon.In this work we extend L¨owdin’s results to the caseof a generic number of fermionic operators in order toobtain compact and easily handled expressions prone toan efficient evaluation in a computer. In addition, theevaluation of hundred of thousands, if not millions of op-erator overlaps, calls for robust evaluation methods ca-pable to handle the cases of zero or nearly zero overlapsof the Slater determinants where L¨owdin’s method be-comes ill-defined. In our derivation we will use a secondquantization formalism from the beginning, which makescalculations more transparent.This article is organized as follows. Our generalizedversion of L¨owdin’s theorem is described and proved inSec. II. The case of zero overlap between both Slaterdeterminants is discussed in Sec. III. A numerical appli-cation is provided in Sec. IV, where our approach is usedto estimate the entanglement of a block for a linear com-bination of Slater determinants. The article finishes withsome conclusions and our proposals for further work. II. GENERALIZED L ¨OWDIN’S THEOREM
To start with, let us consider two generic Slater deter-minants | A (cid:105) = a † i · · · a † i N |−(cid:105) , (1) | B (cid:105) = b † j · · · b † j N |−(cid:105) . (2)of N particles. The a † i and b † j are arbitrary creation oper-ators with quantum numbers denoted by i and j , respec-tively. Their hermitian conjugate annihilates the trueFock vacuum |−(cid:105) : a i |−(cid:105) = b j |−(cid:105) = 0 , (3) a r X i v : . [ c ond - m a t . s t r- e l ] S e p and such that { a † i , b j } = S ∗ ij , (4) { a i , b † j } = S ij , (5)as well as { a † i , b † j } = 0 . (6)The overlap matrix is defined by S ij = (cid:104) a i | b j (cid:105) = (cid:104)−| a i b † j |−(cid:105) . The overlap between both states, (cid:104) A | B (cid:105) isevaluated in a recursive way: (cid:104) A | B (cid:105) = (cid:104)−| a N · · · a b † · · · b † N |−(cid:105) (7)= − (cid:104)−| a N · · · a b † a b † · · · b † N |−(cid:105) (8)+ S (cid:104)−| a N · · · a b † · · · b † N |−(cid:105) (9)by jumping with the b † creation operator over the a annihilation one. The notation has been also simplifiedby replacing indexes i , . . . by 1 , . . . . Let us now introducethe quantity (cid:104) A | B (cid:105) [11] = (cid:104)−| a N · · · a b † · · · b † N |−(cid:105) , (10)that corresponds to the overlap of the two Slater determi-nants, but “removing the a b † pair from (cid:104) A | B (cid:105) ”. Then, (cid:104) A | B (cid:105) = S (cid:104) A | B (cid:105) [11] − S (cid:104) A | B (cid:105) [21] + · · · + ( − N +1 S N (cid:104) A | B (cid:105) [ N , (11)and the expansion ends after N jumps because (cid:104)−| b † = 0.We easily recognize in (cid:104) A | B (cid:105) [11] the minor of S with re-spect to the matrix element (1 , S . Viewed fromthis perspective, the expression of (cid:104) A | B (cid:105) given in Eq.(11) becomes the minor expansion of the determinant of S by the first row, i.e. (cid:104) A | B (cid:105) = det( S ) . (12)Let us now expand a † i and b † j in terms of a common basis { c † k , k = 1 , . . . , N B } a † i = N B (cid:88) k =1 A ki c † k , (13) b † j = N B (cid:88) k =1 B kj c † k , (14)then the N × N overlap matrix S becomes the productof the two expansion matrices, A and B , of dimension N B × N S ij = (cid:88) k A ∗ ki B kj = ( A † B ) ij . (15)The previous result can be easily generalized to the cal-culation of a general overlap (cid:104) A | f M · · · f g † · · · g † M | B (cid:105) , (16) where the f l and g † p are arbitrary annihilation and cre-ation operators expressed in the c † basis as f l = N B (cid:88) k =1 F kl c k , (17) g † p = N B (cid:88) k =1 G kp c † k , (18)in terms of the F and G matrices of dimension N B × N .This kind of overlaps appear when considering a systemof N + M particles where M of them play a differentrole than the remaining N ones and therefore require ofa different set of orbitals. As the f ’s anti-conmute withthe a ’s and the g † ’s with the b † , we can repeat verbatimthe previous considerations for (cid:104) A | B (cid:105) . We only have tobe careful and define four partial overlap matrices:( S fg ) ij = (cid:104)−| f i g † j |−(cid:105) ( M × M ) , (19)( S ag ) ip = (cid:104)−| a l g † j |−(cid:105) ( N × M ) , (20)( S fb ) lj = (cid:104)−| f i b † p |−(cid:105) ( M × N ) , (21)( S ab ) lp = (cid:104)−| a l b † p |−(cid:105) ( N × N ) , (22)to arrive to the formula (cid:104) A | f M · · · f g † · · · g † M | B (cid:105) = det (cid:18) S fg S fb S ag S ab (cid:19) (23)which is the general result for the overlap of Eq. (16). Inorder to disentangle the contributions from each set oforbitals it is convenient to use the well known formulafor the determinant of a partitioned matrixdet (cid:18) P QR S (cid:19) = det P det( S − RP − Q )= det S det( P − QS − R ) . (24)in order to obtaindet (cid:18) S fg S fb S ag S ab (cid:19) = det S ab det (cid:0) S fg − S fb S − S ag (cid:1) (25)which we can call generalized L¨owdin’s theorem (GLT). Itrequires the evaluation of the determinant of one M × M matrix and the determinant and inverse of a N × N ma-trix. This formula is also advantageous over Eq. (23)when N (cid:29) M (cid:29) g or f orbitals are required as only one costly matrixinversion is required. A similar expression has been ob-tained for the more general kind of product wave func-tions of the Hartree-Fock-Bogoliubov (HFB) type [21] us-ing pfaffians [22].In the right hand side of Eq. (25) a potential source ofproblems is identified in the inverse of S ab . If the inverseexists, then det S ab (cid:54) = 0 and it is possible to write (cid:104) A | f M · · · f g † · · · g † M | B (cid:105)(cid:104) A | B (cid:105) = det (cid:0) S fg − S fb S − S ag (cid:1) (26)which is the canonical form of the GLT where the sumof N ! contractions is replaced by the evaluation of thedeterminant of a M × M matrix. On the other hand,the result of Eq. (25) is required to resolve the implicitindeterminacy when det S ab = 0 and S ab is not invertible(see below).The above derivation assumes that the f and g † are innormal order. If this is not the case, operators can alwaysbe brought to normal order using commutation relationsof fermion operators. To illustrate the procedure and toobtain a compact expression we evaluate now the overlapof a one-body operator ˆ Q (cid:104) A | f M · · · f ˆ Qg † · · · g † M | B (cid:105) (27)where ˆ Q is written in terms of fermion operators r † m and t n as ˆ Q = (cid:88) m,n Q mn r † m t n . (28)The matrix element (cid:104) A | f M · · · f r † m t n g † · · · g † M | B (cid:105) isevaluated by using the commutation relation r † m t n = − t n r † m + ( S tr ) nm as (cid:104) A | f M · · · f r † m t n g † · · · g † M | B (cid:105) =( S tr ) mn det (cid:18) S fg S fb S ag S ab (cid:19) − det S tr S tg S tb S fr S fg S fb S ar S ag S ab (29)With obvious notation, we introduce the row S t,gb andcolumn S fa,r vectors as well as the matrix S fa,gb to be ableto use property (24). Straightforward manipulations leadto the final expression (cid:104) A | f M · · · f r † m t n g † · · · g † M | B (cid:105) = S t,gb S − S fa,r det ( S fa,gb ) (30)For the evaluation of the overlap of a two body operatorthe matrix element (cid:104) A | f M · · · f r † m r † m t n t n g † · · · g † M | B (cid:105) (31)is required. Using the same procedure as before and aftera few manipulations we obtain (cid:104) A | f M · · · f r † m r † m t n t n g † · · · g † M | B (cid:105) =det (cid:16) S t,gb S − S fa,r (cid:17) det ( S fa,gb ) (32)where S t,gb and S fa,r have dimensions 2 × ( M + N ) and( M + N ) ×
2, respectively. The matrix element is givenby the product of det( S fa,gb ) times the determinant of a2 × k -particle, k -hole matrix elements is straightforward andleads to the determinant of a k × k matrix of contractions.The combinatorial increase in the number of terms as k increases is thus hidden in the form of a determinant of low dimensionality. This result is the generalization ofEq. (51) of [23].To finish this section let us consider a common sit-uation concerning to symmetry restoration where theoverlap includes a multi-particle unitary operator ˆ T inthe form of an exponentiated one body-operator. Typ-ical examples are the rotation and translation operator.Then the operator b † m generating the | B (cid:105) configurationare transformed to ˜ b † m given byˆ T b † m ˆ T † = ˜ b † m = (cid:88) n ( T b ) nm b † n (33)The overlap becomes (cid:104) A | f M · · · f g † · · · g † M ˆ T | B (cid:105) = det (cid:18) S fg S f˜b S ag S a˜b (cid:19) (34)where the only modification with respect to Eq. (23) isin the overlaps S f˜b and S a˜b which have to be computedwith the ˜ b † m of Eq. (33). III. CASE OF ZERO OVERLAP
Let us study how to apply the GLT of Eq. (25) whenthe overlap between the states | A (cid:105) and | B (cid:105) is zero.The methodology used can be used straightforwardly forthe other form of L¨owdin’s theorem, Eq. (32). When (cid:104) A | B (cid:105) = 0, S ab is a singular matrix and Eq. (25) be-comes indeterminate. To avoid the problem one canalways use the full determinant in Eq. (23), of order( N + M ) × ( N + M ), but this comes at a higher cost thanjust using (25). In addition, resolving the indeterminacyexplicitly is always beneficial in order to avoid numericalartifacts that could eventually appear. To this end, weintroduce the singular value decomposition (SVD) of SS ab = U Σ V † , (35)where U and V are orthogonal matrices ( U † U = I ) andΣ is diagonal. If S ab is near singular, it can be expressedas Σ = σ . . . σ N − k ε . . . ε k ≡ (cid:18) Σ R E (cid:19) (36)with ε i a set of k small numbers, while E is the k × k diagonal matrix with ε i in the diagonal. Using thisdecomposition we arrive atdet S ab = f det Σ , (37)where f = e iϕ UV ≡ det U det V † . One also has S − = V Σ − U † . (38)Let us define now S V fb ≡ S fb V and S U ag ≡ U † S ag , anddecompose them in a regular (R) and a singular (S) part,according to the decomposition in Eq. (36) S V fb = (cid:0)(cid:124) (cid:123)(cid:122) (cid:125) N − k ¯ S V, Rfb (cid:124)(cid:123)(cid:122)(cid:125) k ¯ S V, Sfb (cid:1)(cid:9) M (39)and S U ag = (cid:32) ¯ S U, Rag (cid:124)(cid:123)(cid:122)(cid:125) M ¯ S U, Sag (cid:33) } N-k } k (40)Using this decomposition we can writedet (cid:18) S fg S fb S ag S ab (cid:19) = f det Σ det (cid:0) S fg − S V fb Σ − S U ag (cid:1) , (41)with S V fb Σ − S U ag = ¯ S V, Rfb (cid:0) Σ R (cid:1) − ¯ S U, Rag + ¯ S V, Sfb E − ¯ S U, Sag , (42)Let us now introduce the M × M matrix C ≡ S fg − ¯ S V, Rfb (cid:0) Σ R (cid:1) − ¯ S U, Rag , (43)and consider the determinantdet (cid:16) C − ¯ S V, Sfb E − ¯ S U, Sag (cid:17) . (44)It can be computed using the Woodbury formula for thedeterminant (see, for instance, Ref [24])det (cid:0) A + U W V † (cid:1) = det (cid:0) W − + V † A − U (cid:1) det W det A (45)to obtaindet (cid:18) S fg S fb S ag S ab (cid:19) = f det Σ R det E ×× det (cid:16) E − ¯ S U, Sag C − ¯ S V, Sfb (cid:17) E det C, (46)which is well defined in the limit ε i → (cid:18) S fg S fb S ag S ab (cid:19) = f det Σ R det (cid:16) − ¯ S U, Sag C − ¯ S V, Sfb (cid:17) det C, (47)which is a finite quantity when ε i →
0. This quantity re-quires the SVD of S ab to get Σ R and the U and V matri-ces. For large values of N this can be a costly operation oforder N . The inverse of the diagonal Σ R is also requiredas well as the construction of the S U ag and S V fb matrices.Once we have all the ingredients, the evaluation of thedeterminants in Eq. (47) require little extra work dueto the low dimensionality of the matrices involved ( k × k and M × M ). Please note that the dimensionality of thedifferent matrices appearing in Eq. (47) is not the same(Σ R is ( N − k ) × ( N − k ), ¯ S U, Sag C − ¯ S V, Sfb is k × k and C is M × M ) and therefore the formula for the product ofa determinant does not apply here.Note that the derivation above can also be extendedto the case where the matrix S ab is ill-conditioned andhas a very small (but non-zero) determinant. A blind useof the traditional formulas may contaminate the final re-sults due to numerical artifacts consequence of the finiterepresentation of floating point numbers in computer’sarithmetic. IV. NUMERICAL EXPERIMENTS
In this section we put the extended L¨owdin approachto the test, characterizing the entanglement behavior ofa linear combination of two Slater determinants.
A. The rainbow system
As our physical system, we have chosen the rainbowsystem [25–27], a 1D inhomogeneous fermionic hoppingsystem which presents volumetric entanglement betweenits left and right halves. It can be described on an openchain through the following Hamiltonian, H = − L − (cid:88) m = − L +1 J m c † m c m +1 + h.c. , (48)with hopping amplitudes given by J m = (cid:40) α m +1 if m (cid:54) = 0,1 otherwise, (49)in terms of an inhomogeneity parameter α ∈ (0 , α = 1, the system reduces to the homogeneous case. Forsmall α , the ground state (GS) of Hamiltonian (48) be-comes approximately a valence bond solid with concentricbonds around the center, see Fig. 1. This GS violatesmaximally the area-law [28], giving rise to a volumetricgrowth of entanglement. B. Entanglement and number fluctuations
The eigenstates of Hamiltonian (48) are Slater deter-minants. Therefore, the exact entanglement entropy forany block can be efficiently evaluated [29]. Yet, this tech-nique of evaluating the entanglement entropy cannot beextended to a linear combination of Slater determinants.Let us consider an arbitrary linear combination of theground state and first excited states of Hamiltonian (48),both within the half-filling sector, i.e. with N/ | Ψ (cid:105) = α | (cid:105) + β | (cid:105) . (50) m : − − − − − Figure 1. Illustration of the rainbow ground state in the small α limit: a valence bond solid with concentric bonds aroundthe center of the chain. Notice that the entanglement en-tropy of a left block is proportional to the block size, i.e. avolumetric growth of the entropy. We are interested in estimates of the entanglement be-tween a region A and its complementary, always withinstate | Ψ (cid:105) . The procedure is straightforward in the caseof a single Slater determinants [29]. In that case, thecomputation of the full spectrum of the reduced densitymatrix can be performed in O ( N k ) operations, with alow value of k . On the other hand, the evaluation of theentanglement properties for a generic state takes, in prin-ciple, O (2 N ) operations. Yet, there are some observableswhich act as entanglement witnesses , i.e.: whose valuesact as an indicator of entanglement. The fluctuation ofthe number of particles is one of these observables [30].Let n A be the number of particles on region A : n A = (cid:88) i ∈ A n i = (cid:88) i ∈ A c † i c i , (51)whose fluctuations are given by σ N = (cid:10) n A (cid:11) − (cid:104) n A (cid:105) . (52)The expectation value of the number of particles on re-gion A is found through the following calculation: (cid:104) n A (cid:105) = (cid:88) i ∈ A (cid:104) Ψ | n i | Ψ (cid:105) = (cid:88) i ∈ A (cid:2) | α | (cid:104) | n i | (cid:105) + | β | (cid:104) | n i | (cid:105) + 2Re ¯ αβ (cid:104) | n i | (cid:105) ] . (53)The first two terms are easily found, because they referto a single Slater determinant. The third, notwithstand-ing, must be found using the (generalized) L¨owdin tricksdescribed: (cid:104) | n i | (cid:105) = (cid:104) | c † i c i | (cid:105) = (cid:104) | (cid:105) − (cid:104) | c i c † i | (cid:105) . (54)The quadratic term is more involved: (cid:10) n A (cid:11) = (cid:88) i,j ∈ A (cid:2) | α | (cid:104) | n i n j | (cid:105) + | β | (cid:104) | n i n j | (cid:105) +2Re ¯ αβ (cid:104) | n i n j | (cid:105) ] . (55)Again, the first two terms are straightforward to ob-tain. If C A, is the submatrix of the correlation matrixcorresponding to block A on state | (cid:105) , then (cid:104) | n i n j | (cid:105) = Tr( C A, ) + Tr( C A, ( I − C A, )) . (56)The last term of Eq. (55) is the most involved one,because we cannot assume Wick’s theorem. We find (cid:104) | n i n j | (cid:105) = (cid:104) | (1 − c i c † i )(1 − c j c † j ) | (cid:105) == (cid:104) | (cid:105) − (cid:104) | c j c † j | (cid:105) − (cid:104) | c i c † i | (cid:105) + δ ij (cid:104) | c i c † j | (cid:105) − (cid:104) | c i c j c † i c † j | (cid:105) . (57)It has been proved that, for a single Slater determinant[30], S A ≥ σ N . (58)Thus, even though Eq. (58) is not proved for a genericlinear combination of Slater determinants, we will employ σ N as our estimate for an entanglement witness . C. Numerical experiments
For concreteness, let us set α = √ x and β = √ − x for x ∈ [0 ,
1] in Eq. (50). Thus, our state will be givenby | Ψ (cid:105) = √ x | (cid:105) + √ − x | (cid:105) , (59)Let us notice that the first excitation is obtained fromthe ground state by performing a parity transformationon the Fermi level [27]. Notice that, by construction (cid:104) | (cid:105) = 0, thus forcing us to make all our computationsin the zero overlap case.In the top panel of Fig. 2 we show the variance σ N of the number of particles in the left half of the rain-bow system with N = 8, for different values of α . No-tice that only in the α → + limit the variance is thesame for x = 0 and x = 1, i.e.: in the GS and thefirst excited. This variance has been computed in twodifferent ways: the dots correspond to the exact calcula-tion, with the full Slater determinant, and the continuousline corresponds to the computation performed with thegeneralized L¨owdin formulas derived in this article. Theagreement is complete, and the computational time isenormously reduced with our tools. . . . . . . .
91 0 0 . . . . . . . . . σ N x α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . α = 1 . . .
53 0 0 . . . . . . . . . S x α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . α = 1 Figure 2. Top: deviation of the number of particles in theleft half of state (59) with N = 8 sites and various valuesof α , computed with the full Slater determinants and withour generalized L¨owdin scheme. The theoretical value in the α → + limit is σ A = 1. Bottom: entanglement of the lefthalf of state (59), computed with the full Slater determinants.The theoretical value in the α → + limit is 4 log 2 ≈ . In the α → + limit, the left-half of the rainbow systembecomes an infinite temperature mixed state. Thus, thefluctuations in the particle number are easy to obtain, fol-lowing a binomial distribution, σ A = (cid:112) N/
8, that we canreadily check in Fig. 2. Also, the entanglement entropywill grow up to the maximal possible value, ( N log 2) / x for various values of α in a chainwith N = 40 sites. D. Generalized L¨owdin C++ code
We have uploaded our C++ libraries to compute ma-trix elements of Slater determinants using the generalizedL¨owdin approach to the repository github as free soft- . . . . . . . . . . . . . . . . . σ N x α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . α = 1 Figure 3. Same observable as Fig. 2 for a rainbow systemwith N = 40. The theoretical value for the deviation of theparticle number of the left half in the α → + limit is (cid:112) N/ √ ≈ . ware [31]. The same material can also be downloadedfrom the Supplementary material section of the electronicversion of the journal [32]. V. CONCLUSIONS AND FURTHER WORK
We have extended the seminal results of L¨owdin to ob-tain efficiently the matrix elements of an arbitrarily largeproduct of fermionic operators between arbitrary Slaterdeterminants. Our results are still applicable when theoverlap matrix between the orbitals of the Slater deter-minants is singular, i.e. when the corresponding statesare orthogonal.Efficient computation of matrix elements in non-orthogonal Slater determinants will open a very interest-ing possibility: the creation of Ans¨atze including Slaterdeterminants obtained from different procedures and,therefore, using different orbitals.As proposals of future work, we would like to remarkthe extension of the previous calculations to the obten-tion of the full reduced density matrix, combining ourresults with those of [29] for the reduced density matrixof a block of a single Slater determinant.
ACKNOWLEDGMENTS
We thank G.F. Bertsch for a careful reading of themanuscript and suggestions. This work has been sup-ported by the Spanish Ministerio de Ciencia, Innovaci´ony Universidades and the European regional develop-ment fund (FEDER), grants Nos FIS2015-69167-C2-1(JRL), FIS2015-63770 (JD and LMR), PGC2018-094180-B-I00 (JD), and FPA2015-65929, PGC2018-094583-B-I00 (LMR). [1] P.-O. L¨owdin, Phys. Rev. , 1490 (1955).[2] F. Prosser and S. Hagstrom, International Journal ofQuantum Chemistry , 89 (1968).[3] S. C. Leasure and G. G. Balint-Kurti, Phys. Rev. A ,2107 (1985).[4] H. Koch and E. Dalgaard, Chemical Physics Letters ,193 (1993).[5] J. Verbeek and J. H. V. Lenthe, Journal of MolecularStructure: THEOCHEM , 115 (1991).[6] N. Tomita, S. Ten-no, and Y. Tanimura, ChemicalPhysics Letters , 687 (1996).[7] G. E. Scuseria, C. A. Jim´enez-Hoyos, T. M.Henderson, K. Samanta, and J. K. Ellis, TheJournal of Chemical Physics , 124108 (2011),http://dx.doi.org/10.1063/1.3643338.[8] C. A. Jim´enez-Hoyos, R. Rodr´ıguez-Guzm´an, and G. E.Scuseria, The Journal of Chemical Physics , 204102(2013), https://doi.org/10.1063/1.4832476.[9] R. Rodr´ıguez-Guzm´an, C. A. Jim´enez-Hoyos, and G. E.Scuseria, Phys. Rev. B , 195110 (2014).[10] C. Brouder, Phys. Rev. A , 032720 (2005).[11] B. A. Bernevig and N. Regnault, Phys. Rev. Lett. ,206801 (2009).[12] W. Detmold and K. Orginos, Phys. Rev. D , 114512(2013).[13] T. Yamazaki, Y. Kuramashi, and A. Ukawa (PACS-CSCollaboration), Phys. Rev. D , 111504 (2010).[14] S. Beane, W. Detmold, K. Orginos, and M. Savage,Progress in Particle and Nuclear Physics , 1 (2011).[15] M. Savage, Progress in Particle and Nuclear Physics ,140 (2012), from Quarks and Gluons to Hadrons andNuclei. [16] T. Doi and M. G. Endres, Computer Physics Communi-cations , 117 (2013).[17] T. Otsuka, M. Honma, T. Mizusaki, N. Shimizu, andY. Utsuno, Progress in Particle and Nuclear Physics ,319 (2001).[18] Y. Utsuno, N. Shimizu, T. Otsuka, and T. Abe, Com-puter Physics Communications , 102 (2013).[19] K. Schmid, Progress in Particle and Nuclear Physics ,565 (2004).[20] G. Puddu, Journal of Physics G: Nuclear and ParticlePhysics , 321 (2006).[21] G. F. Bertsch and L. M. Robledo, Phys. Rev. Lett. ,042505 (2012).[22] L. M. Robledo, Phys. Rev. C , 021302 (2009).[23] P.-O. L¨owdin, Phys. Rev. , 1474 (1955).[24] D. Harville, Matrix Algebra From a Statistician’s Per-spective (Springer New York, 2008).[25] G. Vitagliano, A. Riera, and J. I. Latorre, New Journalof Physics , 113049 (2010).[26] G. Ram´ırez, J. Rodr´ıguez-Laguna, and G. Sierra, Jour-nal of Statistical Mechanics: Theory and Experiment , P10004 (2014).[27] G. Ram´ırez, J. Rodr´ıguez-Laguna, and G. Sierra, Jour-nal of Statistical Mechanics: Theory and Experiment , P06002 (2015).[28] J. Eisert, M. Cramer, and M. B. Plenio, Rev. Mod. Phys. , 277 (2010).[29] I. Peschel, Journal of Physics A: Mathematical and Gen-eral , L205 (2003).[30] I. Klich, Journal of Physics A: Mathematical and General39