††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Low-depth simulations of fermionic systems on square-grid quantum hardware

Manuel G. Algaba  IQM, Nymphenburgerstr. 86, 80636 Munich, Germany  P.V.
Sriluckshmy  IQM, Nymphenburgerstr. 86, 80636 Munich, Germany  Martin Leib
IQM, Nymphenburgerstr. 86, 80636 Munich, Germany  Fedor Šimkovic IV  IQM,
Nymphenburgerstr. 86, 80636 Munich, Germany

###### Abstract

We present a general strategy for mapping fermionic systems to quantum
hardware with square qubit connectivity which yields low-depth quantum
circuits, counted in the number of native two-qubit fSIM gates. We achieve
this by leveraging novel operator decomposition and circuit compression
techniques paired with specifically chosen _low-depth fermion-to-qubit
mappings_ and allow for a high degree of gate cancellations and parallelism.
Our mappings retain the flexibility to simultaneously optimize for qubit
counts or qubit operator weights and can be used to investigate arbitrary
fermionic lattice geometries. We showcase our approach by investigating the
tight-binding model, the Fermi-Hubbard model as well as the multi-orbital
Hubbard-Kanamori model. We report unprecedentedly low circuit depths per
single Trotter layer with up to a $70\%$ improvement upon previous state-of-
the-art. Our compression technique also results in significant reduction of
two-qubit gates. We find the lowest gate-counts when applying the XYZ-
formalism to the DK mapping. Additionally, we show that our decomposition and
compression formalism produces favourable circuits even when no native
parameterized two-qubit gates are available.

##  1 Introduction

One of the most promising applications of future quantum computers is the
simulation of fermionic quantum systems. These systems are of major relevance
to various scientific fields, such as quantum chemistry, condensed matter
physics, lattice gauge theories and nuclear physics. Simultaneously, they are
related to many unsolved scientific problems with potential industrial
applications, such as high-temperature superconductivity [1], battery design
[2], chemical reaction rates [3], and nitrogen fixation in fertilizers [4].

Classically, despite extensive scientific efforts in developing and
benchmarking computational tools [5, 6, 7], most such systems remain too
difficult to study with current state-of-the-art techniques. This is mainly
due to their exponentially growing computational space [8] as well as due to
the anti-commuting nature of fermions, which leads to the infamous sign
problem for many computational methods based on quantum Monte Carlo [9, 10].

Whilst it is generally believed that quantum computers will eventually be able
to circumvent these difficulties and provide at least polynomial speedups [11]
for problems of relevance, finding suitable quantum, or hybrid quantum-
classical, algorithms is an active field of research. The need for developing
highly optimized quantum algorithms is further amplified by the shortcomings
of currently available noisy quantum hardware, which is severely limited in
their qubits counts, coherence times and gate fidelities [12]. Small scale
proof-of-principle implementations of algorithms on such NISQ hardware have
nevertheless been performed, encouraging further research in the field [13,
14, 15].

A necessary first step for any quantum computing approach to fermionic systems
is to find an adequate transformation between the fermionic Fock space and the
computational space of a multi-qubit quantum device. This transformation is
not unique and _fermion-to-qubit mappings_ can have vastly different
properties from one another. The right choice of mappings is thus fundamental
to successfully simulating fermionic systems.

The Jordan-Wigner transformation (JW) [16] is not only historically the first
fermion-to-qubit mapping, but also the most popular choice today, especially
due to its simplicity and favourable performance for small systems. It maps
fermionic modes one-to-one to qubits arranged on a line at the expense of
generating qubit operators whose weights scale linearly with the total system
size.

Numerous mappings have been since proposed that improve the JW transformation
in terms of specific optimisation criteria, especially in cases where one can
exploit some known underlying structure in the fermionic Hamiltonian to be
simulated. In particular, some mappings aim at reducing the total number of
qubits needed to encode a fermionic system by using symmetries [17, 18, 19,
20, 21, 22], by encoding fermionic modes through tree-graph structures [23,
24], or by considering specific hardware topology and connectivity structures
[20, 25]. Other mappings are designed with partial quantum error correction in
mind [26], and allow to either identify [27, 20] or correct [28, 29, 30]
errors that occur during circuit execution.

Recently, the focus in the study of fermion-to-qubit mappings has increasingly
shifted towards efficient implementation of mappings on near-term (NISQ)
quantum devices [31, 32, 33]. Indeed, current NISQ systems are mostly limited
by gate count and coherence times. This is reflected in the fact that even
though processors with hundreds of qubits are currently available, the largest
reported fermionic simulations to date typically use less than thirty [34].

It is often implicitly assumed that the mappings which yield the lowest
circuit depths are the ones with the lowest operator weights, corresponding to
the number of qubits they act upon. These weights can be reduced by
introducing additional ancilla qubits used to resolve fermionic commutation
relations [35, 36, 37, 38, 39, 40, 41]. In particular, many (geometrically)
_local_ mappings, with circuit depths which are independent of fermionic
lattice size, have been proposed. Alternatively, one can dynamically alter the
association of fermionic modes with hardware qubits using fermionic swap
(fSWAP) networks [35, 42, 43, 44, 32].

For fermionic systems with some degree of sparsity it turns out that a
combination of local mappings and fSWAP networks yields the lowest-depth
circuits to date [31, 32]. This approach has reduced the necessary circuit
depths per Trotter layer by multiple orders of magnitude to the range of
roughly $10^{2}-10^{4}$, depending on the fermionic model to be simulated. The
lingering question is whether further improvements to those depths are
possible, especially since reducing the circuit depths by even a small
constant factor can potentially accelerate the advent of useful fermionic
simulations on quantum computers by years. Another issue with many of the
existing fermion-to-qubit mappings is that they assume qubit layouts that are
tailored towards one singular model or demand connectivity graphs that are not
compatible with the restricted qubit topologies provided by some of the
leading quantum computing hardware platforms. It is therefore important to
investigate mappings that bring future fermionic simulation away from the
hypothetical and closer to the practical.

In this work we focus on designing a universal strategy of simulating multi-
orbital fermionic systems on an arbitrary two-dimensional lattice by means of
a quantum processor with a square qubit connectivity graph, which is one of
the standard topologies found in the superconducting quantum computer industry
[45, 46]. Our particular focus is on generating minimal-depth quantum
circuits, counted in the number of parallelizable native two-qubit gate layers
per single Trotter step. To this end, we present operator decomposition and
compression techniques which we apply to a number of fermion-to-qubit mappings
which we introduce in this paper. We show that our formalism generally
improves the circuit depths for the investigated models and that mappings with
the lowest operator weights are not necessarily optimal in terms of circuit
depth because of their respective lack in parallelism.

The paper is structured as follows: In Section 2 we introduce our
decomposition and compression techniques (which we collectively name XYZ). In
Section 3 we introduce the fermionic models of interest and in Section 4 we
discuss the general formalism we use to map them to square qubit layouts. In
Section 5 we present our analysis of circuit depths obtained by using the XYZ
formalism in combination with various local fermion-to-qubit mappings for
various fermionic lattice geometries. We inspect the tight-binding model (TB)
in Sec. 5.1, the Fermi-Hubbard model (FH) in Sec. 5.2, and the Hubbard-
Kanamori model (HK) in Sec. 5.3. In Sec. 6 we investigate how the advantage of
using the XYZ formalism depends on the type of available native gates. We
finish with a discussion of our main findings as well as future directions in
Section 7.

##  2 XYZ operator decomposition and circuit compression techniques

Most quantum algorithms involve multi-qubit operators, yet current quantum
processors only allow for one and two-qubit gates to be executed at any given
point. As a consequence, it is necessary to decompose such multi-qubit
operators into circuits involving executable native gates. Depending on which
gates are considered native, there are multiple ways to decompose operators
which yield different levels of efficiency in terms of resulting gate counts
and circuit depths. Let us consider as example the operator $e^{i\alpha
Z_{1}Z_{2}Z_{3}Z_{4}}$ where $\alpha\in\mathbb{R}$ and $Z_{i}$ is a Pauli
operator acting on qubit $q_{i}$ 111In what follows we label circuit qubits
with ascending integers from top to bottom.:

$\raisebox{-0.45pt}{\includegraphics[scale={0.08}]{cnotlad1.png}}=\raisebox{-0.42pt}{\includegraphics[scale={0.08}]{cnotlad2.png}}=\raisebox{-0.42pt}{\includegraphics[scale={0.08}]{cnotlad3.png}}$
(1)

We can decompose the operator (l.h.s. of Eq. 1) using the standard approach in
the field [47], which generates a _V-shaped_ circuit involving CNOT gates and
single qubit rotations (left circuit in Eq. 1). However, it is more efficient
in terms of circuit depth to decompose the operator into an _X-shaped_ circuit
instead (right circuit in Eq. 1) [48]. For multiple operators acting on
overlapping sets of qubits one can further compress the resulting circuits by
combining their outermost CNOT gates using known identities (we list these in
Appendix A).

Here, we use an alternative, recently introduced technique, named XYZ-
decomposition [49] which, as we will proceed to show, yields considerably
shallower circuits compared to the decomposition into CNOTs, especially when
parameterized two-qubit gates are available. The XYZ technique is based on the
decomposition of a multi-qubit gate of the form $e^{i\alpha\mathcal{O}}$,
where $\mathcal{O}$ is a tensor product of Pauli operators on a given number
of qubits, into three operators of lower or equal weight:

$e^{i\alpha\mathcal{O}}=e^{i\frac{\pi}{4}\mathcal{O}_{1}}e^{i\alpha\mathcal{O}_{2}}e^{-i\frac{\pi}{4}\mathcal{O}_{1}}$
(2)

where $\mathcal{O},\mathcal{O}_{1},\mathcal{O}_{2}$ are unitary and fulfil the
relations $\mathcal{O}^{2}=\mathbb{1}$ and
$\mathcal{O}=\frac{i}{2}[\mathcal{O}_{1},\mathcal{O}_{2}]$. All three
operators must act on a connected set of qubits and $\mathcal{O}_{1}$,
$\mathcal{O}_{2}$ are acting on a subset of the qubits that $\mathcal{O}$ is
acting on. For simplicity we assume a linear connectivity, meaning qubits can
interact with up to two of their neighbors. This approach can be recursively
applied to obtain a complete decomposition of the original multi-qubit
operator into two-qubit gates. As an example, let us consider the following
operator acting on four qubits:

$\displaystyle e^{i\alpha
X_{1}Z_{2}Z_{3}Y_{4}}=e^{i\frac{\pi}{4}X_{1}Y_{2}}e^{i\alpha
X_{2}Z_{3}Y_{4}}e^{-i\frac{\pi}{4}X_{1}Y_{2}}$ (3)
$\displaystyle=e^{i\frac{\pi}{4}X_{1}Y_{2}}e^{i\frac{\pi}{4}X_{3}Y_{4}}e^{-i\alpha
X_{2}Y_{3}}e^{-i\frac{\pi}{4}X_{1}Y_{2}}e^{-i\frac{\pi}{4}X_{3}Y_{4}}$

Here, we first decompose the operator at the _central_ qubit $q_{2}$ to
generate two- and three-qubit operators. Subsequently, we further decomposed
the three-qubit operators by breaking it up at $q_{3}$. Depending on the
prefactor in the exponential we can introduce a graphical notation for the
three different types of two-qubit gates occurring in such decompositions:

$\raisebox{-0.45pt}{\includegraphics[scale={0.08}]{xyzdecomp1.png}}\equiv
e^{i\frac{\pi}{4}\sigma_{1}\rho_{2}},\
\raisebox{-0.44pt}{\includegraphics[scale={0.08}]{xyzdecomp2.png}}\equiv
e^{i\alpha\sigma_{1}\rho_{2}},\
\raisebox{-0.44pt}{\includegraphics[scale={0.08}]{xyzdecomp3.png}}\equiv
e^{-i\frac{\pi}{4}\sigma_{1}\rho_{2}}$ (4)

These are the only three types of two-qubit gates generated by the XYZ-
decomposition (single-qubit gate equivalents of the above can also appear).
The two gates with prefactors $\pm\pi/4$ are simply special cases of the gate
with prefactor $\alpha$. Here, $\sigma$ and $\rho$ are Pauli operators acting
on qubits $q_{1}$ and $q_{2}$, respectively. Note that the square two-qubit
operator does not indicate the value of $\alpha$, which has to be tracked
separately. Only a single such two-qubit gate is generated per multi-qubit
operator. If parameterized two-qubit gates are not available it can be further
decomposed into a parameterized single-qubit gate and two two-qubit gates
using Eq. 2.

Generally, one has a choice of which qubit to designate as central for the
next decomposition as well as the choice of how to distribute the two new
Paulis of the new operators acting on this qubit. Let’s consider the operator
$e^{i\alpha Y_{1}Z_{2}Z_{3}X_{4}}$ for which we show below four of the
possible decompositions:

$\begin{matrix}\raisebox{-0.46pt}{\includegraphics[scale={0.08}]{legs1.png}}\
=\ \raisebox{-0.46pt}{\includegraphics[scale={0.08}]{legs3.png}}\ =\
\raisebox{-0.46pt}{\includegraphics[scale={0.08}]{legs5.png}}\ =\\\ \\\
\quad\quad=\ \raisebox{-0.46pt}{\includegraphics[scale={0.08}]{legs2.png}}\ =\
\raisebox{-0.46pt}{\includegraphics[scale={0.08}]{legs4.png}}\end{matrix}$ (5)

Starting with four-qubit operators one has the choice of decomposing operators
into a V-shaped or an X-shaped circuit. As shown in Ref. [49], the latter
always produces lower-depth circuits and we found that this statement
generally holds true even when multiple operators are considered
simultaneously. It is also possible to construct _asymmetric_ X-shaped
circuits for operators acting on five or more qubits, but this rarely leads to
improvements in overall circuit depth. One can also always turn around the
shape of the decomposition, vertically, in some sense reminiscent of using the
Yang-Baxter equation to decrease the depths of standard circuits [50].

Let us now turn to compressing circuits resulting from applying multiple
multi-qubit operators acting on overlapping sets of qubits. Our main focus is
on treating operations between consecutive two-qubit gates, as all single-
qubit gates appearing in between them can be transported through to the side
as long as they have a $\pi/4$ prefactor. Consider an arbitrary two qubit gate
of the form $e^{i\alpha\sigma^{a}_{1}\rho_{2}}$ and a single qubit gate of the
form $e^{\pm i\frac{\pi}{4}\sigma^{b}}$. Then we have:

$\raisebox{-0.4pt}{\includegraphics[scale={0.08}]{transporting_a.png}}=\raisebox{-0.4pt}{\includegraphics[scale={0.08}]{transporting_b.png}}$
(6)

where
$\bar{\sigma}=\pm\frac{i}{2}[\sigma^{a},\sigma^{b}]+\delta_{ab}\sigma^{a}$.
This works for any single qubit gate with a prefactor of $\pm\pi/4$ but in the
negative case we have to reverse the appropriate signs for $\bar{\sigma}$. If
there is an arbitrary prefactor $\alpha$, this transportation identity does
not work. However, in practice, such situations will not appear in our
circuits.

Given a commuting pair of two-qubit gates acting on the same qubits they can
be either commuted through, annihilated if they differ only in their sign
prefactor in the exponential or combined into a single two-qubit gate if this
gate is considered native. Graphically we have:

$\begin{matrix}\ \ \ \ \
\raisebox{-0.4pt}{\includegraphics[scale={0.08}]{xyzdecomp6a.png}}=\raisebox{-0.4pt}{\includegraphics[scale={0.08}]{xyzdecomp6b.png}}=\raisebox{-0.4pt}{\includegraphics[scale={0.08}]{xyzdecomp6c.png}},\\\
\\\ \ \
\,\,\raisebox{-0.36pt}{\includegraphics[scale={0.08}]{xyzdecomp4a.png}}=\raisebox{-0.36pt}{\includegraphics[scale={0.08}]{xyzdecomp4b.png}}\
=\,\ \
\raisebox{-0.36pt}{\includegraphics[scale={0.08}]{xyzdecomp4c.png}}\end{matrix}$
(7)

where we grouped the two last gates into
$e^{i\alpha_{1}X_{1}Y_{2}+\alpha_{2}Y_{1}X_{2}}$, with
$\alpha_{1}=\pm\alpha_{2}$, which up to single qubit rotations corresponds to
a parameterised iSWAP gate ($e^{i\alpha(X_{1}X_{2}+Y_{1}Y_{2})}$). This way,
two commuting Pauli strings with the same prefactor (independently of their
sign) can be implemented using one two-qubit and some single-qubit gates. We
thus assume that the squared ’native’ two-qubit gate on the lower r.h.s. of
Eq. 7 can admit any combination of signs. All commuting operators fulfilling
this can be combined into one two-qubit gate if a set of three native gates is
available, i.e. the set $\\{e^{i\phi
Z_{1}Z_{2}},e^{i\phi(X_{1}Y_{2}+Y_{1}X_{2})},e^{i\phi(X_{1}Y_{2}+Z_{1}Z_{2})}\\}$.
More generally, all such gates can be implemented, up to single qubit
rotations, as a fermionic simulation (fSIM) gate [42]:

$\text{fSIM}_{ij}(\theta,\phi)=e^{i\frac{\theta}{2}(X_{i}X_{j}+Y_{i}Y_{j})+i\frac{\phi}{4}(Z_{i}+Z_{j}-Z_{i}Z_{j})}$
(8)

This gate also allows to merge three commuting Pauli strings if two of them
have a similar absolute prefactor. Such a two-qubit gate has been natively
implemented on existing superconducting platforms [51] and has been shown to
be useful for compressing fermionic operators in the context of the JW mapping
combined with a fermionic swap (fSWAP) gate network [42]. In principle, we
found that nearly all two-qubit gate pairs appearing within our circuits can
be transformed into $e^{i\alpha(X_{1}Y_{2}+Y_{1}X_{2})}$ (with gates coming
from some quartic fermionic interaction terms being notable exceptions which
generate $e^{i\alpha Z_{1}Z_{2}}$ type gates). This is why we will use its
graphical form from Eq. 7 throughout later sections instead of fSIM gates.

For anti-commuting two-qubit gates acting on the same qubits it is possible to
combine them into a single two-qubit gate if at least one of them has a
$\pm\pi/4$ prefactor in the exponential, which corresponds to the common case
encountered:

$\raisebox{-0.42pt}{\includegraphics[scale={0.08}]{xyzdecomp8a.png}}=\raisebox{-0.42pt}{\includegraphics[scale={0.08}]{xyzdecomp8b.png}}$
(9)

where $\sigma^{a}=\frac{i}{2}[\sigma^{b},\bar{\sigma}]$. Let us now show how
one can use these circuit compression identities in practice. Consider the two
operators $e^{i\alpha Y_{1}Z_{2}Z_{3}X_{4}}e^{i\alpha X_{1}Z_{2}Z_{3}Y_{4}}$:

$\begin{matrix}\raisebox{-0.46pt}{\includegraphics[scale={0.08}]{exhca.png}}=\raisebox{-0.46pt}{\includegraphics[scale={0.08}]{exhcb.png}}=\raisebox{-0.46pt}{\includegraphics[scale={0.08}]{exhcc.png}}\end{matrix}$
(10)

We have decomposed the two four-qubit operators (l.h.s.) in such a way that
most gates commute (center) and subsequently paired them into native two-qubit
gates (r.h.s.), which has reduced the depth by half. As we will see, this is a
common theme for fermionic Hamiltonians and can be proven to be optimal.
Specifically, for $\alpha=\pi/8$:

$\displaystyle e^{i\frac{\pi}{8}(X_{1}X_{2}...X_{n-1}Y_{n}+}$
${}^{Y_{1}X_{2}...X_{n-1}X_{n})}|0\rangle^{\otimes n}=$ (11)
$\displaystyle\frac{|0\rangle^{\otimes n}+|1\rangle^{\otimes n}}{\sqrt{2}}$

The implementation of two such operators together therefore gives rise to a
maximally entangled GHZ state. Following the steps in [49], this implies that
for every pair of qubits there has to be a chain of non-commuting two-qubit
gates sequentially connecting them, so a minimal depth of $n-1$ is needed for
operators with even weight and $n$ for operators with odd weight, proving the
XYZ decomposition for fermionic hoppings to be optimal in terms of depth and
in the number of two-qubit gates.

Now let us add another two operators to the example of Eq. 10, which have a
partial qubit overlap with the first two: $e^{i\alpha
Y_{1}Z_{2}Z_{3}X_{4}}e^{i\alpha X_{1}Z_{2}Z_{3}Y_{4}}e^{i\alpha
Y_{3}Z_{4}Z_{5}X_{6}}e^{i\alpha X_{3}Z_{4}Z_{5}Y_{6}}$. We obtain the circuit:

$\begin{matrix}\raisebox{-0.47pt}{\includegraphics[scale={0.08}]{exdisja.png}}=\raisebox{-0.47pt}{\includegraphics[scale={0.08}]{exdisjb.png}}=\raisebox{-0.47pt}{\includegraphics[scale={0.08}]{exdisjc.png}}\end{matrix}$
(12)

where we have been able to mutually annihilate two-qubit gates between the two
pairs, which allowed us to further parallelize gates and reduce the total
circuit depth by a further factor two.

These examples clearly demonstrate the power of our approach. However, the
overall efficiency of the formalism will differ for different fermion-to-qubit
mappings. In the following section, we will present a general strategy for
generating mappings for a square qubit layout. These mappings will then be
compared against each other in the context of these XYZ decomposition and
compression techniques.

##  3 Fermionic models on regular lattices

The most general model used to describe fermionic systems in quantum chemistry
and material science is the electronic structure Hamiltonian. Its second-
quantized form is given by:

$\mathcal{H}_{\text{ES}}=\sum_{pq}h^{pq}c^{\dagger}_{p}c_{q}^{\phantom{\dagger}}+\sum_{pqrs}h^{pqrs}c^{\dagger}_{p}c^{\dagger}_{q}c_{r}^{\phantom{\dagger}}c_{s}^{\phantom{\dagger}}$
(13)

where $p,q,r,s$ are indices of a given basis set, $c^{\dagger}$ and $c$ are
the fermionic creation and annihilation operators and $h^{pq}$, $h^{pqrs}$ are
constants called the one- and two-electron integrals, respectively. In
general, the electronic structure Hamiltonian is maximally dense, thus
containing all $N^{4}$ possible quadratic and quartic terms, where $N$ is the
total number of fermionic modes in the system. For such dense Hamiltonian, it
has been shown that that the optimal strategy is to use the Jordan-Wigner
mapping in combination with fSWAP networks [42, 43]. However, this optimality
no longer holds true for structured Hamiltonians with some degree of sparsity,
like those defined on regular lattices. Indeed, these are the systems for
which we expect our approach to be the most fruitful.

One of the most complex models falling into this class is the Hubbard-Kanamori
model (HK), designed to accurately describe the competing spin- and orbital
degrees of freedom in d- and f- electron materials. Notably, the strongly
correlated transition-metal oxides [52, 53, 54], which are of fundamental
relevance to the design of batteries [2] fall into this category. The HK
Hamiltonian reads:

$\displaystyle\mathcal{H}_{\text{HK}}=\sum_{i,j,m,\sigma}t^{ijm\sigma}c^{\dagger}_{im\sigma}\;c_{jm\sigma}^{\phantom{\dagger}}$
$\displaystyle+\sum_{i,m}U^{im}\;n_{im\uparrow}^{\phantom{\dagger}}n_{im\downarrow}^{\phantom{\dagger}}$
(14)
$\displaystyle+\\!\\!\\!\sum_{i,m<\bar{m}}\\!\\!\\!U_{1}^{im\bar{m}}\left(n_{im\uparrow}^{\phantom{\dagger}}n_{i\bar{m}\downarrow}^{\phantom{\dagger}}+n_{im\downarrow}^{\phantom{\dagger}}n_{i\bar{m}\uparrow}^{\phantom{\dagger}}\right)$
$\displaystyle+\\!\\!\\!\sum_{i,m<\bar{m}}\\!\\!\\!U_{2}^{im\bar{m}}\left(n_{im\uparrow}^{\phantom{\dagger}}n_{i\bar{m}\uparrow}^{\phantom{\dagger}}+n_{im\downarrow}^{\phantom{\dagger}}n_{i\bar{m}\downarrow}^{\phantom{\dagger}}\right)$
$\displaystyle+\\!\\!\\!\sum_{i,m<\bar{m}}\\!\\!\\!J^{im\bar{m}}\left(c^{\dagger}_{im\uparrow}c^{\dagger}_{im\downarrow}c_{i\bar{m}\downarrow}^{\phantom{\dagger}}c_{i\bar{m}\uparrow}^{\phantom{\dagger}}+c^{\dagger}_{i\bar{m}\uparrow}c^{\dagger}_{i\bar{m}\downarrow}c_{im\downarrow}^{\phantom{\dagger}}c_{im\uparrow}^{\phantom{\dagger}}\right.$
$\displaystyle\left.\quad\quad\quad\quad\quad\\!+c^{\dagger}_{im\uparrow}c^{\dagger}_{i\bar{m}\downarrow}c_{im\downarrow}^{\phantom{\dagger}}c_{i\bar{m}\uparrow}^{\phantom{\dagger}}+c^{\dagger}_{i\bar{m}\uparrow}c^{\dagger}_{im\downarrow}c_{i\bar{m}\downarrow}^{\phantom{\dagger}}c_{im\uparrow}^{\phantom{\dagger}}\right)$

where $n$ is the number operator, the indices $i$, $j$ denote lattice sites
($i\neq j$), the indices $m$, $\bar{m}$ label orbitals and there are two spins
$\sigma\in\\{\uparrow,\downarrow\\}$ per orbital. The terms with prefactor $t$
correspond to fermionic hopping, $U$, $U_{1}$ and $U_{2}$ terms are intra- and
inter-band density-density interactions, and $J$ terms represent pair-hopping
and spin-exchange interaction terms. Following common practice [54], we only
allow for hopping to occur between neighboring lattice sites.

One can trivially reduce this Hamiltonian to the single-orbital Fermi-Hubbard
model (FH) which is believed to capture the physics of copper-based high-
temperature superconductors. The model is obtained from Eq. 3 by setting the
number of orbitals to one and therefore suppressing all inter-orbital
interaction terms, i.e. $U_{1}^{im\bar{m}}=U_{2}^{im\bar{m}}=J^{im\bar{m}}=0$.
Finally, we set $t^{ijm\sigma}\equiv-t^{ij}$ with $t^{ij}\neq 0$ only if
$i$,$j$ are considered neighbors on the fermionic lattice. The resulting
Hamiltonian has the form:

$\displaystyle\mathcal{H}_{\text{FH}}$
$\displaystyle=-\sum_{i,j,\sigma}t^{ij}c^{\dagger}_{i\sigma}c_{j\sigma}^{\phantom{\dagger}}+U\sum_{i}n_{i\uparrow}^{\phantom{\dagger}}n_{i\downarrow}^{\phantom{\dagger}}$
(15)

Note, that in what follows we omit any chemical potential terms
($\mu^{\sigma}=t^{iim\sigma}$) of the Fermi-Hubbard Hamiltonian for simplicity
reasons. Such terms could be trivially added to the Hamiltonian without
requiring any changes to our approach. The FH model is extremely challenging
to classical computational methods and, with the exception of some specific
regimes, remains largely unsolved.

If we further allow for only one spin type, and consequently also set $U=0$,
we recover the one-spin tight-binding model Hamiltonian (TB):

$\displaystyle\mathcal{H}_{\text{TB}}$
$\displaystyle=-\sum_{i,j}t^{ij}c^{\dagger}_{i}c_{j}^{\phantom{\dagger}}$ (16)

This model is relatively trivial to solve, but is nevertheless useful to
investigate it for demonstration purposes, as will be done in section 5.1.

In principle, one can define all of the aforementioned models on any fermionic
connectivity graph, but here we only focus on regular, two-dimensional
lattices. The most common choice is the square lattice with nearest-neighbor
(NN) and optionally next-nearest-neighbor (NNN) connectivity. However, other
geometries have also been extensively studied in literature, notably the
triangular [55, 56], honeycomb [57, 58] and Kagome [59] lattices. In this
paper, we consider eight different geometries as shown in Fig. 1. Rather than
investigating optimal mappings for each of these lattices individually, we
embed each of them into a square lattice layout whilst allowing for higher-
neighbor connectivity (and the corresponding inter-site hopping terms). This
way it is sufficient to only consider strategies of mapping a square fermionic
lattice to a square qubit layout. We do not consider one- or three-dimensional
fermionic models here, as for one-dimensional systems the Jordan-Wigner
mapping is already optimal and in the case of three-dimensional systems one
faces the additional difficulty of accommodating the third dimension on a two-
dimensional lattice. The best approach in this case is to treat it on equal
footing with orbital degrees of freedom [32], as we will show in detail in
further sections.

Figure 1: Different fermionic lattices embedded into a square lattice geometry
by allowing for higher-neighbour connectivity. Nodes represent fermionic
lattice sites and links indicate the existence of hopping terms between two
sites.

##  4 Local fermion-to-qubit mappings for square qubit layouts

In this section we will investigate how to efficiently transform fermionic
Hamiltonians into spin Hamiltonians containing operators that act on a set of
qubits $Q$ which can be implemented on a quantum device:

$\displaystyle\mathcal{H}_{\text{Q}}=\sum_{i}a_{i}\bigotimes_{j\in
Q}\,\tilde{\sigma}_{i}^{j},$ (17)

where, $a_{i}$ are constants and $\tilde{\sigma}_{i}^{j}\in\\{I,X,Y,Z\\}$. The
antisymmetric nature of fermionic systems is encoded in the state wavefunction
when working in first quantization. In second quantization this is already
included in the fermionic operators used to describe the Hamiltonian of the
system. This poses a challenge for transforming the antisymmetric operators of
fermionic Hamiltonians into Pauli qubit operators that have no native
antisymmetry. This transformation should preserve the locality of the
fermionic interactions to avoid increases in the scaling of the number of
gates or the circuit depth.

The most convenient approach for creating local fermion-to-qubit mappings
(i.e. those where the operator weight is constant with respect to the system
size) consists of defining the vertex $V_{i}$ and edge $E_{ij}$ operators for
every fermionic mode $i$ and pair of modes $(i,j)$. To facilitate the
conversion between the two types of operators, we can additionally introduce
the Majorana fermionic operators
$\gamma_{i}=c_{i}^{\phantom{\dagger}}+c_{i}^{\dagger}$ and
$\bar{\gamma}_{i}=\frac{1}{i}(c_{i}^{\phantom{\dagger}}-c_{i}^{\dagger})$
[35]. This allows us to compose the edge and vertex operators as
$E_{ij}=-i\gamma_{i}\gamma_{j}=-i(c_{i}c_{j}+c^{\dagger}_{i}c^{\dagger}_{j}+c_{i}c^{\dagger}_{j}+c^{\dagger}_{i}c_{j})$,
$V_{i}=-i\gamma_{i}\bar{\gamma}_{i}=c_{i},c^{\dagger}_{i}-c^{\dagger}_{i}c_{i}$.
In order to be compatible with the fermionic anticommutation relations, the
edge and vertex operators must themselves satisfy the following relations:

$\displaystyle\quad\quad\\{E_{ij},V_{i}\\}=\\{E_{ij},E_{jk}\\}=0$ (18)
$\displaystyle[E_{ij},E_{kl}]=[E_{ij},V_{k}]=[V_{i},V_{j}]=0$

for indices $i\neq j\neq k\neq l$. This means that edges must anticommute with
vertices which they are incident on, two edges must anticommute if they share
a vertex, and all other combinations of two operators must commute.
Additionally, the following relation must be fulfilled for any closed loop of
edges:

$i^{(|p|-1)}\prod_{j}^{|p|-1}E_{p_{j},p_{j+1}}=\mathbb{1}$ (19)

where $p=\\{p_{1},p_{2},...\\}$ forms a closed path. The condition in Eq. 19
is not generally fulfilled since most combinations for edge and vertex
operators will create a non-trivial operator $S_{p}$ instead of the identity
$\mathbb{1}$. However, one can solve this issue by initializing the quantum
device in the common eigenspace $U=\bigcap_{p}U_{p}$ of each product’s $+1$
eigenspace $U_{p}$ and thus $S_{p}\sim\mathbb{1}$.

As an extra condition for the mapping to be correct, one has to define inverse
edge operators as $E_{ji}=-E_{ij}$, which we do by assigning one sign
prefactor to the ordering $(i,j)$ and the opposite sign to the inverse one
$(j,i)$. Edge operators may be constructed from other existing edges by using
the _composite_ rule:

$E_{ik}=iE_{ij}E_{jk}$ (20)

This means that it is not necessary to define a unique edge for every
connected fermionic mode pair, as long as a composite edge path exists between
them. An alternative strategy to composing edge operators is to instead swap
fermionic modes between adjacent vertices. This can be done using an fSWAP
operator between two modes $(i,j)$, which is defined as [35]:

$\operatorname{fSWAP}_{i,j}=e^{i\frac{\pi}{4}V_{i}}e^{i\frac{\pi}{4}V_{j}}e^{\frac{\pi}{4}(E_{ij}V_{j}+V_{i}E_{ij})}$
(21)

Given these definitions we can now translate the quadratic and quartic
fermionic operators occurring in the models of Sec. 3 into products of edge
and vertex operators:

$n_{j}^{\phantom{\dagger}}\rightarrow\frac{1}{2}(1-V_{j})$ (22)
$n_{j}^{\phantom{\dagger}}n_{k}^{\phantom{\dagger}}\rightarrow\frac{1}{4}(1-V_{j}-V_{k}+V_{j}V_{k})$
(23)
$c^{\dagger}_{j}c_{k}^{\phantom{\dagger}}\rightarrow\frac{i}{4}(1+V_{k}-V_{j}-V_{j}V_{k})E_{jk}$
(24) $\displaystyle
c^{\dagger}_{j}c^{\dagger}_{k}c_{l}^{\phantom{\dagger}}c_{m}^{\phantom{\dagger}}\to-\frac{1}{16}(1+V_{m}-V_{k}+V_{l}-V_{j}$
(25) $\displaystyle-
V_{k}V_{m}+V_{l}V_{m}+V_{j}V_{k}-V_{k}V_{l}-V_{j}V_{l}-V_{j}V_{m}$
$\displaystyle-
V_{k}V_{l}V_{m}+V_{j}V_{k}V_{m}-V_{j}V_{l}V_{m}+V_{j}V_{k}V_{l}$
$\displaystyle\quad\quad\quad\quad+V_{j}V_{k}V_{l}V_{m})E_{jl}E_{km}$

where $j\neq k\neq l\neq m$. One can notice that the resulting formulation
involves significantly more terms than the original fermionic operator.
However, most of these cancel out once multiple symmetry-related fermionic
operators are combined:

$c^{\dagger}_{j}c_{k}^{\phantom{\dagger}}+c^{\dagger}_{k}c_{j}^{\phantom{\dagger}}\rightarrow\frac{i}{2}(V_{k}-V_{j})E_{jk}$
(26) $\displaystyle
c^{\dagger}_{j}c^{\dagger}_{k}c_{l}^{\phantom{\dagger}}c_{m}^{\phantom{\dagger}}+c^{\dagger}_{j}c^{\dagger}_{l}c_{k}^{\phantom{\dagger}}c_{m}^{\phantom{\dagger}}+c^{\dagger}_{l}c^{\dagger}_{m}c_{j}^{\phantom{\dagger}}c_{k}^{\phantom{\dagger}}+c^{\dagger}_{m}c^{\dagger}_{k}c_{l}^{\phantom{\dagger}}c_{j}^{\phantom{\dagger}}\to$
$\displaystyle-\frac{1}{4}(V_{j}V_{k}-V_{j}V_{l}+V_{l}V_{m}-V_{m}V_{k})E_{jl}E_{km}$
(27)

which can be expressed using other edges since $E_{jl}E_{km}=-E_{jm}E_{kl}$.
Given a set of edges and vertices from a fermionic connectivity graph, one has
to define their matching (Pauli string) operators acting on hardware qubits
which satisfy all the commutation relations from Eq. 18. The choice of these
operators is not unique and will influence the performance of the underlying
algorithm, quantified either in terms of the qubit-to-mode ratio, the gate
count, the circuit depth or the error detecting and correcting properties of
the mapping.

(a)

(b)

Figure 2: Fermion-to-qubit mappings for square qubit layout: a) Unit cells
showing different horizontal sequences of physical qubits (squares, P) and
ancilla qubits (circles, A) b) Vertex operators and all possible edge
operators defined in the horizontal and vertical directions. The X/Y symbols
indicate two possible Pauli strings for a given edge, chosen depending on the
particular mapping instance.

The family of mappings proposed in this paper uses common structures for the
edge and vertex operators (see Fig. 2). For any given mapping we define a unit
cell, which consists of a sequence of physical (P) and ancilla (A) qubits.
Each unit cell is repeated horizontally and vertically to cover the qubit
lattice. We define physical qubits as those on which a fermionic mode exists,
meaning that a vertex operator will be acting on it with a Pauli Z operator.
Additionally, we introduce ancilla qubits, whose purpose is to resolve
(anti-)commutation relations between edge operators. In cases where the qubit
lattice is incommensurate with the unit cell one can partially cover
additional qubits with the given unit cell pattern. In what follows, we will
not consider the effects of finite lattices as ultimately they will not
influence the relative performance, in terms of circuit depth, of the mappings
considered.

The qubit-to-mode ratio of a mapping is given by $r=(N_{P}+N_{A})/N_{P}$,
where $N_{P}$ and $N_{A}$ are the number of physical and ancilla qubits,
respectively. We allow a unit cell to contain any number of physical qubits
followed by up to two ancilla qubits. We found that whilst two consecutive
ancilla qubits can improve the performance of the mapping by allowing for
additional degrees of parallelism, a third ancilla does not yield further
improvements in this respect but simultaneously worsens other metrics of the
mapping. Similarly, we found that allowing for more than two consecutive
physical qubits only leads to an improved qubit-to-mode ratio whilst
increasing maximum weights of operators as well as the resulting circuit
depths. Nevertheless, it is important to note that such mappings might be a
viable option if the number of qubits is the limiting factor, as one can, in
principle, come arbitrarily close to a ratio of $r=1$. In this limit we would
recover the JW mapping embedded in a two-dimensional lattice using a snake
pattern consisting solely of two-qubit horizontal edge string and two-qubit
vertical edges to connecting them at the ends.

The above physical-ancillary pattern restrictions leave us with four possible
choices of unit cells: PPA, PA, PPAA and PAA (Fig. 2.a)), which we will
consider for mapping particular models in the next sections. We can now define
a common set of edge operators for these four unit cells (Fig. 2.b)). This
includes horizontal edges, which can act on between two, three and four
qubits, depending on the number of ancilla qubits connecting two physical
ones. Vertical edges act on four qubits and can only be defined when a pair of
ancilla qubits is adjacent to a pair of physical qubits. In cases where a
vertical edge exists on the boundary of the lattice and no other vertical edge
shares physical qubits with it (i.e. in the Honeycomb lattice of Fig. 1), one
can omit the ancilla qubits in the vertical edge operator and reduce it to the
two-qubit form shown in Fig. 2.b). All additional higher-order edges occurring
in the lattices of Fig. 1 are generated using the composite rule of Eq. 20.
Further, consecutive edges in the horizontal direction must act on their
common physical qubit with different Pauli operators, which means that we
alternate between the set of operators {XX, XZX, XZZX, XXYX, XYXX} and the set
{YY, YZY, YZZY, YXYX, YYXY}, resulting in a pool of operators of weights 2-4.
This, together with the definitions of the edge operators in Fig. 2,
guarantees that all commutation relations of Eq. 18 are satisfied.

Note that the PA mapping is equivalent to the Verstraete-Cirac mapping [37]
embedded in a square qubit topology (the authors of the original paper
considered two superimposed square grids of qubits, whilst the mapping has
been embedded into a non-square topology in Ref. [31]). The remaining three
mappings have, to our knowledge, not been previously considered.

Figure 3: a) Fermionic connectivity of the 3-by-3 NNN square lattice, arranged
into parallel groups. b) The corresponding numbered qubits in a 6-by-3 square
qubit layout using the PA mapping. Square shapes represent physical qubits,
circles are ancilla qubits. c) Vertex (black) and edge operators for this
mapping. Red and blue strings represent horizontal edge operators while green
and purple strings depict vertical ones. The letters X,Y,Z indicate the Pauli
acting on a given qubit. Dashed qubit connections are not required. d) All
operators obtained from applying the mapping to the hopping terms of the
tight-binding model, grouped according to a). $X/Y$ and $Y/X$ terms indicate
two operators of the form $X\dots Y$ and $Y\dots X$. Diagonal terms are
obtained from the composite edge rule. e) Resulting quantum circuit for the
execution of a single Trotter step for the Hamiltonian, obtained using the
XYZ-decomposition, and assuming the two-qubit gate $e^{i\alpha(XY+YX)}$ is
native. Further circuit compression can be applied within grey shaded areas of
the circuit.

As an example of how to use the XYZ formalism efficiently in conjunction with
these mappings, let us consider PAA, which contains four-qubit edge operators
both in the horizontal and vertical directions. The corresponding nearest-
neighbor fermionic hopping terms can then be efficiently transformed in
Hermitian conjugate pairs into four-qubit operators using Eq. 26. Applying the
tricks we introduced in the context of the XYZ formalism we can compress the
circuit for the horizontal hopping operators,
$t_{h}=e^{i\alpha(X_{1}Z_{2}Z_{3}Y_{4}-Y_{1}Z_{2}Z_{3}X_{4})}$ and the
vertical hopping operators
$t_{v}=e^{i\alpha(X_{1}Y_{2}X_{3}Y_{4}-Y_{1}Y_{2}X_{3}X_{4})}$, into depth
three:

$\raisebox{-0.47pt}{\includegraphics[scale={0.08}]{spinful_operatorsb1.png}}=\raisebox{-0.47pt}{\includegraphics[scale={0.08}]{spinful_operatorsb2.png}}\
,\quad\raisebox{-0.47pt}{\includegraphics[scale={0.08}]{spinful_operatorsc1.png}}=\raisebox{-0.47pt}{\includegraphics[scale={0.08}]{spinful_operatorsc2.png}}$
(28)

where the two outside qubits are physical ones and the two inside qubits are
ancillas. Similarly, all density-density fermionic terms of the form of Eq. 23
can be straightforwardly transformed into circuits with a single two-qubit
gate:

$\raisebox{-0.43pt}{\includegraphics[scale={0.085}]{spinful_operatorsa1.png}}=\raisebox{-0.43pt}{\includegraphics[scale={0.085}]{spinful_operatorsa2.png}}=\raisebox{-0.43pt}{\includegraphics[scale={0.085}]{spinful_operatorsa3.png}}=\raisebox{-0.43pt}{\includegraphics[scale={0.085}]{spinful_operatorsa4.png}}$
(29)

The density-density term and the hopping terms between two modes can, in fact,
be implemented as efficiently as the hopping terms alone, even if the modes
are connected through a number of ancillas qubits:

$\begin{matrix}\ \ \ \
\raisebox{-0.43pt}{\includegraphics[scale={0.105}]{hopwithint1.png}}=\raisebox{-0.43pt}{\includegraphics[scale={0.105}]{hopwithint2.png}}=\quad\quad\\\
\\\
=\raisebox{-0.43pt}{\includegraphics[scale={0.105}]{hopwithint3.png}}=\raisebox{-0.43pt}{\includegraphics[scale={0.105}]{hopwithint4.png}}\end{matrix}$
(30)

Here, the central gate on the r.h.s. represents an interaction of the form
$e^{i\theta(X_{1}Y_{2}+Y_{1}X_{2})+i\phi Z_{1}Z_{2}}$, which is an fSIM gate
up to single qubit rotations. This shows that using the XYZ decomposition for
jointly implementing hopping and density-density operators between two
physical modes is optimal.

One can also efficiently compress the four additional quartic pair-hopping and
spin-exchange terms which appear in the HK model of Eq. 3 into a depth-five
circuit by using Eq. 4:

$\raisebox{-0.48pt}{\includegraphics[scale={0.08}]{spinful_operatorse1.png}}=\raisebox{-0.49pt}{\includegraphics[scale={0.08}]{spinful_operatorse2.png}}$
(31)

Here, we assume that the four modes correspond to pairs of spins
$\\{\uparrow,\downarrow\\}$ on two arbitrary orbitals $\\{1,2\\}$ within this
operator are ordered (up to spin and orbital relabeling) as:
$1\\!\\!\uparrow$, $2\\!\\!\downarrow$, $1\\!\\!\downarrow$,
$2\\!\\!\uparrow$. The shaded area in Eq. 31 indicates that additional
compression is possible within contained gates. Note that whilst for hopping
terms it was more efficient to combine Hermitian conjugates, in the case of
the quartic terms we first combine operators between the two groups (pair
hopping and spin exchange) as these differ by exactly two Paulis acting on the
first and last qubit, which leads to maximal cancellations. Beyond the HK
model, the density-assisted hoppings operators [60] of the form
$(c^{\dagger}_{i,\downarrow}c_{j,\downarrow}+c^{\dagger}_{j,\downarrow}c_{i,\downarrow})n_{i,\uparrow}$
can also be straightforwardly implemented with the presented techniques.

We can also consider the computational cost of implementing an instance of the
fSWAP operator within our mappings. From Eq. 21 it is evident that the first
two exponential terms correspond to single-qubit gates and the last has
similar structure to a hopping term together with its Hermitian conjugate (see
Eq. 26). This term can be implemented especially efficiently, in depth one, if
we swap modes between vertices connected by a two-qubit edge (meaning two
consecutive P qubits). We obtain the following circuit for such fermionic swap
operators:

$\raisebox{-0.42pt}{\includegraphics[scale={0.08}]{spinful_operatorsd1.png}}=\raisebox{-0.42pt}{\includegraphics[scale={0.08}]{spinful_operatorsd2.png}}=\raisebox{-0.3pt}{\includegraphics[scale={0.12}]{spinful_operatorsd3.png}}$
(32)

Note that fSWAPs are reversible and sometimes this can be useful when
combining it with other operators within a specific ordering. Fermionic swap
operators can also be implemented between modes connected through ancilla
qubits, but at a higher cost and we thus refrain from using them in this work.
Nonetheless, it can be noted that through Eq. 30, they can be optimally
combined with a pair of hopping operators together with the density-density
operator by shifting the angle of the hoppings by $\pi/4$. This extends the
work in [42], where they used the fSIM gate for executing together these three
types of operators between neighboring modes on JW chain.

fSWAPs are especially practical for models containing multiple modes per
fermionic lattice site. They allow us to keep the same general unit cell
structure of our mappings derived for one mode per site by replacing every
physical qubit P by a chain of physical qubits corresponding to all modes
defined on the given lattice site (this approach has been first introduced in
Ref.[32]). This chain will locally have the properties of a JW mapping on each
individual site and all modes except for those at the ends will only be
connected by horizontal edges. These modes can then be transported through the
chain using a fermionic SWAP network, which due to the JW nature of the chain
leads to much lower circuit depths compared to applying the composite rule.

As an example, consider the Fermi-Hubbard model with two spin-modes per site.
As explained, we arrange the modes next to each other horizontally
($\textit{P}_{\uparrow}$,$\textit{P}_{\downarrow}$) and connect them through a
two-qubit edge operator $E_{\uparrow\downarrow}=X_{\uparrow}X_{\downarrow}$
and the vertex operators of both modes remain single-qubit Z’s. In cases where
one has more than one mode per lattice site it is necessary to implement an
fSWAP network, which alternates fSWAP operations between even ($2k$, $2k+1$)
and odd ($2k+1$, $2k+2$) pairs of modes associated to the site, as done for
normal qubit swap networks [61]. This ensures that every mode travels to both
ends of the chain and then interacts with similar modes of neighboring lattice
sites after applying only linear number of swap layers.

In order to additionally ensure that all intra-site Hamiltonian terms from Eq.
3 can be implemented at some point in the network, the relative position of
modes is equally important. The described swap network only ensures that every
pair of modes is adjacent at some point, but this requirement cannot be
satisfied for every group of four modes. Luckily, the quartic spin-exchange
and pair-hopping terms always act on two pairs of spin-modes within two
orbitals. With this restriction it is indeed possible to make all such groups
adjacent within the swap network if we ensure that pairs of spins of an
orbital always move in the same direction. If we consider the example of four
orbitals, then this can be ensured by arranging the modes in the pattern:
$1\\!\\!\uparrow$, $2\\!\\!\downarrow$, $1\\!\\!\downarrow$,
$2\\!\\!\uparrow$, $3\\!\\!\uparrow$, $4\\!\\!\downarrow$,
$3\\!\\!\downarrow$, $4\\!\\!\uparrow$. It is easy to see that any two spin
modes of an orbital will at every step of the swap network be at most one
qubit apart. The reverse ordering of the spins between right-moving (1,3) and
left-moving (2,4) pairs further guarantees that all $J$-terms can be
implemented in depth five, as shown using Eq. 31.

Finally, we would like to point out that the edge operators in the presented
mappings generate a group of stabilizers $S_{p}$ with weights ranging from 6
to 8, depending on the particular choice of mapping. These, similarly to the
mappings in Ref. [40], allow for partial single-qubit error detection, but not
for error correction.

##  5 Simulating fermionic systems

Many different algorithms for studying the properties of fermionic models on
quantum computers can be found in literature. The most straightforward example
consists of time-evolving a system under a Hamiltonian of the form
$H=\textstyle{\sum_{j}h_{j}}$ using the Suzuki-Trotter decomposition formula
[62]:

$\displaystyle e^{-iHt}=\lim_{\delta_{t}\to
0}\prod_{k}^{t/\delta_{t}}\prod_{j}e^{-ih_{j}\delta_{t}}$ (33)

where $t/\delta_{t}\in\mathbb{Z}$. For any non-zero value of $\delta_{t}$ this
relation becomes approximate and the quality of the result will depend on the
smallness of $\delta_{t}$. Following this approach, there have been proposed
additional techniques involving more sophisticated summation formulae [63, 64,
65, 66, 67]. Many other quantum algorithms developed for fermionic models also
require implementing sequential evolution operators [68, 69, 70, 71, 72, 73,
74, 75], typically equivalent to a single Trotter step or some of it’s parts.

In this section we aim investigate the performance of our four fermion-to-
qubit mappings (PPA, PA, PPAA, PAA) used together with the XYZ-decomposition
formalism, and evaluated by the depth of the quantum circuits obtained for a
single Trotter step ($k=1$ in Eq. 33). When counting gates (and depth) we
consider single-qubit gates as a free resource and fSIM two-qubit gates as
native. We study three different system classes with varying complexity: the
tight-binding model (TB), the single-band Fermi-Hubbard model (FH) and the
multi-orbital Hubbard-Kanamori model (HK). In the case of TB and FH models we
consider all eight lattices shown in Fig. 1. For geometries which are not
symmetric with respect to 90 degree rotations, such as the honeycomb and
Kagome lattices, we consider both orientations. Note that the Kagome fermionic
lattice can also be embedded into a square qubit layout by without third-
neighbor connections but at the cost of including idling qubits. We therefore
consider both options for its embedding and only report the minimal depth
found. To obtain a better idea of the performance of our fermion-to-qubit
mappings, we benchmark them against the state-of-the-art mapping from Derby
and Klassen (DK) [40] and its later modification for multi-orbital models
[32]. The details for this mapping and details on how one can embed it into a
square qubit layout can be found in Appendix B.

###  5.1 Tight-binding model

Let us start by considering the TB model. In Fig. 3 we show a worked-out
example for this model on an NNN square lattice using the PA mapping. We
investigate a 3-by-3 fermionic NNN square lattice with a single mode per
lattice site and open boundary conditions processed using the PA mapping. The
total circuit depth, however, is independent of system size. Fig. 3.a) shows
the fermionic lattice sites together with all existing connections between
them. Using the PA mapping we can port this lattice onto a square 18-qubit
layout (6-by-3, see Fig. 3.b)). We then proceed by defining all vertices,
horizontal and vertical edges in Fig. 3.c) and, using Eq. 26 the operators
corresponding to hopping terms of the NNN square lattice combined with their
Hermitian conjugates (Fig. 3.d)). Diagonal hopping operators can be
constructed using the composite rule from a vertex and two edges (one
horizontal and one vertical). In the PA mapping, all diagonal hopping terms
act on four qubits. For a particular quantum algorithm one can now efficiently
decompose these qubit operators using the formalism which we introduced in
Sec. 2. The resulting circuit of one step in the Suzuki-Trotter decomposition
of the time-evolution for the TB Hamiltonian on this lattice is shown in Fig.
3.e).

As previously discussed, we find that the strategy generating the shallowest
circuits is to first combine the operator pairs generated from mapping the sum
of Hermitian conjugate hopping terms, in the spirit of Eq. 10, which
immediately reduces the total depth by a factor two. In Fig. 3.e) we have also
rotated (with single-qubit gates) all operators to be of the form $\sigma
Z\dots Z\bar{\sigma}$, guaranteeing that all two-qubit operators are of the
form $e^{i\alpha(XY+YX)}$, which in this case would suffice as a native gate.
Parallel operators acting on distinct sets of qubits can then be collected
into layers and applied simultaneously.

The order in which these layers are applied will also affect the circuit
depth, as additional compression can be achieved between the layers. In Fig.
3.e) we have not performed such compression between layers, but rather
indicated, by shaded areas, where this is possible. The two cases allowing
additional compression that occur in Fig. 3.e) are of one of the two forms:

$\begin{matrix}\raisebox{-0.42pt}{\includegraphics[scale={0.08}]{xyzdecomp10a.png}}=\raisebox{-0.42pt}{\includegraphics[scale={0.08}]{xyzdecomp10b.png}},\\\
\\\
\raisebox{-0.42pt}{\includegraphics[scale={0.08}]{xyzdecomp11a.png}}=\raisebox{-0.42pt}{\includegraphics[scale={0.08}]{xyzdecomp11b.png}}\end{matrix}$
(34)

In such cases one will, through the application of Eq. 9, produce two-qubit
gates of the form $e^{i\alpha ZZ}$, which same as the $e^{i\alpha(XY+YX)}$
gate can be implemented using a single fSIM gate.

Lattice type | PPA | PA | PPAA | PAA | DK
---|---|---|---|---|---
Honeycomb | 8 | 7 | 6 | 8 | 8
Square | 12 | 9 | 8 | 10 | 10
Shastry-Sutherland | 19 | 12 | 16 | 16 | 13
Kagome | 14 | 14 | 13 | 13 | 17
Triangular | 26 | 14 | 26 | 19 | 23
Checkerboard | 25 | 15 | 24 | 22 | 15
Tetrakis | 23 | 18 | 20 | 23 | 22
NNN square | 33 | 18 | 31 | 25 | 29
Table 1: Single Trotter step circuit depths for the tight-binding model on
different lattice geometries, counted in native two-qubit gates (fSIM), is
shown for multiple fermion-to-qubit mappings. Circuits were obtained using the
XYZ formalism. Lowest depths per geometry are indicated in bold.

For circuits which consist of multiple Trotter steps one can compress even
further between steps. The most straightforward way is by inverting the order
of the operators, which is equivalent to the second-order Suzuki-Trotter
formula. This allows to merge two parallel layers of hopping operators on the
boundary between two Trotter steps (i.e.
$e^{-ih_{1}\delta_{t}}e^{-ih_{2}\delta_{t}}e^{-ih_{2}\delta_{t}}e^{-ih_{1}\delta_{t}}=e^{-ih_{1}\delta_{t}}e^{-i2h_{2}\delta_{t}}e^{-ih_{1}\delta_{t}}$)
and save one such layer for each additional Trotter step. In this case it can
pay off to strategically position the most computationally expensive layer at
the boundary.

Our results for the TB model are summarized in Table 1. We find that a single
Trotter step can be simulated with circuit depths ranging from 6 for the
honeycomb lattice to 17 for the NNN square lattice which contains all possible
next-nearest neighbor hopping terms. For lattices without NNN connectivity
(honeycomb and square) we found the PPAA mapping to perform fine, whilst for
all the lattices the PA mapping proved optimal. Surprisingly, the state-of-
the-art DK mapping performed worse for all geometries we have considered,
despite generating lower-weight qubit operators for the vertical hopping
operators (see a comparison of operator weights between mappings in Table. 4
of Appendix B). here, for every further Trotter step, one can reduce the
circuits by an additional two (honeycomb) to three (NNN square) depth.

One might be tempted to argue that for small system sizes, the JW mapping
performs better than the results we obtained for our local mappings. However,
in Appendix C we show that even for a 3-by-3 lattice, arguably the smallest
truly two-dimensional case, the JW mapping already performs worse than the
local mappings considered here.

By itself, the tight-binding Hamiltonian does not present a challenge to
classical methods. Other than for benchmarking purposes, it is however also
interesting with respect to mapping the Fermi-Hubbard model. Namely, if one
would allow for a quantum processor with two superimposed square qubit
lattices, one for each spin, all of the hoppings in the FH model could be
applied in parallel. Then, one only needs to treat the on-site interaction
term between pairs of spins on the same sites which, within this qubit layout,
could be immediately implemented with one layer of two-qubit gates (this
approach was previously used in Ref. [31]). This means that such circuits for
the FH model are always only one deeper than the results from Table 1. We
found these circuit depths to always be slightly superior to using mappings
involving only one square qubit layout, which will be discussed in the next
section. On the other hand, such qubit topologies are rather challenging to be
implemented in most quantum computing platforms (with trapped-ion QPUs with
all-to-all connectivity being the notable exception).

###  5.2 Fermi-Hubbard model

Let us now consider how to map the single-band FH model with two spin-types
per lattice site to a square qubit layout. For the case of two spin modes per
lattice site, we can split all physical qubits
$\textit{P}\rightarrow\textit{P}_{\uparrow}\textit{P}_{\downarrow}$. This way,
the PA pattern becomes $\textit{P}_{\downarrow}\textit{P}_{\uparrow}A$ and PAA
becomes ${P}_{\uparrow}\textit{P}_{\downarrow}AA$ (we will keep using PA and
PAA nomenclature for simplicity). Further, we keep the ordering for vertically
neighboring lattice sites the same, but switch it between horizontally
neighboring ones, so that a site with pattern
$\textit{P}_{\uparrow}\textit{P}_{\downarrow}$ is followed by a site with the
pattern $\textit{P}_{\downarrow}\textit{P}_{\uparrow}$. This allows for
similar modes to interact across neighboring lattice sites. For this model, we
will not consider the PPA and PPAA mappings, as we found they generally
perform worse in comparison to PA and PAA due to their lack of vertical
connectivity.

The general strategy for all mappings involving multiple modes per lattice
site is to devise a fermionic swap network within each lattice site
individually and perform all applicable inter-site interactions at each step
of the network much as was shown for the tight-binding Hamiltonian. In the
case of two modes per lattice site this network consists of a single fSWAP
operation.

Figure 4: a) Fermionic connectivity graph for the Fermi-Hubbard model on a
3-by-4 square lattice. b) Hopping operators (blue, purple, green and red) and
the density-density operator (yellow) for the model using the PAA mapping
where circles are ancilla qubits and filled (empty) squares represent spin up
(down) physical qubits. The X/Y and Y/X symbolize two distinct terms. Black
arrows indicate qubits between which fSWAPs operations are required. c)
quantum circuit for a single Trotter step of the model, restricted to a single
orbital with two spin modes ($1\\!\\!\uparrow$, $1\\!\\!\downarrow$). Color
scheme corresponds to a). Further circuit compression can be applied within
grey shaded areas of the circuit. Dashed qubit connections are not required.
The total circuit depth is 11.

Let us consider the example of the FH model on a NN square lattice and using
the PAA mapping, as shown in Fig. 4. In a) we see the qubit layout containing
one spin up and down physical qubit per lattice site. In b) we show the
correspond qubit layout together with examples of hopping and density-density
interaction operators. Spin-up modes correspond to full squares and spin-down
modes to empty squares. In c) we show the full circuit from the point of view
of a single fermionic site. Given translational invariance, the circuits for
all other lattice sites will be identical. We first apply one horizontal
hopping operator pair for each mode and combine it with both of the vertical
hoppings pairs in order to reduce their combined depth to seven. Next, we
apply the density-density operator, which can be combine with an fSWAP between
the two modes into a single fSIM gate. Finally, we can implement the remaining
horizontal hopping terms in the second direction which were not immediately
available before we exchanged the two spin modes in depth three. The total
circuit depth is therefore 11 per Trotter step. Note, that if only
$e^{i\phi(XY+YX)}$ and $e^{i\phi ZZ}$ native gates were available, the depth
would only increase by one, due to the combined density-density and fSWAP
operators. Generally, this is the highest additional gain from using native
fSIM gates instead of $e^{i\phi(XY+YX)}$ and $e^{i\phi ZZ}$ that can be
obtained for this model.

In Table 2 we provide the single Trotter layer depths for the FH Hamiltonian
on various lattices and for the PA and PAA mappings, as well as the DK mapping
serving as a benchmark. We use the XYZ formalism for all mappings and optimize
over all possible lattice orientations, decompositions and orderings of
operators with the Trotter later. The PAA mapping produces the lowest possible
depths for any lattice considered, the PA mapping is slightly worse throughout
and, surprisingly, the DK mapping produces clearly inferior results with up to
2.5 times longer circuits. We believe that this is due to the lack of
parallelism, which is a consequence of the lower number of ancilla qubits. The
additional parallelism and cancellations within PAA bring the depths for the
FH model always below double that of the TB model, which one would naively
expect as a lower bound. Remarkably, we find that for the honeycomb and square
lattices the resulting depths are 9 and 11, respectively, which is only worse
by tree compared to the TB model. Additional reductions of up to three can be
achieved from cancellations between Trotter layers.

It is common in fermion-to-qubit mapping literature to assume two square
layers of qubits (one for each spin) one on top of each other connected one-
to-one in the vertical axis as the qubit topology [31]. Assumming that
topology, we can implement the hoppings for spin up and down at the same time
as in our TB model and then the density-density terms in depth one. This
implies a depth for the FH model equal to the TB model plus one, i.e. depth
$7$ for a honeycomb lattice and $9$ for a square lattice.

We can compare the latter results with a trivial lower bound for a FHM on a
fermionic lattice with a connectivity graph of maximum degree $k$. In that
case, the depth of any realization of a Trotter step cannot be less than
$k+1$, since it will be limited by the maximally connected site, whose
physical qubit has to interact through mutually anticommuting gates with $k$
other physical qubits of the same spin (through hopping operators) and $1$
physical qubit with the opposite spin (through the density-density term). This
yields a naive minimal depth of $4$ for the honeycomb lattice and $5$ for the
square grid, suggesting that our approach is at least close-to optimal.

Lattice type | PA | PAA | DK
---|---|---|---
Honeycomb | 11 | 9 | 17
Square | 15 | 11 | 21
Shastry-Sutherland | 21 | 19 | 24
Kagome | 25 | 22 | 36
Triangular | 27 | 22 | 46
Checkerboard | 27 | 22 | 30
Tetrakis | 25 | 23 | 46
NNN square | 35 | 30 | 74
Table 2: Single Trotter step circuit depths for the Fermi-Hubbard model on
different lattice geometries, counted in native two-qubit gates (fSIM), is
shown for multiple fermion-to-qubit mappings. Circuits were obtained using the
XYZ formalism. Lowest depths per geometry are indicated in bold.

###  5.3 Hubbard-Kanamori model

Now, we focus on investigating the circuit depths for the Hubbard-Kanamori
model on the square lattice and as a function of the number of orbitals $M$
per lattice site. Compared to the FH model, the HK model has additional
density-density terms between different orbitals ($U_{1}$ and $U_{2}$) and the
Hund’s coupling terms related spin-exchange and pair-hopping terms ($J$) that
need to be implemented.

We pursue a similar strategy to the FH model and group all physical modes of a
given lattice site into horizontal blocks on which we perform an even-odd
fermionic swap network [42]. The modes on a site are ordered, as described
earlier, in a way that ensures that all $J$ terms can be implemented within
adjacent qubits and decomposed into depth-five circuits. This order is
inverted between horizontally subsequent blocks of modes of neighboring
lattice sites. We only consider the PAA mapping for the HK model, as it has
unambiguously demonstrated the best performance already at the level of the FH
model. We also do not consider the model on any other lattice geometries as we
do not expect this to yield additional insights.

Figure 5: a) Fermionic connectivity graph for a 3-by-3 square NN lattice with
4 orbitals per lattice site with two spin-modes each. b) Vertex (black) and
edge (blue, green, red, purple, pink) operators mapping this fermionic
connectivity to a square qubit layout using the PAA mapping. Each lattice site
is encoded into a linear chain consisting of four spin-up (full) and four
spin-down (empty) physical qubits (squares) connected to the neighboring sites
through ancilla qubits (circles). Dashed qubit connections are not required.
c) Minimal-depth circuit for the simulation of one Trotter step of the
Hubbard-Kanamori model in 55 parallel steps. Size of blocks is not to scale:
$t_{h}$-$t_{v}$-$t_{v}$ blocks are depth-seven, $t_{h}$-terms are depth-three,
J terms are depth-five, $U$, $U_{1}$, $U_{2}$ terms and fSWAPs can be
performed with a single two-qubit gate. Further circuit compression can be
applied within grey shaded areas of the circuit.

In Fig. 5 we study the example of the four-band HK model on the square lattice
using PAA in more detail (see Fig. 5.a)). In Fig. 5.b) we show a subset of the
vertex and edge operators for this mapping. Note that all edges between
adjacent modes of the same lattice site are given by either XX or YY (pink).
We also differentiate between vertices for the two different spin types (empty
vs full squares). In Fig. 5.c) we provide the shallowest, 55-depth, circuit we
were able to find for this model. It can be seen that, with the exception of
some fSWAP and $J$ layers, the depth is fully governed by hopping operations.
Similarly to the FH model we split the hoppings into circuit blocks with depth
7 ($t_{h}$-$t_{v}$-$t_{v}$) and depth 3 ($t_{h}$). All of the density-density
terms, those layers of the fSWAP network and those Hund’s terms which don’t
involve any modes on the external edges of the chain can be effectively
parallelized with the hoppings terms. In terms of the swap network, we can
omit one fSWAP layer at the start and at the end of the circuit whilst still
ensuring that every mode has travelled to both edges. In some cases, we delay
fSWAP operations within a layer in order to be able to efficiently implement
density-density terms. One can achieve some additional reduction in circuit
depth by combining fSWAP operations with Hund’s or density-density terms.

This approach can be generally applied also for the HK model with any other
number of orbitals $M$ with only small modifications. In Appendix D we show
the circuits for $M=2$ and $M=3$, with circuit depths of 27 and 35,
respectively. For $M\geq 4$ we can provide a general expression for the depth
per single Trotter layer, $28\lceil\frac{M}{2}\rceil-1$. For an even number of
orbitals this is obtained by counting $M$ layers of hopping terms with depth
of $3+7$, $2M-1$ non-parallelizable layers of fSWAPs of which $\frac{M}{2}$
can be combined with the $J$ terms and $\frac{M}{2}$ non-parallelizable layers
of $J$ terms with depth 5. For an odd number of orbitals (with the exception
of $M=3$) we did not find any fSWAP network that guarantees that all Hund’s
terms can be implemented. Our most efficient strategy is to fill up the chain
of each lattice site with an additional _dummy_ orbital, corresponding to a
pair of qubits, and to proceed as in the case of an even number of orbitals
(hence the ceiling function in the expression). Similarly to the HK model, we
can provide an expression for the depth of a multi-orbital FH model, which is
equivalent to removing all $U_{1}$, $U_{2}$ and $J$ terms from the circuit. In
this case we obtain a depth of $12M-1$ for both even and odd numbers of
orbitals $M$.

Some extra terms in a hypothetical Hamiltonian could, in principle, be
simulated without any additional resource costs. For example, this includes
intra-orbital hoppings, as they can be always be compressed together with
fSWAP operations. Another class of terms that may be relatively cheap to add
are three-mode terms acting on two orbitals.

##  6 Comparison to standard decomposition techniques

In this section we attempt to disentangle the improvement over state-of-the-
art due to our choice of mappings from that coming from the XYZ decomposition
and compression techniques. To this end we investigate the circuit depths for
four examples: the TB model on an NNN square lattice using the PA mapping, the
FH model on an NNN square lattice using PAA and DK mappings, and the HK model
on the square lattice using PAA. For these models we attempt to maximaly
reduce the circuit depth both using XYZ and standard decomposition techniques
(see X-shaped circuit in Eq. 1) and consider either CNOT or the parameterized
fSIM gates to be native. Generally, any two-qubit gate can be transformed into
at most three CNOTs plus single-qubit gates. However, in the case of $e^{i\phi
Z_{1}Z_{2}}$ and $e^{i\phi(X_{1}Y_{2}+Y_{1}X_{2})}$, we only need one a two
CNOTs, respectively:

$\raisebox{-0.4pt}{\includegraphics[scale={0.085}]{zztocnota.png}}=\raisebox{-0.4pt}{\includegraphics[scale={0.085}]{zztocnotb.png}},\
\
\raisebox{-0.4pt}{\includegraphics[scale={0.085}]{xytocnota.png}}=\raisebox{-0.4pt}{\includegraphics[scale={0.085}]{xytocnotb.png}}$
(35)

where the central Y-rotation has the opposite angle of the XY+YX rotation.
When comparing the XYZ formalism to the standard techniques, it is clear that
compression identities also exist for sequences of two-qubit CNOT and single-
qubit gates. We give an overview over compression identities we have
considered for circuits from the standard decomposition techniques in Appendix
A. Equally, we also take into account that it is possible to combine two CNOT
gates with a parameterized single-qubit gate between them (as shown in Eq. 1)
into a single fSIM gate.

In Table 3 we present the summary of our findings for circuit depth of various
configurations. As expected, we find that the combination of XYZ and fSIM
gates yields the lowest circuit depths throughout. Further, both the choice of
native gate set and the decomposition technique have a strong influence on the
result. If we compare the XYZ decomposition with the standard one counted in
fSIM gates, we see an improvement of up to a factor $2.3$ in the case of the
FH model using PAA. Interestingly, even if we only allow for CNOT gates there
is an improvement of up to a factor $1.8$, for the same model and mapping.
This holds true for any of the examples in the Table and suggests that using
XYZ decomposition generally leads to lower-depth circuits. The number of gates
is also reduced to a new minimal bound with respect to previous state-of-the-
art when using XYZ-decomposition and DK mapping, since hoppings can be
implemented with half of the two-qubit gates needed before.

We can also attempt to quantify the total improvement we have achieved related
to the previous state-of-art, which we identify as the combination of using
the standard decomposition, parameterized gates and the DK mapping. For the FH
model on the NNN square lattice this combination yields a circuit depth of 97
per single Trotter layer, compared to 30 for using XYZ together with the PAA
mapping. This corresponds to an improvement of 70 %.

Decomp. |  Native TQGs |  TB NNN PA |  FH NNN PAA |  FH NNN DK |  HK NN PAA
---|---|---|---|---|---
XYZ | fSIM | 17 | 30 | 74 | 55
Standard | fSIM | 31 | 68 | 97 | 84
XYZ | CNOT | 30 | 53 | 125 | 100
Standard | CNOT | 47 | 93 | 130 | 132
Table 3: Depths of a single Trotter step for different models using either the
XYZ or the standard decomposition and expressed in terms of different native
two-qubit gates.

##  7 Discussion

We have presented an efficient operator decomposition and circuit compression
formalism. Whilst generally leading to improvements in circuit depth for the
three fermionic systems studied here, it can, in principle, be applied to any
other systems of interest, including non-fermionic ones [49]. We found that
the XYZ formalism is more powerful when native fSIM gates are available,
however, it can still lead to improvements over standard decomposition
techniques when parameterized gates are not available since it can optimally
combine the fermionic operators acting on each pair of modes. We have also
introduced a set of local fermion-to-qubit mappings which are restricted to a
realistic, square, qubit connectivity graph. For systems with multiple modes
per lattice site we have chosen the strategy of using the composite edge rule
to implement inter-site operators. At the same time we have placed modes
belonging to the same site onto local Jordan-Wigner chains and performed a
fermionic swap network to allow for all inter-site operators to be
implemented. This idea is much in the spirit of Ref. [32].

Especially the PAA mapping allows for an increased degree of parallelism which
leads to significantly lower circuit depths for scientifically relevant
systems compared to previously introduced mappings [36, 37, 40]. In
particular, we found one single Trotter step of the FH model on a honeycomb
lattice, the FH model on an NNN square lattice and the three-band HK model on
a NN square lattice can be simulated on a square QPU in depths 9, 30 and 35,
respectively. These, to the best of our knowledge, correspond to either the
first or the lowest number reported for these examples to date. It is
important to state that the chosen FH and HK models constitute unsolved
actively researched strongly correlated systems in condensed matter physics
with actively debated phase diagrams [76] and direct relevance to the physics
of cuprates, nickelates and transitions-metal oxides. Our formalism could also
be applied to ab-initio models from quantum chemistry, especially when the
degree of geometric locality allows to neglect transfer integrals below a
given threshold. Since the here developed mappings allow for low-weight
interactions between far-neighbors in one lattice dimension (see Table 4),
they might be particularly useful for materials with weak couplings in the
second and third dimensions [77, 78]. Concluding, we believe that the circuit
depths reported in this work certainly improve the chances of fermionic
simulation algorithms successfully being implemented on future NISQ quantum
hardware.

We note that the parallelism of the PAA mapping comes at the expense of
introducing additional ancilla qubits which don’t necessarily have to be
readily available in current processors. The use of techniques meant to reduce
the number of qubits by using the symmetries of the Hamiltonian [79] could
potentially improve this issue. Equally, if the goal is to optimize for the
minimal total gate count, rather than the circuit depth, we expect the DK
mapping to yield favourable results compared to the other mappings introduced
here.

We would like to stress that in most cases one single Trotter layer is not
sufficient to reliably simulate a system and one would potentially need
(polynomially) many such layers instead. This is due to the Trotter error
introduced from approximating the limit in Eq. 33 with a finite number of
steps. It would be of interest to study the actual necessary number of Trotter
steps for these fermionic model in various regimes, for various system sizes,
as well as the impact on the Trotter error that the ordering of operators in
the decomposition has. It has been shown that this effect can be significant
[80] and could therefore easily out-weight the minimal additional
cancellations achieved by our optimal orderings of operators. On the other
hand, it was also shown in Ref. [80] that most orderings produce close to
optimal Trotter errors, whereas finding the optimal ordering is a factorially
hard problem.

We have not discussed efficient state preparation and measurement strategies
in this work and point the interested reader to Refs [81, 82] for more details
on these subjects. We have also not considered the effect of errors on our
circuits, as has been done in Ref. [31]. We note that the results of such
analyses are highly error-model dependent and leave this to future studies.

## Acknowledgments

The authors would like to thank S. Cheylan, S. Ikäläinen, F. R. Fernandes
Pereira, I. de Vega and K. Q. Dog for useful discussions.

## References

  * McArdle  _et al._ [2020] S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Rev. Mod. Phys. 92, 015003 (2020).
  * Delgado  _et al._ [2022] A. Delgado, P. A. M. Casares, R. d. Reis, M. S. Zini, R. Campos, N. Cruz-Hernández, A.-C. Voigt, A. Lowe, S. Jahangiri, M. A. Martin-Delgado, J. E. Mueller, and J. M. Arrazola, How to simulate key properties of lithium-ion batteries with a fault-tolerant quantum computer (2022).
  * Andersson  _et al._ [2022] M. P. Andersson, M. N. Jones, K. V. Mikkelsen, F. You, and S. S. Mansouri, Current Opinion in Chemical Engineering 36, 100754 (2022).
  * Reiher  _et al._ [2017] M. Reiher, N. Wiebe, K. M. Svore, D. Wecker, and M. Troyer, Proceedings of the National Academy of Sciences 114, 7555 (2017).
  * LeBlanc  _et al._ [2015] J. P. F. LeBlanc, A. E. Antipov, F. Becca, I. W. Bulik, G. K.-L. Chan, C.-M. Chung, Y. Deng, M. Ferrero, T. M. Henderson, C. A. Jiménez-Hoyos, E. Kozik, X.-W. Liu, A. J. Millis, N. V. Prokof’ev, M. Qin, G. E. Scuseria, H. Shi, B. V. Svistunov, L. F. Tocchio, I. S. Tupitsyn, S. R. White, S. Zhang, B.-X. Zheng, Z. Zhu, and E. Gull (Simons Collaboration on the Many-Electron Problem), Phys. Rev. X 5, 041041 (2015).
  * Motta  _et al._ [2017] M. Motta, D. M. Ceperley, G. K.-L. Chan, J. A. Gomez, E. Gull, S. Guo, C. A. Jiménez-Hoyos, T. N. Lan, J. Li, F. Ma, A. J. Millis, N. V. Prokof’ev, U. Ray, G. E. Scuseria, S. Sorella, E. M. Stoudenmire, Q. Sun, I. S. Tupitsyn, S. R. White, D. Zgid, and S. Zhang (Simons Collaboration on the Many-Electron Problem), Phys. Rev. X 7, 031059 (2017).
  * Schäfer  _et al._ [2021] T. Schäfer, N. Wentzell, F. Šimkovic, Y.-Y. He, C. Hille, M. Klett, C. J. Eckhardt, B. Arzhang, V. Harkov, F. m. c.-M. Le Régent, A. Kirsch, Y. Wang, A. J. Kim, E. Kozik, E. A. Stepanov, A. Kauch, S. Andergassen, P. Hansmann, D. Rohe, Y. M. Vilk, J. P. F. LeBlanc, S. Zhang, A.-M. S. Tremblay, M. Ferrero, O. Parcollet, and A. Georges, Phys. Rev. X 11, 011058 (2021).
  * Georgescu  _et al._ [2014] I. Georgescu, S. Ashhab, and F. Nori, Reviews of Modern Physics 86, 153 (2014).
  * Troyer and Wiese [2005] M. Troyer and U.-J. Wiese, Phys. Rev. Lett. 94, 170201 (2005).
  * Rossi  _et al._ [2017] R. Rossi, N. Prokof’ev, B. Svistunov, K. Van Houcke, and F. Werner, Europhysics Letters 118, 10004 (2017).
  * Lee  _et al._ [2022] S. Lee, J. Lee, H. Zhai, Y. Tong, A. M. Dalzell, A. Kumar, P. Helms, J. Gray, Z.-H. Cui, W. Liu, _et al._ , arXiv preprint arXiv:2208.02199 (2022).
  * Koch  _et al._ [2020] D. Koch, B. Martin, S. Patel, L. Wessing, and P. M. Alsing, AIP Advances 10, 095101 (2020).
  * Arute  _et al._ [2020] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, A. Bengtsson, S. Boixo, M. Broughton, B. B. Buckley, D. A. Buell, B. Burkett, N. Bushnell, Y. Chen, Z. Chen, Y.-A. Chen, B. Chiaro, R. Collins, S. J. Cotton, W. Courtney, S. Demura, A. Derk, A. Dunsworth, D. Eppens, T. Eckl, C. Erickson, E. Farhi, A. Fowler, B. Foxen, C. Gidney, M. Giustina, R. Graff, J. A. Gross, S. Habegger, M. P. Harrigan, A. Ho, S. Hong, T. Huang, W. Huggins, L. B. Ioffe, S. V. Isakov, E. Jeffrey, Z. Jiang, C. Jones, D. Kafri, K. Kechedzhi, J. Kelly, S. Kim, P. V. Klimov, A. N. Korotkov, F. Kostritsa, D. Landhuis, P. Laptev, M. Lindmark, E. Lucero, M. Marthaler, O. Martin, J. M. Martinis, A. Marusczyk, S. McArdle, J. R. McClean, T. McCourt, M. McEwen, A. Megrant, C. Mejuto-Zaera, X. Mi, M. Mohseni, W. Mruczkiewicz, J. Mutus, O. Naaman, M. Neeley, C. Neill, H. Neven, M. Newman, M. Y. Niu, T. E. O’Brien, E. Ostby, B. Pató, A. Petukhov, H. Putterman, C. Quintana, J.-M. Reiner, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, D. Strain, K. J. Sung, P. Schmitteckert, M. Szalay, N. M. Tubman, A. Vainsencher, T. White, N. Vogt, Z. J. Yao, P. Yeh, A. Zalcman, and S. Zanker, Observation of separated dynamics of charge and spin in the Fermi-Hubbard model (2020).
  * Tazhigulov  _et al._ [2022] R. N. Tazhigulov, S.-N. Sun, R. Haghshenas, H. Zhai, A. T. K. Tan, N. C. Rubin, R. Babbush, A. J. Minnich, and G. K.-L. Chan, Simulating challenging correlated molecules and materials on the Sycamore quantum processor (2022).
  * Stanisic  _et al._ [2022] S. Stanisic, J. L. Bosse, F. M. Gambetta, R. A. Santos, W. Mruczkiewicz, T. E. O’Brien, E. Ostby, and A. Montanaro, Nature Communications 13, 10.1038/s41467-022-33335-4 (2022).
  * Jordan and Wigner [1928] P. Jordan and E. Wigner, Zeitschrift für Physik 47, 631 (1928).
  * Bravyi  _et al._ [2017] S. Bravyi, J. M. Gambetta, A. Mezzacapo, and K. Temme, Tapering off qubits to simulate fermionic Hamiltonians (2017).
  * Steudtner and Wehner [2018] M. Steudtner and S. Wehner, New Journal of Physics 20, 063010 (2018).
  * Chiew and Strelchuk [2021] M. Chiew and S. Strelchuk, Optimal fermion-qubit mappings (2021).
  * Chien and Klassen [2022] R. W. Chien and J. Klassen, Optimizing fermionic encodings for both Hamiltonian and hardware (2022).
  * Kirby  _et al._ [2022] W. Kirby, B. Fuller, C. Hadfield, and A. Mezzacapo, PRX Quantum 3, 10.1103/prxquantum.3.020351 (2022).
  * Harrison  _et al._ [2022] B. Harrison, D. Nelson, D. Adamiak, and J. Whitfield, Reducing the qubit requirement of Jordan-Wigner encodings of $n$-mode, $k$-fermion systems from $n$ to $\lceil\log_{2}{N\choose K}\rceil$ (2022).
  * Havlíček  _et al._ [2017] V. Havlíček, M. Troyer, and J. D. Whitfield, Physical Review A 95, 10.1103/physreva.95.032332 (2017).
  * Vlasov [2019] A. Y. Vlasov, Clifford algebras, spin groups and qubit trees (2019).
  * Miller  _et al._ [2022] A. Miller, Z. Zimborás, S. Knecht, S. Maniscalco, and G. García-Pérez, The Bonsai algorithm: grow your own fermion-to-qubit mapping (2022).
  * Haah [2017] J. Haah, Revista Colombiana de Matemáticas 50, 299 (2017).
  * Bausch  _et al._ [2020] J. Bausch, T. Cubitt, C. Derby, and J. Klassen, Mitigating errors in local fermionic encodings (2020).
  * Setia  _et al._ [2019] K. Setia, S. Bravyi, A. Mezzacapo, and J. D. Whitfield, Physical Review Research 1, 10.1103/physrevresearch.1.033033 (2019).
  * Tantivasadakarn [2020] N. Tantivasadakarn, Physical Review Research 2, 10.1103/physrevresearch.2.023353 (2020).
  * Chen  _et al._ [2022] Y.-A. Chen, A. V. Gorshkov, and Y. Xu, Error-correcting codes for fermionic quantum simulation (2022).
  * Clinton  _et al._ [2021] L. Clinton, J. Bausch, and T. Cubitt, Nature Communications 12, 10.1038/s41467-021-25196-0 (2021).
  * Clinton  _et al._ [2022] L. Clinton, T. Cubitt, B. Flynn, F. M. Gambetta, J. Klassen, A. Montanaro, S. Piddock, R. A. Santos, and E. Sheridan, Towards near-term quantum simulation of materials (2022).
  * Nys and Carleo [2022] J. Nys and G. Carleo, Quantum 6, 833 (2022).
  * Huggins  _et al._ [2022] W. J. Huggins, B. A. O’Gorman, N. C. Rubin, D. R. Reichman, R. Babbush, and J. Lee, Nature 603, 416 (2022).
  * Bravyi and Kitaev [2002] S. B. Bravyi and A. Y. Kitaev, Annals of Physics 298, 210 (2002).
  * Ball [2005] R. C. Ball, Physical Review Letters 95, 10.1103/physrevlett.95.176407 (2005).
  * Verstraete and Cirac [2005] F. Verstraete and J. I. Cirac, Journal of Statistical Mechanics: Theory and Experiment 2005, P09012 (2005).
  * Chen  _et al._ [2018] Y.-A. Chen, A. Kapustin, and D. Radičević, Annals of Physics 393, 234 (2018).
  * Setia and Whitfield [2018] K. Setia and J. D. Whitfield, The Journal of Chemical Physics 148, 164104 (2018).
  * Derby  _et al._ [2021] C. Derby, J. Klassen, J. Bausch, and T. Cubitt, Physical Review B 104, 10.1103/physrevb.104.035118 (2021).
  * Chen and Xu [2022] Y.-A. Chen and Y. Xu, Equivalence between fermion-to-qubit mappings in two spatial dimensions (2022).
  * Kivlichan  _et al._ [2018] I. D. Kivlichan, J. McClean, N. Wiebe, C. Gidney, A. Aspuru-Guzik, G. K.-L. Chan, and R. Babbush, Physical Review Letters 120, 10.1103/physrevlett.120.110501 (2018).
  * O’Gorman  _et al._ [2019] B. O’Gorman, W. J. Huggins, E. G. Rieffel, and K. B. Whaley, Generalized swap networks for near-term quantum computing (2019).
  * Cade  _et al._ [2020] C. Cade, L. Mineh, A. Montanaro, and S. Stanisic, Physical Review B 102, 10.1103/physrevb.102.235122 (2020).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, B. Burkett, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, A. Dunsworth, E. Farhi, B. Foxen, A. Fowler, C. Gidney, M. Giustina, R. Graff, K. Guerin, S. Habegger, M. P. Harrigan, M. J. Hartmann, A. Ho, M. Hoffmann, T. Huang, T. S. Humble, S. V. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, J. Kelly, P. V. Klimov, S. Knysh, A. Korotkov, F. Kostritsa, D. Landhuis, M. Lindmark, E. Lucero, D. Lyakh, S. Mandrà, J. R. McClean, M. McEwen, A. Megrant, X. Mi, K. Michielsen, M. Mohseni, J. Mutus, O. Naaman, M. Neeley, C. Neill, M. Y. Niu, E. Ostby, A. Petukhov, J. C. Platt, C. Quintana, E. G. Rieffel, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, K. J. Sung, M. D. Trevithick, A. Vainsencher, B. Villalonga, T. White, Z. J. Yao, P. Yeh, A. Zalcman, H. Neven, and J. M. Martinis, Nature 574, 505 (2019).
  * Valery  _et al._ [2022] J. A. Valery, S. Chowdhury, G. Jones, and N. Didier, PRX Quantum 3, 10.1103/prxquantum.3.020337 (2022).
  * Nielsen and Chuang [2000] M. Nielsen and I. Chuang,  _Quantum computation and quantum information_ (Cambridge University Press, 2000).
  * Cowtan  _et al._ [2020] A. Cowtan, S. Dilkes, R. Duncan, W. Simmons, and S. Sivarajah, Electronic Proceedings in Theoretical Computer Science 318, 213 (2020).
  * Sriluckshmy  _et al._ [2023] P. V. Sriluckshmy, V. Pina-Canelles, M. Ponce, M. G. Algaba, F. Šimkovic, and M. Leib, Optimal, hardware native decomposition of parameterized multi-qubit pauli gates (2023).
  * Peng  _et al._ [2022] B. Peng, S. Gulania, Y. Alexeev, and N. Govind, Physical Review A 106, 10.1103/physreva.106.012412 (2022).
  * Foxen  _et al._ [2020] B. Foxen, C. Neill, A. Dunsworth, P. Roushan, B. Chiaro, A. Megrant, J. Kelly, Z. Chen, K. Satzinger, R. Barends, _et al._ , Physical Review Letters 125, 120504 (2020).
  * Kanamori [1963] J. Kanamori, Progress of Theoretical Physics 30, 275 (1963).
  * Georges  _et al._ [2013] A. Georges, L. de' Medici, and J. Mravlje, Annual Review of Condensed Matter Physics 4, 137 (2013).
  * Hao  _et al._ [2019] H. Hao, B. M. Rubenstein, and H. Shi, Physical Review B 99, 10.1103/physrevb.99.235142 (2019).
  * Szasz  _et al._ [2020] A. Szasz, J. Motruk, M. P. Zaletel, and J. E. Moore, Phys. Rev. X 10, 021042 (2020).
  * Wietek  _et al._ [2021] A. Wietek, R. Rossi, F. Šimkovic, M. Klett, P. Hansmann, M. Ferrero, E. M. Stoudenmire, T. Schäfer, and A. Georges, Phys. Rev. X 11, 041013 (2021).
  * Raczkowski  _et al._ [2020] M. Raczkowski, R. Peters, T. T. Phùng, N. Takemori, F. F. Assaad, A. Honecker, and J. Vahedi, Phys. Rev. B 101, 125103 (2020).
  * Yang  _et al._ [2021] X. Yang, H. Zheng, and M. Qin, Phys. Rev. B 103, 155110 (2021).
  * Medeiros-Silva  _et al._ [2023] A. R. Medeiros-Silva, N. C. Costa, and T. Paiva, Phys. Rev. B 107, 035134 (2023).
  * Jiang  _et al._ [2023] S. Jiang, D. J. Scalapino, and S. R. White, Density-matrix-renormalization-group-based downfolding of the three-band hubbard model: the importance of density-assisted hopping (2023).
  * Algaba  _et al._ [2022] M. G. Algaba, M. Ponce-Martinez, C. Munuera-Javaloy, V. Pina-Canelles, M. J. Thapa, B. G. Taketani, M. Leib, I. de Vega, J. Casanova, and H. Heimonen, Phys. Rev. Res. 4, 043089 (2022).
  * Lloyd [1996] S. Lloyd, Science 273, 1073 (1996).
  * Suzuki [1976] M. Suzuki, Communications in Mathematical Physics 51, 183 (1976).
  * Childs and Wiebe [2012] A. M. Childs and N. Wiebe, arXiv 10.48550/ARXIV.1202.5822 (2012).
  * Berry  _et al._ [2014] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, arXiv 10.48550/ARXIV.1412.4687 (2014).
  * Childs  _et al._ [2019] A. M. Childs, A. Ostrander, and Y. Su, Quantum 3, 182 (2019).
  * Campbell [2019] E. Campbell, Physical Review Letters 123, 10.1103/physrevlett.123.070503 (2019).
  * Kitaev [1995] A. Y. Kitaev, Quantum measurements and the Abelian stabilizer problem (1995).
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, Nature Communications 5, 10.1038/ncomms5213 (2014).
  * Wecker  _et al._ [2015] D. Wecker, M. B. Hastings, and M. Troyer, Physical Review A 92, 10.1103/physreva.92.042303 (2015).
  * Yuan  _et al._ [2019] X. Yuan, S. Endo, Q. Zhao, Y. Li, and S. Benjamin, Quantum 3, 191 (2019).
  * Grimsley  _et al._ [2019] H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, Nature Communications 10, 10.1038/s41467-019-10988-2 (2019).
  * Somma [2019] R. D. Somma, New Journal of Physics 21, 123025 (2019).
  * Bharti and Haug [2021a] K. Bharti and T. Haug, Physical Review A 104, 10.1103/physreva.104.l050401 (2021a).
  * Bharti and Haug [2021b] K. Bharti and T. Haug, Physical Review A 104, 042418 (2021b).
  * Huang  _et al._ [2022] X.-X. Huang, B. Moritz, M. Claassen, and T. P. Devereaux, Physical Review B 105, 10.1103/physrevb.105.165124 (2022).
  * Helmholz [1947] L. Helmholz, Journal of the American Chemical Society 69, 886 (1947).
  * Amsler  _et al._ [2023] M. Amsler, P. Deglmann, M. Degroote, M. P. Kaicher, M. Kiser, M. Kühn, C. Kumar, A. Maier, G. Samsonidze, A. Schroeder, M. Streif, D. Vodola, and C. Wever, Quantum-enhanced quantum monte carlo: an industrial view (2023).
  * Ralli  _et al._ [2022] A. Ralli, T. Weaving, A. Tranter, W. M. Kirby, P. J. Love, and P. V. Coveney, Unitary partitioning and the contextual subspace variational quantum eigensolver (2022).
  * Tranter  _et al._ [2019] A. Tranter, P. J. Love, F. Mintert, N. Wiebe, and P. V. Coveney, Entropy 21, 1218 (2019).
  * Bharti  _et al._ [2021] K. Bharti, A. Cervera-Lierta, T. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. Kottmann, T. Menke, _et al._ , arXiv preprint arXiv:2101.08448 (2021).
  * Tilly  _et al._ [2022] J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth, and J. Tennyson, Physics Reports 986, 1 (2022).
  * Abraham  _et al._ [2021] H. Abraham, AduOffei, R. Agarwal, I. Y. Akhalwaya, and G. A. et al., Qiskit: An open-source framework for quantum computing (2021).

##  Appendix A CNOT identities

In order to obtain the circuit depths in Table 3, some CNOT identities were
used with the aim of compressing the circuit and making a fair comparison of
the XYZ formalism against the standard decomposition techniques. Most of these
identities have been applied on the boundaries between different fermionic
operators. The identities, obtained from Qiskit [83], use CNOTs, H as Hadamard
gate, S as phase gate and any Pauli rotation
$R_{\sigma}(\alpha)=e^{-i\frac{\alpha}{2}\sigma}$ where
$\sigma\in\\{X,Y,Z\\}$, as well as their Hermitian conjugate:

$\begin{matrix}\raisebox{-0.4pt}{\includegraphics[scale={0.400}]{cnotcomp1a.pdf}}=\raisebox{-0.45pt}{\includegraphics[scale={0.400}]{cnotcomp1b.pdf}},\\\
\\\
\raisebox{-0.4pt}{\includegraphics[scale={0.400}]{cnotcomp2a.pdf}}=\raisebox{-0.43pt}{\includegraphics[scale={0.400}]{cnotcomp2b.pdf}}\end{matrix}$
(36)
$\begin{matrix}\raisebox{-0.4pt}{\includegraphics[scale={0.400}]{cnotcomp3a.pdf}}=\raisebox{-0.47pt}{\includegraphics[scale={0.400}]{cnotcomp3b.pdf}},\\\
\\\
\raisebox{-0.4pt}{\includegraphics[scale={0.400}]{cnotcomp4a.pdf}}=\raisebox{-0.45pt}{\includegraphics[scale={0.400}]{cnotcomp4b.pdf}}\end{matrix}$
(37)
$\begin{matrix}\raisebox{-0.4pt}{\includegraphics[scale={0.400}]{cnotcomp5a.pdf}}=\raisebox{-0.4pt}{\includegraphics[scale={0.400}]{cnotcomp5b.pdf}}\end{matrix}$
(38)
$\begin{matrix}\raisebox{-0.4pt}{\includegraphics[scale={0.400}]{cnotcomp6a.pdf}}=\raisebox{-0.4pt}{\includegraphics[scale={0.400}]{cnotcomp6b.pdf}},\\\
\\\
\raisebox{-0.4pt}{\includegraphics[scale={0.400}]{cnotcomp7a.pdf}}=\raisebox{-0.4pt}{\includegraphics[scale={0.400}]{cnotcomp7b.pdf}}\end{matrix}$
(39)

We note that, in principle, CNOT (and XYZ) decompositions exist that yield a
logarithmic depth at the cost of having an all-to-all connectivity structure
[32]. However, the would not improve the results obtained in this paper, since
for weights lower than 6 using a binary tree does not make ant difference
compared to using the standard approach. Hence, higher connectivities will not
give extra any advantage for the cases considered using our mappings.

##  Appendix B Embedding the DK mapping into a square qubit layout

For the sake of completeness, we give a brief overview on the DK mapping
proposed in Ref. [40]. It is a local fermion-to-qubit mapping constructed
through the definition of edge and vertex operators, much like the mappings
presented in the main text. It has a fermion-to-qubit ratio of $r=1.5$ and
weight operators $w=3$.

The original topology can be understood as a square grid of physical qubits
with ancilla qubits placed on alternating faces of the lattice, connected to
the four neighboring qubits. See the dashed lattice from Fig. 6.a), in which
some faces have no ancilla qubits in the middle. Each edge $(i,j)$ on a face
with no ancilla qubits is assigned an orientation following a clockwise or
counterclockwise direction over the face edges, alternating on every row of
faces. The vertex and edge operators are:

$V_{i}\equiv Z_{i}$ (40)
$E_{i,j}\equiv\left\\{\begin{array}[]{cl}X_{i}Y_{j}X_{f(i,j)}&\mathrm{if}\
(i,j)\ \textrm{is oriented downwards}\\\ -X_{i}Y_{j}X_{f(i,j)}&\mathrm{if}\
(i,j)\ \textrm{is oriented upwards}\\\ X_{i}Y_{j}Y_{f(i,j)}&\mathrm{if}\
(i,j)\ \textrm{is horizontal}\end{array}\right.$ (41)

where $i,j$ are physical qubits and $f$ is an ancilla qubit on the face
corresponding to (i,j). This mapping yields weight-3 hopping operators,
weight-2 density-density operators and weight-8 stabilizers.

(a)

(b)

Figure 6: Edge operators for a) the original DK mapping [40] with one spin and
for b) the square grid-embedded version. Squares represent physical qubits and
circles are ancilla qubits. The dashed lines show the edge operators
connectivity while the solid lines are qubit connections.

(a)

(b)

Figure 7: a) Spinful DK mapping from [32] and b) its embedding on a square grid topology. Full (empty) squares represent physical qubits with spin up (down) qubits and circles are ancilla qubits. The dashed lines show the edge operators connectivity while the solid lines are qubit connections. e  Axis | Weight (DK) | Weight (PA) | Weight (PAA)
---|---|---|---
$\mathrel{\text{\rotatebox[origin={c}]{0.0}{$\vrule height=0.0pt,width=0.0pt\shortrightarrow$}}}$ | $2+s+M(s-1)$ | $2+s+M(s-1)$ | $2+2s+M(s-1)$
$\mathrel{\text{\rotatebox[origin={c}]{90.0}{$\vrule height=0.0pt,width=0.0pt\shortrightarrow$}}}$ | $2+s+M(s-1)$ | $3+s$ | $3+s$
$\mathrel{\text{\rotatebox[origin={c}]{45.0}{$\vrule height=0.0pt,width=0.0pt\shortrightarrow$}}}$ | $2+2s+Ms$ | $2+2s+M(s-1)$ | $2+3s+M(s-1)$
$\mathrel{\text{\rotatebox[origin={c}]{135.0}{$\vrule height=0.0pt,width=0.0pt\shortrightarrow$}}}$ | $2+2s+Ms$ | $2+2s+M(s-1)$ | $2+3s+M(s-1)$
Table 4: Worst-case operator weight of the edge operators between lattice
sites separated by $s\geq 1$ steps on the respective axis for a fermionic
model with $M$ modes.

While it is obvious that the DK mapping can be embedded into a 4-connected
square grid topology by adding idling qubits in the faces with no ancilla
qubits, an embedding without idling qubits can also be found as shown in Fig.
6.b). We found that the circuit depths obtained in Table 1 do not change if we
use this implementation of the DK mapping.

An extension to multi-orbital models was proposed in [32], substituting every
physical qubit from the DK mapping by a chain of physical qubits representing
different fermionic modes on each site and using JW chains for defining the
edges, identical to our treatment in the main text. The resulting topology for
two spins can be seen in Fig. 7.a) and its (more straight-forward) embedding
into a square grid topology is shown in Fig. 7.b).

The interactions between fermions on non-neighboring sites can also be
included in the DK mapping, however the operator weights will scale
differently with the distance and number of modes per lattice site. In Table 4
we show the edge operator weights of the DK, PA and PAA mappings for a model
with $M$ modes per lattice site between far-neighbors along different axes,
and separated by $s\geq 1$ steps. For the DK mapping, the unit cell consists
of two lattice sites for which this edge operator weight in some cases scales
differently. We consider the worst-case operator weight in such case. While in
DK, the edges along all axes scale with $M$, in the case of PA and PAA lattice
sites can be connected vertically through ancilla qubits, and thus do not
scale with $M$. For some multi-orbital fermionic systems with higher-connected
geometries this can scaling can yield an advantage in the operator weights and
hence also the circuit depths. Since the linear dependence on $M$ is usually
translated into quadratic dependence of the circuit depth, we can therefore
expect to reduce the depth from quadratic to linear in some cases.

Figure 8:  a) Qubit layout representing the JW snake-like pattern of a TB
model on a 3-by-3 fermionic lattice and b) quantum circuit of the
corresponding Trotter step, showing in grey the combinations of gates and
through red lines the cancellations of gates. The two-qubit gates represent
fSIM gates.  (a)  Figure 9: Single-Trotter step simulation circuits for
Hubbard-Kanamori model with a) 2 orbitals and b) 3 orbitals. Size of blocks is
not to scale: $t_{h}$-$t_{v}$-$t_{v}$ blocks are depth-seven, $t_{h}$-terms
are depth-three, J terms are depth-five, $U$, $U_{1}$, $U_{2}$ terms and
fSWAPs can be performed with a single two-qubit gate. Further circuit
compression can be applied within grey shaded areas of the circuit.

##  Appendix C Jordan-Wigner vs PA on a 3-by-3 lattice

Due to the relatively local nature of hopping operators when using the Jordan
Wigner mapping (JW) on a small lattices, one could argue that this will yield
smaller circuit depths than local mappings. To address this question , we
investigate in Fig. 8 the minimal JW example of a TB model on a 3-by-3 square
grid and obtain a circuit depth of 10, which is worse than the depth 9
obtained for PA or depth 8 for PPAA in Table 1. Using the same approach, one
can also show that the number of two-qubit gates for this system size is
exactly the same between the Derby-Klassen mapping and JW (here the JW will ).
This result leads us to the conclusion that it is never favourable to use JW
for the simulation of the such fermionic models on two-dimensional lattices
with a system size larger than 3-by-3.

It should, however, be noted, that we have not included any overhead due to
state preparation for local mappings into this calculation. This overhead can
be either negligible, when stabilizers can be measured at the start of the
circuit, or otherwise add a one-time overhead at the beginning of the circuit,
which scales linearly with the system size. In the second case, it is expected
that local mappings still become favourable to the JW mapping, albeit at
somewhat larger system sizes.

##  Appendix D Hubbard-Kanamori simulation

The simulation of a HK model for 4 orbitals has been shown in the main text.
In this appendix, we additionally provide the circuits for the HK model with
two and three orbitals in Fig. 5.

In the case of two orbitals (Fig. 5.a) we follow the standard fSWAP network
strategy described in the main text. In order to minimize the circuit depth we
combine fSWAP operators with interaction terms from the Hamiltonian. Note that
in one case an fSWAP operator can be sandwiched between a $U_{1}$ and a $J$
interaction term. The final depth for this combination of operators is
therefore still five. As a consequence, we can implement the circuit in total
depth of 27.

In the case of the three-orbital model we pursue a somewhat different ordering
of modes within a site:
$\\{1\\!\\!\uparrow$,$2\\!\\!\uparrow$,$3\\!\\!\downarrow$,$2\\!\\!\downarrow$,$3\\!\\!\uparrow$,$1\\!\\!\downarrow\\}$
together with a modified fSWAP network (see the circuit in Fig. 5.b)), which
allows for fSWAP operators to be merged with $J$-terms from both the left and
the right sides. The total circuit depth for this model is 35, remarkably only
eight longer than for the two-orbital model. This is mainly due to the fact
that we can implement all quartic $J$-terms in parallel with the hoppings
between lattice sites. This system is also of scientific interest as it is
considered to be the minimal number of orbitals that is necessary for Hund’s
physics appears in the HK model [53].