# Recursive relations and quantum eigensolver algorithms within modified
Schrieffer–Wolff transformations for the Hubbard dimer

Quentin Marécat  ICGM, Université de Montpellier, CNRS, ENSCM, Montpellier,
France  Bruno Senjean  ICGM, Université de Montpellier, CNRS, ENSCM,
Montpellier, France  Matthieu Saubanère<EMAIL_ADDRESS>ICGM, Université de Montpellier, CNRS, ENSCM, Montpellier, France

###### Abstract

We derive recursive relations for the Schrieffer–Wolff (SW) transformation
applied to the half-filled Hubbard dimer. While the standard SW transformation
is set to block-diagonalize the transformed Hamiltonian solely at the first
order of perturbation, we infer from recursive relations two types of
modifications, variational or iterative, that approximate, or even enforce for
the homogeneous case, the desired block-diagonalization at infinite order of
perturbation. The modified SW unitary transformations are then used to design
an test quantum algorithms adapted to the noisy and fault-tolerant era. This
work paves the way toward the design of alternative quantum algorithms for the
general Hubbard Hamiltonian.

##  I Introduction

By describing the competition between kinetically induced electron
delocalization and electron localization due to the Coulomb repulsion, the
non-trivial Hubbard model remains one of the most challenging systems in
condensed matter physics Hubbard (1963). Indeed, despite its simplicity, no
general and analytic solution exists. Besides exact results at certain limits
such as the Nagaoka theorem Nagaoka (1965) close to half-band filling or the
Bethe Ansatz Lieb and Wu (1968) in one dimension, different approximations,
strategies and numerical algorithms have been designed to solve this
cornerstone problem on classical computers. More precisely, one could mention
density functional López-Sandoval and Pastor (2002); Lima  _et al._ (2003) or
Green’s functions Georges  _et al._ (1996); Sénéchal  _et al._ (2000);
Potthoff (2003) based theories, renormalization methods White (1992) or more
recently divide and conquer strategies Knizia and Chan (2012); Sekaran  _et
al._ (2021), to cite but a few.

In that context, the emergence of quantum computers has revived the hope of
obtaining accurate physically relevant quantities for any dimension, size,
regime and filling. Indeed, a growing interest on developing quantum
algorithms to solve the Hubbard model emerges from the literature Wecker  _et
al._ (2015a, b); Kivlichan  _et al._ (2018); Reiner  _et al._ (2019);
Montanaro and Stanisic (2020); Cai (2020); Cade  _et al._ (2020); Mineh and
Montanaro (2022); Martin  _et al._ (2022); Stanisic  _et al._ (2022);
Dallaire-Demers  _et al._ (2019, 2020); Suchsland  _et al._ (2022); Gard and
Meier (2022); Kivlichan  _et al._ (2020); Campbell (2022); Clinton  _et al._
(2021). On the one hand, most of the proposed algorithms targets Noisy
Intermediate Scale Quantum (NISQ) devices and relies mainly on hybrid
classical/quantum strategies such as the Variational Quantum Eigensolver (VQE)
Peruzzo  _et al._ (2014); Bharti  _et al._ (2022). Roughly speaking, it
consists in applying a parameterized unitary transformation on an easy-to-
prepare initial state, generally the Hartree–Fock state, on the quantum device
while the variational parameters are optimized on a classical computer.
Several type of Ansatz have been proposed to design this unitary
transformation, either physically motivated such as the variational
Hamiltonian Ansatz Wecker  _et al._ (2015b); Kivlichan  _et al._ (2018);
Reiner  _et al._ (2019); Montanaro and Stanisic (2020); Cai (2020); Cade  _et
al._ (2020); Mineh and Montanaro (2022); Martin  _et al._ (2022); Stanisic
_et al._ (2022) and the unitary coupled cluster Ansatz Dallaire-Demers  _et
al._ (2019), or hardware efficient ones Dallaire-Demers  _et al._ (2020);
Suchsland  _et al._ (2022); Gard and Meier (2022). Most of these approaches,
as they are based on an initial Hartree–Fock state, are particularly relevant
for the weakly correlated regime. Cade  _et al._ (2020); Martin  _et al._
(2022) In any case, a compromise between the desired accuracy and the
computational cost has to be reached. It depends in particular on the Ansatz
circuit depth, the number of CNOT gates and the number of variational
parameters, for which the development of improved or new types of Ansatz is
needed. On the second hand, some algorithms target long-term expected fault-
tolerant devices Kivlichan  _et al._ (2020); Campbell (2022); Clinton  _et
al._ (2021), and rely for instance on Hamiltonian propagation for which the
associated quantum circuits are much deeper than those devoted to the NISQ
era.

Concerning the application of a unitary transformation onto a easy-to-prepare
known state, the unitary Van–Vleck (VV) similarity transformation, developed
in the framework of many-body perturbation theory Van Vleck (1929); Jordahl
(1934); Foldy and Wouthuysen (1950); Primas (1963); Brandow (1979); Shavitt
and Redmon (1980); Bravyi  _et al._ (2011), appears relevant to serve as a
basis for new quantum algorithms. In few words, given an Hamiltonian
$\hat{H}=\hat{H}_{0}+\hat{V}$ where $\hat{H}_{0}$ is called the unperturbed
Hamiltonian whose eigenstates are known, and $\hat{V}$ is a perturbation, the
VV similarity transformation aims to design perturbatively a unitary
transformation $\hat{U}=e^{\hat{S}^{\rm VV}}$, where $\hat{S}^{\rm VV}$ is
called the generator, that leads to an effective Hamiltonian $\bar{H}_{\rm
eff}$ in the low-energy subspace of $\hat{H}_{0}$. Ultimately, at infinite
order of perturbation, the transformed Hamiltonian
$\bar{H}=\hat{U}\hat{H}\hat{U}^{\dagger}$ is block-diagonal and is reduced to
$\bar{H}_{\rm eff}$ in the low-energy subspace of $\hat{H}_{0}$, such that the
eigenvalues of $\bar{H}_{\rm eff}$ strictly match the lowest eigenvalues of
$\hat{H}$. It follows a straightforward quantum algorithm for which the ground
state (or excited states) of a given Hamiltonian can be prepared on a quantum
computer by applying $e^{\hat{S}^{\rm VV}}$ on the known ground state (or
excited states) of the unperturbed Hamiltonian $\hat{H}_{0}$. However, an
explicit expression for $\hat{S}^{\rm VV}$ is in general unknown and
truncation of the perturbative order or approximations are mandatory.
Considering the non-interacting Hamiltonian as $\hat{H}_{0}$ and the electron-
electron Coulomb repulsion as the perturbation, the VV similarity
transformation is closely related to the unitary coupled cluster Ansatz
Shavitt and Bartlett (2009). On the other limit where $\hat{H}_{0}$ is the
Coulomb repulsion operator and $\hat{V}$ is the non-interacting Hamiltonian,
Schrieffer and Wolff (SW) derived an analytic form of $\hat{S}^{\rm VV}$, such
that $\bar{H}$ is block-diagonal at the first-order of perturbation Schrieffer
and Wolff (1966). Moreover, they showed that at the limit of small
perturbation, the Kondo model corresponds to the effective low-energy
approximation of the Anderson model. Following the work of SW, the Heisenberg
model was also shown to be the effective Hamiltonian of the Hubbard model at
half-band filling for large Coulomb repulsion strength Harris and Lange
(1967); Chao  _et al._ (1977). Yet, improvements of the SW approximation can
fairly serve as a basis to approximate $\hat{S}^{\rm VV}$ and construct an
efficient and hopefully accurate quantum algorithm for the Hubbard model. In
that context, Zhang et al. proposed two quantum algorithms devoted to finding
the VV unitary transformation in the context of spin chains. Zhang  _et al._
(2022) The first one is a quantum phase estimation based algorithm that
provides the exact transformation, but which is only realizable in the fault-
tolerant era. The second one, more adapted to the NISQ era, is an hybrid
quantum-classical algorithm based on a variational approach where the unitary
transformation (Ansatz) is built from the exponentiation of the commutator
$[\hat{H}_{0},\hat{V}]$, expressed as a linear combination of Pauli operators.

In this contribution, we derive recursive relations to the perturbative
expansion of $\bar{H}$ within the standard SW generator for the Hubbard dimer.
Following these relations, we propose two modifications of this generator, one
variational with a single parameter thanks to the recursive relations, and the
other iterative in the spirit of the Foldy–Wouthuysen transformation Foldy and
Wouthuysen (1950). Both modified SW transformations are shown to approximate,
or even perform for the homogeneous case, the desired block-diagonalization at
infinite order of perturbation, as the VV generator would provide. As a proof
of concept, we introduce two quantum algorithms associated to the modified SW
transformations on the Hubbard dimer. Finally, in light of our findings, we
discuss the perspective of generalizing our approach to larger Hubbard systems
that is left for future investigations. In particular, we show that in
contrast to most of the currently proposed Ansatz, our strategies are relevant
close to the strongly interacting regime.

##  II Van–Vleck similarity and standard Schrieffer–Wolff transformations

Let us first recall the Van–Vleck canonical perturbation theory following
Shavitt and Redmon Shavitt and Redmon (1980). Consider a Hamiltonian $\hat{H}$
with (unknown) orthonormal eigenvectors $\\{\ket{\Psi_{i}}\\}$ such that

$\hat{H}|\Psi_{i}\rangle=E_{i}|\Psi_{i}\rangle,$ (1)

and another orthonormal basis set $\\{\ket{\Phi_{i}}\\}$, eigenvectors of
another Hamiltonian $\hat{H}^{0}$ with the same dimension than $\hat{H}$, that
is related to $\\{\ket{\Psi_{i}}\\}$ by a unitary transformation,

$|\Psi_{i}\rangle=\hat{U}^{\dagger}|\Phi_{i}\rangle=\sum_{j}|\Phi_{j}\rangle\langle\Phi_{j}|\hat{U}^{\dagger}|\Phi_{i}\rangle=\sum_{j}|\Phi_{j}\rangle
U_{ij}^{\dagger}.$ (2)

The eigenvalues of $\hat{H}$ can be inferred as the elements of the diagonal
representation of the similar Hamiltonian,

$\bar{H}^{\rm VV}=\hat{U}\hat{H}\hat{U}^{\dagger},$ (3)

in the orthonormal basis $\\{\ket{\Phi_{i}}\\}$. Thus, solving the eigenvalue
problem in Eq. (1) is equivalent to searching for a unitary transformation
$\hat{U}$ such that $\bar{H}^{\rm VV}=\hat{U}\hat{H}\hat{U}^{\dagger}$ is
diagonal in a given basis set $\\{\ket{\Phi_{i}}\\}$. The reasoning remains
equivalent, though less restrictive, if solely a block-diagonalization in a
target subspace is desired. In other words, we are looking for an unknown
Hamiltonian $\bar{H}^{\rm VV}$ with eigenvectors $|\Phi_{i}\rangle$ that
shares the same eigenvalues than $\hat{H}$. If one focuses on the ground state
$\ket{\Psi_{0}}$, it is enough to only block-diagonalize $\bar{H}$,

$\displaystyle\langle\Phi_{i}|\hat{U}\hat{H}\hat{U}^{\dagger}|\Phi_{0}\rangle$
$\displaystyle=\langle\Phi_{0}|\hat{U}^{\dagger}\hat{H}\hat{U}|\Phi_{i}\rangle=0\quad\forall
i\neq 0,$ (4)
$\displaystyle\langle\Phi_{0}|\hat{U}\hat{H}\hat{U}^{\dagger}|\Phi_{0}\rangle$
$\displaystyle=E_{0}.$ (5)

Many $\hat{U}$ fulfill these conditions up to a unitary transformation acting
only on the subspace of $\\{\ket{\Phi_{i}}\\}$ with $i\neq 0$.
Let us now consider the following decomposition of the Hamiltonian,

$\hat{H}=\hat{H}^{0}+\hat{V},$ (6)

where $\hat{H}^{0}$ is diagonal in the $\\{\ket{\Phi_{i}}\\}$ basis set, i.e.
$\langle\Phi_{i}|\hat{H}^{0}|\Phi_{j}\rangle=E_{i}^{0}\delta_{ij}$. If one
wants to block-diagonalize $\bar{H}$ with respect to a given subspace
$\Omega$, for instance the one that contains all degenerate ground states of
$\hat{H}^{0}$, one can define the operator

$\hat{P}=\sum_{i\in\Omega}|\Phi_{i}\rangle\langle\Phi_{i}|$ (7)

that projects onto $\Omega$, and its complementary projector

$\hat{Q}=\hat{1}-\hat{P}=\sum_{i\notin\Omega}|\Phi_{i}\rangle\langle\Phi_{i}|.$
(8)

We note $\hat{O}_{D}=\hat{P}\hat{O}\hat{P}+\hat{Q}\hat{O}\hat{Q}$ the block-
diagonal projection of an operator $\hat{O}$ and its complementary off-block-
diagonal part $\hat{O}_{X}=\hat{P}\hat{O}\hat{Q}+\hat{Q}\hat{O}\hat{P}$.
Adopting the exponential form of the unitary transformation
$\hat{U}=e^{\hat{G}}$, $\hat{G}$ being an anti-Hermitian generator with
$\hat{G}=\hat{G}_{X}$ and $\hat{G}_{D}=0$, we seek conditions for $\hat{G}$
such that $\bar{H}^{\rm VV}$ is block-diagonal, i.e. $\bar{H}_{X}^{\rm
VV}=\hat{0}$. Within the super-operator formalism Primas (1963), $\bar{H}^{\rm
VV}$ reads:

$\displaystyle\bar{H}^{\rm VV}$
$\displaystyle=e^{\hat{G}}\hat{H}e^{-\hat{G}}=\hat{H}+[\hat{G},\hat{H}]+\frac{1}{2}[\hat{G},[\hat{G},\hat{H}]]+\dots$
$\displaystyle=\sum_{n=0}^{\infty}\frac{1}{n!}\mathcal{G}^{n}(\hat{H})=e^{\mathcal{G}}(\hat{H}),$
(9)

where $\mathcal{G}(\hat{X})=[\hat{G},\hat{X}]$. By decomposing
$e^{\mathcal{G}}(\hat{H})={\rm
cosh}\mathcal{G}(\hat{H})+\rm{sinh}\mathcal{G}(\hat{H})$, it follows that the
condition $\bar{H}_{X}^{\rm VV}=\hat{0}$ is fulfilled for

$[\hat{G},\hat{H}^{0}]=-[\hat{G},\hat{V}_{D}]-\sum_{n=0}^{\infty}c_{n}\mathcal{G}^{2n}(\hat{V}_{X}),$
(10)

where $c_{n}=2^{2n}B_{2n}/(2n)!$ are functions of Bernoulli numbers $B_{2n}$.
Eq. (10) is the central equation of the VV canonical perturbation theory that
defines the generator $\hat{G}$ such that $\bar{H}^{\rm VV}$ is block-
diagonal, thus expressing the eigenstates of $\hat{H}$ in terms of eigenstates
of $\hat{H}^{0}$ through $\hat{G}$. Using an order by order expansion of
$\hat{G}$, i.e. $\hat{G}=\sum_{n=1}\hat{G}^{(n)}$, conditions to cancel
$\bar{H}_{X}$ can be obtained at each order as,

$\displaystyle[\hat{G}^{(1)},\hat{H}^{0}]=-\hat{V}_{X},$ (11)
$\displaystyle[\hat{G}^{(2)},\hat{H}^{0}]=-[\hat{G}^{(1)},\hat{V}_{D}],$ (12)
$\displaystyle[\hat{G}^{(3)},\hat{H}^{0}]=-[\hat{G}^{(2)},\hat{V}_{D}]-\frac{1}{3}[\hat{G}^{(1)},[\hat{G}^{(1)},\hat{V}_{X}]],$
(13) $\displaystyle\dots$

It follows that $\bar{H}^{\rm VV}$ can also be expressed order by order as

$\bar{H}^{\rm
VV}=\hat{H}^{0}+\hat{V}_{D}+\sum_{n=0}^{\infty}t_{n}\mathcal{G}^{2n+1}(\hat{V}_{X}),$
(14)

with $t_{n}=2(2^{2n+2}-1)B_{2n+2}/(2n+2)!$. As mentioned in Ref. [Shavitt and
Redmon, 1980], the Van–Vleck perturbation theory equations (10) and (14) are
expressed in the domain of a Lie algebra, thus allowing an equivalent
diagrammatic expansion. Note that the convergence of perturbative series and
the diagrammatic expansion has been thoroughly investigated in Ref. [Bravyi
_et al._ , 2011], which also provides recursive relations to obtain the $n$-th
order term $\hat{G}^{(n)}$ of the VV generator $\hat{G}$ as a function of the
previous $n-1$ terms.

In practice, finding both an analytic and a numerical form of $\hat{G}$ for a
given $\hat{H}^{0}$ and $\hat{V}$ remains challenging, at least equivalent as
the explicit diagonalisation of $\hat{H}$. From the perspective of developing
quantum algorithms based on the VV formalism, one realizes that the number of
terms in the generator drastically increases order by order, thus leading to
deeper circuits and, consequently, to an increase in complexity and
sensibility to noise of quantum algorithms. To overcome this issue, we explore
an alternative approach which consists in using a truncated generator, the
Schrieffer–Wolff generator, that is later modified by adding a variational
parameter or by using an iterative process to compensate the resulting
truncation error.

First of all, following Ref. [Schrieffer and Wolff, 1966], let us recall the
Schrieffer–Wolff transformation in the context of the half-filled Hubbard
model that we decompose as in Eq. (6) into a local part,

$\displaystyle\hat{H}^{0}$ $\displaystyle=$
$\displaystyle\sum_{i\sigma}\mu_{i}\hat{n}_{i\sigma}+\sum_{i}U_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow},$
(15)

and a non-local (kinetic) part,

$\displaystyle\hat{V}$ $\displaystyle=$ $\displaystyle-\dfrac{1}{2}\sum_{i\neq
j,\sigma}t_{ij}\left(\hat{\gamma}_{ij\sigma}+\hat{\gamma}_{ji\sigma}\right),$
(16)

with $\hat{n}_{i\sigma}=\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma}$ and
$\hat{\gamma}_{ij\sigma}=\hat{c}_{i\sigma}^{\dagger}\hat{c}_{j\sigma}$, and
$\hat{c}^{\dagger}_{i\sigma}$ ($\hat{c}_{i\sigma}$) the creation
(annihilation) operator of an electron of spin
$\sigma=\\{\uparrow,\downarrow\\}$ in site $i$. This decomposition contrasts
with the usual decomposition between the non-interacting part for which the
solution is easily accessible and the non-trivial canonical (interacting)
part. Indeed, the ground state of $\hat{H}^{0}$ is degenerate at half filling
for $U>0$, and consists in a superposition of all states having no double
occupation (spanning the so-called Heisenberg subspace in this paper).
Starting from the atomic limit ($U/t\rightarrow\infty$), Schrieffer and Wolff
have proposed, in the original context of an Anderson Hamiltonian, to use the
unitary transformation $\hat{U}=e^{\hat{S}}$ such that

$\bar{H}^{\rm SW}=e^{\hat{S}}He^{-\hat{S}},$ (17)

that we denote simply $\bar{H}$ in the following to simplify notations, is
block-diagonalized at first order of perturbation, i.e.

$[\hat{S},\hat{H}^{0}]=-\hat{V}.$ (18)

We highlight that the above equation corresponds to the first order of
perturbation of the VV relations, Eq. (11), i.e that the SW generator
$\hat{S}$ block-diagonalizes the Hamiltonian only at first order, contrary to
the VV generator $\hat{G}$.
Note that $\hat{V}_{X}=\hat{V}$ when the operator $\hat{P}$ projects onto the
Heisenberg subspace. It can be shown that under the SW condition (18),
$\hat{S}$ takes the following form,

$\hat{S}=\dfrac{1}{2}\sum_{i\neq
j,\sigma}\hat{p}_{ij\bar{\sigma}}\left(\hat{\gamma}_{ij\sigma}-\hat{\gamma}_{ji\sigma}\right),$
(19)

with $\hat{p}_{ij\sigma}$ defined as

$\displaystyle\hat{p}_{ij\sigma}=\sum_{x=0}^{3}\lambda_{ij\sigma,x}\,\hat{p}_{ij\sigma,x},$
(20)

where
$\hat{p}_{ij\sigma,0}=\left(1-\hat{n}_{i\sigma}\right)\left(1-\hat{n}_{j\sigma}\right)$,
$\hat{p}_{ij\sigma,1}=\hat{n}_{i\sigma}\left(1-\hat{n}_{j\sigma}\right)$,
$\hat{p}_{ij\sigma,2}=\left(1-\hat{n}_{i\sigma}\right)\hat{n}_{j\sigma}$,
$\hat{p}_{ij\sigma,3}=\hat{n}_{i\sigma}\hat{n}_{j\sigma}$ are projectors, i.e.
$\sum_{x=0}^{3}\hat{p}_{ij\sigma,x}=\hat{1}$, and

$\displaystyle\lambda_{ij\sigma,0}$
$\displaystyle=-\frac{t_{ij}}{\Delta\mu_{ij}}\;\;{\rm if}\;\Delta\mu_{ij}\neq
0;\;\;\lambda_{ij\sigma,0}=0\;\;{\rm else},$ (21)
$\displaystyle\lambda_{ij\sigma,1}$
$\displaystyle=-\frac{t_{ij}}{\Delta\mu_{ij}+U_{i}}\;\;{\rm
if}\;\Delta\mu_{ij}+U_{i}\neq 0;\;\;\lambda_{ij\sigma,1}=0\;\;{\rm else},$
(22) $\displaystyle\lambda_{ij\sigma,2}$
$\displaystyle=-\frac{t_{ij}}{\Delta\mu_{ij}-U_{j}}\;\;{\rm
if}\;\Delta\mu_{ij}-U_{j}\neq 0;\;\;\lambda_{ij\sigma,2}=0\;\;{\rm else},$
(23) $\displaystyle\lambda_{ij\sigma,3}$
$\displaystyle=-\frac{t_{ij}}{\Delta\mu_{ij}+\Delta U_{ij}}\;\;{\rm
if}\;\Delta\mu_{ij}+\Delta U_{ij}\neq 0;\;\;\lambda_{ij\sigma,3}=0\;\;{\rm
else},$ (24)

with $\Delta\mu_{ij}=\mu_{i}-\mu_{j}$ and $\Delta U_{ij}=U_{i}-U_{j}$. Within
the SW transformation, we obtain

$\bar{H}=\hat{H}^{0}+\sum_{n=2}^{\infty}\frac{n-1}{n!}\mathcal{S}^{n-1}(\hat{V})=\hat{H}^{0}+\sum_{n=2}^{\infty}\bar{H}^{(n)},$
(25)

where $\mathcal{S}$ is the super-operator defined as
$\mathcal{S}(\hat{X})=[\hat{S},\hat{X}]$, and
$\bar{H}^{(n)}=\dfrac{n-1}{n!}\mathcal{S}^{n-1}(\hat{V})$.

Since $\hat{S}$ consists in a truncated form of $\hat{G}$, $\bar{H}$ is not
expected to be block-diagonal anymore, i.e. $\bar{H}_{X}\neq\hat{0}$. In the
following, we propose recursive relations between $\mathcal{S}^{n}(\hat{V})$
and $\mathcal{S}^{n-1}(\hat{V})$, derived for the Hubbard dimer, that provide
an explicit expression for $\bar{H}$ and in particular for $\bar{H}_{X}$ in
terms of two- and three-body operators. These relations are further exploited
to develop two modifications of $\hat{S}$, one variational and the other
iterative, designed to minimize or even cancel $\bar{H}_{X}$ while conserving
the same complexity as $\hat{S}$.

##  III Modified Schrieffer–Wolff transformations

###  III.1 Recursive relations

We establish recursive relations for each order of Eq. (25), which details are
provided in Appendix A. More precisely, we find that even orders are block-
diagonal, i.e. $\bar{H}^{(2n)}_{D}=\bar{H}^{(2n)}$ and reads

$\displaystyle\bar{H}^{(2n)}$ $\displaystyle=$
$\displaystyle\dfrac{2n-1}{2(2n)!}\sum_{\begin{subarray}{c}i\neq j\\\
\sigma\end{subarray}}\sum_{x=0}^{3}K^{(2n-1)}_{ij\sigma,x}\hat{p}_{ij\bar{\sigma},x}\left(\hat{n}_{i\sigma}-\hat{n}_{j\sigma}\right)$
(26) $\displaystyle+\dfrac{2n-1}{2(2n)!}\sum_{\begin{subarray}{c}i\neq j\\\
\sigma\end{subarray}}J^{(2n-1)}_{ij\sigma}\left(\hat{\gamma}_{ij\sigma}\hat{\gamma}_{ji\bar{\sigma}}+\hat{\gamma}_{ji\sigma}\hat{\gamma}_{ij\bar{\sigma}}\right)$
$\displaystyle+\dfrac{2n-1}{2(2n)!}\sum_{\begin{subarray}{c}i\neq j\\\
\sigma\end{subarray}}L^{(2n-1)}_{ij\sigma}\left(\hat{\gamma}_{ij\sigma}\hat{\gamma}_{ij\bar{\sigma}}+\hat{\gamma}_{ji\sigma}\hat{\gamma}_{ji\bar{\sigma}}\right),$

while odd orders are found to be off-block-diagonal, i.e.
$\bar{H}^{(2n+1)}_{X}=\bar{H}^{(2n+1)}$ and take the following expression,

$\bar{H}^{(2n+1)}=\dfrac{2n}{2(2n+1)!}\sum_{i\neq
j\sigma}\sum_{x=0}^{3}T_{ij\sigma,x}^{(2n)}\hat{p}_{ij\bar{\sigma},x}\left(\hat{\gamma}_{ij\sigma}+\hat{\gamma}_{ji\sigma}\right).$
(27)

and where only the expression of interaction integrals $I_{ij\sigma,x}^{(k)}$
($I=J,K,L$ or $T$) depend on the order $k$. Explicit formulas for the
interaction integrals $I_{ij\sigma,x}^{(k)}$ are given in Appendix A. By
summing over all orders, $\bar{H}$ takes exactly the following form,

$\displaystyle\bar{H}$
$\displaystyle=H^{0}+\sum_{n=2}^{\infty}\bar{H}^{(n)}=\bar{H}_{D}+\bar{H}_{X},$
(28)

with $\bar{H}_{D}=\hat{H}^{0}+\sum_{n=1}^{\infty}\bar{H}^{(2n)}$ and
$\bar{H}_{X}=\sum_{n=1}^{\infty}\bar{H}^{(2n+1)}$. Specifically, the non
block-diagonal contribution $\bar{H}_{X}$ couples states from the Heisenberg
subspace to states from its complementary subspace, and reads explicitely

$\displaystyle\bar{H}_{X}=\dfrac{1}{2}\sum_{\begin{subarray}{c}i\neq j\\\
\sigma\end{subarray}}\sum_{x=0}^{3}T_{ij\sigma,x}\hat{p}_{ij\bar{\sigma},x}\left(\hat{\gamma}_{ij\sigma}+\hat{\gamma}_{ji\sigma}\right),$
(29)

where the interaction integrals $T_{ij\sigma,x}$ are obtained by summing over
all orders as

$T_{ij\sigma,x}=\sum_{n=1}^{\infty}\dfrac{2n}{(2n+1)!}T_{ij\sigma,x}^{(2n)}.$
(30)

For the homogeneous case, the associated integrals are simply given by

$\displaystyle T_{ij\sigma,0}=T_{ij\sigma,3}=0,$ (31)

and

$\displaystyle T_{ij\sigma,1}=T_{ij\sigma,2}=-t_{ij}\left(\cos{(4t/U)}-{\rm
sinc}\,(4t/U)\right).$ (32)

Explicit form of the block-diagonal contributions and relations for the
inhomogeneous Hubbard dimer are derived in Appendix A.

Figure 1: Schematic representation of the action of the different operators
$\bar{H}$ in the Hilbert space of the half-filled Hubbard dimer.

At this stage, we have established recursive relations to obtain the similar
Hamiltonian $\bar{H}$ within the standard SW transformation (SWT) at infinite
order of perturbation. However, given the definition of $\hat{S}$ in Eq. (18),
the standard SW transformation at infinite order does not lead to a block-
diagonal representation of $\bar{H}$ with respect to the Heisenberg subspace
since $\bar{H}_{X}\neq\hat{0}$, see Fig 1. Based on the previous recursive
relations, in the following subsections we present two strategies, denoted as
modified SW (MSW) transformations, one variational and the other iterative, to
fully perform the desired block-diagonalization.

###  III.2 Variational Schrieffer–Wolff transformation

We propose to introduce a variational scaling parameter $\theta$ to the
unitary transformation,

$\hat{U}(\theta)=e^{\theta\hat{S}},$ (33)

such that for $\theta=0$, $\hat{U}=\hat{1}$, and for $\theta=1$ one recovers
the standard SW transformation $\hat{U}=e^{\hat{S}}$. Within this unitary
transformation, the similar Hamiltonian $\bar{H}(\theta)$ reads

$\displaystyle\bar{H}(\theta)$
$\displaystyle=e^{\theta\hat{S}}\hat{H}e^{-\theta\hat{S}},$
$\displaystyle=\hat{H}^{0}+\hat{V}\left(1-\theta\right)+\sum_{n=2}^{\infty}\frac{\theta^{n-1}(n-\theta)}{n!}\mathcal{S}^{n-1}(\hat{V}).$
(34)

Using the previously established recursive relations and after summation till
the infinite order, see Appendix B, it can be decomposed as follows, similarly
as in Eq. (28),

$\displaystyle\bar{H}(\theta)=\bar{H}_{D}(\theta)+\bar{H}_{X}(\theta),$ (35)

where the $\theta$-dependence lies in the renormalized interaction integrals
that read for the non block-diagonal contribution in the homogeneous case,

$\displaystyle T_{ij\sigma,0}(\theta)=T_{ij\sigma,3}(\theta)=0,$ (36)

and

$\displaystyle
T_{ij\sigma,1}(\theta)=T_{ij\sigma,2}(\theta)=-t\left(\cos{(4t\theta/U)}-\theta{\rm
sinc}\,(4t\theta/U)\right),$ (37)

The scaling parameter $\theta$ can be optimized to minimize the contributions
from the coupling operator $\bar{H}_{X}(\theta)$, which is shown to cancel out
for the homogeneous case at

$\displaystyle\theta=(U/4t)\tan^{-1}{(4t/U)},$ (38)

leading to an exact block-diagonalization of $\bar{H}(\theta)$. More
precisely, we demonstrate in Appendix B that in this case and at the saddle
point, the variational generator fulfill the VV condition in Eq. (10).
Relations become more complex for the inhomogeneous cases and the variational
process has to be done numerically, see Appendix B. In this case, the
cancellation of $\bar{H}_{X}(\theta)$ cannot always be reached. Alternatively,
one can minimize the energy of $\bar{H}(\theta)$ restricted to the Heisenberg
subspace, which is equivalent to maximize the overlap between the minimizing
state and the exact ground state $\ket{\Psi_{0}}$. The difference between the
two optimization schemes is discussed in Appendix D.

###  III.3 Iterative Schrieffer–Wolff transformation

Alternatively to the variational approach, one can take advantage of the
similarity between the coupling terms in Eq. (76) and the perturbation
$\hat{V}$ in Eq. (50). In the spirit of the Foldy–Wouthuysen transformation
Foldy and Wouthuysen (1950), we propose the following iterative scheme:

  1. 1.

Initialize the iterative process by applying the standard SW transformation on
the original problem to obtain $\bar{H}^{(s=0)}$.

  2. 2.

At the iteration $s=s+1$, define the new problem
$\bar{H}^{0(s)}=\bar{H}^{(s)}-\bar{V}^{(s)}$ and
$\bar{V}^{(s)}=\bar{H}_{X}^{(s-1)}$.

  3. 3.

Find the corresponding generator $\hat{S}^{(s)}$ such that
$[\bar{H}^{0(s)},\hat{S}^{(s)}]=\bar{V}^{(s)}$.

  4. 4.

Use the recursive relations derived in Appendix C to obtain the new
$\bar{H}^{(s)}$.

  5. 5.

Repeat steps 2 to 4 until convergence is reached, i.e.
$\bar{H}_{X}^{(s)}\rightarrow\hat{0}$.

After $N_{s}$ iterations, the iterative unitary transformation and the similar
Hamiltonian are given by

$\hat{U}^{(N_{s})}=\left(\prod_{s=0}^{N_{s}-1}e^{\hat{S}^{(s)}}\right),$ (39)

and

$\bar{H}^{(N_{s})}=\hat{U}^{(N_{s})}\hat{H}\hat{U}^{(N_{s})^{\dagger}},$ (40)

respectively. The amplitudes of the resulting coupling terms for large $U/t$
behave asymptotically in $(t^{2}/U)^{N_{s}}$ for $N_{s}$ iterations, such that
the iterative algorithm converges exponentially to a precise decoupling.

###  III.4 Perspectives for larger Hubbard rings

The iterative and variational MSW transformations are shown to perform exact
(quasi) block-diagonalization for the homogeneous (inhomogeneous) Hubbard
dimer, respectively, thanks to the recursive properties in Eqs. (26) and (27).
Before investigating the quantum algorithms associated to the presented MSW
transformation applied to the Hubbard dimer, we discuss possible extensions to
larger systems. First of all, note that the recursive relations obtained in
Eqs. (26) and (27) are not valid for larger systems, where additional terms
giving rise to interactions among more than two sites emerge. Nonetheless, for
the purposes of this study, we neglect these terms, meaning that VV
perturbation condition in Eq. (10) is not satisfied, but the SW condition in
Eq. (18) (i.e first order) still is. In this section, the truncation error is
assessed on a classical computer for homogeneous nearest neighbour (NN)
Hubbard rings of up to $N=10$ sites. To do so we apply the unitary
transformation $\hat{U}^{\dagger}(\theta)=e^{-\theta\hat{S}}$ to the ground
state $\ket{\Phi_{\rm Heis}}$ of the NN antiferromagnetic Heisenberg model
$J\sum_{ij}\hat{s}_{i}\hat{s}_{j}$, where $\hat{s}_{i}$ denotes the spin
operator on site $i$ and $J>0$ is the spin-coupling element, which corresponds
to the strongly correlated limit of the NN Hubbard model Schrieffer and Wolff
(1966). We follow the variational scheme presented in Sec. III.2, where
$\theta$ is optimized to minimize the expectation value
$E(\theta)=\langle\Phi_{\rm Heis}|\bar{H}(\theta)|\Phi_{\rm Heis}\rangle$ with
$\bar{H}(\theta)=U(\theta)\hat{H}U^{\dagger}(\theta)$, which is equivalent to
maximize the overlap of $U^{\dagger}(\theta)|\Phi_{\rm Heis}\rangle$ with the
exact ground state $|\Psi_{0}\rangle$.

Figure 2: Relative errors in the ground-state energy calculated for
homogeneous half-filled Hubbard rings with respect to the number of sites.
Results are given for the variational MSW transformation (solid lines) and the
standard SWT ($\theta=1$, dashed lines). Lines are guide for the eye.

In Fig. 2 we show, as a function of the number of sites $N$ and for different
values of the Coulomb repulsion strenght $U/t$, the relative error (in $\%$)
of $E(\theta)$ with respect to the exact ground-state energy. Results are
provided for $\theta=1$, corresponding to the standard SWT, and for the
optimal value $\theta^{*}$. As $N$ increases the relative error increases and
appears to converge to what would correspond to the truncation error. As
expected, the truncation error increases as $U/t$ decreases, i.e. $\sim
1\%\,(1\%),\sim 5\%\,(7\%)$ and $\sim 11\%\,(30\%)$ for $U/t=20$, 8, 4 and
$\theta=\theta^{*}$ ($\theta=1)$, respectively. The introduction of a single
and variational parameter systematically and drastically improves over the
standard SWT. Consequently, the variational extension to the SW approximation,
that is exact for the homogeneous half-filled Hubbard dimer, remains a good
approximation for larger system sizes in the intermediate to strongly
correlated regime. Straightforward improvements can be envisioned by
considering higher-order contributions from the generator, following the
recursive relations proposed in Ref. [Bravyi  _et al._ , 2011], for instance.
However, the implementation of the MSW transformations on classical computers
is computationally intractable for systems beyond $\sim 16$ orbitals, in
analogy with the unitary coupled cluster ansatz Romero  _et al._ (2018). This
is also the case for the construction of the trial state $|\Phi_{\rm
Heis}\rangle$ for large system’s size, which we disregard in the following by
investigating quantum algorithms applied to the Hubbard dimer, for which
$|\Phi_{\rm Heis}\rangle$ is easy to prepare.

##  IV Modified SW transformations applied on quantum computers

At this stage, we investigate the relevance of the aforementioned MSW
transformations, $\bar{H}^{\rm MSW}=\bar{H}(\theta)$ and $\bar{H}^{\rm
MSW}=\bar{H}^{(N_{s})}$ in Secs. III.2 and III.3, respectively, for the design
of new quantum algorithms.

In both cases, $\bar{H}^{\rm MSW}$ is block-diagonal for the homogeneous case,
so that ground or excited states can easily be constructed as linear
combination of two basis vectors for the Hubbard dimer, see Fig. 1. For the
homogeneous half-filled Hubbard dimer, relevant trial eigenstates of
$\bar{H}^{\rm MSW}$ consist in the Heisenberg state $|\Phi_{\rm
Heis}\rangle=(1/\sqrt{2})\left(|\uparrow\;\downarrow\,\rangle+|\downarrow\;\uparrow\,\rangle\right)$
and the ionic state $|\Phi_{\rm
Ionic}^{\alpha}\rangle=\cos(\alpha)|\uparrow\downarrow\;\cdot\,\rangle+\sin(\alpha)|\,\cdot\;\uparrow\downarrow\,\rangle$.
Indeed, $e^{\hat{S}}$ preserving the spin symmetry, triplet states
$|\uparrow\;\uparrow\,\rangle$ and $|\downarrow\;\downarrow\,\rangle$ are
discarded. The eigenstates of $\hat{H}$ can then be constructed from the trial
eigenstates of $\bar{H}^{\rm MSW}$ by applying the transformation
$\hat{U}^{\rm MSW}$, which refers to the variational [see Eq. (33)] or to the
iterative [see Eq. (39)] MSW transformation. It appears clear that both the
variational or iterative MSW approaches are adapted to the design of quantum
algorithms, as they are both formulated as a unitary transformation applied to
an easy-to-prepare initial state.

As a proof of concept, we have implemented both quantum algorithms to treat
the homogeneous and inhomogeneous Hubbard dimer, using Qiskit tA v  _et al._
(2021) to construct the quantum circuits. We use the one-to-one correspondence
between the states of the qubits and the occupation of the spin-orbitals of
the Hubbard dimer to map our states onto qubits, with even-numbered qubits
corresponding to spin-up orbitals and odd-numbered qubits to spin-down
orbitals. The fermionic creation and annihilation operators are mapped onto
Pauli strings $\hat{\mathcal{P}}_{i}$ using the Jordan–Wigner (JW)
transformation Jordan and Wigner (1928). To implement the unitary
transformation on quantum circuits, the first-order Trotter–Suzuki
approximation is used, i.e. the exponential of the sum of Pauli strings is
decomposed into a product of exponential of a single Pauli string,

$\displaystyle e^{\theta\hat{S}}\xrightarrow{\rm
JW}e^{\theta\sum_{i}\xi_{i}\hat{\mathcal{P}}_{i}}\approx\prod_{i}e^{\theta\xi_{i}\hat{\mathcal{P}}_{i}},$
(41)

for which the associated circuit is known (see panel (c) of Fig. 3), and
$\\{\xi_{i}\\}$ are the coefficients that are functions of the SW generator
parameters $\\{\lambda\\}$ obtained after the JW transformation. The trial
Heisenberg and ionic states can be easily prepared on the quantum computer, as
shown in panels (a) and (b) of Fig. 3. Finally, we simulate our variational
MSW transformation using a noise model built on Qiskit. This noise model
consists in a depolarizing quantum error channel applied on every one- and
two-qubit gates, with depolarizing error parameters of $\lambda_{1}=0.0001$
and $\lambda_{2}=0.001$, respectively. Note that the 4-qubit circuit resulting
from the variational MSW transformation is composed of 32 one-qubit gates and
35 CNOT gates, and that no readout error is considered. Sampling noise is also
added to this noise model by considering $n_{\rm shots}=8192$ for the
estimation of the expectation value of each Pauli string. The variational
parameter was optimized by using the SPSA optimizer with a maximum of 1000
iterations. Then, the optimal parameter $\theta^{*}$ is calculated as the mean
of the last 25 iterations, and the expectation values of $\bar{H}(\theta^{*})$
with respect to the Heisenberg and ionic states are estimated as the mean of
another 100 noisy simulations (with fixed parameter $\theta^{*}$). The noisy
results are then compared to the exact references obtained by exact
diagonalization, as well as to the noiseless state-vector simulation, without
considering any quantum or sampling noise and for which the L-BFGS-B optimizer
was used to update the variational parameter. For the iterative MSW
transformation, only state-vector simulations is performed.

Figure 3: a) Quantum circuit corresponding to the Heisenberg state
$\ket{\Phi_{\rm
Heis}}=\left(\ket{\,\uparrow\;\downarrow\,}-\ket{\,\downarrow\;\uparrow\,}\right)/\sqrt{2}=\left(\ket{1001}-\ket{0110}\right)/\sqrt{2}$,
b) Quantum circuit corresponding to the linear combination of the ionic states
$\ket{\Phi_{\rm
Ionic}^{\alpha}}=\cos(\alpha)\ket{\,\uparrow\downarrow\;\cdot\,}+\sin(\alpha)\ket{\,\cdot\;\uparrow\downarrow\,}=\cos(\alpha)\ket{1100}+\sin(\alpha)\ket{0011}$,
c) Quantum circuit corresponding to the implementation of $e^{\xi
X_{0}Z_{1}Y_{2}Z_{3}}$.

###  IV.1 Variational approach

The variational approach described in Sec. III.2 consists in finding the
optimal parameter $\theta_{X}$ of the unitary in Eq. (33) such that the
couplings $\bar{H}_{X}(\theta_{X})$ are minimized, thus enforcing the block-
diagonalization of $\bar{H}(\theta_{X})$. Compared to the strategy of Zhang
and coworkers Zhang  _et al._ (2022), our minimization process implies only a
single variational parameter, rather than a number of parameters that would
correspond to the number of Pauli strings composing the generator [i.e. each
$\xi_{i}$ in Eq. (41)]. Minimizing the couplings $\bar{H}_{X}(\theta_{X})$
requires the estimation of the off-diagonal matrix elements of $\bar{H}$ such
as done in Ref. [Zhang  _et al._ , 2022]. However, this is not straightforward
on quantum computers in contrast to the measurement of expectation values,
although one can note some improvements in the literature Huggins  _et al._
(2020); Stair and Evangelista (2021). Consequently in this work, in analogy
with the VQE algorithm, we minimize the energy
$\langle\Phi\hat{U}(\theta)|\hat{H}|\hat{U}^{\dagger}(\theta)\Phi\rangle$
measured on the quantum device, rather than minimizing the norm of
$\bar{H}_{X}$. Note that the two strategies are equivalent when the initial
trial state $|\Phi\rangle$ is indeed the ground state of $\bar{H}(\theta)$.
Otherwise, it does not lead to the expected block-diagonalization of the
Hamiltonian, as discussed in more details in Appendix D.

Figure 4: Energies of the half-filled Hubbard dimer for $\Delta\mu=0$ with
respect to the repulsion strength, using the variational SW method to minimize
the energy $\bra{\Phi_{\rm Heis}}\bar{H}(\theta)\ket{\Phi_{\rm Heis}}$.

Let us start with the homogeneous Hubbard dimer ($\Delta\mu=0$). In this case,
analytical expressions for the variational SW transformation can be derived as
shown in Sec. III.2. As readily seen in Fig. 4, the energy obtained by
minimizing $\bra{\Phi_{\rm Heis}}\bar{H}(\theta)\ket{\Phi_{\rm Heis}}$ matches
exactly the ground-state energy of $\hat{H}$, and the minimizing parameter,
denoted by $\theta_{\rm Heis}$, is exactly the same as the analytical
expression in Eq. (38) (not shown). In addition, using the exact same unitary
$\hat{U}(\theta_{\rm Heis})$ but on the equi-weighted ionic state $|\Phi_{\rm
Ionic}^{\alpha=\pi/4}\rangle$, one recovers the first-excited singlet energy
of $\hat{H}$. Thus, our variationally optimized SW transformation has indeed
block-diagonalized $\hat{H}$ exactly for any repulsion strength $U/t$, with
the Heisenberg subspace containing the singlet ground state and the triplet
states. Looking at the energies obtained from the noisy simulation, they
follow closely the noiseless results, especially for the first-excited state
energy for which the relative error does not exceed 1.5%. We note also an
increase of around 0.03 in the expectation value of the spin operator
$\hat{S}^{2}$ due to the noise, showing that the final state is not a pure
singlet state anymore.

Figure 5: Energies of the half-filled Hubbard dimer for $\Delta\mu/t=2$ (top
panel) with respect to the repulsion strength, using the variational SW method
to minimize the energies $\bra{\Phi_{\rm Heis}}\bar{H}(\theta)\ket{\Phi_{\rm
Heis}}$ (orange markers, shown for $U>\Delta\mu$) and $\bra{\Phi_{\rm
Ionic}^{\alpha=0}}\bar{H}(\theta)\ket{\Phi_{\rm Ionic}^{\alpha=0}}$ (blue
markers, shown for $U<\Delta\mu$). The associated minimizing parameters
$\theta_{\rm Heis}$ and $\theta_{\rm Ionic}$ are shown on the bottom panel,
respectively. The vertical dotted line corresponds to $U=\Delta\mu$.

Turning to the inhomogeneous Hubbard dimer with $\Delta\mu/t=2$, no analytical
expressions are known for the optimal parameter $\theta$. In contrast to the
homogeneous case, minimizing the energy $\bra{\Phi_{\rm
Heis}}\bar{H}(\theta)\ket{\Phi_{\rm Heis}}$ doesn’t lead to a block-diagonal
$\bar{H}(\theta_{\rm Heis})$ in the entire range of interaction, but only for
$U\gg\Delta\mu$ as shown in Fig. 5. In the other case, the ground state
doesn’t belong to the Heisenberg subspace such that $\bar{H}(\theta_{\rm
Heis})$ is not block-diagonal. Hence, the Heisenberg state is not an
eigenstate of $\bar{H}(\theta_{\rm Heis})$, neither is the ionic state (see
Appendix D for more details). However, rather than minimizing the energy with
respect to the Heisenberg state, one can prepare a different initial trial
state corresponding to the ground state (or a good approximation of it) that
belongs to the other subspace. In the case of the Hubbard dimer, this is the
ionic subspace which ground state is a linear combination of the ionic states
(see Fig. 1 and panel (b) of Fig. 3). As $\Delta\mu/t=2$, the optimal $\alpha$
value of $\ket{\Phi_{\rm Ionic}^{\alpha}}$ is not trivial and we approximate
it as 0, i.e. $|\Phi_{\rm
Ionic}^{\alpha=0}\rangle=\ket{\,\uparrow\downarrow\;\cdot\,}$. Minimizing
$\langle\Phi_{\rm Ionic}^{\alpha=0}|\bar{H}(\theta)|\Phi_{\rm
Ionic}^{\alpha=0}\rangle$ now leads to an optimal $\theta_{\rm Ionic}$ that
approximately block-diagonalizes $\bar{H}(\theta_{\rm Ionic})$. More
precisely, one recovers the correct ground-state and first-excited-state
singlet energies for $U\ll\Delta\mu$ by measuring the expectation values
$\langle\Phi_{\rm Ionic}^{\alpha=0}|\bar{H}(\theta_{\rm Ionic})|\Phi_{\rm
Ionic}^{\alpha=0}\rangle$ and $\bra{\Phi_{\rm Heis}}\bar{H}(\theta_{\rm
Ionic})\ket{\Phi_{\rm Heis}}$, where $\theta_{\rm Ionic}$ is defined as the
optimal parameter that minimizes $\langle\Phi_{\rm
Ionic}^{\alpha=0}|\bar{H}(\theta)|\Phi_{\rm Ionic}^{\alpha=0}\rangle$.
Interestingly, the Heisenberg state now belongs to the subspace which contains
the first-excited singlet state, as opposed to the correlation regime
$U\gg\Delta\mu$. Note that in the strictly correlated (or atomic) limit
$U/t\rightarrow\infty$, $\theta_{\rm Heis}$ tends to 1 (see bottom panel of
Fig. 5), which is expected as the variational SW transformation tends to the
standard SW transformation that is exact in this limit. Moving from this
limit, the value of the optimal parameter $\theta_{\rm Heis}$ decreases to
compensate the error from applying the MSW transformation in the non-atomic
limit. Finally, the noisy simulations show a relatively good agreement with
the noiseless results. In analogy with the homogeneous model on Fig. 4, the
expectation value of $\hat{S}^{2}$ also increases from 0 to around 0.05, and
the deviation in energy is more significant on the ground-state energy and
when $U/t$ increases. According to the bottom panel of Fig. 5, it seems that
the optimized parameter obtained from the noisy simulation deviates
significantly from the exact one for large $U/t$ values (last blue circle on
the curve), thus indicating that the classical optimization for large $U/t$
values is more challenging in the noisy environment. This could be mitigated
by employing error mitigation strategies that are outside of the scope of this
manuscript Cai  _et al._ (2022).

###  IV.2 Iterative approach

Figure 6: Energies of the half-filled Hubbard dimer for $\Delta\mu=0$ with
respect to the repulsion strength, using the iterative SW transformation
applied on the Heisenberg state (orange triangles) and the equi-weighted ionic
state (blue crosses). The exchange integrals $J_{01\sigma}^{(N_{s})}$ obtained
at convergence are also represented.

The steps 1 to 5 of the iterative approach described in Sec. III.3 can all be
performed on a classical computer by using the recursive relations derived in
Appendix B, such that only the preparation of the final state
$U^{(N_{s})\dagger}\ket{\Phi}$ and the measurement of $\hat{H}$ are done on
the quantum device. Let us start with the homogeneous dimer in Fig. 6.
Interestingly, and in contrast to the variational approach, the ground state
doesn’t always belong to the Heisenberg subspace. Indeed, the ground and
second-excited singlet states of $\bar{H}^{(N_{s})}$ oscillate between the
Heisenberg and the equi-weighted ($\alpha=\pi/4$) ionic states. This can be
rationalized by analyzing the behaviour of the exchange integrals. Indeed, the
analytical function in Eq. (A) at iteration 0 (corresponding to the standard
SW transformation) shows that the exchange integrals oscillate and change sign
for different correlation strength (not shown). The iterative process strongly
sharpens these oscillations (though the function remains continuous and
infinitely differentiable for all $U/t>0$), as shown by the solid blue lines
in Fig. 6. The change of sign of the exchange integrals indicates a change in
the ground state of $\bar{H}^{(N_{s})}$. If they are negative, the ground
state belongs to the Heisenberg subspace, while it belongs to the ionic
subspace if they are positive. One can also verify that the energy required to
go from the Heisenberg state to the equi-weighted ionic state is of
$4J_{01\sigma}^{(N_{s})}$, where $J_{01\sigma}^{(N_{s})}$ are the couplings
terms obtained after $N_{s}$ iterations. Finally, note that the first-excited
singlet state energy of $\hat{H}$ is actually exactly recovered from the
$|\Phi_{\rm Ionic}^{\alpha=-\pi/4}\rangle$ state that is an eigenstate of
$\bar{H}^{(N_{s})}$ (not shown).

Figure 7: Energies of the half-filled Hubbard dimer for $\Delta\mu/t=2$ with
respect to the repulsion strength. The iterative SW transformation is applied
on the Heisenberg state (orange triangles) and on the pure ionic state (blue
crosses), with (top panel) and without (bottom panel) trotterization error.
The vertical dotted line corresponds to $U=\Delta\mu$.

Turning to the inhomogeneous dimer with $\Delta\mu/t=2$ in Fig. 7, one
observes a similar behaviour than for the variational approach in Fig. 5, i.e.
the ground state belongs to the Heisenberg subspace for $U\gg\Delta\mu$ and to
the ionic subspace for $U\ll\Delta\mu$. However, the values around the
transition $U\sim\Delta\mu$ (top panel of Fig. 7) are much less accurate than
for the variational method. This can be rationalized by comparing the energies
obtained with and without trotterizing the SW transformation (top and bottom
panels of Fig. 7, respectively). Indeed, $\bar{H}^{(N_{s})}$ obtained without
trotterization is block-diagonal (appart from some small deviation for a few
points), although there are some interchanges between the nature of the ground
state around the transition $U\sim\Delta\mu$ in contrast to the variational
approach. Trotterizing the iterative SW transformation does lead to
significant errors and to a non-block-diagonalized $\bar{H}^{(N_{s})}$, such
that the Heisenberg and the ionic states are not eigenstates of
$\bar{H}^{(N_{s})}$ anymore. Such trotterization errors are much more
pronounced within the iterative method than the variational one for two
reasons. On the one hand, the successive applications of more than one (most
of the time, 3 iterations in this work) unitary transformation does multiply
the number of operators that have to be trotterized. On the other hand, it is
known that the variational optimization of the parameters in VQE-based
algorithms does compensate the trotter errors Grimsley  _et al._ (2019).

###  IV.3 Variational versus iterative approach: numerical efficiency

In contrast to the variational approach, the iterative approach leads to a
fully quantum (parameter-free) algorithm as it simply consists in applying the
unitary transformation of Eq. (39) on a prepared eigenstate of $\bar{H}^{\rm
MSW}=\bar{H}^{(N_{s})}$. However, the associated quantum circuit is much
deeper than for the variational approach, which applies a single unitary
transformation only. In terms of gate complexity, the number of CNOT required
to implement a single SW unitary transformation scales with the number of
Pauli terms in the SW generator, as well as with the number of qubits as shown
by the cascade of CNOT in panel (c) of Fig. 3. To evaluate the relevance of
our approach to more complex systems, we extrapolate the computational scaling
for a $N$-sites Hubbard model. As the operators of the SW generator only act
on nearest-neighbor sites, the number of Pauli terms scales as
$\mathcal{O}(N)$. For the iterative approach, one has to multiply by the
number of iterations $N_{s}$. Hence, the number of CNOT scales as
$\mathcal{O}(N^{2})$ and $\mathcal{O}(N_{s}N^{2})$ for the variational and the
iterative approach, respectively. Although the variational approach is more
attractive in the NISQ era due to its shallower circuit depth, it is at the
expense of much more measurements as it has to be multiplied at least by the
number of iterations dictated by the type of cost function and the method used
for the classical optimization of the circuit parameters. Note that while the
iterative approach appears less adapted to NISQ computers, its associated
circuit depth still remains far shallower than quantum phase estimation based
approaches, such as the fault-tolerant one proposed in Ref. [Zhang  _et al._ ,
2022]. Which method is the most efficient will depend on the ability of the
considered quantum computer to afford deep quantum circuit. Within noisy
quantum computer, the variational approach appears more adapted, while the
iterative method can be used on fault tolerant devices.

##  V Conclusions and perspectives

In this paper, we derived recursive relations for the Schrieffer–Wolff
transformation applied on the half-filled Hubbard dimer. Based on these
findings, we proposed a variational and an iterative modification of the
standard SW transformation to approach, or even to perform for homogeneous
case, a block-diagonalization of the Hamiltonian. These modified
Schrieffer–Wolff transformations have been used to design two quantum
algorithms that have been implemented and compared on the half-filled Hubbard
dimer. Regarding the extension of this work to design efficient and
alternative quantum algorithms for the general Hubbard model, or even for
other models or ab-initio Hamiltonian, several challenges have to be
addressed.

At this stage, one could directly, and without modification, use the
variational SW Ansatz [Eq. (33)], the iterative Ansatz [Eq. (39)] or a
combination of both to evaluate the ground-state energy of a given Hubbard
model. Beside the fact that it consists in a serious approximation, as
additional terms in the perturbative expansion will implicitly be neglected
for Hubbard models larger than two sites, it also requires to prepare a
relevant trial state $\ket{\Phi}$ that generalizes the Heisenberg state used
for the Hubbard dimer. If no trivial and easy-to-prepare trial eigenstates of
$\hat{H}^{0}$ are known, this step could be performed variationally using the
VQE algorithm, for instance. Within this strategy, one can expect valuable
results for the regime of large $U/t$ values and close to half-filling.
Alternatively, one could apply the modified SW transformations on a few
relevant and easy-to-prepare states that we know belong to the low-energy
subspace we are interested in, thus forming a basis on which all the
Hamiltonian matrix elements are measured on the quantum computer, followed by
a classically diagonalization in the same spirit of the quantum subspace
diagonalization methods McClean  _et al._ (2017); Motta  _et al._ (2020);
Stair  _et al._ (2020).

Finally, the generalization of this work to any filled Hubbard model or to the
Quantum Chemistry Hamiltonian probably requires to improve the generator. This
could be done for instance, in the spirit of coupled cluster approaches, by
introducing more complex terms or more variational parameters. All the
aforementioned developments are beyond the scope of this manuscript and are
left for future work.

###### Acknowledgements.

The authors would like to thank the ANR (Grant No. ANR-19-CE29-0002 DESCARTES
project) for funding.

##  Appendix A Recursive relations for the standard Schrieffer–Wolff
transformation

As originally presented by SW, defining $\hat{S}$ with the SW condition (18)
is relevant to construct low energy effective Hamiltonian close to the atomic
limit $U/t\rightarrow\infty$. However, as shown in the following, the use of
the generator $\hat{S}$ rapidly leads to off-block-diagonal terms in the
perturbative expansion that are far from negligible, already at intermediate
value of $U/t$. More precisely, second-order ($n=2$) contributions read

$\displaystyle\bar{H}^{(2)}=\dfrac{1}{2}\mathcal{S}(\hat{V})=\dfrac{1}{2}\left[\hat{S},\hat{V}\right]=\bar{H}^{(2),{\rm
dia}}+\bar{H}^{(2),{\rm ex}}+\bar{H}^{(2),{\rm de}},$ (42)

where

$\displaystyle\bar{H}^{(2),{\rm
dia}}=\dfrac{1}{2}\sum_{\begin{subarray}{c}i\neq j\\\
\sigma\end{subarray}}\sum_{x=0}^{3}K^{(1)}_{ij\sigma,x}\hat{p}_{ij\bar{\sigma},x}\left(\hat{n}_{i\sigma}-\hat{n}_{j\sigma}\right)$
(43)

corresponds to diagonal contributions with
$K^{(1)}_{ij\sigma,x}=-2t_{ij}\lambda_{ij\sigma,x}$, and

$\displaystyle\bar{H}^{(2),{\rm
ex}}=\dfrac{1}{2}\sum_{\begin{subarray}{c}i\neq j\\\
\sigma\end{subarray}}J^{(1)}_{ij\sigma}\left(\hat{\gamma}_{ij\sigma}\hat{\gamma}_{ji\bar{\sigma}}+\hat{\gamma}_{ji\sigma}\hat{\gamma}_{ij\bar{\sigma}}\right)$
(44)

corresponds to exchange terms, couples spins of different sites and acts
solely in the Heisenberg subspace with
$J^{(1)}_{ij\sigma}=t_{ij}(\lambda_{ij\sigma,1}-\lambda_{ij\sigma,2})$.
Finally,

$\displaystyle\bar{H}^{(2),{\rm
de}}=\dfrac{1}{2}\sum_{\begin{subarray}{c}i\neq j\\\
\sigma\end{subarray}}L^{(1)}_{ij\sigma}\left(\hat{\gamma}_{ij\sigma}\hat{\gamma}_{ij\bar{\sigma}}+\hat{\gamma}_{ji\sigma}\hat{\gamma}_{ji\bar{\sigma}}\right)$
(45)

creates and annihilates double occupations, with
$L^{(1)}_{ij\sigma}=-t_{ij}(\lambda_{ij\sigma,1}-\lambda_{ij\sigma,2})$, and
thus does not act on the Heisenberg subspace at half band filling. At this
stage, the perturbed Hamiltonian is stable through the Heisenberg space. The
third-order contributions lead to

$\displaystyle\bar{H}^{(3)}=\dfrac{1}{3}\mathcal{S}^{2}(\hat{V})=\dfrac{1}{3}\left[\hat{S},\left[\hat{S},\hat{V}\right]\right]=\bar{H}^{(3),\rm
cpl},$ (46)

where

$\displaystyle\bar{H}^{(3),{\rm
cpl}}=\dfrac{1}{2}\sum_{\begin{subarray}{c}i\neq j\\\
\sigma\end{subarray}}\sum_{x=0}^{3}T^{(2)}_{ij\sigma,x}\hat{p}_{ij\bar{\sigma},x}\left(\hat{\gamma}_{ij\sigma}+\hat{\gamma}_{ji\sigma}\right),$
(47)

with

$\displaystyle T^{(2)}_{ij\sigma,0}$ $\displaystyle=$ $\displaystyle
4t_{ij}\lambda^{2}_{ij\sigma,0},$ $\displaystyle T^{(2)}_{ij\sigma,1}$
$\displaystyle=$ $\displaystyle
2t_{ij}[\lambda_{ij\sigma,1}(2\lambda_{ij\sigma,1}+(\lambda_{ij\sigma,1}-\lambda_{ij\sigma,2}))$
$\displaystyle+(\lambda_{ij\sigma,1}-\lambda_{ij\sigma,2})^{2}],$
$\displaystyle T^{(2)}_{ij\sigma,2}$ $\displaystyle=$ $\displaystyle
2t_{ij}[\lambda_{ij\sigma,2}(2\lambda_{ij\sigma,2}+(\lambda_{ij\sigma,2}-\lambda_{ij\sigma,1}))$
(48) $\displaystyle+(\lambda_{ij\sigma,2}-\lambda_{ij\sigma,1})^{2}],$

and

$\displaystyle T^{(2)}_{ij\sigma,3}$ $\displaystyle=$ $\displaystyle
4t_{ij}\lambda^{2}_{ij\sigma,3}.$ (49)

Obviously, $\bar{H}^{(3),\rm cpl}$ couples states from the Heisenberg subspace
to the other states belonging to the complementary subspace. Hence, already by
truncating at the third order, the similar Hamiltonian
$\bar{H}=H^{0}+\sum_{n=2}^{3}\bar{H}^{(n)}$ is not block-diagonal (i.e stable)
anymore with respect to the Heisenberg subspace, as one can expect since the
SW generator $\hat{S}$ is designed to keep $\bar{H}$ block-diagonal at first
order only, see Fig 1.
Interestingly, the form of $\bar{H}^{(3),\rm cpl}$ in Eq. (47) is analogous to
the original $\hat{V}$ which can be rewritten as

$\hat{V}=\dfrac{1}{2}\sum_{ij}\sum_{x=0}^{3}t_{ij\sigma,x}\hat{p}_{ij\bar{\sigma},x}\left(\hat{\gamma}_{ij\sigma}+\hat{\gamma}_{ji\sigma}\right),$
(50)

with $t_{ij\sigma,x}=-t_{ij}$, $\forall 0\leq x\leq 3$. Given the SW
generator, see in Eq. (19),

$\hat{S}=\dfrac{1}{2}\sum_{i\neq
j,\sigma}\sum_{x=0}^{3}\lambda_{ij\sigma,x}\hat{p}_{ij\bar{\sigma},x}\left(\hat{\gamma}_{ij\sigma}-\hat{\gamma}_{ji\sigma}\right),$
(51)

where the coefficients $\lambda_{ij\sigma,x}$ are defined in Eqs. (21)–(24),
the definition $\hat{H}=\hat{H}^{0}+\hat{V}$ and the unitary transformation
$e^{\hat{S}}$, one obtains the transformed Hamiltonian

$\displaystyle\bar{H}$ $\displaystyle=e^{\hat{S}}\hat{H}e^{-\hat{S}}$
$\displaystyle=\left(\sum_{n}\frac{\hat{S}^{n}}{n!}\right)\hat{H}\left(\sum_{n}\frac{(-1)^{n}\hat{S}^{n}}{n!}\right)$
$\displaystyle=\hat{H}^{0}+\sum_{n=2}^{\infty}\frac{n-1}{n!}\mathcal{S}^{n-1}(\hat{V}),$
(52)

where $\mathcal{S}(\hat{X})=[\hat{S},\hat{X}]$ is a super-operator that
transforms $\hat{X}$ into another operator acting in the same Hilbert space
Primas (1963). Within the previous definition of $\hat{S}$ and $\hat{V}$, it
can be shown that even orders of Eq. (52) take the following form,

$\mathcal{S}^{2n}(\hat{V})=\dfrac{1}{2}\sum_{i\neq
j\sigma}\sum_{x=0}^{3}T_{ij\sigma,x}^{(2n)}\hat{p}_{ij\bar{\sigma},x}\left(\hat{\gamma}_{ij\sigma}+\hat{\gamma}_{ji\sigma}\right),$
(53)

while odd orders are written as

$\displaystyle\mathcal{S}^{2n+1}(\hat{V})$ $\displaystyle=$
$\displaystyle\dfrac{1}{2}\sum_{\begin{subarray}{c}i\neq j\\\
\sigma\end{subarray}}\sum_{x=0}^{3}K^{(2n+1)}_{ij\sigma,x}\hat{p}_{ij\bar{\sigma},x}\left(\hat{n}_{i\sigma}-\hat{n}_{j\sigma}\right)$
(54) $\displaystyle+\frac{1}{2}\sum_{\begin{subarray}{c}i\neq j\\\
\sigma\end{subarray}}J^{(2n+1)}_{ij\sigma}\left(\hat{\gamma}_{ij\sigma}\hat{\gamma}_{ji\bar{\sigma}}+\hat{\gamma}_{ji\sigma}\hat{\gamma}_{ij\bar{\sigma}}\right)$
$\displaystyle+\dfrac{1}{2}\sum_{\begin{subarray}{c}i\neq j\\\
\sigma\end{subarray}}L^{(2n+1)}_{ij\sigma}\left(\hat{\gamma}_{ij\sigma}\hat{\gamma}_{ij\bar{\sigma}}+\hat{\gamma}_{ji\sigma}\hat{\gamma}_{ji\bar{\sigma}}\right).$

Moreover, interaction integrals in Eqs. (53) and (54) can be recursively
obtained as follows,

$\displaystyle K^{(2n+1)}_{ij\sigma,x}=$ $\displaystyle
2\lambda_{ij\sigma,x}T_{ij\sigma,x}^{(2n)},$ (55) $\displaystyle
J^{(2n+1)}_{ij\sigma}=$
$\displaystyle\dfrac{1}{2}\left(K^{2n+1}_{ij\sigma,2}-K^{(2n+1)}_{ij\sigma,1}\right),$
(56) $\displaystyle L^{(2n+1)}_{ij\sigma}=$
$\displaystyle\dfrac{1}{2}\left(K^{(2n+1)}_{ij\sigma,2}\beta-K^{(2n+1)}_{ij\sigma,1}/\beta\right),$
(57) $\displaystyle T_{ij\sigma,0}^{(2n)}=$
$\displaystyle-2\lambda_{ij\sigma,0}K^{(2n-1)}_{ij\sigma,0}=-4\lambda_{ij\sigma,0}^{2}T_{ij\sigma,0}^{(2n-2)}$
$\displaystyle T_{ij\sigma,0}^{(2n)}=$
$\displaystyle(-1)^{n}t_{ij\sigma,0}(2\lambda_{ij\sigma,0})^{2n},$ (58)
$\displaystyle T_{ij\sigma,1}^{(2n)}=$
$\displaystyle-\lambda_{ij\sigma,1}\left(3K^{(2n-1)}_{ij\sigma,1}-K^{(2n-1)}_{ij\sigma,2}-2J^{(2n-1)}_{ij\sigma}\right)$
$\displaystyle+2\lambda_{ij\sigma,2}L^{(2n-1)}_{ij\sigma},$ (59)
$\displaystyle T_{ij\sigma,2}^{(2n)}=$
$\displaystyle-\lambda_{ij\sigma,2}\left(3K^{(2n-1)}_{ij\sigma,2}-K^{(2n-1)}_{ij\sigma,1}+2J^{(2n-1)}_{ij\sigma}\right)$
$\displaystyle-2\lambda_{ij\sigma,1}L^{(2n-1)}_{ij\sigma},$ (60)

and, finally,

$\displaystyle T_{ij\sigma,3}^{(2n)}=$
$\displaystyle-2\lambda_{ij\sigma,3}K^{(2n-1)}_{ij\sigma,3}=-4\lambda_{ij\sigma,3}^{2}T_{ij\sigma,3}^{(2n-2)}$
$\displaystyle=$
$\displaystyle(-1)^{n}t_{ij\sigma,3}(2\lambda_{ij\sigma,3})^{2n},$ (61)

with $\beta_{ij\sigma}=\lambda_{ij\sigma,1}/\lambda_{ij\sigma,2}$ and the
initial condition $T_{ij\sigma,x}^{(0)}=t_{ij\sigma,x}$.
We uncouple Eqs. (59) and (60) by introducing the following geometric series,

$\displaystyle W_{1}^{(2n)}$ $\displaystyle=$
$\displaystyle\lambda_{ij\sigma,2}T_{ij\sigma,1}^{(2n)}+\lambda_{ij\sigma,1}T_{ij\sigma,2}^{(2n)}$
(62) $\displaystyle=$ $\displaystyle W_{1}^{(0)}(-1)^{n}(2\alpha)^{2n},$

and

$\displaystyle W_{2}^{(2n)}$ $\displaystyle=$
$\displaystyle\lambda_{ij\sigma,1}T_{ij\sigma,1}^{(2n)}-\lambda_{ij\sigma,2}T_{ij\sigma,2}^{(2n)}$
(63) $\displaystyle=$ $\displaystyle W_{2}^{(0)}(-1)^{n}(4\alpha)^{2n},$

with
$\alpha^{2}=\left(\lambda_{ij\sigma,1}^{2}+\lambda_{ij\sigma,2}^{2}\right)/2$
and the initial conditions
$W_{1}^{(0)}=\lambda_{ij\sigma,2}t_{ij\sigma,1}+\lambda_{ij\sigma,1}t_{ij\sigma,2}$
and
$W_{2}^{(0)}=\lambda_{ij\sigma,1}t_{ij\sigma,1}-\lambda_{ij\sigma,2}t_{ij\sigma,2}$,
thus leading to

$\displaystyle
T_{ij\sigma,1}^{(2n)}=\frac{(-1)^{n}}{\lambda_{ij\sigma,2}\left(1+\beta^{2}\right)}\left(W_{1}^{(0)}(2\alpha)^{2n}+W_{2}^{(0)}(4\alpha)^{2n}\beta\right),$
(64)

and

$\displaystyle
T_{ij\sigma,2}^{(2n)}=\frac{(-1)^{n}}{\lambda_{ij\sigma,2}\left(1+1/\beta^{2}\right)}\left(W_{1}^{(0)}(2\alpha)^{2n}-W_{2}^{(0)}(4\alpha)^{2n}/\beta\right).$
(65)

Summing all odd and even contributions at the infinite order allows to recover
Eq. (28) to Eq. (29) of the manuscript. The values of the different
interaction integrals at infinite order can be obtained using the recursive
relations, for instance for $T_{ij\sigma,0}$,

$\displaystyle T_{ij\sigma,0}$
$\displaystyle=\sum_{n=0}^{\infty}\frac{2n}{(2n+1)!}T^{(2n)}_{ij\sigma,0}$
$\displaystyle=t_{ij\sigma,0}\sum_{n=0}^{\infty}\frac{2n(-1)^{n}}{(2n+1)!}(2\lambda_{ij\sigma,0})^{2n}$
$\displaystyle=t_{ij\sigma,0}\sum_{n=0}^{\infty}\frac{(-1)^{n}}{(2n)!}(2\lambda_{ij\sigma,0})^{2n}$
$\displaystyle\quad-\frac{t_{ij\sigma,0}}{\lambda_{ij\sigma,0}}\sum_{n=0}^{\infty}\frac{(-1)^{n}}{(2n+1)!}(2\lambda_{ij\sigma,0})^{2n+1}$
$\displaystyle=t_{ij\sigma,0}\left(\cos(2\lambda_{ij\sigma,0})-{\rm
sinc}(2\lambda_{ij\sigma,0})\right).$ (66)

The other integrals are similarly obtained and reads, for the electronic
integrals $T$ corresponding to the coupling between the Heisenberg subspace
and its complementary subspace,

$\displaystyle T_{ij\sigma,1}=$
$\displaystyle\frac{W_{1}^{(0)}}{\lambda_{ij\sigma,2}(1+\beta_{ij\sigma}^{2})}\left(\cos(2\alpha)-{\rm
sinc}(2\alpha)\right)$
$\displaystyle+\frac{W_{2}^{(0)}}{\lambda_{ij\sigma,1}(1+1/\beta_{ij\sigma}^{2})}\left(\cos(4\alpha)-{\rm
sinc}(4\alpha)\right),$ (67) $\displaystyle T_{ij\sigma,2}=$
$\displaystyle\frac{W_{1}^{(0)}}{\lambda_{ij\sigma,1}(1+1/\beta_{ij\sigma}^{2})}\left(\cos(2\alpha)-{\rm
sinc}(2\alpha)\right)$
$\displaystyle-\frac{W_{2}^{(0)}}{\lambda_{ij\sigma,2}(1+\beta_{ij\sigma}^{2})}\left(\cos(4\alpha)-{\rm
sinc}(4\alpha)\right),$ (68)
$T_{ij\sigma,3}=t_{ij\sigma,3}\left(\cos(2\lambda_{ij\sigma,3})-{\rm
sinc}(2\lambda_{ij\sigma,3})\right),$ (69)

for the electronic integrals $K$ associated to operators that acts diagonally
on each subspace,

$K_{ij\sigma,0}=t_{ij\sigma,0}\left(\sin(2\lambda_{ij\sigma,0})-\lambda_{ij\sigma,0}\,{\rm
sinc}^{2}(\lambda_{ij\sigma,0})\right),$ (70) $\displaystyle K_{ij\sigma,1}=$
$\displaystyle\frac{W_{1}^{(0)}}{\beta_{ij\sigma}+1/\beta_{ij\sigma}}\left(2\,{\rm
sinc}(2\alpha)-{\rm sinc}^{2}(\alpha)\right)$
$\displaystyle\frac{W_{2}^{(0)}}{1+1/\beta_{ij\sigma}^{2}}\left(2\,{\rm
sinc}(4\alpha)-{\rm sinc}^{2}(2\alpha)\right),$ (71) $\displaystyle
K_{ij\sigma,2}=$
$\displaystyle\frac{W_{1}^{(0)}}{\beta_{ij\sigma}+1/\beta_{ij\sigma}}\left(2\,{\rm
sinc}(2\alpha)-{\rm sinc}^{2}(\alpha)\right)$
$\displaystyle-\frac{W_{2}^{(0)}}{1+\beta_{ij\sigma}^{2}}\left(2\,{\rm
sinc}(4\alpha)-{\rm sinc}^{2}(2\alpha)\right),$ (72)
$K_{ij\sigma,3}=t_{ij\sigma,3}\left(\sin(2\lambda_{ij\sigma,3})-\lambda_{ij\sigma,3}\,{\rm
sinc}^{2}(\lambda_{ij\sigma,3})\right),$ (73)

and finally, for the electronic integrals $J$ and $L$ associated to spin-spin
operators and doubly occupied- empty sites operators,

$J_{ij\sigma}=(1/2)\left(K_{ij\sigma,1}-K_{ij\sigma,2}\right),$ (74)

and

$L_{ij\sigma}=(1/2)\left(K_{ij\sigma,1}/\beta_{ij\sigma}-\beta_{ij\sigma}K_{ij\sigma,2}\right),$
(75)

respectively. For the homogeneous case, the associated integrals are simply
given by

$\displaystyle T_{ij\sigma,0}=T_{ij\sigma,3}=K_{ij\sigma,0}=K_{ij\sigma,3}=0,$
(76) $\displaystyle T_{ij\sigma,1}=T_{ij\sigma,2}=-t\left(\cos{(4t/U)}-{\rm
sinc}\,(4t/U)\right),$ (77)

and

$\displaystyle J_{ij\sigma}=K_{ij\sigma,1}$
$\displaystyle=-K_{ij\sigma,2}=-L_{ij\sigma}$
$\displaystyle=-t\left(\sin{(4t/U)}-(2t/U){\rm sinc}^{2}\,(2t/U)\right)/2.$
(78)

##  Appendix B Recursive relations for the variational Schrieffer–Wolff
transformation

The variational SW transformation given in Eq. (33) leads to the following
transformed Hamiltonian,

$\displaystyle\bar{H}(\theta)$
$\displaystyle=e^{\theta\hat{S}}\hat{H}e^{-\theta\hat{S}}$
$\displaystyle=\hat{H}^{0}+\hat{V}\left(1-\theta\right)+\sum_{n=2}^{\infty}\frac{\theta^{n-1}(n-\theta)}{n!}\mathcal{S}^{n-1}(\hat{V}).$
(79)

Interestingly, Eqs. (53) and (54) still hold in this case, such that a
strategy similar to the one introduced in Appendix A can be used to obtain
interaction integrals at each order, recursively. The summation to the
infinite order is then possible, thus leading to Eq. (35) with the electronic
integrals $T(\theta)$ corresponding to the coupling between the Heisenberg
subspace and its complementary subspace,

$T_{ij\sigma,0}(\theta)=t_{ij\sigma,0}\left(\cos(2\theta\lambda_{ij\sigma,0})-\theta{\rm
sinc}(2\theta\lambda_{ij\sigma,0})\right),$ (80) $\displaystyle
T_{ij\sigma,1}(\theta)=$
$\displaystyle\frac{W_{1}^{(0)}}{\lambda_{ij\sigma,2}(1+\beta_{ij\sigma}^{2})}\left(\cos(2\theta\alpha)-\theta{\rm
sinc}(2\theta\alpha)\right)$
$\displaystyle+\frac{W_{2}^{(0)}}{\lambda_{ij\sigma,1}(1+1/\beta_{ij\sigma}^{2})}\left(\cos(4\theta\alpha)-\theta{\rm
sinc}(4\theta\alpha)\right),$ (81) $\displaystyle T_{ij\sigma,2}(\theta)=$
$\displaystyle\frac{W_{1}^{(0)}}{\lambda_{ij\sigma,1}(1+1/\beta_{ij\sigma}^{2})}\left(\cos(2\theta\alpha)-\theta{\rm
sinc}(2\theta\alpha)\right)$
$\displaystyle-\frac{W_{2}^{(0)}}{\lambda_{ij\sigma,2}(1+\beta_{ij\sigma}^{2})}\left(\cos(4\theta\alpha)-\theta{\rm
sinc}(4\theta\alpha)\right),$ (82)
$T_{ij\sigma,3}(\theta)=t_{ij\sigma,3}\left(\cos(2\theta\lambda_{ij\sigma,3})-\theta{\rm
sinc}(2\theta\lambda_{ij\sigma,3})\right),$ (83)

the electronic integrals $K(\theta)$ associated to operators that acts
diagonally on each subspace,

$K_{ij\sigma,0}(\theta)=t_{ij\sigma,0}\left(\sin(2\theta\lambda_{ij\sigma,0})-\lambda_{ij\sigma,0}\theta^{2}\,{\rm
sinc}^{2}(\theta\lambda_{ij\sigma,0})\right),$ (84) $\displaystyle
K_{ij\sigma,1}(\theta)=$
$\displaystyle\frac{W_{1}^{(0)}\theta}{\beta_{ij\sigma}+1/\beta_{ij\sigma}}\left(2\,{\rm
sinc}(2\theta\alpha)-\theta{\rm sinc}^{2}(\theta\alpha)\right)$
$\displaystyle\frac{W_{2}^{(0)}\theta}{1+1/\beta_{ij\sigma}^{2}}\left(2\,{\rm
sinc}(4\theta\alpha)-\theta{\rm sinc}^{2}(2\theta\alpha)\right)$ (85)
$\displaystyle K_{ij\sigma,2}(\theta)=$
$\displaystyle\frac{W_{1}^{(0)}\theta}{\beta_{ij\sigma}+1/\beta_{ij\sigma}}\left(2\,{\rm
sinc}(2\theta\alpha)-\theta{\rm sinc}^{2}(\theta\alpha)\right)$
$\displaystyle-\frac{W_{2}^{(0)}\theta}{1+\beta_{ij\sigma}^{2}}\left(2\,{\rm
sinc}(4\theta\alpha)-\theta{\rm sinc}^{2}(2\theta\alpha)\right)$ (86)
$K_{ij\sigma,3}(\theta)=t_{ij\sigma,3}\left(\sin(2\theta\lambda_{ij\sigma,3})-\lambda_{ij\sigma,3}\theta^{2}\,{\rm
sinc}^{2}(\theta\lambda_{ij\sigma,3})\right),$ (87)

and finally, the electronic integrals $J(\theta)$ and $L(\theta)$ associated
to spin-spin operators and doubly occupied-empty sites operators,

$J_{ij\sigma}(\theta)=(1/2)\left(K_{ij\sigma,1}(\theta)-K_{ij\sigma,2}(\theta)\right),$
(88)

and

$L_{ij\sigma}(\theta)=(1/2)\left(K_{ij\sigma,1}(\theta)/\beta_{ij\sigma}-\beta_{ij\sigma}K_{ij\sigma,2}(\theta)\right),$
(89)

respectively.

In the following we show that for the homogenous case and at the saddle point,
the variational SW transformation with generator $\hat{G}=\theta\hat{S}$
fulfills the VV condition in Eq. (10). More precisely, in this case, Eq. (10)
reads

$\displaystyle[\theta\hat{S},\hat{H}^{0}]=-\hat{V}-\sum_{n=1}^{\infty}c_{n}\theta^{2n}\mathcal{S}^{2n}(\hat{V}),$
$\displaystyle\theta\hat{V}=\hat{V}+\sum_{n=1}^{\infty}c_{n}\theta^{2n}\mathcal{S}^{2n}(\hat{V}).$
(90)

Inserting Eq. (53) on the right handside of Eq. (90) with
$T_{ij\sigma,0}^{(2n)}=T_{ij\sigma,3}^{(2n)}=0$ and
$T_{ij\sigma,1}^{(2n)}=-T_{ij\sigma,2}^{(2n)}=(-1)^{n}4t/U$ for the
homogeneous case [see Eqs. (58), (A), (64) and (65)] leads to

$\displaystyle\theta\hat{V}=\left(1+\sum_{n=1}^{\infty}c_{n}(-1)^{n}\left(\dfrac{4t\theta}{U}\right)^{2n}\right)\hat{V}$
(91)

Using $c_{n}=B_{2n}2^{2n}/(2n)!$ $=(-1)^{n-1}|B_{2n}|2^{2n}/(2n)!$, we obtain

$\displaystyle\theta=1-\sum_{n=1}^{\infty}\dfrac{|B_{2n}|2^{2n}}{(2n)!}\left(\dfrac{4t\theta}{U}\right)^{2n}$
(92) $\displaystyle\theta=\left(\dfrac{4t\theta}{U}\right){\rm
cotan}\left(\dfrac{4t\theta}{U}\right)$ (93)

or, equivalently, $\theta=(U/4t)\tan^{-1}{(4t/U)}$, which is the analytical
expression of the minimizing variational parameter in Eq. (38).

##  Appendix C Iterative generator

In this section, we detail the construction of the iterative transformation,
in the spirit of the Foldy–Wouthuysen transformation Foldy and Wouthuysen
(1950):

$U^{(N_{s})}=\left(\prod_{s=0}^{N_{s}-1}e^{\hat{S}^{(s)}}\right).$ (94)

At each iteration, the generator takes the form

$\hat{S}^{(s)}=\frac{1}{2}\sum_{i\neq
j,\sigma}\sum_{x}\lambda_{ij\sigma,x}^{(s)}\hat{p}_{ij\bar{\sigma},x}\left(\hat{\gamma}_{ij\sigma}-\hat{\gamma}_{ji\sigma}\right),$
(95)

and the iteration $s=0$ refers to the standard SW transformation at infinite
order, for which $\lambda_{ij\sigma,x}^{(0)}$ parameters are set to satisfy
Eq. (11), i.e.,

$\displaystyle[\hat{H}^{0},\hat{S}^{(0)}]=\hat{V}_{X}.$ (96)

This equation can be rewritten as follows,

$\displaystyle\sum_{i\neq
j,\sigma}\sum_{x=0}^{3}\left(f_{ij\sigma,x}^{(0)}(\lambda_{ij\sigma,x}^{(0)})-t_{ij\sigma,x}\right)\hat{p}_{ij\bar{\sigma},x}\left(\hat{\gamma}_{ij\sigma}+\hat{\gamma}_{ji\sigma}\right)=\hat{0},$
(97)

with

$\displaystyle f_{ij\sigma,0}^{(0)}$ $\displaystyle=$
$\displaystyle\lambda_{ij\sigma,0}^{(0)}\Delta\mu_{ij}/2,$ $\displaystyle
f_{ij\sigma,1}^{(0)}$ $\displaystyle=$
$\displaystyle\lambda_{ij\sigma,1}^{(0)}(\Delta\mu_{ij}+U_{i})/2,$
$\displaystyle f_{ij\sigma,2}^{(0)}$ $\displaystyle=$
$\displaystyle\lambda_{ij\sigma,2}^{(0)}(\Delta\mu_{ij}-U_{j})/2,$
$\displaystyle f_{ij\sigma,3}^{(0)}$ $\displaystyle=$
$\displaystyle\lambda_{ij\sigma,x}^{(0)}(\Delta\mu_{ij}+\Delta U_{ij})/2$ (98)

thus leading to Eqs. (21)–(24).

Through the recursive relations, we established that
$\bar{H}^{0}=e^{\hat{S}^{(0)}}\hat{H}e^{-\hat{S}^{(0)}}$ is given by Eq. (28)
where $\bar{H}^{\rm cpl}$ is off-block-diagonal. The iterative process
consists in repeating SW-type unitary transformation, providing that
$\hat{H}^{0}+\bar{H}^{{\rm dia}(s)}+\bar{H}^{{\rm ex}(s)}+\bar{H}^{{\rm
de}(s)}\rightarrow\hat{H}^{0(s+1)}$ and $\bar{H}^{{\rm
cpl}(s)}\rightarrow\hat{V}^{(s+1)}$. Eq. (96) then defines the $(s+1)$ SW
generator,

$\displaystyle[\hat{H}^{0(s+1)},\hat{S}^{(s+1)}]=\hat{V}^{(s+1)},$
$\displaystyle[\hat{H}^{0}+\bar{H}^{{\rm dia}(s)}+\bar{H}^{{\rm
ex}(s)}+\bar{H}^{{\rm de}(s)},\hat{S}^{(s+1)}]=\bar{H}^{{\rm cpl}(s)},$ (99)

that, similarly to Eq. (97), reduces a set of linearly coupled equations,

$\displaystyle f_{ij\sigma,x}^{(s+1)}-T^{(s)}_{ij\sigma,x}=0,$ (100)

or, with more details,

$\displaystyle
f_{ij\sigma,0}^{(s+1)}=\lambda_{ij\sigma,0}^{(s+1)}\left(\Delta\mu_{ij}+2K^{(s)}_{ij\sigma,0}\right),$
$\displaystyle f_{ij\sigma,1}^{(s+1)}$
$\displaystyle=2\lambda_{ij\sigma,2}^{(s+1)}L^{(s)}_{ij}+\lambda_{ij\sigma,1}^{(s+1)}\times$
$\displaystyle\left(\Delta\mu_{ij}+U_{i}+3K^{(s)}_{ij\sigma,1}-K^{(s)}_{ij\sigma,2}+2J^{(s)}_{ij}\right),$
$\displaystyle f_{ij\sigma,2}^{(s+1)}$
$\displaystyle=-2\lambda_{ij\sigma,1}^{(s+1)}L^{(s)}_{ij\sigma}+\lambda_{ij\sigma,2}^{(s+1)}\times$
$\displaystyle\left(\Delta\mu_{ij}-U_{j}+3K^{(s)}_{ij\sigma,2}-K^{(s)}_{ij\sigma,1}-2J^{(s)}_{ij}\right),$

and, finally,

$\displaystyle
f_{ij\sigma,3}^{(s+1)}=\lambda_{ij\sigma,3}^{(s+1)}\left(\Delta\mu_{ij}+\Delta
U_{ij}+2K^{(s)}_{ij\sigma,3}\right).$

It follows straightforwardly that

$\displaystyle\lambda_{ij\sigma,0}^{(s+1)}=T^{(s)}_{ij\sigma,0}/(\Delta\mu_{ij}+2K^{(s)}_{ij\sigma,0}),$
(101)
$\displaystyle\lambda_{ij\sigma,1}^{(s+1)}=\frac{2T^{(s)}_{ij\sigma,2}L^{(s)}_{ij\sigma}+T^{(s)}_{ij\sigma,1}B_{2,ij\sigma}^{(s)}}{4L^{(s)}_{ij\sigma}L^{(s)}_{ij\sigma}+B_{1,ij\sigma}^{(s)}B_{2,ij\sigma}^{(s)}},$
(102)
$\displaystyle\lambda_{ij\sigma,2}^{(s+1)}=\frac{-2T^{(s)}_{ij\sigma,1}L^{(s)}_{ij\sigma}+T^{(s)}_{ij\sigma,2}B_{1,ij\sigma}^{(s)}}{4L^{(s)}_{ij\sigma}L^{(s)}_{ij\sigma}+B_{1,ij\sigma}^{(s)}B_{2,ij\sigma}^{(s)}},$
(103)

and

$\displaystyle\lambda_{ij\sigma,3}^{(s+1)}=T^{(s)}_{ij\sigma,3}/(\Delta\mu_{ij}+\Delta
U_{ij}+2K^{(s)}_{ij\sigma,3}),$ (104)

with
$B_{1,ij\sigma}^{(s)}=(\Delta\mu_{ij}+U_{i})+3K^{(s)}_{ij\sigma,1}-K^{(s)}_{ij\sigma,2}+2J^{(s)}_{ij\sigma}$
and
$B_{2,ij\sigma}^{(s)}=(\Delta\mu_{ij}-U_{j})+3K^{(s)}_{ij\sigma,2}-K^{(s)}_{ij\sigma,1}-2J^{(s)}_{ij\sigma}$
.

##  Appendix D Cost functions

Figure 8: Energies of the half-filled Hubbard dimer for $\Delta\mu/t=2$ (top
panel) with respect to the repulsion strength, using the variational MSW
transformation method to minimize the matrix elements of
$|\bar{H}_{X}(\theta)|$. The minimizing parameter $\theta_{X}$ is shown in the
bottom panel. The vertical dotted line corresponds to $U=\Delta\mu$.

Figure 9: Top panel: Energies of the half-filled Hubbard dimer for
$\Delta\mu/t=2$ with respect to the repulsion strength, using the variational
SW method to minimize the energy $\bra{\Phi_{\rm
Heis}}\bar{H}(\theta)\ket{\Phi_{\rm Heis}}$ (orange triangles) or
$\langle\Phi_{\rm Ionic}^{\alpha=0}|\bar{H}(\theta)|\Phi_{\rm
Ionic}^{\alpha=0}\rangle$ (blue crosses). Bottom panel: Minimizing parameter
$\theta_{\rm Heis}$ (yellow crosses) and $\theta_{\rm Ionic}$ (blue dots). The
vertical dotted line corresponds to $U=\Delta\mu$.

In the variational approach described in Sec. III.2, one can optimize the
scaling parameter $\theta$ to minimize the contributions from the coupling
operator $\bar{H}_{X}(\theta)$ , or to minimize the energy of
$\bar{H}(\theta)$ restricted to the Heisenberg subspace. In this section, we
investigate the difference of the two strategies.

The first one requires to minimize the Frobenius norm
$||\bar{H}_{X}(\theta)||_{F}$, for which the saddle point $\theta_{X}$ gives
the transformation that maximally decouples the Heisenberg subspace $\Omega$
from the complementary subspace, without any warranty that the low-energy
ground state belongs to $\Omega$. On a quantum computer, it requires the
estimation of off-diagonal elements
$\bra{\Phi_{i}}\bar{H}(\theta)\ket{\Phi_{j}}$, $i$ belonging to $\Omega$ and
$j$ to the complementary subspace, which is generally non trivial, although
one can note some recent improvements made in Refs. [Huggins  _et al._ , 2020]
and [Stair and Evangelista, 2021]. Zhang and coworkers Zhang  _et al._ (2022)
used a similar cost function, evaluated using only the estimation of off-
diagonal terms over states belonging to $\Omega$. In Fig. 8, we plot the
minimizing $\theta_{X}$ and the energies associated to the rotated states
$U^{\dagger}(\theta_{X})\ket{\Phi_{\rm Heis}}$ and
$U^{\dagger}(\theta_{X})\ket{\Phi_{\rm Ionic}}$.

The other method consists in minimizing the energy $\bra{\Phi_{\rm
Heis}}U(\theta)HU^{\dagger}(\theta)\ket{\Phi_{\rm Heis}}$ ($\bra{\Phi_{\rm
Ionic}}U(\theta)HU^{\dagger}(\theta)\ket{\Phi_{\rm Ionic}}$) with respect to
$\theta$ and is analog to the VQE algorithm. The saddle point obtained is
denoted as $\theta_{\rm Heis}$ ($\theta_{\rm Ionic}$) and gives the
transformation that maximally overlaps the rotated Heisenberg (Ionic) state
with the exact ground state. Results are shown in Fig. 9

Comparing Fig. 8 to Fig. 9, one can directly see that minimizing the energy
with respect to the Heisenberg state does work for $U\gg\Delta\mu$ only, where
$\theta_{\rm Heis}\simeq\theta_{X}$, while minimizing the energy with respect
to the ionic state works only for $U\ll\Delta\mu$, where $\theta_{\rm
ionic}\simeq\theta_{X}$. Note that in virtue of the variational principle, the
energies of the states $U^{\dagger}(\theta_{\rm Heis})\ket{\Phi_{\rm Heis}}$
and $U^{\dagger}(\theta_{\rm Ionic})\ket{\Phi_{\rm Ionic}}$ are always below
$U^{\dagger}(\theta_{X})\ket{\Phi_{\rm Heis}}$ and
$U^{\dagger}(\theta_{X})\ket{\Phi_{\rm Ionic}}$, respectively, for the entire
range of correlation regime. Therefore, minimizing
$||\bar{H}_{X}(\theta)||_{F}$ gives the rotation that maximally satisfies the
VV conditions [Eq. (10)] but, in contrast to the VQE-like algorithm, it does
not ensure that the Heisenberg space $\Omega$ maximally overlaps with the low-
effective sub-space.

## References

  * Hubbard (1963) J. Hubbard, Proc. Math. Phys. Eng. Sci. P ROY SOC A-MATH PHY 276, 238 (1963).
  * Nagaoka (1965) Y. Nagaoka, Solid State Commun. 3, 409 (1965).
  * Lieb and Wu (1968) E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. 20, 1445 (1968).
  * López-Sandoval and Pastor (2002) R. López-Sandoval and G. M. Pastor, Phys. Rev. B 66 (2002).
  * Lima  _et al._ (2003) N. A. Lima, M. F. Silva, L. N. Oliveira, and K. Capelle, Phys. Rev. Lett. 90 (2003).
  * Georges  _et al._ (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
  * Sénéchal  _et al._ (2000) D. Sénéchal, D. Perez, and M. Pioro-Ladriere, Phys. Rev. Lett. 84, 522 (2000).
  * Potthoff (2003) M. Potthoff, Eur. Phys. J. B 32, 429 (2003).
  * White (1992) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
  * Knizia and Chan (2012) G. Knizia and G. K.-L. Chan, Phys. Rev. Lett. 109 (2012).
  * Sekaran  _et al._ (2021) S. Sekaran, M. Tsuchiizu, M. Saubanère, and E. Fromager, “Householder transformed density matrix functional embedding theory,”  (2021).
  * Wecker  _et al._ (2015a) D. Wecker, M. B. Hastings, N. Wiebe, B. K. Clark, C. Nayak, and M. Troyer, Phys. Rev. A 92, 062318 (2015a).
  * Wecker  _et al._ (2015b) D. Wecker, M. B. Hastings, and M. Troyer, Phys. Rev. A 92, 042303 (2015b).
  * Kivlichan  _et al._ (2018) I. D. Kivlichan, J. McClean, N. Wiebe, C. Gidney, A. Aspuru-Guzik, G. K.-L. Chan, and R. Babbush, Phys. Rev. Lett. 120, 110501 (2018).
  * Reiner  _et al._ (2019) J.-M. Reiner, F. Wilhelm-Mauch, G. Schön, and M. Marthaler, Quantum Sci. Technol. 4, 035005 (2019).
  * Montanaro and Stanisic (2020) A. Montanaro and S. Stanisic, “Compressed variational quantum eigensolver for the Fermi-Hubbard model,”  (2020).
  * Cai (2020) Z. Cai, Phys. Rev. Applied 14, 014059 (2020).
  * Cade  _et al._ (2020) C. Cade, L. Mineh, A. Montanaro, and S. Stanisic, Phys. Rev. B 102, 235122 (2020).
  * Mineh and Montanaro (2022) L. Mineh and A. Montanaro, Phys. Rev. B 105, 125117 (2022).
  * Martin  _et al._ (2022) B. A. Martin, P. Simon, and M. J. Rančić, Phys. Rev. Research 4, 023190 (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, Nat. Commun. 13, 5743 (2022).
  * Dallaire-Demers  _et al._ (2019) P.-L. Dallaire-Demers, J. Romero, L. Veis, S. Sim, and A. Aspuru-Guzik, Quantum Sci. Technol. 4, 045005 (2019).
  * Dallaire-Demers  _et al._ (2020) P.-L. Dallaire-Demers, M. Stȩchły, J. F. Gonthier, N. T. Bashige, J. Romero, and Y. Cao, arXiv:2003.01862 (2020).
  * Suchsland  _et al._ (2022) P. Suchsland, P. K. Barkoutsos, I. Tavernelli, M. H. Fischer, and T. Neupert, Phys. Rev. Research 4, 013165 (2022).
  * Gard and Meier (2022) B. T. Gard and A. M. Meier, Phys. Rev. A 105, 042602 (2022).
  * Kivlichan  _et al._ (2020) I. D. Kivlichan, C. Gidney, D. W. Berry, N. Wiebe, J. McClean, W. Sun, Z. Jiang, N. Rubin, A. Fowler, A. Aspuru-Guzik, H. Neven, and R. Babbush, Quantum 4, 296 (2020).
  * Campbell (2022) E. T. Campbell, Quantum Sci. Technol. 7, 015007 (2022).
  * Clinton  _et al._ (2021) L. Clinton, J. Bausch, and T. Cubitt, Nat. Commun. 12, 4989 (2021).
  * 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, Nat. Commun. 5, 4213 (2014).
  * Bharti  _et al._ (2022) K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C. Kwek, and A. Aspuru-Guzik, Rev. Mod. Phys. 94, 015004 (2022).
  * Van Vleck (1929) J. H. Van Vleck, Phys. Rev. 33, 467 (1929).
  * Jordahl (1934) O. M. Jordahl, Phys. Rev. 45, 87 (1934).
  * Foldy and Wouthuysen (1950) L. L. Foldy and S. A. Wouthuysen, Phys. Rev. 78, 29 (1950).
  * Primas (1963) H. Primas, Rev. Mod. Phys. 35, 710 (1963).
  * Brandow (1979) B. H. Brandow, Int. J. Quantum Chem. 15, 207 (1979).
  * Shavitt and Redmon (1980) I. Shavitt and L. T. Redmon, J. Chem. Phys. 73, 5711 (1980).
  * Bravyi  _et al._ (2011) S. Bravyi, D. P. DiVincenzo, and D. Loss, Ann. Physik. 326, 2793 (2011).
  * Shavitt and Bartlett (2009) I. Shavitt and R. J. Bartlett,  _Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory_, 1st ed. (Cambridge University Press, 2009).
  * Schrieffer and Wolff (1966) J. R. Schrieffer and P. A. Wolff, Phys. Rev. 149, 491 (1966).
  * Harris and Lange (1967) A. B. Harris and R. V. Lange, Phys. Rev. 157, 295 (1967).
  * Chao  _et al._ (1977) K. A. Chao, J. Spalek, and A. M. Oles, J. Phys. C: Solid State Phys. 10, L271 (1977).
  * Zhang  _et al._ (2022) Z. Zhang, Y. Yang, X. Xu, and Y. Li, Phys. Rev. Research 4, 043023 (2022).
  * Romero  _et al._ (2018) J. Romero, R. Babbush, J. R. McClean, C. Hempel, P. J. Love, and A. Aspuru-Guzik, Quantum Sci. Technol. 4, 014008 (2018).
  * tA v  _et al._ (2021) A. tA v, M. S. ANIS, Abby-Mitchell, H. Abraham, AduOffei, R. Agarwal, G. Agliardi, M. Aharoni, V. Ajith, I. Y. Akhalwaya, G. Aleksandrowicz, T. Alexander, M. Amy, S. Anagolum, Anthony-Gandon, I. F. Araujo, E. Arbel, A. Asfaw, I. E. Ashimine, A. Athalye, A. Avkhadiev, C. Azaustre, P. BHOLE, V. Bajpe, A. Banerjee, S. Banerjee, W. Bang, A. Bansal, P. Barkoutsos, A. Barnawal, G. Barron, G. S. Barron, L. Bello, Y. Ben-Haim, M. C. Bennett, D. Bevenius, D. Bhatnagar, P. Bhatnagar, A. Bhobe, P. Bianchini, L. S. Bishop, C. Blank, S. Bolos, S. Bopardikar, S. Bosch, S. Brandhofer, Brandon, S. Bravyi, Bryce-Fuller, D. Bucher, L. Burgholzer, A. Burov, F. Cabrera, P. Calpin, L. Capelluto, J. Carballo, G. Carrascal, A. Carriker, I. Carvalho, R. Chakrabarti, A. Chen, C.-F. Chen, E. Chen, J. C. Chen, R. Chen, F. Chevallier, K. Chinda, R. Cholarajan, J. M. Chow, S. Churchill, CisterMoke, C. Claus, C. Clauss, C. Clothier, R. Cocking, R. Cocuzzo, J. Connor, F. Correa, Z. Crockett, A. J. Cross, A. W. Cross, S. Cross, J. Cruz-Benito, C. Culver, A. D. Córcoles-Gonzales, N. D, S. Dague, T. E. Dandachi, A. N. Dangwal, J. Daniel, DanielAja, M. Daniels, M. Dartiailh, A. R. Davila, F. Debouni, A. Dekusar, A. Deshmukh, M. Deshpande, D. Ding, J. Doi, E. M. Dow, P. Downing, E. Drechsler, M. S. Drudis, E. Dumitrescu, K. Dumon, I. Duran, K. EL-Safty, E. Eastman, G. Eberle, A. Ebrahimi, P. Eendebak, D. Egger, EgrettaThula, ElePT, I. Elsayed, Emilio, A. Espiricueta, M. Everitt, D. Facoetti, Farida, P. M. Fernández, S. Ferracin, D. Ferrari, A. H. Ferrera, R. Fouilland, A. Frisch, A. Fuhrer, B. Fuller, M. GEORGE, J. Gacon, B. G. Gago, C. Gambella, J. M. Gambetta, A. Gammanpila, L. Garcia, T. Garg, S. Garion, J. R. Garrison, J. Garrison, T. Gates, N. Gavrielov, G. Gentinetta, H. Georgiev, L. Gil, A. Gilliam, A. Giridharan, Glen, J. Gomez-Mosquera, Gonzalo, S. de la Puente González, J. Gorzinski, I. Gould, D. Greenberg, D. Grinko, W. Guan, D. Guijo, Guillermo-Mijares-Vilarino, J. A. Gunnels, H. Gupta, N. Gupta, J. M. Günther, M. Haglund, I. Haide, I. Hamamura, O. C. Hamido, F. Harkins, K. Hartman, A. Hasan, V. Havlicek, J. Hellmers, Ł. Herok, R. Hill, S. Hillmich, I. Hincks, C. Hong, H. Horii, C. Howington, S. Hu, W. Hu, C.-H. Huang, J. Huang, R. Huisman, H. Imai, T. Imamichi, K. Ishizaki, Ishwor, R. Iten, T. Itoko, A. Ivrii, A. Javadi, A. Javadi-Abhari, W. Javed, Q. Jianhua, M. Jivrajani, K. Johns, S. Johnstun, Jonathan-Shoemaker, JosDenmark, JoshDumo, J. Judge, T. Kachmann, A. Kale, N. Kanazawa, J. Kane, Kang-Bae, A. Kapila, A. Karazeev, P. Kassebaum, T. Kato, T. Kehrer, J. Kelso, S. Kelso, H. van Kemenade, V. Khanderao, S. King, Y. Kobayashi, Kovi11Day, A. Kovyrshin, J. Krishna, R. Krishnakumar, P. Krishnamurthy, V. Krishnan, K. Krsulich, P. Kumkar, G. Kus, LNoorl, R. LaRose, E. Lacal, R. Lambert, H. Landa, J. Lapeyre, D. Lasecki, J. Latone, S. Lawrence, C. Lee, G. Li, T. J. Liang, J. Lishman, D. Liu, P. Liu, Lolcroc, A. K. M, L. Madden, Y. Maeng, S. Maheshkar, K. Majmudar, A. Malyshev, M. E. Mandouh, J. Manela, Manjula, J. Marecek, M. Marques, K. Marwaha, D. Maslov, P. Maszota, D. Mathews, A. Matsuo, F. Mazhandu, D. McClure, M. McElaney, J. McElroy, C. McGarry, D. McKay, D. McPherson, S. Meesala, D. Meirom, C. Mendell, T. Metcalfe, M. Mevissen, A. Meyer, A. Mezzacapo, R. Midha, D. Millar, D. Miller, H. Miller, Z. Minev, A. Mitchell, A. Mohammad, N. Moll, A. Montanez, G. Monteiro, M. D. Mooring, R. Morales, N. Moran, D. Morcuende, S. Mostafa, M. Motta, R. Moyard, P. Murali, D. Murata, J. Müggenburg, T. NEMOZ, D. Nadlinger, K. Nakanishi, G. Nannicini, P. Nation, E. Navarro, Y. Naveh, S. W. Neagle, P. Neuweiler, A. Ngoueya, T. Nguyen, J. Nicander, Nick-Singstock, P. Niroula, H. Norlen, NuoWenLei, L. J. O’Riordan, O. Ogunbayo, P. Ollitrault, T. Onodera, R. Otaolea, S. Oud, D. Padilha, H. Paik, S. Pal, Y. Pang, A. Panigrahi, V. R. Pascuzzi, S. Perriello, E. Peterson, A. Phan, K. Pilch, F. Piro, M. Pistoia, C. Piveteau, J. Plewa, P. Pocreau, C. Possel, A. Pozas-Kerstjens, R. Pracht, M. Prokop, V. Prutyanov, S. Puri, D. Puzzuoli, Pythonix, J. Pérez, Quant02, Quintiii, R. I. Rahman, A. Raja, R. Rajeev, I. Rajput, N. Ramagiri, A. Rao, R. Raymond, O. Reardon-Smith, R. M.-C. Redondo, M. Reuter, J. Rice, M. Riedemann, Rietesh, D. Risinger, P. Rivero, M. L. Rocca, D. M. Rodríguez, RohithKarur, B. Rosand, M. Rossmannek, M. Ryu, T. SAPV, N. R. C. Sa, A. Saha, A. Ash-Saki, A. Salman, S. Sanand, M. Sandberg, H. Sandesara, R. Sapra, H. Sargsyan, A. Sarkar, N. Sathaye, N. Savola, B. Schmitt, C. Schnabel, Z. Schoenfeld, T. L. Scholten, E. Schoute, J. Schuhmacher, M. Schulterbrandt, J. Schwarm, P. Schweigert, J. Seaward, Sergi, D. E. Serrano, I. F. Sertage, K. Setia, F. Shah, N. Shammah, W. Shanks, R. Sharma, P. Shaw, Y. Shi, J. Shoemaker, A. Silva, A. Simonetto, D. Singh, D. Singh, P. Singh, P. Singkanipa, Y. Siraichi, Siri, J. Sistos, J. Sistos, I. Sitdikov, S. Sivarajah, Slavikmew, M. B. Sletfjerding, J. A. Smolin, M. Soeken, I. O. Sokolov, I. Sokolov, V. P. Soloviev, SooluThomas, Starfish, D. Steenken, M. Stypulkoski, A. Suau, S. Sun, K. J. Sung, M. Suwama, O. Słowik, R. Taeja, H. Takahashi, T. Takawale, I. Tavernelli, C. Taylor, P. Taylour, S. Thomas, K. Tian, M. Tillet, M. Tod, M. Tomasik, C. Tornow, E. de la Torre, J. L. S. Toural, K. Trabing, M. Treinish, D. Trenev, TrishaPe, F. Truger, TsafrirA, G. Tsilimigkounakis, K. Tsuoka, D. Tulsi, D. Tuna, W. Turner, Y. Vaknin, C. R. Valcarce, F. Varchon, A. Vartak, A. C. Vazquez, P. Vijaywargiya, V. Villar, B. Vishnu, D. Vogt-Lee, C. Vuillot, WQ, J. Weaver, J. Weidenfeller, R. Wieczorek, J. A. Wildstrom, J. Wilson, E. Winston, WinterSoldier, J. J. Woehr, S. Woerner, R. Woo, C. J. Wood, R. Wood, S. Wood, J. Wootton, M. Wright, L. Xing, J. YU, Yaiza, B. Yang, U. Yang, J. Yao, D. Yeralin, R. Yonekura, D. Yonge-Mallo, R. Yoshida, R. Young, J. Yu, L. Yu, Yuma-Nakamura, C. Zachow, L. Zdanski, H. Zhang, E. Zheltonozhskii, I. Zidaru, B. Zimmermann, B. Zindorf, C. Zoufal, a matsuo, aeddins ibm, alexzhang13, b63, bartek bartlomiej, bcamorrison, brandhsn, nick bronn, chetmurthy, choerst ibm, comet, dalin27, deeplokhande, dekel.meirom, derwind, dime10, ehchen, ewinston, fanizzamarco, fs1132429, gadial, galeinston, georgezhou20, georgios ts, gruu, hhorii, hhyap, hykavitha, itoko, jeppevinkel, jessica angel7, jezerjojo14, jliu45, johannesgreiner, jscott2, kUmezawa, klinvill, krutik2966, ma5x, michelle4654, msuwama, nico lgrs, nrhawkins, ntgiwsvp, ordmoj, sagar pahwa, pritamsinha2304, rithikaadiga, ryancocuzzo, saktar unr, saswati qiskit, sebastian mair, septembrr, sethmerkel, sg495, shaashwat, smturro2, sternparky, strickroman, tigerjack, tsura crisaldo, upsideon, vadebayo49, welien, willhbang, wmurphy collabstar, yang.luh, yuri@FreeBSD, and M. Čepulkovskis, “Qiskit: An open-source framework for quantum computing,”  (2021).
  * Jordan and Wigner (1928) P. Jordan and E. Wigner, Z. Phys. 47, 631 (1928).
  * Huggins  _et al._ (2020) W. J. Huggins, J. Lee, U. Baek, B. O’Gorman, and K. B. Whaley, New J. Phys. (2020).
  * Stair and Evangelista (2021) N. H. Stair and F. A. Evangelista, PRX Quantum 2, 030301 (2021).
  * Cai  _et al._ (2022) Z. Cai, R. Babbush, S. C. Benjamin, S. Endo, W. J. Huggins, Y. Li, J. R. McClean, and T. E. O’Brien, arXiv:2210.00921 (2022).
  * Grimsley  _et al._ (2019) H. R. Grimsley, D. Claudino, S. E. Economou, E. Barnes, and N. J. Mayhall, J. Comp. Theor. Chem. 16, 1 (2019).
  * McClean  _et al._ (2017) J. R. McClean, M. E. Kimchi-Schwartz, J. Carter, and W. A. de Jong, Phys. Rev. A 95, 042308 (2017).
  * Motta  _et al._ (2020) M. Motta, C. Sun, A. T. Tan, M. J. O’Rourke, E. Ye, A. J. Minnich, F. G. Brandão, and G. K.-L. Chan, Nat. Phys. 16, 205 (2020).
  * Stair  _et al._ (2020) N. H. Stair, R. Huang, and F. A. Evangelista, J. Chem. Theory Comput. 16, 2236 (2020).