Operator Locality in Quantum Simulation of Fermionic Models

TL;DR

This paper reviews and analyzes various fermion-to-qubit mappings for quantum simulation, focusing on locality improvements and encoding efficiency, especially for lattice models like the Hubbard model.

Contribution

It reformulates the Bravyi-Kitaev transform with a classical data structure and introduces locality improvements for fermionic models, comparing different mappings' efficiencies.

Findings

Bravyi-Kitaev transform can be reformulated for better locality.

Auxiliary Fermion scheme achieves the most compact encoding with ancillas.

A variant of Bravyi-Kitaev is optimal without ancillas for Hubbard models.

Abstract

Simulating fermionic lattice models with qubits requires mapping fermionic degrees of freedom to qubits. The simplest method for this task, the Jordan-Wigner transformation, yields strings of Pauli operators acting on an extensive number of qubits. This overhead can be a hindrance to implementation of qubit-based quantum simulators, especially in the analog context. Here we thus review and analyze alternative fermion-to-qubit mappings, including the two approaches by Bravyi and Kitaev and the Auxiliary Fermion transformation. The Bravyi-Kitaev transform is reformulated in terms of a classical data structure and generalized to achieve a further locality improvement for local fermionic models on a rectangular lattice. We conclude that the most compact encoding of the fermionic operators can be done using ancilla qubits with the Auxiliary Fermion scheme. Without introducing ancillas, a variant of the Bravyi-Kitaev transform provides the most compact fermion-to-qubit mapping for Hubbard-like models.

Tables

Table 1. Table 1: Operator locality and qubit resource overheads for transformation of a 2 2 2 D rectangular lattice Hubbard model.

Method	Density-density	Horizontal	Vertical	Qubits
JW	2	2	w+1	$2wh$
BK	$2\lfloor {\log }_2\left(wh\right)\rfloor +2$	$\lfloor {\log }_2\left(wh\right)\rfloor +\lceil {\log }_2\left(wh\right)\rceil$	$\lfloor {\log }_2\left(wh\right)\rfloor +\lceil {\log }_2\left(wh\right)\rceil$	$2wh$
SBK	$2\lfloor {\log }_{2}w\rfloor +2$	$\lfloor {\log }_{2}w\rfloor +\lceil {\log }_{2}w\rceil$	$2\lfloor {\log }_{2}w\rfloor +1$	$2wh$
AF	2	2	4	$4\left(wh-1\right)$
LSFS	8	7	7	$4wh-2h-2w$

Table 2. Table 2: Worst-case locality of the hopping term as a function of the hypercubic lattice dimension.

Method	Worst-case locality	Qubits
JW	$w^{D-1}+1$	$2w^D$
BK	$2\lfloor {\log }_2{w}^D\rfloor$	$2w^D$
SBK	$2\lfloor {\log }_2{w}^{D-1}\rfloor +1$	$2w^D$
AF	$2D$	$2D{w}^D$
LSFS	$4D$	$2D\left(w-1\right)w^{D-1}$

Equations

$$H=-t\sum _{(i,j)\in E}\sum _{\sigma =\uparrow ,\downarrow }(a_{i\sigma }^{†}{a}_{j\sigma }+a_{j\sigma }^{†}{a}_{i\sigma })+U\sum _{i\in V}{n}_{i\uparrow }{n}_{i\downarrow },$$

$$a_k$$

$$\prod _{j=0}^{k-1}{Z}_j=\underset{k}{\underbrace{Z\otimes Z\dots Z}}\otimes \underset{N-k}{\underbrace{1\otimes \dots 1}},$$

$$x_0$$

$$x_3$$

$$a_2^{†}$$

$$a_3^{†}\to$$

$$+\left(|0⟩{⟨0|}_1|0⟩{⟨0|}_2+|1⟩{⟨1|}_1|1⟩{⟨1|}_2\right)|1⟩{⟨0|}_3{X}_6.$$

$$c_3=a_3^{†}+a_3$$

$$c_j=a_j+a_j^{†}$$

$$d_j=i\left(a_j^{†}-a_j\right)$$

$$P(j)$$

$$B_k$$

$$A_{(jk)}$$

$$a_k^{†}{a}_j+a_j^{†}{a}_k$$

$$A_{(jk)}{B}_l$$

$$A_{(jk)}{A}_{(ls)}$$

$$(i{)}^p{A}_{(j_0{j}_1)}{A}_{(j_1{j}_2)}\dots A_{(j_p{j}_0)}$$

$$\prod _{k\in V}{B}_k$$

$${\tilde{B}}_k=\prod _{j\in n(k)}{Z}_{(jk)},$$

$${\tilde{A}}_{(jk)}$$

$${\tilde{A}}_{(jk)}={ϵ}_{jk}{X}_{(jk)}\prod _{l<k}^{n(j)}{Z}_{(lj)}\prod _{s<j}^{n(k)}{Z}_{(sk)}.$$

$${\tilde{B}}_k$$

$${\tilde{A}}_{(jk)}$$

$$A_{(\alpha \beta )}{A}_{(\beta \gamma )}{A}_{(\gamma \delta )}{A}_{(\delta \alpha )}$$

$$c_{\alpha }=\left(a_{\alpha }+a_{\alpha }^{†}\right),$$

$$C_{(\alpha \beta \gamma \delta )}$$

$$a_{k+1}^{†}{a}_k+a_k^{†}{a}_{k+1}$$

$$a_j^{†}{a}_k+a_k^{†}{a}_j$$

$$H_{\uparrow }$$

$$H_{\downarrow }$$

$$H$$

$$n_{k\uparrow }{n}_{k\downarrow }$$

$$n_j$$

Full text

I Introduction

Among the various applications of quantum computing, quantum simulation has long stood out as a primary motivation Feynman ; Lloyd . Classical computers can often perform rapid electronic structure calculations without explicit electron-electron interaction and obtain relatively accurate results Ramakrishnan14 . However, systems where the electron-electron interaction cannot be integrated out are called strongly correlated and represent a new frontier for electronic structure in both theoretical chemistry Jiang15 and strongly correlated materials, such as high-temperature superconductors Altland2010Condensed . It is in this regime that quantum simulation is a promising route forward Wecker ; Nori ; Georgescu ; Dallaire16 ; Kreula16 .

Quantum simulation comes in two distinct flavors: digital and analog, are each subject to different mindset and constraints. In the digital context, the hardware is thought of as a universal quantum computer where an arbitrary quantum circuit can be implemented and used to approximate the system of interest Lloyd . Since high-quality qubits are required in this context, the simulation qubit count can be thought of as an important constraint.

An analog quantum simulator on the other hand approximates the system with another, easier to implement, control and measure. Such a simulator or emulator is usually tailored to a specific problem and it is therefore argued to be technologically more viable to build such a chip rather than a general purpose quantum computer Georgescu . These analog simulators are typically restricted to 2-qubit couplings and a limited set of global operations - examples being both the trapped ions Blatt or the superconducting qubits Georgescu .

Quantum simulation of strongly correlated fermionic systems has recently been a focus of algorithmic developments Bela16 ; Kreula16 ; Dallarie16 ; Wecker . Besides direct simulation Wecker , it has been pointed out that quantum simulation of a strongly correlated region can act as an impurity solver for dynamical mean field theories Bela16 ; Dallaire16 ; Kreula16 . In these recent investigations, the authors have chosen to use the Jordan-Wigner fermionic encoding scheme JordanWigner ; Whitfield11 ; Somma02 for their specific simulations. However, under this transformation, local fermionic operators become spin operators acting on an extensive number of qubits, which may be problematic especially in the analog context. This can be avoided using other encoding schemes and we contribute to the ongoing line of research by investigating various fermion-to-qubit mappings.

The Hubbard model has served as a paradigmatic example for strongly correlated problems Wecker ; Altland2010Condensed ; Georgescu . We will continue this trend and use the Hubbard model as testbed for our ideas. Its Hamiltonian, on a graph with edges $E$ and vertices $V$, is given by:

[TABLE]

where $t$ and $U$ are parameters of the model, $n_{j\sigma}=a_{j\sigma}^{\dagger}a_{j\sigma}$, and the fermion creation operators $\{a_{i\sigma}\}$ satisfy $a_{i\sigma}^{\dagger}a_{j\tau}+a_{j\tau}a^{\dagger}_{i\sigma}=\delta_{ij}\delta_{\sigma\tau}$. We will consider $t$ and $U$ to be fixed and assume, for now, that we are on a square lattice.

The paper is organized into two parts. The first part reviews and extends mappings from fermionic Hamiltonians to qubits and the second part studies operator locality of the various fermion encoding methods.

II Mapping Fermionic Hamiltonians to Qubits

The following section reviews and expands on a set of locality improving transformations for mapping fermionic Hamiltonians to qubits. Subsection II.1 briefly summarizes the Jordan-Wigner transformation which will be used as a baseline for locality overhead comparison. In subsection II.2, we review the first method originally outlined in Ref. Kitaev which has been referred to as Bravyi-Kitaev transformation in the literature Seeley ; Tranter . We reformulate the transform in terms of a classical data structure (different perspective on the construction can be found in Seeley ; Tranter ). This allows for its generalization outlined in subsection II.3. The generalized Bravyi-Kitaev transformation corresponds to a whole class of fermion-to-spin transformations characterized by a transition from linear to logarithmic operator locality. In part II.5, we review the second method outlined in Ref. Kitaev and provide an example of $2$D Hubbard model mapping. Lastly, subsection II.6 reviews the Auxiliary Fermion method introduced in Ball ; Verstraete with focus on operator locality analysis. We have previously worked out construction details of this transformation in Ref. us .

II.1 Jordan-Wigner Transform

The usual way to map fermionic operators to qubits is the Jordan-Wigner (JW) transformation JordanWigner . This encoding stores information about the occupancy of $N$ fermionic sites in $N$ qubits. The fermion raising/lowering operators on $k$-th site are mapped to qubit operators by:

[TABLE]

where $Z_{j}$ stands for an $N$-qubit operator corresponding to a Pauli $Z$ operator applied to the $j$-th qubit and $\mathbf{1}$ to the rest of the qubit register. The above operators obey fermionic anti-commutation relations and therefore generate fermionic algebra on qubits. The

[TABLE]

string of Pauli $Z$ operators counts the excitation parity.

The action of raising/lowering operators can be therefore thought of as a composition of two operations on qubit states: (1) counting the parity and (2) updating the fermionic site occupancy. The number of single-qubit Pauli $Z$ operators used for parity counting scales asymptotically as $O(N)$, while the update is implemented with a single qubit operator $\ket{0}\bra{1}_{k}$ or its conjugate. The composite operation hence costs $O(N)$ in operator locality.

The fermionic raising/lowering operators occur only in pairs in any physical Hamiltonian. The Pauli strings could therefore cancel, as is the case in $1$D Hubbard model. In a general case however (specifically for a Hubbard model on higher dimensional lattices), the hopping operator locality scales with the size of the lattice. We therefore proceed by introducing an alternative scheme which improves locality of the resulting qubit Hamiltonian.

II.2 The Bravyi-Kitaev Transform

The first of the two fermionic transformations introduced in Kitaev , the Bravyi-Kitaev (BK) transform, can be described by a classical data structure, the Fenwick tree Fenwick94anew , which we will introduce below. The BK transform has been previously reviewed in Tranter ; Seeley and formulated in terms of recursive prescription for transformation matrices. This carried an implicit constraint on the number of qubits being a power of $2$. Our approach defines the scheme for an arbitrary number of qubits.

II.2.1 Fenwick Trees

In context of classical computation, a Fenwick tree can be used to map binary strings $n_{0}\,n_{1}\ldots n_{N}$, $n_{i}\in\{0,1\}$ to binary strings $x_{0}\,x_{1}\ldots x_{N}$, $x_{i}\in\{0,1\}$ such that both the prefix sum $\left(\sum_{m=0}^{k-1}n_{m}\right)$ and bit-flip operations have $O(\log N)$ access costs in the encoded representation. This optimization is achieved by storing partial occupancy sums ($x_{i}$) rather than occupancies/bits ($n_{i}$) in a way we now describe. The partial occupancy sums $x_{i}$ are dictated by the tree constructed using Algorithm 1:

A tree generated by $\textbf{Fenwick}(0,\,N-1)$ has depth $d=\lceil\log_{2}N\rceil$ and number of root-children equal to $n=\lfloor\log_{2}N\rfloor$. An example of a Fenwick tree for $N=7,\,d=3$ is shown in Fig. 1. The partial sums $x_{j}$ of the encoded representation are given by a (mod 2) sum of $j$-th fermionic occupancy $n_{j}$ with the descendants of $j$ in the Fenwick tree. For example the zeroth bit encoded by a Fenwick tree in Fig. 1 stores only the occupancy of the zeroth fermionic site as it has no descendants, while the first bit stores $x_{1}=n_{1}+x_{0}=n_{1}+n_{0}$. Likewise, the sixth bit has $\{3,5\}$ as its children and therefore stores $x_{6}=n_{6}+(x_{3}+x_{5})=n_{0}+n_{1}+n_{2}+n_{3}+n_{4}+n_{5}+n_{6}$. The remaining bits are given by:

[TABLE]

As a specific example, $n_{0}\,n_{1}\ldots n_{6}=0111010$ is encoded as $x_{0}\,x_{1}\ldots x_{6}=0111010$.

II.2.2 Bravyi-Kitaev Transformation

The BK transform uses Fenwick trees to improve qubit operator locality of the fermionic parity counting string to $O(\log N)$, while increasing fermionic occupancy update cost to $O(\log N)$. The raising/lowering operators are hence mapped with $O(\log N)$ operator locality overhead, which is substantially better than $O(N)$ for JW.

Starting with the simplest example, consider $a^{\dagger}_{2}$ applied to the second fermionic site in a qubit register $\ket{x_{0}\,x_{1}\ldots x_{6}}$ encoding an occupancy state $\ket{n_{0}\,n_{1}\ldots n_{6}}$ of 7 fermionic sites as in Fig. 1. This operator acts as:

[TABLE]

on the encoded states, as one needs to count the excitation parity of [math] and $1$ by applying $Z_{1}$, change the occupancy of the second node by applying $\ket{1}\bra{0}_{2}$ and ensure consistency of the encoding by updating sites 3 and 6 (ancestors of 2) with $X_{3}$ and $X_{6}$.

Mapping $a_{j}^{\dagger}$ for a general $j$ is more complicated, as one has to condition application of $\ket{0}\bra{1}_{j}$ or $\ket{1}\bra{0}_{j}$ on content of children of $j$ in the Fenwick tree. This is the case for $j=3$ in Fig. 1 for example. If the third fermionic site is initially unoccupied ($n_{3}=0$), the raising operator changes $n_{3}$ from 0 to 1. In the encoded representation, the third qubit stores $x_{3}=n_{0}+n_{1}+n_{2}+n_{3}=(x_{1}+x_{2})+n_{3}$. So if $(x_{1}+x_{2})=1$, a qubit lowering operator $\ket{0}\bra{1}_{3}$ should be applied in the encoded representation instead of $\ket{1}\bra{0}_{3}$. It follows that one has to condition this operation on the children’s parity. The operator hence maps to:

[TABLE]

The above description considerably simplifies by working in the Majorana basis:

[TABLE]

Consider now the set of children with indices less than $j$ of all ancestors of $j$. We label this set as $C(j)$. For example, the set of children of all ancestors of qubit $9$ in Fig. 3 is given by $\{7,\,9,\,10,\,11,\,13,\,14\}$ and out of this, only $7$ is less than $9$ and hence $C(9)=\{7\}$. For consistency with refs. Seeley ; Tranter , we denote the set of children of the $j$-th site by $F(j)$ and work with a set $P(j)=C(j)\cup F(j)$. If $U(j)$ labels the set of all ancestors of $j$, then:

[TABLE]

where $Z_{P(j)}$ implies Pauli $Z$ operators applied to qubits in a set $P(j)$.

Note that $P(j)\cap U(j)=\emptyset$, since all nodes in $U(j)$ have indices greater than $j$ while $P(j)$ have all indices less than $j$. Also note that the $d_{j}$ operator acts trivially on the $F(j)$ qubits. Locality of the $c_{j}$ Majorana on qubits is hence never better than $d_{j}$, since no operators are applied to children of the $j$-th node (Eq. 3). The worst-case locality for the $c_{j}$ Majorana operator is therefore given by $|U(j)\cup P(j)|+1=|U(j)|+|P(j)|+1$. In fact, for a Fenwick tree of $N=2^{d}$ sites, the locality of $c_{j}$ becomes exactly $\log_{2}N+1$ as we now show:

Proof.

Let $d=0$. Then $N=1$ and the $c_{j}$ locality is $\log_{2}1+1=1$ as there is only single node in the tree. Now suppose the locality is $\log_{2}N+1$ for a tree with $N=2^{d}$ nodes. The $2N=2^{d+1}$ tree is constructed by connecting roots of two $N$ trees - compare for example the descendants of $7$ to the tree of the remaining nodes in Fig. 10. Any node in the right subtree has to update the new root, which worsens the locality by $1$ - hence the terms are now $\log_{2}2N+1$ local. Operators on the rest of the tree will have to lookup an additional node to obtain the parity, which also implies $\log_{2}2N+1$ locality. So every $c_{j}$ on a tree with $N=2^{d}$ nodes is $\log_{2}N+1$ local. ∎

We now provide a unifying framework for JW and BK encodings and introduce an optimized variant of the method suitable for rectangular lattice geometries.

II.3 Fenwick Trees as a Class of Fermionic Encodings

The recursive Algorithm 1 from the previous section gives rise to a class of fermionic encodings with JW and BK schemes as limiting cases. Instead of using a single Fenwick tree to encode all fermionic modes, we partition the fermionic sites into Fenwick trees of varying depth and include the set of all roots of segmented trees less than $j$ to $P(j)$. The definition of $P(j)$ in this context becomes:

[TABLE]

In particular, we could choose to put each node in its own Fenwick tree of depth 0, which would correspond directly to the JW transformation (Fig. 4).

Every iteration of the algorithm defines new representation of fermionic algebra with qubit operators. With increasing recursion depth, the occupancy update cost worsens, while the parity counting costs improve. Asymptotically, the representation transitions from operator locality overheads of $O(1)\rightarrow O(\log N)$ for occupancy update and $O(N)\rightarrow O(\log N)$ for parity computation.

II.4 Segmented Bravyi-Kitaev (SBK) transform

It is possible to use the segmented transform for operator locality improvement of specific fermionic Hamiltonians on rectangular lattices. As already discussed, the BK transform optimizes update and lookup costs simultaneously. Composition of these operations describes the action of fermionic raising/lowering operators, but does not strictly correspond to operators occurring in the Hubbard model - or in fact any other physical fermionic Hamiltonian. For closed systems, physical fermionic Hamiltonians only contain pairs of raising and lowering operators Kitaev , which places further constraints on the set of operators we need to map to qubits. If we additionally restrict our attention to operators which are local (as is the Hubbard model), there is a lot of redundancy one can exploit for further qubit operator locality optimization.

We focus on the case of a $w\times h,\;w\leq h$ rectangular lattice and build a Fenwick tree for every row - this is a specific example of segmented Fenwick tree as defined in the previous section. At first sight, this appears to worsen locality of the qubit operators as the locality of single raising/lowering operators now scales asymptotically as $O(h\log w)$. If we however restrict the set of possible operations to on-site and nearest-neighbor terms, the single-qubit operators on segmented tree roots cancel and the locality becomes $O(\log w)$ - a substantial improvement, as the operator locality is now independent of the lattice height.

It turns out to be slightly more optimal to store two disconnected trees per lattice row, as shown in Fig. 5. This is because the full parity of the row is not necessary for the vertical nor horizontal hopping terms.

II.5 Loop-Stabilized Fermion Simulation (LSFS)

We now shift our attention to the second method introduced in Ref. Kitaev which finds an alternative representation of the local fermionic Hamiltonians on qubits in a line graph (i.e. on graph edges - see Fig. 6). The method, referred to as a “Superfast Simulation of Fermions” in Kitaev improves qubit operator locality to a constant for any bounded-degree graph. The newly proposed name “Loop-Stabilized Fermion Simulation” is motivated by the fermionic algebra being represented in a subspace defined by a stabilizer condition on the set of all possible loops in the line graph.

II.5.1 Loop-Stabilized Fermionic Simulation

As discussed in the construction of the SBK method, physical fermionic Hamiltonians are sums or products of fermionic raising/lowering operator pairs. Equivalently, all fermionic Hamiltonians can be obtained by combining the following operators Kitaev :

[TABLE]

Here, the subscript on $B$ corresponds to a vertex and the subscript on $A$ labels a graph edge. Specifically, the fermionic site to site hopping is expressed by:

[TABLE]

More formally, the $A_{(jk)},\,B_{j}$ operators generate the algebra of physical fermionic Hamiltonians. The algebra defining rules are Kitaev :

[TABLE]

In other words, $A_{(jk)}$ anticommutes with any other generator as long as they share precisely one vertex. Furthermore note that $B_{k}^{\dagger}=B_{k}$, $A_{(ij)}^{\dagger}=A_{(ij)}$, $B_{k}^{2}=A_{(ij)}^{2}=\mathbf{1}$ and $A_{(ij)}=-A_{(ji)}$. Additionally:

[TABLE]

for any closed path $j_{0}\,j_{1}\,j_{2}\ldots\,j_{p}$.

Bravyi and Kitaev found in Kitaev qubit operators $\tilde{B}_{l},\,\tilde{A}_{(jk)}$ obeying the above algebra, subject to further constraints on excitation parity sector that we now impose. For an even number $N_{f}$ of fermions, one has that:

[TABLE]

as an additional rule to the above set. This motivates the following choice for the qubit operator $\tilde{B}_{k}$:

[TABLE]

where $n(k)$ is the set of nearest neighbors of $k$. The $\prod_{k\in V}\tilde{B}_{k}~{}=~{}\mathbf{1}$ condition is then trivially satisfied for the operator since each edge shares precisely two vertices (this is colloquially known as the “handshaking lemma”).

We now derive the form of $\tilde{A}_{(jk)}$. Firstly assume that $\tilde{A}_{(jk)}$ is a tensor product of Pauli operators and/or identity on edges adjacent to vertices $j,\,k$. In order to satisfy $A_{(jk)}B_{l}=(-1)^{\delta_{jl}+\delta_{kl}}B_{l}A_{(jk)}$, we first address the case for which $k\neq l,\,j\neq l$, so that the fermionic operators commute. If $\tilde{A}$ and $\tilde{B}$ are to obey the same relation, all single qubit operators of $A_{(jk)}$ on edges adjacent to $j,\,k$, except for the edge $(j,k)$, must be Hermitian and identity-squaring operators in a subspace spanned by $\{\mathbf{1},Z\}$. There are only two options - either $Z$ or $\mathbf{1}$. Focusing on the case $l\in\{j,k\}$, the $A_{(jk)}B_{l}$ operators anti-commute, implying that the operator on the $(j,k)$ edge qubit is a Hermitian and identity-squaring operator in the subspace spanned by $\{X,Y\}$ - again $X$ or $Y$ are the only possibilities.

It remains to satisfy the condition imposed by the generalized commutator of two $A_{(jk)},A_{(lm)}$ operators. These anticommute iff they share a vertex. An example of qubit operators $\tilde{A}_{(jk)}$ satisfying this and the previous constraints on a square lattice is given in Fig. 7.

An example of operators satisfying these constraints on a general lattice is given by:

[TABLE]

and is a specific case of the prescription found in Ref. Kitaev . Lastly, we impose anticommutation of $A_{(jk)}$ using antisymmetric tensor $\epsilon_{jk}$ which is +1 when $j>k$ and $-1$ when $j<k$. Thus,

[TABLE]

Even though this implies the generator itself is multiply defined for each edge, we will see that all interesting physical operators will be independent of $\epsilon_{jk}$.

The loop condition of Eq. 4 will be imposed by a set of stabilizer operators which can be concisely presented with a specific lattice geometry in mind. We therefore postpone its discussion until after the following example.

II.5.2 Example: 2D Hubbard Model

We illustrate this method by mapping the Hubbard Hamiltonian on rectangular lattice to qubits. Let $Z^{\uparrow}_{k}$ denote a Pauli $Z$ operator applied to the qubit on the vertical edge adjacent to the vertex $k$ - if there is no such edge, substitute the term with an identity operator (for example, $Z_{1}^{\uparrow}=\mathbf{1}$ in the diagram in Fig. 6). Operators on other adjacent edges are defined analogously by using $\{\rightarrow,\leftarrow,\uparrow,\downarrow\}$ superscripts.

The $\tilde{B}_{k}$ operator is represented with a “cross” of Pauli $Z$ operators:

[TABLE]

The form of $\tilde{A}_{(jk)}$ differs for horizontal and vertical edges, which we denote by $E_{H}$ and $E_{V}$ respectively. It also depends on a specific lattice indexing - we choose the simplest one shown in Fig. 6 and leave it an open question whether a better ordering exists. The operator $\tilde{A}_{(jk)}$ is then given by:

[TABLE]

It remains to satisfy the loop condition of Eq. 4, since in order to represent the fermionic algebra, the qubit operators must obey the same relation. Because the fermionic operators obey:

[TABLE]

for any plaquette $(\alpha\beta\gamma\delta)$ and an arbitrary closed loop of Eq. 4 can be obtained by taking product of such plaquette operators, encoded physical fermionic states correspond to qubit states which are $+1$ eigenstates of $\tilde{A}_{(\alpha\beta)}\tilde{A}_{(\beta\gamma)}\tilde{A}_{(\gamma\delta)}\tilde{A}_{(\delta\alpha)}$. Formally, we restrict the qubit states by a set of stabilizer operators:

[TABLE]

where $(\alpha\beta\gamma\delta)$ labels vertices of the plaquette. An example of such stabilizer is worked out in Figure 6 for a trivial lattice ordering. Since $C_{(\alpha\beta\gamma\delta)}$ is a Krönecker product of Pauli matrices, its eigenvalues are $\pm 1$, implying that violation of a single stabilizer condition costs $2$ units. Note that if the simulation starts in the subspace encoding physical fermionic states, it stays in it, as the generating operators $B_{k},\,A_{(jk)}$ commute with the stabilizers Kitaev .

Inspired by the toric code construction Toric , we note that one can restrict the system to the physical codespace by including a penalty term $H_{penalty}=-\frac{\Delta}{2}\sum_{k}C_{k}$, where $\Delta\gg t,U,\epsilon$ corresponds to the “energy gap” of the system and $k$ runs over all plaquettes.

The above representation is now applied to the fermionic nearest-neighbor hopping operator. The nearest-neighbor couplings for horizontal edges maps to $5$-qubit-local operators:

[TABLE]

Analogously, the vertical nearest-neighbor couplings are encoded by 7-local operators:

[TABLE]

It remains to account for the Hubbard repulsion term. The simplest way to implement it is by simulating the above on two lattices labeled with $\uparrow,\downarrow$, coupled by the density-density interaction term. We obtain that:

[TABLE]

which leads to the following expression for the full Hamiltonian:

[TABLE]

The $H_{\uparrow}$ and $H_{\downarrow}$ Hamiltonian terms are decoupled and have been already mapped to spins. The density-density interaction term maps to:

[TABLE]

where the primed indices correspond to fermions in spin $\downarrow$ lattice and the unprimed ones to the sites in the spin $\uparrow$ lattice. It hence follows that in order to implement this inter-lattice coupling on spins, the density-density Hamiltonian term acts on 8 qubits simultaneously in the worst case. This presents an upper bound on the spin Hamiltonian locality in this setting.

II.6 The Auxiliary Fermion Scheme

The auxiliary fermion (AF) scheme us ; Ball ; Verstraete was introduced in Ball and also (independently) by Verstraete . The scheme uses auxiliary fermionic modes to allow fermionic models on general lattice geometries to be simulated locally. We recall the relevant details from us to allow comparison against the other fermionic encoding schemes.

A given set of site-to-site hopping terms can be characterized by a graph $G=(E,\,V)$, where vertices correspond to sites and the edges to pairs of sites participating in the hopping. In one dimension, this graph is a linear path, $G_{1}$, where the degree of the endpoints is one and all other sites have degree two. In higher dimensional settings such as the $2$D Hubbard model studied here, the degree of the sites may be greater than two. The non-local degree $d_{nl}$ of a site is the number of edges that are not included in the linear path $G_{1}$. Therefore, it is important to choose a path which overlaps maximally with the desired interaction graph $G$. In the case of $2$D Hubbard, this is accomplished by a snake-like pattern (Fig. 8).

Next, we must account for spin in the Hubbard model. One requires $2wh$ spinless fermionic sites if the model is defined on $w\times h$ rectangular lattice. If we take the first half of these sites to be spin down and the second half to be spin up, there is no need to track phase factors between the spin up and spin down subsets. This is because the Hubbard model preserves spin; thus, there are no hopping terms between the first half of the sites and the second half. The density-density terms are each a product of two one-point coupling terms which does not require the tracking of antisymmetric phase factors.

The nonlocal degree of each node (when $G_{1}$ is a subgraph of $G$) is $d_{nl}=d(G)-d(G_{1})$ and, as we’ve previously shown, each auxiliary fermionic site can facilitate up to two non-local couplings us . Hence the number of auxiliary sites required per fermionic site is given by $\lceil d_{nl}/2\rceil$. Since each non-local degree in the $2$D Hubbard model is less than or equal to two, only a single auxiliary mode is needed for each site with non-local degree greater than zero. The $(1,w)$ and $(h,1)$ corners (for each spin) require no auxiliary sites since the edge set of $G_{1}$ is sufficient for coupling them to their nearest neighbors.

As an example, we illustrate these ideas using a square lattice of width $w=3$ and height $h=3$. In the general case, there are $(w-2)(h-2)$ sites with degree four in the interior of the graph $G$. Each of these sites has non-local degree two and hence each needs just a single auxiliary fermionic site. Then there are $2(h-2)+2(w-2)$ sites along the boundary with degree three. Each of these sites have non-local degree of 1 and will also require one auxiliary fermionic site each. Finally, the four sites occupying the corners of the lattice to be simulated have degree two. Presuming that the snake-like pattern is used, two of the corners will have non-local degree zero while the other two each have non-local degree 1. Summing all the cost together we find that the total number of qubits is $4wh-4$.

The locality of simulation operators was detailed in our previous publication us . For a hopping term from site $i$ to a non-consecutive site $j$, the qubit operator acts on site $i$ and on site $j$ as well as auxiliary modes associated with site $j$ and with site $i$. Hence, the qubit operator is either $4$-local or $3$-local in the Hubbard model. The density-density operators are $2$-local us .

III Comparison

In the following we assess the various methods in context of the $2$D Hubbard Hamiltonian. We analyze locality of the nearest-neighbor fermionic hopping and density operators. The first subsection goes through locality analysis for $2$D Hubbard model. We then present results in higher dimensions.

III.1 2D Hubbard Model

We consider a Hubbard model defined on a $w\times h,\;w\,<\,h$ rectangular lattice. For JW, we order the fermions as in Fig. 8. The longest string of $Z$ operators introduced by the mapping has then length $(w+1)$.

For LSFS, we have already shown that the density-density terms in the Hubbard Hamiltonian are all $8$-local for lattices with $w,\,h\geq 3$. It is therefore only sensible to use LSFS for locality reduction if $w\geq 8$. Likewise AF only improves locality compared to JW if $w\geq 4$, since the most nonlocal term of the qubit Hamiltonian is the 4-local vertical hopping term.

The AF transform is hence superior to LSFS for the Hubbard model defined on a rectangular lattice. Both methods use more qubits than fermionic sites of the original model. For a $w\times h$ rectangular lattice, LSFS requires $4(w-1)(h-1)$ qubits, while the auxiliary fermion method needs $4wh-4$, compared to a minimum of $2wh$.

We now compare these methods to the BK transform and its optimized SBK variant introduced in subsection II.4. The density operator can be written in Majorana basis and converted by BK transform as:

[TABLE]

Its locality therefore only depends on $|F(j)|$, i.e. the number of children of a node $j$ in the Fenwick tree. The root has the largest number of children, which implies that the worst-case locality is $\left(\lfloor\log_{2}N\rfloor+1\right)$. The density-density term is hence $\left(2\lfloor\log_{2}N\rfloor+2\right)$-local.

Using the expression for Majorana operators in the Bravyi-Kitaev mapping (Eq. 3), we have for the hopping operator that:

[TABLE]

If all-to-all couplings were allowed, the worst-case locality would be $\lfloor\log_{2}N\rfloor+\lceil\log_{2}N\rceil$ - the sum of tree depth $d$ and the number of root children $n$. This corresponds to hopping from the first to the last site of the lattice (see Fig. 10).

Focusing our attention to the nearest-neighbor hopping in the SBK transformation, we have to analyze the vertical and horizontal hoppings independently. For hopping along horizontal edges, a loose bound on worst case locality is given by $2\lceil\log_{2}w\rceil-1=2d-1$ and corresponds to hopping to the root child with smallest index $j$ from the $(j+1)$-th node.

The worst case locality for hopping along the verticals is given by $2\lceil\log_{2}w\rceil+2=2d+2$ (where $d$ is the tree depth) and corresponds to hopping between two deepest leaf nodes. This can be reduced to $2\lceil\log_{2}w\rceil+1=2d+1$ by splitting the segmented trees in half, which however increases the worst case locality for horizontal hoppings to $2\lceil\log_{2}w\rceil$. The scheme hence provides locality advantage compared to JW for lattices of width $w>2$. The worst case localities for $2$D are sumarized in Tab. 1 and in Fig. 11.

III.2 Hypercubic Lattices

We now extend the previous analysis to higher dimensions and consider mapping of the hopping term in a general hypercubic lattice of $D>0$ dimensions with side $w$. With the simplest ordering, the hopping term in Jordan-Wigner transformation becomes $w^{d-1}+1$ local, where $w$ is the number of sites along a side of the hypercubic lattice.

In the case of AF, the number of auxiliaries per site scales as $D-1$, since in $D$-dimensional hypercubic lattice, each vertex has $2D$ nearest neighbors. This translates to $d_{nl}=2D-2$ non-local degree and $D-1$ auxiliary fermion sites per each site of the original lattice. The worst case hopping term locality hence goes as $2D-2$, whereas the number of qubits scales as $D\times w^{D}$.

Since each vertex has $2D$ neighbors, locality of the LSFS hopping term goes as $4D-1$ in a bulk of the hypercube (one neighbor is shared), while locality of the density-density term becomes $4D$. The number of edges in a hypercubic lattice of side $w$ in dimension $D$ is given by $E(D,w)=D(w-1)w^{D-1}$, which is also the number of qubits required for the mapping.

Proof.

One has that $E(2,w)=2(w-1)w$, as there are $2$ edges per vertex with the exception of the boundary, where we over-count by $2w$. In $3$ dimensions, we construct a $w\times w\times w$ cube by connecting vertices of $w\times w$ squares by $(w-1)w^{2}$ edges. It follows that $E(3,w)=wE(2,w)+w^{2}(w-1)=4w^{2}(w-1)$. To count the number of edges in a $w^{\times D}$ hypercube, one can take $w$ $(D-1)$-dimensional hypercubes and connect them by $(w-1)w^{D-1}$ edges. This implies that $E(D,w)=wE(D-1,w)+w^{D-1}(w-1)$. It follows that $E(D,w)=D(w-1)w^{D-1}$. ∎

The BK method requires $N=w^{D}$ qubits in general, while its worst-case locality scales as $\lceil\log_{2}(w^{D})\rceil+\lfloor\log_{2}(w^{D})\rfloor=d+n$. For nearest neighbour hopping, this can be further optimized to $2\lceil\log_{2}(w^{D-1})\rceil+1=2d+1$ by using the segmentation trick of the SBK method at the highest level of the lattice. Hopping term locality as a function of hypercubic lattice dimension is summarized in Tab. 2.

IV Conclusion

Which scheme is superior depends on the type of quantum simulation. For digital quantum simulation on a universal quantum computer with logical qubits, $N$-qubit operators can be simulated with just linear overhead in $N$. On such devices the main limitation will be the number of available logical qubits and then our proposed modification of the Bravyi-Kitaev transform leads to the best locality improvement of $O(\log_{2}w)$-local spin operators, where $w$ is the lattice width.

For analog simulations, on the other hand, operator locality will be the decisive factor. On such quantum simulators natively only few-qubit couplings are available, typically only two-qubit terms. In that setting multi-qubit terms have to be generated using perturbative gadgets KKR ; Oliveira . These, however, require large energy penalties to be sufficiently deep inside the perturbative regime where the effective higher-order interactions appear. One hence wants to optimize the locality of the terms and aim for transformations with the most local terms. The presented analysis shows that it is possible to map the $2$D Hubbard to a $4$-local qubit Hamiltonian by the Auxiliary Fermion method, at the expense of using number of ancillary qubits. This is the optimal fermion representation for 2D lattices of width $w\geq 4$.

V Acknowledgements

We would like to thank to Peter D. Johnson, Peter Winkler, Amit Chakrabarti, Thomas H. Cormen and Alexey Soluyanov for useful discussions. V. H. would like to thank Dartmouth College for support and hospitality while finishing the work and the Clarendon and Keble de Breyne scholarships for support. This project was supported by the Swiss National Science Foundation through the National Competence Center in Research NCCR QSIT and by the European Research Council through ERC Advanced Grant SIMCOFE.