Efficient computation of matrix elements of generic Slater determinants
Javier Rodríguez-Laguna Luis Miguel Robledo Jorge Dukelsky












apsrev4-1


<
>
||⟩#1>|#1>
⟨<|#1|<#1|
<#1|#2|#3><#1|#2|#3>
⟨<|#⟩1|#2>⟨#1|#2⟩
()

    [
]

∙ 
ℤ
ℂ
Tr
eff

    #1,#2,#3,#4.[ #1 #2; #3 #4 ]#1,#2,#3,#4,#5,#6,#7,#8,#9.[ #1 #2 #3; #4 #5 #6; #7 #8 #9 ]fgagfbabShort Title]Efficient computation of matrix elements of
      generic Slater determinantsDpto. Física Fundamental, Universidad Nacional de
      Educación a Distancia (UNED), Madrid, Spain.Dpto. Física Teórica, Universidad Autónoma de
      Madrid (UAM), Madrid, Spain.Center for Computational Simulation, Universidad Politécnica de Madrid,
      Campus de Montegancedo, 
      Boadilla del Monte, 28660-Madrid, SpainInstituto de Estructura de la Materia (IEM), CSIC,
      Madrid, Spain.
    We present an extension of the Löwdin strategy to find arbitrary
    matrix elements of generic Slater determinants. 
    The new method applies to arbitrary number of fermionic operators, even in the case of a singular overlap matrix. 
    [
        Jorge Dukelsky
        September 9, 2019
    =====================
    
    
    
    § INTRODUCTION
    
    Many developments in quantum many body physics require the efficient
    computation of matrix elements of fermionic operators between Slater
    determinant states.  In many relevant cases, the matrix elements
    should be computed between Slater determinant states which are not
    based on the same set of single-body fermionic states. An early
    example of this kind of calculations can be found in the seminal paper
    by the Swedish physicist Per-Olov Löwdin <cit.> who developed in 1955 a
    very smart strategy, based on a careful application of the properties
    of the determinants. The strategy was developed in full for the case
    of two and four fermionic operators, and required the overlap matrix
    between the two Slater determinants to be non-singular
    <cit.>.  Further developments, aimed to simplify the
    original complexity of the formulas and to facilitate their use in the
    framework of the valence bond theory, can be found in the literature
    <cit.>. There are many applications in
    quantum chemistry requiring overlaps of operators between Slater
    determinant including Configuration Interaction (CI) and symmetry
    restoration methods (see Refs
    <cit.>
    as an example).  More recently Brouder <cit.>
    proposed a method to reduce the combinatorial complexity of Wick's
    theorem to a more manageable algebraic complexity. As applications of
    the method, formulas for the overlap of a general product of creation
    and annihilation operators between arbitrary Slater determinants were
    proposed and used to compute, for instance, the generating function of
    the Green function or k-density correlation operators. The method
    uses general ideas coming from the world of quantum groups and Hopf
    algebras but leads to rather involved expressions. Slater determinants
    are also used to expand the wave functions of the fractional quantum
    Hall effect (FQHE) as discussed in Ref <cit.>.
    In lattice QCD the study of physical systems involving several hadrons
    <cit.> as it is the case in the description of collisions
    <cit.> or when the hadrons aggregate to form atomic
    nuclei <cit.>, require the evaluation
    of matrix elements involving the product of (3N)^2 (N is the
    number of hadrons) creation or annihilation operators coming from the
    quark fields. If the matrix elements are computed in terms of Wick
    contractions (see below), the required number of terms grows
    exponentially fast with N<cit.>. General overlaps of
    Slater determinants are also required in nuclear physics, in the
    framework of the MonteCarlo Shell Model (MCSM)
    <cit.>, or in the field of symmetry
    restored quasiparticle excitations <cit.>.
    
    In all the above situations the matrix elements can be computed with the
    help of Wick's theorem or its generalizations, but the number of contractions
    to consider grows with the factorial of the number of operators involved and 
    therefore it becomes unmanageable very soon.  
    
    In this work we extend Löwdin's results to the case of a generic
    number of fermionic operators in order to obtain compact and easily 
    handled expressions prone to an efficient evaluation in a computer. In
    addition, the evaluation of hundred of thousands, if not millions of
    operator overlaps, calls for robust evaluation methods capable to
    handle the cases of zero or nearly zero overlaps of the Slater
    determinants where Löwdin's method becomes ill-defined. In our derivation we will use
    a second quantization formalism from the beginning, which makes
    calculations more transparent.
    
    This article is organized as follows. Our generalized version of
    Löwdin's theorem is described and proved in Sec. <ref>. The
    case of zero overlap between both Slater determinants is discussed in
    Sec. <ref>. A numerical application is provided in
    Sec. <ref>, where our approach is used to estimate the
    entanglement of a block for a linear combination of Slater
    determinants. The article finishes with some conclusions and our
    proposals for further work.
    
    
    
    
    § GENERALIZED LÖWDIN'S THEOREM
    
    
    To start with, let us consider two generic Slater determinants
    
    
        ||⟩A> =    a^†_i_1⋯ a^†_i_N||⟩->, 
        ||⟩B> =    b^†_j_1⋯ b^†_j_N||⟩->.
    
    of N particles. The a^†_i and b^†_j are arbitrary
    creation operators with quantum numbers denoted by i and j, respectively.
    Their hermitian conjugate annihilates the true
    Fock vacuum ||⟩->:
    
    
        a_i||⟩->=b_j ||⟩-> = 0,
    
    and such that
    
    
        { a^†_i, b_j}   = S^*_ij, 
        { a_i, b^†_j}   = S_ij,
    
    as well as
    
    
        {a^†_i,b^†_j}=0.
    
    The overlap matrix is defined by S_ij=a_i|b_j=⟨<|-|a_i
    b^†_j||⟩->. The overlap between both states, A|B is
    evaluated in a recursive way:
    
    
        A|B =    ⟨<|-| a_N ⋯ a_1 b^†_1 ⋯ b^†_N ||⟩-> 
        
          =    - ⟨<|-| a_N ⋯ a_2 b^†_1 a_1
                       b^†_2 ⋯ b^†_N ||⟩-> 
             + S_11⟨<|-| a_N ⋯ a_2 b^†_2 ⋯ b^†_N ||⟩->
    
    by  jumping with the b^†_1 creation operator over the
    a_1 annihilation one. The notation has been also simplified by
    replacing indexes i_1, … by 1,….  Let us now introduce
    the quantity
    
    
        A|B_[11]=⟨<|-| a_N ⋯ a_2 b^†_2⋯ b^†_N
          ||⟩->,
    
    that corresponds to the overlap of the two Slater determinants, but
    “removing the a_1b^†_1 pair from A|B”. Then,
    
    
        A|B   = S_11A|B_[11] -S_21A|B_[21] +⋯
            +(-1)^N+1S_N1A|B_[N1],
    
    and the expansion ends after N jumps because ⟨<|-|b^†_1 =0.
    We easily recognize in A|B_[11] the minor of S with respect 
    to the matrix element (1,1), i.e. S_11. Viewed from this perspective, 
    the expression of A|B given in Eq. (<ref>)  
    becomes the minor expansion of the determinant of S by the first row, i.e.
    
    
        A|B=(S).
    
    Let us now expand a^†_i and b^†_j in terms of a common
    basis {c^†_k,k=1,…,N_B}
    
        a^†_i    = ∑_k=1^N_B A_ki c^†_k,  
        
          b^†_j    = ∑_k=1^N_B B_kj c^†_k,
    
    then the N × N overlap matrix S becomes the product of the
    two expansion matrices, A and B, of dimension N_B × N
    
        S_ij=∑_k A^*_ki B_kj = (A^† B)_ij.
    
    The previous result can be easily generalized to the calculation of a general
    overlap
    
    
        ⟨<|A| f_M ⋯ f_1 g^†_1 ⋯ g^†_M ||⟩B>,
    
    where the f_l and g^†_p are arbitrary annihilation and
    creation operators expressed in the c^† basis as
    
    
        f_l           = ∑_k=1^N_B F_kl c_k,  
        
          g^†_p    = ∑_k=1^N_B G_kp c^†_k,
    
    in terms of the F and G matrices of dimension N_B × N. This kind
    of overlaps appear when considering a system of N+M particles where
    M of them play a different role than the remaining N ones and
    therefore require of a different set of orbitals. 
    
    As the f's anti-conmute with the a's and the g^†'s with the
    b^†, we can repeat verbatim the previous considerations for
    A|B. We only have to be careful and define four partial overlap
    matrices: 
    
    
    
    
        (S_)_ij = ⟨<|-| f_i g^†_j ||⟩->         (M× M), 
        
          (S_)_ip = ⟨<|-| a_l g^†_j ||⟩->         (N× M), 
        
          (S_)_lj = ⟨<|-| f_i b^†_p ||⟩->         (M× N), 
        
          (S_)_lp = ⟨<|-| a_l b^†_p ||⟩->         (N× N),
    
    to arrive to the formula
    
    
        ⟨<|A| f_M ⋯ f_1 g^†_1 ⋯ g^†_M ||⟩B> =
           S_,S_,S_,S_.
    
    which is the general result for the overlap of Eq. (<ref>). 
    In order to disentangle the contributions from each set of orbitals it 
    is convenient to use the
    well known formula for the determinant of a partitioned matrix
    
    
        P,Q,R,S.    =  P ( S-RP^-1Q ) 
           =  S  (P-QS^-1R).
    
    in order to obtain
    
    
        S_,S_,S_,S_. =  S_S_-S_ S_^-1 S_
    
    which we can call  generalized Löwdin's theorem (GLT). It requires the
    evaluation of the determinant of one M× M matrix and the determinant
    and inverse of a N× N matrix. This formula is also advantageous 
    over Eq. (<ref>) when N≫ M≫ 1
    and many matrix elements for different g or f orbitals are required as
    only one costly matrix inversion is required.
    
    
    
    
    
    
    A similar expression has been obtained for the more general kind of product wave functions of the Hartree-Fock-Bogoliubov
    (HFB) type <cit.> using pfaffians <cit.>. 
    
    In the right hand side of 
    Eq. (<ref>) a potential source of problems is identified in the
    inverse of S_. If the inverse exists, then S_≠ 0 and 
    it is possible to write
    
    
        ⟨<|A| f_M ⋯ f_1 g^†_1 ⋯ g^†_M ||⟩B>/A|B = S_-S_ S_^-1 S_
    
    which is the canonical form of the GLT where the sum of N!
    contractions is replaced by the evaluation of the determinant of a
    M× M matrix. On the other hand, the result of
    Eq. (<ref>) is required to resolve the implicit indeterminacy
    when S_ = 0 and S_ is not invertible (see below).
    
    
    The above derivation assumes that the f and g^† are in
    normal order. If this is not the case, operators can always be brought
    to normal order using commutation relations of fermion operators. To
    illustrate the procedure and to obtain a compact expression we
    evaluate now the overlap of a one-body operator Q̂
    
        ⟨<|A| f_M ⋯ f_1 Q̂ g^†_1 ⋯ g^†_M ||⟩B>
    
    where Q̂ is written in terms of fermion operators
    r^†_m and t_n as
    
    
        Q̂ = ∑_m,n Q_mn r^†_m t_n.
    
    The matrix element ⟨<|A| f_M ⋯ f_1 r^†_m t_n
    g^†_1 ⋯ g^†_M ||⟩B> is evaluated by using the
    commutation relation r^†_m t_n= - t_n r^†_m +
    S_tr_nm as
    
    
    
        ⟨<|A| f_M ⋯ f_1 r^†_m t_n g^†_1 ⋯ g^†_M ||⟩B> = 
        S_tr_mn S_,S_,S_,S_. - 
           S_tr,S_tg,S_tb,S_fr,S_,S_,S_ar,S_,S_.
    
    With obvious notation, we introduce the row  S_t,gb and column
    S_fa,r vectors as well as the matrix S_fa,gb to be able to use
    property (<ref>). Straightforward 
    manipulations lead to the final expression
    
    
        ⟨<|A| f_M ⋯ f_1 r^†_m t_n g^†_1 ⋯ g^†_M ||⟩B> = 
        
          S_t,gb S_fa,gb^-1 S_fa,rS_fa,gb
    
    For the evaluation of the overlap of a two body operator the matrix element
    
    
        ⟨<|A| f_M ⋯ f_1 r^†_m_1 r^†_m_2 t_n_1t_n_2 g^†_1 ⋯ g^†_M ||⟩B>
    
    is required. Using the same procedure as before and after a few
    manipulations we obtain 
    
    
        ⟨<|A| f_M ⋯ f_1 r^†_m_1 r^†_m_2 t_n_1t_n_2 g^†_1 ⋯ g^†_M ||⟩B> = 
        S_t,gb S_fa,gb^-1 S_fa,rS_fa,gb
    
    where S_t,gb and S_fa,r have dimensions 2× 
    (M+N)  and (M+N)× 2, respectively. The matrix element is given 
    by the product of ( S_fa,gb ) times the determinant 
    of a 2× 2 matrix with entries corresponding to the “elementary 
    contractions". The generalization to more general k-particle, 
    k-hole matrix elements is straightforward and leads 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 of Eq. (51) of <cit.>. 
    
    
    
    To finish this section let us consider a common situation concerning
    to symmetry restoration where the overlap includes a multi-particle
    unitary operator 𝒯̂ in the form of an exponentiated one
    body-operator. Typical examples are the rotation and translation operator. Then the
    operator b^†_m generating the |B⟩ configuration are transformed to
    b̃^†_m given by
    
    
    
        𝒯̂ b^†_m 𝒯̂^† =
          b̃^†_m = ∑_n T_b_nm b^†_n
    
    The overlap becomes
    
    
        ⟨<|A| f_M ⋯ f_1  g^†_1 ⋯ g^†_M 𝒯̂||⟩B> =
           S_,S_fb̃,S_,S_ab̃.
    
    where the only modification with respect to Eq. (<ref>) is in the 
    overlaps S_fb̃ and S_ab̃ which have to be computed
    with the b̃^†_m of Eq. (<ref>).
    
    
    
    
    
    § CASE OF ZERO OVERLAP
    
    
    Let us study how to apply the GLT of Eq. (<ref>) when the overlap
    between the states ||⟩A> and ||⟩B> is zero. The methodology
    used can be used straightforwardly for the other form of Löwdin's
    theorem, Eq. (<ref>). When A|B=0, S_ab is a
    singular matrix and Eq. (<ref>) becomes indeterminate. To avoid
    the problem one can always use the full determinant in Eq.
    (<ref>), of order (N+M)×(N+M), but this comes at a
    higher cost than just using (<ref>). In addition, resolving the
    indeterminacy explicitly is always beneficial in order to avoid
    numerical artifacts that could eventually appear. To this end, we introduce
    the  singular value decomposition (SVD) of S
    
        S_=U Σ V^†,
    
    where U and V are orthogonal matrices (U^† U=I) and
    Σ is diagonal. If S_ is near singular, it can be expressed as
    
    ε
    
        Σ = [   σ_1      ;           ⋱      ;             σ_N-k;                      _1;                             ⋱      ;                                  _k ]≡[ Σ^R    ;       E ]
    
    with _i a set of k small numbers, 
    while E is the k× k diagonal matrix with _i in the
    diagonal. Using this decomposition we arrive at
    
    
        S_= f Σ,
    
    where  f=e^iφ_UV≡ U  V^†.
    
    
    
    
    
    
    One also has
    
    
        S_^-1=V Σ^-1 U^†.
    
    Let us define now S^V_≡ S_ V and S^U_≡
    U^† S_, and decompose them in a  regular (R) and a
     singular (S) part, according to the decomposition in
    Eq. (<ref>)
        S^V_= .[ N-kS̅^V,R_   kS̅^V,S_ ]} M
    
    and 
    
    
        S^U_ = (
          [  S̅^U,R_; MS̅^U,S_;          ])
            [ .} N-k;   .} k;        ]
    
    
    
    
    Using this decomposition we can write
    
    
        S_,S_,S_,S_. =
          f ΣS_ - S^V_Σ^-1 S^U_,
    
    with
    
    
        S^V_Σ^-1 S^U_ =    S̅^V,R_Σ^R^-1S̅^U,R_
         
           +   S̅^V,S_ E^-1S̅^U,S_,
    
    Let us now introduce the M× M matrix
    
    
        C ≡ S_ - S̅^V,R_Σ^R^-1S̅^U,R_,
    
    and consider the determinant
    
    
        C-S̅^V,S_ E^-1S̅^U,S_.
    
    It can be computed using the Woodbury formula for the determinant 
    (see, for instance, Ref <cit.>)
    
    
        A+UWV^† = W^-1+V^† A^-1 U W
           A
    
    to obtain
    
    
        S_,S_,S_,S_. =
          f Σ^R E ×
           ×E - S̅^U,S_ C^-1S̅^V,S_1/ E C,
    
    which is well defined in the limit _i→ 0
    
        S_,S_,S_,S_. = fΣ^R-S̅^U,S_ C^-1S̅^V,S_ C,
    
    which is a finite quantity when _i→ 0. This quantity requires the 
    SVD of S_ to get Σ^R and the U and V matrices.
    For large values of N this can be a costly operation of order N^3.
    The inverse of the diagonal Σ^R is also required
    as well as the construction of the  S^U_ and
    S^V_ matrices. Once we have all the ingredients, the evaluation
    of the determinants in Eq. (<ref>) require little extra work
    due to the low dimensionality of the matrices involved (k × k and
    M× M). Please note that the dimensionality of the different matrices
    appearing in Eq. (<ref>) is not the same (Σ^R
    is (N-k)×(N-k), S̅^U,S_ C^-1S̅^V,S_
    is k× k and C is M× M) and therefore the formula for the product
    of a determinant does not apply here. 
        
    Note that the derivation above can also be extended to the case where
    the matrix S_ is ill-conditioned and has a very small (but
    non-zero) determinant. A blind use of the traditional
    formulas may contaminate the final results due to numerical artifacts
    consequence of the finite representation of floating point numbers in
    computer's arithmetic.
    
    
    
    § NUMERICAL EXPERIMENTS
    
    
    In this section we put the extended Löwdin approach to the test,
    characterizing the entanglement behavior of a linear combination of
    two Slater determinants.
    
    
    
     §.§ The rainbow system
    
    
    As our physical system, we have chosen the  rainbow system<cit.>, a 1D inhomogeneous
    fermionic hopping system which presents volumetric entanglement
    between its left and right halves. It can be described on an open
    chain through the following Hamiltonian,
    
    
    H=-∑_m=-L+1^L-1 J_m c^†_m c_m+1 + h.c.,
    
    with hopping amplitudes given by
    
    
    J_m = α^2m+1   if m≠ 0,
    
      1   otherwise,
    
    in terms of an inhomogeneity parameter α∈ (0,1].
    For α=1, the system reduces to the homogeneous case. For small
    α, the ground state (GS) of Hamiltonian (<ref>)
    becomes approximately a  valence bond solid with concentric bonds
    around the center, see Fig. <ref>. This GS violates
    maximally the area-law <cit.>, giving rise to a volumetric
    growth of entanglement.
    
    
    
     §.§ Entanglement and number fluctuations
    
    
    The eigenstates of Hamiltonian (<ref>) are Slater
    determinants. Therefore, the exact entanglement entropy for any block
    can be efficiently evaluated <cit.>. Yet, this technique of
    evaluating the entanglement entropy cannot be extended to a linear
    combination of Slater determinants.
    
    Let us consider an arbitrary linear combination of the ground state
    and first excited states of Hamiltonian (<ref>), both
    within the half-filling sector, i.e. with N/2 particles.
    
    ||⟩Ψ>=α||⟩0>+β||⟩1>.
    
    We are interested in estimates of the entanglement between a region
    A and its complementary, always within state ||⟩Ψ>. The
    procedure is straightforward in the case of a single Slater
    determinants <cit.>. In that case, the computation of the
    full spectrum of the reduced density matrix can be performed in
    O(N^k) operations, with a low value of k. On the other hand, the
    evaluation of the entanglement properties for a generic state takes,
    in principle, O(2^N) operations. Yet, there are some observables
    which act as  entanglement witnesses, i.e.: whose values act as
    an indicator of entanglement. The fluctuation of the number of
    particles is one of these observables <cit.>. 
    
    Let n_A be the number of particles on region A:
    
    
    n_A=∑_i∈ A n_i = ∑_i∈ A c^†_i c_i,
    
    whose fluctuations are given by
    
    σ_N^2 =  n_A^2  - n_A^2.
    
    The expectation value of the number of particles on region A is
    found through the following calculation:
    
    
        n_A   =∑_i∈ A<Ψ|n_i|Ψ> 
        
          =   ∑_i∈ A
            |α|^2 <0|n_i|0> + |β|^2 <1|n_i|1> .
            
            
             +   . 2Re α̅β<0|n_i|1>
        .
    
    The first two terms are easily found, because they refer to a single
    Slater determinant. The third, notwithstanding, must be found using
    the (generalized) Löwdin tricks described:
    
    <0|n_i|1>=<0|c^†_i c_i|1>=0|1 -
    <0|c_ic^†_i|1>.
    
    
    The quadratic term is more involved:
    
    
        n_A^2=    ∑_i,j∈ A
            |α|^2<0|n_in_j|0>
            +|β|^2<1|n_in_j|1> . 
               . + 2 Re α̅β<0|n_in_j|1>
        .
    
    
    Again, the first two terms are straightforward to obtain. If C_A,0 is
    the submatrix of the correlation matrix corresponding to block A on
    state ||⟩0>, then
    
    <0|n_in_j|0>=(C_A,0)^2+(C_A,0(I-C_A,0)).
    
    
    The last term of Eq. (<ref>) is the most involved one, because
    we cannot assume Wick's theorem. We find
    
    
        <0|n_in_j|1>   = <0|(1-c_ic^†_i)(1-c_jc^†_j)|1>= 
           = 0|1 - <0|c_jc^†_j|1> - <0|c_ic^†_i|1> 
           + δ_ij<0|c_ic^†_j|1>
        -<0|c_ic_jc^†_ic^†_j|1>.
    
    
    It has been proved that, for a single Slater determinant
    <cit.>,
    
    
    S_A≥ 4log 2σ_N^2.
    
    Thus, even though Eq. (<ref>) is not proved for a generic
    linear combination of Slater determinants, we will employ σ_N
    as our estimate for an  entanglement witness.
    
    
    
    
     §.§ Numerical experiments
    
    
    For concreteness, let us set α=√(x) and β=√(1-x) for
    x∈ [0,1] in Eq. (<ref>). Thus, our state will be
    given by 
    
    ||⟩Ψ> = √(x)||⟩0> + √(1-x)||⟩1>,
    
    Let us notice that the first excitation is obtained from the ground
    state by performing a parity transformation on the Fermi level
    <cit.>. Notice that, by construction 0|1=0, thus
    forcing us to make all our computations in the  zero overlap case.
    
    In the top panel of Fig. <ref> we show the variance σ_N of the
    number of particles in the left half of the rainbow system with N=8,
    for different values of α. Notice that only in the
    α→ 0^+ limit the variance is the same for x=0 and x=1,
    i.e.: in the GS and the first excited. This variance has been computed
    in two different ways: the dots correspond to the exact calculation,
    with the full Slater determinant, and the continuous line corresponds
    to the computation performed with the generalized Löwdin formulas derived
    in this article. The agreement is complete, and the
    computational time is enormously reduced with our tools.
    
    
    
    In the α→ 0^+ limit, the left-half of the rainbow system
    becomes an infinite temperature mixed state. Thus, the fluctuations in
    the particle number are easy to obtain, following a binomial
    distribution, σ_A = √(N/8), that we can readily check in
    Fig. <ref>. Also, the entanglement entropy will grow up to the
    maximal possible value, (Nlog 2)/2.
    
    The lower panel of Fig. <ref> shows the Von Neumann entropy for
    the same blocks as the top panel. Notice that, in this case, our
    calculations cannot be easily extended. Our generalized Löwdin
    calculation scheme is much more efficient, and can be extended to
    larger system sizes. Fig. <ref> shows the deviation in the
    number of particles of state (<ref>) as a function of x for
    various values of α in a chain with N=40 sites.
    
    
    
     §.§ Generalized Löwdin C++ code
    
    
    We have uploaded our C++ libraries to compute matrix elements of
    Slater determinants using the generalized Löwdin approach to the
    repository  github as  free software<cit.>. The same
    material can also be downloaded from the Supplementary material section of
    the electronic version of the journal <cit.>.
    
    
    
    § CONCLUSIONS AND FURTHER WORK
    
    
    We have extended the seminal results of Löwdin to obtain efficiently
    the matrix elements of an arbitrarily large product of fermionic
    operators between arbitrary Slater determinants. Our results are still
    applicable when the overlap matrix between the orbitals of the Slater
    determinants is singular, i.e. when the corresponding states are
    orthogonal.
    
    Efficient computation of matrix elements in non-orthogonal Slater
    determinants will open a very interesting possibility: the creation of
    Ansätze including Slater determinants obtained from different
    procedures and, therefore, using different orbitals.
    
    As proposals of future work, we would like to remark the extension of
    the previous calculations to the obtention of the full reduced density
    matrix, combining our results with those of <cit.> for the
    reduced density matrix of a block of a single Slater determinant.
    
    
    
    We thank G.F. Bertsch for a careful reading of the manuscript and 
    suggestions. This work has been supported by the Spanish Ministerio de 
    Ciencia, Innovación y Universidades and the European regional 
    development 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).
    
    apsrev4-1
    99
