Zombie states are a recently introduced formalism to describe coupled coherent fermionic states that address the fermionic sign problem in a computationally tractable manner. Previously, it has been shown that Zombie states with fractional occupations of spin orbitals obeyed the correct fermionic creation and annihilation algebra and presented results for real-time evolution [D. V. Shalashilin, J. Chem. Phys. **148**, 194109 (2018)]. In this work, we extend and build on this formalism by developing efficient algorithms for evaluating the Hamiltonian and other operators between Zombie states and address their normalization. We also show how imaginary time propagation can be used to find the ground state of a system. We also present a biasing method, for setting up a basis set of random Zombie states, that allows much smaller basis sizes to be used while still accurately describing the electronic structure Hamiltonian and its ground state and describe a technique of wave function “cleaning” that removes the contributions of configurations with the wrong number of electrons, improving the accuracy further. We also show how low-lying excited states can be calculated efficiently using a Gram–Schmidt orthogonalization procedure. The proposed algorithm of imaginary time propagation on biased random grids of Zombie states may present an alternative to the existing quantum Monte Carlo methods.

## I. INTRODUCTION

The wavefunction of a multi-electron system is usually accurately described by a linear combination of Slater determinants.^{1} Each Slater determinant represents a specific configuration—how the available electrons are distributed across spin orbitals. However, in a system with many electrons and many orbitals, the number of possible Slater determinants increases rapidly. Even with sensible approximations and truncation, accurately describing the electronic structure of a system requires a large number of Slater determinants and with that comes computational cost. Introducing an active space of a few most important orbitals and few active electrons is the main method of treating this problem. Active space removes abruptly all configurations outside of it, but even though the contribution of individual configurations outside of the active space can be small, their total contribution can be significant. A smooth transition between active and virtual orbitals in active space methods can potentially be beneficial. Such a transition can be perhaps achieved with quantum Monte Carlo (MC) methods. Diffusion MC^{2} and Green’s function MC^{3} propagate walkers in continuum space, which removes the need for a one-electron basis set, but they both fall foul of the fermionic sign problem caused by the antisymmetry of the electronic wavefunction. It has been possible to solve this problem using a fixed node approximation,^{4} but this has had limited further work. A more robust approach is to introduce a Monte Carlo way of selecting the configurations via a random walk in the manifold of Slater Determinants. Full Configuration Interaction Quantum Monte Carlo (FCIQMC), developed by Alavi and co-workers, works by using a long-time propagation in imaginary time and random walkers to stochastically describe the full configuration interaction wavefunction (FCI).^{5,6} The method does not converge to the bosonic solution^{5} and obtains the fermionic ground state without a fixed node approximation. FCIQMC has been successful at obtaining FCI results for large systems that were not previously possible, including the neutral and cationic elements from Li to Mg.^{5–9} A similar method is Monte Carlo Configuration Interaction (MCCI), developed by Greer,^{10,11} which originally was applied to the single point energy of water^{10} and later used to find the dissociation energy of water and Hartree–Fock (HF).^{12} Like FCIQMC, MCCI uses a Monte Carlo procedure to build a compact wavefunction containing the important configurations intended to give accuracy close to FCI results. Coe *et al.* continued to develop MCCI and showed that it can generate potential energy curves, approaching levels of chemical accuracy, for a range of small molecules, including N_{2} and CH_{4}.^{12} MCCI was then extended to include natural orbitals and second-order perturbation theory, which, particularity, at longer bond lengths was shown to improve the accuracy and convergence time of calculated energies.^{13} MCCI has also been used to simulate multipole moments achieving good accuracy when compared to FCI values.^{14} MCCI also showed good accuracy when compared to FCIQMC when calculating ionization energies while only requiring a relatively small number of basis functions when compared to the FCI space.^{14} The method has also been extended to include state-averaging to allow computation of excited states, which was shown to be effective for both H_{3} and LiF.^{15}

However, when used as part of a molecular dynamics simulation, the computational cost of electronic structure calculations, especially Monte Carlo methods, is still significant. Hence, finding methods with acceptable cost of electronic structure calculations, which at the same time would not neglect configurations that normally would be outside of the active space, is of great benefit. In a recent paper,^{16} we introduced the idea of Zombie States (ZSs), which potentially can combine the ideas of randomness and active space. Zombie states have the possibility of an electron in every orbital in a superposition of “dead” and “alive” states^{16} so that the occupation of an orbital is fractional. This allows multiple Slater determinant configurations to be described by a single ZS. However, ZSs may not have a well-defined number of electrons such that there is a nonzero probability of the state containing a wrong number of electrons. But, a small superposition of ZSs might be able to describe a wave function with the required number of electrons accurately, and the hope is that a smaller basis set size could be used, reducing the computational cost while maintaining acceptable accuracy.

ZSs can be viewed as fermionic Coherent States (CSs) of a two-level system. Previously, coherent states of the harmonic oscillator have been used to describe bosonic systems in second quantization with the help of Herman-Kluck^{17} and coupled coherent states propagation methods^{18} and also with generalized coherent states.^{19} Second quantization Hamiltonians look like those of coupled oscillators with the difference that the oscillators represent the amplitudes and populations of the orbitals occupied by Bosons. Following the idea of these methods, it seems reasonable to try to treat fermions on a similar level. The complication here is that the elements of Grassmann algebra enter into the standard definition of a fermionic coherent state:^{20} $\eta j=a(0j+\eta j1j)$, where $1j$ and $0j$ are the *j*th orbital in its occupied and unoccupied vacuum states, respectively; *a* is a normalization factor; and *η*_{j} is an element of Grassmann algebra. The necessity of Grassmann algebra makes computation difficult but is required to maintain the correct permutation antisymmetry of the multi-electron fermionic coherent states as well as the anticommutation of creation and annihilation operators. Zombie states were shown^{16} to be capable of describing fermionic coherent states while removing the need for Grassmann algebra and the use of Wick’s theorem, which are normally required to evaluate matrix elements between fermionic coherent states.^{21}

Earlier work^{16} gave a mathematical treatment of the ZSs’ second quantization for fermions, and creation and annihilation operators were defined with the use of a simple sign change rule, which replaces Grassmann algebra. These results were verified by the reproduction of full configuration interaction^{1} electronic energies for Li_{2} and LiH via diagonalization of the electronic structure Hamiltonian in the complete basis of randomly selected ZSs. It was also shown that Zombie states can be used for real time propagation of the electronic wave function. A Fourier transform of real-time Zombie state evolution also has reproduced the exact electronic energy levels in the above-mentioned molecules. Both methods are, however, not very efficient and could be used only as a demonstration that Zombie states’ mathematics was correct. It is therefore necessary to develop better ZS-based methods for finding low-lying state electronic energies in molecules, which are most important in chemistry. In this vein, it is necessary to improve upon the naïve algorithms initially presented to increase efficiency; computation of the Hamiltonian in the original article^{16} had $O(M5)$ scaling, which can be reduced. This also follows for other key operators. Bosonic coherent states can be defined as an eigenstate of the annihilation operator, whereas the only eigenstate of the annihilation operator for a fermionic coherent state is the vacuum. Thus, we also seek to provide an alternative formalism for Zombie states in the language of second quantization that allows Zombie states to be created from a vacuum state. In this paper, we will present the theoretical advancement in the formalism of Zombie states, deriving a stronger normalization condition. Efficient algorithms for key operators are then given utilizing scaled algorithms and sensible manipulations of the program to greatly reduce the computation time of the Hamiltonian as well as other properties, such as spin and the number of electrons.

Furthermore, imaginary time propagation is validated as an effective method for finding the ground-state energies of states of real fermionic systems, capable of reproducing full CI (in a truncated basis),^{1} in this case for Li_{2}. We also show how a biasing method can be used to create a small basis of Zombie states, which with imaginary time evolution can find low-lying state energies. To further improve the accuracy, we show how the procedure called “cleaning” removes the contribution of unwanted numbers of electrons and improves the accuracy of electronic energy. We also make an attempt to create a method for finding excited states. It has previously been shown that a Gram–Schmidt procedure can be applied to orthogonalize higher energy states against lower lying ones while adding minimal computational cost.^{22} We apply a similar Gram–Schmidt procedure to imaginary time propagation of ZSs too, showing that low-lying excited states can be found in this manner.

This article is structured as follows: a brief recap of the original Zombie state formulation is given followed by the algebraic developments; the algorithmic improvements for the two-electron Hamiltonian are detailed; results for imaginary time evolution of Li_{2} for various basis sets and finally results for imaginary time evolution for excited states of Li_{2} using Gram–Schmidt orthogonalization are given.

## II. THEORY

### A. Original formulation

Zombie states^{16} are coherent antisymmetrized superpositions of “dead” and “alive” electronic states. Considering the single *m*th spin orbital, a coherent state is

where $1m$ corresponds to there being an electron in spin orbital *m* and $0m$ corresponds to the *m*th spin orbital being empty, consistent with the conventional electronic structure notation.

We can generalize this to a coherent state, which is a Slater Determinant of one-electron Zombie states,

and which can be summarized by 2*M* coefficients,

The notation used here is synonymous with that in Ref. 16, where in $amjj$, *m*_{j} = 0 refers to the *m*th spin orbital being “dead” (unoccupied) and *m*_{j} = 1 refers to the *m*th spin orbital being “alive” (occupied) and *j* refers to the spin orbital number. Note that the Zombie state contains 2*M* complex coefficients, where *M* is the total number of spin orbitals.

The overlap of two states is

Calculation of this is $O(M)$, where there are *M* spin orbitals, whether or not they are occupied.

The action of the creation and annihilation operators on a ZS is given by

The operators $b\u0302m\u2020$ and $b\u0302m$ not only act on the *m*th orbital but also change sign of the amplitudes *a*_{1n} of alive coefficients for all orbitals with *n* < *m*. If *a*_{1m} and *a*_{0n} are given by ones and zeros so that Zombie states represent standard Slater determinants, it is easy to see that the sign changing rule is equivalent to the so called Wigner–Jordan rule. As previously shown, all anticommutation relations are as expected.^{16} Most straightforwardly, but not very efficiently, matrix elements $\zeta (a)H\u0302\zeta (b)$ of the second quantized electronic structure Hamiltonian can be computed by sequential application of the creation and annihilation operators,^{16}

to $\zeta (b)$ and overlapping the result with $\zeta (a)$,

where $\zeta klmn(b)=b\u0302k\u2020b\u0302l\u2020b\u0302mb\u0302n\zeta (b)$ and $\zeta mn(b)=bm\u2020b\u0302n\zeta mn(b)$; the overlaps are calculated using Eq. (4).

An electronic wavefunction can be represented as a superposition of *K* basis Zombie states,

and the matrix elements described above now allow usage of the ansatz equation (8) for propagation or finding quantum states and their energies. Individual Zombie states generally are not restricted to a particular number of electrons, but for a wavefunction such as Eq. (8), a sufficiently large *K* number of Zombie states with coefficients *d*_{k} usually can be chosen such that contributions from unwanted numbers of electrons cancel out. This is very different from standard electronic structure methods, where the problem is to include as many configurations with the right number of electrons as possible.

A standard Hartree–Fock configuration $\phi me(j)$ that corresponds to *m*_{e} electrons can be written as a Zombie state with “binary” amplitudes of dead and alive states and *m*_{e} ones in the upper row,

and it can be treated just like any other ZS.

In Subsections II B and II C, algebraic and algorithmic improvements are presented, which allow much more efficient calculation of matrix elements of operators between ZSs and faster propagation of the wave function equation (8). These sections are technical: an alternative derivation of Zombie states from the vacuum state is presented as is the full method we used to reduce Hamiltonian scaling. These are important to the ZS methodology as a whole, but a reader less interested in technical aspects of the method can go directly to the Sec. II D, where it is shown that the right electronic energies can be recovered from imaginary time propagation of the wave function equation (8) with an unrestricted number of electrons.

### B. Algebraic developments

The previous article^{16} gave the action of creation and annihilation operators for Zombie states and how the Hamiltonian could be computed for them. Here, we extend the algebraic formalism by showing how a Zombie state can be created from a vacuum state. We show how creation and annihilation operators act in this formalism, as well as computation of the overlap of two Zombie states. Finally, we derive a stronger condition for normalization of a Zombie state.

From the standard electronic structure theory,^{1} a given electronic occupancy can be written as an antisymmetrized Slater determinant as

where {*m*_{i}} = 0 if the *i*th orbital is empty and 1 if occupied, as detailed in Ref. 23 or, equivalently, in the form of Eq. (9). This is equivalent to stating in second quantization notation,

where $$ is the vacuum state, ensuring here and in what follows that the creation operators are applied in the reverse order they appear in the Slater determinant. We could also trivially rewrite Eq. (11) as

We now consider a generalized form of Eq. (12),

where *a*_{0j} and *a*_{1j} are complex scalar coefficients. Clearly, if *a*_{1j} = *m*_{j} and *a*_{0j} = (1 − *m*_{j}), then the same Slater determinant electron state is produced as in Eq. (12).

By inference from Eq. (13), we can define a Zombie operator

whose adjoint is

and commutators and anti-commutators

such that Eq. (13) can be written more concisely as

Note that $z\u0302j$ is a function of two complex numbers, so we can write $z\u0302j\u2261z\u0302j(a0j,a1j)$. This shorthand will become useful later.

The idea of defining new operators from linear combinations of creation and annihilation operators is, of course, nothing new, the most famous example perhaps being Majorana Fermions, which are their own antiparticle, defined as

However, Zombie operators are not their own antiparticles. We now consider the action of creation and annihilation operators on Zombie operators.

#### 1. Creation and annihilation operators

The action of creation and annihilation operators on a single Zombie operator is given, for *j* ≠ *k*, by

and if *j* = *k*,

where we use the standard relations^{1} ${b\u0302j,b\u0302k}=0$, ${b\u0302j\u2020,b\u0302k\u2020}=0$, and ${b\u0302j,b\u0302k\u2020}=\delta jk$.

Using these results, we can compute the action of the creation and annihilation operators on a Zombie state, which is defined by the product of Zombie operators as in Eq. (17),

where the $z\u0302j(a0j,\u2212a1j)b\u0302j$ term vanishes because $b\u0302j$ tries to destroy an electron that is not there.

It also follows that

However, Eqs. (22c) and (23b) themselves describe Zombie states with changed coefficients to those of $\zeta $. The Zombie states they refer to agree exactly with the sign-changing rules presented in the original Zombie states paper, though through more convoluted algebra.

These results can equivalently be stated in terms of the ${anjj}$ coefficients. As $\zeta $ is a function of the Zombie coefficients ${anjj}=a$, we can write $\zeta \u2261\zeta (a)$. Then, let $\zeta (j)(a(j))=b\u0302j\zeta (a)$. We then see that

If $\zeta (j)(a(j))=b\u0302j\u2020\zeta (a)$,

again in exact agreement with the earlier results.

#### 2. Overlap of two Zombie states

We can use this formalism to compute the overlap of two Zombie states, $\zeta (a)\zeta (b)$, where both the bra and ket are formed from the action of *M* Zombie operators as in Eq. (17). We consider a computation recursively, noting that the vacuum state is normalized, namely, $\u3009=1$, and defining

such that

We, therefore, have

The second term in Eq. (28b) is of the form $\zeta j+1(a)b\u0302j\u2020\zeta j+1(b)$, and since neither $\zeta j+1(a)$ nor $\zeta j+1(b)$ contain any contributions from electron *j*, $b\u0302j\u2020$ acting to the left will try to destroy an electron that is not there, giving zero. Similar arguments show why the third term in Eq. (28b) vanishes and why the fourth term survives.

#### 3. Normalization

Clearly, if we wish a Zombie state (the product of *M* Zombie operators) to be normalized, then the RHS of Eq. (29) must equal unity when $\zeta (a)=\zeta (b)$. However, we may want to consider a more general case, where we may not necessarily have wavefunctions formed from a fixed *M* number of Zombie operators but a variable number of Zombie operators.

Algebraically, if we require

this will ensure that $\zeta (a)\zeta (a)\u2261\zeta 1(a)\zeta 1(a)=1$, but what if we further require

for all *j*, 1 ≤ *j* ≤ *M*.

This problem can also be solved recursively, since trivially $\zeta M+1(a)\zeta M+1(a)\u2261\u3009=1$. For $\zeta M(a)\zeta M(a)=1$, from Eq. (28), it is required that $|a0M(a)|2+|a1M(a)|2=1$. However, from Eq. (26), $\zeta j(a)\zeta j(a)=(|a0j(a)|2+|a1j(a)|2)\zeta j+1(a)\zeta j+1(a)$, and so for Eq. (31) to be true for all *j*, we require

Note that Eq. (32) is a stronger criterion than Eq. (30), i.e., satisfying Eq. (32) will ensure that Eq. (30) is satisfied, but not the other way around.

Equation (32) might correspond to normalization of a “one dimensional” ZS. This normalization is equivalent to the normalization used when introducing quantum superposition sampling to the multiconfigurational Ehrenfest method presented in Ref. 24.

In the event that the Zombie state is real (useful for imaginary time propagation to find stationary states), then a simple choice to satisfy Eq. (32) is

where *θ*_{j} is a real number such that 0 ≤ *θ*_{j} < 2*π*. This means that a real, normalized Zombie state can be completely described by *M* real {*θ*_{j}} values. This can be generalized to complex states^{24} with

where *ϕ*_{j} is drawn from the interval 0 ≤ *ϕ*_{j} < 2*π*.

### C. Algorithmic developments

Here, we consider how to efficiently compute the action of the Hamiltonian and other operators on a Zombie state. We show how the electronic Hamiltonian matrix elements with Zombie states $\zeta (a)H\u0302\zeta (b)$ can be calculated more efficiently. First, we show how a more efficient algorithm with a smaller prefactor can calculate matrix elements faster, and second, we will present an algorithm with $O(M4)$ scaling as opposed to the $O(M5)$ scaling of the naïve algorithm used in the initial paper^{16} described above.

#### 1. Reduced prefactor Hamiltonian algorithms

The bottleneck in the computation of FCI states and energies, in general, is evaluation of the two-electron Hamiltonian

where $\zeta ijlk(b)=bi\u2020b\u0302j\u2020b\u0302lb\u0302k\zeta (b)$. This can be solved with iterative applications of the creation and annihilation operators, meaning that the original algorithm, therefore, requires *M*^{4} terms to be evaluated to solve Eq. (36). Each term has four creation/annihilation operators [each of which is $\u223cO(M)$ due to the sign changing rule] and one overlap computation, which is $O(M)$ as

Consequently, naïve evaluation of Eq. (6) requires *M*^{4} creation/annihilation operations and is $O(M5)$ overall.

Many of the two-electron integrals are zero by spin symmetry, i.e., $ijkl\u3009=0$ unless *σ*(*i*) = *σ*(*k*) and *σ*(*j*) = *σ*(*l*), where *σ*(*i*) is the spin of spin orbital *i*. One method for accelerating the naïve algorithm is by evaluating only $\zeta ijlk(b)=bi\u2020b\u0302j\u2020b\u0302lb\u0302k\zeta (b)$ if $ijkl\u3009\u22600$ on spin symmetry grounds. When applied to the restricted Hartree–Fock (RHF) calculation, this simple code modification gave nearly a four times speed up when calculating a matrix element between two ten-orbital ZSs. This is roughly to be expected since three quarters of the two-electron integrals are zero from spin symmetry.

Furthermore, it is also possible to loop only over combinations of {*ijkl*} already known to be nonzero by altering the loops over orbitals such that only nonzero (by spin) $\zeta ijlk(b)$ are considered in the first place. This produces a comparable, albeit slightly larger time improvement, ignoring zero two-electron integrals.

In an attempt, to accelerate the code further, we note that computing the action of *M*^{4} creation and annihilation operators is expensive and look to evaluate fewer of them. Re-examining Eqs. (35) and (36), we see

where

We can then compute ${\zeta ij(a)}\u2200i,j$ and ${\zeta lk(b)}\u2200l,k$, which requires *M*^{2} creation and annihilation operator applications, while overall the Hamiltonian evaluation is $O(M5)$. Having precomputed ${\zeta ij(a)}$ and ${\zeta lk(b)}$, we can then evaluate Eq. (38).

From the definition of the overlap integral in Eq. (4), we see that if *both* the dead and alive Zombie amplitudes are zero for a given state, i.e., for the *m*th spin orbital,

then the overlap of this state with any other state will be zero. Furthermore, since the Zombie creation and annihilation operators only move the dead and alive amplitudes within a spin orbital, and not between spin orbitals, $\zeta (c)$ will continue to have zero overlap with any other state *even after application of creation and annihilation operators*. Consequently, if a state such as $\zeta (c)$ is itself generated by application of creation and annihilation operators, it will not contribute to the two-electron Hamiltonian and so does not need to be considered in further calculations.

We know that^{16}

Consequently, if *a*_{1} = 0, application of $b\u0302$ returns a vanishing state that will have no overlap with any other. This is a rephrasing, in Zombie terms, of the well-known result from electronic structure theory that trying to annihilate an electron, which is not there will return zero. Applied to calculation of the two-electron Hamiltonian, let us consider $\zeta ij(a)=b\u0302jb\u0302i\zeta (a)$. If $a1i(a)=0$, then $b\u0302i\zeta (a)$ returns a vanishing state and will contribute nothing to the two-electron calculation, even after application of $b\u0302j$. If $b\u0302i\zeta (a)\u2255\zeta i(a)$ is nonzero, then if $a1j(a)=0$, the Zombie state vanishes. Consequently, when looping over orbital indices in Eq. (38), if the application of a given annihilation operator would return a vanishing state (which can be determined from $a1i(a)=0$) and similarly for *j*, *l*, *k*, then the all further calculation for that creation/annihilation operator can be discarded.

In addition, considering the overlap calculation equation (4) and rewriting it as

we see that if *any* of the $a0m(a)*a0m(b)+a1m(a)*a1m(b)$ terms are zero, then Ω_{ab} = 0. This means that calculation can be terminated as soon as one of the terms is zero and the overlap returned as zero, without having to compute any remaining terms. These improvements gave a similar time improvement to previous changes when working with ten orbitals; however, the time improvement was approximately half as good when calculating a matrix element with 50 orbitals. Full details of the algorithmic code racing can be found in Sec. 1.1 of the supplementary material. All of these algorithms still have (approximate) fifth order scaling, i.e., $O(M5)$, although with a much lower prefactor than earlier versions.

#### 2. Lower-scaling Hamiltonian algorithm

It is, however, possible to reduce the scaling of the two-electron Hamiltonian from $O(M5)$ to $O(M4)$. We first note that evaluation of the two-electron Hamiltonian could be written as

We have previously shown that $\zeta ij(a)$ can be calculated in $O(M3)$ and $\zeta k(b)$ in $O(M2)$, and these can be calculated separately. Using previous methods, even with the calculation improvements, to evaluate $\zeta ij(a)b\u0302l\zeta k(b)\u2200i,j,k,l$ requires $O(M3)$ for looping over *i*, *j*, *k*, then an $O(M)$ for summing over *l*, and then another $O(M)$ for evaluating the overlap, making $O(M5)$ overall. We, therefore, consider how to evaluate $\zeta ij(a)b\u0302l\zeta k(b)$ for a specific set of *i*, *j*, *k*, but over all values of *l*, in $O(M)$ operations, such that the overall algorithm would be $O(M4)$.

We start by simplifying our notation to $\zeta (a)b\u0302l\zeta (b)$, which yields

Note the minus sign in the first product caused by the action of the annihilation operator. We then define

These can all trivially be calculated in $O(M)$ steps. We also define

which can be calculated recursively in $O(M)$ steps. Then,

which can also be found in $O(M)$ steps.

The reduced scaling improves calculation of two-electron matrix elements over 20 times for ten orbitals calculation and over 100 times with 50 orbital calculations compared to the naïve calculation. Further details are given in Table 4 in the supplementary material. In the case of 50 orbitals, the scaled algorithm does incur a larger peak memory usage, requiring ∼414% more memory for the calculation. However, the fastest algorithm still only needs 2 MB of memory and so the speed increase is worth the relatively larger memory requirement.

We have also developed an algorithm to calculate the number operator, which scales linearly with the number of orbitals in the system. Using a similar method, we also show how to reduce the scaling of the algorithms used to calculate *S*_{z}, $S\u0302z2$ and total spin. Details of these two algorithms are presented in Appendixes A and B.

### D. Imaginary time evolution, wave function cleaning, biasing, and excited state calculation

Imaginary time propagation, i.e., solution of the Schrödinger equation in imaginary time,

is commonly used to find the ground state of a quantum mechanical system.^{5,6,22} Many methods of quantum Monte Carlo exploit it in electronic structure theory. For example, a random walk in the Fock space of Slater determinants has been shown to be a beautiful technique to solve Eq. (49) and to avoid the wave function node problem, which otherwise makes quantum Monte Carlo for fermions difficult. In this section, we will show that a small randomly selected set of Zombie states can provide an accurate description of the ground state of a molecule.

#### 1. Imaginary time algorithm

Consider a wavefunction $\Psi $ that can be expanded in the basis of *K* Zombie states $\zeta (k)$, where *K* ≤ 2^{M}, and which are normalized but not necessarily orthogonal,

Then, applying the identity matrix from Ref. 25 for a nonorthogonal basis set $\u2211l,k\zeta l)\Omega lk\u22121\zeta (k)$ and using the reasoning from Ref. 26, we yield

which can be rearranged in matrix notation as

Note how the inverse overlap matrix with the elements $\Omega lk=\zeta (l)\zeta (k)$ is required and how the algebra developed for bosonic coherent states is being applied to fermionic ZSs.^{26}

#### 2. Cleaning

The propagation of Eq. (52) relaxes the wave function to the lowest energy state of the system.^{27} As occupations of ZS spin orbitals are fractional Zombie states, their linear combinations are not restricted to a particular number of electrons. However, they can be projected onto a Fock Space with a given number of electrons. The identity operator can be written as a sum of *m*_{e} electron identities $I\u0302me$,

where the identity covering the Fock Space with *m*_{e} electrons is

where the sum in *j* is over all Fock configurations $\phi me(j)$, with *m*_{e} being number of electrons, which can be viewed as Zombie states with “binary” amplitudes 1 and 0 and *m*_{e} being unit amplitudes of alive electrons on *m*_{e} occupied spin orbitals.

Then, the superposition of Zombie states in Eq. (8) or Eq. (50) can be projected onto Fock Space of *m*_{e} electrons as

where

The energy of the *m*_{e} contribution and its norm can then be calculated as

and

respectively. The sum of energies of all *m*_{e} contributions is equal to the energy of the whole Zombie wave function Eq. (8) and all *m*_{e} norms add up to 1, i.e.,

and

For a completely converged Zombie wave function, the energy and norm of all *m*_{e} other than the right number of electrons *n*_{e} will be zero, but if convergence is incomplete, the energy of the converged *m*_{e} wave function $\Psi c$ can be estimated as

#### 3. Biased bases

In this section, we will consider various ways of how a basis set of ZSs can be chosen to solve Eq. (52). The best choice (in theory) would be a complete set of 2^{M} Slater determinants $\phi me(j)$ covering all possible electron occupancies from 0 to *M* electrons in *M* spin orbitals. Diagonalization of the Hamiltonian matrix in this basis gives all possible energy levels as it contains all possible combinations, from 0 to *M* electrons, placed in *M* spin orbitals. A basis of randomly selected 2^{M} ZSs $\zeta (k)$ is also complete and is capable of yielding correct quantum energies and wave functions. However, for efficient calculation, the basis set must be as small as possible. We employ a biased basis that sits between the two extremes of complete order and randomness. In complete active space self-consistent field (CASSCF) calculations,^{30} electrons and orbitals are split into three groups: inactive, active, and virtual. From this, we have designed a biasing method to set up a ZS basis. Inactive, or core electrons, are low lying and are always (at least in a CASSCF calculation) to be occupied; the active orbitals can either be occupied or unoccupied and virtual orbitals are always empty. We take this idea of splitting electrons into different groups and then exploit the fractional occupation of a Zombie state. This allows orbitals to be set as completely “alive” (or “dead”) while also biasing orbitals to be more or less “dead” or “alive.” It should be noted that though we are using the ideas of splitting electrons into different types from CASSCF, this does not make ZS an SCF method any more than using HF energies as a starting point. Amplitudes for Zombie states can be randomly chosen from normal distributions centered around either completely “alive” or completely “dead” values and the distributions’ width set according to its activity level. To satisfy the normalization condition given by Eq. (32), in this paper, only one Zombie amplitude (dead or alive) is generated using the normal distribution and the other is set using Eq. (33). Hence, the low-lying core electrons have “alive” amplitudes set to one, although a very narrow Gaussian could also be employed, which corresponds to a normal distribution width of zero giving a *δ* function with *a*_{1j} = 1 and *a*_{0j} = 0. The active electron normal distributions can be centered to favor either completely alive or dead amplitudes; the point where centering changes from “alive” to “dead” naturally changes with each system. Thinner distributions are used for orbitals at either end of the active electron/orbital group and wider distributions are used for orbitals where the chance, or not, of occupation is relatively equal. Virtual orbitals would have “dead” coefficients set to 1. However, again, narrow distribution near completely “dead” amplitudes can be employed.

#### 4. Excited states

As previously stated, the aim of developing Zombie states is so it can be implemented in no-adiabatic molecular dynamics simulations and so it would be useful to compute low-lying excited states with the low computational expense that Zombie states could provide. It has previously been shown that a Gram–Schmidt orthogonalization procedure can be combined with full configuration interaction quantum Monte Carlo^{22} and was found to add little extra computational cost while allowing multiple low-energy states to be studied. We will apply the same process here combining imaginary time propagation and Gram–Schmidt. Generally, Gram–Schmidt orthogonalization, for any set of vectors, is defined by

with the projection operator

In this implementation, the *k* lowest states of interest are arbitrarily set, creating **d**_{k} vectors that are not necessarily orthogonal. The vectors are then each orthogonalized,

where $dk\u2032$ denotes an orthogonalized vector, and they are then normalized. Each state $dk\u2032$ is then put into the differential equation equivalent to Eq. (52),

A single time step is taken, and a new set of $dk(\beta +\Delta \beta )=dk\u2032(\beta )+d\u0307k(\beta )\Delta \beta $ is then found. The process then repeats: vectors are orthogonalized using the Gram–Schmidt process and normalized and a time step is taken.

## III. RESULTS

Here, we will demonstrate that imaginary time propagation is an effective and efficient method for finding ground state energies using a Slater determinant basis as a point of reference. Earlier work^{16} used Li_{2} to verify the theoretical basis of Zombie states, and so here we too will use Li_{2} as an example system. PyScf^{31} (version 1.7.6) using a 6-31G** basis and the *D*_{∞h} symmetry group has been used to calculate one and two-electron integrals. These are then converted from spatial to spin orbitals, in this case, five spatial molecular orbitals (and therefore ten spin orbitals). A Hamiltonian for the system is computed, using the Hamiltonian detailed in Sec. II C 2, with the specified ZSs. The full code used can be found in Ref. 32 All calculations are carried out in atomic units so energies are in Hartrees.

### A. Ground-state imaginary time propagation

The first trivial case considered is a complete basis of Slater determinant $\phi me(j)$ Zombie states. Li_{2} consists of ten spin orbitals, which gives a full basis a size of 2^{10} = 1024 basis functions. In Fig. 1, the six-electron restricted Hartree–Fock (RHF) determinant is used as a starting point and imaginary time evolution gives the neutral ground state energy of Li_{2}. A full sized basis of 2^{10} = 1024 randomly generated basis functions $\zeta (k)$ was also tested. Random Zombie states are generated by randomly generating *θ*, between 0 and 2*π*, for each orbital, and the dead and alive coefficients are then found by calculating *a*_{1} = cos(*θ*) and *a*_{0} = sin(*θ*), respectively. The random basis is made up of Zombie states consisting of superposition of “dead” and “alive” electrons compared to the Slater determinant basis, which is entirely binary. The initial vector is a superposition of random Zombie states, which is chosen to be equal to the RHF determinant, and that propagation of this using imaginary time evolution (solid line in Fig. 1) matches the result obtained by using the Slater Determinant basis. The final energy of each state obtained through imaginary time propagation is verified by comparison to the eigenvalues found by diagonalizing the complete Slater determinant Hamiltonian.

When a reduced basis of 200 random Zombie states is used as seen in Fig. 2, the basis set is no longer capable of reproducing the ground state energy. The reduced basis is created by selecting 200 basis functions $\zeta (k)$ at random from the complete random basis previously used. The energy obtained via imaginary time propagation, in this smaller basis, no longer matches the value found by diagonalizing the Slater determinant basis. However, the cleaned energy provides a much better estimate of the ground state energy. For a poorly selected basis, the distribution of energy and norm over possible electronic numbers *m*_{e} is shown in Fig. 3, which demonstrates that only a fraction (≈17.5%) of the final wave function belongs to the right Fock space on *n*_{e} = 6 electrons. The highest fraction of the norm (≈40%) belongs to the 7*e* Fock space. The cleaned energy *E*(6*e*)/*N*(6*e*) = −14.598 267 is somewhat above the FCI energy of the 6*e* ground state (−14.871 914), but the cleaned energy of the 7*e* Fock Space *E*(7*e*)/*N*(7*e*) = −14.750 162 is approaching the 7*e* FCI (−14.858 062) energy. Thus, the randomly selected basis has, by chance, been set up in favor of the anion.

The electronic structure of the Li_{2} is not completely random and unknown, so it is sensible to exploit this. A biased random basis set of the Zombie states was generated using the biasing regime previously described. Equation (34) was used to ensure that the “dead” and “alive” coefficients were normalized, values *θ*_{j} were randomly generated using a normal distribution centered at $12\pi $ for electrons “alive” in the restricted Hartree–Fock Slater determinant [*j* = 1…6] and centered around 0 for the “dead” electrons. The width, *σ*, of the normal distributions used to bias each electron is summarized in Table I, and these values are used to create the biased bases in all subsequent results. There are six electrons in the neutral Li_{2} molecule and four core inactive electrons that are set as always “alive.” Figure 4 schematically shows how this biasing method is applied: the first two spatial orbitals (first four spin orbitals) are set to always be occupied with the final three spatial orbitals (last six spin orbitals) having fractional occupancy.

. | Figure . | |
---|---|---|

j-th electron
. | μ\2π
. | σ\2π
. |

1, 2, 3, 4 | 0.25 | 0 |

5, 6 | 0.25 | 0.175 |

7, 8 | 0 | 0.351 |

9, 10 | 0 | 0.120 |

. | Figure . | |
---|---|---|

j-th electron
. | μ\2π
. | σ\2π
. |

1, 2, 3, 4 | 0.25 | 0 |

5, 6 | 0.25 | 0.175 |

7, 8 | 0 | 0.351 |

9, 10 | 0 | 0.120 |

Figure 5 also demonstrates that ZS propagation can yield the ground state of the $Li2\u2212$ anion. A 64 biased Zombie state basis was generated, but the first basis function was set as the seven electron open-shell restricted Hartree–Fock Zombie state and the remaining 63 Zombie states generated using the same biasing method as before. Similarly, the first Zombie state can be set to the six-electron restricted Hartree–Fock Zombie state. It can be seen in Fig. 5 that depending on the choice of the initial condition, the same ZS propagation can lead the ground state of neutral molecule and its anion. Figure 6 shows the result of the propagation using the basis of 10, 30, and 50 random basis ZS functions together with their improved energies obtained by cleaning. The first basis function in each basis set is the six electron restricted Hartree–Fock determinant to ensure each evolution starts at the same energy. As expected, the energies for these smaller biased bases are not as accurate as obtained using the “complete” 64 basis function Hamiltonian. However, the biasing method is considerably more accurate when a directly comparing random and biased bases of the same size.

Figure 7 shows the distribution of energy and norm over the Fock spaces with *m*_{e} electrons for the case of 30 BFs. Almost all norm is in the 6*e* Fock space, but even in this case, cleaning improves the energy and brings it closer to the full CI result. Figure 5 shows that a basis of 64 randomly selected basis functions is complete and yields the exact result even without cleaning.

### B. Low-lying excited states

Gram–Schmidt orthogonalization was then applied to the complete random basis and is shown in Fig. 8. Four separate states were computed, all of which were initially set to arbitrary vectors. The Gram–Schmidt process was applied and then each state was propagated in imaginary time with Gram–Schmidt orthogonalization applied after every time step. It can be seen in Fig. 8 that the first state is the neutral ground state, the second and third states are the degenerate anion ground states (degenerate as *M*_{s} can be ±1/2), and the fourth is the first neutral excited state for Li_{2}. All of these values match the corresponding eigenvalues found by diagonalizing the complete Slater determinant basis.

The same Gram–Schmidt process was then applied to four states using a 64 Zombie state biased Hamiltonian. The imaginary time propagation of these four states can be seen in Fig. 9. As with complete random basis propagating in imaginary time with Gram–Schmidt produces the lowest neutral state, the degenerate anion ground states, and the first excited energy, which are also shown on both figures as dashed lines.

## IV. CONCLUSIONS

The findings of this paper are broad but substantially advance Zombie states as a potential method for simulating electronic structure and dynamics. While previously shown to be numerically consistent, the alternative formalism of Zombie states, presented here, builds from the vacuum state using creation and annihilation operators to create Zombie states. This forms a clearer link between Zombie states and more traditional electronic structure theory. This new formulation also leads to a stricter normalization condition, which not only demonstrated how ideas from bosonic systems can be translated to fermionic systems but was key in the implementation of the biasing method used. A key facet of a practicable modern methods is efficient algorithms for the calculation of the Hamiltonian and other operators. In this paper, we have presented algorithms that greatly improve on the naïve algorithms for the number operator, spin operators ($S\u0302z$ and $S\u03022$), and the two-electron Hamiltonian. A scaled algorithm for the two-electron Hamiltonian is presented and reduces computational cost significantly. Computationally inexpensive calculation of important system properties is obviously very important for any subsequent work. However, this also sets up an easily replicable method for reducing the scaling of any algorithms that may be developed in subsequent work. All of the algorithms follow a similar method using the action of creation and annihilation operators to obtain recursion relations, which lead to lower-scaling algorithms.

We have also shown that using imaginary time propagation on a system of Zombie states, it is possible to find the lowest-lying state of that system. A complete basis of Slater determinants and a basis of complete random Zombie states were shown to give the exact same ground state energies for Li_{2}. Imaginary time propagation removes the need for long real-time evolution and Fourier transforms. Reducing the size of the random basis resulted in the ground state energy obtained no longer matching the benchmark set by the complete Slater determinant basis. However, this problem is overcome by the use of the biasing method, where we use the ideas of CASSCF^{30} (or CASCI) by splitting orbitals into inactive, active, and virtual types and then biasing individual orbitals in the active space. The lowest four spin orbitals were set as core orbitals and so always “alive,” leaving two electrons to fill the remaining six spin orbitals, such that there were 2^{6} = 64 Slater Determinants in a complete basis set for this active space and any more basis functions would lead to linear dependencies and a singularities in the overlap matrix. However, with 64 biased Zombie state basis functions, it was possible to produce, with imaginary time evolution, energy levels that matched the 16 times larger Slater determinant basis. We also showed that with a single Zombie state changed to a specified restricted Hartree–Fock determinant, it was possible to find a specific energy level with imaginary time propagation. It should be noted that the result of converged imaginary time propagation does not give the correct number of electrons or spin. However, applying cleaning, by construction, does recover this information. Reducing the number of biased Zombie states gives less accurate energy levels; however, when compared to imaginary time propagation for random bases of the same size, the biased bases are considerably better. In some respects, the accuracy of calculations with 64 basis functions is not a great surprise from a chemistry viewpoint, since the lowest four electrons are in 1s orbitals and, therefore, tightly bound to the lithium nuclei. It would be energetically very costly to excite these electrons (compared to the 2s electrons), and they are therefore expected to remain core 1s in any low-lying states. Clearly, the use of Zombie states is not giving new scientific insight into the electronic structure of Li_{2}, but we have demonstrated that ZS can reproduce results with a reduced basis in an efficient manner.

It is worth noting that 64 basis functions is less than the 210 configurations possible with six electrons in ten spin orbitals but far more than the 15 possible configurations possible when the first four electrons are fixed, and the remaining two can be arranged across six spin orbitals. Thus, introduction of Zombie states with non-integer occupations of orbitals for a small system has not decreased the required number of basis functions but, in fact, increased them. This is not surprising as the basis of Zombie States describes multi-electronic states with all possible numbers of electrons from no electrons at all to electrons occupying all spin orbitals. For systems with small active space, such as Li_{2}, the price of not sticking to the right number of electrons outweighs the benefit of bringing in other excited electronic configurations, which may be disregarded by active space methods. However, as can be seen from the Fig. 6, a small ZS basis accounts for some correlation energy. For larger systems, where the active space is large, this may be an affordable way of taking electron correlation into account. The exponential scaling of required configurations is a common problem in electronic structure calculations, but ZS should reduce this, with adequate biasing, with each ZS containing many configurations. As ZSs are made up of superpositions of many electronic configurations, we expect that a manageable number of them will be required to describe larger systems. Therefore, future work will focus on applying Zombie states to much larger systems that traditionally require extremely large basis sets, such as N_{2},^{6} to obtain FCI values while using significantly smaller ZS basis sets. We envisage that this work will require the application of sampling techniques previously developed for bosonic systems,^{33} which when applied to Zombie states should aid in converging results and address the possibility of choosing a linear combination of Zombie states such that the overall wavefunction is an eigenstate of the number operator.

## SUPPLEMENTARY MATERIAL

See the supplementary material for details of the code racing for all of the algorithms presented. There is also detail of the code used to generate the electron integrals using PySCF.

## ACKNOWLEDGMENTS

The authors would like to thank George Booth for discussions on PySCF and Dmitry Makhov for help with Molpro. The support of the EPSRC via Grant No. EP/P021123/1 is acknowledged. T.J.H.H. also acknowledges a Royal Society University Research Fellowship (Grant No. URF\R1\201502).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

O.A.B. and T.J.H.H. contributed equally to this work.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: NUMBER OPERATOR

Although the observable properties of a Zombie state can be formally computed exactly by iterative application of creation and annihilation operators followed by computation of the overlap of the resultant states,^{16} this usually is not an efficient calculation. Here, we present algorithms for the commonly used properties of electronic states, such as the number of electrons, and their spin properties.

If we consider the action of the number operator,

where $b\u0302m$ is the annihilation operator and $b\u0302m\u2020$ is the creation operator for Zombie state *m*, we see

such that $n\u0302m$ effectively “deletes” the coefficient $a0m(b)$. Note that the sign switching cancels itself out. It, therefore, appears that acting with the number operator $N\u0302$ is $O(M)$ as the summation is over *M* terms.

A comparison of Eqs. (4) and (A1) suggests that computing the number of electrons represented by any Zombie state, i.e., $\zeta (a)N\u0302\zeta (b)$, would be $O(M2)$ as computing the overlap of two Zombie states is $O(M)$.

However, it is possible to construct an algorithm, which is only $O(M)$ by careful consideration of the summation and the product by adapting an algorithm from Ref. 34. Mathematically, calculation of the action of the number operator is

We can define

and use Eq. (46c) to define

This can be used along with Eq. (47b) to find the following recursion relations:

Computation of Eqs. (A4) and (46c) is clearly $O(M)$ and using the recursion relations and so is computation of {*g*_{l}} and {*h*_{m}}. Inserting the aforementioned equations into Eq. (A3),

which is a summation that can also be computed in $O(M)$ steps. When applied, this new algorithm is significantly faster than the previous method; increasing the number of orbitals by a factor of 10, the new algorithm becomes roughly another 10 times faster than the old one as to be expected from the scaling arguments above. Full details of this code racing can be found in Table 5 of the supplementary material.

Since Zombie states are not usually eigenstates of the number operator, unlike most Slater determinants, the uncertainty in the in number of electrons can be defined as a standard deviation,

For $N\u03022$,

Computing $n\u0302k\zeta (b)$ is trivial (simply delete the “dead” coefficient for orbital *k*). We can therefore define $\zeta k(b)=n\u0302k\zeta (b)$ and compute

using the known algorithm for computing the number operator, so this is $O(M2)$ instead of $O(M3)$.

The general idea presented trivially extends to the “Ghost” operator,

which counts how many of the Zombie states are “dead” (unoccupied), as opposed to the number operator that counts how many of them are alive.

### APPENDIX B: SPIN OPERATORS

#### 1. $S\u0302z$ operator

Orbitals are defined with the assumption that any used in the Zombie state calculation are from a *restricted* calculation such that for each orbital with spin $\alpha $ (spin up, *m*_{s} = +1/2), there exists an orbital with the same spatial wavefunction but spin component $\beta $ (spin down, *m*_{s} = −1/2). For a system with *M* spin orbitals, there will consequently be *K* spatial orbitals where $K=M/2,K\u2208N$. We further define that all spin orbitals with up spin have an odd index *m* = 1, 3, …, *M* − 1, and all spin orbitals with down spin have an even index *m* = 2, 4, …, *M* such that the *k*th spatial orbital has $\alpha $ spin orbital *m* = 2*k* − 1 and $\beta $ spin orbital *m* = 2*k*, *k* = 1, 2, …, *K*.

Using this numbering convention, in second quantization notation,

The optimized number operator is then adapted introducing a simple sign change rule to calculate *S*_{z} in $O(M)$ steps. A comparison of the naïve algorithm, using creation and annihilation operators, and the optimized algorithm, co-opting the optimized number operator equations, can be found in Table 5 of the supplementary material.

#### 2. Faster $S\u0302z2$ computation

Prima facie, computation of $S\u0302z2$ would be $O(M3)$ using Eq. (B1),

so there is $O(M)$ for the summation over *k*, another over *j*, and another to evaluate overlap.

However, extending the adapted number operator algorithm used for $S\u0302z$ in $O(M)$ calculation to $S\u0302z2$ is straightforward,

but $n\u0302m\zeta (b)$ gives another Zombie state that we call $\zeta (b,m)$ such that

The bra-ket can be evaluated in $O(M)$ and, therefore, the action of the $S\u0302z2$ operator in $O(M2)$.

#### 3. Total spin

Usually, in electronic structure calculation, one also wants to know the total spin,^{1}

We can apply the faster $S\u0302z$ and $S\u0302z2$ algorithms as detailed above. $S\u0302+$ and $S\u0302\u2212$ are raising and lowering operators,

where *K* is the number of *spatial orbitals*, viz., *K* = *M*/2. We are numbering the orbitals 1, 2, …, *K* such that the *k*th spatial orbital has an alpha spin orbital number 2*k* − 1 and a beta spin orbital number 2*k*.

We then observe the effect of $s\u0302+,k$ and $s\u0302\u2212,k$ on the wavefunction (where *m* = 2*k*),

and

In both cases, the sign changes cancel each other out exactly. This means that no sign change is required for evaluating $S\u0302+$ or $S\u0302\u2212$, saving computational cost. However, there are clearly reductions to be made in the scaling that should reduce the time needed to compute $\zeta (a)s\u0302+,ls\u0302\u2212,k\zeta (b)$. The action of $s\u0302\u2212,k\zeta (b)$ defined by Eq. (B8) can be defined as $s\u0302\u2212,k\zeta (b)=\zeta k(c)$ (where *K* = *M*/2 and *m* = 2*l*), which has $O(M)$ scaling, and so it is possible to compute

In a similar manner to the scaled Hamiltonian with Eq. (46c) and the number operator with Eqs. (A5) and (A4), we define

which gives

Overall, this reduces the scaling from $O(M3)$ to $O(M2)$. This new algorithm was then used alongside the original speed improvements giving a large time speedup (see Table 8 in the supplementary material for full details). However, it should be possible to scale the $\zeta (a)s\u0302+s\u0302\u2212\zeta (b)$, removing the need to calculate $\zeta (c)$. We can set *m* = 2*l* and *n* = 2*k* and then define

which can have the following outcomes:

which gives

We then have to define the following recursion relations:

These can all trivially be calculated in $O(M)$ steps. We, then, also define

Equations (B16a) and (B16b) can be calculated recursively in $O(M)$ steps and Eq. (B16c) in $O(M2)$. Therefore,

The naïve total spin algorithm was made up of three separate algorithms $S\u0302+S\u0302\u2212$, $S\u0302z$, and $S\u0302z2$ that scaled $O(M3)$, $O(M2)$, and $O(M3)$, respectively. The scaled algorithm reduced the scaling of each algorithm to $O(M2)$, $O(M)$ and $O(M2)$ respectively. This scaling improvement resulted in $S\u03022$ for 1000 orbitals being calculated over 460 times faster than the naïve algorithm; full code racing details are given in Tables 8 and 9 of the supplementary material.

## REFERENCES

Providing that the starting wavefunction contains contributions from wavefunctions of the same symmetry and number of electrons as the lowest-energy state.