#  Spatial, spin, and charge symmetry projections for a Fermi-Hubbard model on
a quantum computer

Kazuhiro Seki  Quantum Computational Science Research Team, RIKEN Center for
Quantum Computing (RQC), Saitama 351-0198, Japan  Seiji Yunoki  Quantum
Computational Science Research Team, RIKEN Center for Quantum Computing (RQC),
Saitama 351-0198, Japan  Computational Quantum Matter Research Team, RIKEN
Center for Emergent Matter Science (CEMS), Saitama 351-0198, Japan
Computational Materials Science Research Team, RIKEN Center for Computational
Science (R-CCS), Hyogo 650-0047, Japan  Computational Condensed Matter Physics
Laboratory, RIKEN Cluster for Pioneering Research (CPR), Saitama 351-0198,
Japan

###### Abstract

We propose an extended version of the symmetry-adapted variational-quantum-
eigensolver (VQE) and apply it to a two-component Fermi-Hubbard model on a
bipartite lattice. In the extended symmetry-adapted VQE method, the Rayleigh
quotient for the Hamiltonian and a parametrized quantum state in a properly
chosen subspace is minimized within the subspace and is optimized among the
variational parameters implemented on a quantum circuit to obtain
variationally the ground state and the ground-state energy. The corresponding
energy derivative with respect to a variational parameter is expressed as a
Hellmann-Feynman-type formula of a generalized eigenvalue problem in the
subspace, which thus allows us to use the parameter-shift rules for its
evaluation. The natural-gradient-descent method is also generalized to
optimize variational parameters in a quantum-subspace-expansion approach. As a
subspace for approximating the ground state of the Hamiltonian, we consider a
Krylov subspace generated by the Hamiltonian and a symmetry-projected
variational state, and therefore the approximated ground state can restore the
Hamiltonian symmetry that is broken in the parametrized variational state
prepared on a quantum circuit. We show that spatial symmetry operations for
fermions in an occupation basis can be expressed as a product of the nearest-
neighbor fermionic swap operations on a quantum circuit. We also describe how
the spin and charge symmetry operations, i.e., rotations, can be implemented
on a quantum circuit. By numerical simulations, we demonstrate that the
spatial, spin, and charge symmetry projections can improve the accuracy of the
parametrized variational state, which can be further improved systematically
by expanding the Krylov subspace without increasing the number of variational
parameters.

##  I Introduction

Recent technological advances in quantum devices Nakamura  _et al._ (1999);
Kok  _et al._ (2007); Ladd  _et al._ (2010); Xiang  _et al._ (2013); Chow  _et
al._ (2014); Barends  _et al._ (2014); Ristè  _et al._ (2015); Kelly  _et al._
(2015); Arute  _et al._ (2019); Zhong  _et al._ (2020); Blais  _et al._ (2021)
have suggested that quantum computation of quantum physics and chemistry
systems Feynman (1982) is becoming a reality in the not-so-distant future
Arute  _et al._ (2020); Stanisic  _et al._ (2021). Currently available quantum
computers are, however, prone to noise and hence the size of a quantum circuit
to be reliably executed is limited. Such quantum computers are called noisy
intermediate-scale quantum (NISQ) computers Preskill (2018). Despite the
limitation, NISQ computers with 50-100 qubits are anticipated to show
advantage over classical computers with the best known algorithm for
particular tasks (for example, see Refs. Arute  _et al._ (2019); Zhong  _et
al._ (2020); Huang  _et al._ (2020)). Thus, in parallel with research to
accomplish fault-tolerant quantum computers, for which an increasing number of
experimental developments towards realization of logical qubit operations has
been reported recently Nigg  _et al._ (2014); Erhard  _et al._ (2021);
Satzinger  _et al._ (2021); Marques  _et al._ (2021); Hilder  _et al._ (2021),
it is of importance to find practical applications of NISQ computers for
further stimulating progress in the field of quantum computing. To this and,
several quantum-classical hybrid algorithms and error mitigation schemes have
been developed Li and Benjamin (2017); Endo  _et al._ (2018); Bonet-Monroig
_et al._ (2018); Endo  _et al._ (2019).

The variational-quantum eigensolver (VQE) Yung  _et al._ (2014); Peruzzo  _et
al._ (2014); McClean  _et al._ (2016); Kandala  _et al._ (2017) is one of the
potentially promising quantum-classical hybrid schemes for solving eigenvalue
problems in quantum chemistry and quantum physics with NISQ computers. For
recent reviews on variational quantum algorithms, see for example Refs. Endo
_et al._ (2021); Cerezo  _et al._ (2021); Tilly  _et al._ (2021). In the VQE,
the expectation value of a Hamiltonian of interest with respect to a
variational state represented on a parametrized quantum circuit, i.e., the
variational energy, is evaluated with a quantum computer, while variational
parameters are optimized by minimizing the variational energy on a classical
computer. Depending on the form of variational states and the assignment of
tasks for quantum and classical computers, several variants of the VQE scheme
have been proposed. For example, VQE-type approaches based on the quantum-
subspace expansion (QSE) Colless  _et al._ (2018) performs subspace
diagonalization in an appropriately chosen subspace to approximate the target
state(s) (e.g. the ground state) better than the conventional VQE scheme, at a
cost of polynomial numbers of additional measurements Nakanishi  _et al._
(2019); Heya  _et al._ (2019). Aiming at a systematic construction of the
subspace, a Krylov subspace Chatelin (2012) generated by properly chosen
initial states with a real-time evolution operator Parrish and McMahon (2019);
Huggins  _et al._ (2020); Stair  _et al._ (2020), an imaginary-time evolution
operator Motta  _et al._ (2019); Yeter-Aydeniz  _et al._ (2020), or a
Hamiltonian Seki and Yunoki (2021) is often adopted.

In addition, to exploit the Hamiltonian symmetry, several variants or
extensions of the VQE approach have been proposed. For example, we have
proposed the symmetry-adapted VQE (SAVQE) method to encompass the Hamiltonian
symmetry in the VQE scheme Seki  _et al._ (2020). In the SAVQE method, the
symmetry that is broken in the variational state generated on a parametrized
quantum circuit is restored with a symmetry-projection operator. The
nonunitarity of the projection operator is treated classically as a post
processing with increasing the number of measurements, a similar idea for
treating Jastrow-type correlators reported earlier in Ref. Mazzola  _et al._
(2019). Moreover, related schemes have been employed to restore symmetry in a
continuous group such as the SU(2) total-spin conservation Tsuchimochi  _et
al._ (2020); Siwach and Lacroix (2021) and the U(1) particle-number
conservation Khamoshi  _et al._ (2020); Guzman and Lacroix (2021), where the
integral over continuous parameters of the group in the projection operator is
properly discretized. Concerning the translational symmetry of periodic
systems, another VQE-type scheme has been proposed to formulate directly in
the reciprocal space Liu  _et al._ (2020); Manrique  _et al._ (2020); Yoshioka
_et al._ (2022). Furthermore, in a different VQE variant, an appropriate
penalty term is introduced into a cost function to obtain an eigenstate of the
Hamiltonian in a desired symmetry sector Ryabinkin  _et al._ (2019); Kuroiwa
and Nakagawa (2021). It should also be noted that the Hamiltonian symmetry can
also be utilized to mitigate errors due to different noise channels Bonet-
Monroig  _et al._ (2018). Remarkably, recent experiments have demonstrated a
significant improvement for mitigating errors in a VQE simulation of a Fermi-
Hubbard model by an error mitigation technique based on the Hamiltonian
symmetry, including spin- and particle-number conservations, time-reversal
symmetry, particle-hole symmetry, and spatial symmetry Stanisic  _et al._
(2021).

In this paper, we propose a QSE-based VQE method that incorporates the
symmetry of Hamiltonian. The main idea is based on the SAVQE method Seki  _et
al._ (2020) and the quantum power method (QPM) Seki and Yunoki (2021), both of
which have been developed recently by the present authors. In the proposed
method, we construct a Krylov subspace by multiplying Hamiltonian power onto a
quantum state that is obtained by applying the symmetry-projection operators
to a single variational state prepared on a parametrized quantum circuit. We
then perform a subspace diagonalization by minimizing the Rayleigh quotient
for the Hamiltonian and the quantum states in the Krylov subspace, and obtain
the lowest eigenvalue of the generalized eigenvalue problem that depends on
the variational parameters. The variational parameters are optimized so as to
minimize the lowest eigenvalue of the generalized eigenvalue problem, which
gives us the variational ground-state energy with the optimal set of
variational parameters that represents the variational ground state.
Considering a two-component Fermi-Hubbard model in a ladder lattice structure,
we numerically demonstrate the proposed method by showing that the estimated
ground-state energy as well as the ground-state fidelity can be improved by
the spatial, spin, and charge symmetry projections, and they are further
improved systematically with expanding the subspace without increasing the
number of variational parameters.

The rest of this paper is organized as follows. In Sec. II, we briefly
summarize a general formalism of the QSE method that is relevant for our
purpose. In Sec. III, we formulate the natural-gradient-descent (NGD) method
for the QSE scheme, and we describe how the corresponding energy gradient and
the Fubini-Study metric tensor can be evaluated with a quantum computer using
the parameter-shift rule. We then describe in Sec. IV the Krylov-extended
SAVQE by introducing the Krylov subspace generated by the Hamiltonian and a
symmetry-projected variational state. In Sec. V, we define the Fermi-Hubbard
model on a bipartite lattice, and we briefly review its spatial, spin, and
charge symmetry. We then describe the corresponding symmetry-projection
operators. In Sec. VI, we first show how the symmetry operations for a fermion
model in general can be implemented on a quantum circuit, and then we explain
the case of the spin and charge symmetry operations. In Sec. VII, we
numerically demonstrate the proposed method for the Fermi-Hubbard model. A
conclusion and a discussion are given in Sec. VIII. Additional details on the
NGD method and the fermionic symmetry operations are provided in Appendixes A
and B, respectively. Matrix representations of typical two-qubit two-level
unitary gates for quantum many-body systems, such as the Givens-rotation gate
and the Bogoliubov-transformation gate, are provided in Appendix C. A simple
parallelization scheme for classical simulation of the VQE method is described
in Appendix D. A remark on the normalization factor of the symmetry-projected
quantum state is made in Appendix E. Further numerical results for a different
type of variational states are discussed in Appendix F.

##  II Formalism of quantum-subspace expansion

###  II.1 Trial state and variational principle

Consider a subspace

${\cal U}={\rm span}\left(|u_{0}\rangle,|u_{1}\rangle,\cdots,|u_{d_{\cal
U}-1}\rangle\right)$ (1)

with $\dim{\cal U}\equiv d_{\cal U}$. The basis states
$\\{|u_{i}\rangle\\}_{i=0}^{d_{\cal U}-1}$ should be linearly independent of
each other but they are not necessarily orthonormalized. We assume that
$|u_{i}\rangle$ has a form of

$|u_{i}(\boldsymbol{\theta})\rangle=\hat{O}_{i}|\psi(\boldsymbol{\theta})\rangle,$
(2)

where $|\psi(\boldsymbol{\theta})\rangle$ is a variational state (ansatz)
parametrized by a set of $N_{\rm v}$ variational parameters
$\boldsymbol{\theta}=\\{\theta_{k}\\}_{k=1}^{N_{\rm v}}$, assuming they are
real, and $\hat{O}_{i}$ is an operator independent of the variational
parameters. $\hat{O}_{i}$ is not necessarily unitary but it is assumed to be
given as a linear combination of unitary operators. Note that by definition,
the parametrized part $|\psi(\boldsymbol{\theta})\rangle$ is common to all
$\\{|u_{i}(\boldsymbol{\theta})\rangle\\}_{i=0}^{d_{\cal U}-1}$ in Eq. (2).

We now intend to approximate the exact ground state $|\Psi_{0}\rangle$ of the
Hamiltonian $\hat{\cal H}$ within the subspace ${\cal U}$. This can be done
simply by assuming that a trial state $|\Psi_{\cal
U}(\boldsymbol{\theta})\rangle$ for approximating the ground state
$|\Psi_{0}\rangle$ is given by

$|\Psi_{\cal U}(\boldsymbol{\theta})\rangle=\sum_{i=0}^{d_{\cal
U}-1}v_{i}(\boldsymbol{\theta})|u_{i}(\boldsymbol{\theta})\rangle$ (3)

in the subspace ${\cal U}$, where
$\boldsymbol{v}(\boldsymbol{\theta})=\\{v_{i}(\boldsymbol{\theta})\\}_{i=0}^{d_{\cal
U}-1}$ are coefficients to be determined. According to the variational
principle, the optimal variational parameters $\boldsymbol{\theta}_{\rm opt}$
are obtained by minimizing the expectation value of the Hamiltonian $\hat{\cal
H}$, i.e., the variational energy, with respect to both $\boldsymbol{\theta}$
and $\boldsymbol{v}(\boldsymbol{\theta})$:

$\displaystyle\boldsymbol{\theta}_{\rm
opt}\equiv\underset{\boldsymbol{\theta}}{\rm
arg~{}min}\left\\{\min_{\boldsymbol{v}(\boldsymbol{\theta})}\left\\{\frac{\langle\Psi_{\cal
U}(\boldsymbol{\theta})|\hat{\cal H}|\Psi_{\cal
U}(\boldsymbol{\theta})\rangle}{\langle\Psi_{\cal
U}(\boldsymbol{\theta})|\Psi_{\cal
U}(\boldsymbol{\theta})\rangle}\right\\}\right\\}.$ (4)

###  II.2 Subspace diagonalization

Let us consider the optimization with respect to
$\boldsymbol{v}(\boldsymbol{\theta})$ for a given set of variational
parameters $\boldsymbol{\theta}$ in Eq. (4) and denote the optimal
coefficients as
$\boldsymbol{v}_{0}(\boldsymbol{\theta})\equiv\\{v_{i,0}(\boldsymbol{\theta})\\}_{i=0}^{d_{\cal
U}-1}$, i.e.,

$\boldsymbol{v}_{0}(\boldsymbol{\theta})=\underset{\boldsymbol{v}(\boldsymbol{\theta})}{\rm
arg~{}min}\left\\{\frac{\langle\Psi_{\cal U}(\boldsymbol{\theta})|\hat{\cal
H}|\Psi_{\cal U}(\boldsymbol{\theta})\rangle}{\langle\Psi_{\cal
U}(\boldsymbol{\theta})|\Psi_{\cal U}(\boldsymbol{\theta})\rangle}\right\\}.$
(5)

The optimal coefficients $\boldsymbol{v}_{0}(\boldsymbol{\theta})$ subject to
a normalization condition, e.g., $\langle\Psi_{\cal
U}(\boldsymbol{\theta})|\Psi_{\cal U}(\boldsymbol{\theta})\rangle=1$, can be
determined by solving the generalized eigenvalue problem Seki and Yunoki
(2021)

$\boldsymbol{H}(\boldsymbol{\theta})\boldsymbol{v}_{0}(\boldsymbol{\theta})=E_{0}(\boldsymbol{\theta})\boldsymbol{S}(\boldsymbol{\theta})\boldsymbol{v}_{0}(\boldsymbol{\theta}),$
(6)

where $\boldsymbol{H}(\boldsymbol{\theta})$ is the Hamiltonian matrix with its
element

$[\boldsymbol{H}(\boldsymbol{\theta})]_{ij}=\langle
u_{i}(\boldsymbol{\theta})|\hat{\cal
H}|u_{j}(\boldsymbol{\theta})\rangle=\langle\psi(\boldsymbol{\theta})|\hat{O}_{i}^{\dagger}\hat{\cal
H}\hat{O}_{j}|\psi(\boldsymbol{\theta})\rangle,$ (7)

$\boldsymbol{S}(\boldsymbol{\theta})$ is the overlap matrix with its element

$[\boldsymbol{S}(\boldsymbol{\theta})]_{ij}=\langle
u_{i}(\boldsymbol{\theta})|u_{j}(\boldsymbol{\theta})\rangle=\langle\psi(\boldsymbol{\theta})|\hat{O}_{i}^{\dagger}\hat{O}_{j}|\psi(\boldsymbol{\theta})\rangle,$
(8)

and $E_{0}(\boldsymbol{\theta})$ is the smallest eigenvalue. The trial state

$|\Psi_{\cal U}^{(0)}(\boldsymbol{\theta})\rangle\equiv\sum_{i=0}^{d_{\cal
U}-1}v_{i,0}(\boldsymbol{\theta})|u_{i}(\boldsymbol{\theta})\rangle$ (9)

with $[\boldsymbol{v}_{0}]_{i}=v_{i,0}$ being the $i$th entry of the
corresponding eigenvector in Eq. (6) is normalized as long as the eigenvector
$\boldsymbol{v}_{0}(\boldsymbol{\theta})$ is normalized with respect to
$\boldsymbol{S}(\boldsymbol{\theta})$, i.e.,

$\langle\Psi_{\cal U}^{(0)}(\boldsymbol{\theta})|\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta})\rangle=\boldsymbol{v}_{0}^{\dagger}(\boldsymbol{\theta})\boldsymbol{S}(\boldsymbol{\theta})\boldsymbol{v}_{0}(\boldsymbol{\theta})=1,$
(10)

where
$\boldsymbol{v}_{0}^{\dagger}(\boldsymbol{\theta})=[\boldsymbol{v}_{0}^{\rm
T}(\boldsymbol{\theta})]^{*}$.

###  II.3 Energy and other expectation values

The smallest eigenvalue $E_{0}(\boldsymbol{\theta})$ of the generalized
eigenvalue problem in Eq. (6) corresponds to the minimum value of the Rayleigh
quotient, i.e.,

$E_{0}(\boldsymbol{\theta})=\min_{\boldsymbol{v}(\boldsymbol{\theta})}\left\\{\frac{\langle\Psi_{\cal
U}(\boldsymbol{\theta})|\hat{\cal H}|\Psi_{\cal
U}(\boldsymbol{\theta})\rangle}{\langle\Psi_{\cal
U}(\boldsymbol{\theta})|\Psi_{\cal
U}(\boldsymbol{\theta})\rangle}\right\\}=\boldsymbol{v}_{0}^{\dagger}(\boldsymbol{\theta})\boldsymbol{H}(\boldsymbol{\theta})\boldsymbol{v}_{0}(\boldsymbol{\theta}),$
(11)

where the normalization condition for $|\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta})\rangle$ in Eq. (10) is used in the second
equality. Therefore, the approximated ground-state energy of the Hamiltonian
$\hat{\cal H}$ is obtained by minimizing $E_{0}(\boldsymbol{\theta})$ with
respect to the variational parameters $\boldsymbol{\theta}$,

$\displaystyle E_{0}(\boldsymbol{\theta}_{\rm opt})$
$\displaystyle=\min_{\boldsymbol{\theta}}\left\\{E_{0}(\boldsymbol{\theta})\right\\}=\min_{\boldsymbol{\theta}}\left\\{\frac{\langle\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta})|\hat{\cal H}|\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta})\rangle}{\langle\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta})|\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta})\rangle}\right\\}$
$\displaystyle=\min_{\boldsymbol{\theta}}\left\\{\min_{\boldsymbol{v}(\boldsymbol{\theta})}\left\\{\frac{\langle\Psi_{\cal
U}(\boldsymbol{\theta})|\hat{\cal H}|\Psi_{\cal
U}(\boldsymbol{\theta})\rangle}{\langle\Psi_{\cal
U}(\boldsymbol{\theta})|\Psi_{\cal
U}(\boldsymbol{\theta})\rangle}\right\\}\right\\},$ (12)

and the optimized variational parameters $\boldsymbol{\theta}_{\rm opt}$ give
us the approximated ground state $|\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta}_{\rm opt})\rangle$. Once a set of optimal
variational parameters $\boldsymbol{\theta}_{\rm opt}$ is obtained, the
ground-state expectation value of an observable $\hat{A}$ can be approximated
as

$\langle\Psi_{0}|\hat{A}|\Psi_{0}\rangle\approx\langle\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta}_{\rm opt})|\hat{A}|\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta}_{\rm
opt})\rangle=\boldsymbol{v}_{0}^{\dagger}(\boldsymbol{\theta}_{\rm
opt})\boldsymbol{A}(\boldsymbol{\theta}_{\rm
opt})\boldsymbol{v}_{0}(\boldsymbol{\theta}_{\rm opt}),$ (13)

where $[\boldsymbol{A}(\boldsymbol{\theta})]_{ij}=\langle
u_{i}(\boldsymbol{\theta})|\hat{A}|u_{j}(\boldsymbol{\theta})\rangle=\langle\psi(\boldsymbol{\theta})|\hat{O}_{i}^{\dagger}\hat{A}\hat{O}_{j}|\psi(\boldsymbol{\theta})\rangle$
and the normalization condition for $|\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta}_{\rm opt})\rangle$ in Eq. (10) is assumed.

##  III Parameter-optimization method

###  III.1 Natural-gradient-descent method

The NGD method optimizes $\boldsymbol{\theta}$ by minimizing
$E_{0}(\boldsymbol{\theta})$ with the following iteration Amari (1998):

$\boldsymbol{\theta}^{(x+1)}=\boldsymbol{\theta}^{(x)}-\tau[\boldsymbol{G}(\boldsymbol{\theta}^{(x)})]^{-1}\nabla
E_{0}(\boldsymbol{\theta}^{(x)}),$ (14)

where $\boldsymbol{\theta}^{(x)}=\\{\theta_{k}^{(x)}\\}_{k=1}^{N_{\rm v}}$ are
the variational parameters at the $x$th iteration and $\tau>0$ is a parameter
called learning rate. $\boldsymbol{G}(\boldsymbol{\theta})$ is the Fubini-
Study metric tensor Provost and Vallee (1980); Kolodrubetz  _et al._ (2017);
Stokes  _et al._ (2020) defined below and $\left[\nabla
E_{0}(\boldsymbol{\theta})\right]_{k}=\partial
E_{0}(\boldsymbol{\theta})/\partial\theta_{k}$ is the energy gradient with
respect to the variational parameter $\theta_{k}$. The gradient descent (GD)
method can be obtained by setting $\boldsymbol{G}=\boldsymbol{I}$ in Eq. (14).
The initial variational parameters $\boldsymbol{\theta}^{(1)}$ can be set
arbitrary. In the following, we shall derive explicit forms of the energy
gradient and the Fubini-Study metric tensor for a QSE-based approach.

From Eqs. (6) and (10) and the Hermiticity of
$\boldsymbol{H}(\boldsymbol{\theta})$ and
$\boldsymbol{S}(\boldsymbol{\theta})$, the energy derivative is given by

$\frac{\partial
E_{0}(\boldsymbol{\theta})}{\partial\theta_{k}}=\boldsymbol{v}_{0}^{\dagger}(\boldsymbol{\theta})\left(\frac{\partial\boldsymbol{H}(\boldsymbol{\theta})}{\partial\theta_{k}}-E_{0}(\boldsymbol{\theta})\frac{\partial\boldsymbol{S}(\boldsymbol{\theta})}{\partial\theta_{k}}\right)\boldsymbol{v}_{0}(\boldsymbol{\theta}).$
(15)

Equation (15) can be considered as a variant of the Hellmann-Feynman theorem
for a generalized eigenvalue problem in the sense that the energy derivative
does not require the derivative of the eigenvector
$\boldsymbol{v}_{0}(\boldsymbol{\theta})$. The difference from the Hellmann-
Feynman theorem for a standard eigenvalue problem is that Eq. (15) involves
the derivative of the overlap matrix $\boldsymbol{S}(\boldsymbol{\theta})$.

The Fubini-Study metric tensor $\boldsymbol{G}(\boldsymbol{\theta})$ is given
by

$[\boldsymbol{G}(\boldsymbol{\theta})]_{kl}={\rm
Re}\left[\gamma_{kl}(\boldsymbol{\theta})-\beta_{k}^{*}(\boldsymbol{\theta})\beta_{l}(\boldsymbol{\theta})\right],$
(16)

where

$\gamma_{kl}(\boldsymbol{\theta})=\langle\partial_{k}\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta})|\partial_{l}\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta})\rangle$ (17)

and

$\beta_{k}(\boldsymbol{\theta})=\langle\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta})|\partial_{k}\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta})\rangle$ (18)

with $\partial_{k}\equiv\partial/\partial{\theta_{k}}$ (see Appendix A for
details). Note that our definition of $\beta_{k}$ in Eq. (18) differs from
that in Ref. Provost and Vallee (1980) by a multiplicative factor $-{\rm i}$.
The derivative of the normalization condition $\partial_{k}\langle\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta})|\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta})\rangle=0$ implies that
$\beta_{k}(\boldsymbol{\theta})$ is pure imaginary and hence
$\boldsymbol{G}(\boldsymbol{\theta})$ is a real symmetric matrix. By
substituting

$\displaystyle|\partial_{k}\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta})\rangle=\sum_{i}\partial_{k}v_{i,0}(\boldsymbol{\theta})|u_{i}(\boldsymbol{\theta})\rangle+\sum_{i}v_{i,0}(\boldsymbol{\theta})|{\partial_{k}}u_{i}(\boldsymbol{\theta})\rangle$
(19)

into Eqs. (17) and (18), we obtain that

$\displaystyle\gamma_{kl}(\boldsymbol{\theta})$
$\displaystyle=\sum_{ij}(\partial_{k}v_{i,0}^{*})(\partial_{l}v_{j,0})\langle
u_{i}|u_{j}\rangle+\sum_{ij}v_{i,0}^{*}v_{j,0}\langle\partial_{k}u_{i}|\partial_{l}u_{j}\rangle$
$\displaystyle+\sum_{ij}(\partial_{k}v_{i,0}^{*})v_{j,0}\langle
u_{i}|\partial_{l}u_{j}\rangle+\sum_{ij}v_{i,0}^{*}(\partial_{l}v_{j,0})\langle\partial_{k}u_{i}|u_{j}\rangle$
$\displaystyle=(\partial_{k}\boldsymbol{v}_{0})^{\dagger}\boldsymbol{S}(\partial_{l}\boldsymbol{v}_{0})+\boldsymbol{v}_{0}^{\dagger}\boldsymbol{S}_{(k,l)}\boldsymbol{v}_{0}$
$\displaystyle+(\partial_{k}\boldsymbol{v}_{0})^{\dagger}\boldsymbol{S}_{(,l)}\boldsymbol{v}_{0}+\boldsymbol{v}_{0}^{\dagger}[\boldsymbol{S}_{(,k)}]^{\dagger}(\partial_{l}\boldsymbol{v}_{0})$
(20)

and

$\displaystyle\beta_{k}(\boldsymbol{\theta})$
$\displaystyle=\sum_{ij}v_{i,0}^{*}(\partial_{k}v_{j,0})\langle
u_{i}|u_{j}\rangle+\sum_{ij}v_{i,0}^{*}v_{j,0}\langle
u_{i}|\partial_{k}u_{j}\rangle$
$\displaystyle=\boldsymbol{v}_{0}^{\dagger}\boldsymbol{S}(\partial_{k}\boldsymbol{v}_{0})+\boldsymbol{v}_{0}^{\dagger}\boldsymbol{S}_{(,k)}\boldsymbol{v}_{0},$
(21)

where $[\boldsymbol{S}_{(,k)}(\boldsymbol{\theta})]_{ij}\equiv\langle
u_{i}(\boldsymbol{\theta})|\partial_{k}u_{j}(\boldsymbol{\theta})\rangle$ and
$[\boldsymbol{S}_{(k,l)}(\boldsymbol{\theta})]_{ij}\equiv\langle\partial_{k}u_{i}(\boldsymbol{\theta})|\partial_{l}u_{j}(\boldsymbol{\theta})\rangle$.
Note that the $\boldsymbol{\theta}$ dependence of $|u_{i}\rangle$,
$\boldsymbol{v}_{0}$, and $\boldsymbol{S}$ is implicitly assumed in the right-
hand sides of Eqs. (20) and (21). We should also note that
$[\boldsymbol{S}_{(k,l)}(\boldsymbol{\theta})]^{\dagger}=\boldsymbol{S}_{(l,k)}(\boldsymbol{\theta})$
and hence not all matrices $\boldsymbol{S}_{(k,l)}(\boldsymbol{\theta})$ are
independent.

To obtain the derivative of the eigenvector
$\boldsymbol{v}_{0}(\boldsymbol{\theta})$ with respect to the variational
parameters $\boldsymbol{\theta}$, let us first expand it in the subspace
${\cal U}$ as

$\frac{\partial\boldsymbol{v}_{0}(\boldsymbol{\theta})}{\partial\theta_{k}}=\sum_{n=0}^{d_{\cal
U}-1}c_{k,n}\boldsymbol{v}_{n}(\boldsymbol{\theta}),$ (22)

where $\boldsymbol{v}_{n}$ is the $n$th eigenvector satisfying
$\boldsymbol{H}(\boldsymbol{\theta})\boldsymbol{v}_{n}(\boldsymbol{\theta})=E_{n}(\boldsymbol{\theta})\boldsymbol{S}(\boldsymbol{\theta})\boldsymbol{v}_{n}(\boldsymbol{\theta})$
with
$\boldsymbol{v}_{m}^{\dagger}(\boldsymbol{\theta})\boldsymbol{S}(\boldsymbol{\theta})\boldsymbol{v}_{n}(\boldsymbol{\theta})=\delta_{mn}$
for $0\leqslant m,n\leqslant d_{\cal U}-1$, and $\\{c_{k,n}\\}_{n=0}^{d_{\cal
U}-1}$ are complex numbers to be determined. By taking the derivative of
$\boldsymbol{H}(\boldsymbol{\theta})\boldsymbol{v}_{n}(\boldsymbol{\theta})=E_{n}(\boldsymbol{\theta})\boldsymbol{S}(\boldsymbol{\theta})\boldsymbol{v}_{n}(\boldsymbol{\theta})$,
multiplying $\boldsymbol{v}_{n}^{\dagger}(\boldsymbol{\theta})$ from the left
on it, and using
$\boldsymbol{v}_{m}^{\dagger}(\boldsymbol{\theta})\boldsymbol{S}(\boldsymbol{\theta})\boldsymbol{v}_{n}(\boldsymbol{\theta})=\delta_{mn}$
for $n\not=0$, we find that
$c_{k,n}=\frac{\boldsymbol{v}_{n}^{\dagger}\left(\partial_{k}\boldsymbol{H}-E_{0}\partial_{k}\boldsymbol{S}\right)\boldsymbol{v}_{0}}{E_{0}-E_{n}}\boldsymbol{v}_{n}$
for $n\not=0$ and hence

$\frac{\partial\boldsymbol{v}_{0}(\boldsymbol{\theta})}{\partial\theta_{k}}=c_{k,0}\boldsymbol{v}_{0}(\boldsymbol{\theta})+\sum_{n\not=0}^{d_{\cal
U}-1}\frac{\boldsymbol{v}_{n}^{\dagger}(\boldsymbol{\theta})\left(\tfrac{\partial\boldsymbol{H}(\boldsymbol{\theta})}{\partial\theta_{k}}-E_{0}(\boldsymbol{\theta})\tfrac{\partial\boldsymbol{S}(\boldsymbol{\theta})}{\partial\theta_{k}}\right)\boldsymbol{v}_{0}(\boldsymbol{\theta})}{E_{0}(\boldsymbol{\theta})-E_{n}(\boldsymbol{\theta})}\boldsymbol{v}_{n}(\boldsymbol{\theta}),$
(23)

where $E_{0}(\boldsymbol{\theta})\not=E_{n}(\boldsymbol{\theta})$ is assumed.
For the remaining coefficient $c_{k,0}$, the derivative of the normalization
condition $\partial_{k}\langle\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta})|\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta})\rangle=\partial_{k}(\boldsymbol{v}_{0}^{\dagger}(\boldsymbol{\theta})\boldsymbol{S}(\boldsymbol{\theta})\boldsymbol{v}_{0}(\boldsymbol{\theta}))=0$
implies that

$\displaystyle{\rm
Re}c_{k,0}=-\frac{1}{2}\boldsymbol{v}_{0}^{{\dagger}}(\boldsymbol{\theta})\frac{\partial\boldsymbol{S}(\boldsymbol{\theta})}{\partial\theta_{k}}\boldsymbol{v}_{0}(\boldsymbol{\theta}),$
(24) $\displaystyle{\rm Im}c_{k,0}\text{:\,undetermined}.$ (25)

As we shall discuss below, the coefficient $c_{k,0}$ can be chosen arbitrarily
as far as the Fubini-Study metric tensor $\boldsymbol{G}(\boldsymbol{\theta})$
is concerned. For numerical simulations in Sec. VII, we set that ${\rm
Im}c_{k,0}=0$.

Here are two remarks regarding the coefficient $c_{k,0}$. First, ${\rm
Re}c_{k,0}$ in Eq. (24) ensures that $\beta_{k}(\boldsymbol{\theta})$ is pure
imaginary [see the discussion below Eq. (18)], which follows from Eq. (21)
because one can easily show that
$\boldsymbol{v}_{0}^{\dagger}\boldsymbol{S}(\partial_{k}\boldsymbol{v}_{0})=c_{k,0}$
and ${\rm
Re}\left[\boldsymbol{v}_{0}^{\dagger}\boldsymbol{S}_{(,k)}\boldsymbol{v}_{0}\right]=\frac{1}{2}\boldsymbol{v}_{0}^{{\dagger}}(\partial_{k}\boldsymbol{S})\boldsymbol{v}_{0}$.
This is a natural generalization of the corresponding result for real
symmetric $\boldsymbol{H}(\boldsymbol{\theta})$ and
$\boldsymbol{S}(\boldsymbol{\theta})$ in a real inner product space reported
in Ref. Ruymbeek  _et al._ (2019) to a Hermitian case in a complex inner
product space. Second, the Fubini-Study metric tensor
$\boldsymbol{G}(\boldsymbol{\theta})$ is independent of $c_{k,0}$. To prove
this statement, let us consider a shift of the coefficient $c_{k,0}\mapsto
c_{k,0}+\Delta c_{k,0}$, i.e.,

$\frac{\partial\boldsymbol{v}_{0}(\boldsymbol{\theta})}{\partial\theta_{k}}\mapsto\frac{\partial\boldsymbol{v}_{0}(\boldsymbol{\theta})}{\partial\theta_{k}}+\Delta
c_{k,0}\boldsymbol{v}_{0}(\boldsymbol{\theta}).$ (26)

Then, $\gamma_{kl}(\boldsymbol{\theta})$ and
$\beta_{k}^{*}(\boldsymbol{\theta})\beta_{l}(\boldsymbol{\theta})$ are
transformed accordingly as

$\displaystyle\gamma_{kl}(\boldsymbol{\theta})$
$\displaystyle\mapsto\gamma_{kl}(\boldsymbol{\theta})+\Delta c_{k,0}^{*}\Delta
c_{l,0}+\Delta c_{k,0}^{*}\beta_{l}(\boldsymbol{\theta})+\Delta
c_{l,0}\beta_{k}^{*},(\boldsymbol{\theta})$ (27)

and

$\displaystyle\beta_{k}^{*}(\boldsymbol{\theta})\beta_{l}(\boldsymbol{\theta})$
$\displaystyle\mapsto(\beta_{k}^{*}(\boldsymbol{\theta})+\Delta
c_{k,0}^{*})(\beta_{l}(\boldsymbol{\theta})+\Delta c_{l,0})$
$\displaystyle=\beta_{k}^{*}(\boldsymbol{\theta})\beta_{l}(\boldsymbol{\theta})+\Delta
c_{k,0}^{*}\Delta c_{l,0}+\Delta
c_{k,0}^{*}\beta_{l}(\boldsymbol{\theta})+\Delta
c_{l,0}\beta_{k}^{*}(\boldsymbol{\theta}).$ (28)

Therefore, the Fubini-Study metric tensor
$\boldsymbol{G}(\boldsymbol{\theta})$ is invariant under the shift in Eq.
(26), i.e.,

$\boldsymbol{G}(\boldsymbol{\theta})\mapsto\boldsymbol{G}(\boldsymbol{\theta}).$
(29)

Notice that the transformations in Eqs. (27) and (28) are similar to those
considered in Ref. Provost and Vallee (1980), where the Fubini-Study metric
tensor $\boldsymbol{G}(\boldsymbol{\theta})$ is constructed so as to be
invariant under multiplication of a global phase factor to the state
$|\Psi_{\cal U}^{(0)}(\boldsymbol{\theta})\rangle$. Indeed, the former
transformations are exactly the same as the latter ones if $\Delta c_{k,0}$
and $\Delta c_{l,0}$ in Eqs. (27) and (28) are pure imaginary. This also
implies that the undetermined and in principle arbitrarily chosen ${\rm
Im}c_{k,0}$ in Eq. (25) can be absorbed into the global phase factor of the
state $|\Psi_{\cal U}^{(0)}(\boldsymbol{\theta})\rangle$, which is analogous
to the case of the first-order perturbation theory for a standard (i.e., not a
generalized) Hermitian eigenvalue problem in quantum mechanics, where the
imaginary part of the coefficient corresponding to $c_{k,0}$ is also
undetermined Schiff (1949). Since
$\partial_{k}\boldsymbol{v}_{0}(\boldsymbol{\theta})$ appears only in the
Fubini-Study metric tensor $\boldsymbol{G}(\boldsymbol{\theta})$ in the
present study, $\Delta c_{k,0}$ and hence $c_{k,0}$ can be chosen arbitrarily.

###  III.2 Derivatives

We now derive analytical expressions of the derivatives, assuming that
$|\psi(\boldsymbol{\theta})\rangle$ in Eq. (2) has a particular form of

$|\psi(\boldsymbol{\theta})\rangle\equiv\prod_{k=N_{\rm
v}}^{1}\hat{U}_{k}(\theta_{k})\hat{W}|0\rangle^{\otimes N},$ (30)

where $\hat{U}_{k}(\theta_{k})$ is a unitary operator parametrized by
$\theta_{k}$ and $\hat{W}$ is a unitary operator independent of the
variational parameters. Since $\\{\hat{U}_{k}(\theta_{k})\\}_{k=1}^{N_{\rm
v}}$ and $\hat{W}$ are unitary, $|\psi(\boldsymbol{\theta})\rangle$ can be
prepared efficiently on a quantum computer. Moreover, we assume that
$\hat{U}_{k}(\theta_{k})$ in Eq. (30) is an exponential of an involutory
operator $\hat{P}_{k}$ of the form

$\hat{U}_{k}(\theta_{k})={\rm e}^{-{\rm
i}\hat{P}_{k}\theta_{k}/2}=\hat{1}\cos{\tfrac{\theta_{k}}{2}}-{\rm
i}\hat{P}_{k}\sin{\tfrac{\theta_{k}}{2}},$ (31)

where $\hat{1}$ is the identity operator and the second equality is because of
$\hat{P}_{k}^{2}=\hat{1}$.

According to the parameter-shift rules Li  _et al._ (2017); Mitarai  _et al._
(2018); Schuld  _et al._ (2019), the derivatives appearing in Eqs. (15) and
(23) can be evaluated as

$\frac{\partial\boldsymbol{H}(\boldsymbol{\theta})}{\partial\theta_{k}}=\frac{1}{2}\left[\boldsymbol{H}(\boldsymbol{\theta}+\tfrac{\pi}{2}\boldsymbol{e}_{k})-\boldsymbol{H}(\boldsymbol{\theta}-\tfrac{\pi}{2}\boldsymbol{e}_{k})\right]$
(32)

and

$\frac{\partial\boldsymbol{S}(\boldsymbol{\theta})}{\partial\theta_{k}}=\frac{1}{2}\left[\boldsymbol{S}(\boldsymbol{\theta}+\tfrac{\pi}{2}\boldsymbol{e}_{k})-\boldsymbol{S}(\boldsymbol{\theta}-\tfrac{\pi}{2}\boldsymbol{e}_{k})\right],$
(33)

where $[\boldsymbol{e}_{k}]_{k^{\prime}}=\delta_{kk^{\prime}}$ is the $N_{\rm
v}$-dimensional basis vector. Note that the derivative of the parametrized
state $|u_{i}(\boldsymbol{\theta})\rangle$ in Eq. (2) itself is also evaluated
simply as

$|\partial_{k}u_{i}(\boldsymbol{\theta})\rangle=\frac{1}{2}|u_{i}(\boldsymbol{\theta}+\pi\boldsymbol{e}_{k})\rangle.$
(34)

Instead of Eqs. (32) and (33), the derivatives of
$\boldsymbol{H}(\boldsymbol{\theta})$ and
$\boldsymbol{S}(\boldsymbol{\theta})$ can also be written as

$\frac{\partial\boldsymbol{H}(\boldsymbol{\theta})}{\partial\theta_{k}}=\boldsymbol{H}_{(,k)}(\boldsymbol{\theta})+[\boldsymbol{H}_{(,k)}(\boldsymbol{\theta})]^{\dagger}$
(35)

and

$\frac{\partial\boldsymbol{S}(\boldsymbol{\theta})}{\partial\theta_{k}}=\boldsymbol{S}_{(,k)}(\boldsymbol{\theta})+[\boldsymbol{S}_{(,k)}(\boldsymbol{\theta})]^{\dagger},$
(36)

where $[\boldsymbol{H}_{(,k)}(\boldsymbol{\theta})]_{ij}\equiv\langle
u_{i}(\boldsymbol{\theta})|\hat{\cal
H}|\partial_{k}u_{j}(\boldsymbol{\theta})\rangle$. By substituting Eq. (34)
into Eqs. (35) and (36), different analytical expressions for the derivatives
can be obtained. As far as numerical simulations are concerned, Eqs. (35) and
(36) are more preferable than Eqs. (32) and (33). This is because the number
of the parametrized states required for Eqs. (35) and (36) is $N_{\rm v}+1$,
i.e., $|\psi(\boldsymbol{\theta})\rangle$ and
$\\{|\psi(\boldsymbol{\theta}+\pi\boldsymbol{e}_{k})\rangle\\}_{k=1}^{N_{\rm
v}}$, while Eqs. (32) and (33) require $2N_{\rm v}$ states
$\\{|\psi(\boldsymbol{\theta}\pm\tfrac{\pi}{2}\boldsymbol{e}_{k})\rangle\\}_{k=1}^{N_{\rm
v}}$, and thus the numerical complexity is approximately half.

##  IV Krylov-extended symmetry-adapted VQE

Equation (9) implies that the selection of the subspace $\cal U$ is crucial
for approximating the ground state (or any target eigenstate) of a Hamiltonian
$\hat{\cal H}$. In the Krylov-subspace SAVQE for the ground state calculation,
the subspace $\cal U$ is constructed as

$\displaystyle{\cal U}$ $\displaystyle={\cal K}_{d_{\cal U}}(\hat{\cal
H},\hat{\cal P}|\psi(\boldsymbol{\theta})\rangle)$ $\displaystyle={\rm
span}\left(\hat{\cal P}|\psi(\boldsymbol{\theta})\rangle,\hat{\cal H}\hat{\cal
P}|\psi(\boldsymbol{\theta})\rangle,\cdots,\hat{\cal H}^{d_{\cal
U}-1}\hat{\cal P}|\psi(\boldsymbol{\theta})\rangle\right),$ (37)

where $\hat{\cal P}$ is a symmetry-projection operator and by definition
$\hat{\cal H}$ and $\hat{\cal P}$ commute to each other, i.e., $[\hat{\cal
H},\hat{\cal P}]=0$. The explicit form of $\hat{\cal P}$ is determined by
considering the Hamiltonian symmetry and the details are described in Sec.
V.4. $|\psi(\boldsymbol{\theta})\rangle$ is a quantum circuit ansatz
parametrized by the variational parameters
$\boldsymbol{\theta}=\\{\theta_{k}\\}_{k=1}^{N_{\rm v}}$ as in Eq. (30) and
generally breaks the symmetry of the Hamiltonian $\cal{H}$. With the
projection operator $\hat{\cal P}$, the state $\hat{\cal
P}|\psi(\boldsymbol{\theta})\rangle$ is projected onto the symmetry sector
relevant to the ground state. It is obvious that the Krylov subspace in Eq.
(37) is a special case of the more general subspace in Eq. (1) with the
operator $\hat{O}_{n}$ in Eq. (2) given now by

$\hat{O}_{n}=\hat{\cal H}^{n}\hat{\cal P}.$ (38)

In the VQE scheme, each term in the matrix elements of
$\boldsymbol{H}(\boldsymbol{\theta})$ and
$\boldsymbol{S}(\boldsymbol{\theta})$ for a given set of the variational
parameters $\boldsymbol{\theta}$ is evaluated on a quantum computer, while a
classical computer is employed to solve the generalized eigenvalue problem in
Eq. (6) and to optimize the variational parameters, here using the NGD method
described in Sec. III, for which each term in the matrix elements of the
Fubini-Study tensor $\boldsymbol{G}(\boldsymbol{\theta})$ is also evaluated on
a quantum computer. As explained in the next section, the nonunitary
projection operator $\hat{\cal P}$ can be expressed by a linear combination of
unitary operators and hence can be implemented on a quantum computer with the
appropriate postprocessing Seki  _et al._ (2020).

To evaluate the expectation value of $\hat{\cal H}^{n}$ on a quantum computer,
several direct Kowalski and Peng (2020); Vallury  _et al._ (2020); Claudino
_et al._ (2021) and approximate Bespalova and Kyriienko (2021); González  _et
al._ (2021) methods have been proposed recently. For a recent review, see for
example Ref. Aulicino  _et al._ . Here, we employ the QPM, a scheme introduced
in Ref. Seki and Yunoki (2021), where the Hamiltonian power $\hat{\cal H}^{n}$
is directly treated by approximating it with a linear combination of unitary
time-evolution operators at $n+1$ different times. Using the second-order
symmetric Suzuki-Trotter (ST) decomposition Trotter (1959); Suzuki (1976), the
Hamiltonian power $\hat{\cal H}^{n}$ is approximated as Seki and Yunoki (2021)

$\hat{\cal H}^{n}=\hat{\cal H}^{n}_{\rm ST}(\Delta)+O(\Delta^{2}),$ (39)

where

$\hat{\cal H}^{n}_{\rm ST}(\Delta)=\frac{{\rm
i}^{n}}{\Delta^{n}}\sum_{k=0}^{n}(-1)^{k}\binom{n}{k}\
\left[\hat{S}_{2}(\Delta/2)\right]^{n-2k},$ (40)

and $\Delta$ is a parameter and should be in general a small real number.
$\hat{S}_{2}(\Delta/2)$ is the second-order symmetric ST decomposition of the
time-evolution operator ${\rm e}^{-{\rm i}\hat{\cal H}\Delta/2}$ and
$O(\Delta^{2})$ in Eq. (39) is the leading systematic error. The explicit form
of $\hat{S}_{2}(\Delta/2)$ depends on the Hamiltonian $\hat{\cal H}$ and it is
given in Sec. VII.2 for the two-component Fermi-Hubbard model in a ladder
lattice structure.

As described in detail in Ref. Seki and Yunoki (2021), the order of the
leading systematic error can be eliminated systematically with the Richardson
extrapolation. For instance, the first-order Richardson extrapolation
$\hat{\cal H}_{\rm ST(1)}^{n}(\Delta)$ of $\hat{\cal H}_{\rm ST}^{n}(\Delta)$
can be given by

$\hat{\cal H}_{\rm ST(1)}^{n}(\Delta)=\frac{1}{3}\left[{4\hat{\cal H}_{\rm
ST}^{n}(\Delta/2)-\hat{\cal H}_{\rm ST}^{n}(\Delta)}\right],$ (41)

which approximates $\hat{\cal H}^{n}$ to $O(\Delta^{4})$, i.e.,

$\hat{\cal H}^{n}=\hat{\cal H}^{n}_{\rm ST(1)}(\Delta)+O(\Delta^{4}).$ (42)

By using Eq. (41), we can approximate $\hat{O}_{n}$ in Eq. (38) as

$\hat{O}_{n}\approx\hat{\cal H}_{\rm ST(1)}^{n}(\Delta)\hat{\cal P}.$ (43)

Accordingly, the Hamiltonian matrix and the overlap matrix are given
respectively as

$[\boldsymbol{H}(\boldsymbol{\theta})]_{ij}\approx\langle\psi(\boldsymbol{\theta})|\hat{\cal
P}\hat{\cal H}_{\rm ST(1)}^{i}(\Delta)\hat{\cal H}\hat{\cal H}_{\rm
ST(1)}^{j}(\Delta)\hat{\cal P}|\psi(\boldsymbol{\theta})\rangle$ (44)

and

$[\boldsymbol{S}(\boldsymbol{\theta})]_{ij}\approx\langle\psi(\boldsymbol{\theta})|\hat{\cal
P}\hat{\cal H}_{\rm ST(1)}^{i}(\Delta)\hat{\cal H}_{\rm
ST(1)}^{j}(\Delta)\hat{\cal P}|\psi(\boldsymbol{\theta})\rangle,$ (45)

where $\hat{\cal P}^{\dagger}=\hat{\cal P}$ and $[\hat{\cal H}_{\rm
ST(1)}^{n}]^{\dagger}=\hat{\cal H}_{\rm ST(1)}^{n}$ are used. On the other
hand, because of $[\hat{\cal H},\hat{\cal P}]=0$, we can also approximate
$\hat{O}_{n}$ in Eq. (38) as

$\hat{O}_{n}\approx\hat{\cal P}\hat{\cal H}_{\rm ST(1)}^{n}(\Delta).$ (46)

Correspondingly, the Hamiltonian matrix and the overlap matrix can be given
respectively as

$[\boldsymbol{H}(\boldsymbol{\theta})]_{ij}\approx\langle\psi(\boldsymbol{\theta})|\hat{\cal
H}_{\rm ST(1)}^{i}(\Delta)\hat{\cal H}\hat{\cal P}\hat{\cal H}_{\rm
ST(1)}^{j}(\Delta)|\psi(\boldsymbol{\theta})\rangle$ (47)

and

$[\boldsymbol{S}(\boldsymbol{\theta})]_{ij}\approx\langle\psi(\boldsymbol{\theta})|\hat{\cal
H}_{\rm ST(1)}^{i}(\Delta)\hat{\cal P}\hat{\cal H}_{\rm
ST(1)}^{j}(\Delta)|\psi(\boldsymbol{\theta})\rangle,$ (48)

where $\hat{\cal P}^{2}=\hat{\cal P}$ is also used. Although the apparent
forms of $[\boldsymbol{H}(\boldsymbol{\theta})]_{ij}$
($[\boldsymbol{S}(\boldsymbol{\theta})]_{ij}$) in Eqs. (44) and (47) [Eqs.
(45) and (48)] are different, they are exactly the same up to $O(\Delta^{4})$
because of $[\hat{\cal P},\hat{\cal H}_{\rm ST(1)}^{n}]=O(\Delta^{4})$.

An important difference between Eqs. (44) and (45) and Eqs. (47) and (48) is
that the latter formulas contain the projection operator $\cal{\hat{P}}$ only
once for each matrix element. This is advantageous for quantum-classical
hybrid computing because the number of measurements required for Eqs. (47) and
(48) is smaller than that for Eqs. (44) and (45) by a factor of the number of
terms in $\hat{\cal P}$, assuming that $\hat{\cal P}$ is given by a linear
combination of unitary operators. For this reason, we adopt Eqs. (47) and (48)
for numerical simulations in Sec. VII. It should be noted that, due to the
presence of the projection operator $\hat{\cal P}$, the number of measurements
required for evaluating the Hamiltonian and overlap matrices in Eqs. (47) and
(48) is increased by a factor of the number of terms in $\hat{\cal P}$ from
that discussed in Ref. Seki and Yunoki (2021).

After optimizing the variational parameters
$\boldsymbol{\theta}=\\{\theta_{k}\\}_{k=1}^{N_{\rm v}}$ in Eq. (12), we
obtain the approximated ground state $|\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta}_{\rm opt})\rangle$ in the Krylov subspace $\cal
U$ defined in Eq. (37), i.e.,

$|\Psi_{0}\rangle\approx|\Psi_{\cal U}^{(0)}(\boldsymbol{\theta}_{\rm
opt})\rangle=\sum_{i=0}^{d_{\cal U}-1}v_{i,0}(\boldsymbol{\theta}_{\rm
opt})\hat{\cal P}\hat{\cal H}_{\rm
ST(1)}^{i}(\Delta)|\psi(\boldsymbol{\theta}_{\rm opt})\rangle.$ (49)

The expectation value of any quantity $\hat{\cal O}$ with respect to
$|\Psi_{\cal U}^{(0)}(\boldsymbol{\theta}_{\rm opt})\rangle$ is then evaluated
simply as

$\displaystyle\langle\Psi_{0}|\hat{\cal O}|\Psi_{0}\rangle$
$\displaystyle\approx$ $\displaystyle\frac{\langle\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta}_{\rm opt})|\hat{\cal O}|\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta}_{\rm opt})\rangle}{\langle\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta}_{\rm opt})|\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta}_{\rm opt})\rangle}$ (50) $\displaystyle=$
$\displaystyle\frac{\boldsymbol{v}_{0}^{\dagger}(\boldsymbol{\theta}_{\rm
opt})\boldsymbol{O}(\boldsymbol{\theta}_{\rm
opt})\boldsymbol{v}_{0}(\boldsymbol{\theta}_{\rm
opt})}{\boldsymbol{v}_{0}^{\dagger}(\boldsymbol{\theta}_{\rm
opt})\boldsymbol{S}(\boldsymbol{\theta}_{\rm
opt})\boldsymbol{v}_{0}(\boldsymbol{\theta}_{\rm opt})},$

where

$[\boldsymbol{O}(\boldsymbol{\theta})]_{ij}=\langle\psi(\boldsymbol{\theta})|\hat{\cal
H}_{\rm ST(1)}^{i}(\Delta)\hat{\cal O}\hat{\cal P}\hat{\cal H}_{\rm
ST(1)}^{j}(\Delta)|\psi(\boldsymbol{\theta})\rangle$ (51)

and $[\hat{\cal O},\hat{\cal P}]=0$ is assumed. This is easily generalized for
a quantity $\hat{\cal O}$ with $[\hat{\cal O},\hat{\cal P}]\neq 0$.

##  V Model, symmetry, and projection operator

###  V.1 Fermi-Hubbard model

As a demonstration, we apply the extended SAVQE method to the two-component
Fermi-Hubbard model. The Hamiltonian $\hat{\cal H}$ of the Fermi-Hubbard model
is given by

$\hat{\cal H}=-t\sum_{\langle
i,j\rangle,\sigma}\left(\hat{c}_{i\sigma}^{\dagger}\hat{c}_{j\sigma}+{\rm
H.c.}\right)+U_{\rm
H}\sum_{i}\left(\hat{n}_{i\uparrow}-\frac{1}{2}\right)\left(\hat{n}_{i\downarrow}-\frac{1}{2}\right),$
(52)

where $\hat{c}_{i\sigma}^{\dagger}$ ($\hat{c}_{i\sigma}$) is the creation
(annihilation) operator of a fermion at site $i\,(=1,2,\dots,L)$ with spin
$\sigma\,(=\uparrow,\downarrow)$ on a lattice having $L$ sites, and
$\hat{n}_{i\sigma}=\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma}$ is the
fermion density operator. $t$ and $U_{\rm H}$ are the hopping and the on-site
interaction parameters, respectively, and $\langle i,j\rangle$ runs over all
nearest-neighbor pairs of sites $i$ and $j$. We assume that $t>0$ and $U_{\rm
H}>0$.

We adopt the Jordan-Wigner transformation Jordan and Wigner (1928); Rodriguez
(1959) of the form

$\hat{c}_{i\sigma}^{\dagger}\overset{\rm
JWT}{=}\frac{1}{2}(\hat{X}_{i_{\sigma}}-{\rm
i}\hat{Y}_{i_{\sigma}})\prod_{k<{i_{\sigma}}}\hat{Z}_{k}$ (53)

and

$\hat{c}_{i\sigma}\overset{\rm JWT}{=}\frac{1}{2}(\hat{X}_{i_{\sigma}}+{\rm
i}\hat{Y}_{i_{\sigma}})\prod_{k<{i_{\sigma}}}\hat{Z}_{k},$ (54)

where “$\overset{\rm JWT}{=}$” implies that the Jordan-Wigner transformation
transforms fermion operators on the left-hand side to a qubit representation
on the right-hand side, and ${\hat{X}}_{i_{\sigma}}$,
${\hat{Y}}_{i_{\sigma}}$, and ${\hat{Z}}_{i_{\sigma}}$ are the Pauli operators
at qubit $i_{\sigma}$. We use the convention that the qubits are ordered as
$i_{\uparrow}=i\,(=1,2,\dots,L)$ and $i_{\downarrow}=i+L\,(=L+1,L+2,\dots,2L)$
with the total number of qubits being $N=2L$ in Eqs. (53) and (54). The
fermion density operator is given by

$\hat{n}_{i\sigma}=\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma}\overset{\rm
JWT}{=}\frac{1}{2}(1-\hat{Z}_{i_{\sigma}}),$ (55)

and thus it is plausible to say that the single-particle state $(i,\sigma)$ is
occupied (unoccupied) if the $i_{\sigma}$th qubit state is
$|1\rangle_{i_{\sigma}}$ ($|0\rangle_{i_{\sigma}}$) because
$\hat{n}_{i\sigma}|1\rangle_{i_{\sigma}}=1|1\rangle_{i_{\sigma}}$
($\hat{n}_{i\sigma}|0\rangle_{i_{\sigma}}=0|0\rangle_{i_{\sigma}}$), where
$|1\rangle_{i_{\sigma}}$ and $|0\rangle_{i_{\sigma}}$ are the eigenstates of
Pauli operator $\hat{Z}_{i_{\sigma}}$ at qubit $i_{\sigma}$, i.e.,
$\hat{Z}_{i_{\sigma}}|0\rangle_{i_{\sigma}}=|0\rangle_{i_{\sigma}}$ and
$\hat{Z}_{i_{\sigma}}|1\rangle_{i_{\sigma}}=-|1\rangle_{i_{\sigma}}$. The
Fermi-Hubbard Hamiltonian in the qubit representation reads

$\hat{\mathcal{H}}\overset{\rm JWT}{=}-\frac{t}{2}\sum_{\sigma}\sum_{\langle
i_{\sigma},j_{\sigma}\rangle}\left(\hat{X}_{i_{\sigma}}\hat{X}_{j_{\sigma}}+\hat{Y}_{i_{\sigma}}\hat{Y}_{j_{\sigma}}\right)\hat{Z}_{{\rm
JW},i_{\sigma}j_{\sigma}}+\frac{U_{\rm
H}}{4}\sum_{i=1}^{L}\hat{Z}_{i}\hat{Z}_{i+L},$ (56)

where $\hat{Z}_{{\rm JW},ij}=\prod_{i\lessgtr k\lessgtr j}\hat{Z}_{k}$ is the
Jordan-Wigner string for $i\lessgtr k\lessgtr j$ and $\hat{Z}_{{\rm
JW},ij}=\hat{1}$ for $i=j\pm 1$.

Figure 1:  Geometry of a $4\times 2$ lattice in which the Fermi-Hubbard model
is defined. Circles represent sites that are labeled by the numbers shown
inside the circles. The hopping term is present only between the nearest-
neighboring sites connected by black solid lines. The lattice sites belonging
to $A$ ($B$) sublattice are denoted with red (blue) circles. The two-fold
rotational axis ($C_{2}$) and the mirror planes ($\sigma_{1},\sigma_{2}$) for
the $C_{2v}$ symmetry are also indicated in the figure.

In this paper, we consider the Fermi-Hubbard model in a $4\times 2$ lattice
under open-boundary conditions (see Fig. 1) at half filling. Therefore, to
represent the variational state, we use $N=16$ qubits, where the first 8
qubits are assigned to the single-particle states labeled by site $i$ and spin
$\uparrow$, and the remaining 8 qubits are assigned to the single-particle
states labeled by site $i$ and spin $\downarrow$ (also see Fig. 8).

###  V.2 Hamiltonian symmetry

We now consider the Hamiltonian symmetry. From the particular geometry of the
cluster shown in Fig. 1, it is obvious that the model has the $C_{2v}$ point-
group symmetry, i.e.,

$\left[\hat{\cal H},\hat{g}_{m}\right]=0,$ (57)

where $\hat{g}_{m}$ is a group element of $C_{2v}$. We note that the $4\times
2$ cluster itself possesses the $D_{2h}$ point-group symmetry that has 8
symmetry operations (i.e., symmetry elements), and the $C_{2v}$ group is a
subgroup of the $D_{2h}$ group. Indeed, the $D_{2h}$ group can be seen as a
direct-product group of the $C_{2v}$ group and the $C_{i}$ group, i.e.,
$D_{2h}=C_{2v}\times C_{i}$. However, in a pure two-dimensional space, the
inversion and the two-fold rotation around the principal axis are equivalent.
Therefore, among the 8 symmetry operations in the $D_{2h}$ group, only the 4
symmetry operations, which are all the symmetry operations of the $C_{2v}$
group, are independent. This is the reason why we consider the $C_{2v}$ point-
group symmetry here.

It can also be shown that the Fermi-Hubbard model considered here is spin
symmetric, i.e.,

$\left[\hat{\cal H},\hat{S}^{2}\right]=0,\quad\left[\hat{\cal
H},\hat{S}_{z}\right]=0,$ (58)

where $\hat{S}^{2}=\hat{S}_{x}^{2}+\hat{S}_{y}^{2}+\hat{S}_{z}^{2}$,
$\hat{S}_{x}=\sum_{i}\hat{S}^{x}_{i}$, $\hat{S}_{y}=\sum_{i}\hat{S}^{y}_{i}$,
and $\hat{S}_{z}=\sum_{i}\hat{S}^{z}_{i}$ with

$\displaystyle\hat{S}_{i}^{x}$
$\displaystyle=\hat{\mathbf{c}}_{i}^{\dagger}\mathbf{s}_{x}\hat{\mathbf{c}}_{i}=\frac{1}{2}\left(\hat{c}_{i\uparrow}^{\dagger}\hat{c}_{i\downarrow}+\hat{c}_{i\downarrow}^{\dagger}\hat{c}_{i\uparrow}\right),$
(59) $\displaystyle\hat{S}_{i}^{y}$
$\displaystyle=\hat{\mathbf{c}}_{i}^{\dagger}\mathbf{s}_{y}\hat{\mathbf{c}}_{i}=\frac{1}{2{\rm
i}}\left(\hat{c}_{i\uparrow}^{\dagger}\hat{c}_{i\downarrow}-\hat{c}_{i\downarrow}^{\dagger}\hat{c}_{i\uparrow}\right),$
(60) $\displaystyle\hat{S}_{i}^{z}$
$\displaystyle=\hat{\mathbf{c}}_{i}^{\dagger}\mathbf{s}_{z}\hat{\mathbf{c}}_{i}=\frac{1}{2}\left(\hat{n}_{i\uparrow}-\hat{n}_{i\downarrow}\right).$
(61)

Here, we have introduced that
$\hat{\mathbf{c}}_{i}^{\dagger}=\begin{bmatrix}\hat{c}_{i\uparrow}^{\dagger}&\hat{c}_{i\downarrow}^{\dagger}\end{bmatrix}$
and

$\mathbf{s}_{x}=\frac{1}{2}\begin{bmatrix}0&1\\\
1&0\end{bmatrix},\quad\mathbf{s}_{y}=\frac{1}{2}\begin{bmatrix}0&-{\rm i}\\\
{\rm i}&0\end{bmatrix},\quad\mathbf{s}_{z}=\frac{1}{2}\begin{bmatrix}1&0\\\
0&-1\end{bmatrix}.$ (62)

Hereafter, we refer to Eq. (58) as the $S$ symmetry.

Due to the bipartite structure of the hopping term (see Fig. 1), the Fermi-
Hubbard model has another symmetry Yang (1989)

$\left[\hat{\cal H},\hat{\eta}^{2}\right]=0,\quad\left[\hat{\cal
H},\hat{\eta}_{z}\right]=0,$ (63)

where
$\hat{\eta}^{2}=\hat{\eta}_{x}^{2}+\hat{\eta}_{y}^{2}+\hat{\eta}_{z}^{2}$,
$\hat{\eta}_{x}=\sum_{i}\hat{\eta}^{x}_{i}$,
$\hat{\eta}_{y}=\sum_{i}\hat{\eta}^{y}_{i}$, and
$\hat{\eta}_{z}=\sum_{i}\hat{\eta}^{z}_{i}$ with

$\displaystyle\hat{\eta}_{i}^{x}$
$\displaystyle=\hat{\mathbf{b}}_{i}^{\dagger}\mathbf{s}_{x}\hat{\mathbf{b}}_{i}=\frac{{\rm
e}^{{\rm
i}\phi_{i}}}{2}\left(\hat{c}_{i\uparrow}^{\dagger}\hat{c}_{i\downarrow}^{\dagger}+\hat{c}_{i\downarrow}\hat{c}_{i\uparrow}\right),$
(64) $\displaystyle\hat{\eta}_{i}^{y}$
$\displaystyle=\hat{\mathbf{b}}_{i}^{\dagger}\mathbf{s}_{y}\hat{\mathbf{b}}_{i}=\frac{{\rm
e}^{{\rm i}\phi_{i}}}{2{\rm
i}}\left(\hat{c}_{i\uparrow}^{\dagger}\hat{c}_{i\downarrow}^{\dagger}-\hat{c}_{i\downarrow}\hat{c}_{i\uparrow}\right),$
(65) $\displaystyle\hat{\eta}_{i}^{z}$
$\displaystyle=\hat{\mathbf{b}}_{i}^{\dagger}\mathbf{s}_{z}\hat{\mathbf{b}}_{i}=\frac{1}{2}\left(\hat{n}_{i\uparrow}+\hat{n}_{i\downarrow}-1\right).$
(66)

Here, we have introduced that
$\hat{\mathbf{b}}_{i}^{\dagger}=\begin{bmatrix}\hat{c}_{i\uparrow}^{\dagger}&{\rm
e}^{{\rm i}\phi_{i}}\hat{c}_{i\downarrow}\end{bmatrix}$ and

${\rm e}^{{\rm i}\phi_{i}}={\rm e}^{-{\rm
i}\phi_{i}}=\begin{cases}+1&\text{for $i\in A$ sublattice}\\\ -1&\text{for
$i\in B$ sublattice}\end{cases}$ (67)

with $A$ and $B$ sublattices being indicated in Fig. 1. Hereafter, we refer to
Eq. (63) as the $\eta$ symmetry. Note that, as in the case of the spin
operators ${\hat{S}}_{i}^{\alpha}$ with $\alpha=x,y,z$, the $\eta$ operators
${\hat{\eta}}_{i}^{\alpha}$ in Eqs. (64)-(66) satisfy the commutation
relations of the angular momentum. It is also important to note that
$[{\hat{S}}^{\alpha}_{i},{\hat{\eta}}^{\beta}_{j}]=0$ with
$\alpha,\beta=x,y,z$. Therefore, any eigenstate of $\hat{\cal H}$ can be
simultaneously an eigenstate of ${\hat{S}}^{2}$, ${\hat{S}}_{z}$,
${\hat{\eta}}^{2}$, and ${\hat{\eta}}_{z}$, and thus it is characterized with
these eigenvalues.

###  V.3 Spin- and $\eta$-singlet ground state

Via the particle-hole transformation for spin-down fermions

$\displaystyle\hat{c}_{i\downarrow}$ $\displaystyle\mapsto{\rm e}^{-{\rm
i}\phi_{i}}\hat{c}_{i\downarrow}^{\dagger}$ (68)

with spin-up fermions unaltered, the spin and the $\eta$ operators are
transformed into each other because
$\hat{\mathbf{c}}_{i}^{\dagger}\mapsto\hat{\mathbf{b}}_{i}^{\dagger}$ and
$\hat{\mathbf{b}}_{i}^{\dagger}\mapsto\hat{\mathbf{c}}_{i}^{\dagger}$,
assuming that ${\rm e}^{{\rm i}\phi_{i}}$ is given in Eq. (67). As stated
above, this observation also implies that the $\eta$ operators satisfy the
same commutation relations as the spin operators. Note that a unitary operator
$\hat{U}_{i}^{\rm ph}=\hat{c}_{i\downarrow}+{\rm e}^{-{\rm
i}\phi_{i}}\hat{c}_{i\downarrow}^{\dagger}$ Tasaki (2020) satisfies the
relation $(\hat{U}_{i}^{\rm
ph})^{\dagger}\hat{c}_{i\downarrow}\hat{U}_{i}^{\rm ph}={\rm e}^{-{\rm
i}\phi_{i}}\hat{c}_{i\downarrow}^{\dagger}$, which corresponds to the mapping
in Eq. (68).

With the same particle-hole transformation, the repulsive Fermi-Hubbard model
on a bipartite lattice maps to the attractive Fermi-Hubbard model on the same
lattice Shiba (1972). The particle-hole transformation is useful to gain
insight into the correspondence between the repulsive and attractive Fermi-
Hubbard models as well as that between the spin and charge degrees of freedom
(see for example Refs. You  _et al._ (2017); Karakuzu  _et al._ (2018);
Fujiuchi  _et al._ (2020); Li (2020)). Indeed, the fact Lieb (1989) that the
ground state of the attractive Fermi-Hubbard model for even number of fermions
in the absence of magnetic field has $S=0$ (i.e., spin singlet) implies that
the ground state of the repulsive Fermi-Hubbard model on a bipartite lattice
at half filling has $\eta=0$ (i.e., $\eta$ singlet). Furthermore, the ground
state of the repulsive Fermi-Hubbard model on a bipartite lattice at half
filling is proved to have $S=\frac{1}{2}||B|-|A||$, where $|B|$ ($|A|$) is the
number of sites in the $B$ ($A$) sublattice Lieb (1989). Since $|B|=|A|$ in
the present case, the ground state has $S=0$.

###  V.4 Projection operators

To take into account the spatial, $S$, and $\eta$ symmetry, we adopt the
projection operator $\hat{\cal P}$ of the form

$\hat{\cal P}=\hat{\cal P}^{(\eta)}\hat{\cal P}^{(S)}\hat{\cal P}^{(\alpha)},$
(69)

where $\hat{\cal P}^{(\alpha)}$, $\hat{\cal P}^{(S)}$, and $\hat{\cal
P}^{(\eta)}$ are the projection operators of the spatial, $S$, and $\eta$
symmetry, respectively. In this section, we specify the form of these
projection operators. Note that these three projection operators commute with
each other.

####  V.4.1 Spatial symmetry

The projection operator of the spatial symmetry for an irreducible
representation $\alpha$ of group $\cal G$ is given by Inui  _et al._ (1990)

$\hat{\cal P}^{(\alpha)}_{\mu\nu}=\frac{d_{\alpha}}{|\cal
G|}\sum_{\hat{g}_{m}\in{\cal
G}}\left[D^{(\alpha)}_{\mu\nu}(\hat{g}_{m})\right]^{*}\hat{g}_{m},$ (70)

where $|{\cal G}|$ is the number of group elements in ${\cal G}$ (i.e., order
of group $\cal G$), $d_{\alpha}$ is the dimension of the irreducible
representation $\alpha$, $D_{\mu\nu}^{(\alpha)}(\hat{g}_{m})$ is the
$(\mu,\nu)$ entry of the representation matrix in the irreducible
representation $\alpha$ for the symmetry operation $\hat{g}_{m}$, and the sum
runs over all the group elements. The character of the representation matrix
$D^{(\alpha)}_{\mu\nu}(\hat{g}_{m})$ is defined by
$\chi^{(\alpha)}(\hat{g}_{m})=\sum_{\mu=1}^{d_{\alpha}}D^{(\alpha)}_{\mu\mu}(\hat{g}_{m})$.
Note that $D_{\mu\nu}^{(\alpha)}$ satisfies the orthogonality relation

$\sum_{\hat{g}_{m}\in{\cal
G}}\left[D_{\mu\nu}^{(\alpha)}(\hat{g}_{m})\right]^{*}D_{\mu^{\prime}\nu^{\prime}}^{(\alpha^{\prime})}(\hat{g}_{m})=\frac{|{\cal
G}|}{d_{\alpha}}\delta_{\alpha\alpha^{\prime}}\delta_{\mu\mu^{\prime}}\delta_{\nu\nu^{\prime}}.$
(71)

$\hat{\cal P}_{\mu\nu}^{(\alpha)}$ extracts a basis state $|\alpha\mu\rangle$
from $|\alpha\nu\rangle$, i.e., $\hat{\cal
P}_{\mu\nu}^{(\alpha)}|\alpha\nu\rangle=|\alpha\mu\rangle$, where
$|\alpha\nu\rangle$ is a state such that
$\hat{g}_{m}|\alpha\nu\rangle=\sum_{\mu}D_{\mu\nu}^{(\alpha)}(\hat{g}_{m})|\alpha\mu\rangle$.
Strictly speaking, $\hat{\cal P}_{\mu\nu}^{(\alpha)}$ is a projection operator
only when $\mu=\nu$ because $\hat{\cal P}_{\mu\nu}^{(\alpha)}\hat{\cal
P}_{\nu^{\prime}\mu^{\prime}}^{(\alpha^{\prime})}=\delta_{\alpha\alpha^{\prime}}\delta_{\nu\nu^{\prime}}\hat{\cal
P}_{\mu\mu^{\prime}}^{(\alpha)}$, but here we loosely use the term “projection
operator” for $\hat{\cal P}_{\mu\nu}^{(\alpha)}$.

For the point group $C_{2v}$, $|{\cal G}|=4$, $d_{\alpha}=1$, and
$\\{\hat{g}_{m}\\}_{m=1}^{|{\cal
G}|}=\\{\hat{I},\hat{C}_{2},\hat{\sigma}_{1},\hat{\sigma_{2}}\\}$, where
$\hat{C}_{2}$ is the $\pi$ rotation about the center of the $4\times 2$
cluster, and $\sigma_{1}$ and $\sigma_{2}$ are reflections with respect to the
corresponding planes (see Fig. 1). Omitting the subscript $\mu\nu$ in Eq. (70)
because of $d_{\alpha}=1$, we have

$\hat{\cal
P}^{(\alpha)}=\frac{1}{4}\sum_{m=1}^{4}\left[\chi^{(\alpha)}(\hat{g}_{m})\right]^{*}\hat{g}_{m}.$
(72)

The ground state of the Fermi-Hubbard model studied here belongs to the
irreducible representation $\alpha=A_{1}$, in which
$\chi^{(\alpha)}(\hat{g}_{m})=1$ for all $\hat{g}_{m}$.

####  V.4.2 $S$ and $\eta$ symmetry

Since both $S$ and $\eta$ operators satisfy the same commutation relations of
the angular momentum, we denote by $J$ either $S$ or $\eta$ for convenience.
According to the theory of rotation group, the projection operator of the $J$
symmetry ($J=S,\eta$) is given by Ring and Schuck (1980); Schmid (2004);
Mizusaki and Imada (2004); Bally and Bender (2021)

$\hat{\cal P}^{(J)}_{MK}=\frac{2J+1}{\Omega}\int{\rm
d}\omega\left[D^{(J)}_{MK}(\omega)\right]^{*}\hat{R}(\omega),$ (73)

where

$\hat{R}(\omega)={\rm e}^{-{\rm i}\alpha\hat{J}_{z}}{\rm e}^{-{\rm
i}\beta\hat{J}_{y}}{\rm e}^{-{\rm i}\gamma\hat{J}_{z}}$ (74)

is the rotation operator with $\omega=(\alpha,\beta,\gamma)$, i.e., the
parameters of rotation group (Euler angles) specifying the group element, and

$\Omega=\int{\rm d}\omega=\int_{0}^{2\pi}{\rm d}\alpha\int_{0}^{\pi}{\rm
d}\beta\sin{\beta}\int_{0}^{2\pi}{\rm d}\gamma=8\pi^{2}$ (75)

defines the volume of the parameter $\omega$ space.

$D_{MK}^{(J)}(\omega)$ in Eq. (73) is the Wigner’s $D$ function defined by

$D^{(J)}_{MK}(\omega)=\langle JM|\hat{R}(\omega)|JK\rangle={\rm e}^{-{\rm
i}\alpha M}d^{(J)}_{MK}(\beta){\rm e}^{-{\rm i}\gamma K}$ (76)

with $|JM\rangle$ being an eigenstate of $\hat{J}^{2}$ and $\hat{J}_{z}$ such
that $\hat{J}^{2}|JM\rangle=J(J+1)|JM\rangle$ and
$\hat{J}_{z}|JM\rangle=M|JM\rangle$, and

$d^{(J)}_{MK}(\beta)=\langle JM|{\rm e}^{-{\rm i}\beta\hat{J}_{y}}|JK\rangle$
(77)

is the Wigner’s small $d$ function that is real Rose (1957); Wigner (1959);
Varshalovich  _et al._ (1988). Note that the Wigner’s $D$ function
$D_{MK}^{(J)}(\omega)$ satisfies the orthogonality relation

$\int{\rm
d}\omega\left[D_{MK}^{(J)}(\omega)\right]^{*}D_{M^{\prime}K^{\prime}}^{(J^{\prime})}(\omega)=\frac{\Omega}{2J+1}\delta_{JJ^{\prime}}\delta_{MM^{\prime}}\delta_{KK^{\prime}}.$
(78)

Therefore, $\hat{\cal P}_{MK}^{(J)}$ extracts a basis state $|JM\rangle$ from
$|JK\rangle$, i.e., $\hat{\cal P}_{MK}^{(J)}|JK\rangle=|JM\rangle$ because
$\hat{R}(\omega)|JK\rangle=\sum_{M}D_{MK}^{(J)}(\omega)|JM\rangle$. In Eqs.
(75) and (78), $J$ is assumed to be integer. For a half integer $J$, the range
of either $\alpha$ or $\gamma$ integration should be $[0,4\pi]$ and hence
$\Omega=16\pi^{2}$ Inui  _et al._ (1990); Bally and Bender (2021);
Varshalovich  _et al._ (1988). Details on $D_{MK}^{(J)}(\omega)$ and
$d_{MK}^{(J)}(\beta)$ can be found, for example, in Ref. Varshalovich  _et
al._ (1988).

Since our target state is the ground state of the Fermi-Hubbard model in the
sector of $\eta_{z}=0$ (i.e., half filling) and $S_{z}=0$ (i.e., zero
magnetization), we can set $M=0$ in Eq. (73) for both $J=\eta$ and $S$. In
addition, the variational state $|\psi(\boldsymbol{\theta})\rangle$, whose
concrete form is described in Sec. VII.1, is constructed to satisfy
$\hat{J}_{z}|\psi(\boldsymbol{\theta})\rangle=0$, implying that we can set
$K=0$ in Eq. (73) for $J=\eta$ and $S$. For $K=M=0$, the Wigner’s $D$ function
is given by

$D^{J}_{00}(\omega)=d^{J}_{00}(\beta)=P_{J}(\cos{\beta}),$ (79)

where $P_{J}$ is the $J$th order Legendre polynomial. By integrating out
$\gamma$, the state $\hat{\cal
P}^{(J)}|\psi(\boldsymbol{\theta})\rangle\equiv\hat{\cal
P}^{(J)}_{00}|\psi(\boldsymbol{\theta})\rangle$ is given by Tahara and Imada
(2008)

$\displaystyle\hat{\cal P}^{(J)}|\psi(\boldsymbol{\theta})\rangle$
$\displaystyle=\frac{2J+1}{4\pi}\int_{0}^{2\pi}{\rm d}\alpha\int_{0}^{\pi}{\rm
d}\beta\sin{\beta}P_{J}(\cos{\beta})$ $\displaystyle\times{\rm e}^{-{\rm
i}\alpha\hat{J}_{z}}{\rm e}^{-{\rm
i}\beta\hat{J}_{y}}|\psi(\boldsymbol{\theta})\rangle.$ (80)

Furthermore, when a matrix element $\langle\psi^{\prime}|\hat{\cal
P}^{(J)}|\psi(\boldsymbol{\theta})\rangle$ with
$\hat{J}_{z}|\psi^{\prime}\rangle=0$ is considered, one can also integrate out
$\alpha$ as Tahara and Imada (2008)

$\displaystyle\langle\psi^{\prime}|\hat{\cal
P}^{(J)}|\psi(\boldsymbol{\theta})\rangle$
$\displaystyle=\frac{2J+1}{2}\int_{0}^{\pi}{\rm
d}\beta\sin{\beta}P_{J}(\cos{\beta})$
$\displaystyle\times\langle\psi^{\prime}|{\rm e}^{-{\rm
i}\beta\hat{J}_{y}}|\psi(\boldsymbol{\theta})\rangle.$ (81)

Examples of $|\psi^{\prime}\rangle$ particularly relevant for the present
study include $|\psi(\boldsymbol{\theta})\rangle$, $\hat{\cal
H}|\psi(\boldsymbol{\theta})\rangle$, and
$\hat{J}^{2}|\psi(\boldsymbol{\theta})\rangle$. Moreover, Eqs. (80) and (81)
still hold even when $|\psi(\boldsymbol{\theta})\rangle$ in these equations is
replaced with $\hat{\cal H}_{\rm ST(1)}^{n}|\psi(\boldsymbol{\theta})\rangle$
because $\hat{J}_{z}\hat{\cal H}_{\rm
ST(1)}^{n}|\psi(\boldsymbol{\theta})\rangle=0$ for the ST decomposition
employed (see Sec. VII.2).

###  V.5 Projection operator as a linear combination of unitaries

Since the rotation group is a continuous group, $\hat{\cal P}^{(J)}$ involves
integration over continuous variables, i.e., Euler angles, as in Eqs. (80) and
(81). However, for numerical simulations as well as quantum-classical hybrid
calculations, a proper discretization of the integration is necessary. The
integration over $\alpha$ and $\beta$ in Eq. (80) can be discretized with the
trapezoidal rule and the Gauss-Legendre quadrature, respectively. By omitting
irrelevant normalization factor, $\hat{\cal
P}^{(J)}|\psi(\boldsymbol{\theta})\rangle$ is now approximated as

$\displaystyle\hat{\cal P}^{(J)}|\psi(\boldsymbol{\theta})\rangle$
$\displaystyle\approx$ $\displaystyle\frac{2J+1}{2N_{J,\rm
azimuth}}\sum_{i=1}^{N_{J,\rm azimuth}}\sum_{j=1}^{N_{J,\rm
polar}}w_{j,J}P_{J}(\cos{\beta}_{j,J}){\rm e}^{-{\rm
i}\alpha_{i,J}\hat{J}_{z}}{\rm e}^{-{\rm
i}\beta_{j,J}\hat{J}_{y}}|\psi(\boldsymbol{\theta})\rangle,$ (82)

where $\\{\alpha_{i,J}\\}_{i=1}^{N_{J,\rm azimuth}}$ are $N_{J,\rm azimuth}$
azimuth angles equally spaced in $[0,2\pi]$,
$\\{\beta_{j,J}\\}_{j=1}^{N_{J,\rm polar}}$ are polar angles such that
$P_{N_{J,\rm polar}}(\cos{\beta_{j,J}})=0$, and $\\{w_{j,J}\\}_{j=1}^{N_{J,\rm
polar}}$ are the corresponding integration weight of the Gauss-Legendre
quadrature.

The full symmetry-projected state is now given approximately as

$\displaystyle\hat{\cal P}^{(\eta)}\hat{\cal P}^{(S)}\hat{\cal
P}^{(\alpha)}|\psi(\boldsymbol{\theta})\rangle$ $\displaystyle\approx$
$\displaystyle\frac{2\eta+1}{2N_{\eta,\rm azimuth}}\frac{2S+1}{2N_{S,\rm
azimuth}}\frac{d_{\alpha}}{|{\cal G}|}\sum_{i=1}^{N_{\eta,\rm
azimuth}}\sum_{j=1}^{N_{\eta,\rm polar}}\sum_{k=1}^{N_{S,\rm
azimuth}}\sum_{l=1}^{N_{S,\rm polar}}\sum_{m=1}^{|\cal G|}$
$\displaystyle\times$ $\displaystyle
w_{j,\eta}w_{l,S}P_{\eta}(\cos{\beta}_{j,\eta})P_{S}(\cos{\beta}_{l,S})\left[\chi^{(\alpha)}(\hat{g}_{m})\right]^{*}$
$\displaystyle\times$ $\displaystyle{\rm e}^{-{\rm
i}\alpha_{i,\eta}\hat{\eta}_{z}}{\rm e}^{-{\rm
i}\beta_{j,\eta}\hat{\eta}_{y}}{\rm e}^{-{\rm i}\alpha_{k,S}\hat{S}_{z}}{\rm
e}^{-{\rm
i}\beta_{l,S}\hat{S}_{y}}\hat{g}_{m}|\psi(\boldsymbol{\theta})\rangle.$ (83)

Note that the operators in the last line are unitary and thus the symmetry-
projected state is evaluated by applying a linear combination of unitary
operators to the state. The symmetry-projected state given in Eq. (83) is in
general not normalized even though the state
$|\psi(\boldsymbol{\theta})\rangle$ is normalized. However, this is not a
problem because the normalization of states spanning the Krylov subspace $\cal
U$ is not required (see Sec II and Sec. IV) but the normalization of the
approximated ground state is guaranteed by Eq. (10).

A matrix element similar to that in Eq. (81) but for the full projection
operators can be evaluated simply by setting $N_{J,\rm azimuth}=1$ and ${\rm
e}^{-{\rm i}\alpha_{1,J}\hat{J}_{z}}=\hat{I}$ in Eq. (83) and taking the
overlap with $\langle\psi^{\prime}|$, i.e.,

$\displaystyle\langle\psi^{\prime}|\hat{\cal P}^{(\eta)}\hat{\cal
P}^{(S)}\hat{\cal P}^{(\alpha)}|\psi(\boldsymbol{\theta})\rangle$
$\displaystyle\approx$
$\displaystyle\frac{2\eta+1}{2}\frac{2S+1}{2}\frac{d_{\alpha}}{|{\cal
G}|}\sum_{j=1}^{N_{\eta,\rm polar}}\sum_{l=1}^{N_{S,\rm
polar}}\sum_{m=1}^{|\cal G|}$ $\displaystyle\times$ $\displaystyle
w_{j,\eta}w_{l,S}P_{\eta}(\cos{\beta}_{j,\eta})P_{S}(\cos{\beta}_{l,S})\left[\chi^{(\alpha)}(\hat{g}_{m})\right]^{*}$
$\displaystyle\times$ $\displaystyle\langle\psi^{\prime}|{\rm e}^{-{\rm
i}\beta_{j,\eta}\hat{\eta}_{y}}{\rm e}^{-{\rm
i}\beta_{l,S}\hat{S}_{y}}\hat{g}_{m}|\psi(\boldsymbol{\theta})\rangle.$ (84)

Concerning quantum-classical hybrid calculations, Eq. (84) implies that the
matrix element on the left-hand side can be estimated by evaluating the matrix
elements ${\Big{\\{}}{\Big{\\{}}{\Big{\\{}}\langle\psi^{\prime}|{\rm e}^{-{\rm
i}\beta_{j,\eta}\hat{\eta}_{y}}{\rm e}^{-{\rm
i}\beta_{l,S}\hat{S}_{y}}\hat{g}_{m}|\psi(\boldsymbol{\theta})\rangle{\Big{\\}}}_{m=1}^{|{\cal
G}|}{\Big{\\}}}_{l=1}^{N_{S,\rm polar}}{\Big{\\}}}_{j=1}^{N_{\eta,\rm polar}}$
with $N_{\eta,\rm polar}N_{S,\rm polar}|\cal G|$ different circuit structures
using quantum computers, and then adding all of them with the proper weights
on classical computers. We should recall that $d_{\alpha}=1$ and $|{\cal
G}|=4$ for the $C_{2v}$ point group, $P_{\eta}(\cos{\beta}_{j,\eta})=1$ for
$\eta=0$, $P_{S}(\cos{\beta}_{l,S})=1$ for $S=0$, and
$\chi^{(\alpha)}(\hat{g}_{m})=1$ for $\alpha=A_{1}$.

While Eq. (84) suffices for our purpose, we briefly discuss how many of the
integration points $N_{J,\rm polar}$ and $N_{J,\rm azimuth}$ in Eq. (83) are
required for the variational state $|\psi(\boldsymbol{\theta})\rangle$ that
will be described in Sec. VII.1. We find numerically that typically $N_{S,\rm
polar}=2$ and $N_{S,\rm azimuth}=4$ ($N_{\eta,\rm polar}=3$ and $N_{\eta,\rm
azimuth}=5$) are enough to ensure that the calculated value of $S$ ($\eta$)
with respect to the state $\hat{\cal P}|\psi(\boldsymbol{\theta})\rangle$ can
be regarded as integer within the double-precision arithmetic. We also find
that typically the lowest-order Lebedev quadrature $\hat{\cal
P}^{(J)}\approx\sum_{i=1}^{6}{\rm e}^{-{\rm i}\alpha_{i,J}\hat{J}_{z}}{\rm
e}^{-{\rm i}\beta_{i,J}\hat{J}_{y}}$, which uses six integration points
$(\alpha_{i,J},\beta_{i,J})=(0,0),(0,\pi/2),(\pi/2,\pi/2),(\pi,\pi/2),(3\pi/2,\pi/2),$
and $(0,\pi)$ with equal integration weights, is adequate at least for the
$S$-symmetry projection. On the other hand, when we use Eq. (84), where
$N_{J,\rm azimuth}=1$ by construction, we find that typically $N_{S,\rm
polar}=2$ ($N_{\eta,\rm polar}=3$) is sufficient to evaluate $S$ ($\eta$)
within the double-precision arithmetic. However, to ensure that the symmetry
projection is made essentially exactly, we set $N_{J,\rm polar}=4$ for both
$J=S$ and $\eta$ projections in our numerical simulations.

##  VI Symmetry operation on quantum circuit

In this section, we describe how the spatial symmetry operators $\hat{g}_{m}$
and the rotation operators ${\rm e}^{-{\rm i}\alpha\hat{J}_{y}}$ and ${\rm
e}^{-{\rm i}\beta\hat{J}_{z}}$ with $J=S,\eta$ for fermions can be implemented
on a quantum circuit, assuming the qubit representation of fermion operators
by the Jordan-Wigner transformation in Eqs. (53) and (54).

###  VI.1 Quantum circuit for spatial symmetry operations

Let $\hat{g}_{m}$ be a spatial-symmetry operation that transfers the local
state at the $i$th qubit to the $m(i)$th qubit. Such operation can be
represented as a permutation

${\sigma}_{m}\equiv\underbrace{\left(\begin{matrix}N&N-1&\cdots&2&1\\\
m(N)&m(N-1)&\cdots&m(2)&m(1)\\\ \end{matrix}\right)}_{N\text{\ columns}},$
(85)

where $N$ is the number of qubits. We first briefly review an implementation
of the spatial-symmetry operations for spins, and then proceed to the spatial-
symmetry operations for fermions.

####  VI.1.1 Spatial symmetry operations for spins

For a system composed of spins, each having spin 1/2, the spatial-symmetry
operator $\hat{g}_{m}$ can be expressed as a product of nearest-neighbor swap
operators $\hat{\cal S}_{i\delta(i)}$

$\hat{g}_{m}=\prod_{{\sigma}_{m}}\hat{\cal S}_{i\delta(i)},$ (86)

because any permutation can be expressed as a product of transpositions (i.e.,
nearest-neighbor swap operations). Here, $\delta(i)(=i\pm 1)$ denotes a
neighboring qubit of qubit $i$, and the product $\prod_{{\sigma}_{m}}$ should
contain sequences of the swap operators of the form $\hat{\cal
S}_{\delta(\cdots\delta(\delta(k))),m(k)}\cdots\hat{\cal
S}_{\delta(k),\delta(\delta(k))}\hat{\cal S}_{k,\delta(k)}$ for every $k$.
Such a product in Eq. (86) can be constructed accordingly to the “Amida
lottery” construction Seki  _et al._ (2020), as shown in Fig. 2(a) for the
case of

${\sigma}_{m}=\left(\begin{matrix}6&5&4&3&2&1\\\ 4&3&2&1&6&5\\\
\end{matrix}\right)$ (87)

with
$|b_{6}\,b_{5}\,b_{4}\,b_{3}\,b_{2}\,b_{1}\rangle=|b_{6}\rangle_{6}|b_{5}\rangle_{5}|b_{4}\rangle_{4}|b_{3}\rangle_{3}|b_{2}\rangle_{2}|b_{1}\rangle_{1}$
being transferred to
$|b_{2}\,b_{1}\,b_{6}\,b_{5}\,b_{4}\,b_{3}\rangle=|b_{2}\rangle_{6}|b_{1}\rangle_{5}|b_{6}\rangle_{4}|b_{5}\rangle_{3}|b_{4}\rangle_{2}|b_{3}\rangle_{1}$,
where $|b\rangle_{i}$ is the local state $b$ ($=0$ or $1$) at the $i$th qubit.

Figure 2:  Amida-lottery construction of a symmetry operator $\hat{g}_{m}$ on
a quantum circuit for (a) spin models in the computational basis Seki  _et
al._ (2020) and (b) fermion models in an occupation basis [see Eq. (94)]. The
figure refers to the case of $N=6$ qubits. The vertical lines represent qubits
and ${\rm sgn}(\cal R)$ in (b) denotes an appropriate sign factor for fermions
[see Eq. (95) and also Eq. (163)]. In the middle panel of (a), swap gates
$\hat{\cal S}_{i\delta(i)}$ are highlighted with thick red lines. In (b), a
fermionic-swap gate $\hat{\cal F}_{i\delta(i)}$ is represented by a two-qubit
gate with symbol $\cal F$. Assuming that the circuit evolves from top to
bottom, the middle panels in (a) and (b) show that the symmetry operation
$\hat{g}_{m}$ in this example can be expressed as $\hat{g}_{m}=\hat{\cal
S}_{45}\hat{\cal S}_{56}\hat{\cal S}_{34}\hat{\cal S}_{45}\hat{\cal
S}_{23}\hat{\cal S}_{34}\hat{\cal S}_{12}\hat{\cal S}_{23}$ for spin models
and $\hat{g}_{m}=\hat{\cal F}_{45}\hat{\cal F}_{56}\hat{\cal F}_{34}\hat{\cal
F}_{45}\hat{\cal F}_{23}\hat{\cal F}_{34}\hat{\cal F}_{12}\hat{\cal F}_{23}$
for fermion models, respectively.

We note that swap operator $\hat{\cal S}_{ij}$ acts on qubits $i$ and $j$ as

$\displaystyle\hat{\cal
S}_{ij}|0\rangle_{i}|0\rangle_{j}=|0\rangle_{i}|0\rangle_{j},$ (88)
$\displaystyle\hat{\cal
S}_{ij}|0\rangle_{i}|1\rangle_{j}=|1\rangle_{i}|0\rangle_{j},$ (89)
$\displaystyle\hat{\cal
S}_{ij}|1\rangle_{i}|0\rangle_{j}=|0\rangle_{i}|1\rangle_{j},$ (90)
$\displaystyle\hat{\cal
S}_{ij}|1\rangle_{i}|1\rangle_{j}=|1\rangle_{i}|1\rangle_{j}.$ (91)

The matrix representation of $\hat{\cal S}_{ij}$ in these four basis states
$|0\rangle_{i}|0\rangle_{j}$, $|0\rangle_{i}|1\rangle_{j}$,
$|1\rangle_{i}|0\rangle_{j}$, and $|1\rangle_{i}|1\rangle_{j}$ is thus given
by

$\hat{\cal S}_{ij}\overset{\cdot}{=}\begin{bmatrix}1&0&0&0\\\ 0&0&1&0\\\
0&1&0&0\\\ 0&0&0&1\end{bmatrix},$ (92)

where $\overset{\cdot}{=}$ implies that the operator on the left-hand side can
be represented in the matrix form on the right-hand side with the basis states
given above.

####  VI.1.2 Spatial symmetry operations for fermions

For a system composed of fermions, $\hat{g}_{m}$ has to take into account the
anticommutation relations of fermion operators, in addition to the permutation
of local states. This implies that symmetry operators can no longer be given
as a product of the swap operators for fermions in general. As it is explained
in Appendix B, the symmetry operation for fermions can be implemented simply
by replacing the swap operators $\hat{\cal S}_{i\delta(i)}$ in Eq. (86) with
the fermionic-swap operators $\hat{\cal F}_{i\delta(i)}$:

$\hat{g}_{m}=\prod_{{\sigma}_{m}}\hat{\cal F}_{i\delta(i)}.$ (93)

Here, $\hat{\cal F}_{i\delta(i)}$ transforms a fermion creation operator
$\hat{c}_{i}^{\dagger}$ in such a way that $\hat{\cal
F}_{i\delta(i)}\hat{c}_{i}^{\dagger}\hat{\cal
F}_{i\delta(i)}^{-1}=\hat{c}_{\delta(i)}^{\dagger}$, implying that
$\hat{g}_{m}\hat{c}_{i}^{\dagger}\hat{g}_{m}^{-1}=\hat{c}_{m(i)}^{\dagger}$.
Note that the subscript $i\ (=1,2,\cdots,N)$ here labels all the single-
particle local states including sites and spins.

Before discussing quantum circuits for fermionic symmetry operations, let us
first explain how a state of the form

$\displaystyle|b_{N}b_{N-1}\cdots b_{1}\rangle_{\rm F}$
$\displaystyle\equiv\prod_{i=N}^{1}(\hat{c}_{i}^{\dagger})^{b_{i}}|0\rangle_{\rm
F}$
$\displaystyle=(\hat{c}_{N}^{\dagger})^{b_{N}}(\hat{c}_{N-1}^{\dagger})^{b_{N-1}}\cdots(\hat{c}_{1}^{\dagger})^{b_{1}}|0\rangle_{\rm
F}$ (94)

is transformed by $\hat{g}_{m}$, where $b_{i}=0$ or $1$ with
$(\hat{c}_{i}^{\dagger})^{0}\equiv\hat{1}$ and $|0\rangle_{\rm F}$ denotes the
fermion vacuum such that $\hat{c}_{i}|0\rangle_{\rm F}=0$ for any $i$. We call
a state of the form in Eq. (94) an occupation-basis state. Notice that here we
have adopted the convention for the order of the single-particle local states
$i\ (=1,2,\cdots,N)$ such that the smaller $i$ resides on the further right,
indicating the product $\prod_{i=N}^{1}\cdots$, although any convention for
the order can be used. Then, the occupation-basis state is transformed by
$\hat{g}_{m}$ as

$\displaystyle\hat{g}_{m}|b_{N}b_{N-1}\cdots b_{1}\rangle_{\rm F}$
$\displaystyle=$
$\displaystyle\hat{g}_{m}\prod_{i=N}^{1}(\hat{c}_{i}^{\dagger})^{b_{i}}|0\rangle_{\rm
F}$ $\displaystyle=$
$\displaystyle\prod_{i=N}^{1}(\hat{c}_{m(i)}^{\dagger})^{b_{i}}|0\rangle_{\rm
F}$ $\displaystyle=$
$\displaystyle\prod_{m=N}^{1}(\hat{c}_{m(i)}^{\dagger})^{b_{i}}|0\rangle_{\rm
F}\times{\rm sgn}({\cal R})$ $\displaystyle=$
$\displaystyle\prod_{i=N}^{1}(\hat{c}_{i}^{\dagger})^{b_{m^{-1}(i)}}|0\rangle_{\rm
F}\times{\rm sgn}({\cal R})$ $\displaystyle=$
$\displaystyle\,\,|b_{m^{-1}(N)}b_{m^{-1}(N-1)}\cdots
b_{m^{-1}(1)}\rangle_{\rm F}\times{\rm sgn}({\cal R}),$ (95)

where ${\rm sgn}(\cal R)$ denotes an appropriate sign factor due to the
reordering of the fermion operators, and $m^{-1}(i)$ denotes the number
obtained by applying the inverse permutation $\sigma_{m}^{-1}$ to $i$. More
details on fermionic symmetry operations are briefly reviewed in Appendix B.
An example shown in Fig. 2(b) corresponds to a symmetry operation

$\displaystyle\hat{g}_{m}$
$\displaystyle|b_{6}\,b_{5}\,b_{4}\,b_{3}\,b_{2}\,b_{1}\rangle_{\rm F}$
$\displaystyle=$
$\displaystyle\hat{g}_{m}(\hat{c}_{6}^{\dagger})^{b_{6}}(\hat{c}_{5}^{\dagger})^{b_{5}}(\hat{c}_{4}^{\dagger})^{b_{4}}(\hat{c}_{3}^{\dagger})^{b_{3}}(\hat{c}_{2}^{\dagger})^{b_{2}}(\hat{c}_{1}^{\dagger})^{b_{1}}|0\rangle_{\rm
F}$ $\displaystyle=$
$\displaystyle(\hat{c}_{4}^{\dagger})^{b_{6}}(\hat{c}_{3}^{\dagger})^{b_{5}}(\hat{c}_{2}^{\dagger})^{b_{4}}(\hat{c}_{1}^{\dagger})^{b_{3}}(\hat{c}_{6}^{\dagger})^{b_{2}}(\hat{c}_{5}^{\dagger})^{b_{1}}|0\rangle_{\rm
F}$ $\displaystyle=$
$\displaystyle(\hat{c}_{6}^{\dagger})^{b_{2}}(\hat{c}_{5}^{\dagger})^{b_{1}}(\hat{c}_{4}^{\dagger})^{b_{6}}(\hat{c}_{3}^{\dagger})^{b_{5}}(\hat{c}_{2}^{\dagger})^{b_{4}}(\hat{c}_{1}^{\dagger})^{b_{3}}|0\rangle_{\rm
F}\times{\rm sgn}({\cal R})$ $\displaystyle=$
$\displaystyle|b_{2}\,b_{1}\,b_{6}\,b_{5}\,b_{4}\,b_{3}\rangle_{\rm
F}\times{\rm sgn}({\cal R})$ (96)

with the permutation $\sigma_{m}$ given in Eq. (87). We should note that the
${\rm sgn}({\cal R})$ is not the parity of the permutation $\sigma_{m}$ but
depends on how fermions occupy the single-particle local states.

Let us now discuss quantum circuits for fermionic symmetry operations.
Following the convention of the Jordan-Wigner transformation in Eqs. (53) and
(54), here we assume the correspondence that the $i$th fermionic single-
particle state is occupied (unoccupied) if the state of the $i$th qubit is
$|1\rangle_{i}$ ($|0\rangle_{i}$). The fermion vacuum in the occupation basis
can then be represented in the computational basis as

$|0\rangle_{\rm F}\overset{\rm
JWT}{=}|0\rangle_{N}|0\rangle_{N-1}\cdots|0\rangle_{1}.$ (97)

Note however that, in general, the occupation-basis state $|b_{N}b_{N-1}\cdots
b_{1}\rangle_{\rm F}$ under the Jordan-Wigner transformation is not identical
to the computational basis state $|b_{N}b_{N-1}\cdots
b_{1}\rangle=|b_{N}\rangle_{N}|b_{N-1}\rangle_{N-1}\cdots|b_{1}\rangle_{1}$,
i.e.,

$|b_{N}b_{N-1}\cdots b_{1}\rangle_{\rm F}\overset{\rm
JWT}{\not=}|b_{N}b_{N-1}\cdots b_{1}\rangle$ (98)

because of the sign factor due to the anticommutation relation of fermions,
e.g., $|11\rangle_{\rm
F}=\hat{c}_{2}^{\dagger}\hat{c}_{1}^{\dagger}|0\rangle_{\rm F}\overset{\rm
JWT}{=}-|1\rangle_{2}|1\rangle_{1}=-|11\rangle$. The sign factor due to the
anticommutation relations of fermions on a quantum state can be tracked by
using the fermionic-swap operator Bravyi and Kitaev (2002); Essler  _et al._
(2005); Verstraete  _et al._ (2009); Barthel  _et al._ (2009); Wecker  _et
al._ (2015a); Kivlichan  _et al._ (2018) on a quantum circuit, which does not
depend on the convention adopted for the occupation-basis state. The explicit
form of the fermionic-swap operator in terms of the fermion operators is given
in Appendix B.

The nearest-neighbor fermionic-swap operator $\hat{\cal F}_{i\delta(i)}$
acting on neighboring qubits $i$ and $\delta(i)$ in the computational basis is
defined as

$\displaystyle\hat{\cal
F}_{i\delta(i)}|0\rangle_{i}|0\rangle_{\delta(i)}=|0\rangle_{i}|0\rangle_{\delta(i)},$
(99) $\displaystyle\hat{\cal
F}_{i\delta(i)}|0\rangle_{i}|1\rangle_{\delta(i)}=|1\rangle_{i}|0\rangle_{\delta(i)},$
(100) $\displaystyle\hat{\cal
F}_{i\delta(i)}|1\rangle_{i}|0\rangle_{\delta(i)}=|0\rangle_{i}|1\rangle_{\delta(i)},$
(101) $\displaystyle\hat{\cal
F}_{i\delta(i)}|1\rangle_{i}|1\rangle_{\delta(i)}=-|1\rangle_{i}|1\rangle_{\delta(i)}.$
(102)

The matrix representation of $\hat{\cal F}_{i\delta(i)}$ in the above four
basis states is thus given by

$\hat{\cal F}_{i\delta(i)}\overset{\cdot}{=}\begin{bmatrix}1&0&0&0\\\
0&0&1&0\\\ 0&1&0&0\\\ 0&0&0&-1\end{bmatrix}.$ (103)

Since the controlled-$Z$ (CZ) operator acting on qubits $i$ and $j$ is
represented as

$\widehat{\rm CZ}_{ij}\overset{\cdot}{=}\begin{bmatrix}1&0&0&0\\\ 0&1&0&0\\\
0&0&1&0\\\ 0&0&0&-1\end{bmatrix},$ (104)

the nearest-neighbor fermionic-swap operator can be written as [also see Fig.
3(a)]

$\hat{\cal F}_{i\delta(i)}\overset{\rm JWT}{=}\hat{\cal
S}_{i\delta(i)}\widehat{\rm CZ}_{i\delta(i)}=\widehat{\rm
CZ}_{i\delta(i)}\hat{\cal S}_{i\delta(i)}.$ (105)

Figure 3:  Fermionic-swap gate acting on a pair of (a) nearest-neighbor qubits
$i$ and $\delta(i)=i\pm 1$ and (b) distant qubits $i$ and $j$ (the figure
refers to the case of $|i-j|=5$) is decomposed into a product of the swap and
CZ gates . In (b), the third and fourth equalities follow from the gate
identities shown in the boxes.

The decomposition of the fermionic-swap operator as in Eq. (105) is valid only
when the qubits are nearest neighbored, i.e., $\delta(i)=i\pm 1$. As shown in
Fig. 3(b), the long-range fermionic-swap operator $\hat{\cal F}_{ij}$ for
$j\not=\delta(i)$ is inherently nonlocal and can be implemented as a product
of nearest-neighbor fermionic-swap operators [also see Fig. 2(b)].
Furthermore, using the identities shown in Fig. 3(b), $\hat{\cal F}_{ij}$ can
be expressed by $\hat{\cal S}_{ij}\widehat{\rm CZ}_{ij}$ sandwiched between
two sequences of $|i-j|-1$ CZ gates. The example shown in Fig. 3(b) can be
written as

$\displaystyle\hat{\cal F}_{ij}$ $\displaystyle\overset{\rm
JWT}{=}\left[\prod_{i\lessgtr k\lessgtr j}\widehat{\rm
CZ}_{jk}\right]\hat{f}_{ij}\left[\prod_{i\lessgtr k\lessgtr j}\widehat{\rm
CZ}_{jk}\right],$ (106)

where

$\hat{f}_{ij}\equiv\hat{\cal S}_{ij}\widehat{\rm CZ}_{ij}=\widehat{\rm
CZ}_{ij}\hat{\cal S}_{ij}\overset{\cdot}{=}\begin{bmatrix}1&0&0&0\\\
0&0&1&0\\\ 0&1&0&0\\\ 0&0&0&-1\end{bmatrix}.$ (107)

Note that this matrix is exactly the same as that in Eq. (103), but now qubits
$i$ and $j$ are not necessarily nearest neighbored.

Let us now consider how one can implement a unitary operator of the form
$\exp{\left[-{\rm i}\hat{q}_{ij}\right]}\overset{\rm JWT}{=}\exp{\left[-{\rm
i}\hat{h}_{ij}\hat{Z}_{{\rm JW},ij}\right]}$ in a quantum circuit, where
$\hat{q}_{ij}\overset{\rm JWT}{\equiv}\hat{h}_{ij}\hat{Z}_{{\rm JW},ij}$ is
Hermitian and quadratic in terms of fermion operators, satisfying that
$\hat{\cal F}_{kj}\hat{q}_{ij}\hat{\cal F}_{kj}=\hat{q}_{ik}$, and
$\hat{h}_{ij}$ is a Hermitian operator in the qubit representation acting on
qubits $i$ and $j$, i.e., $[\hat{h}_{ij},\hat{Z}_{{\rm JW},ij}]=0$. Examples
of $\hat{q}_{ij}$ include
$\hat{q}_{ij}=\theta(\hat{c}_{i}^{\dagger}\hat{c}_{j}+{\rm H.c.})$ and
$\hat{q}_{ij}={\rm i}\theta(\hat{c}_{i}^{\dagger}\hat{c}_{j}^{\dagger}-{\rm
H.c.})$ with real $\theta$. First, the equalities $\hat{\cal F}_{kj}\exp[-{\rm
i}\hat{q}_{ij}]\hat{\cal F}_{kj}=\exp[-{\rm i}\hat{q}_{ik}]\overset{\rm
JWT}{=}\exp[-{\rm i}\hat{h}_{ik}\hat{Z}_{{\rm JW},ik}]$ and $\hat{Z}_{{\rm
JW},i\delta(i)}=\hat{I}$ imply that, by using two long-range fermionic-swap
gates ${\cal F}_{\delta(i)j}$, where qubit $\delta(i)$ is located between
qubits $i$ and $j$, we can remove the Jordan-Wigner string $\hat{Z}_{{\rm
JW},ij}$ from the exponent (first equality in Fig. 4):

$\exp\left[-{\rm i}\hat{h}_{ij}\hat{Z}_{{\rm JW},ij}\right]=\hat{\cal
F}_{\delta(i)j}\exp\left[-{\rm i}\hat{h}_{i\delta(i)}\right]\hat{\cal
F}_{\delta(i)j}.$ (108)

Then, by representing the long-range fermionic-swap gates $\hat{\cal
F}_{\delta(i)j}$ as a product of the CZ and the swap gates (second equality in
Fig. 4) and using the identity $\widehat{\rm CZ}_{kl}^{2}=\hat{I}$, we can
cancel the redundant $\widehat{\rm CZ}_{kl}$ gates and obtain (third equality
in Fig. 4)

$\exp\left[-{\rm i}\hat{h}_{ij}\hat{Z}_{{\rm
JW},ij}\right]=\left[\prod_{i\lessgtr k\lessgtr j}\widehat{\rm
CZ}_{jk}\right]\exp\left[-{\rm i}\hat{h}_{ij}\right]\left[\prod_{i\lessgtr
k\lessgtr j}\widehat{\rm CZ}_{jk}\right].$ (109)

A similar strategy for eliminating the redundancy in consecutive Jordan-Wigner
strings has been reported in Ref. Hastings  _et al._ (2015).

Figure 4:  Implementation of the ($|i-j|+1$)-qubit unitary gate $\exp[-{\rm
i}\hat{h}_{ij}\hat{Z}_{{\rm JW},ij}]$ (top left) by a product of two-qubit
unitary gates (bottom). The figure refers to the case of $|i-j|=5$. The second
equality follows from the decomposition of the long-range fermionic-swap gate
shown in Fig. 3(b), and the third equality follows from the identity
$\widehat{\rm CZ}_{kl}^{2}=\hat{I}$. The last equality makes use of the
decomposition of the long-range fermionic-swap gate into a product of the
nearest-neighbor fermionic-swap gates [for a variant of this decomposition,
see the first equality in Fig. 3(b)].

Two remarks are in order. First, by substituting
$\hat{h}_{ij}\propto\hat{X}_{i}\hat{X}_{j}+\hat{Y}_{i}\hat{Y}_{j}$ in Eq.
(109), we can reproduce the quantum circuit for the exponentiated hopping term
of fermions reported previously in Ref. Reiner  _et al._ (2019), i.e., an
exchange-type gate McKay  _et al._ (2016) sandwiched between two sequences of
$|i-j|-1$ CNOT gates. Second, the ($|i-j|+1$)-qubit unitary gate
$\exp\left[-{\rm i}\hat{h}_{ij}\hat{Z}_{{\rm JW},ij}\right]$ can also be
implemented by a product of neighboring two-qubit gates at a cost of
additional $2(|i-j|-1)$ fermionic-swap gates (last equality in Fig. 4). For
example, assuming that $i<\delta(i)=i+1<j$, $\exp\left[-{\rm
i}\hat{h}_{ij}\hat{Z}_{{\rm JW},ij}\right]$ can also be implemented as

$\exp\left[-{\rm i}\hat{h}_{ij}\hat{Z}_{{\rm
JW},ij}\right]=\left[\prod_{k=\delta(i)+1}^{j}\hat{\cal
F}_{k,k-1}\right]\exp\left[-{\rm
i}\hat{h}_{i\delta(i)}\right]\left[\prod_{k=j}^{\delta(i)+1}\hat{\cal
F}_{k,k-1}\right].$ (110)

The resulting quantum circuit first moves the local state at the $j$th qubit
to the qubit next to qubit $i$ by applying the nearest-neighbor fermionic-swap
gates successively, thus keeping track of the antisymmetric nature of fermion
exchange, then it applies the gate $\exp[-{\rm i}\hat{h}_{i\delta(i)}]$ on
qubits $i$ and $\delta(i)$, and finally it moves the local state at the
$\delta(i)$th qubit to the original location, i.e., the $j$th qubit, by
undoing the fermionic-swap operations successively Stoudenmire and White
(2010). We should note here that such a way of exploiting the nearest-neighbor
fermionic-swap gates has been discussed extensively for simulating fermions
under the Jordan-Wigner transformation in Refs. Kivlichan  _et al._ (2018);
Cai (2020); Cade  _et al._ (2020). While Eq. (110) might be preferable for a
real quantum device, depending on the connectivity of qubits and its native
gate set, we employ Eq. (109) for our numerical simulations with classical
computers because of the smaller number of operations.

###  VI.2 Quantum circuit for spin rotation

Figure 5:  (a) A quantum circuit to implement spin rotation ${\rm e}^{-{\rm
i}\alpha\hat{S}_{z}}{\rm e}^{-{\rm i}\beta\hat{S}_{y}}$ for the $4\times 2$
site Fermi-Hubbard model. $f{\cal G}(\beta)$ denotes a fermionic Givens-
rotation gate for distant qubits $i_{\uparrow}$ and $i_{\downarrow}$ defined
in Eq. (116). For evaluation of Eq. (84), ${\rm e}^{-{\rm
i}\alpha\hat{S}_{z}}$ is not required and hence the single-qubit rotation
$R_{z}(\theta)$ should be replaced with identity. (b) A decomposition of a
fermionic Givens-rotation gate $f{\cal G}(\beta)$ into a Givens-rotation gate
${\cal G}(\beta)$ sandwiched with CZ gates, as in Eq. (115). (c) A
decomposition of the Givens-rotation gate ${\cal G}(\beta)$ defined in Eq.
(117). Here, $\hat{R}_{Y}(\theta)={\rm e}^{-{\rm i}\theta\hat{Y}_{i}/2}$
acting on qubit $i$. The matrix representation of ${\cal G}(\beta)$ with the
basis states
$\\{|0\rangle_{i}|0\rangle_{j},|0\rangle_{i}|1\rangle_{j},|1\rangle_{i}|0\rangle_{j},|1\rangle_{i}|1\rangle_{j}\\}$
is also shown.

With the Jordan-Wigner transformation, the local spin operators are
represented as

$\displaystyle\hat{S}_{i}^{x}$ $\displaystyle\overset{\rm
JWT}{=}\frac{1}{4}\left(\hat{X}_{i_{\uparrow}}\hat{X}_{i_{\downarrow}}+\hat{Y}_{i_{\uparrow}}\hat{Y}_{i_{\downarrow}}\right)\hat{Z}_{{\rm
JW},i_{\uparrow}i_{\downarrow}},$ (111) $\displaystyle\hat{S}_{i}^{y}$
$\displaystyle\overset{\rm
JWT}{=}\frac{1}{4}\left(\hat{X}_{i_{\uparrow}}\hat{Y}_{i_{\downarrow}}-\hat{Y}_{i_{\uparrow}}\hat{X}_{i_{\downarrow}}\right)\hat{Z}_{{\rm
JW},i_{\uparrow}i_{\downarrow}},$ (112) $\displaystyle\hat{S}_{i}^{z}$
$\displaystyle\overset{\rm
JWT}{=}\frac{1}{4}\left(\hat{Z}_{i_{\downarrow}}-\hat{Z}_{i_{\uparrow}}\right).$
(113)

Therefore, the rotation operator ${\rm e}^{-{\rm
i}\beta\hat{S}_{y}}=\prod_{i=1}^{L}{\rm e}^{-{\rm i}\beta\hat{S}^{y}_{i}}$ can
be given by a product of

$\displaystyle{\rm e}^{-{\rm i}\beta\hat{S}^{y}_{i}}$
$\displaystyle\overset{\rm JWT}{=}\exp\left[-{\rm
i}\frac{\beta}{4}\left(\hat{X}_{i_{\uparrow}}\hat{Y}_{i_{\downarrow}}-\hat{Y}_{i_{\uparrow}}\hat{X}_{i_{\downarrow}}\right)\hat{Z}_{{\rm
JW},i_{\uparrow},i_{\downarrow}}\right]$ (114)
$\displaystyle=\left[\prod_{i_{\uparrow}<k<i_{\downarrow}}\widehat{{\rm
CZ}}_{i_{\downarrow}k}\right]\hat{\cal
G}_{i_{\uparrow}i_{\downarrow}}\left(\beta\right)\left[\prod_{i_{\uparrow}<k<i_{\downarrow}}\widehat{{\rm
CZ}}_{i_{\downarrow}k}\right]$ (115) $\displaystyle\equiv\widehat{f{\cal
G}}_{i_{\uparrow}i_{\downarrow}}(\beta),$ (116)

where the ${\rm CZ}$ gates in Eq. (115) account for the Jordan-Wigner string
in Eq. (114), as shown for the more general case in Eq. (109), and

$\hat{\cal G}_{ij}(\theta)=\exp\left[{-}{\rm
i}\frac{\theta}{4}\left(\hat{X}_{i}\hat{Y}_{j}-\hat{Y}_{i}\hat{X}_{j}\right)\right]$
(117)

is the Givens-rotation gate for $i\neq j$ Wecker  _et al._ (2015a) whose
matrix representation in the computational basis is given by

$\hat{\cal G}_{ij}(\theta)\overset{\cdot}{=}\begin{bmatrix}1&0&0&0\\\
0&\cos{\tfrac{\theta}{2}}&\sin{\tfrac{\theta}{2}}&0\\\
0&-\sin{\tfrac{\theta}{2}}&\cos{\tfrac{\theta}{2}}&0\\\ 0&0&0&1\end{bmatrix}.$
(118)

For deriving this matrix representation, it is useful to notice that
$(\hat{X}_{i}\hat{Y}_{j}-\hat{Y}_{i}\hat{X}_{j})^{2}=2(1-\hat{Z}_{i}\hat{Z}_{j})$
and
$(1-\hat{Z}_{i}\hat{Z}_{j})(\hat{X}_{i}\hat{Y}_{j}-\hat{Y}_{i}\hat{X}_{j})=2(\hat{X}_{i}\hat{Y}_{j}-\hat{Y}_{i}\hat{X}_{j})$
when $i\neq j$. A more detailed description for a general single-particle
fermion operator is found in Appendix C.2. In Eq. (116), we have defined the
fermionic Givens-rotation gate $\widehat{f{\cal
G}}_{i_{\uparrow}i_{\downarrow}}(\beta)$, which is an extension of the Givens-
rotation gate for fermions.

Similarly, the rotation around the $z$ direction by $\alpha$ can be given by

${\rm e}^{-{\rm i}\alpha\hat{S}^{z}_{i}}\overset{\rm
JWT}{=}\hat{R}_{Z,i_{\uparrow}}\left(-\frac{\alpha}{2}\right)\hat{R}_{Z,i_{\downarrow}}\left(\frac{\alpha}{2}\right),$
(119)

where

$\hat{R}_{Z,i}(\theta)={\rm e}^{-{\rm i}\theta\hat{Z}_{i}/2}.$ (120)

Figure 5 shows a quantum circuit corresponding to the product of rotations
${\rm e}^{-{\rm i}\alpha\hat{S}_{z}}{\rm e}^{-{\rm
i}\beta\hat{S}_{y}}=\prod_{i}{\rm e}^{-{\rm
i}\alpha\hat{S}_{i}^{z}}\prod_{i}{\rm e}^{-{\rm i}\beta\hat{S}_{i}^{y}}$. For
the Givens rotation, we adopt the gate decomposition given in Ref. Jiang  _et
al._ (2018). One can also find another way of implementing the rotation ${\rm
e}^{-{\rm i}\beta\hat{S}_{y}}$ in Refs. Wecker  _et al._ (2015a); Tsuchimochi
_et al._ (2020)

###  VI.3 Quantum circuit for $\eta$ rotation

With the Jordan-Wigner transformation, the local $\eta$ operators are
represented as

$\displaystyle\hat{\eta}_{i}^{x}$ $\displaystyle\overset{\rm JWT}{=}\frac{{\rm
e}^{{\rm
i}\phi_{i}}}{4}\left(\hat{X}_{i_{\uparrow}}\hat{X}_{i_{\downarrow}}-\hat{Y}_{i_{\uparrow}}\hat{Y}_{i_{\downarrow}}\right)\hat{Z}_{{\rm
JW},i_{\uparrow}i_{\downarrow}},$ (121) $\displaystyle\hat{\eta}_{i}^{y}$
$\displaystyle\overset{\rm JWT}{=}-\frac{{\rm e}^{{\rm
i}\phi_{i}}}{4}\left(\hat{X}_{i_{\uparrow}}\hat{Y}_{i_{\downarrow}}+\hat{Y}_{i_{\uparrow}}\hat{X}_{i_{\downarrow}}\right)\hat{Z}_{{\rm
JW},i_{\uparrow}i_{\downarrow}},$ (122) $\displaystyle\hat{\eta}_{i}^{z}$
$\displaystyle\overset{\rm
JWT}{=}-\frac{1}{4}\left(\hat{Z}_{i_{\uparrow}}+\hat{Z}_{i_{\downarrow}}\right).$
(123)

Therefore, the rotation operator ${\rm e}^{-{\rm
i}\beta\hat{\eta}_{y}}=\prod_{i=1}^{N}{\rm e}^{-{\rm
i}\beta\hat{\eta}^{y}_{i}}$ can be given by a product of

$\displaystyle{\rm e}^{-{\rm i}\beta\hat{\eta}^{y}_{i}}$
$\displaystyle\overset{\rm JWT}{=}\exp\left[{\rm i}\frac{{\rm e}^{{\rm
i}\phi_{i}}\beta}{4}\left(\hat{X}_{i}\hat{Y}_{j}+\hat{Y}_{i}\hat{X}_{j}\right)\hat{Z}_{{\rm
JW},i_{\uparrow},i_{\downarrow}}\right]$ (124)
$\displaystyle=\left[\prod_{i_{\uparrow}<k<i_{\downarrow}}\widehat{{\rm
CZ}}_{i_{\downarrow}k}\right]\hat{\cal
B}_{i_{\uparrow}i_{\downarrow}}\left(-{\rm e}^{{\rm
i}\phi_{i}}\beta\right)\left[\prod_{i_{\uparrow}<k<i_{\downarrow}}\widehat{{\rm
CZ}}_{i_{\downarrow}k}\right]$ (125) $\displaystyle\equiv\widehat{f{\cal
B}}_{i_{\uparrow}i_{\downarrow}}(-{\rm e}^{{\rm i}\phi_{i}}\beta),$ (126)

where the ${\rm CZ}$ gates in Eq. (125) account for the Jordan-Wigner string
in Eq. (124), as shown for the more general case in Eq. (109), and

$\hat{\cal B}_{ij}(\theta)=\exp\left[{-}{\rm
i}\frac{\theta}{4}\left(\hat{X}_{i}\hat{Y}_{j}+\hat{Y}_{i}\hat{X}_{j}\right)\right].$
(127)

is the Bogoliubov-transformation gate Wecker  _et al._ (2015a); Jiang  _et
al._ (2018) whose matrix representation in the computational basis is given by

$\hat{\cal
B}_{ij}(\theta)\overset{\cdot}{=}\begin{bmatrix}\cos{\tfrac{\theta}{2}}&0&0&-\sin{\tfrac{\theta}{2}}\\\
0&1&0&0\\\ 0&0&1&0\\\
\sin{\tfrac{\theta}{2}}&0&0&\cos{\tfrac{\theta}{2}}\end{bmatrix}.$ (128)

For deriving this matrix representation, it is useful to notice that
$(\hat{X}_{i}\hat{Y}_{j}+\hat{Y}_{i}\hat{X}_{j})^{2}=2(1+\hat{Z}_{i}\hat{Z}_{j})$
and
$(1+\hat{Z}_{i}\hat{Z}_{j})(\hat{X}_{i}\hat{Y}_{j}+\hat{Y}_{i}\hat{X}_{j})=2(\hat{X}_{i}\hat{Y}_{j}+\hat{Y}_{i}\hat{X}_{j})$
when $i\neq j$. A more detailed description for a general anomalous single-
particle fermion operator is found in Appendix C.3. In Eq. (126), we have
defined the fermionic Bogoliubov-transformation gate $\widehat{f{\cal
B}}_{i_{\uparrow}i_{\downarrow}}(\beta)$, which is an extension of the
Bogoliubov-transformation gate for fermions.

Figure 6:  (a) A quantum circuit to implement $\eta$ rotation ${\rm e}^{-{\rm
i}\beta\hat{\eta}_{y}}{\rm e}^{-{\rm i}\alpha\hat{\eta}_{z}}$ for the $4\times
2$ site Fermi-Hubbard model. $f{\cal B}(\beta)$ denotes a fermionic
Bogoliubov-transformation gate for distant qubits $i_{\uparrow}$ and
$i_{\downarrow}$ defined in Eq. (126). For evaluation of Eq. (84), ${\rm
e}^{-{\rm i}\alpha\hat{\eta}_{z}}$ is not required and hence the single-qubit
rotation $R_{z}(\theta)$ should be replaced with identity. (b) A decomposition
of a fermionic Bogoliubov-transformation gate $f{\cal B}(\beta)$ into a
Bogoliubov-transformation gate ${\cal B}(\beta)$ sandwiched with CZ gates, as
in Eq. (125), assuming that site $i$ represented by qubits $i_{\uparrow}$ and
$i_{\downarrow}$ belongs to $B$ sublattice. When site $i$ represented by
qubits $i_{\uparrow}$ and $i_{\downarrow}$ belongs to $A$ sublattice, $f{\cal
B}(\beta)$ and ${\cal B}(\beta)$ should be replaced with $f{\cal B}(-\beta)$
and ${\cal B}(-\beta)$, respectively. (c) A decomposition of the Bogoliubov-
transformation gate $f{\cal B}(\beta)$ defined in Eq. (127). Here,
$\hat{R}_{Y}(\theta)={\rm e}^{-{\rm i}\theta\hat{Y}_{i}/2}$ acting on qubit
$i$. The matrix representation of ${\cal B}(\beta)$ with the basis states
$\\{|0\rangle_{i}|0\rangle_{j},|0\rangle_{i}|1\rangle_{j},|1\rangle_{i}|0\rangle_{j},|1\rangle_{i}|1\rangle_{j}\\}$
is also shown.

Similarly, the rotation around the $z$ direction by $\alpha$ can be given as

${\rm e}^{-{\rm i}\alpha\hat{\eta}^{z}_{i}}\overset{\rm
JWT}{=}\hat{R}_{Z,i_{\uparrow}}\left(-\frac{\alpha}{2}\right)\hat{R}_{Z,i_{\downarrow}}\left(-\frac{\alpha}{2}\right).$
(129)

Figure 6 shows a quantum circuit corresponding to the product of rotations
${\rm e}^{-{\rm i}\alpha\hat{\eta}_{z}}{\rm e}^{-{\rm
i}\beta\hat{\eta}_{y}}=\prod_{i}{\rm e}^{-{\rm
i}\alpha\hat{\eta}_{i}^{z}}\prod_{i}{\rm e}^{-{\rm
i}\beta\hat{\eta}_{i}^{y}}$. For the Bogoliubov transformation, we adopt the
gate decomposition given in Ref. Jiang  _et al._ (2018). One can also find
another way of implementing the rotation ${\rm e}^{-{\rm
i}\beta\hat{\eta}_{y}}$ in Ref. Wecker  _et al._ (2015a).

##  VII Numerical simulations

In this section, we demonstrate the Krylov-extended SAVQE method by
numerically simulating the two-component Fermi-Hubbard model on the $4\times
2$ cluster under open boundary conditions at half filling, i.e., one fermion
per site. To speed up the calculations, a simple strategy of parallelizing
numerical simulations is employed (see Appendix D).

###  VII.1 Variational state

As a variational state $|\psi(\boldsymbol{\theta})\rangle$ in Eq. (30) for the
ground state of the Fermi-Hubbard model, we first construct a product state of
two-site bonding orbitals for each spin of fermions by applying Hadamard,
Pauli $X$, and CNOT gates to $|0\rangle^{\otimes N}$ in an appropriate manner
[see Fig. 7(a)]. Note that the bonding state is equivalent to one of the Bell
states, $\tfrac{1}{\sqrt{2}}(|0\rangle|1\rangle+|1\rangle|0\rangle)$ (for its
preparation, see Ref. Nielsen and Chuang (2000)). The preparation of the
product state of bonding orbitals is independent of the variational parameters
and hence corresponds to $\hat{W}|0\rangle^{\otimes N}$ in Eq. (30), also
indicated in Fig. 7(a). Note that the total number of fermions for the state
represented by this first part of the quantum circuit is $N/2=L$ (i.e., half
filling) with the same number of up and down fermions, assuming that $L$ is
even. This implies that the expectation values of $\hat{S}_{z}$ and
$\hat{\eta}_{z}$ are both zero. As described below, we will construct a
parametrized part of the quantum circuit for the variational state
$|\psi(\boldsymbol{\theta})\rangle$ that preserves these features, i.e.,
$\langle\psi(\boldsymbol{\theta})|\hat{S}_{z}|\psi(\boldsymbol{\theta})\rangle=\langle\psi(\boldsymbol{\theta})|\hat{\eta}_{z}|\psi(\boldsymbol{\theta})\rangle=0$
for an arbitrary set of variational parameters $\boldsymbol{\theta}$. However,
this does not necessarily imply that the expectation values of ${\hat{S}}^{2}$
and ${\hat{\eta}}^{2}$ are zero. Instead, these expectation values are
generally nonzero, as shown later in Figs. 12(c) and 12(d), for example.

Figure 7:  (a) A quantum circuit for preparing the variational state
$|\psi(\boldsymbol{\theta})\rangle=\prod_{k=N_{\rm
v}}^{1}\hat{U}_{k}(\theta_{k})\hat{W}|0\rangle^{\otimes N}$. The parametrized
gates corresponding to $\hat{U}_{k}(\theta_{k})=\exp(-{\rm i}\hat{\cal
F}_{i_{k}j_{k}}\theta_{k}/2)$ and $\hat{U}_{k}(\theta_{k})=\exp(-{\rm
i}\hat{Z}_{i_{k}}\hat{Z}_{j_{k}}\theta_{k}/2)$ are denoted as e$f$swap and
e$ZZ$, respectively. The qubit numbers indicated in the left most side
correspond to the numbering of qubits in Fig. 8 for the two-component Fermi-
Hubbard model on the $4\times 2$ cluster. (b) A decomposition of the e$f$swap
gate $\hat{U}_{k}(\theta_{k})=\exp(-{\rm i}\hat{\cal
F}_{i_{k}j_{k}}\theta_{k}/2)$, as given in Eq. (130). The lower part of the
panel shows that the parametrized gate $\hat{f}_{ij}(\theta)=\exp(-{\rm
i}\hat{f}_{ij}\theta/2)$ can be expressed as a product of the exponential-swap
gate $\exp(-{\rm i}\hat{\cal S}_{ij}\theta/2)$ and the controlled-phase gate
${\rm CPHASE}_{ij}(\theta)$ [also see Eq. (132)]. (c) A decomposition of the
e$ZZ$ gate $\hat{U}_{k}(\theta_{k})=\exp(-{\rm
i}\hat{Z}_{i_{k}}\hat{Z}_{j_{k}}\theta_{k}/2)$.

In order to construct a parametrized part of the quantum circuit for the
variational state $|\psi(\boldsymbol{\theta})\rangle$ in Eq. (30), let $N_{\rm
hop}$ ($N_{\rm int}$) be the number of hopping (interaction) terms in the
Hamiltonian $\hat{\cal H}$ in Eq. (52). For $1\leqslant k\leqslant N_{\rm
hop}$ in Eq. (30), we apply

$\displaystyle\hat{U}_{k}(\theta_{k})$ $\displaystyle=\exp{(-{\rm i}\hat{\cal
F}_{i_{k}j_{k}}\theta_{k}/2)}$
$\displaystyle=\hat{I}\cos\frac{\theta_{k}}{2}-{\rm i}\hat{\cal
F}_{i_{k}j_{k}}\sin{\frac{\theta_{k}}{2}}$ $\displaystyle\overset{\rm
JWT}{=}\left[\prod_{i_{k}\lessgtr m\lessgtr j_{k}}\widehat{{\rm
CZ}}_{i_{k}m}\right]\hat{f}_{i_{k}j_{k}}\left(\theta_{k}\right)\left[\prod_{i_{k}\lessgtr
m\lessgtr j_{k}}\widehat{{\rm CZ}}_{i_{k}m}\right]$ (130)

to every pair of qubits $i_{k}$ and $j_{k}$ between which the hopping ($t$)
term is present in the Hamiltonian $\hat{\cal H}$. Since the long-range
fermionic-swap operator $\hat{\cal F}_{i_{k}j_{k}}$ is nonlocal, the operator
$\hat{U}_{k}(\theta_{k})$ above operates also onto all qubits between qubits
$i_{k}$ and $j_{k}$ [see Fig. 7(b)]. We refer to this gate as an e$f$swap
gate. In Eq. (130), the parametrized two-qubit gate
$\hat{f}_{ij}\left(\theta\right)$ is defined as

$\hat{f}_{ij}\left(\theta\right)\equiv\exp(-{\rm i}\hat{f}_{ij}\theta/2)$
(131)

and $\hat{f}_{ij}$ is given in Eq. (107). Note that the last equality in Eq.
(130) can be proved simply from Eq. (106) and $\hat{f}_{ij}^{2}=\hat{I}$.
Since $\hat{f}_{ij}=\hat{\cal S}_{ij}\widehat{{\rm CZ}}_{ij}=\hat{\cal
S}_{ij}+(\widehat{{\rm CZ}}_{ij}-\hat{I})$ and $\hat{\cal S}_{ij}$ commutes
with $\widehat{{\rm CZ}}_{ij}$, the parametrized two-qubit gate
$\hat{f}_{ij}\left(\theta\right)$ in Eq. (131) can be given by a product of
the exponential-swap gate $\exp(-{\rm i}\hat{\cal S}_{ij}\theta/2)$ and the
controlled-phase gate ${\rm CPHASE}_{ij}(\theta)\equiv\exp{[-{\rm
i}(\widehat{{\rm CZ}}_{ij}-\hat{I})\theta/2]}$, i.e.,

$\displaystyle\hat{f}_{ij}\left(\theta\right)=$ $\displaystyle\exp(-{\rm
i}\hat{\cal S}_{ij}\theta/2)\ {\rm CPHASE}_{ij}(\theta)$
$\displaystyle\overset{\cdot}{=}$ $\displaystyle\begin{bmatrix}{\rm e}^{-{\rm
i}\theta/2}&0&0&0\\\ 0&\cos{\tfrac{\theta}{2}}&-{\rm
i}\sin{\tfrac{\theta}{2}}&0\\\ 0&-{\rm
i}\sin{\tfrac{\theta}{2}}&\cos{\tfrac{\theta}{2}}&0\\\ 0&0&0&{\rm e}^{-{\rm
i}\theta/2}\end{bmatrix}\begin{bmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\
0&0&0&{\rm e}^{{\rm i}\theta}\end{bmatrix}$ $\displaystyle=$
$\displaystyle\begin{bmatrix}{\rm e}^{-{\rm i}\theta/2}&0&0&0\\\
0&\cos{\tfrac{\theta}{2}}&-{\rm i}\sin{\tfrac{\theta}{2}}&0\\\ 0&-{\rm
i}\sin{\tfrac{\theta}{2}}&\cos{\tfrac{\theta}{2}}&0\\\ 0&0&0&{\rm e}^{+{\rm
i}\theta/2}\end{bmatrix}$ (132)

in the computational basis [also see Fig. 7(b)]. From this matrix
representation, we can indeed confirm directly that
$\hat{f}_{ij}\left(\theta\right)=\hat{I}\cos\frac{\theta}{2}-{\rm
i}\hat{f}_{ij}\sin{\frac{\theta}{2}}$.

For $N_{\rm hop}+1\leqslant k\leqslant N_{\rm hop}+N_{\rm int}$, we apply

$\hat{U}_{k}(\theta_{k})=\exp{(-{\rm
i}\hat{Z}_{i_{k}}\hat{Z}_{j_{k}}\theta_{k}/2)}$ (133)

to every pair of qubits $i_{k}$ and $j_{k}$ between which the interaction
($U_{\rm H}$) term is present in the Hamiltonian $\hat{\cal H}$. As shown in
Fig. 7(c), this gate (eZZ gate) is easily implemented in the circuit.

We define the depth $D$ of the entire quantum circuit for the variational
state $|\psi(\boldsymbol{\theta})\rangle$ in such a way that a single layer of
gates contains the sequences of $N_{\rm hop}+N_{\rm int}$ gates,
$\prod_{k=1}^{N_{\rm hop}+N_{\rm int}}\hat{U}_{k}(\theta_{k})$, and hence the
total number $N_{\rm v}$ of the variational parameters is $N_{\rm v}=D(N_{\rm
hop}+N_{\rm int})$ [see Fig. 7(a)]. Because $\hat{\cal F}_{ij}^{2}=\hat{1}$
and $(\hat{Z}_{i}\hat{Z}_{j})^{2}=\hat{1}$ are satisfied, this parametrized
variational state allows us to use the parameter-shift rules for derivatives
described in Sec. III.2.

###  VII.2 Suzuki-Trotter decomposition

As described in Sec. IV, we employ the QPM to generate the Hamiltonian power
$\hat{\cal H}^{n}$ for the Krylov subspace $\cal U$ in Eq. (37). In the QPM,
the Hamiltonian power $\hat{\cal H}^{n}$ is approximated with controlled
accuracy by a linear combination of ST decomposed time-evolution operators. In
particular, as shown in Eqs. (39) and (40), we employ the second-order
symmetric ST decomposition of the time-evolution operator ${\rm e}^{-{\rm
i}\hat{\cal H}\Delta}$:

${\rm e}^{-{\rm i}\hat{\cal H}\Delta}=\hat{S}_{2}(\Delta)+O(\Delta^{3})$ (134)

with

$\displaystyle\hat{S}_{2}(\Delta)$ $\displaystyle=$ $\displaystyle{\rm
e}^{-{\rm i}\hat{\cal H}_{t}^{A}\Delta/2}{\rm e}^{-{\rm i}\hat{\cal
H}_{t}^{B}\Delta/2}{\rm e}^{-{\rm i}\hat{\cal H}_{t}^{C}\Delta/2}$ (135)
$\displaystyle\quad\quad\times{\rm e}^{-{\rm i}\hat{\cal H}_{U}\Delta}{\rm
e}^{-{\rm i}\hat{\cal H}_{t}^{C}\Delta/2}{\rm e}^{-{\rm i}\hat{\cal
H}_{t}^{B}\Delta/2}{\rm e}^{-{\rm i}\hat{\cal H}_{t}^{A}\Delta/2},$

where $\hat{\cal H}_{t}^{A}$, $\hat{\cal H}_{t}^{B}$, and $\hat{\cal
H}_{t}^{C}$ are the hopping ($t$) terms between sites connected by different
types of bonds and $\hat{\cal H}_{U}$ is the interaction ($U_{\rm H}$) term in
the Fermi-Hubbard Hamiltonian $\hat{\cal H}$, i.e., $\hat{\cal H}=\hat{\cal
H}_{t}^{A}+\hat{\cal H}_{t}^{B}+\hat{\cal H}_{t}^{C}+\hat{\cal H}_{U}$ (see
Fig. 8). This decomposition scheme is used for $\hat{H}_{\rm ST}^{n}(\Delta)$
in Eq. (40) and thus also for $\hat{H}_{\rm ST(1)}^{n}(\Delta)$ in Eq. (41).

Figure 8:  In the ST decomposition, the Fermi-Hubbard Hamiltonian $\hat{\cal
H}$ on the $4\times 2$ cluster is divided into four parts $\hat{\cal
H}_{t}^{A}$, $\hat{\cal H}_{t}^{B}$, $\hat{\cal H}_{t}^{C}$, and $\hat{\cal
H}_{U}$, where $\hat{\cal H}_{t}^{A}$, $\hat{\cal H}_{t}^{B}$, and $\hat{\cal
H}_{t}^{C}$ are the hopping ($t$) terms between sites connected by different
types of bonds (indicated by blue, orange, and green lines, respectively), and
$\hat{\cal H}_{U}$ is the interaction ($U_{\rm H}$) term (indicated by red
curvy lines). Circles represent qubits that are numbered from 1 to 8 for
single-particle states at site $i\,(=1,2,\dots,8)$ with spin up and from 9
to16 for single-particle states at site $i\,(=1,2,\dots,8)$ with spin down.
This numbering of qubits follows the spin-uniform labeling (see Fig. 11).

###  VII.3 Numerical results

In order to assess the accuracy of the Krylov-extended SAVQE, we evaluate the
variational energy

$E_{0}(\boldsymbol{\theta}^{(x)})=\langle\Psi_{\cal
U}(\boldsymbol{\theta}^{(x)})|\hat{\cal{H}}|\Psi_{\cal
U}(\boldsymbol{\theta}^{(x)})\rangle,$ (136)

the fidelity of the ground state

$F(\boldsymbol{\theta}^{(x)})=\left|\langle\Psi_{0}|\Psi_{\cal
U}(\boldsymbol{\theta}^{(x)})\rangle\right|^{2},$ (137)

the expectation value of the total spin squared

$\langle\hat{S}^{2}\rangle_{\boldsymbol{\theta}^{(x)}}=\langle\Psi_{\cal
U}(\boldsymbol{\theta}^{(x)})|\hat{S}^{2}|\Psi_{\cal
U}(\boldsymbol{\theta}^{(x)})\rangle,$ (138)

and the expectation value of the total $\eta$ squared

$\langle\hat{\eta}^{2}\rangle_{\boldsymbol{\theta}^{(x)}}=\langle\Psi_{\cal
U}(\boldsymbol{\theta}^{(x)})|\hat{\eta}^{2}|\Psi_{\cal
U}(\boldsymbol{\theta}^{(x)})\rangle,$ (139)

as a function of the number of the optimization iteration $x$ in Eq. (14).
Here, $|\Psi_{0}\rangle$ is the exact ground state obtained by the Lanczos
exact diagonalization method,

$|\Psi_{\cal U}(\boldsymbol{\theta}^{(x)})\rangle=\frac{|\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta}^{(x)})\rangle}{\sqrt{\langle\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta}^{(x)})|\Psi_{\cal
U}^{(0)}(\boldsymbol{\theta}^{(x)})\rangle}}$ (140)

is the approximated ground state with the variational parameters
$\boldsymbol{\theta}^{(x)}$ at the $x$th optimization iteration, and
$|\Psi_{\cal U}^{(0)}(\boldsymbol{\theta}^{(x)})\rangle$ is given in Eq. (9).
All these quantities except for the fidelity $F(\boldsymbol{\theta}^{(x)})$
can be calculated as in Eq. (50) along with Eq. (84). For the fidelity
calculation, a careful treatment for the normalization factor of the
projection operator is required (see Appendix E). Note also that the
expectation values of $\hat{S}_{z}$ and $\hat{\eta}_{z}$ are both zero by
construction, i.e., $\langle\Psi_{\cal
U}(\boldsymbol{\theta}^{(x)})|\hat{S}_{z}|\Psi_{\cal
U}(\boldsymbol{\theta}^{(x)})\rangle=\langle\Psi_{\cal
U}(\boldsymbol{\theta}^{(x)})|\hat{\eta}_{z}|\Psi_{\cal
U}(\boldsymbol{\theta}^{(x)})\rangle=0$, independently of
$\boldsymbol{\theta}^{(x)}$.

The numerical results shown in this section are all for the Hubbard
interaction $U_{\rm H}/t=4$. Although the learning rate $\tau$ of the
parameter optimization in Eq. (14) can be varied during the optimization
iteration in general, here we fix $\tau=0.025/t$ throughout the optimization
iteration. The parameter $\Delta$ appearing in Eq. (40) for the QPM is set to
be $\Delta=0.05/t$. When the first-order Richardson extrapolation is employed,
the expected systematic error for approximating the Hamiltonian power
$\hat{\cal H}^{n}$ is $O(\Delta^{4})$ [see Eq. (42)]. Indeed, we find that the
errors in $\boldsymbol{H}$ and $\boldsymbol{S}$ are negligible for this
$\Delta$ value.

Let us first examine the vanishing-gradient or the barren-plateau problem
McClean  _et al._ (2018) for our particular study. While the barren-plateau
problem in the context of variational quantum algorithms has been discussed
often with its dependence on the number of qubits, our focus here is on the
comparison between the natural gradient and the gradient of the cost function,
i.e., the expectation value of energy $E_{0}(\boldsymbol{\theta})$, for a
fixed number of qubits. Figure 9 shows the variance of the natural gradient
defined as

$\sigma_{\rm NG}^{2}\equiv\frac{1}{R}\sum_{r=1}^{R}\left[\frac{1}{N_{\rm
v}}\sum_{k=1}^{N_{\rm
v}}[\tilde{\nabla}E_{0}(\boldsymbol{\theta}_{r})]_{k}^{2}\right],$ (141)

where $\boldsymbol{\theta}_{r}=\\{\theta_{r,k}\\}_{k=1}^{N_{\rm v}}$ is a set
of the randomly generated variational parameters, $R$ is the number of the
random instances, and

$\tilde{\nabla}E_{0}(\boldsymbol{\theta}_{r})\equiv[\boldsymbol{G}(\boldsymbol{\theta}_{r})]^{-1}\nabla
E_{0}(\boldsymbol{\theta}_{r})$ (142)

is the natural gradient Amari (1998) of the expectation value of energy
$E_{0}(\boldsymbol{\theta}_{r})$ [see Eq. (14)]. In the calculations, we set
$d_{\cal U}=1$ and the projection operator of the form $\hat{\cal{P}}=\hat{I}$
(i.e., without any symmetry-projection operators), and we estimate
$\sigma_{\rm NG}^{2}$ from $R=64$ random instances. The error bar in Fig. 9
indicates the standard error of the mean. For comparison, we also calculate
the variance of the gradient

$\sigma_{\rm G}^{2}\equiv\frac{1}{R}\sum_{r=1}^{R}\left[\frac{1}{N_{\rm
v}}\sum_{k=1}^{N_{\rm v}}[\nabla
E_{0}(\boldsymbol{\theta}_{r})]_{k}^{2}\right],$ (143)

in the same manner with the same sets of the randomly generated variational
parameters. In Eqs. (141) and (143), it is assumed that the averages of the
natural gradient and the gradient over the random instances vanish.

Figure 9:  Variance of the natural gradient $\sigma_{\rm NG}^{2}$ (squares)
and that of the gradient $\sigma_{\rm G}^{2}$ (circles) as a function of the
number $D$ of layers in a quantum circuit shown in Fig. 7 for the parametrized
variational state $|\psi(\boldsymbol{\theta})\rangle$ without any symmetry
projections. The Fermi-Hubbard model on the $4\times 2$ cluster is considered
and hence there are $N=16$ qubits. The number $N_{\nu}$ of variational
parameters increase linearly in $D$, i.e., $N_{\nu}=28D$ in this case.

In Euclidean space, the natural gradient should coincide with the gradient.
However, as shown in Fig. 9, this is not the case here. Namely, the
variational-parameter space is certainly non Euclidean. Indeed, the steepest-
descent direction of $E_{0}(\boldsymbol{\theta})$ in the variational-parameter
space, which can be seen as a Riemannian manifold where the Fubini-Study
metric tensor $\boldsymbol{G}(\boldsymbol{\theta})$ is attributed at each
point $\boldsymbol{\theta}$, is in general not along $-\nabla
E_{0}(\boldsymbol{\theta})$ but along
$-\tilde{\nabla}E_{0}(\boldsymbol{\theta})=-[\boldsymbol{G}(\boldsymbol{\theta})]^{-1}\nabla
E_{0}(\boldsymbol{\theta})$ (see Appendix A). Even when the natural gradient
is used, the variance decreases exponentially in $D$ until it is saturated for
large $D$, which is similar to the results reported in Ref. McClean  _et al._
(2018). This indicates that the natural gradient cannot solve the barren-
plateau problem. However, it is remarkable to notice in Fig. 9 that the
variance of the natural gradient is more than one order of magnitude larger
than that of the gradient. The larger variance of the natural gradient than
that of the gradient suggests that the natural gradient can alleviate the
barren-plateau problem by capturing the correct steepest-descent direction at
each point $\boldsymbol{\theta}$ in the variational-parameter space, at the
expense of computing the Fubini-Study metric tensor
$\boldsymbol{G}(\boldsymbol{\theta})$.

Figure 10 shows the ground-state energy and the ground-state fidelity as a
function of the number $D$ of the layers of the parametrized gates (see Fig.
7) calculated for the variational states with the full projection operator
$\hat{\cal P}=\hat{\cal P}^{(\eta)}\hat{\cal P}^{(S)}\hat{\cal P}^{(\alpha)}$.
We perform 64 independent calculations with different sets of random initial
parameters $\boldsymbol{\theta}^{(1)}$, and we evaluate these quantities for
the well optimized variational parameters $\boldsymbol{\theta}^{(1000)}$.
Figure 10 shows the results averaged over the 64 independent calculations and
also the best results in terms of the fidelity among the 64 independent
calculations. For comparison, we also show the results obtained for a
Hamiltonian variational ansatz (HVA) Wecker  _et al._ (2015b); Wiersema  _et
al._ (2020), starting also with the initial state $\hat{W}|0\rangle^{\otimes
N}$, which also preserves the spatial, $S$, and $\eta$ symmetry without using
the projection operators. Each layer of the HVA consists of $O(L^{2})$ two
qubit gates due to the Jordan-Wigner string, as in
$|\psi(\boldsymbol{\theta})\rangle$. Further details of the HVA are found in
Appendix F.

Figure 10:  (a) The ground-state energy $E_{0}$ and (b) the ground-state
fidelity $F$ as a function of the number $D$ of layers for various variational
states with the Krylov subspace dimension $d_{\mathcal{U}}=1$ and 2. The
horizontal line in (a) indicates the exact ground-state energy $E_{\rm
exact}(=-13.01250315t)$ obtained by the Lanczos exact-diagonalization method.
For comparison, the results obtained for a HVA are also shown. Here, we
perform 64 different independent calculations with different sets of random
initial parameters $\boldsymbol{\theta}^{(1)}$ and show the results averaged
over these 64 independent calculations (open symbols) and the best results in
terms of the fidelity among these 64 independent calculations (solid symbols).

As shown in Fig. 10(b), the fidelity of $F\approx 0.949$ is achieved at $D=2$
and $d_{\cal U}=1$ on average, and $F\approx 0.991$ is achieved at best. By
increasing the dimension of the subspace to $d_{\cal U}=2$, $F\approx 0.922$
on average and $F\approx 0.980$ at best at $D=1$, and $F\approx 0.994$ on
average and $F\approx 0.999$ at best at $D=2$ are achieved. Accordingly, as
shown in Fig. 10(a), the ground-state energy $E_{0}$ is also systematically
improved with increasing the number $D$ of layers or the Krylov subspace
dimension $d_{\mathcal{U}}$ to attain almost the exact energy already at $D=1$
and $d_{\mathcal{U}}=2$ at best and the exact energy at $D=2$ on average. On
the other hand, the HVA, which also respects all the Hamiltonian symmetry,
requires six layers on average and four layers at best to achieve high
fidelity of 0.9. These results are consistent with the empirical observation
that the more symmetry is broken to restore it by projection, the better
quality a variational wave function acquires Shi  _et al._ (2014).

It should be noted that applying the projection operators and the Hamiltonian
powers in the Krylov-extended SAVQE requires more gates than the HVA, as
summarized in Table 1. Moreover, the number of CZ gates required for the
Jordan-Wigner string depends on the way of mapping fermion indexes (i.e.,
local single-particle states) to qubits. Figure 11 shows two particular
labeling schemes, which we refer to as spin-uniform and spin-alternating
labelings. Generally, the full projection operator requires $O(L^{2})$ more
two-qubit gates, and the Hamiltonian power requires $O(L)$ two-qubit gates
(assuming the Fermi-Hubbard model on the ladder lattice as in Fig. 8), which
are comparable to one to two more layers in the HVA (see Appendix F). As shown
in Table 1, CZ gates for the spin and $\eta$ rotations can be eliminated by
using the spin-alternating labeling, while the number of CZ gates for the
spatial-symmetry operation is double. Furthermore, the number of CZ gates for
the time-evolution operator by the hopping term scales similarly but with
different prefactors in the two different labelings for the Fermi-Hubbard
model on the ladder lattice. Overall, the two labeling schemes are comparable
in terms of the number of CZ gates required for the Fermi-Hubbard model on the
ladder lattice. However, the spin-alternating labeling would be preferred for
the Fermi-Hubbard model on a two-dimensional square lattice.

Table 1:  Order estimation and counting of the number of two-qubit gates required for the symmetry and the time-evolution operations of the two-component Fermi-Hubbard model in a ladder lattice structure under open-boundary conditions such as the one shown in Fig. 1. Two labeling schemes for qubits, spin-uniform and spin-alternating labelings, are considered (see Fig. 11). “$\times 2$” (“$\times 4$”) in the right-most column indicates that the number of the CZ gates required for the Jordan-Wigner (JW) string in the spin-alternating labeling is approximately doubled (quadrupled) as compared with that in the spin-uniform labeling. “Givens”, “Bogoliubov”, and “Exchange” are three different types of two-qubit gates defined in Eqs. (117), (127), and (211), respectively. $L$ is the number of sites and $\hat{\cal{H}}_{t}\equiv\hat{\cal H}-\hat{\cal H}_{U}$ is the hopping term of the Hamiltonian.  unitary operators | two-qubit gates
---|---
| spin-uniform labeling |  | spin-alternating labeling
spatial symmetry $\hat{g}$ | CZ and swap | CZ for JW string |  | CZ and swap | CZ for JW string
| $O(L)$ | $O(L^{2})$ |  | $O(L)$ | $O(L^{2})\times 2$
spin rotation ${\rm e}^{-{\rm i}\hat{S}_{y}\beta}$ | Givens | CZ for JW string |  | Givens | CZ for JW string
| $L$ | $2L(L-1)$ |  | $L$ | $0$
$\eta$ rotation ${\rm e}^{-{\rm i}\hat{\eta}_{y}\beta}$ | Bogoliubov | CZ for JW string |  | Bogoliubov | CZ for JW string
| $L$ | $2L(L-1)$ |  | $L$ | $0$
time evolution by hopping ${\rm e}^{-{\rm i}\hat{\cal H}_{t}\Delta}$ 111The same counting of the number of the CZ gates is applied for the product of exponentiated fermionic-swap operators [defined in Eq. (130)] in the variational state $|\psi(\boldsymbol{\theta})\rangle$ in Eq. (30). | Exchange | CZ for JW string |  | Exchange | CZ for JW string
| $O(L)$ |  $O(L)$ 222 For the Fermi-Hubbard model on a one-dimensional (two-dimensional square) lattice under open-boundary conditions, the number of the CZ gates required is 0 [$O(L^{3/2})$]. |  | $O(L)$ |  $O(L)\times 4$ 333 For the Fermi-Hubbard model on a one-dimensional (two-dimensional square) lattice under open boundary conditions, the number of the CZ gates required is $O(L)$ [$O(L^{3/2})$].
time evolution by interaction ${\rm e}^{-{\rm i}\hat{\cal H}_{U}\Delta}$ | CNOT | CZ for JW string |  | CNOT | CZ for JW string
| $2L$ | 0 |  | $2L$ | 0
Figure 11:  Labeling schemes for qubits. (a) spin-uniform labeling and (b)
spin-alternating labeling. Here, $L=N/2$ is the number of sites.

Now we show the results not only for the variational states with the full
projection operator $\hat{\cal P}=\hat{\cal P}^{(\eta)}\hat{\cal
P}^{(S)}\hat{\cal P}^{(\alpha)}$, but also with various combinations of the
projection operators, $\hat{\cal P}^{(\alpha)}$, $\hat{\cal P}^{(S)}$, and
$\hat{\cal P}^{(\eta)}$, to resolve the efficiency of each projection
operator. Figure 12 shows the results for the quantum circuit of depth $D=1$
and the Krylov-subspace dimension $d_{\cal U}=1$. Each of the results is
average over 64 independent calculations with different sets of random initial
parameters $\boldsymbol{\theta}^{(1)}$ distributed in $[-0.05,0.05]$. As shown
in Fig. 12(b), the fidelity at $x=1$ is essentially zero for all cases with
these 64 different sets of initial parameters $\boldsymbol{\theta}^{(1)}$. As
compared to the results for the bare variational state
$|\psi(\boldsymbol{\theta})\rangle$ without any projection operators, a slight
improvement in the ground-state energy and a substantial improvement in the
ground-state fidelity are found for $\hat{\cal
P}^{(S)}|\psi(\boldsymbol{\theta})\rangle$ and also for $\hat{\cal
P}^{(\alpha)}|\psi(\boldsymbol{\theta})\rangle$. On the other hand, a slight
improvement in the energy and no improvement in the fidelity are found for
$\hat{\cal P}^{(\eta)}|\psi(\boldsymbol{\theta})\rangle$. Among the
variational states with double projection operators, those with the spatial-
symmetry projection, $\hat{\cal P}^{(S)}\hat{\cal
P}^{(\alpha)}|\psi(\boldsymbol{\theta})\rangle$ and $\hat{\cal
P}^{(\eta)}\hat{\cal P}^{(\alpha)}|\psi(\boldsymbol{\theta})\rangle$, give the
better energy than that without the spatial-symmetry projection, $\hat{\cal
P}^{(\eta)}\hat{\cal P}^{(S)}|\psi(\boldsymbol{\theta})\rangle$. However, the
fidelity for $\hat{\cal P}^{(S)}\hat{\cal
P}^{(\alpha)}|\psi(\boldsymbol{\theta})\rangle$ is substantially larger than
that for $\hat{\cal P}^{(\eta)}\hat{\cal
P}^{(\alpha)}|\psi(\boldsymbol{\theta})\rangle$. It is not surprising that the
variational state with the full projection operator, i.e., $\hat{\cal
P}^{(\eta)}\hat{P}^{(S)}\hat{\cal
P}^{(\alpha)}|\psi(\boldsymbol{\theta})\rangle$, achieves the best ground-
state energy and fidelity.

Figure 12:  (a) The ground-state energy $E_{0}(\boldsymbol{\theta}^{(x)})$,
(b) the ground-state fidelity $F(\boldsymbol{\theta}^{(x)})$, (c) the
expectation value of total spin squared
$\langle\hat{S}^{2}\rangle_{\boldsymbol{\theta}^{(x)}}$, and (d) the
expectation value of total $\eta$ squared
$\langle\hat{\eta}^{2}\rangle_{\boldsymbol{\theta}^{(x)}}$ as a function of
the optimization iteration $x$ with various combinations of projection
operators in $|\Psi_{\cal U}(\boldsymbol{\theta})\rangle$ for the Fermi-
Hubbard model on a $4\times 2$ cluster with open boundary conditions at
$U_{\rm H}/J=4$. $E_{\rm exact}$ in (a) is the exact ground-state energy. The
quantum circuit of depth $D=1$ is used for $|\psi(\boldsymbol{\theta})\rangle$
and hence the number $N_{\rm v}$ of variational parameters is 28. The
dimension of the Krylov subspace is $d_{\cal U}=1$. Each result for various
combinations of projection operators in $|\Psi_{\cal
U}(\boldsymbol{\theta})\rangle$ is averaged over 64 independent calculations
started with 64 different sets of random initial parameters
$\boldsymbol{\theta}^{(1)}$, which are common to all calculations.

As expected, $\langle\hat{S}^{2}\rangle_{\boldsymbol{\theta}^{(x)}}$ and
$\langle\hat{\eta}^{2}\rangle_{\boldsymbol{\theta}^{(x)}}$ are exactly zero,
regardless of the values of variational parameters
$\boldsymbol{\theta}^{(x)}$, for the variational states containing the spin-
symmetry projection $\hat{\cal P}^{(S)}$ and the $\eta$-symmetry projection
$\hat{\cal P}^{(\eta)}$, respectively, shown in Figs. 12(c) and 12(d). It is
also found that while
$\langle\hat{\eta}^{2}\rangle_{\boldsymbol{\theta}^{(x)}}$ for the variational
states without $\hat{\cal P}^{(\eta)}$ tends to decrease toward zero with
increasing the optimization iteration $x$,
$\langle\hat{S}^{2}\rangle_{\boldsymbol{\theta}^{(x)}}$ for the variational
states without $\hat{\cal P}^{(S)}$ remains finite. Moreover,
$\langle\hat{S}^{2}\rangle_{\boldsymbol{\theta}^{(x)}}$ for the variational
state $\hat{\cal P}^{(\eta)}|\psi(\boldsymbol{\theta})\rangle$ converges to a
larger value than that for the bare variational state
$|\psi(\boldsymbol{\theta})\rangle$. This can explain the less improved
fidelity for $\hat{\cal P}^{(\eta)}|\psi(\boldsymbol{\theta})\rangle$ in Fig.
12(b), although the energy is improved better than that for the bare
variational state $|\psi(\boldsymbol{\theta})\rangle$ [see Fig. 12(a)], in
accordance with the variational principle.

Figure 13 shows the same results as in Fig. 12 but with the subspace dimension
$d_{\cal U}=2$. Here, the initial variational parameters are set to be the
optimal values obtained with $d_{\cal U}=1$ in Fig. 12. As shown in Figs.
13(a) and 13(b), by expanding the Krylov subspace, both the ground-state
energy and the ground-state fidelity are significantly improved from $d_{\cal
U}=1$, without increasing the number of the variational parameters. Note that
the most noticeable improvement of the energy and the fidelity is achieved
already at the moment when the Krylov subspace is expanded without further
optimizing the variational parameters. The further improvement of these
quantities is also observed with increasing the optimization iteration $x$.
The relative accuracy among the eight different variational states follows
essentially the same trend as in the case with $d_{\cal U}=1$ shown in Fig.
12. However, we should note that the great improvement of
$\langle\hat{S}^{2}\rangle_{\boldsymbol{\theta}^{(x)}}$ and
$\langle\hat{\eta}^{2}\rangle_{\boldsymbol{\theta}^{(x)}}$ is not found even
when the Krylov subspace is expanded to $d_{\cal U}=2$, as shown in Figs.
13(c) and 13(d).

Figure 13:  Same as Fig. 12 but with the Krylov-subspace dimension $d_{\cal
U}=2$. The 64 different sets of optimized variational parameters
$\boldsymbol{\theta}^{(x)}$ at the last optimization iteration $x=10^{3}$
obtained for each variational state with $d_{\cal U}=1$ in Fig. 12 are used as
the initial sets of variational parameters $\boldsymbol{\theta}^{(1)}$ for the
corresponding variational states here. For comparison, the results obtained at
the last optimization iteration $x=10^{3}$ in Fig. 12 are also plotted by
solid symbols at $x=0.7$. The spiky behavior is occasionally observed in the
ground state energy $E_{0}(\boldsymbol{\theta}^{(x)})$ during the parameter
optimization iteration. This occurs for particular sets of calculations when
the learning rate $\lambda$ in Eq. (14) is too large to guarantee the
monotonic decrease of $E_{0}(\boldsymbol{\theta}^{(x)})$ with the optimization
iteration $x$ (for more details, see Appendix A.2).

Figure 14 shows the same results as in Fig. 12 but with the circuit of depth
$D=2$. As shown in Fig. 14(a), the relative error of the ground-state energy
is decreased by an order of magnitude for all the variational states as
compared with that for the variational states with $D=1$ in Fig. 12(a).
Similarly, the ground-state fidelity in Fig. 14(b) is increased almost two
times larger that that for the variational states with $D=1$ in Fig. 12(b). We
can also notice in Fig. 14(c) that no improvement of
$\langle\hat{S}^{2}\rangle_{\boldsymbol{\theta}^{(x)}}$ is made for the
variational state $\hat{\cal P}^{(\eta)}|\psi(\boldsymbol{\theta})\rangle$
against the bare variational state $|\psi(\boldsymbol{\theta})\rangle$, as in
the case of $D=1$ in Fig. 12(c), even though the ground-state energy is indeed
improved. This is consistent with the moderate improvement of the ground-state
fidelity for $\hat{\cal P}^{(\eta)}|\psi(\boldsymbol{\theta})\rangle$. As
shown in Fig. 14(d),
$\langle\hat{\eta}^{2}\rangle_{\boldsymbol{\theta}^{(x)}}$ converges
essentially to zero for all the variational states even without containing the
$\eta$-symmetry projection $\hat{\cal P}^{(\eta)}$. The spatial-symmetry-
projected state $\hat{\cal P}^{(\alpha)}|\psi(\boldsymbol{\theta})\rangle$ is
the best among the variational states with a single projection operator, while
the variational state with the spatial-symmetry and the spin-symmetry
projections, $\hat{\cal P}^{(S)}\hat{\cal
P}^{(\alpha)}|\psi(\boldsymbol{\theta})\rangle$, is the best among those with
double projection operators.

Figure 14:  Same as Fig. 12 but with the quantum circuit of depth $D=2$ and
hence the number $N_{\rm v}$ of variational parameters is 56.