Tapering off qubits to simulate fermionic Hamiltonians

# Tapering off qubits to simulate fermionic Hamiltonians

Sergey Bravyi IBM T.J. Watson Research Center, Yorktown Heights, NY 10598, USA    Jay M. Gambetta IBM T.J. Watson Research Center, Yorktown Heights, NY 10598, USA    Antonio Mezzacapo IBM T.J. Watson Research Center, Yorktown Heights, NY 10598, USA    Kristan Temme IBM T.J. Watson Research Center, Yorktown Heights, NY 10598, USA
July 18, 2019
###### Abstract

We discuss encodings of fermionic many-body systems by qubits in the presence of symmetries. Such encodings eliminate redundant degrees of freedom in a way that preserves a simple structure of the system Hamiltonian enabling quantum simulations with fewer qubits. First we consider symmetry describing the particle number conservation. Using a previously known encoding based on the first quantization method a system of fermi modes with particles can be simulated on a quantum computer with qubits. We propose a new version of this encoding tailored to variational quantum algorithms. Also we show how to improve sparsity of the simulator Hamiltonian using orthogonal arrays. Next we consider encodings based on the second quantization method. It is shown that encodings with a given filling fraction and a qubit-per-mode ratio can be constructed from efficiently decodable classical LDPC codes with the relative distance and the encoding rate . A family of codes based on high-girth bipartite graphs is discussed. Graph-based encodings eliminate roughly qubits. Finally we consider symmetries, and show how to eliminate qubits using previously known encodings, illustrating the technique for simple molecular-type Hamiltonians.

## I Introduction

Quantum information processing holds the promise of solving some of the problems that are deemed too challenging for conventional computers. One important problem in this category is the simulation of strongly interacting fermionic systems, in the context of quantum chemistry or material science. A natural application for a quantum computer would be preparing low-energy states and estimating the ground energy of a fermionic Hamiltonian. Several methods have been proposed in the literature to accomplish this task, for example, preparation of a good trial state followed by the quantum phase estimation AbramsLloyd1 (); aspuru2005simulated () or state preparation by the adiabatic evolution farhi2000quantum (); babbush2014adiabatic (). These methods however require a universal quantum computer capable of implementing very long circuits, exceeding the state-of-the-art demonstrations by many orders of magnitude wecker2015solving (). Alternative methods that could be more viable in the near future are variational quantum algorithms peruzzo2014variational (); wecker2015progress (); mcclean2016theory (); Li2016efficient (). Such algorithms minimize the energy of a target fermionic Hamiltonian over a class of trial states that can be prepared on a given quantum hardware by varying control parameters.

Since the basic units of a quantum computer are qubits rather than fermi modes, any simulation method relies on a certain encoding of fermionic degrees of freedom by qubits Wigner1928 (); BK2002 (); verstraete2005mapping (); seeley2012bravyi (); Tranter15 (); Moll2015 (); whitfield2016local (). A choice of a good encoding is important as it may affect both the number of qubits and the running time of a simulation algorithm. Here we propose encodings tailored to variational quantum algorithms and fermionic Hamiltonians that possess symmetries such as the particle number conservation. The presence of symmetries allows one to restrict the simulation to an eigenspace of the symmetry operator whereby reducing the number of qubits that encode a fermionic system. In several specific examples, such as the hydrogen molecule and the Fermi-Hubbdard model it has been observed that some qubits can indeed be removed from the simulation without loss of information o2015scalable (); Moll2015 (). Importantly, the removal of qubits in these examples preserves a simple structure of the encoded Hamiltonian enabling efficient simulations with fewer qubits. This motivates the question of how to generalize these methods and how to eliminate redundant qubits in a computationally efficient manner for larger systems.

To address these questions let us first give a more precise notion of a simulation. We shall describe a fermionic system to be simulated by a target Hamiltonian composed of one- and two-body operators such as those describing hopping, chemical potential, and two-particle interactions,

 Htgt=M∑α,β=1tαβa†αaβ+M∑α,β,γ,δ=1uαβγδa†αa†βaγaδ. (1)

Here is the number of fermi modes, and are creation and annihilation operators for a mode , and are complex coefficients such that and . Leaving aside superconductivity and relativistic effects, all natural fermionic Hamiltonians have the above form. Since each term in has equal number of creation and annihilation operators, commutes with the particle number operator . We shall assume that the system contains a fixed number of particles . For example, could be the number of valence electrons in a molecule. Define a target Hilbert space as the -particle subspace spanned by all states of fermi modes satisfying . Without loss of generality, (otherwise consider holes instead of particles). Our goal is to estimate the minimum energy of restricted to the -particle subspace . Below we shall often identify and the restriction of onto the subspace . We shall choose the energy scale such that all coefficients in Eq. (1) have magnitude at most one.

Let us now formally define an encoding of fermionic degrees of freedom by qubits. Following Ref. TIM2014 (), we shall describe such encoding as an isometry , where is the Hilbert space of qubits. A state of the target system is identified with a state of the simulator system. Encoded states span a codespace .

A Hamiltonian describing a system of qubits is called a simulator of if it satisfies two conditions. First, we require that

 HsimE=EHtgt. (2)

In words, must preserve the codespace and the restriction of onto the codespace must be unitarily equivalent to . The action of on the orthogonal complement of the codespace may be arbitrary. Secondly, we require that the codespace contains a ground state of . This guarantees that and have the same ground energy and their ground states coincide modulo the encoding .

Our goal is to construct encodings that require qubits and, at the same time, all target Hamiltonians Eq. (1) have sufficiently simple simulators. We shall say that a simulator Hamiltonian is -sparse if

 Hsim=r∑i=1Di, (3)

where each operator is diagonal in some tensor product basis of qubits. This basis may depend on . We require that matrix elements of are efficiently computable. A family of encodings as above is called sparse if there exist small constants such that any target Hamiltonian Eq. (1) has an -sparse simulator Eq. (3) with and for all . Sparse encodings are well-suited for applications in variational quantum algorithms peruzzo2014variational (); wecker2015progress (); mcclean2016theory (); Li2016efficient (). Indeed, a basic subroutine of variational algorithms is estimating energy of a given trial state that can be prepared on the available quantum hardware. One can estimate the expectation value by preparing the trial state , performing a local change of basis in each qubit such that becomes diagonal in the standard , basis, and then measuring each qubit. Performing a sequence of such measurements with a freshly prepared trial state for each term gives an estimate of the energy .

### Summary of results

First assume that the number of particles is small such that . We expect that this regime may be relevant for high-accuracy quantum chemistry calculations with large basis sets. We construct a sparse encoding with qubits such that any target Hamiltonian Eq. (1) has a simulator Eq. (3) with sparsity

 r≤9M3.17. (4)

We also give a upper bound on the norm of the terms although we do not expect this bound to be tight. This encoding mostly follows ideas of Refs. Toloui2013 (); Babbush2015exponentially () and relies on the first-quantization method. We extend the results of Refs. Toloui2013 (); Babbush2015exponentially () in two respects. First we show how to improve the sparsity of using orthogonal arrays Hedayat2012orthogonal (). Such arrays have been previously used for quantum simulations and dynamical decoupling Wocjan2002 (); Kern2005 (); Rotteler2006 () but their application in the context of variational quantum algorithms appears to be new. Secondly, we introduce a sparse Hamiltonian enforcing the anti-symmetric structure of encoded states and compute the spectral gap of this Hamiltonian using arguments based on the Schur duality Goodman2000representations (); bacon2005quantum (); Bacon2006efficient (). This allows us to bound the norm .

Next consider encodings based on the second quantization method. Let be the filling fraction and be the desired qubit-per-mode ratio. We show that sparse encodings with prescribed , can be constructed from classical error correcting codes with certain properties. Namely, we need binary linear codes that have a column-sparse parity check matrix, relative distance , and the encoding rate . It is known that the requisite codes can be constructed whenever

 ν<1/4andη>h(2ν), (5)

where is the binary Shannon entropy. For example, one can use good LDPC codes Gallager1962 () whose parameters approach the Gilbert-Varshamov bound Varshamov1957 (); MacWilliamsBook (). Thus a constant fraction of qubits can be eliminated if the target system has a filling fraction . The simulator Hamiltonian Eq. (3) has sparsity proportional to the number of non-zero coefficients , in the target Hamiltonian. Furthermore, for all .

The bound Eq. (5) is worse than what one could expect naively. Indeed, since has dimension , the information-theoretic bound on the qubit-per-mode ratio is . We leave as an open question whether sparse encodings can achieve this bound.

The above result has one important caveat. Namely, we shall see that matrix elements of the simulator Hamiltonian can be computed efficiently only if the chosen code is efficiently decodable. In Appendix C we describe a brute force implementation of the decoding algorithm that may be practical for small number of modes . Simulating larger systems may require LDPC codes that are both good and efficiently decodable. Designing such codes is an active research area, see Refs. sipser1996expander (); zemor2001expander (); guruswami2005linear ().

To enable efficient decoding and improve sparseness of the simulator Hamiltonian we consider a special class of LDPC codes associated with high-girth bipartite graphs. For such encodings any two-body operator has a -sparse simulator, while any four-body operator has -sparse simulator. Furthermore, matrix elements of the simulators can be computed in time . Graph-based encodings can eliminate qubits for . Figure 1 shows a numerically computed lower bound on the number of modes that can be simulated for fixed values of using the graph-based encodings.

Our third result concerns -symmetries such as the fermionic parity conservation. For concreteness, we consider a symmetry group describing parity conservation for electrons with a fixed spin orientation. Such symmetry is present, for example, in the standard molecular electron Hamiltonians that neglect spin-spin and spin-orbit interactions. The target Hilbert space is chosen as a subspace in which the number of electrons with a given spin is fixed modulo two. Using the second quantization method and both the parity and the binary tree encoding of Ref. BK2002 () we construct a simulator Hamiltonian Eq. (3) that acts on qubits. More generally, we give a systematic method of detecting symmetries in a given target Hamiltonian and show how to construct sparse encodings that eliminate one qubit for each symmetry. Our techniques are illustrated using quantum chemistry Hamiltonians describing simple molecules.

To summarize, we observed that the encoding based on the first quantization method achieves the best performance in terms of the number of qubits that can be eliminated. However, it is applicable only if the number of particles is relatively small. Furthermore, the encoding does not take advantage of any structure present in the coefficients of . For example, might have sparseness even if has only non-zero coefficients. In contrast, encodings based on the second quantization method eliminate fewer qubits but have broader applicability and produce more sparse simulator Hamiltonian such that the number of terms in and are roughly the same (up to a constant factor). Which encoding should be preferred may depend on details of the target system.

We expect our results to find applications in the near-future experimental demonstrations of variational quantum algorithms. In this context the number of qubits is fixed by the hardware constraints and may not be large enough to simulate interesting molecules directly. Combining the standard variational algorithms peruzzo2014variational (); wecker2015progress (); mcclean2016theory () with the encodings described in this paper may extend the range of molecules that can be simulated on a given quantum hardware. We leave as open questions whether sparse encoding can be constructed for the filling fraction , what is the tradeoff between the parameters and the encoding sparseness, and how to generalize our techniques to other types of symmetries.

The rest of the paper is organized as follows. Section II summarizes our notations. Encodings based on the first quantization method are described in Section III. These encoding are applicable if the number of particles is sufficiently small. Section IV shows how to construct sparse encodings for a constant filling fraction using classical LDPC codes. Encodings based on high-girth bipartite graphs are described in Sections V,VI. Discrete symmetries and applications of our techniques to small molecular-type Hamiltonians are discussed in Sections VII,VIII. Appendix A summarizes the previously known encodings of fermions by qubits. Appendix B illustrates our methods using the hydrogen molecule as an example. Appendix C shows how to compute matrix elements of simulator Hamiltonians constructed from LDPC codes.

## Ii Notations

A system of fermi modes is described by the Fock space of dimension equipped with the standard basis , where is the occupation number of the mode such that . Our target Hilbert space is defined as the -particle subspace of ,

 Htgt=span(|x⟩∈FM:|x|=N). (6)

Here denotes the Hamming weight of . The simulator system consists of qubits, where satisfies

 dim(Htgt)=(MN)≤2Q≤2M.

The Hilbert space is equipped with the standard basis , where . We shall reserve letters for qubit basis vectors and letters for the Fock basis vectors. For any integer let .

Suppose is a two-body or four-body fermionic observable (hermitian operator) listed below

 iϵ(a†αaβ±a†βaα),iϵ(a†αa†βaγaδ±a†δa†γaβaα), (7)

where is chosen to make the operator hermitian. We shall say that a qubit observable acting on is a simulator of if

 OsimE=EOtgt. (8)

A direct consequence of Eq. (8) is that preserves the codespace and the restriction of onto the codespace is unitarily equivalent to . The action of on the orthogonal complement of the codespace may be arbitrary. Let us say that the simulator is -sparse if it can be written as

 Osim=r∑i=1Di, (9)

where are hermitian operators such that each operator is diagonal in some tensor product basis of qubits. This basis may depend on . The maximum sparsity of two-body and four-body simulators will be denoted and respectively.

## Iii Sparse encodings for small number of particles

Here we discuss encodings based on the first quantization method. Assume for simplicity that the number of modes is a power of two, . Otherwise, round up to the nearest power of two. Given a fermi mode , let be the binary representation of the integer .

The simulator system consists of qubits partitioned into consecutive registers of qubits each. For any -tuple of modes let be a basis vector such that the register is in the state .

Consider a basis vector and let

 1≤α1<α2<…<αN≤M

be the subset of modes that are occupied in the state , that is, iff . Define the encoding as

 E|x⟩=1√N!∑π∈SN(−1)πPπ|¯¯¯¯α1,¯¯¯¯α2,⋯,¯¯¯¯αN⟩, (10)

where is the group of permutations of objects, is the sign of a permutation , and is a unitary operator that applies a permutation to the registers such that

 Pπ|¯¯¯¯α1,¯¯¯¯α2,⋯,¯¯¯¯αN⟩=|¯¯¯¯απ(1),¯¯¯¯απ(2),⋯,¯¯¯¯απ(N)⟩.

The righthand side of Eq. (10) can be viewed as the first-quantized version of the state . The codespace is spanned by antisymmetric states such that

 Pπ|ψ⟩=(−1)π|ψ⟩,for all π∈SN.

We choose the simulator Hamiltonian as

 Hsim=T+U+gH⊥, (11)

where is the first-quantized version of , namely

 T=N∑i=1M∑α,β=1tαβ|α⟩⟨β|i (12)
 U=−∑1≤i≠j≤NM∑α,β,γ,δ=1uαβγδ|α,β⟩⟨γ,δ|i,j (13)

Here and below the subscripts indicate the registers acted upon by an operator. Finally, penalizes states orthogonal to the codespace,

 H⊥=∑1≤i

Here is the SWAP operator that exchanges and . Note that has zero ground energy and its ground subspace coincides with the codespace . The coefficient in Eq. (11) will be chosen large enough so that the ground state of belongs to the codespace.

Note that for any permutation . Thus preserves the codespace . The standard correspondence between the first and the second quantized Hamiltonians implies that the restriction of onto the codespace is unitarily equivalent to , so that Eq. (2) is satisfied. Thus is indeed a simulator of (for large enough to be chosen later).

Let us show that is -sparse, where

 r=9m≤M3.17. (15)

Furthermore, , where each term is diagonal in a tensor product of Pauli bases

 X≡{(|0⟩±|1⟩)/√2},Y≡{(|0⟩±i|1⟩)√2},

and . Let

 Pauli(m)={σ1,σ2,…,σM2}

be the set of all -qubit Pauli operators (ignoring the overall phase). By definition, each operator is a tensor product of single-qubit Pauli operators . Note that there are such operators. We note that the SWAP operator on two qubits can be written as

 12(I⊗I+σx⊗σx+σy⊗σy+σz⊗σz).

Since the SWAP operator exchanging -qubit registers , is a tensor product of two-qubit SWAPS, one gets

 (↔)i,j=M−1M2∑a=1(σa⊗σa)i,j. (16)

Expanding each term in Eqs. (12,13) in the basis of Pauli operators and using Eq. (16) to expand one gets

 Hsim=∑1≤i

for some real coefficients . We shall group Pauli operators that appear in Eq. (17) into bins such that Pauli operators from the same bin are diagonal in the same tensor product basis. First consider a single register . Obviously, any Pauli operator acting on is diagonal in a tensor product of the bases . Such tensor product bases can be labeled by letters in the alphabet

 A≡{X,Y,Z}m. (18)

Recall that an orthogonal array Hedayat2012orthogonal () over an alphabet is a matrix of size with entries from such that any pair of columns of contains each two-letter word in the alphabet the same number of times. (More precisely, the above defines an orthogonal array with strength two). Orthogonal arrays have been previously used for quantum simulations and dynamical decoupling Wocjan2002 (); Kern2005 (); Rotteler2006 (). We shall use a family of orthogonal arrays based on the Galois field known as Rao-Hamming construction Hedayat2012orthogonal (). It gives an orthogonal array over an alphabet of size with and . Note that the equality is possible only if any pair of columns of contains each two-letter word in the alphabet exactly one time. Also note that the number of particles obeys . We shall use only the first columns of . Let be some row of . It defines a tensor product of Pauli bases

 Rf≡Rf,1⊗Rf,2⊗⋯⊗Rf,N

for the system of qubits. By construction, each Pauli term that appears in Eq. (17) is diagonal in at least one basis . Thus we can choose a decomposition where is diagonal in the basis . This shows that any target Hamiltonian Eq. (1) has a simulator with sparsity . Rounding up to the nearest power of two gives Eq. (5).

The coefficient in Eq. (11) must be large enough so that the ground state of belongs to the codespace. This is always the case if , where is the smallest non-zero eigenvalue of , see Eq. (14). Indeed, suppose is an eigenvector of orthogonal to the codespace. Then the corresponding eigenvalue is

 ⟨ψ|Hsim|ψ⟩≥⟨ψ|U+T|ψ⟩+gΔ⊥>∥U+T∥.

Such eigenvalue cannot be the ground energy of since the latter is unitarily equivalent to a submatrix of . Thus the ground state of must belong to the codespace. Recall that we assume the coefficients to have magnitude at most one. This gives a conservative estimate . Below we show that

 Δ⊥=N2for all N≤M. (19)

Thus it suffices to choose and all terms in the simulator Hamiltonian have norm . We do not expect the bound to be tight. In practical implementation of variational quantum algorithms one can start from and gradually increase until the best variational state satisfies .

Let us now prove Eq. (19). Recall that is the smallest non-zero eigenvalue of the Hamiltonain defined in Eq. (14). We shall use a symmetry-based argument to compute all eigenvalues of . Consider the permutation group and the unitary group . The Hilbert space defines a unitary representation of these groups such that a permutation acts on as and a unitary matrix acts on as . The standard result from the group representation theory known as Schur duality Goodman2000representations () gives a decomposition

 (\mathbbmCM)⊗N=⨁λPλ⊗Qλ, (20)

where the sum runs over all Young diagrams with boxes, is the irreducible representation (irrep) of the permutation group and is the irrep of the unitary group . The operators and are block-diagonal with respect to Schur decomposition Eq. (20). Furthermore, within each sector the operator acts non-trivially only on the subsystem , while acts non-trivially only on subsystem .

Importantly, the Hamiltonian commutes with the action of both groups and , that is,

 [H⊥,Pπ]=[H⊥,V⊗N]=0

for all and for all . Since the groups and generate the full operator algebra in each sector in the decomposition Eq. (20), we conclude that

 H⊥=⨁λeλΠλ, (21)

where is the projector onto the sector in Eq. (20) and are eigenvalues of . Thus we can compute by picking an arbitrary state from the sector and computing .

We shall consider a Young diagram with columns as a partition of the integer , that is,

 N=λ1+…+λd,λ1≥…≥λd≥1.

Namely, is the number of boxes in the -th column of . For any integer define a state such that

 |ϕ(u)⟩=1√u!∑π∈Su(−1)π|π(1),π(2),…,π(u)⟩.

For example, , ,

 |ϕ(3)⟩=1√6(|1,2,3⟩+|3,1,2⟩+|2,3,1⟩−|2,1,3⟩−|3,2,1⟩−|1,3,2⟩)

etc. We claim that a state

 |ψλ⟩≡|ϕ(λ1)⟩⊗⋯⊗|ϕ(λd)⟩ (22)

belongs to the sector of the decomposition Eq. (20). Namely, such state can be obtained by applying a suitable Young symmetrizer bacon2005quantum () to a basis vector. Indeed, define a Young tableau obtained by filling columns of one by one with consecutive integers starting from the first column, see Figure 2 for an example. Let and be the subgroups that permute integers from the same column and from the same row of respectively. The Young symmetrizer corresponding to tableau is defined as

 Πλ,T∼⎛⎝∑π∈Scol(−1)πPπ⎞⎠⋅(∑τ∈SrowPτ). (23)

It is well-known Goodman2000representations (); bacon2005quantum () that is proportional to a (non-orthogonal) projector onto a subspace of the sector in the Schur decomposition. In particular, maps any state to some state that belongs to the sector . Let be a sequence of integers obtained by filling columns of one by one with consecutive integers such that the -th column is filled with integers , see Figure 2 for an example. Let be the basis vector corresponding to .

We observe that the second factor in Eq. (23) has trivial action on since for all . It follows that

 Πλ,T|s(λ)⟩∼|ψλ⟩ (24)

and thus indeed belongs to the sector of the Schur decomposition. One can easily check that is an eigenvector of with the eigenvalue

 eλ=12[(N2)−d∑a=1(λa2)+∑1≤a

From Eq. (21) one infers that any eigenvalue of must have a form . Note that iff is a single column, that is, , . One can check that the smallest non-zero value of is achieved when has two columns with length and , that is, , , and . In this case which proves Eq. (19).

## Iv Sparse encoding based on LDPC codes

In this section we use the second quantization method and classical LDPC codes to construct sparse encodings for the case when the target system has a constant filling fraction .

Let be a binary matrix with rows and columns. Given a binary vector of length , we shall write for the matrix-vector multiplication modulo two, that is,

 (Ax)i=M∑α=1Ai,αxα(mod2) (26)

We consider encodings defined by

 E|x⟩=|Ax⟩ (27)

where and . Let us say that a matrix is -injective if it maps distinct -bit vectors with the Hamming weight to distinct -bit vectors . It follows directly from the definitions that is an isometry iff is -injective. The -injectivity condition is satisfied if is chosen as a parity check matrix describing a binary linear code of length with the minimum distance . Indeed, in this case all errors of weight up to must have different syndromes and thus a syndrome uniquely identifies a weight- error .

How sparse is the encoding defined in Eq. (27)? Let columns of be and let

 c(A)=maxα|Aα|

be the maximum column weight. We claim that fermionic observables defined in Eq. (7) have simulators Eq. (9) with sparsity

 r2≤22c(A)−1andr4≤24c(A)−1 (28)

for two-body and four-body observables respectively. Furthermore, for all . Thus the encoding defined by Eq. (27) is sparse whenever is a column-sparse matrix.

First, let us introduce some notations. Let be a set of all weight- vectors satisfying ,

 A−1s≡{x∈{0,1}M:Ax=sand|x|=N}. (29)

The set may be empty for some . By definition, is -injective iff the set contains at most one element for any . Below denotes a string with a single non-zero at the position . We use the notation for the bitwise XOR.

For concreteness, consider a pair of modes and a target observable

 Otgt=a†αaβ+a†βaα. (30)

We have if and

 Otgt|x⟩=Sαβ(x)|x⊕eα⊕eβ⟩

if where is the parity of all bits of located between and , that is,

 Sαβ(x)=β−1∏γ=α+1(−1)xγ.

Since , we have

 EOtgt|x⟩ = Sαβ(x)|Ax⊕Aα⊕Aβ⟩if xαxβ=01,10, EOtgt|x⟩ = 0if xαxβ=00,11. (31)

Let us say that a basis vector is -flippable if for some weight- string such that or . Note that is a single string whenever is -flippable. Define an operator acting on such that

 Γαβ|s⟩ = Sαβ(A−1s)|s⟩if s is αβ-flippable Γαβ|s⟩ = 0otherwise. (32)

Given a bit string , let be the product of Pauli operators over all qubits such that . We claim that the observable has a simulator

 Osim=X(Aα⊕Aβ)Γαβ=ΓαβX(Aα⊕Aβ). (33)

First let us check that commutes with . Suppose is -flippable and let . By assumption, for some such that and, say, . It follows that has weight and . Furthermore, . Thus is -flippable and . Since , we have shown that and thus

 X(Aα⊕Aβ)Γαβ|s⟩=Sαβ(x)|t⟩=ΓαβX(Aα⊕Aβ)|s⟩.

If is not -flippable then so is , so that

 X(Aα⊕Aβ)Γαβ|s⟩=ΓαβX(Aα⊕Aβ)|s⟩=0.

We have shown that commutes with .

Next, let us check the simulation condition Eq. (2). Suppose has weight and let . Using the first equality in Eq. (33) one infers that

 OsimE|x⟩=X(Aα⊕Aβ)Γαβ|s⟩=Sαβ(x)|s⊕Aα⊕Aβ⟩

if and otherwise. Comparing this and Eq. (31) shows that .

Let us show that is -sparse with . By construction, is diagonal in the standard basis, that is,

 Γαβ=∑s∈{0,1}Mg(s)|s⟩⟨s|,g(s)=0,±1. (34)

Let . To simplify notations, let us reorder the qubits such that acts on the first qubits. Decompose , where and . Define a function obtained from by applying the Walsh-Hadamard transform with respect to the first argument:

 h(t,s′)≡2−k∑u∈{0,1}k(−1)t⋅ug(u,s′). (35)

Here . Substituting the identity

 |s⟩⟨s|≡|u,s′⟩⟨u,s′|=2−k∑t∈{0,1}k(−1)t⋅uZ(t)⊗|s′⟩⟨s′|

into Eq. (34) gives

 Γαβ=∑t∈{0,1}k∑s′∈{0,1}Q−kh(t,s′)Z(t)⊗|s′⟩⟨s′|. (36)

For each define an operator

 Γαβ(t)=∑s∈{0,1}Q−kh(t,s)|s⟩⟨s|. (37)

acting on the last qubits. Then

 Γαβ=∑t∈{0,1}kZ(t)⊗Γαβ(t). (38)

As was shown above, commutes with . Since commutes (anti-commutes) with for even (odd) , we infer from Eq. (38) that whenever has odd weight. Combining Eqs. (33,38) we arrive at

 Osim=∑t∈{0,1}k|t|evenX⊗kZ(t)⊗Γαβ(t), (39)

This gives a -sparse decomposition of as defined in Eq. (9) where

 Di≡X⊗kZ(t)⊗Γαβ(t). (40)

This operator can be made diagonal in the standard basis by applying a Clifford operator exchanging Pauli and for all qubits such that . It remains to note that

 k=|Aα⊕Aβ|≤|Aα|+|Aβ|≤2c(A).

Thus the simulator Eq. (39) has sparsity . Furthermore, all matrix elements of are contained in the interval , see Eqs. (34,35,37,40). We omit the derivation of simulators for other observables defined in Eq. (7) since it follows exactly the same steps as above.

Consider a target Hamiltonian defined in Eq. (1). Decompose as a linear combination of two-body and four-body observables defined in Eq. (7). Replacing each observable by a qubit simulator constructed above gives a simulator Hamiltonian

 Hsim=gEE†+r∑i=1Di. (41)

Here we combined the terms from each simulator into a single sum. The term penalizes states orthogonal to the codespace. Note that is upper bounded by a constant times the number of non-zero coefficients in the target Hamiltonian. Since , we can guarantee that the ground state of belongs to the codespace provided that .

Let us choose as a parity check matrix of a binary linear code that encodes bits into bits with the minimum distance . As was argued above, such matrix is -injective. It is known that certain families of codes described by sparse parity check matrices can approach the Gilbert-Varshamov bound Varshamov1957 (); MacWilliamsBook (), namely,

 K=M(1−h(2N/M)−ϵ), (42)

where is the binary Shannon entropy function and can be made arbitrarily small by choosing large enough . This claim follows from the existence of good LDPC codes Gallager1962journal (), see for instance Theorem A.3 of Gallager1962 (). We can assume wlog that all rows of are linearly independent in which case has rows. Then a family of good LDPC codes as above gives a family of sparse encodings with the filling fraction and the qubit-per-mode ratio that satisfy , as claimed in Eq. (5). Unfortunately, the constant grows quickly as approaches , see Gallager1962 (). Since the sparsity of simulators constructed for few-body fermionic observables is exponential in , see Eq. (28), encodings based on good LDPC codes are not quite practical. We show how to overcome this problem using “bad” LDPC codes in Sections V,VI.

Next let us discuss how to compute matrix elements of the simulator Hamiltonian Eq. (41). Note that all steps in the definition of are computationally efficient except for inverting the action of , that is, computing the set defined in Eq. (29). Define a function

 xmin(s)=argminx∈{0,1}MAx=s|x|. (43)

It returns an error of minimum weight consistent with a given syndrome . Suppose is a parity check matrix of a linear code with the minimum distance . It follows easily from the definitions that if has weight and otherwise. Thus it suffices to give an efficient algorithm for computing . The latter is known as a minimum weight decoding problem. Although in general this problem is NP-hard Berlekamp1978 (), there are special classes of LDPC codes that admit a linear time decoder sipser1996expander (); zemor2001expander (). These codes have a non-zero encoding rate and relative distance, but they are not good in the sense of Eq. (42). In Section VI we discuss a special class of LDPC codes based on high-girth bipartite graphs that can be decoded in time . Appendix C gives a simple algorithm that computes the set for any -injective matrix. Although this algorithm is not efficient asymptotically, it can be implemented for small system sizes .

## V Improving the sparsity

Here we show how to improve the sparsity bounds in Eq. (28) if the parity check matrix has a certain additional structure. At this point we shall exploit the fact that simulators only need to reproduce the action of target observables within the codespace and can act arbitrarily on the orthogonal complement to the codespace.

Let be a binary matrix of size with columns . We shall say that is bipartite if the set of rows can be partitioned into two disjoint subsets, , such that each column intersects both and on odd number of rows,

 |Aα∩L|(mod2)=|Aα∩R|(mod2)