2021 Generalised Nonorthogonal Matrix Elements: Unifying Wick’s Theorem and the Slater–Condon Rules

Evaluating matrix elements between nonorthogonal Slater determinants represents a key component of several modern electronic structure methods. One approach, the ﬁrst-quantised generalised Slater–Condon rules, is applicable to any pair of determinants with mutually nonorthogonal orbitals, but involves biorthogonalising the occupied orbitals. Alternatively, the nonorthogonal Wick’s theorem allows matrix elements to be evaluated using second-quantisation, but fails if the two determinants have a zero overlap overall. In this contribution, we unify these two approaches by deriving an entirely generalised variant of the nonorthogonal Wick’s theorem that is applicable to pairs of determinants with both nonorthogonal orbitals and a total zero overlap. Our approach is therefore applicable to any pair of nonorthogonal Slater determinants, and can be used to evaluate matrix elements between excited conﬁgurations while only biorthogonalising the occupied orbitals in the reference determinants.

Matrix elements between nonorthogonal Slater determinants represent an essential component of many emerging electronic structure methods. However, evaluating nonorthogonal matrix elements is conceptually and computationally harder then their orthogonal counterparts. While several different approaches have been developed, these are predominantly derived from the first-quantised generalised Slater-Condon rules and usually require biorthogonal occupied orbitals to be computed for each matrix element. For coupling terms between nonorthogonal excited configurations, a second-quantised approach such as the nonorthogonal Wick's theorem is more desirable, but this fails when the two reference determinants have a zero many-body overlap. In this contribution, we derive an entirely generalised extension to the nonorthogonal Wick's theorem that is applicable to all pairs of determinants with nonorthogonal orbitals. Our approach creates a universal methodology for evaluating any nonorthogonal matrix element and allows Wick's theorem and the generalised Slater-Condon rules to be unified for the first time. Furthermore, we present a simple well-defined protocol for deriving arbitrary coupling terms between nonorthogonal excited configurations. In the case of overlap and one-body operators, this protocol recovers efficient formulae with reduced scaling, promising significant computational acceleration for methods that rely on such terms.

I. INTRODUCTION
Matrix elements between nonorthogonal Slater determinants are increasingly common in emerging electronic structure methods. For example, capturing strong correlation using a linear combination of nonorthogonal Slater determinants is a relatively old idea 1-4 that has seen a renaissance in the past decade. [5][6][7][8][9][10][11][12][13][14] . Similarly, nonorthogonal matrix elements arise in projected Hartree-Fock methods 15,16 while the combination of geminal-based nonorthogonal functions is an area of ongoing research. 17 In each method, the nonorthogonality of different determinants can capture strong static correlation effects by breaking and restoring symmetries of the Hamiltonian, 15 or it can provide a quasi-diabatic representation of dominant electronic configurations. 5,9 Alternatively, multiple wave functions built from different orbitals arise in orbital-optimised excited states identified through methods such as ∆SCF, [18][19][20][21] excited-state mean-field theory, 22,23 or the complete active space self-consistent field. 24,25 In these cases, orbital optimisation can significantly improve predictions of charge transfer excitations, but nonorthogonal matrix elements are required for inter-state coupling terms such as oscillator strengths.
A variety of different approaches have been developed for the efficient evaluation of nonorthogonal matrix elements, [26][27][28][29][30] which are predominantly derived from Löwdin's general formula. 31 e most popular framework in quantum chemistry is the generalised Slater-Condon rules, 27,32 where biorthogonal occupied orbitals are constructed 33,34 and a modified form of the Slater-Condon rules 35 is applied depending on the number of zero-overlap orbital pairs in the biorthogonal basis. is approach is applia) Electronic mail: hugh.burton@chem.ox.ac.uk cable to any pair of determinants, but requires the diagonalisation of the occupied orbital overlap matrix each time. In contrast, the development of many-body correlation methods using orthogonal determinants has greatly benefited from the second-quantised Wick's theorem. 36 While a nonorthogonal variant of Wick's theorem exists, it is limited to determinants that have a non-zero many-body overlap and is not applicable if there are zero-overlap orbital pairs in the biorthogonal basis. 37,38 is limitation arises because the ouless theorem, 39 used to relate two nonorthogonal determinants via an exponential transformation, breaks down when the two determinants have a zero many-body overlap. As a result, the nonorthogonal Wick's theorem has seen only limited use in quantum chemistry. 16,40,41 Computationally efficient nonorthogonal matrix elements become increasingly important in methods that use orthogonally-excited configurations from nonorthogonal reference determinants. For example, including post-NOCI dynamic correlation in methods such as perturbative NOCI-MP2 [42][43][44] and NOCI-PT2, 14 or nonorthogonal multireference CI, 3,16,41,45 requires overlap, one-body, or two-body coupling terms between excitations from nonorthogonal determinants. e number of nonorthogonal matrix elements therefore grows rapidly, and repeated biorthogonalisation of the occupied orbitals becomes prohibitively expensive. In principle, the nonorthogonal Wick's theorem could allow these matrix elements to be evaluated using only biorthogonal reference orbitals but, until now, this requires the reference determinants to have a strictly non-zero overlap.
In this contribution, we derive an entirely generalised nonorthogonal form of Wick's theorem that applies to any pair of determinants with nonorthogonal orbitals, even if the overall determinants have a zero overlap.
is new framework, which we call the "Extended Nonorthogonal Wick's eorem", provides the most general approach for deriving matrix elements using second-quantisation, allowing Wick's theorem and the generalised Slater-Condon rules to be unified for the first time.
One particular advantage of our approach is that it allows all matrix elements between excited configurations from a pair of nonorthogonal determinants to be derived using a single well-defined protocol. While some of the resulting expressions have previously been derived for various bespoke applications (see e.g. Refs. 6, 42, and 45), these have o en relied on the properties of matrix determinants to account for orbital excitations. Instead, we present a unifying theory that can recover all of these results and is automatically applicable in cases where the reference determinants have a zero overlap. Furthermore, we show how evaluating intermediates for a given pair of determinants can reduce the scaling of overlap and one-body coupling terms between excited configurations to O(1). ese particular non-orthogonal matrix elements then become almost as straightforward as the orthogonal Slater-Condon rules or Wick's theorem, promising considerable acceleration for methods that rely on such terms.
To derive our generalised nonorthogonal matrix elements, we first define our notation in Section II. In Section III, we extend ouless' theorem 39 to the case where the two determinants have nonorthogonal orbitals and a zero many-body overlap. Section IV combines this extended ouless transformation with Wick's theorem to create a generalised protocol for evaluating nonorthogonal matrix elements using secondquantisation. We then illustrate the application of this approach by re-deriving the generalised Slater-Condon rules for one-and two-body operators in Section V. Finally, in Section VI, we extend our framework to the matrix elements between excited configurations and show how O(1) scaling can be achieved for overlap and one-body operators.

II. NOTATION
We will consider matrix elements between the two determinants | w Φ and | x Φ . Each determinant is constructed from a bespoke set of molecular orbitals (MOs), represented in terms of the atomic spin-orbital basis functions |χ µ as Here, we employ the nonorthogonal tensor notation of Head-Gordon et al. 46 to explicitly keep track of any required overlap matrices. Occupied MOs are indexed as (i, j, k, etc), virtual MOs as (a, b, c, etc), and any general MO as (p, q, r, etc). We emphasise that the MOs are orthogonal within each Slater determinant, but are nonorthogonal between the different determinants.
In second-quantisation, the N -electron determinant | w Φ is defined as where |− is the physical vacuum and the molecular orbital creation operators w b † i satisfy the standard fermionic anticommutation rules. 47 Using the expansion Eq. (1), the MO creation and annihilation operators can be represented in terms of the covariant atomic spin-orbital creation a † µ and annihilation a µ operators as e covariant atomic spin-orbital operators have only one non-zero anticommutator 47 where g µν = χ µ |χ ν defines the corresponding covariant metric tensor (overlap matrix). 46 roughout this paper, we will need to express the atomic spin-orbital creation and annihilation operators in terms of the MO creation and annihilation operators using the inverse of Eq. (3) To avoid the introduction of overlap matrices throughout our expressions, we will o en use the contravariant atomic spin-orbital operators (a µ ) † and a µ defined as where, g µν is the contravariant metric tensor corresponding to the inverse covariant overlap matrix, 46 i.e.
If the AO basis is overcomplete, this contravariant metric tensor becomes the pseudo-inverse of the covariant overlap matrix. Note that the anticommutator of these contravariant atomic spin-orbital operators is

III. EXTENDED THOULESS TRANSFORMATION
A. Conventional Thouless Transformation e conventional form of ouless' theorem allows two nonorthogonal determinants to be related by an exponential operator of single excitations as 39 To derive the single excitation operator Z, the occupied orbitals can be be transformed to a biorthogonal basis using Löwdin pairing 33,34 such that e virtual-occupied and virtual-virtual blocks of the biorthogonalised overlap matrix become µν ( xC * ) ·µ a· g µν wCν· ·p = xwS ap , (11) and the transformed molecular orbital creation and annihilation operators are given as e single excitation operator in Eq. (9) is then defined as with the xw Z ai matrix elements given by A brief derivation of this result can be found in Appendix A. Unfortunately, this exponential representation relies on the strict nonorthogonality of the two determinants x Φ| w Φ = 0; in other words, it is not applicable to a pair of determinants that are orthogonal but contain mutually nonorthogonal orbitals. Our first step is therefore a generalisation of the ouless transformation to the case where x Φ| w Φ = 0.

B. Introducing Zero-Overlap Orbitals
We begin in the biorthogonal basis identified through Löwdin pairing, with orbital coefficients satisfying Eq. (10). For a general pair of nonorthogonal orbitals, it is possible for orbital pairs to have a zero-overlap in the biorthogonal basis, where xwS i = 0. Taking the case with m zero overlaps between orbitals k 1 , · · · , k m , we construct "reduced" determinants by removing the electrons in these zero-overlap orbitals to give ese reduced determinants are strictly nonorthogonal with the non-zero reduced overlap defined as erefore, the ouless transformation can now be applied to these reduced determinants to give Here, we have introduced the reduced single excitation oper-atorZ that only contains excitations from occupied orbitals with a non-zero overlap as e full N -electron determinants are then related through second-quantisation as Equation (19) can be further simplified by exploiting the commutativity relation [Z, xb k ] = 0 to shi the exp Z operator to the far right-hand side, giving We can then introduce single-electron excitation operatorŝ z k for the zero-overlap orbitals aŝ where we have exploited the biorthogonality and zerooverlap of the occupied orbitals such that xwS ik = 0 for all i. e commutativity of these single excitation operators with the xb k annihilation operators leads to the simplified relationship Introducing the relationshipẑ k = exp(ẑ k ) − 1 allows Eq. (22) to be expanded as Expanding the product of (exp(ẑ k ) − 1) terms then leads to a sum of exponential transformations where every combination of the zero-overlap single excitation operatorsẑ k is either included or excluded with an appropriate phase factor, giving Here, we have introduced the compound index C n to denote a particular combination of n out of m zero-overlap orbitals, while the superscript notationZ Cn indicates which particular zero-overlap excitations are included in the corresponding operator, i.e.Z We refer to this transformation in its various forms (22)- (24) as the "Extended ouless Transformation". To explicitly illustrate its application, the case of two zero-overlaps in orbital pairs k 1 and k 2 leads to whereZ (k1,k2) =Z +ẑ k1 +ẑ k2 andZ (k1) =Z +ẑ k1 .

A. Conventional Wick's Theorem
Efficiently deriving matrix elements using the conventional Wick's theorem requires the introduction of contractions, defined for two creation or annihilation operators x b p and x b q as where { x b p x b q } represents a normal-ordered operator string with respect to the reference Fermi vacuum x Φ| · · · | x Φ . 36 e only non-zero contractions between creation and annihilation operators with respect to this symmetric Fermi vacuum are rough Wick's theorem, the Fermi vacuum expectation of an operator product is given by the sum over all fully contracted products of operators, e.g.

B. Zero-Overlap Transformed Operators
Using the extended ouless transformation, we can now extend the nonorthogonal Wick's theorem 37,38,48 to derive matrix elements between any pair of determinants with mutually nonorthogonal orbitals. In what follows, we will consider the contravariant atomic spin-orbital operators (see Section II) to avoid large numbers of overlap matrices in our expressions. e matrix elements for general operators expressed in the atomic spin-orbital basis requires the evaluation of terms containing a string of creation and annihilation operators, such as x Φ|(a µ ) † (a ν ) † · · · a σ a τ | w Φ . Applying the extended ouless transformation leads to the linear combination To evaluate each constituent matrix element for the combinations C n , we follow the approach described in Refs. 38 and 40 and introduce a similarity-transformed set of spin-orbital creation and annihilation operators as ese operators clearly depend on the particular combination C n of included zero-overlap single excitation operators. Expanding the similarity transformation as and similarly for (a µ ) † , leads to the explicit forms An explicit derivation of these relationships can be found in Appendix B. Exploiting the relationship and the resolution of the identity then allows the constituent matrix elements within Eq. (30) to be expressed as C. The Fundamental Contraction e extended ouless transformation essentially converts the nonorthogonal matrix element with an asymmetric Fermi vacuum x Φ| · · · | w Φ to a transformed matrix element with respect to symmetric Fermi vacuum x Φ| · · · | x Φ . Since the transformed operators d[C n ] µ and (d[C n ] µ ) † are expressed purely in terms of the xb p creation and xb † p annihilation operators, with respect to the x Φ| · · · | x Φ vacuum, their nonzero contractions with respect to x Φ| · · · | x Φ can be derived by combining Eqs. (28) and (33) to give where we have introduced the general notation ese expressions closely resemble co-density matrices, 5 and their derivation can be found in Appendix C. e overall matrix element x Φ|(a µ ) † (a ν ) † · · · a σ a τ | w Φ requires the derivation of contractions between the a µ and (a µ ) † atomic spin-orbital operators rather than the transformed d[C n ] µ and (d[C n ] µ ) † operators. To show how a general string of operators can be evaluated using the contractions defined Eqs (37a) and (37b), we first demonstrate the derivation of the one-body co-density matrix element Assuming that some form of Wick's theorem can be derived, we expect this matrix element to be represented by the single contraction Taking the most general case with m zero-overlap orbitals, the extended ouless transformation leads to Reversing the order of summation over k and C n yields e first term in square brackets C n 1 can be recognised as the total number of ways to pick n orbitals from the m zerooverlap orbitals, given by m n . Similarly, the second term in square brackets C n / ∋ k 1 is the total number of ways to pick n orbitals from the m − 1 zero-overlap orbitals that remain when orbital k is excluded, given by m−1 n . ese combinatorial expansions allow Eq. (42) to be expressed as Here, we note that there are no ways to exclude an orbital k when all zero-orbital overlaps are included in the complete combination C m . Exploiting the binomial expansion then leads to the closed-form expression e reduced overlap xwS will be a prefactor for every matrix element between these nonorthogonal determinants. e remaining terms can then be used to define "fundamental contractions" for second-quantisation operators with respect to the asymmetric Fermi vacuum x Φ| · · · | w Φ . e form of these contractions depends on the number of zero-overlap orbitals m and we define the first fundamental contraction as Similarly, the second fundamental contraction can be identified as Crucially, we emphasise that these fundamental contractions are defined with respect to the asymmetric Fermi vacuum x Φ| · · · | w Φ .

D. Combining Several Contractions
Next, we show how the fundamental contractions can be combined to derive matrix elements for longer products of creation and annihilation operators. As an example, consider the two-body reduced co-density matrix element, defined as Applying Wick's theorem, this matrix element should be given by the sum of the two contractions Note that the second term in this expression carries a phase of −1 from the intersection of the contraction lines, representing the fundamental parity of fermionic operators. 36 Taking the first contraction as an example, we apply the extended ouless transformation and the transformed contractions defined in Eqs. (37a) and (37b) to give Once again, the reduced overlap xwS appears as an overall prefactor. e order of summation over the k 1 , k 2 indices and C n can then be swapped to give 1 .
e third term in square brackets C n / ∋k1,k2 1 is simply the number of ways to pick n orbitals from the m − 2 zero-overlap orbitals that remain when orbitals k 1 and k 2 are removed, given by m−2 n . Applying the binomial expansion Eq. (44) in an analogous way to the single contraction leads to the closed form A similar expression can be derived for the second contraction as which, analogously with the orthogonal case, will carry a −1 phase factor from the intersection of the contraction lines. Combining these two equations, with their associated phase factors, yields the full expression for the two-body reduced codensity matrix elements with m zero-overlap orbitals as

E. General Rules for Constructing Matrix Elements
To simplify the derivation of even longer operator strings, we note that the two-body matrix elements in Eq. (52) can be factorised into the product of two fundamental contractions with individual m 1 and m 2 values under the constraint m 1 +m 2 = m, giving is factorisation can be extended for a general product of contractions with each m i restricted to the values 0 or 1 for the overall product to be non-zero. We can therefore define an intuitive approach for extending Wick's theorem to generalised nonorthogonal matrix elements: 1. Construct all fully contracted combinations of the operator product with the associated phase factors.
2. For each term, sum every possible way to distribute m zeros among the contractions such that i m i = m.
3. For every set of {m i } in each term, construct the relevant contribution as a product of fundamental contractions depending on whether each contraction has m i = 0 or 1.

4.
Multiply the final combined expression with the reduced overlap xwS . ese rules and contractions therefore allow any matrix element to be evaluated with respect to the asymmetric Fermi vacuum x Φ| · · · | w Φ for nonorthogonal orbitals, regardless of whether the many-body determinants are orthogonal or not. As a result, our formulation is the most flexible form of Wick's theorem, and reduces to the previous nonorthogonal 37,38,48 or orthogonal 36 variants under suitable restrictions on the MO coefficients. We refer to this approach as the "Extended Nonorthogonal Wick's eorem". In the following Sections, we will show how these steps can be applied to recover the generalised Slater-Condon rules, 32 and to derive matrix elements between excited configurations with respect to the nonorthogonal reference determinants.
V. GENERALISED SLATER-CONDON RULES e generalised Slater-Condon rules provide a firstquantised approach for evaluating matrix elements between nonorthogonal determinants. A detailed description of these rules, their derivation, and their application, can be found in Ref. 32. However, to the best of our knowledge, the generalised Slater-Condon rules have never previously been derived though a fully second-quantised framework. In this Section, we will show how these rules can be recovered using the extended nonorthogonal Wick's theorem.

A. One-Body Operators
Consider first the one-body operator with the corresponding matrix element Applying the extended nonorthogonal Wick's theorem yields only one non-zero contraction to give Substituting the fundamental contraction [Eq. (46)] and considering the possible values of m immediately yields the onebody generalised Slater-Condon rules 32 in the formed pre-sented in Ref. 5 as Here, we note that xw M νµ = xw W νµ when there are no zero-overlap orbitals, while xw P νµ = xw P νµ k for one-zero overlap in orbital k. e matrices xw M and xw P k correspond respectively to the weighted and unweighted co-density matrices discussed in Ref. 5.

B. Two-Body Operators
Next, consider a general two-body operatorv defined in second-quantisation aŝ with the two-electron integrals in the spin-orbital basis defined as v µντ σ = χ µ χ ν |v|χ τ χ σ .
(61) e corresponding matrix element is given by Recognising the x Φ|(a µ ) † (a ν ) † a τ a σ | w Φ term as the twobody reduced co-density matrix derived in Eq. (54), we immediately recover Here, note that the m = 1 terms each contain two terms with the single zero-overlap in either the first or second contraction respectively. Exploiting the symmetry v µντ σ = v νµστ and the identity xw P σµ xw P τ ν − xw P τ µ xw P σν = k1 =k2 then allows the two-body generalised Slater-Condon rules 32 to be recovered as Note that, for the m = 1 case with one zero-overlap in orbital k, we have exploited the identities xw P τ µ k xx P σν k = xw P σµ k xx P µν k and xw P τ µ k xw P σν k = xw P σµ k xw P µν k to introduce the simplification xw M → xw W .

VI. MATRIX ELEMENTS FOR EXCITED CONFIGURATIONS
While re-deriving the generalised Slater-Condon rules provides an important verification of the extended nonorthogonal Wick's theorem, it does not provide any computational advantage over the original framework. On the contrary, the primary focus of our new framework involves deriving matrix elements between excited configurations from a pair of nonorthogonal determinants e.g., x Φ ab... ij... | · · · | w Φ cd... kl... . Terms of this form arise in perturbative corrections to NOCI, 14,42-44 the NOCI-CIS approach for core excitations, 8,49 and the evaluation of S 2 coupling terms in NOCI expansions. 6 Furthermore, these nonorthogonal matrix elements will be required to evaluate inter-state coupling elements between orbital-optimised excited-state wave functions identified using excited-state mean-field theory 22 or state-specific complete active space SCF. 24 Until now, evaluating these matrix elements has required the direct application of the generalised Slater-Condon rules to each pair of excitations.
is approach leads to significant computational costs associated with the biorthogonalisation of the excited determinants, which scales as O N 3 each time. In this Section, we show how the extended nonorthogonal Wick's theorem allows these matrix elements to be evaluated using only biorthogonalised reference determinants. When large numbers of coupling elements are required, avoiding the biorthogonalisation of these excited configurations significantly reduces the overall computational cost. Furthermore, the cost of certain nonorthogonal matrix elements (including all one-body operators) becomes independent of the number of electrons or basis functions if additional intermediate matrix elements are computed for each pair of reference determinants.
In this Section, we describe how the extended nonorthogonal Wick's theorem can be applied to matrix elements between excited configurations. First, we discuss how the fundamental contractions described in Section IV C can be modified to an MO-based form that makes excited configurations easier to handle. We then illustrate the derivation of certain overlap, one-body, and two-body matrix elements between nonorthogonal excited configurations. e resulting expressions are entirely generalised for any pair of nonorthogonal reference determinants and can significantly reduce the computational scaling compared to a naïve application of the generalised Slater-Condon rules.

A. Asymmetric Representation
Evaluating matrix elements between two excited determinants will require the evaluation of the asymmetric contractions x b † p w b q with respect to the asymmetric Fermi vacuum x Φ| · · · | w Φ . Expanding the molecular orbital creation and annihilation operators using Eq. (3) yields Introducing the fundamental contractions for (a µ ) † a ν and (a µ ) † a ν defined in Eqs. (46) and (47) respectively then reduces these expressions to different forms depending on the number of zero-overlaps m k associated with the contraction, Here, we have defined the "screened" overlap terms and with xw M and xw P defined in Eqs. (38b) and (38d) respectively, and the overlap element Although the indexing notation used in Eq. (68a) may seem counterintuitive, we find that it helps to keep track of the bra and ket orbital coefficients in the screened overlap terms. Crucially, the orbital coefficients used to evaluate these contractions do not need to be the same as those used to evaluate the xw M and xw P matrices. is feature is particularly advantageous as the excited configurations can be represented in terms of the original orbital basis while the xw M and xw P matrices are evaluated in the biorthogonal basis. As a result, only the reference determinants need to be biorthogonalised, and the remaining matrix elements are evaluated in terms of these screened overlap terms. Furthermore, the screened overlap elements are themselves one-body matrix elements that can be computed once for a given pair of determinants and stored, before being combined to evaluate more complicated matrix elements. With N ref determinants, the total cost of computing these intermediates therefore scales as O N 2 ref n 3 . To take full advantage of these asymmetric contractions, the one-and two-body operators can also be represented in terms of one set of molecular orbitals aŝ where we have defined the transformed matrix elements and Evaluating nonorthogonal matrix elements through the extended nonorthogonal Wick's theorem then proceeds using similar steps outlined in Section IV E: 1. Assemble all fully contracted combinations of the asymmetric operator and excitation operator product and compute the corresponding phase factors.
2. For each term, sum every possible way to distribute m zeros among the contractions such that i m i = m.
3. For every set of {m i } in each term, construct the relevant contribution as a product of fundamental contractions defined in Eqs. (67a) and (67b) depending on whether each contraction has m i = 0 or 1.

4.
Multiply the combined expression with the reduced overlap xwS of the reference determinants. Note that the number of zero-overlap orbitals in these expressions corresponds to the biorthogonalised reference determinants, not the excited configurations.
In the following Sections, we will illustrate this process through a series of typical nonorthogonal matrix elements. As we shall see below, evaluating the sum of every combination of zero-overlap indices assigned to each contraction quickly leads to complicated equations. However, we can derive the general structure of a matrix element for the case with no zero-overlap orbitals m = 0, and then recover the forms for different values of m by distributing the zerooverlap indices over each fundamental contraction. We will therefore focus on the parent equation with m = 0, which we refer to as the "canonical form" and denote using the notation · · · 0 . Furthermore, any matrix element with more zero-overlap orbitals than the total number of contractions must be strictly zero.

B. Overlap Terms
First, we consider the overlap element between excited determinants. e simplest overlap matrix element involves only a single excitation x Φ| w Φ a i . For m zero-overlap orbitals in the reference determinants, this single excitation matrix element can be identified using one contraction as Here, the canonical form is simply and there is only one way to assign one zero-overlap orbital to the single contraction. Note that the reduced overlap between the reference determinants remains a prefactor for the overall matrix elements. Next, the overlap of two single excitations x Φ a i | w Φ b j can be evaluated as (76) e canonical form for this element is the m = 0 term e m = 1 term is recovered by taking the sum of the two different ways to distribute one zero-overlap orbital to the two contractions, while there is only one way to distribute two zero-overlap orbitals for the m = 2 term.
As a third example, consider the double excitation overlap x Φ| w Φ ab ij , which can be evaluated as In this case, the minus sign arises from the intrinsic −1 phase that results from the intersection of the contraction lines in from which the terms for m = 1 and m = 2 can be derived.
Notably, this form of x Φ| w Φ ab ij has previously be derived for m = 0 in Refs. 42 and 45, but our derivation generalises for any number of zero-overlap orbitals.

C. One-Body Operators
We now consider matrix elements for one-body operators of the form given in Eq. (71a). To further simplify the subsequent expressions, we can introduce intermediate matrices that account for partial contraction with the one-body operator. In particular, we introduce the partially contracted inter-mediate terms Note that the individual contraction in the F 0 term, indicated by the X matrix, may correspond to a zero-overlap orbital pair, as indicated by the notationF 0 . e yw [XF X] rs (and similar) terms correspond to two contractions, and can thus be assigned two zero overlap orbitals. e different possibilities of assigning these zero-overlap contractions is denoted as yw [XF X] rs or yw [XFX] rs , and similarly for terms involving the Y contraction. Crucially, these intermediate values also correspond to orbital pairs, and can be precomputed once for each pair of reference determinants, leading to a oneoff computational cost that scales as O N 2 ref n 3 . As a result, the summation over the p, q indices is avoided for the subsequent evaluation of matrix elements between excited determinants, and the computational scaling of these one-body terms is the same as the overlap matrix elements.
To illustrate the application of this approach, we take the simplest one-body matrix element x Φ|f | w Φ a i . For m zero-overlap orbitals in the reference determinants, this matrix element can be expanded as Using the contractions defined in Eqs. (67a) and (67b), and the intermediate terms defined in Eq. (80), the canonical form for m = 0 is given as If the F 0 , ww X ia , and ww [XF Y ] ia are precomputed and stored once for the reference determinants, the subsequent cost of evaluating x Φ|f | w Φ a i is independent of the number of electrons or basis functions as scales as O (1). From this canonical form, expressions for m = 1 and 2 (or higher) can be identified as Next, consider the coupling of two single excitations x Φ a i |f | w Φ b j , corresponding to the fully contracted terms from which expressions for m > 0 can be obtained as described above.

D. Two-body Operators
Finally, we consider matrix elements for two-body operators, with the general form given in Eq. (71b). Again, we can define intermediate matrices in the orbital basis that can be pre-computed for a given pair of nonorthogonal determinants. First, we introduce analogues of the Coulomb and exchange matrices, respectively defined as Partially contracted intermediate matrices can then be defined as e.g., with the constant value Like the one-body intermediate matrices, when the zerooverlap orbitals are distributed among the contractions, these zeros may be independently assigned to any of the X, Y , or (J − K) terms, as each of these contain one contraction. For example, with m = 1, one must consider the sum of three terms given as e.g., On the contrary, the V 0 constant term corresponds to two contractions and can be assigned two zero-overlap orbitals.
Noting the symmetry x v pqrs = x v qpsr , we denote these possibilities as Here, the factor of two in Eq. (93a) arises because the zerooverlap orbital can be assigned to either the first or second contraction to give the same result.
To illustrate the application of this approach for excited configurations, we first consider the two-body matrix element Combining each contraction, and exploiting the intermediates in Eqs. (89) and (90) yields the canonical form (m = 0) for this matrix element as Distributing the zero-overlap orbitals among each contraction then leads to explicit expressions for all m values as Next, we consider the two-body coupling of two singly excited determinants as With a total of 24 possible ways to fully contract the corresponding matrix element, we maintain brevity by directly presenting the canonical form for m = 0 as From here, expressions for the cases with m > 0 can be recovered by distributing the zero-overlap orbitals over the contractions associated with the V , X, (J − K), or Y matrices. While partially contracted intermediate expressions could also be evaluated to avoid the nested summation on the last line in Eq. (98), this will generally incur an unacceptably large storage cost. Finally, we consider the two-body coupling of a reference determinant and a double excitation Again, there are a total of 24 possible ways to fully contract the corresponding matrix element, so we advance directly to the canonical form for m = 0, given as We note that this formula has previously been identified in Ref. 45 for m = 0, but our approach now allows this result to be extended for m > 0.

E. Illustration of Scaling
rough the use of intermediate matrices in Sections VI A-VI C, the cost of evaluating overlap and one-body matrix elements between excited configurations can be made independent of the number of electrons N or basis functions n. In both cases, the computational cost of evaluating the intermediate matrices scales as O N 2 ref n 3 , which is determined by the cost of evaluating the screened overlap terms in Eqs. (68a) and (68b). e subsequent cost of evaluating matrix elements between excited configurations then scales as O(1). In contrast, the cost of applying the generalised Slater-Condon rules for overlap or one-body operators is dominated by the evaluation of the occupied orbital overlap matrix, giving a scaling of O n 2 N for every excited coupling term.
For two-body operators, a similar application of two-body intermediates requires transformations of the two-electron integrals into an asymmetric MO representation for each pair of determinants, and for every possible combination of four contractions in the four possible forms given in Eq. (67a) or (67b). ese two-body intermediates therefore carry a stor-age overhead that scales as O 4 4 N 2 ref n 4 , which will generally be unacceptably large. We imagine there may be certain applications that require only a subset of these two-electron intermediates, which could then be stored to achieve the subsequent O(1) scaling. In general, however, the n 4 scaling for two-body matrix elements cannot overcome, although the overall cost is still reduced by avoiding biorthogonalisation of the occupied orbitals in the excited configurations.
To explicitly illustrate the computational speed-up that might be expected using the extended nonorthogonal Wick's theorem, consider the evaluation of matrix elements between two configuration interaction singles (CIS) wave functions built from different reference Slater determinants In particular, consider the evaluation of the overlap a one-body operatorf , e.g. the transition dipole moment, and a two-body operatorv such as the two-electron repulsion e relevant overlap, one-body, and two-body nonorthogonal matrix elements are given in Eqs (76), (84), and (98) respectively. In each case, only one set of intermediates need to be evaluated using the biorthogonalised reference determinants, followed by a total of N 2 n 2 nonorthogonal matrix elements between singly excited configurations (assuming n ≫ N ). Here, we only consider one pairing of the reference determinants.
e computational scaling associated with first evaluating the intermediate terms from the reference determinants, and then evaluating all the nonorthogonal coupling terms between the excited configurations is summarised in Table I. For the overlap and one-body terms, the computational scaling is dominated by the evaluation of the intermediates, which scales as O n 3 . In comparison, the total cost of biorthogonalising the occupied orbitals for every pair of excited configurations scales as O N 3 n 4 . For two-body operators, the cost for both approaches is dominated by the contraction with the two-body integrals for each pair of excited configurations.
is leads to a total scaling O N 2 n 6 for both cases, but the extended nonorthogonal Wick's theorem may still be faster as the prefactor is reduced for each element.

VII. CONCLUDING REMARKS
Intuitive derivations and efficient implementations of nonorthogonal matrix elements are becoming increasingly important for the development of nonorthogonal configuration interaction methods and inter-state coupling terms for state-specific excited state wave functions. However, until now, the evaluation of these matrix elements has relied on the generalised Slater-Condon rules, while the second-quantised nonorthogonal Wick's theorem fails when there are any zerooverlap orbital pairs between the reference determinants. In this work, we have extended the nonorthogonal Wick's theorem to the case where two determinants have nonorthogonal orbitals, but have a zero many-electron overlap. is new theory, which we call the Extended Nonorthogonal Wick's theorem, provides the most generalised framework for evaluating matrix elements between two nonorthogonal determinants using second quantisation.
Among the primary advantages of our new approach is the ease of deriving and evaluating matrix elements between excited configurations from nonorthogonal reference determinants. To illustrate this feature, we have derived a series of overlap, one-body, and two-body coupling terms between nonorthogonal excited configurations. For overlap terms and one-body operators, these excited nonorthogonal matrix elements can be expressed in terms of one-body intermediate terms that can be precomputed and stored for a given pair of reference determinants. As a result, the subsequent cost of evaluating the coupling of excited configurations scales as O (1), which is the same as the conventional Slater-Condon rules or Wick's theorem for orthogonal reference determinants. However, in the current form of the theory, the cost of evaluating two-body coupling terms scales as O n 4 , which is the same as applying the generalised Slater-Condon rules for each pair of nonorthogonal excited configurations. From a computational perspective, the extended nonorthogonal Wick's theorem will therefore provide the greatest acceleration in tasks that require at most one-body coupling terms, such as the evaluation of transition dipole moments for orbital-optimised excited states, [22][23][24][25] or the coupling elements of a one-body reference Hamiltonian in nonorthogonal perturbation theories. 14,[42][43][44] Looking forwards, we hope that this work will encourage and accelerate future developments in nonorthogonal electronic structure theory by creating a unifying theory for deriving challenging nonorthogonal matrix elements. As part of this vision, we are working on an open-source C++ library for evaluating typical nonorthogonal matrix elements, and we intend to report on this project soon.