# Expanding variational quantum eigensolvers to larger systems by dividing the
calculations between classical and quantum hardware

John P. T. Stenger  NRC Postdoctoral Associate, U.S. Naval Research
Laboratory, Washington, DC 20375, United States  Daniel Gunlycke  U.S. Naval
Research Laboratory, Washington, DC 20375, United States  C. Stephen Hellberg
U.S. Naval Research Laboratory, Washington, DC 20375, United States

###### Abstract

We present a hybrid classical/quantum algorithm for efficiently solving the
eigenvalue problem of many-particle Hamiltonians on quantum computers with
limited resources by splitting the workload between classical and quantum
processors. This algorithm reduces the needed number of qubits at the expense
of an increased number of quantum evaluations. We demonstrate the method for
the Hubbard model and show how the conservation of the $z$-component of the
total spin allows the spin-up and spin-down configurations to be computed on
classical and quantum hardware, respectively. Other symmetries can be
exploited in a similar manner.

##  I Introduction

When quantum computers meet their full potential, they will be able to
drastically out perform classical computers for a certain class of problems
[1, 2, 3, 4, 5, 6]. One example is determining the ground state of many-
particle interacting quantum systems [7, 8, 9, 10, 11, 12, 13, 14, 15, 16].
With the current state of Noisy Intermediate Scale Quantum (NISQ) devices, it
is impractical to have only algorithms that put the burden of computation
either totally on the quantum computer or totally on the classical computer.
Algorithms that split the workload between the classical and quantum devices
have the potential to exceed the performance of either device alone. In this
way, algorithms can evolve alongside the hardware, shifting more of the work
onto the quantum machines as their performance become more robust.

The variational quantum eigensolver (VQE) is anticipated to be one the first
algorithms that will show quantum superiority [17, 18, 19]. In the VQE
algorithm, the quantum state is prepared by applying a set of parameterized
quantum gates. Expectation values are calculated on the quantum device using
the quantum state and then the parameters are optimized on a classical
computer using the output from the quantum computer. VQE is already a nice
example of an algorithm that splits the work between the classical and quantum
devices. However, the limitation of current quantum hardware allows for only
very short quantum circuits. Recent works have proposed shifting more of the
algorithm to the classical computer by approximating parts of the simulated
system using tensor networks [20, 21], and by approximating the whole system
using classical shadows [22, 23]. We propose algorithms which reduce the
number of qubits by shifting the calculation of part of the simulated system
onto the classical computer without making approximations.

In this work, we propose hybrid quantum-classical algorithms to compute the
ground state of many body, fermionic Hamiltonians. First, we present a
generalized approach for splitting a quantum many-particle system into two
sections and a VQE type algorithm which can deal with one section classically
and the other quantum mechanically. In this generalized approach, there are no
restrictions on how the system is split. If the available classical resources
have been maxed out, then each additional one-particle state would have to be
included in the set treated by quantum hardware. Our algorithm scales in this
case polynomially with an increasing number of one-particle states. On the
other hand, if the available quantum resources have been maxed out, then each
additional one-particle state would have to be included in the set treated
exclusively by classical hardware. Our algorithm would then scale
exponentially with an increasing number of one-particle states. Next, we focus
on the Hubbard model and split the system into spin-up and spin-down sectors.
The spin up sector is solved on the classical computer while the spin-down
sector is solved on the quantum computer. We will describe three different
approaches to solving this spin-split model: In the first approach, the number
of spin-up electrons is fixed while the spin-down space includes all possible
electron numbers; in the second, we fix the parity of the spin-down sector; in
the third, the total number of electrons in both spin sectors are fixed.
Solving the spin-up sector on the classical computer cuts the required number
of qubits in half. Putting restrictions on the spin-down sector can further
reduce the number of qubits [24] but can also increase the classical workload.

##  II General Method

The interactions in a physical system of identical fermions can generally be
described by the Hamiltonian

$\hat{H}=\sum_{\mu\mu^{\prime}}t_{\mu\mu^{\prime}}c^{\dagger}_{\mu}c_{\mu^{\prime}}+\sum_{\mu\mu^{\prime}\nu\nu^{\prime}}v_{\mu\nu\mu^{\prime}\nu^{\prime}}c^{\dagger}_{\mu}c^{\dagger}_{\nu}c_{\nu^{\prime}}c_{\mu^{\prime}},$
(1)

where $c^{\dagger}_{\mu}$ and $c_{\mu}$ are fermionic creation and
annihilation operators acting on the one-particle index $\mu$ in the set
$\boldsymbol{\mathrm{M}}$ and $t_{\mu\mu^{\prime}}$ and
$v_{\mu\nu\mu^{\prime}\nu^{\prime}}$ are coefficients specifying the one- and
two-particle interactions, respectively. See the appendix for expressions that
can be used to compute the coefficients.

###  II.1 System Splitting

Our goal is to separate the computational problem into parts such that one
part is performed classically and the other part is performed using a modified
VQE approach. To that end, let us begin by defining two subsystems $A$ and $B$
and splitting the one-particle index set $\boldsymbol{\mathrm{M}}$ into two
sets $\boldsymbol{\mathrm{M}}_{A}$ and $\boldsymbol{\mathrm{M}}_{B}$. Next, we
define for each $\Upsilon\in\\{A,B\\}$, the subsystem creation operator

$C^{\dagger}_{\vec{n}_{\Upsilon}}=\prod_{\mu\in\boldsymbol{\mathrm{M}}_{\Upsilon}}\big{(}c^{\dagger}_{\mu}\big{)}^{n_{\mu}},$
(2)

where
$\vec{n}_{\Upsilon}=(n_{\mu})_{\mu\in\boldsymbol{\mathrm{M}}_{\Upsilon}}$ are
particle configurations comprising particle occupation numbers
$n_{\mu}\in\\{0,1\\}$. With these operators, we can introduce the Fock states

$\ket{\vec{n}}=C^{\dagger}_{\vec{n}_{A}}C^{\dagger}_{\vec{n}_{B}}\ket{0},$ (3)

where $\ket{0}$ is the vacuum state, which forms a basis for the antisymmetric
Fock space $\mathcal{F}$ for our system. Rather than mapping the entire
$\mathcal{F}$ to the Hilbert space of the quantum register, however, we only
plan to map the antisymmetric Fock space $\mathcal{F}_{B}$ for subsystem $B$.
A standard basis for the latter $\mathcal{F}_{B}$ can be formed by the Fock
states

$\ket{\vec{n}_{B}}=C^{\dagger}_{\vec{n}_{B}}\ket{0},$ (4)

for this subsystem. For any given particle configuration $\vec{n}_{A}$, we can
then expand an arbitrary state in $\mathcal{F}_{B}$ as

$\ket{\Psi_{\vec{n}_{A}}}=\sum_{\vec{n}_{B}}\beta_{\vec{n}_{A}\vec{n}_{B}}\ket{\vec{n}_{B}},$
(5)

where the coefficients $\beta_{\vec{n}_{A}\vec{n}_{B}}$ are normalized such
that $\braket{\Psi_{\vec{n}_{A}}}{\Psi_{\vec{n}_{A}}}=1$. An arbitrary state
in $\mathcal{F}$ can then be expressed as the linear combination

$\displaystyle\ket{\Psi}$
$\displaystyle=\sum_{\vec{n}_{A}}\alpha_{\vec{n}_{A}}C^{\dagger}_{\vec{n}_{A}}\ket{\Psi_{\vec{n}_{A}}}$
$\displaystyle=\sum_{\vec{n}_{A}\vec{n}_{B}}\alpha_{\vec{n}_{A}}\beta_{\vec{n}_{A}\vec{n}_{B}}C^{\dagger}_{\vec{n}_{A}}\ket{\vec{n}_{B}}$
(6)

with the normalization $\braket{\Psi}{\Psi}=1$. This form allows us to find
any particular state $\ket{\Psi}$, e.g., the ground state of the system, by
determining $\alpha_{\vec{n}_{A}}$ classically and
$\beta_{\vec{n}_{A}\vec{n}_{B}}$ using our modified VQE approach described
below.

VQE requires the calculation of expectation values of $\hat{H}$ for various
states in $\mathcal{F}$. The expectation value for the arbitrary state
$\ket{\Psi}$ is

$\braket{\Psi}{\hat{H}}{\Psi}=\sum_{\vec{n}_{A}\vec{n}_{A}^{\prime}}\alpha^{*}_{\vec{n}_{A}^{\prime}}\alpha_{\vec{n}_{A}}\braket{\Psi_{\vec{n}_{A}^{\prime}}}{C_{\vec{n}_{A}^{\prime}}\hat{H}C^{\dagger}_{\vec{n}_{A}}}{\Psi_{\vec{n}_{A}}}.$
(7)

To evaluate the expectation values on the right-hand side of this equation, we
start by grouping the terms in the Hamiltonian

$\hat{H}=\hat{H}_{A}+\hat{H}_{B}+\hat{H}_{AB}^{(t)}+\hat{H}_{AB}^{(v1)}+\hat{H}_{AB}^{(v2)}+\hat{H}_{AB}^{(v3)}$
(8)

such that $\hat{H}_{A}$ and $\hat{H}_{B}$ describe interactions exclusively in
subsystems $A$ and $B$, respectively, $\hat{H}_{AB}^{(t)}$ describe one-
particle interactions involving both subsystems, and $\hat{H}_{AB}^{(v1)}$,
$\hat{H}_{AB}^{(v2)}$, and $\hat{H}_{AB}^{(v3)}$ describes two-particle
interactions with one, two, and three fermionic operators acting on subsystem
$B$, respectively. The contributions to the expectation value from
$\hat{H}_{A}$ and $\hat{H}_{B}$, defined by

$\hat{H}_{\Upsilon}=\sum_{\mu\mu^{\prime}\in\boldsymbol{\mathrm{M}}_{\Upsilon}}t_{\mu\mu^{\prime}}c^{\dagger}_{\mu}c_{\mu^{\prime}}+\sum_{\mu\mu^{\prime}\nu\nu^{\prime}\in\boldsymbol{\mathrm{M}}_{\Upsilon}}v_{\mu\nu\mu^{\prime}\nu^{\prime}}c^{\dagger}_{\mu}c^{\dagger}_{\nu}c_{\nu^{\prime}}c_{\mu^{\prime}},$
(9)

can be obtained respectively from

$\braket{\Psi_{\vec{n}_{A}^{\prime}}}{C_{\vec{n}_{A}^{\prime}}\hat{H}_{A}C^{\dagger}_{\vec{n}_{A}}}{\Psi_{\vec{n}_{A}}}=\braket{\vec{n}_{A}^{\prime}}{\hat{H}_{A}}{\vec{n}_{A}}\braket{\Psi_{\vec{n}_{A}^{\prime}}}{\Psi_{\vec{n}_{A}}},$
(10)

where $\ket{\vec{n}_{A}}=C^{\dagger}_{\vec{n}_{A}}\ket{0}$, and

$\braket{\Psi_{\vec{n}_{A}^{\prime}}}{C_{\vec{n}_{A}^{\prime}}\hat{H}_{B}C^{\dagger}_{\vec{n}_{A}}}{\Psi_{\vec{n}_{A}}}=\delta_{{\vec{n}_{A}}\vec{n}_{A}^{\prime}}\braket{\Psi_{\vec{n}_{A}^{\prime}}}{\hat{H}_{B}}{\Psi_{\vec{n}_{A}}},$
(11)

where $\delta_{{\vec{n}_{A}}\vec{n}_{A}^{\prime}}$ is the Kronecker delta.
Similarly, the contributions from terms involving interactions in both
subspaces can be obtained from

$\displaystyle\braket{\Psi_{\vec{n}_{A}^{\prime}}}{C_{\vec{n}_{A}^{\prime}}\hat{H}_{AB}^{(t)}C^{\dagger}_{\vec{n}_{A}}}{\Psi_{\vec{n}_{A}}}=\sum_{\mu\in\boldsymbol{\mathrm{M}}_{A}}\sum_{\mu^{\prime}\in\boldsymbol{\mathrm{M}}_{B}}(-1)^{N_{A}}\Big{[}t_{\mu\mu^{\prime}}\braket{\vec{n}_{A}^{\prime}}{c^{\dagger}_{\mu}}{\vec{n}_{A}}\braket{\Psi_{\vec{n}_{A}^{\prime}}}{c_{\mu^{\prime}}}{\Psi_{\vec{n}_{A}}}-t_{\mu^{\prime}\mu}\braket{\vec{n}_{A}^{\prime}}{c_{\mu}}{\vec{n}_{A}}\braket{\Psi_{\vec{n}_{A}^{\prime}}}{c^{\dagger}_{\mu^{\prime}}}{\Psi_{\vec{n}_{A}}}\Big{]}$
(12)

where $N_{A}$ is the number of particles $\sum_{\mu}n_{\mu}$ in the
configuration $\vec{n}_{A}$, and

$\displaystyle\braket{\Psi_{\vec{n}_{A}^{\prime}}}{C_{\vec{n}_{A}^{\prime}}\hat{H}_{AB}^{(v1)}C^{\dagger}_{\vec{n}_{A}}}{\Psi_{\vec{n}_{A}}}=$
$\displaystyle\sum_{\mu\mu^{\prime}\nu\in\boldsymbol{\mathrm{M}}_{A}}\sum_{\nu^{\prime}\in\boldsymbol{\mathrm{M}}_{B}}(-1)^{N_{A}}\Big{[}\big{(}v_{\mu\nu\nu^{\prime}\mu^{\prime}}-v_{\mu\nu\mu^{\prime}\nu^{\prime}}\big{)}\braket{\vec{n}_{A}^{\prime}}{c^{\dagger}_{\mu}c^{\dagger}_{\nu}c_{\mu^{\prime}}}{\vec{n}_{A}}\braket{\Psi_{\vec{n}_{A}^{\prime}}}{c_{\nu^{\prime}}}{\Psi_{\vec{n}_{A}}}$
$\displaystyle+\big{(}v_{\mu\nu^{\prime}\mu^{\prime}\nu}-v_{\nu^{\prime}\mu\mu^{\prime}\nu}\big{)}\braket{\vec{n}_{A}^{\prime}}{c^{\dagger}_{\mu}c_{\nu}c_{\mu^{\prime}}}{\vec{n}_{A}}\braket{\Psi_{\vec{n}_{A}^{\prime}}}{c^{\dagger}_{\nu^{\prime}}}{\Psi_{\vec{n}_{A}}}\Big{]},$
(13)
$\displaystyle\braket{\Psi_{\vec{n}_{A}^{\prime}}}{C_{\vec{n}_{A}^{\prime}}\hat{H}_{AB}^{(v2)}C^{\dagger}_{\vec{n}_{A}}}{\Psi_{\vec{n}_{A}}}=$
$\displaystyle\sum_{\mu\mu^{\prime}\in\boldsymbol{\mathrm{M}}_{A}}\sum_{\nu\nu^{\prime}\in\boldsymbol{\mathrm{M}}_{B}}\Big{[}\big{(}v_{\mu\nu\mu^{\prime}\nu^{\prime}}-v_{\mu\nu\nu^{\prime}\mu^{\prime}}-v_{\nu\mu\mu^{\prime}\nu^{\prime}}+v_{\nu\mu\nu^{\prime}\mu^{\prime}}\big{)}\braket{\vec{n}_{A}^{\prime}}{c^{\dagger}_{\mu}c_{\mu^{\prime}}}{\vec{n}_{A}}\braket{\Psi_{\vec{n}_{A}^{\prime}}}{c^{\dagger}_{\nu}c_{\nu^{\prime}}}{\Psi_{\vec{n}_{A}}}$
$\displaystyle+v_{\mu\mu^{\prime}\nu\nu^{\prime}}\braket{\vec{n}_{A}^{\prime}}{c^{\dagger}_{\mu}c^{\dagger}_{\mu^{\prime}}}{\vec{n}_{A}}\braket{\Psi_{\vec{n}_{A}^{\prime}}}{c_{\nu^{\prime}}c_{\nu}}{\Psi_{\vec{n}_{A}}}+v_{\nu^{\prime}\nu\mu^{\prime}\mu}\braket{\vec{n}_{A}^{\prime}}{c_{\mu}c_{\mu^{\prime}}}{\vec{n}_{A}}\braket{\Psi_{\vec{n}_{A}^{\prime}}}{c^{\dagger}_{\nu^{\prime}}c^{\dagger}_{\nu}}{\Psi_{\vec{n}_{A}}}\Big{]},$
(14)
$\displaystyle\braket{\Psi_{\vec{n}_{A}^{\prime}}}{C_{\vec{n}_{A}^{\prime}}\hat{H}_{AB}^{(v3)}C^{\dagger}_{\vec{n}_{A}}}{\Psi_{\vec{n}_{A}}}=$
$\displaystyle\sum_{\mu\in\boldsymbol{\mathrm{M}}_{A}}\sum_{\mu^{\prime}\nu\nu^{\prime}\in\boldsymbol{\mathrm{M}}_{B}}(-1)^{N_{A}}\Big{[}\big{(}v_{\mu\nu\mu^{\prime}\nu^{\prime}}-v_{\nu\mu\mu^{\prime}\nu^{\prime}}\big{)}\braket{\vec{n}_{A}^{\prime}}{c^{\dagger}_{\mu}}{\vec{n}_{A}}\braket{\Psi_{\vec{n}_{A}^{\prime}}}{c^{\dagger}_{\nu}c_{\nu^{\prime}}c_{\mu^{\prime}}}{\Psi_{\vec{n}_{A}}}$
$\displaystyle+\big{(}v_{\nu^{\prime}\nu\mu^{\prime}\mu}-v_{\nu^{\prime}\nu\mu\mu^{\prime}}\big{)}\braket{\vec{n}_{A}^{\prime}}{c_{\mu}}{\vec{n}_{A}}\braket{\Psi_{\vec{n}_{A}^{\prime}}}{c^{\dagger}_{\nu^{\prime}}c^{\dagger}_{\nu}c_{\mu^{\prime}}}{\Psi_{\vec{n}_{A}}}\Big{]}.$
(15)

Every factor of the form
$\braket{\vec{n}^{\prime}_{A}}{\hat{O}_{A}}{\vec{n}_{A}}$, where $\hat{O}_{A}$
is some operator acting on states in subspace $A$, can be expressed as a
Kronecker delta function eliminating all but one term in the sum over
$\vec{n}_{A}^{\prime}$ in Eq. (7). These delta functions will be found
classically. The other type of factors of the form
$\braket{\Psi_{\vec{n}_{A}^{\prime}}}{\hat{O}_{B}}{\Psi_{\vec{n}_{A}}}$, where
$\hat{O}_{B}$ is some operator acting on states in subspace $B$, is computed
using our modified VQE algorithm.

First, let us require that the mapping from $\mathcal{F}_{B}$ to the Hilbert
space for the quantum register takes $\ket{\Psi_{\vec{n}_{A}}}$ to
$\ket{\Phi_{\vec{n}_{A}}}$ for all configurations $\vec{n}_{A}$. Moreover, if
we define $\hat{\mathcal{O}}$ such that $\hat{O}_{B}\mapsto\hat{\mathcal{O}}$,
we find that the factors
$\braket{\Psi_{\vec{n}_{A}^{\prime}}}{\hat{O}_{B}}{\Psi_{\vec{n}_{A}}}$ map to
$\braket{\Phi_{\vec{n}^{\prime}_{A}}}{\hat{\mathcal{O}}}{\Phi_{\vec{n}_{A}}}$.
In the case, $\vec{n}_{A}^{\prime}=\vec{n}_{A}$, the latter become expectation
values that can be determined using the standard VQE algorithm by preparing
the state $\ket{\Phi_{\vec{n}_{A}}}$ on the quantum register and performing
appropriate measurements. For $\vec{n}_{A}^{\prime}\neq\vec{n}_{A}$, this
approach is not possible, so we instead introduce an ancilla qubit and prepare
the state

$\ket{\Phi_{\vec{n}_{A}\vec{n}^{\prime}_{A}}}=\frac{1}{\sqrt{2}}\Big{(}\ket{0;\Phi_{\vec{n}_{A}}}+\ket{1;\Phi_{\vec{n}^{\prime}_{A}}}\Big{)}$
(16)

on the quantum register. This allows us to determine

$\braket{\Phi_{\vec{n}^{\prime}_{A}}}{\hat{\mathcal{O}}}{\Phi_{\vec{n}_{A}}}=\braket{\Phi_{\vec{n}_{A}\vec{n}^{\prime}_{A}}}{(X-iY)\otimes\hat{\mathcal{O}}}{\Phi_{\vec{n}_{A}\vec{n}^{\prime}_{A}}},$
(17)

where $X$ and $Y$ are Pauli operators, by computing the latter expectation
value using the same technique as in the standard VQE algorithm for the non-
ancilla qubits.

###  II.2 VQE on blocks

In standard VQE, the ansatz is a parameterized unitary operator
$\hat{U}(\vec{\theta})$, which when acting on the initial state
$\ket{\vec{0}}=\ket{00\ldots 0}$ of the quantum register produce $\ket{\Phi}$
representing the many-particle state $\ket{\Psi}\in\mathcal{F}$. Herein, we
instead want to represent $\ket{\Psi_{\vec{n}_{A}}}\in\mathcal{F}_{B}$. As the
states on subsystem $B$ are dependent on the particle configuration of
subsystem $A$, we are interested in the $\vec{n}_{A}$-dependent ansatz states
$\hat{U}(\vec{\theta}_{\vec{n}_{A}})\ket{\vec{0}}$ to describe the quantum
states

$\ket{\Phi_{\vec{n}_{A}}}\approx\hat{U}(\vec{\theta}_{\vec{n}_{A}})\ket{\vec{0}}$
(18)

on the quantum register corresponding to $\ket{\Psi_{\vec{n}_{A}}}$. The
energy expectation values can be calculated as described in sec II.1. These
values are then used to optimize, not only each $\vec{\theta}_{\vec{n}_{A}}$
but also each $\alpha_{\vec{n}_{A}}$ subject to the normalization constraint
$\sum_{\vec{n}_{A}}|\alpha_{\vec{n}_{A}}|^{2}=1$. As noted above, to calculate

$\braket{\Phi_{\vec{n}^{\prime}_{A}}}{\hat{\mathcal{O}}}{\Phi_{\vec{n}_{A}}}\approx\braket{\vec{0}}{\hat{U}^{\dagger}(\vec{\theta}_{\vec{n}_{A}^{\prime}})\hat{\mathcal{O}}\hat{U}(\vec{\theta}_{\vec{n}_{A}})}{\vec{0}}$
(19)

we use an ancilla qubit to form the state
$\ket{\Phi_{\vec{n}_{A}\vec{n}_{A}^{\prime}}}$ in Eq. (16).

Figure 1: Unitary which creates the
$\ket{+;\Phi_{\vec{n}_{A}\vec{n}_{A}^{\prime}}}$ state used for calculating
off diagonal brackets. In this example, each block is described by three
qubits, however, the unitary can easily be generalized to any number of
qubits.

What we have not yet addressed is how to create this state using quantum
gates. Specifically, we need to implement a controlled unitary operator
$\hat{C}_{\hat{U}(\vec{\theta}_{\vec{n}_{A}})}$ that applies
$\hat{U}(\vec{\theta}_{\vec{n}_{A}})$, if and only if the ancilla qubit is in
the state $\ket{1}$. We can then obtain the desired state by using

$\ket{+;\Phi_{\vec{n}_{A}\vec{n}_{A}^{\prime}}}\approx\hat{C}_{\hat{U}(\vec{\theta}_{\vec{n^{\prime}}_{A}})}\hat{C}_{\hat{U}^{\dagger}(\vec{\theta}_{\vec{n}_{A}})}\hat{[}\mathcal{H}_{0}\otimes
U(\vec{\theta}_{\vec{n}_{A}})]\ket{0;\vec{0}}$ (20)

where $\ket{+}=\big{(}\ket{0}+\ket{1}\big{)}/\sqrt{2}$ and
$\hat{\mathcal{H}_{0}}$ is the Hadamard gate acting on the ancilla qubit, see
Fig 1. The $\hat{C}_{\hat{U}(\vec{\theta})}$ gates can be constructed using
$C_{\text{NOT}}$ gates. However, $C_{\text{NOT}}$ gates are a major source of
error. Since often $|\vec{\theta}_{\vec{n}_{A}}|\ll\pi$, we use scaled two-
qubit cross-resonance gates which are significantly more efficient at small
angles [25].

##  III Hubbard Model

To demonstrate our approach, we focus below on the $L$-site Hubbard model of a
system with identical spin-$1/2$ particles described by the Hamiltonian

$\hat{H}=\epsilon\sum_{\sigma}\hat{N}_{\sigma}+t\sum_{\sigma}\hat{T}_{\sigma}+U\sum_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow},$
(21)

where

$\displaystyle\hat{N}_{\sigma}$ $\displaystyle=\sum_{i}\hat{n}_{i\sigma},$
(22a) $\displaystyle\hat{T}_{\sigma}$
$\displaystyle=\sum_{i}\left(c^{\dagger}_{i\sigma}c_{i+1\sigma}+\text{h.c.}\right),$
(22b)

and the fermionic operators, including the number operators
$\hat{n}_{i\sigma}=c^{\dagger}_{i\sigma}c_{i\sigma}$, are identified by
single-particle indices split into site indices $i\in\\{1,2,\ldots,L\\}$ and
spin indices $\sigma\in\\{\uparrow,\downarrow\\}$, $\epsilon$ is the
electrochemical potential, and $t$ and $U$ are the hopping and onsite Coulomb
interaction coefficients, respectively. The sites are arranged in a one-
dimensional ring with the boundary condition $c_{L\sigma}=c_{0\sigma}$.

###  III.1 Splitting the spin sectors

Placing spin-up and spin-down particles in subsystems $A$ and $B$,
respectively, allows us to identify the particle configurations as
$\vec{n}_{\sigma}=\big{(}n_{1\sigma},n_{2\sigma},\ldots,n_{L\sigma}\big{)}$
and write the many-particle state as

$\ket{\Psi}=\sum_{\vec{n}_{\uparrow}\vec{n}_{\downarrow}}\alpha_{\vec{n}_{\uparrow}}\beta_{\vec{n}_{\uparrow}\vec{n}_{\downarrow}}\ket{\vec{n}},$
(23)

where the normalization yields
$\sum_{\vec{n}_{\uparrow}}|\alpha_{\vec{n}_{\uparrow}}|^{2}=1$ and
$\sum_{\vec{n}_{\downarrow}}|\beta_{\vec{n}_{\uparrow}\vec{n}_{\downarrow}}|^{2}=1$,
the latter for all $\vec{n}_{\uparrow}$. Without any restrictions on the
vector spaces for the two subsystems, the number of spin-up states and the
number of spin-down states are both $2^{L}$.

In order to perform VQE we need to calculate the energy expectation value
which is used as the cost function. We will calculate the spin-up sector of
the expectation value classically. This will leave us with a set of
$\vec{n}_{\uparrow}$-dependent spin-down states

$\ket{\Psi_{\vec{n}_{\uparrow}}}=\sum_{\vec{n}_{\downarrow}}\beta_{\vec{n}_{\uparrow}\vec{n}_{\downarrow}}\ket{\vec{n}_{\downarrow}}$
(24)

in the Fock space $\mathcal{F}_{\downarrow}$. We then use the quantum device
to determine the coefficients
$\beta_{\vec{n}_{\uparrow}\vec{n}_{\downarrow}}$. The expectation value of the
Hamiltonian can be expressed as

$\displaystyle\bra{\Psi}\hat{H}\ket{\Psi}=$ $\displaystyle\leavevmode\nobreak\
\epsilon\sum_{\vec{n}_{\uparrow}}|\alpha_{\vec{n}_{\uparrow}}|^{2}\Big{(}N_{\vec{n}_{\uparrow}}+\braket{\Psi_{\vec{n}_{\uparrow}}}{\hat{N}_{\downarrow}}{\Psi_{\vec{n}_{\uparrow}}}\Big{)}$
$\displaystyle+t\sum_{\vec{n}_{\uparrow}\vec{n}_{\uparrow}^{\prime}}\alpha_{\vec{n}_{\uparrow}^{\prime}}^{*}\alpha_{\vec{n}_{\uparrow}}T_{\vec{n}_{\uparrow}^{\prime}\vec{n}_{\uparrow}}\braket{\Psi_{\vec{n}_{\uparrow}^{\prime}}}{\Psi_{\vec{n}_{\uparrow}}}$
$\displaystyle+t\sum_{\vec{n}_{\uparrow}}|\alpha_{\vec{n}_{\uparrow}}|^{2}\braket{\Psi_{\vec{n}_{\uparrow}}}{\hat{T}_{\downarrow}}{\Psi_{\vec{n}_{\uparrow}}}$
$\displaystyle+U\sum_{\vec{n}_{\uparrow}}|\alpha_{\vec{n}_{\uparrow}}|^{2}\sum_{i}n_{i\uparrow}\braket{\Psi_{\vec{n}_{\uparrow}}}{\hat{n}_{i\downarrow}}{\Psi_{\vec{n}_{\uparrow}}},$
(25)

where
$N_{\vec{n}_{\uparrow}}=\bra{\vec{n}_{\uparrow}}\hat{N}_{\uparrow}\ket{\vec{n}_{\uparrow}}$
and
$T_{\vec{n}_{\uparrow}^{\prime}\vec{n}_{\uparrow}}=\bra{\vec{n}_{\uparrow}^{\prime}}\hat{T}_{\uparrow}\ket{\vec{n}_{\uparrow}}$.
The spin-up factors $n_{i\uparrow}$, $N_{\vec{n}_{\uparrow}}$, and
$T_{\vec{n}_{\uparrow}^{\prime}\vec{n}_{\uparrow}}$ will be calculated on the
classical computer and stored for use in our modified VQE. The spin-down
factors
$\braket{\Psi_{\vec{n}_{\uparrow}^{\prime}}}{\Psi_{\vec{n}_{\uparrow}}}$,
$\braket{\Psi_{\vec{n}_{\uparrow}}}{\hat{T}_{\downarrow}}{\Psi_{\vec{n}_{\uparrow}}}$,
and
$\braket{\Psi_{\vec{n}_{\uparrow}}}{\hat{n}_{i\downarrow}}{\Psi_{\vec{n}_{\uparrow}}}$
will be calculated using the quantum computer.

The expectation values
$\braket{\Psi_{\vec{n}_{\uparrow}}}{\hat{T}_{\downarrow}}{\Psi_{\vec{n}_{\uparrow}}}$,
and
$\braket{\Psi_{\vec{n}_{\uparrow}}}{\hat{n}_{i\downarrow}}{\Psi_{\vec{n}_{\uparrow}}}$
are straightforward to calculate using the quantum device. First, map fermonic
states $\ket{\Psi_{\vec{n}_{\uparrow}}}$ onto states of the quantum computer
$\ket{\Phi_{\vec{n}_{\uparrow}}}$ and decompose $\hat{T}_{\downarrow}$ and
$\hat{n}_{i\downarrow}$ into Pauli operators. Then measure the expectation
value of the Pauli matrices using standard techniques. The bracket
$\braket{\Psi_{\vec{n}_{\uparrow}^{\prime}}}{\Psi_{\vec{n}_{\uparrow}}}$ is
harder to calculate since it does not have the form of an expectation value.
For this bracket we add an ancilla qubit as in Eqns (16) and (17).

By calculating the spin up sector on the classical computer, we have reduced
the number of qubits required for the quantum computation from $2L$ to $L+1$.
An additional advantage of splitting between the spin sectors is that it
exploits the sparsity of the Hamiltonian. Everything is diagonal in the spin-
up sector except for the elements
$T_{\vec{n}_{\uparrow}^{\prime}\vec{n}_{\uparrow}}$. Furthermore,
$T_{\vec{n}_{\uparrow}^{\prime}\vec{n}_{\uparrow}}$ is sparse with many blocks
having all elements equal to zero. Because we calculate the
$T_{\vec{n}_{\uparrow}^{\prime}\vec{n}_{\uparrow}}$ terms before running the
spin-down sector on the quantum processor, we can ensure that only blocks with
non-zero elements are run on the quantum device.

The price of calculating the spin-up sector classically is that we must sum
over $2^{L}$ spin-up configurations. This number can be reduced, however, by
restricting the total number spin-up fermions.

###  III.2 Fixing the number of spin-up fermions

The Hubbard model, Eq. (21), conserves the total number of spin up and spin
down fermions. Therefore, we can reduce the size of the classical computation
by fixing the number of spin-up electrons $\hat{N}_{\uparrow}\rightarrow
N_{\uparrow}$. In general there are

$\Lambda={L\choose N_{\uparrow}}$ (26)

spin up states.

This procedure still requires $Q=L+1$ qubits to solve the down spins. However,
the complexity of the classical computation is greatly reduced. There is no
disadvantage to fixing the number of spin-up fermions.

###  III.3 Fixing the spin-down parity

A single qubit can be removed by fixing the total parity of the spin-down
sector [26]. This can be done by changing the qubit encoding of our fermion
operators. A natural choice for for encoding is the Jordan–Wigner
transformation [27] for which

$c_{i\downarrow}\mapsto\frac{1}{2}(X_{i}+iY_{i})\prod_{j=0}^{L-1}Z_{j}$ (27)

where $X_{i}$, $Y_{i}$, and $Z_{i}$ are the three Pauli matrices acting on
qubit $i$. This is convenient because of the simplicity of interpretation, a
fermion operator at site $i$ flips the state of qubit $i$. However, it uses
more qubits than is necessary as fermion parity is always conserved.

We can remove a qubit by changing the encoding to, for example, the parity
encoding [28, 29]. Here the fermion operators are given by

$c_{i\downarrow}\mapsto\frac{1}{2}(Z_{i-1}X_{i}+iY_{i})\prod_{j=i+1}^{L}X_{j}$
(28)

Note that in this encoding the total parity maps to qubit zero. Therefore, we
can remove this qubit from the encoding and set the value of the parity by
hand. The downside is that the Pauli decomposition of $\hat{T}_{\downarrow}$
and $\hat{n}_{i\downarrow}$ contains a greater number of Pauli terms. For the
most part, the Pauli-strings only grow by one qubit, however, the boundary
term in $\hat{T}_{\downarrow}$ (i.e.
$c^{\dagger}_{L-1\downarrow}c_{0\downarrow}+\text{h.c.}$) is mapped onto
Pauli-strings of length $L$.

This situation can be improved by selecting yet another encoding, namely the
Bravyi-Kitaev encoding [29, 28]. The Bravyi-Kitaev transformation balances
information about locality with information about parity so that the length of
the Pauli-strings for local fermion operators is $\mathcal{O}(\log L)$ while
also having the property that the total parity is mapped to the first qubit.
The exact form of the transformation is complicated and can be found in the
literature [29, 28, 26] so we will not write it here. For the one-dimensional
Hubbard model, the Pauli-strings are still larger on average than in the
Jordan–Wigner transformation. One has to make a hardware-based decision for
whether it is better to have fewer qubits or shorter Pauli-strings.

###  III.4 Fixing the number of spin down fermions

Since the Hubbard model conserves the number of spin-down electrons as well as
the number of spin-up electrons, we can also fix the total number of spin down
fermions to $\hat{N}_{\downarrow}\rightarrow N_{\downarrow}$. By fixing the
total number of spin-down electrons, we reduce the Hilbert space of
$\hat{T}_{\downarrow}$ and $\hat{n}_{i\downarrow}$ which means we can write
$\ket{\Phi_{\vec{n}_{\uparrow}}}$ using fewer qubits from $Q=L+1$ to

$Q=\bigg{\lceil}\log_{2}{L\choose N_{\downarrow}}\bigg{\rceil}+1$ (29)

This can be a major gain, however, it does increase the amount of work
required for the classical computer. In fact, the number of Pauli matrices
required to describe $\hat{T}_{\downarrow}$ is no longer guaranteed to be
linear with system size.

To find the Pauli decomposition of $\hat{T}_{\downarrow}$ and
$\hat{n}_{i\downarrow}$ we have to first create the matrix representation of
both in the basis were both spins are fixed. This step is done classically and
can be done efficiently by calling the states using the binary representation
of integers then we can find matrix elements using binary operations, see [30]
for details. Once the matrix has been specified, the Pauli coefficients can be
found in $\mathcal{O}\left(L^{2}\log_{2}L\right)$ operations, or less [31].

Figure 2: Simulation details. (A) depiction of the four site Hubbard Model at
half filling. The orange circles represent sites and the blue lines show
connectivity. The arrows represent spin-up and spin-down electrons. The
placement of the electrons in the figure does not necessarily reflect the
actual location of electrons. (B) two iterations of the ansatz used for our
simulations. The single qubit rotation during iteration $n$ on qubit $q$ is
parameterized by $\theta_{nq}$. In the simulation we use two iterations of the
ansatz. (C) circuit to apply a single iteration of the ansatz to off diagonal
blocks. The set of rotations on block $A$ are $\\{\theta_{Anq}\\}$ and those
on block $B$ are $\\{\theta_{Bnq}\\}$. The difference between two phases from
different blocks is labeled as $\Delta\theta_{nq}=\theta_{Bnq}-\theta_{Anq}$.
This circuit applies the block $A$ rotations if the ancilla qubit is off and
the block $B$ rotations if the ancilla qubit is on. The controlled-$R_{y}$
gates are generated using parameterized $R_{ZX}$ gates [25].

##  IV Simulation of the four site Hubbard model at half filling fixing the
number of both spin sectors

In this section, we present results for a simulation of the four site ($L=4$)
Hubbard model [Eq. (21)] at half filling ($N_{\uparrow}=N_{\downarrow}=2$).
See Fig. 2(A). The number of spin-up states is $[L\leavevmode\nobreak\
\text{choose}\leavevmode\nobreak\ N_{\uparrow}]=6$. The number of qubits we
need is $Q=4$ (see Eq. (29)).

We use the hardware-efficient ansatz [19, 24] where single qubit gates are
sandwiched between CNOT gates. See Fig. 2(B). We use two iterations of the
ansatz for our simulation.

Figure 3: Results from the simulation of our block VQE method on the four site
Hubbard model at half filling, fixing the number of particles in both spin
sectors. The results of block VQE are plotted as green circles. The exact
ground state is plotted as the solid black curve. The mean field result is
plotted as the solid blue curve.

In order to calculate the off diagonal blocks, we add an ancilla qubit. For
the particular ansatz we use, the general form of the off diagonal unitary in
Eq. 20 can be simplified, see Fig 2(C).

The results of the simulations are plotted as the green circles in Fig. 3. The
results are compared against the exact ground state (solid black curve) and
the mean field solution (solid blue curve). For large values of the
interaction strength $U$, our VQE is much better than the mean field
calculation.

##  V Conclusion

We have described a method which splits the work of solving many-particle
Hamiltonians between classical and quantum computers. In our main example,
this is accomplished by splitting the spin-sectors and solving one of them
classically before solving the other on the quantum device. Splitting the
spins in this way immediately reduces the number of qubits needed for the
quantum calculation by half. A further reduction of qubits can be achieved by
evoking parity, and even more qubits can be removed by invoking spin number
conservation. Each time qubits are removed from the quantum calculation, the
classical calculation becomes more involved. Thus, we propose a ladder of
strategies going from most classically intensive and least quantum intensive
to least classically intensive and most quantum intensive. As quantum
processors improve we will be able to climb up this ladder.

*

##  Appendix A System specification

The coefficients describing the system can be computed using the integral
forms:

$t_{\mu\mu^{\prime}}=\int\psi_{\mu}^{*}(x)\hat{T}(x)\psi_{\mu^{\prime}}(x)\mathrm{d}x$
(30)

and

$v_{\mu\nu\mu^{\prime}\nu^{\prime}}=\iint\psi_{\mu}^{*}(x)\psi_{\nu}^{*}(x^{\prime})\hat{V}(x,x^{\prime})\psi_{\mu^{\prime}}(x)\psi_{\nu^{\prime}}(x^{\prime})\mathrm{d}x\mathrm{d}x^{\prime},$
(31)

where $\psi_{\mu}(x)$ are one-particle basis functions on some coordinate
space and $\hat{T}$ and $\hat{V}$ are first-quantized one- and two-particle
operators, respectively.

For an isolated molecular system within the Born–Oppenheimer approximation,
the integral coefficients for the electronic system can be expressed as

$\displaystyle t_{\mu\mu^{\prime}}$
$\displaystyle=\sum_{\sigma}\int\psi_{\mu}^{*}(\vec{r},\sigma)\Bigg{[}-\frac{\hbar^{2}}{2m}\nabla^{2}-\frac{e^{2}}{4\pi\varepsilon_{0}}\sum_{n}\frac{Z_{n}}{\vec{r}-\vec{R}_{n}}\Bigg{]}\psi_{\mu^{\prime}}(\vec{r},\sigma)\,\mathrm{d}\vec{r},$
(32) $\displaystyle v_{\mu\nu\mu^{\prime}\nu^{\prime}}$
$\displaystyle=\sum_{\sigma\sigma^{\prime}}\iint\psi_{\mu}^{*}(\vec{r},\sigma)\psi_{\nu}^{*}(\vec{r}^{\prime},\sigma^{\prime})\frac{e^{2}}{4\pi\varepsilon_{0}|\vec{r}-\vec{r}^{\prime}|}\psi_{\mu^{\prime}}(\vec{r},\sigma)\psi_{\nu^{\prime}}(\vec{r}^{\prime},\sigma^{\prime})\,\mathrm{d}\vec{r}\,\mathrm{d}\vec{r}^{\prime}$
(33)

where $\vec{r}$ and $\sigma$ are spatial and spin coordinates, respectively,
$\hbar$, $m$, $e$, and $\epsilon_{0}$ are the reduced Planck constant, the
electron mass, the electron charge, and the vacuum permittivity, respectively,
and $Z_{n}$ and $\vec{R}_{n}$ are the charges and locations of the nuclei $n$.

###### Acknowledgements.

This work has been supported by the Office of Naval Research (ONR) through the
U.S. Naval Research Laboratory (NRL). J.P.T.S. thanks the National Research
Council Research Associateship Programs. We acknowledge quantum computing
resources from IBM through a collaboration with the Air Force Research
Laboratory (AFRL).

## References

  * Feynman [1982] R. P. Feynman, Int. J. Theor Phys. 21, 467 (1982).
  * Lloyd [1996] S. Lloyd, Science (New York, N.Y.) 273, 1073 (1996).
  * Boixo  _et al._ [2018] S. Boixo, S. V. Isakov, V. N. Smelyanskiy, R. Babbush, N. Ding, Z. Jiang, M. J. Bremner, J. M. Martinis, and H. Neven, Nature Physics 14, 595–600 (2018).
  * Shor [1997] P. W. Shor, SIAM Journal on Computing 26, 1484–1509 (1997).
  * Cirac and Zoller [2012] J. I. Cirac and P. Zoller, nature physics 8, 264–266 (2012).
  * Georgescu  _et al._ [2014] I. Georgescu, S. Ashhab, and F. Nori, Reviews of Modern Physics 86, 153–185 (2014).
  * Tarruell and Sanchez-Palencia [2019] L. Tarruell and L. Sanchez-Palencia, (2019), arXiv:1809.00571 [cond-mat.quant-gas] .
  * Salfi  _et al._ [2016] J. Salfi, J. A. Mol, R. Rahman, G. Klimeck, M. Y. Simmons, L. C. L. Hollenberg, and S. Rogge, Nature Communications 7, 11342 (2016).
  * Hensgens  _et al._ [2017] T. Hensgens, T. Fujita, L. Janssen, X. Li, C. J. Van Diepen, C. Reichl, W. Wegscheider, S. Das Sarma, and L. M. K. Vandersypen, Nature 548, 70–73 (2017).
  * Singha  _et al._ [2011] A. Singha, M. Gibertini, B. Karmakar, S. Yuan, M. Polini, G. Vignale, M. I. Katsnelson, A. Pinczuk, L. N. Pfeiffer, K. W. West, and et al., Science 332, 1176–1179 (2011).
  * Smith  _et al._ [2019] A. Smith, M. S. Kim, F. Pollmann, and J. Knolle, npj Quantum Information 5, 106 (2019).
  * Bakr  _et al._ [2009] W. S. Bakr, J. I. Gillen, A. Peng, S. Fölling, and M. Greiner, Nature 462, 74–77 (2009).
  * Schreiber  _et al._ [2015] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Luschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Science 349, 842–845 (2015).
  * Choi  _et al._ [2016] J.-y. Choi, S. Hild, J. Zeiher, P. Schauss, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch, and C. Gross, Science 352, 1547–1552 (2016).
  * Bordia  _et al._ [2017] P. Bordia, H. Lüschen, S. Scherg, S. Gopalakrishnan, M. Knap, U. Schneider, and I. Bloch, Physical Review X 7, 041047 (2017).
  * Mitra  _et al._ [2017] D. Mitra, P. T. Brown, E. Guardado-Sanchez, S. S. Kondov, T. Devakul, D. A. Huse, P. Schauß, and W. S. Bakr, Nature Physics 14, 173–177 (2017).
  * 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, 4213 (2014).
  * McClean  _et al._ [2016] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, New Journal of Physics 18, 023023 (2016).
  * Kandala  _et al._ [2017] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242–246 (2017).
  * Yuan  _et al._ [2021] X. Yuan, J. Sun, J. Liu, Q. Zhao, and Y. Zhou, Physical Review Letters 127, 10.1103/physrevlett.127.040501 (2021).
  * Peng  _et al._ [2020] T. Peng, A. W. Harrow, M. Ozols, and X. Wu, Physical Review Letters 125, 10.1103/physrevlett.125.150504 (2020).
  * Huang  _et al._ [2020] H.-Y. Huang, R. Kueng, and J. Preskill, Nature Physics 16, 1050–1057 (2020).
  * Paini  _et al._ [2021] M. Paini, A. Kalev, D. Padilha, and B. Ruck, Quantum 5, 413 (2021).
  * Fischer and Gunlycke [2019] S. A. Fischer and D. Gunlycke, (2019), arXiv:1907.01493 [quant-ph] .
  * Stenger  _et al._ [2021] J. P. T. Stenger, N. T. Bronn, D. J. Egger, and D. Pekker, Phys. Rev. Research 3, 033171 (2021).
  * Bravyi  _et al._ [2017] S. Bravyi, J. M. Gambetta, A. Mezzacapo, and K. Temme, (2017), arXiv:1701.08213 [quant-ph] .
  * Jordan and Wigner [1928] P. Jordan and E. Wigner, Zeitschrift fur Physik 47, 631 (1928).
  * Seeley  _et al._ [2012] J. T. Seeley, M. J. Richard, and P. J. Love, The Journal of Chemical Physics 137, 224109 (2012).
  * Bravyi and Kitaev [2002] S. B. Bravyi and A. Y. Kitaev, Annals of Physics 298, 210–226 (2002).
  * Lin  _et al._ [1993] H. Lin, J. Gubernatis, H. Gould, and J. Tobochnik, Computers in Physics 7, 400 (1993).
  * Gunlycke  _et al._ [2020] D. Gunlycke, M. C. Palenik, A. R. Emmert, and S. A. Fischer, (2020), arXiv:2011.08942 [quant-ph] .