The Fermionic Quantum Emulator
Nicholas C. Rubin1 , Klaas Gunst2 , Alec White2 , Leon Freitag2 , Kyle Throssell2 ,Garnet Kin-Lic Chan3 , Ryan Babbush1 , and Toru Shiozaki2
1

Google Quantum AI, Mountain View, CA, 94043

2

Quantum Simulation Technologies, Inc., Cambridge, MA 02139

3

Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena CA 91125

arXiv:2104.13944v2 [quant-ph] 20 Oct 2021

10/15/2021

The fermionic quantum emulator (FQE) is a collection of protocols for emulating quantum dynamics of fermions efficiently taking advantage of common
symmetries present in chemical, materials, and condensed-matter systems. The
library is fully integrated with the OpenFermion software package and serves
as the simulation backend. The FQE reduces memory footprint by exploiting
number and spin symmetry along with custom evolution routines for sparseand dense Hamiltonians, allowing us to study significantly larger quantum circuits at modest computational cost when compared against qubit state vector
simulators. This release paper outlines the technical details of the simulation
methods and key advantages.

1 Introduction
High accuracy simulation of fermionic systems is an important and challenging problem
and is a major motivation behind current attempts to develop quantum computers [1, 4,
12, 31]. There has been significant experimental progress towards realizing the simulationof fermionic systems on current quantum devices [3, 19]. As these experiments scale insize there is a growing need to understand the possibilities for quantum advantage, with
one approach being to characterize the classical emulation complexity of the corresponding
quantum circuits. In addition, the efficient emulation of near-term fermionic simulation
experiments is crucial for experiment design, algorithm design, and testing. In this work, we
describe an implementation of protocols to efficiently emulate quantum circuits describing
time evolution under fermionic generators. We name the library that implements these
protocols the Fermionic Quantum Emulator (FQE).1
There have been many developments in quantum circuit simulation and emulation.
Broadly these advancements can be classified into algorithmic improvements [6, 7, 14, 14,
16–18, 26] and computational implementation improvements [24, 37, 39, 41]. However,despite this progress, there remains potential to improve the emulation of circuits relevant
Nicholas C. Rubin: Corresponding author: nickrubin@google.com
Garnet Kin-Lic Chan: Corresponding author: garnetc@caltech.edu
Ryan Babbush: Corresponding author: babbush@google.com
Toru Shiozaki: Corresponding author: shiozaki@qsimulate.com
1
We make a distinction between simulating quantum circuits directly and emulating the quantum
circuits when the quantum circuit corresponds to dynamics of fermions.

Accepted in Quantum 2021-09-21, click title to verify. Published under CC-BY 4.0.

1

to fermionic simulation. For example, in general circuit emulation it is necessary to include
additional circuit elements to handle the fermion encoding, but these can be eliminated andthe fermionic sign accounted for implicitly in a specialized fermionic emulator. In addition,
many fermionic systems have symmetries that can be used to reduce the resource requirements of the classical emulation. For example, molecular and material problems are often
described by Hamiltonians that commute with a variety of global symmetry operators,
such as the total particle number, total spin, time-reversal, and point-group and crystallographic symmetries. Though there are many open source software packages to carry out
quantum chemistry calculations of fermionic systems using these symmetries [35, 38, 40],these are not designed to support computations using the quantum circuit model. Similarly, existing quantum circuit simulators and the corresponding computational techniques
used within them are not naturally suited to efficiently working within symmetry reduced
Hilbert spaces.
The FQE is a simulator of fermions represented in second quantization. It corresponds
to a statevector simulator in that the wavefunction is explicitly stored. In the first releaseof the emulator, number and spin (Sz ) symmetry are used to minimize the wavefunction
memory footprint by using a generalization of the Knowles–Handy [21] determinant based
configuration interaction wavefunction organization scheme. Time evolution of the state
under a fermionic Hamiltonian is implemented in two ways: (1) via a custom routine similar to the cosine-sine construction for nilpotent operators applicable to sparse fermion
generators, and (2) via a series expansion of the propagator for dense Hamiltonians taking
advantage of efficient intermediate data structures. The library completely integrates withthe fermion quantum simulation library OpenFermion [27] and has built in functions toconvert wavefunction objects in FQE to the format used in OpenFermion as well as for thegeneral purpose circuit simulator Cirq [9]. For example, the time evolution of a wavefunction can be generated using OpenFermion’s FermionOperator objects or Hamiltonians
defined through FQE utility functions. The FQE is an open source library distributed
under the Apache 2.0 license [33]. The library is implemented in Python to facilitate code
extension and reuse. However, some hot spots have been addressed with high-performance
kernels written in the C language. These C extensions can be enabled or disabled by
the user. This allows the library to function as either a high-performance library forcomputationally demanding applications or as a pure-Python library implemented with
understandable and easily modifiable expressions.
This work describes the key technical aspects of the library, demonstrates how it canbe used in quantum algorithm development, and makes single-thread and multi-threaded
timing comparisons for key circuit primitives against Cirq [9] and the highly performant
general purpose quantum circuit simulator Qsim [32]. We close with a perspective on thedevelopment of the library and future directions.

2 Wavefunctions and Hilbert space organization
When simulating spin- 12 fermions in second quantization, the many-particle fermionic
Hilbert space (sometimes called the Fock space) generated by a basis of L single-fermion
states (termed orbitals) {|φ1 i, |φ2 i, . . . |φL i}, is spanned by basis states (termed determinants) labelled by L-digit binary strings {n1 n2 . . . nL } where n ∈ {0, 1}, and ni is theoccupancy of orbital i. It is evident that the fermionic Hilbert space is isomorphic tothe Hilbert space of L qubits, with the qubit computational basis corresponding to thefermion determinant basis. This allows one to use the traditional state vector representation of qubits to encode fermionic states.
Accepted in Quantum 2021-09-21, click title to verify. Published under CC-BY 4.0.

2

For many molecular and materials systems, the electronic Hamiltonian is an operator
that commutes with various global symmetry operators. These symmetries thus provide
useful quantum numbers to label sectors of Hilbert space. Because the Hamiltonian does
not mix sectors, many simulations can be performed either in a single sector, or in a collection of independent sectors. A simple example is simulating a fixed number of electrons
n. In this case, the Hilbert space contains only Ln states. Because the determinants are
eigenstates of the particle number operator, the spanning basis can be chosen to be deterP
minants with total occupancy i ni = n. If we further assume the orbitals are eigenstatesof Sz (i.e. spin-orbitals), then the determinants are also eigenstates of Sz . Then a Hilbert
space sector labelled by n and Sz is spanned by determinants only with the given n andSz .
The FQE uses such a decomposition into symmetry sectors to store the wavefunction
compactly. A user specifies symmetry sectors of fixed particle number and given Sz . The
total wavefunction is then stored in the direct sum of these sectors. Considering all possible
particle sectors and Sz sectors corresponds to working in the full many-particle fermionic
Hilbert space, and thus storing the full set of 2L qubit amplitudes.
An efficient way to represent wavefunctions for a fixed particle number and Sz sector was
first introduced by Siegbahn [36] and refined by Knowles and Handy [21] in the context
of the exact diagonalization of chemistry Hamiltonians. The binary string encoding adeterminant is separated into a string for the spin-up (α) spin-orbitals and the spin-down(β) spin-orbitals. Commonly these two binary strings (integers) are referred to as the αand β-string respectively. A further simplification arises by assuming that the α and β
spin-orbitals share a common set of spatial functions, known as spatial orbitals. Given M
spatial orbitals in total, the total number of spin-orbitals is thus L = 2M .
The following code example initializes a state with four electrons in four spatial orbitals
over a superposition of all possible Sz values.
wfn = fqe . Wavefunction ([[4 , 4 , 4] , [4 , 2 , 4] , [4 , 0 , 4] , [4 , -2 , 4] , [4 , -4 , 4]])Each element of the input list labels a Hilbert space sector as the triple (n-electrons, Sz ,M -spatial-orbitals). The wavefunction in each sector is stored as a matrix of size nMα ×M
nβ where nα and nβ are the total number of α and β electrons. Each row of the matrix
corresponds to a specific α-string, and each column to a β-string. The mapping between
the integer values of the α- and β-strings to the row and column indices is somewhat
arbitrary; we use the lexical ordering proposed by Knowles and Handy [20]. By storing
the wavefunction coefficients as matrices we can leverage vectorized multiplications when
performing updates or contractions on the wavefunction coefficients.
The memory savings from simulating in specific symmetry sectors is substantial whencompared to traditional state-vector simulators which use the full 2L qubit space. In
Figure 1 we show the relative size of half filling (n = L/2) and quarter filling (n = L/4)
subspaces with Sz = 0 along with the wavefunction memory footprint in gigabytes.
To make the FQE interoperable with other simulation tools we provide a number ofutilities for transforming the FQE Wavefunction representation. We also provide humanreadable printing, saving wavefunctions as binary data in numpy’s .npy format, various
wavefunction initialization routines for random, Hartree-Fock, or user defined states, andconversion functions to OpenFermion and Cirq wavefunction representations. The printing
functionality is demonstrated in the code snippet below.

Accepted in Quantum 2021-09-21, click title to verify. Published under CC-BY 4.0.

3

10
10

2

Gigabytes| symm|/| 22M|

104

n/L = 1/2
n/L = 1/2, sz = 0
n/L = 1/4, sz = 0

6

10

10

18

10

14

10

10

10

6

10

2

10

2

20
30
40
Num. Orbitals M [Num. Qubits L / 2]

2L
n/L = 1/2
n/L = 1/2, sz = 0
n/L = 1/4, sz = 0

10

20
30
40
Num. Orbitals M [Num. Qubits L / 2]

Figure 1: left: Relative sizes of Hilbert space sectors at half and quarter filling compared against full
qubit Hilbert space size. right: Gigabytes required to represent spaces of Sz = 0, half and quarter
filling, compared against the full qubit Hilbert space and half filling with no Sz symmetry restriction.

import fqe# two sectors N =2 Sz =[0 , -2] on 4 orbitals
wf = fqe . Wavefunction ([[2 , 0 , 4] , [2 , -2 , 4]])
wf . set_wf ( strategy = ’ random ’)
wf . print_wfn ()Sector N = 2 : S_z = -2
a ’0000 ’ b ’0011 ’ ( 0 . 1 0 8 5 8 0 9 6 1 8 3 7 3 8 3 2 6 - 0 . 3 4 9 6 8 3 2 1 9 2 7 4 3 8 2 3 4 j )a ’0000 ’ b ’0101 ’ ( -0.09372404351124913 -0.19775900492450152 j )... elided output
Sector N = 2 : S_z = 0
a ’0001 ’ b ’0001 ’ ( 0 . 0 9 2 5 1 6 1 1 8 2 0 2 9 7 7 7 7 + 0 . 1 4 5 0 0 9 3 8 2 7 5 1 4 9 1 7 j )a ’0001 ’ b ’0010 ’ ( 0 . 1 5 5 1 6 1 1 1 8 7 5 5 0 0 6 9 3 + 0 . 1 7 7 5 2 4 5 3 7 9 8 3 1 4 6 1 9 j )a ’0001 ’ b ’0100 ’ ( 0 . 1 2 4 9 9 0 2 6 8 8 6 4 8 0 3 1 3 + 0 . 0 2 8 5 4 6 9 0 4 1 7 4 2 5 9 4 3 5 j )... elided output

The conversion to OpenFermion is handled by converting each coefficient in the Wavefunction
object with its associated α, β-string to an ordered string of creation operators acting ona vacuum state, with the string represented by an OpenFermion FermionOperator object.
This allows the user to leverage the normal ordering routines in OpenFermion. By using
the FermionOperator intermediate we can also directly map the string of normal ordered
operators acting on a vacuum to a state vector for integration with Cirq simulations.
of_ops = fqe . fqe _t o_f er mio n_ op ( wf )cirq_wf = fqe . to_cirq ( wf )new_fqe_wf = fqe . from_cirq ( cirq_wf , thresh =1.0 E -12)

# same as original wf

3 Unitary evolution
3.1 Evolution of Sparse Hamiltonians
The FQE can evolve a state by any fermionic excitation Hamiltonian. A fermionic excitation Hamiltonian is of the formĤexcite = ĝ + ĝ †



Nop

ĝ =g

(1)
Nop

Y †

Y

p=1

q=1

âip σpâjq ρq

Accepted in Quantum 2021-09-21, click title to verify. Published under CC-BY 4.0.

(2)

4

where g is a complex number and ĝ is a single product of an arbitrary number (Nop ) of creation and annihilation operators that create/annihilate orbitals with spatial orbital indiceslabelled by i, j and Sz indices labelled by σ, ρ. We refer to this Hamiltonian as the excitation Hamiltonian because it only involves a single “excitation” term ĝ and its Hermitian conjugate. To specify ĝ the FQE can digest a FermionOperator from OpenFermion. If Ĥexcite
is diagonal, i.e. a polynomial of number operators (such as â†1α â†2α â1α â2α = −n̂1α n̂2α , withn̂1α = â†1α a1α and similarly for n̂2α ) evolution of the wave function can be performed using
the techniques described later in Sec. 3.3.1. When Ĥexcite is not diagonal and contains norepeated indices (such as â†4α â†2β â1α â3β ), evolution of the wavefunction is accomplished bye−i(ĝ+ĝ † )

= 1 + cos(|g|) − 1 − iĝ† sin(|g|)P̂ĝĝ†|g|sin(|g|)+ cos(|g|) − 1 − iĝP̂ĝ† ĝ|g|(3)

where we define P̂x as the projector onto the basis of determinants that are not annihilated by x̂. For other recent applications of this relation in quantum chemistry see e.g.
Refs. [13], [11], and Ref. [8] in the context of quantum and quantum-inspired algorithms.
A derivation of Eq. (3) is presented in Appendix A. When Ĥexcite is not diagonal but has
several repeated indices as well (such as â†1α â†2β â1α â3β = −n̂1α â†2β â3β ), a hybrid approachis used with the only additional complication being the fermion parity evaluation thatarises from operator reordering.
Time evolution of a FQE wavefunction is implemented as a method of the Wavefunction
object and can be easily accessed through the Wavefunction object interface. For example,
the code snippetfrom openfermion import FermionOperator
import fqe
wf = fqe . Wavefunction ([[4 , 0 , 6]])
wf . set_wfn ( strategy = ’ random ’)i , j , theta = 0 , 1 , 2 / 3
op = ( FermionOperator (((2 * i , 1) , (2 * j , 0) ) , coefficient = -1 j * theta ) +FermionOperator (((2 * j , 1) , (2 * i , 0) ) , coefficient =1 j * theta ) )new_wfn = fqe_wfn . time_evolve (1.0 , op )performs evolution of a random state ψ0 with a one-particle excitation Hamiltonian acting
on the α-spin sector††|ψf i = e−iθ(âiα âjα +âjα âiα ) |ψ0 i.

(4)

The Wavefunction object correctly evolves this type of sparse Hamiltonian so long as the
Hamiltonian commutes with the underlying wavefunction symmetries. However, the user
is responsible for providing a Hamiltonian of this form, and the FQE does not correctly
handle the evolution of Hamiltonians which break the specified wavefunction symmetry.
Evolving under a single fermionic excitation Hamiltonian is the basis for arbitrary
fermionic quantum circuit evolution and can be used to simulate arbitrary fermionic Hamiltonian evolution. Because of the isomorphism between fermion and qubit spaces, general
qubit Hamiltonians, so long as they are symmetry preserving in the same sense, can be simulated with little overhead within the same scheme either by modifying the code to ignore
signs arising from fermion parity or by explicitly inserting swap operations. This raises
the possibility to simulate large quantum circuits associated with non-fermionic Hamiltonian evolution, that benefit from the memory saving associated with simulating the fixed“particle” or “spin” sectors of the corresponding qubit Hilbert space.
Accepted in Quantum 2021-09-21, click title to verify. Published under CC-BY 4.0.

5Though evolution by single excitation Hamiltonian provides a primitive to implement
many fermionic circuit simulations, the FQE also implements efficient time-evolution routines for Hamiltonians with special forms and for generic quantum chemical Hamiltonians.
These will be discussed in more detail in the following sections.

3.2 Evolution of Dense Hamiltonians
The FQE provides special routines to evolve sums of excitation Hamiltonians through series
expansions. We begin by discussing dense Hamiltonians, which for the purposes of thiswork we define as a weighted sum over all possible excitation Hamiltonians with the sameindex structure. A dense two-particle Hamiltonian, for which the FQE implements special
routines, would take the form
Ĥ =

X X

Viσ,jρ,kσ0 ,lρ0 â†iσ â†jρ âkσ0 âlρ0(5)ijkl σσ 0 ρρ0where Ĥ is Hermitian. In addition, specialized code is implemented for spin-conserving
spin-orbital Hamiltonians (σ = σ 0 and ρ = ρ0 in Eq. (5)),
Ĥ =

XX

Viσ,jρ,kσ,lρ â†iσ â†jρ âkσ âlρ .

(6)ijkl σρand for spin-free Hamiltonians (Vijkl ≡ Viσ,jρ,kσ0 ,lρ0 in Eq. (6)),
Ĥsf =

XX

Vijkl â†iσ â†jρ âkσ âlρ ,(7)ijkl σρwhich frequently appear in molecular and materials systems. Other cases which havespecialized subroutines discussed below include sparse Hamiltonians (arbitrary sums of
excitation Hamiltonians), quadratic Hamiltonians, and diagonal pair Hamiltonians.
Time evolution with such dense Hamiltonians is performed by means of series expansions such as the Taylor and Chebyshev expansions. When the time step t is taken to besmall, the operator exponential can be efficiently computed by a Taylor expansion thatconverges rapidly,exp(−iĤt)|Ψi =

X 1
n

n!

(−iĤt)n |Ψi.

(8)

We evaluate this by recursively computing the action of the operator on the wave functions,|Ψn+1 i = Ĥ|Ψn i(9)with |Ψ0 i ≡ |Ψi, using which the previous expression can be rewritten asexp(−iĤt)|Ψi =X (−it)n
n

n!

|Ψn i

(10)

We evaluate the norm of each term, |tn |||Ψn i|/n!, and when the norm becomes smaller
than the given threshold, we stop the computation. Alternatively, one can specify the
number of terms in the expansion, or both.
When the spectral range of the operator (i.e., extremal eigenvalues [min , max ]) is known
or can be estimated in advance, the expansion can be made more robust by using the

Accepted in Quantum 2021-09-21, click title to verify. Published under CC-BY 4.0.

6

Chebyshev expansion [22]. The Chebyshev expansion is known to be a near optimal approximation of the exponential in a minimax sense. First, we scale the operator such thatthe eigenvalues are in [−1, +1],
Ĥ 0 =

Ĥ + shift
∆

(11)

where ∆ = (max − min )/2w0 and shift = −(max + min )/2. The factor w0 is a numberthat is slightly smaller than 1 to ensure that the eigenvalues do not lie outside the window
due to numerical noise or insufficient accuracy in the estimation of extremal eigenvalues.
In the FQE, we use w0 = 0.9875. Then we expand asexp(−iĤt)|Ψi = exp(ishift t) exp(−i∆Ĥ 0 t)|Ψi= exp(ishift t)

X

2Jn (∆t)(−i)n |Ψ̃n i

(12)

n

Here Jn (x) is the modified Bessel function of the first kind, and |Ψ̃n i is obtained byrecursion as|Ψ̃n+1 i = 2Ĥ 0 |Ψ̃n i − |Ψ̃n−1 i(13)

We stop the expansion when the contribution from rank n becomes smaller than the given
threshold, as in the Taylor expansion above.
The dense evolution routines can be accessed through the Wavefunction interface
which intelligently dispatches to the specialized routines depending on the form of the
Hamiltonian the user provides. As an example, below is a code snippet assuming the user
has defined an OpenFermion MolecularData object called ‘molecule’. The Wavefunction
interface provides access to the apply_generated_unitary routine which has options foreach series expansion evolution described above.
from fqe . o pen fe rmi on _ut il s import i n t e g r a l s _ t o _ f q e _ r e s t r i c t e d
dt = 0.23 # evolution timeoei , tei = molecule . get_integrals ()# A dense fqe . R e s t r i c t e d H a m i l t o n i a n object
fqe_rham = i n t e g r a l s _ t o _ f q e _ r e s t r i c t e d ( oei , tei )wfn = fqe . Wavefunction ([[ molecule . n_electrons , 0 , molecule . n_orbitals ]])wfn . set_wfn ( strategy = ’ hartree - fock ’)
new_wfn1 = wfn . time_evolve ( dt , fqe_rham )# equivalent path
# Taylor expansion algorithm
new_wfn2 = wfn . a p p l y _ g e n e r a t e d _ u n i t a r y ( dt , ’ taylor ’ , fqe_rham )The above algorithms that use the series expansion are also employed for time evolution with sparse Hamiltonians, which are arbitrary weighted sums of multiple excitation
Hamiltonians. In this case, the apply function used in Eqs. 9 and 13 is overloaded by aspecial function for sparse Hamiltonians, in which the action of the operator is evaluated
one term at a time.

3.3 Special routines for other structured Hamiltonians
3.3.1

Diagonal pair Hamiltonians

Hamiltonians that are diagonal in the determinant basis are particularly simple as they
are polynomials of number operators with terms that all commute. Evolution under such
Hamiltonians appears as an important quantum circuit for chemical Hamiltonian Trotter
Accepted in Quantum 2021-09-21, click title to verify. Published under CC-BY 4.0.

7

steps and certain variational ansätze [23, 29]. The FQE leverages the mutual commuting
nature of the terms of Hamiltonians of this form to efficiently evolve under the generated
unitaries. Specifically, the action of the unitary generated by diagonal pair Hamiltoniansof the formĤdiag =X

Wrs n̂r n̂s ,

(14)

rs

where n̂r is the spin-summed number operator for spatial orbital r (i.e., n̂r = σ â†rσ ârσ ),is constructed by iterating over all determinants and generating a phase associated withthe coefficient Wrs for each determinant. A direct wall clock time comparison to othersimulators is made in Figure 2 by translating evolution under a random diagonal pair
Hamiltonian (of the form in Eq. (14)) into a series of CPHASE gates with phases corresponding to 2Wrs . For the half filling circuits on 14 orbitals (i.e., 14 electrons in 28 spin
orbitals to be represented by 28 qubits) we observe that the python reference implementation of FQE is approximately eight times faster than Cirq and five times faster thansingle-threaded Qsim. The C-accelerated FQE single threaded is five hundred times fasterthan Qsim single threaded with default gate fusion. The C-accelerated FQE on twenty
threads is twenty times faster than QSim using twenty threads and gate-fusion level eight
for the same 14 orbital half filled system. FQE performance gains for lower filling fractions
are even more substantial.

10

1

10

1

10

3

Cirq
FQE-Half
FQE-Quarter

2

4

6

Wall Clock Time [s]Wall Clock Time [s]P

8 10 12
Num. Orbitals

14

16

10

1

10

1

10

3

Qsim-1-8
Qsim-20-8
FQE-C-half-1
FQE-C-half-20
FQE-C-quarter-1

2

4

6

8 10 12
Num. Orbitals

14

16

Figure 2: Run time comparisons for implementing time evolution of a random diagonal Hamiltonian:left The Python reference implementation of FQE is compared against Cirq’s Python circuit evolution
routines; right C-accelerated FQE is compared against Qsim, a highly optimized C++ circuit simulator
with and without threading. FQE computes the results in-place. Qsim is executed in single precisionmode with one (green dashed line) and twenty (green solid line) threads with gate gate-fusion set toeight through the qsimcirq python interface. Quantum circuits for the time evolution generated by the
diagonal Coulomb Hamiltonian are constructed using 4M 2 CPHASE gates where M is the number ofspatial orbitals. FQE at half filling with one and twenty threads is shown in blue while quarter (single
thread) filling is in red.

3.3.2

Quadratic Hamiltonians

Evolving states by quadratic Hamiltonians is another important algorithmic primitive. In
the quantum circuit context, it is closely related to matchgate circuit simulation [42]. The
FQE performs this task by transforming the wave functions into the orbital representation that diagonalizes the quadratic Hamiltonian, evolving in the quadratic Hamiltonian
eigenbasis, and then returning to the original basis. The algorithm is based on those in
the quantum chemistry literature [5, 25, 28] that utilize the LU decomposition, with an
Accepted in Quantum 2021-09-21, click title to verify. Published under CC-BY 4.0.

8

improvement in the handling of pivoting in the LU decomposition. Using the unitary X
Pthat diagonalizes the operator matrix in Â = ij (A)ij â†i âj , i.e., X† AX = a, one obtains
exp(−iÂt)|Ψi = T̂ (X† ) exp(−iât)T̂ (X)|Ψi

(15)

10

3

10

2

10

1

10

0

10

1

10

2

Cirq
FQE-Python-Half
FQE-Python-quarter

2

4

6

8 10 12
Num. Orbitals

Wall Clock Time [s]Wall Clock Time [s]where â is the diagonal operator after orbital rotation, and T̂ (X) is the transformationon the wave function due to the change in the orbital basis described by the unitary X.
Here we ignore the pivoting for brevity; see Appendix B for more algorithmic details. The
overall cost of rotating the wave function (i.e., action of T̂ (X)) is roughly equivalent to thecost to apply a single dense one-particle linear operator to a wave function. The diagonal
operator is then applied to the rotated wave function, followed by the back transformation
to the wave function in the original orbital basis.
Figure 3 compares the wall clock time required to evolve a random quadratic Hamiltonian. For Cirq and Qsim, a circuit is generated using the optimal Givens rotation decomposition construction of Clements [10] in OpenFermion which is translated into a sequenceof 2M 2 gates of the form exp (−iθ(XY − Y X)/2) and Rz. For the half filling circuits onfourteen orbitals (twenty eight qubits) we observe that the python-FQE LU decomposition
algorithm is approximately forty times faster than Cirq. The single-thread C-accelerated
FQE is approximately ten times faster than single-threaded Qsim with gate-fusion set tofour. For twenty threads the C-accelerated FQE is five times faster than Qsim. Lower
filling fractions show even more substantial savings. This illustrates the large advantage
physics, symmetry, and algorithmic considerations can provide when targeting specialized
quantum circuit emulation.

14

16

10

2

10

1

10

0

10

1

10

2

10

3

Qsim-1-4
Qsim-20-4
FQE-C-1-half
FQE-C-20-half
FQE-C-1-quarter

2

4

6

8 10 12
Num. Orbitals

14

16

Figure 3: Wall clock time comparison for a random quadratic Hamiltonian evolution at half filling andquarter filling, left between Cirq and the Python implementation in FQE and right between Qsim and Caccelerated FQE. Note that Cirq and Qsim run times are unaffected by filling fraction. Qsim is executed
accessed through the qsimcirq Python interface. Each Qsim trace is labeled by how many threads areused (one or twenty) and the gate-fusion level (four). Quantum circuits for the time evolution generated
by the quadratic Hamiltonian are constructed using 2M 2 exp (−iθ(XY − Y X)/2) gates where M is
the number of spatial orbitals. FQE evolves the quadratic Hamiltonian in Section 3.3.2.

4 Operator Action on the Wavefunction
The FQE provides a variety of methods to apply fermionic operators to states. The general code allows a user to define an arbitrary individual excitation operator and apply it
to a wavefunction. This functionality constructs the action of the user defined operator
by computing the sparse matrix elements of such an operator in the determinant basis,

Accepted in Quantum 2021-09-21, click title to verify. Published under CC-BY 4.0.

9

which are then contracted against the wavefunction coefficients to obtain the action of theoperator on the state.
The FQE also provides efficient methods to apply dense fermionic operators (such as adense Hamiltonian) to wave functions, because these are the fundamental building blocks
for fermionic time evolution under dense Hamiltonians (see, for example, Eqs. (9) and (13)).
It is equally important for reduced density matrix (RDM) construction, efficient access towhich is a crucial aspect of many fermionic simulation algorithms. In the following, we
present the algorithms to apply dense spin conserving, spin-free Hamiltonians (Eq. (7)),
but the FQE implements similar algorithms for Sz conserving spin-orbital Hamiltonians(Eq. (6)) and Sz non-conserving Hamiltonians (Eq. (5)) as well.

4.1 Knowles–Handy algorithm
In the Knowles–Handy algorithm [21], a resolution of the identity (RI) is inserted in theP
n-electron determinant space after operator reordering (1̂ = K |KihK|) as
hI|Ĥsf |Ψi = −

XX



X
hI|â† âkσ |Ki  Vijkl hK|â† âlρ |JiCJ iσjρσρ ik,K

jl,I

!

+

Vikkl hI|â†iσ âlσ |JiCJ ,

XX X
σ

il

(16)

k

where I, J, and K label determinants and CI is the wave function coefficient associatedwith determinant |Ii. The second term is trivially evaluated. The primitive operations forthe first term are
J
Dij
=
J
D̃ij
=

XX

hJ|â†iσ âjσ |IiCI ,

I

σ

X

J
Vikjl Dkl
,

(17)
(18)

kl

RI =

XX
J

J
.
hI|â†iσ âjσ |JiD̃ij

(19)

σ

This form of the RI insertion has the advantage that it can be easily generalized to apply
3- and 4-particle operators, since these primitive operations can be performed recursively,i.e.,J
Eij,kl=(20)

JhI|â†kσ âlσ |JiEij,kl
.

(21)

σ

I
I
Dij
=

I
hJ|â†iσ âjσ |IiDkl
,

XX

XX
J

σ

The step in Eq. (18) is typically a computational hot spot and is performed using efficient
numpy functions.
For efficiency, the expectation values hIσ |â†iσ âjσ |Jσ i are precomputed for the α- andβ-strings and stored in an FQE object FciGraph. To facilitate the evaluation, the interI in Eq. (17)mediate tensors are also stored using the α- and β-strings; for example, Dij
I I

is stored as a four-index tensor Dijα β where I = Iα ⊗ Iβ . When the Hamiltonian breaks
spin symmetry (and, hence, Sz symmetry), the mappings between strings that have Nσ ± 1
electrons are also generated (hIσ0 |â†iσ |Jσ i and hIσ00 |âiσ |Jσ i) which are stored in the FQE
object FciGraphSet. From these quantities, hI|â†iσ âjσ |Ji can be easily computed on thefly.
Accepted in Quantum 2021-09-21, click title to verify. Published under CC-BY 4.0.

10

4.2 Harrison–Zarrabian algorithm for low filling
Low filling cases are computed using the Harrison–Zarrabian algorithm [15], which is closely
related to the algorithm developed by Olsen et al. [30]. This algorithm inserts the RI inP
the n − 2 electron space (1̂ = L |LihL|) as

hI|Ĥsf |Ψi =XX

hI|â†iσ â†jρ |Li σρ ij,L
X

Vijkl hL|âkσ âlρ |JiCJ  ,

(22)

kl,J

where I and J label a determinant with n electrons and L labels determinants with n − 2electrons. This is computed by the following primitive steps,
LFij,σρ
=X

hL|âkσ âlρ |IiCI ,(23)

L
Vijkl Fkl,σρ
,

(24)

I
L
F̃ij,σρ
=

X
kl

RI =

L
hI|â†iσ â†jρ |LiF̃ij,σρ
.

XX
L

(25)

σρ

The advantage of this algorithm, in comparison to the Knowles–Handy algorithm, is thatthe number of determinants for the RI can be far smaller than the number of the original
determinants. This is pronounced when the wave function is at low filling. The disadvantage, however, is that there is certain overhead in computational cost because one cannot
perform spin summation in Eq. (23). Therefore, the FQE switches to this algorithm when
the filling is smaller than 0.3. Note that a similar algorithm can be devised for the highfilling cases by sorting the operators in the opposite order, but this is not implemented in
the FQE. When such a system is considered, one can reformulate the high-filling problem
into a low-filling one by rewriting the problem in terms of the hole operators b̂iσ = 1 − âiσ .

4.3 Olsen’s algorithm
The approach of Knowles and Handy, while being simple to implement with optimized
linear algebra routines in Python, can be further improved by taking advantage of the
J appearing in Equation 17. For example, D J will be zero if determinantsparsity in the Dij
ij
J has orbital i unoccupied or orbital j doubly occupied. This means that extra work is
performed during the most expensive step of the algorithm, Equation 18. Olsen and coworkers were able to leverage this sparsity by avoiding explicit reference to intermediate
states [30]. In this approach, the same replacement lists used by Knowles and Handy in
J are used to loop over only the non-zero Hamiltonian elements
the construction of the Dij
connecting bra and ket. The essence of the algorithm is to structure the loops in sucha way that vector operations can be used. Details of the algorithm are given by Olsenet al. [30] and clearly reviewed by Sherrill and Schaefer [34]. While Olsen’s algorithm
cannot be efficiently implemented in pure Python, it is the basis of the FQE C extension
for applying dense, two-particle Hamiltonians.

4.4 RDM computation
The FQE provides efficient routines using the intermediate tensors in the above algorithms
to efficiently compute 1- through 4-RDMs. For example, a spin-summed 2-RDM element
Γijkl defined asΓijkl =

XhΨ|â†iσ â†jρ âkσ âlρ |Ψi

(26)

σρ

Accepted in Quantum 2021-09-21, click title to verify. Published under CC-BY 4.0.

11can be computed either using the intermediate tensor in the Knowles–Handy algorithm(Eq. (17)),
Γijkl = −

X

I∗ I
Dki
Djl + δjk

I

X

DliI∗ CI ,

(27)

I

or, in low filling cases, using the intermediate in the Harrison–Zarrabian algorithm,
Γijkl = −

XX
σρ

L∗
L
Fij,σρ
Fkl,σρ
.

(28)

L

The corresponding generalization of Olsen’s algorithm is used when the C extensions areenabled. In both Python algorithms, efficient numpy functions for matrix–matrix andmatrix–vector multiplication are used. Spin-orbital RDMs (i.e. without the spin summation in Eq. (26)) are computed similarly. In addition to the standard (i.e., particle) RDMs,
the FQE can compute RDMs that correspond to other operator orders, e.g., hole RDMs.
This is done by performing Wick’s theorem as implemented in fqe.wick and evaluating
the resulting expression using the particle RDMs. For example, the spin-summed 2-bodyhole RDM element (Γhijkl ) defined asΓhijkl =X

hΨ|âiσ âjρ â†kσ â†lρ |Ψi

(29)

σρ

is computed asΓhijkl = Γklij + δjk Γil − 2δik Γjl − 2δjl Γik .

(30)

The FQE implements the spin-orbital counterpart of this feature as well.
For spin conserving Hamiltonians, the FQE also provides up to 4-body spin-summed
RDMs and up to 3-body spin-orbital RDMs in the OpenFermion format. An example is:
wf = fqe . Wavefunction ([[4 , 0 , 6]])
wf . set_wfn ( strategy = ’ random ’)
spin_sum_opdm = wf . ex pectationValue ( ’i ^ j ’)
spin_sum_oqdm = wf . ex pectationValue ( ’i j ^ ’)
spin_sum_tpdm = wf . ex pectationValue ( ’i ^ j ^ k l ’)opdm , tpdm = wf . sector ((4 , 0) ) . g e t _ o p e n f e r m i o n _ r d m s ()d3 = wf . sector ((4 , 0) ) . get_d3 () # returns spinful 3 - RDM in openfermion ordering

The implementation to compute RDMs offers a significant improvement over OpenFermion’s
native RDM generators which map Pauli operator expectations to RDM elements. As an
example of this efficient RDM computation the FQE library provides a base implementation of the energy gradient due to a two-particle generator Ĝiσ,jρ,kσ0 ,lρ0 = â†iσ â†jρ âkσ0 âlρ0 ,which corresponds to evaluating giσ,jρ,kσ0 ,lρ0 = hΨ|[Ĥ, Ĝiσ,jρ,kσ0 ,lρ0 ]|Ψi for all spatial andspin indices.

5 Closing thoughts and future directions
The Fermionic Quantum Emulator (FQE) described above is an open source library underthe Apache 2.0 license [33]. Currently, the library is completely implemented in Python
to facilitate extension and code reuse. Despite not being written in a high performance
programming language, the FQE’s algorithmic advantages allow us to outperform evenheavily optimized quantum circuit simulators. In future releases, current computationalbottlenecks will be addressed with additional performance improvements. The FermionicAccepted in Quantum 2021-09-21, click title to verify. Published under CC-BY 4.0.

12Quantum Emulator (FQE) provides a user-friendly and developer-friendly code stack without sacrificing performance. Algorithmic advantages allow the FQE to clearly outperform
both Python-based and heavily optimized quantum circuit simulators, even at half filling,with our C-accelerated code. Away from half filling, the efficiency of FQE is even more
pronounced. Extending the code to GPUs or other architectures may lead to even better
performance.
Though this first release of the FQE focuses on exact evolution in a statevector representation, our long-term goal is to provide implementations of fermionic circuit emulationthat exploit symmetry within different simulation strategies. For example, approximate
representation using fermionic matrix product states, can also be adapted to this setting.
Ultimately, a variety of such techniques will need to be explored to reach an honest assessment of possible quantum advantage when simulating molecular and materials systems.

6 Acknowledgments
NCR thanks Ryan LaRose for testing and feedback on the code.

References[1] Daniel S Abrams and Seth Lloyd. Simulation of Many-Body Fermi Systems
on a Universal Quantum Computer. Phys. Rev. Lett., 79(13):4, 1997. DOI:
https://doi.org/10.1103/PhysRevLett.79.2586.
[2] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz,
A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen. LAPACK Users’
Guide. Society for Industrial and Applied Mathematics, Philadelphia, PA, third
edition, 1999. ISBN 0-89871-447-8 (paperback). URL https://www.netlib.org/
lapack/lug/.
[3] Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C. Bardin, Rami
Barends, Sergio Boixo, Michael Broughton, Bob B. Buckley, David A. Buell, Brian
Burkett, Nicholas Bushnell, Yu Chen, Zijun Chen, Benjamin Chiaro, Roberto
Collins, William Courtney, Sean Demura, Andrew Dunsworth, Edward Farhi, Austin
Fowler, Brooks Foxen, Craig Gidney, Marissa Giustina, Rob Graff, Steve Habegger, Matthew P. Harrigan, Alan Ho, Sabrina Hong, Trent Huang, William J. Huggins, Lev Ioffe, Sergei V. Isakov, Evan Jeffrey, Zhang Jiang, Cody Jones, Dvir
Kafri, Kostyantyn Kechedzhi, Julian Kelly, Seon Kim, Paul V. Klimov, Alexander
Korotkov, Fedor Kostritsa, David Landhuis, Pavel Laptev, Mike Lindmark, Erik
Lucero, Orion Martin, John M. Martinis, Jarrod R. McClean, Matt McEwen, Anthony Megrant, Xiao Mi, Masoud Mohseni, Wojciech Mruczkiewicz, Josh Mutus,
Ofer Naaman, Matthew Neeley, Charles Neill, Hartmut Neven, Murphy Yuezhen
Niu, Thomas E. O’Brien, Eric Ostby, Andre Petukhov, Harald Putterman, Chris
Quintana, Pedram Roushan, Nicholas C. Rubin, Daniel Sank, Kevin J. Satzinger,
Vadim Smelyanskiy, Doug Strain, Kevin J. Sung, Marco Szalay, Tyler Y. Takeshita,
Amit Vainsencher, Theodore White, Nathan Wiebe, Z. Jamie Yao, Ping Yeh, andAdam Zalcman. Hartree-fock on a superconducting qubit quantum computer. Science, 369(6507):1084–1089, 2020. ISSN 0036-8075. DOI: 10.1126/science.abb9811.
URL https://science.sciencemag.org/content/369/6507/1084.
[4] Alan Aspuru-Guzik, Anthony D Dutoi, Peter J Love, and Martin Head-Gordon. Simulated Quantum Computation of Molecular Energies. Science, 309(5741):1704, 2005.
DOI: 10.1126/science.1113479.
Accepted in Quantum 2021-09-21, click title to verify. Published under CC-BY 4.0.

13

[5] Gregory J. Atchity and Klaus Ruedenberg. Orbital transformations and configurational transformations of electronic wavefunctions. J. Chem. Phys., 111(7):2910–2920,
1999. DOI: 10.1063/1.479573. URL https://doi.org/10.1063/1.479573.
[6] Sergio Boixo, Sergei V Isakov, Vadim N Smelyanskiy, and Hartmut Neven. Simulationof low-depth quantum circuits as complex undirected graphical models. arXiv preprint
arXiv:1712.05384, 2017. URL https://arxiv.org/abs/1712.05384.
[7] Sergey Bravyi, Dan Browne, Padraic Calpin, Earl Campbell, David Gosset, and MarkHoward. Simulation of quantum circuits by low-rank stabilizer decompositions. Quantum, 3:181, 2019. URL https://dx.doi.org/10.22331/q-2019-09-02-181.
[8] Jia Chen, Hai-Ping Cheng, and James K. Freericks. Quantum-inspired algorithm forthe factorized form of unitary coupled cluster theory. J. Chem. Theory Comput., 17(2):841–847, 2021. DOI: 10.1021/acs.jctc.0c01052. URL https://doi.org/10.1021/
acs.jctc.0c01052.
[9] Cirq Developers. Cirq, March 2021. URL https://doi.org/10.5281/zenodo.
4586899. See full list of authors on Github: https://github .com/quantumlib/Cirq/graphs/contributors.
[10] William R Clements, Peter C Humphreys, Benjamin J Metcalf, W Steven Kolthammer, and Ian A Walmsley. Optimal design for universal multiport interferometers. Optica, 3(12):1460–1465, 2016. DOI: 10.1364/OPTICA.3.001460. URL https://www.osapublishing.org/optica/fulltext.cfm?uri=optica-3-12-1460.
[11] Francesco A. Evangelista, Garnet Kin-Lic Chan, and Gustavo E. Scuseria. Exact
parameterization of fermionic wave functions via unitary coupled cluster theory. J.
Chem. Phys., 151(24):244112, 2019. DOI: 10.1063/1.5133059. URL https://doi.
org/10.1063/1.5133059.
[12] Richard P Feynman. Simulating physics with computers. Int. J. Theor. Phys., 21(6-7):467–488, 1982. ISSN 00207748. DOI: 10.1007/BF02650179.
[13] Maria-Andreea Filip and Alex J. W. Thom. A stochastic approach to unitary coupledcluster. J. Chem. Phys., 153(21):214106, 2020. DOI: 10.1063/5.0026141. URL https://doi.org/10.1063/5.0026141.
[14] Johnnie Gray and Stefanos Kourtis. Hyper-optimized tensor network contraction.
Quantum, 5:410, 2021. URL https://dx.doi.org/10.22331/q-2021-03-15-410.
[15] Robert J. Harrison and Sohrab Zarrabian. An efficient implementation of the full-ci
method using an (n–2)-electron projection space. Chem. Phys. Lett., 158(5):393–398,
1989. ISSN 0009-2614. DOI: https://doi.org/10.1016/0009-2614(89)87358-0. URL
https://www.sciencedirect.com/science/article/pii/0009261489873580.
[16] Cupjin Huang, Fang Zhang, Michael Newman, Junjie Cai, Xun Gao, Zhengxiong
Tian, Junyin Wu, Haihong Xu, Huanjun Yu, Bo Yuan, et al. Classical simulation ofquantum supremacy circuits. arXiv preprint arXiv:2005.06787, 2020. URL https://arxiv.org/abs/2005.06787.
[17] Yifei Huang and Peter Love. Approximate stabilizer rank and improved weak simulation of clifford-dominated circuits for qudits. Phys. Rev. A, 99:052307, May 2019.
DOI: 10.1103/PhysRevA.99.052307. URL https://link.aps.org/doi/10.1103/
PhysRevA.99.052307.
[18] Yifei Huang and Peter Love.
Feynman-path-type simulation using stabilizer
projector decomposition of unitaries.
Phys. Rev. A, 103:022428, Feb 2021.
DOI: 10.1103/PhysRevA.103.022428. URL https://link.aps.org/doi/10.1103/
PhysRevA.103.022428.
[19] Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink,
Jerry M Chow, and Jay M Gambetta. Hardware-efficient variational quantum

Accepted in Quantum 2021-09-21, click title to verify. Published under CC-BY 4.0.

14

eigensolver for small molecules and quantum magnets. Nature, 549(7671):242–246,
2017. DOI: https://doi.org/10.1038/nature23879. URL https://www.nature.com/
articles/nature23879.
[20] Peter J. Knowles and Nicholas C. Handy. A determinant based full configurationinteraction program. Comput. Phys. Comm., 54(1):75–83, 1989. ISSN 0010-4655. DOI:
https://doi.org/10.1016/0010-4655(89)90033-7. URL https://www.sciencedirect.
com/science/article/pii/0010465589900337.
[21] P.J. Knowles and N.C. Handy. A new determinant-based full configuration interaction method. Chem. Phys. Lett., 111(4):315–321, 1984. ISSN 0009-2614. DOI:
https://doi.org/10.1016/0009-2614(84)85513-X. URL https://www.sciencedirect.
com/science/article/pii/000926148485513X.
[22] R Kosloff. Propagation methods for quantum molecular dynamics. Annu. Rev.
Phys. Chem., 45(1):145–178, 1994. DOI: 10.1146/annurev.pc.45.100194.001045. URL
https://doi.org/10.1146/annurev.pc.45.100194.001045.
[23] Joonho Lee, William J Huggins, Martin Head-Gordon, and K Birgitta Whaley. Generalized unitary coupled cluster wave functions for quantum computation. J. Chem. Theory Comput., 15(1):311–324, 2018. DOI: https://doi.org/10.1021/acs.jctc.8b01004.
URL https://pubs.acs.org/doi/10.1021/acs.jctc.8b01004.
[24] Xiu-Zhe Luo, Jin-Guo Liu, Pan Zhang, and Lei Wang. Yao. jl: Extensible, efficient framework for quantum algorithm design. Quantum, 4:341, 2020. DOI:
https://doi.org/10.22331/q-2020-10-11-341. URL https://quantum-journal.org/
papers/q-2020-10-11-341/.
[25] Per Åke Malmqvist. Calculation of transition density matrices by nonunitary
orbital transformations. Int. J. Quantum Chem., 30(4):479–494, 1986. DOI:
https://doi.org/10.1002/qua.560300404. URL https://onlinelibrary.wiley.com/
doi/abs/10.1002/qua.560300404.
[26] Igor L Markov and Yaoyun Shi. Simulating Quantum Computation by Contracting
Tensor Networks. SIAM J. Comput., 38(3):963–981, 2008. ISSN 0097-5397. DOI:
10.1137/050644756. URL http://epubs.siam.org/doi/abs/10.1137/050644756.
[27] Jarrod R McClean, Nicholas C Rubin, Kevin J Sung, Ian D Kivlichan, Xavier BonetMonroig, Yudong Cao, Chengyu Dai, E Schuyler Fried, Craig Gidney, Brendan Gimby,
et al. Openfermion: the electronic structure package for quantum computers. Quantum
Science and Technology, 5(3):034014, 2020. DOI: 10.1088/2058-9565/ab8ebc. URL
https://doi.org/10.1088/2058-9565/ab8ebc.
[28] A. Mitrushchenkov and H.-J. Werner. Calculation of transition moments between
internally contracted mrci wave functions with non-orthogonal orbitals. Mol. Phys.,105(9):1239–1249, 2007. DOI: 10.1080/00268970701326978. URL https://doi.org/
10.1080/0026897070132697.
[29] Mario Motta, Erika Ye, Jarrod R McClean, Zhendong Li, Austin J Minnich, Ryan
Babbush, and Garnet Kin-Lic Chan. Low rank representations for quantum simulationof electronic structure. arXiv preprint arXiv:1808.02625, 2018. URL https://arxiv.
org/abs/1808.02625.
[30] Jeppe Olsen, Björn O. Roos, Poul Jørgensen, and Hans J. "Aa". Jensen. Determinant
based configuration interaction algorithms for complete and restricted configurationinteraction spaces. J. Chem. Phys., 89(4):2185–2192, 1988. DOI: 10.1063/1.455063.
URL https://doi.org/10.1063/1.455063.
[31] G Ortiz, J Gubernatis, E Knill, and R Laflamme. Quantum algorithms for fermionicsimulations. Phys. Rev. A, 64(2):22319, 2001. ISSN 1050-2947. DOI: 10.1103/PhysRevA.64.022319. URL http://link.aps.org/doi/10.1103/PhysRevA.64.022319.

Accepted in Quantum 2021-09-21, click title to verify. Published under CC-BY 4.0.

15

[32] Quantum AI team and collaborators. qsim, September 2020. URL https://doi.
org/10.5281/zenodo.4023103.
[33] Nicholas C. Rubin,
Toru Shiozaki,
Kyle Throssell,
Garnet K.-L.
Chan,and Ryan Babbush.
The Fermionic Quantum Emulator:
https://github.com/quantumlib/openfermion-fqe, Aug 2020.
URL https://github.com/quantumlib/OpenFermion-FQE.
[34] C. David Sherrill and Henry F. Schaefer. The Configuration Interaction Method:
Advances in Highly Correlated Approaches. Advances in Quantum Chemistry, 34(C):
143–269, 1999. ISSN 00653276. DOI: 10.1016/S0065-3276(08)60532-8.
[35] Toru Shiozaki.
Bagel: Brilliantly advanced general electronic-structure library.
Wiley Interdiscip. Rev. Comput. Mol. Sci., 8(1):e1331, 2018.
DOI:
https://doi.org/10.1002/wcms.1331.
URL https://onlinelibrary.wiley.com/
doi/abs/10.1002/wcms.1331.
[36] Per E.M. Siegbahn. A new direct ci method for large ci expansions in a small
orbital space. Chem. Phys. Lett., 109(5):417–423, 1984. ISSN 0009-2614. DOI:
https://doi.org/10.1016/0009-2614(84)80336-X. URL https://www.sciencedirect.
com/science/article/pii/000926148480336X.
[37] Mikhail Smelyanskiy, Nicolas PD Sawaya, and Alán Aspuru-Guzik. qhipster:
The quantum high performance software testing environment. arXiv preprint
arXiv:1601.07195, 2016. URL https://arxiv.org/abs/1601.07195.
[38] Daniel GA Smith, Lori A Burns, Andrew C Simmonett, Robert M Parrish, Matthew C
Schieber, Raimondas Galvelis, Peter Kraus, Holger Kruse, Roberto Di Remigio, AsemAlenaizan, et al. Psi4 1.4: Open-source software for high-throughput quantum chemistry. J. Chem. Phys., 152(18):184108, 2020. DOI: https://doi.org/10.1063/5.0006002.
URL https://aip.scitation.org/doi/10.1063/5.0006002.
[39] Nicholas H Stair and Francesco A Evangelista. Qforte: an efficient state simulator
and quantum algorithms library for molecular electronic structure. arXiv preprint
arXiv:2108.04413, 2021. URL https://arxiv.org/abs/2108.04413.
[40] Qiming Sun, Timothy C Berkelbach, Nick S Blunt, George H Booth, Sheng Guo,
Zhendong Li, Junzi Liu, James D McClain, Elvira R Sayfutyarova, Sandeep Sharma,
et al. Pyscf: the python-based simulations of chemistry framework. Wiley Interdiscip.
Rev. Comput. Mol. Sci., 8(1):e1340, 2018. DOI: https://doi.org/10.1063/5.0006002.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/wcms.1340.
[41] Yasunari Suzuki, Yoshiaki Kawase, Yuya Masumura, Yuria Hiraga, Masahiro Nakadai,
Jiabao Chen, Ken M Nakanishi, Kosuke Mitarai, Ryosuke Imai, Shiro Tamiya, et al.
Qulacs: a fast and versatile quantum circuit simulator for research purpose. arXiv
preprint arXiv:2011.13524, 2020. URL https://arxiv.org/abs/2011.13524.
[42] Barbara M. Terhal and David P. DiVincenzo. Classical simulation of noninteractingfermion quantum circuits. Phys. Rev. A, 65:032325, Mar 2002. DOI: 10.1103/PhysRevA.65.032325. URL https://link.aps.org/doi/10.1103/PhysRevA.65.032325.

A Evolution under a Hermitian Hamiltonian generated by a sum of excitation operators
Here we discuss the unitary evolution associated with a fermionic excitation (Eq. (2)) thatis not diagonal and has no repeated orbital indices. In this case, one can trivially showthat ĝĝ = 0 and ĝ † ĝ † = 0. Using these relations, it can be shown that(ĝ + ĝ † )n = ĝĝ † ĝ · · · + ĝ † ĝĝ † · · · .
Accepted in Quantum 2021-09-21, click title to verify. Published under CC-BY 4.0.

(31)
16

In addition, ĝ † ĝ and ĝĝ † are both diagonal in the computational basis, because any of theunpaired operators in ĝ would be paired with its conjugate in ĝ † ĝ and ĝĝ † . For example,
let ĝ = gâ†4 â†2 â3 â1 (we have omitted the spin indices for brevity), thenĝĝ † = |g|2 â†4 â†2 â3 â1 â†1 â†3 â2 â4= |g|2 n̂4 n̂2 (1 − n̂3 )(1 − n̂1 )(32)

It is worth noting that the parity associated with this reordering for normal-ordered g(all creation operators to the left of annihilation operators) is always even. We define thesquare root of this diagonal operator with the following phase convention,qĝĝ † ≡ |g|n̂4 n̂2 (1 − n̂3 )(1 − n̂1 )(33)

where we used the fact that the action of the number operators gives 0 or 1, and therefore,
they are idempotent. Using these expressions, the Taylor expansion of the evolution operator is exactly re-summed (where the summations are performed separately for the odd
and even rank terms) as†

e−i(g+g ) =∞X(−i)n

n!

n=0

=−1+

∞
X

(g + g † )n
h

(−1)n

n=0

−i

∞
X

(34)
i

(ĝĝ † )n + (ĝ † ĝ)n 2n
(2n)!

h

(−1)n

n=0

q

i

ĝ(ĝ † ĝ)n + ĝ † (ĝĝ † )n 2n+1(2n + 1)!
q

= − 1 + cos( ĝĝ † ) + cos( ĝ † ĝ)qq
11
sin( ĝĝ † ) − iĝ psin( ĝ † ĝ).
− ig † pĝĝ †ĝ † ĝDenoting the projection to the set of determinants that are not annihilated by an operatorx̂ as P̂x , together with the convention in Eq. (33), this can be further simplified to Eq. (3).

B Basis change: Evolution under quadratic Hamiltonians
In this work, we make use of the following primitive. Given an orbital basis {|φi} and manyPelectron wave function |Ψi = I(φ) CI(φ) |I(φ)i and a linear transformation X̂|φi → |φ0 i,
P
we wish to re-express |Ψi = I(φ0 ) CI(φ0 ) |I(φ0 )i where |I(φ0 )i is a determinant in the neworbital basis {|φ0 i}. There have been several discussions of how to implement this transformation efficiently [5, 25, 28]. In the original work by Malmqvist [25], it was understood thatthe transformation can be performed by the successive application of one-body operators
to the wave functions, in which the operator was computed from the LU decomposition of
the orbital transformation matrix X (Xij = hφi |φ0j i). Mitrushchenkov and Werner (MW)[28] later reported that the pivoting in the LU decomposition is necessary to make thetransformation stable. In their work, a LAPACK [2] function was used to perform the LU
decomposition with pivoting on the rows of the orbital transformation matrix. This leads
to reordering of the orbitals in the determinants and additional phase factors in CI(φ0 ) thatmake it complicated to, for example, evaluate the overlap and expectation values. Therefore, we have implemented a column-pivoted formulation of the MW algorithm, which is
Accepted in Quantum 2021-09-21, click title to verify. Published under CC-BY 4.0.

17

presented below. We will show the spin-free formulation (i.e. same transformation for α
and β orbitals) for spin-conserving wave functions as an example, but the procedure forthe spin-broken case can be derived in the same way.
First, we obtain the column pivoted LU decomposition. This is done using numpy.linalg.lu(which pivots rows) for the transpose of the orbital transformation matrix X,
XT = P̄L̄Ū,(35)

where L̄ and Ū are lower- and upper-triangular matrices, and the diagonal elements of L̄
are unit. It is then easily seen that, by taking the transpose, one obtains
X = ŪT L̄T P̄T ,(36)

where ŪT and L̄T are lower and upper triangular, respectively. If one scales the rows andcolumns of ŪT and L̄T , respectively, such that the diagonal elements of ŪT become unit,
this can be rewritten as
X = LUP,(37)

in which L and U are lower- and upper-triangular matrices with the diagonal elementsof L being unit, and P = P̄T . These are done in the ludecomp and transpose_matrix
subfunctions in the wavefunction.transform function.
When pivoting is not considered (i.e., P = 1), wave functions are transformed as follows(see Ref. [25] for details). Suppose that the spin-free MO coefficients are rotated asC̃ri =X

Crj (LU)ji

(38)

j

Using L and U, we compute the matrix elements
F = U−1 − L − I(39)which is performed in the process_matrix function. It has been shown [25] that thisoperator can be used to transform the many-body wave functions. Let T̂ be the oneparticle operator associated with this transformation, i.e.,|Ψ0 i = T̂ (LU)|Ψi.

(40)|Ψk+1,σ i = (1 + F̂kσ )|Ψk,σ i(41)

We perform

in which F̂kσ is a one-body, spin-orbital operator whose matrix elements are the kth columnof F, i.e.,F̂kσ =

X

Fik â†iσ âkσ

(42)

i

This operation is performed recursively first for σ = α from k = 0 (|Ψ0,α i = |Ψi) to k =n − 1 where n = dim(X) using a specialized apply function for "one-column" one-particle
operator, followed by the same procedure for σ = β from k = 0 (|Ψ0,β i = |Ψn,α i) to k =n−1. The resulting wave function, |Ψ0 i = |Ψn,β i, then has the same determinant expansion
coefficients CI(φ) as the coefficients of |Ψi (CI(φ0 ) ) when expressed in the determinants of
the new orbital basis.
Accepted in Quantum 2021-09-21, click title to verify. Published under CC-BY 4.0.

18

In the context of exact evolution with a quadratic Hermitian operator, the orbital transformation is performed to diagonalize the quadratic operator. Assume that the operator
PP
is Â = ij Aij Êij with Êij = σ a†iσ ajσ ; we diagonalize it so that X† AX = a where a isa diagonal matrix. Note that X is unitary. Using the column-pivoted LU decomposition(Eq. 37), the propagation can be performed asexp(−iÂt)|Φi = T̂ ((LU)† ) exp(−iâ0 t)T̂ (LU)|Φi(43)

where â0 is a diagonal operator,
â0 =

X
i

(PaPT )ii

X †âiσ âiσ(44)σ

To summarize, in this column-pivoted formulation, the pivot matrix P is used only to reorder the orbitals in the diagonal operator. This is more trivial and efficient than reordering
the orbitals in the wave functions.

Accepted in Quantum 2021-09-21, click title to verify. Published under CC-BY 4.0.

19