1]Hefei National Research Center for Physical Sciences at the Microscale and
School of Physical Sciences, University of Science and Technology of China,
Hefei 230026, China 2]Shanghai Research Center for Quantum Science and CAS
Center for Excellence in Quantum Information and Quantum Physics, University
of Science and Technology of China, Shanghai 201315, China 3]Hefei National
Laboratory, University of Science and Technology of China, Hefei 230088, China
4]Shanghai Artificial Intelligence Laboratory, Shanghai, 200232, China
5]Center on Frontiers of Computing Studies, Peking University, Beijing 100871,
China 6]School of Computer Science, Peking University, Beijing 100871, China

# Boson sampling enhanced quantum chemistry

Zhong-Xia Shang  Han-Sen Zhong  Yu-Kun Zhang  Cheng-Cheng Yu  Xiao Yuan  Chao-
Yang Lu  Jian-Wei Pan  Ming-Cheng Chen  [  [  [  [  [  [

###### Abstract

In this work, we give a new variational quantum algorithm for solving
electronic structure problems of molecules using only linear quantum optical
systems. The variational ansatz we proposed is a hybrid of non-interacting
Boson dynamics and classical computational chemistry methods, specifically,
the Hartree-Fock method and the Configuration Interaction method. The Boson
part is built by a linear optical interferometer which is easier to realize
compared with the well-known Unitary Coupled Cluster (UCC) ansatz composed of
quantum gates in conventional VQE and the classical part is merely classical
processing acting on the Hamiltonian. We called such ansatzes Boson Sampling-
Classic (BS-C). The appearance of permanents in the Boson part has its
physical intuition to provide different kinds of resources from commonly used
single-, double-, and higher-excitations in classical methods and the UCC
ansatz to exploring chemical quantum states. Such resources can help enhance
the accuracy of methods used in the classical parts. We give a scalable hybrid
homodyne and photon number measurement procedure for evaluating the energy
value which has intrinsic abilities to mitigate photon loss errors and discuss
the extra measurement cost induced by the no Pauli exclusion principle for
Bosons with its solutions. To demonstrate our proposal, we run numerical
experiments on several molecules and obtain their potential energy curves
reaching chemical accuracy.

## Main text

Linear quantum optical systems (LQOS) kok2007linear  have the advantage that
the photons are very clean and more robust to the noise compared with other
architectures. Therefore, LQOS have played important roles in the development
of quantum information science from the first quantum teleportation experiment
bouwmeester1997experimental  to recent demonstrations of quantum advantage
zhong2020quantum ; zhong2021phase ; madsen2022quantum . However, this
advantage comes with the lack of natural interactions between photons as a
price which can create non-linearity that is essential for universal quantum
computing. To enable universal quantum computing on LQOS, one can introduce
non-linearity by feedforward measurements on ancilla qubits such as the KLM
proposal knill2001scheme  and Measurement-based quantum computing
raussendorf2003measurement , which however, is hard to implement in reality
and has made LQOS fall behind ion traps cirac1995quantum  and superconducting
circuits blais2004cavity .

Passive LQOS i.e. LQOS without feedforward measurements could not be used to
make a universal quantum computer is a consensus. But can it show advantages
on certain tasks that are classically hard? Aaronson and Arkhipov
aaronson2011computational  gave a yes answer, which is quite a surprise to the
community. Specifically, they showed sampling the outputs of a Linear Optical
Interferometer (LOI) with Fock state inputs can not likely be simulated by a
classical computer known as the Boson Sampling (BS) problem. Due to its
friendliness to the experiments, the BS problem has made passive LQOS a
promising platform to show early-stage quantum advantages. This milestone was
first finished by Zhong et al zhong2020quantum  recently using Gaussian Boson
Sampling, a variant of BS, and was improved by later works zhong2021phase ;
madsen2022quantum . However, the demonstration of quantum advantage won’t
change the fact that passive LQOS is non-universal. Thus, the next outstanding
challenge is to find specific problems of practical interest that passive LQOS
can solve and have potential advantages. While there have been several
proposals on molecular docking banchi2020molecular , molecular vibrational
spectra huh2015boson  and graph theory arrazola2018using , further efforts are
still required for more applications.

When we look at other platforms in the Noisy Intermediate-Scale Quantum (NISQ)
era preskill2018quantum , solving Electronic Structure Problem (ESP)
mcardle2020quantum  using the Variational Quantum Eigensolver (VQE)
mcclean2016theory ; cerezo2021variational  is one of the most anticipated
applications. For ESP, we aim to solve the ground energy of the corresponding
chemical Hamiltonian to help predict the chemical reaction rates. This can be
done by VQE where a classical optimizer is used to optimize the energy
expectation measured from the output of a shallow quantum circuit ansatz. The
workflow of VQE has some degree of resilience to the noise, which makes
solving ESP using NISQ quantum computers possible. Since solving ESP exactly
requires exponential resources (Full Configuration Interaction (FCI)) using
classical computers, quantum computing may manifest potential advantages.
However, VQE is not friendly to passive LQOS since there is no natural
correspondence from multi-mode Boson Fock space to multi-qubit space and the
quantum circuit language can not be used to describe general operations of
linear interferometers.

In this work, we give a passive LQOS-based VQE algorithm for solving ESP. The
circuit ansatz in this algorithm is a hybrid of non-interacting Boson dynamics
and classical computational chemistry methods. In the following, we will
introduce this algorithm from its physical intuition to its workflow and
implementation. Specifically, we will talk about how to measure the cost
function efficiently by a combination of homodyne and photon number
detections. We will also present several numerical results to show the
properties and the performance of this algorithm.

In quantum chemistry mcardle2020quantum , we are interested in the electronic
structures of molecules. Under the Born-Oppenheimer approximation, treating
the nuclei as classical point charges, the Hamiltonian of a molecule only has
the kinetic energy terms of the electrons, the electrons’ Coulomb interaction
with the nuclei, and the electron-electron Coulomb repulsion left. For ESP, we
aim to find the ground energy of this electronic structure Hamiltonian.
However, this Hamiltonian is currently in the first quantization picture and
we have to take care of the Pauli exclusion principle explicitly. To simplify
the situation, we typically take the Hamiltonian to the second quantization
picture where electrons can be created or annihilated from a set of orbitals
(e.g. the STO-3G basis set used in our numerical examples):

$H^{F}=\sum_{pq}h_{pq}f_{p}^{\dagger}f_{q}+\frac{1}{2}\sum_{pqrs}h_{pqrs}f_{p}^{\dagger}f_{q}^{\dagger}f_{r}f_{s}$
(1)

where $f_{p}^{\dagger}$ and $f_{p}$ are the Fermion creation and annihilation
operations respectively and the $h_{pq}$ and the $h_{pqrs}$ correspond to the
one- and two-electron integrals. The form of the Hamiltonian $H^{F}$ naturally
preserves the number of electrons. If we consider $M$ orbitals as the basis
set and the molecule contains $N$ electrons, then solving ESP corresponds to
finding the smallest eigenenergy in the $N$-electron subspace of dimension
$\binom{M}{N}$ which we will denote as the legal cluster.

The Hartree-Fock (HF) mcardle2020quantum  method is a well-known method to
give an approximate solution to the problem. HF method can be understood as a
VQE procedure using the non-interacting Fermion dynamics
$U_{F}(\vec{\beta})=\exp(-i\sum_{pq}\beta_{pq}f_{p}^{\dagger}f_{q})$ as an
ansatz to minimize the energy:
$|HF\rangle=U_{F}(\vec{\beta}_{hf})|\phi_{0}\rangle,\
\vec{\beta}_{hf}=\text{argmin}_{\vec{\beta}}\langle\phi_{0}|U_{F}^{\dagger}(\vec{\beta})H^{F}U_{F}(\vec{\beta})|\phi_{0}\rangle$
where $|\phi_{0}\rangle$ is an initial product state with $N$ electrons.
$U_{F}$ belongs to the Fermion Gaussian operations that are classically
tractable. Since $U_{F}$ only works as a basis rotation of orbitals and
contains no electron-electron interactions, the HF method can only give a
single-determinant solution whose accuracy is limited. To get more accurate
results, one typically uses $|HF\rangle$ as the reference state
$|\varPhi_{ref}\rangle$ for post-HF methods such as Configuration Interaction
(CI) mcardle2020quantum  and Coupled Cluster (CC) mcardle2020quantum  methods
to give multi-determinant solutions brandas2011advances . Specifically, CI
wavefunction is $(1+T)|\varPhi_{ref}\rangle$ and CC wave function is
$e^{T}|\varPhi_{ref}\rangle$ with $T=T_{1}+T_{2}+...$ contains single, double,
and higher order excitations. When only considering
$T=T_{1}+T_{2}=\sum_{pq}\beta_{pq}f_{p}^{\dagger}f_{q}+\sum_{pqrs}\beta_{pqrs}f_{p}^{\dagger}f_{q}^{\dagger}f_{r}f_{s}$,
we have the CI singles and doubles (CISD) and CC singles and doubles (CCSD).
For quantum computing, the famous Unitary Coupled Cluster (UCC) ansatz
romero2018strategies ; anand2022quantum  used for VQE is the unitary version
of CCSD whose dynamics can be expressed as $e^{T-T^{\dagger}}$. To implement
commonly used UCC Singles and Doubles (UCCSD) on quantum computers, we usually
use the Trotter decomposition to approximate the full dynamics in terms of
local unitary operations, which gives rise to deep circuits intractable in
NISQ devices. The idea of these methods is to introduce additional resources
of interacting Fermion terms to explore states beyond a single determinant in
the legal cluster. The more freedom you have to express states, the higher
accuracy is possible.

Let us briefly introduce BS and make a comparison with Fermion’s version. The
general operations of a multi-mode LOI are non-interacting Boson dynamics that
can be expressed as
$U_{B}(\vec{\alpha})=\exp(-i\sum_{pq}\alpha_{pq}b_{p}^{\dagger}b_{q})$ where
$b_{p}^{\dagger}$ and $b_{p}$ are the Boson creation and annihilation
operations respectively. The (anti-)commutation relations of $b^{\dagger}$ and
$b$ are different from $f^{\dagger}$ and $f$, which reflects the different
behaviors between Bosons and Fermions. This also leads to the difference in
the complexity between the BS and the Fermion Sampling (FS) brod2021bosons .
More specifically, the transition amplitude from an initial product state
$|S\rangle=|s_{1}s_{2}...s_{M}\rangle$ to an output product state
$|T\rangle=|t_{1}t_{2}...t_{M}\rangle$ ($s_{p}$ and $t_{p}$ are natural
numbers for Bosons and 0, 1 for Fermions.) under the non-interacting Boson and
the Fermion dynamics respectively can be expressed as:

$\displaystyle\langle
T|U_{B}(\vec{\alpha})|S\rangle=\frac{\text{Per}(U_{B,ST}(\vec{\alpha}))}{\sqrt{\Pi_{p}s_{p}!\Pi_{p}t_{p}!}},\
\langle T|U_{F}(\vec{\beta})|S\rangle=\text{Det}(U_{F,ST}(\vec{\beta}))$ (2)

where Per and Det denote the matrix permanent and the matrix determinant and
$U_{B,ST}(\vec{\alpha})$ ($U_{F,ST}(\vec{\beta})$) is generated from the mode
transformation matrix $e^{i\alpha}$ ($e^{i\beta}$) corresponding to
$U_{B}(\vec{\alpha})$ ($U_{F}(\vec{\beta})$) brod2021bosons . Since the
computation of permanent is in #P-hard complexity and computation of
determinant are in P complexity, we believe BS is classically intractable but
FS is classically tractable aaronson2011computational . This difference is the
key inspiration of our proposal.

To enable passive LQOS for solving ESP, we use the single-rail encoding
kok2007linear  where 0, 1 photon states $|0\rangle_{B}$ and $|1\rangle_{B}$
represents 0, 1 electron states $|0\rangle_{F}$ and $|1\rangle_{F}$. When we
use the Jordan-Wigner (JW) mapping mcardle2020quantum  to transform Fermion
Hamiltonians to qubit Hamiltonians, this encoding also means $|0\rangle_{B}$
and $|1\rangle_{B}$ represents qubit 0, 1 states $|0\rangle_{Q}$ and
$|1\rangle_{Q}$. The conservation of photon number in passive LQOS naturally
corresponds to the conservation of electron number. As mentioned above, to go
beyond a single determinant and create superpositions in the legal cluster, it
is important to have resources different from the HF method based on non-
interacting Fermion dynamics. The comparison between BS and FS tells us that
permanents in non-interacting Boson dynamics under the single-rail encoding
exactly provide such a resource and are different from the high-order
excitations in CI, CC, and UCC. Also, this resource is classically
intractable. This observation helps us create the following ansatz which we
named Boson Sampling-Classic (BS-C) for VQE:

$\displaystyle|\varPhi_{BS-C}\rangle=V_{C}(\vec{\beta})|\varPhi(\vec{\alpha})\rangle=V_{C}(\vec{\beta})U_{B}(\vec{\alpha})|\varPhi_{ref}\rangle$
(3)

The $U_{B}(\vec{\alpha})$ part is the non-interacting Boson dynamics that
serves as the additional resources to enhance the performances of classical
operations $V_{C}(\vec{\beta})$ such as
$V_{HF}(\vec{\beta})=U_{B}(\vec{\beta})$ in HF that leads to BS-HF and
$V_{CISD}(\vec{\beta})=1+T_{1}+T_{2}$ in CISD that leads to BS-CISD.

The realization of BS-C and the implementation of the measurement of energy
require a detailed discussion. The initial single-determinant reference state
$|\varPhi_{ref}\rangle$ directly corresponds to an initial photon state. The
occupied orbitals correspond to optical modes with single photon inputs and
the virtual orbitals correspond to modes with vacuum inputs.
$U_{B}(\vec{\alpha})$ is a real part realized by a parametrized LOI whose
universal realization can be found in Ref.reck1994experimental .
$V_{C}(\vec{\beta})$ is a virtual part that merely corresponds to classical
processing acting on the Hamiltonian $H^{F}$. Similar ideas have recently been
investigated in Ref.shang2021schr ; mizukami2020orbital . Specifically, we
need to calculate the transformed Hamiltonian
$H^{F}_{T}(\vec{\beta})=V_{C}^{\dagger}(\vec{\beta})H^{F}V_{C}(\vec{\beta})$
which is used to replace $H^{F}$ for JW transformation to get a qubit
Hamiltonian $H^{JW}_{T}(\vec{\beta})$ for measurements. (In current encoding,
$H^{JW}_{T}(\vec{\beta})$ has the same matrix elements as
$H^{F}_{T}(\vec{\beta})$. The reason we add JW mapping is that the Pauli
expression of $H^{JW}_{T}(\vec{\beta})$ can be easier to formulate the
measurement procedure.) Calculating $H^{F}_{T}(\vec{\beta})$ is efficient for
proper $V_{C}(\vec{\beta})$ like HF and CISD since $H^{F}$ in Eq.1 is a local
Fermion Hamiltonian and $V_{HF}(\vec{\beta})$ is merely a basis rotation for
HF and $V_{CISD}(\vec{\beta})$ is also a local operator for CISD. The
resulting $H^{F}_{T}(\vec{\beta})$ is also local, which means there is only a
polynomial number of Fermion terms in $H^{F}_{T}(\vec{\beta})$ and also a
polynomial number of Pauli terms in $H^{JW}_{T}(\vec{\beta})$. The polynomial
number of terms in $H^{JW}_{T}(\vec{\beta})$ is crucial for the measurement
procedure to be scalable which will talked about later. Note that other
classical blocks can also be used as the classical part as long as they
satisfy the conditions on efficient measurements shown later.

For the cost function, we should realize that the value
$\langle\varPhi(\vec{\alpha})|H^{JW}_{T}(\vec{\beta})|\varPhi(\vec{\alpha})\rangle$
doesn’t correspond to the true expectation value of energy and two corrections
are needed. One direct correction is inherited from $V_{C}(\vec{\beta})$. When
$V_{C}(\vec{\beta})$ is non-unitary, the norm of $|\varPhi_{BS-C}\rangle$ will
be changed. The other correction is due to the fact that
$|\varPhi(\vec{\alpha})\rangle$ may have populations outside the encoding
space since there is no Pauli exclusion principle for photons, which will also
influence the length of $|\varPhi_{BS-C}\rangle$. Thus, the true energy value
should be:

$\displaystyle
E(\vec{\alpha},\vec{\beta})=\frac{\langle\varPhi(\vec{\alpha})|H^{JW}_{T}(\vec{\beta})|\varPhi(\vec{\alpha})\rangle}{\langle\varPhi(\vec{\alpha})|V_{C}^{\dagger}(\vec{\beta})V_{C}(\vec{\beta})|\varPhi(\vec{\alpha})\rangle}$
(4)

We should remember that both $H^{JW}_{T}(\vec{\beta})$ and
$V_{C}(\vec{\beta})$ are defined in the encoded Fermion $M$-orbital
$N$-electron Hilbert space of dimension $\binom{M}{N}$ which have zero matrix
element values outside this subspace. For BS-HF, $V_{HF}(\vec{\beta})$ is
unitary in this subspace, $E(\vec{\alpha},\vec{\beta})$ is reduced to
$\langle\varPhi(\vec{\alpha})|H^{JW}_{T}(\vec{\beta})|\varPhi(\vec{\alpha})\rangle/\langle\varPhi(\vec{\alpha})|Q|\varPhi(\vec{\alpha})\rangle$
where $Q$ is the projection operator from Boson $M$-mode $N$-photon Hilbert
space of dimension $\binom{M+N-1}{N}$ to the encoded Fermion $M$-orbital
$N$-electron Hilbert space of dimension $\binom{M}{N}$ (Note that
$|\varPhi(\vec{\alpha})\rangle$ is composed of only $N$-photon states). We
name the value of
$\langle\varPhi(\vec{\alpha})|Q|\varPhi(\vec{\alpha})\rangle$ the projection
ratio. When considering the whole measurement cost of
$E(\vec{\alpha},\vec{\beta})$, the value of the projection ratio is crucial
which we will discuss later.

Having the definition of the cost function in Eq. 4, the question is how to
measure it. According to the transformation rule of the JW mapping, the terms
we need to measure will only be composed of $I$, $Z$, $\sigma_{+}$, and
$\sigma_{-}$. Under the single-rail encoding, different quantum states
correspond to different photon numbers, thus, conventional photon number basis
measurement in dual-encoding is not universal kok2007linear  for measuring
$\sigma_{-}$ and $\sigma_{+}$. On the other hand, the homodyne measurement
d2007homodyne  based on continuous variables in phase space which can be done
experimentally by mixing the target quantum mode with a local oscillator by a
balanced beam splitter and measuring two output modes using two photo-
detectors with post-processing weedbrook2012gaussian  can only be efficient
for evaluating local terms which is not compatible with the JW mapping due to
the existence of $Z$ (See Appendix). Here, a term is local means it has $Z$,
$\sigma_{+}$, and $\sigma_{-}$ act on a limited number of qubits. To resolve
this obstacle, we give a hybrid measurement strategy where for qubits with $I$
and $Z$, we use the photon number measurements, and for qubits with
$\sigma_{+}$, and $\sigma_{-}$ , we use the homodyne measurements. We want to
mention that the introduction of Gaussian measurements will not affect the
hardness of BS chakhmakhchyan2017boson .

To illustrate the workflow, we can consider a term with the form
$H_{i}=Z_{1}Z_{2}\sigma_{+,3}\sigma_{-,4}$, the value $Tr(H_{i}\rho)$ with
$\rho$ a 4-mode optical state can be expressed as:

$\displaystyle Tr(H_{i}\rho)$ $\displaystyle=\sum_{jk}q_{jk}\langle
j|Z_{1}|j\rangle\langle k|Z_{2}|k\rangle$
$\displaystyle\int_{0}^{\pi}\frac{d\phi_{3}d\phi_{4}}{\pi^{2}}\int_{-\infty}^{\infty}dx_{3}dx_{4}p_{jk}(\phi_{3},\phi_{4},x_{3},x_{4})K(\phi_{3},\phi_{4},x_{3},x_{4},\sigma_{+,3}\sigma_{-,4})$
(5)

where $q_{jk}p_{jk}(\phi_{3},\phi_{4},x_{3},x_{4})=\langle jk|\langle
x_{\phi_{3}}x_{\phi_{4}}|\rho|jk\rangle|x_{\phi_{3}}x_{\phi_{4}}\rangle$ is
the joint probability of getting $j$ photons and $k$ photons on the first two
modes and obtaining the results $x_{\phi_{3}}$ and $x_{\phi_{4}}$ when
measuring the quadratures $X_{3,\phi_{3}}$ and $X_{4,\phi_{4}}$ where
$X_{p,\phi_{p}}=(a_{p}^{\dagger}e^{i\phi_{p}}+a_{p}e^{-i\phi_{p}})/2$ and
$K(\phi_{1},...,\phi_{M},x_{1},...,x_{M},\vec{\beta})$ the kernel function
associated with $\sigma_{+,3}\sigma_{-,4}$ which can be efficiently
calculated. Eq.Main text is a typical Monte-Carlo integration which means we
can do repeated samplings to estimate the value. From the operational level,
to get a sample, Eq. Main text means we can first do a photon number
measurement on the first two qubits and then do a homodyne measurement on the
last two qubits. For the homodyne measurement, we need to get $\phi_{3}$ and
$\phi_{4}$ from a uniform distribution, and then measure the corresponding
$X_{3,\phi_{3}}$ and $X_{4,\phi_{4}}$ and get the results $x_{\phi_{3}}$ and
$x_{\phi_{4}}$ which are then used to calculate the value of the corresponding
kernel function. For the measurement cost, we prove that the required number
of samples following the procedure of Eq. Main text to estimate Eq. 4 is
proportional to $\\{2.07317^{k},\epsilon^{-2},m_{H},m_{V},\chi^{-2}\\}$ with
$m_{H}$ the number of terms in $H^{JW}_{T}(\vec{\beta})$, $m_{V}$ the number
of terms in $V_{C}^{\dagger}(\vec{\beta})V_{C}(\vec{\beta})$ (Note that when
$V_{C}^{\dagger}(\vec{\beta})V_{C}(\vec{\beta})$ is reduced to $Q$, we can
simply do pure photon number measurements to estimate the projection ratio.),
$k$ the maximum number of qubits with $\sigma_{-}$ and $\sigma_{+}$ in terms,
$\epsilon$ the desired precision, and $\chi$ the lower bound of the projection
ratio. Note that $k$ has restricted values for HF and CISD methods.

Currently, the above discussions treat LQOS to be ideal, which neglects the
photon loss error in real devices. Photon loss as the prominent noise channel
is non-negligible, especially for our proposal where photon loss will lead
BS-C VQE searching outside the legal subspace. Thus, besides the corrections
in Eq. 4, we also need corrections for photon loss. Since we use a hybrid
measurement strategy which seems unable to correct the photon loss error by
commonly used post selections at first glance, however, we now show a
procedure to correct the photon loss error whose reason can be found in the
appendix. To illustrate the concrete correction procedure, without loss of
generality, we can consider an $H_{i}$ only has one $\sigma_{+}$ and one
$\sigma_{-}$ as an example. We have $N$ single-photon sources as input to
match the correct number of electrons. First, we do the mentioned hybrid
measurements. Here, since $\sigma_{+}\sigma_{-}=|10\rangle\langle 01|$, we
require to detect $N-1$ photons on $I$ and $Z$ modes with at most one photon
on each mode to proceed with homodyne detections. We can do repeated such
measurements to obtain a raw expectation value of $H_{i}$ denoted as $\langle
H_{i}\rangle_{raw}$. Next, we need a normalization factor which can be
obtained by doing repeated pure photon number measurements on all modes many
times. Suppose among these measurements, we have $n_{1}$ times of getting $N$
photons in total on all modes with at most one photon on each mode and $n_{2}$
times of getting $N-1$ photons in total on $I$ and $Z$ modes with at most one
photon on each mode. And among those $n_{1}$ events, we have $n_{3}$ times
where there is one photon in two non-diagonal modes and $N-1$ photons in
diagonal modes. The true value is then corrected by
$\frac{n_{3}}{n_{1}}\frac{n_{2}}{n_{3}}\langle H_{i}\rangle_{raw}$. This value
has already taken the projection ratio into consideration. For other $H_{i}$,
this procedure can be easily generalized. The efficiency of this error
mitigation method merely depends on the level of photon loss and has no
additional non-scalable cost.

We have introduced all components of our proposal and the whole workflow of
the algorithm which we named BS-C VQE has been summarized in Fig.1. Note that
the parameters $\vec{\alpha}$ and $\vec{\beta}$ are optimized simultaneously.
We also explain the physical intuitions of BS-C VQE in Fig.2a-b. We want to
emphasize that the additional resources provided by BS are non-trivial. First,
under our encoding, since BS and FS are different, the non-interacting Boson
dynamics in LOI leads to a different transition probability distribution from
non-interacting Fermion dynamics since permanents and determinants in Eq. 2
are fundamentally different. Thus, unlike HF, LOI belongs to the multi-
determinant type. This “interacting” part can be understood as a multi-
determinant reference state brandas2011advances  generator for the followed
classical part and the Boson nature of these resources indicates they have
small overlaps with those classical methods originated from electron
excitations. Thus, when the optimization procedure is converged, we should
expect to obtain more accurate ground energy estimations than using the
classical methods alone. Second, the numerator of
$E(\vec{\alpha},\vec{\beta})$ can be expanded as

$\displaystyle\langle\varPhi(\vec{\alpha})|H^{JW}_{T}(\vec{\beta})|\varPhi(\vec{\alpha})\rangle=\sum_{ij}H^{JW}_{T}(\vec{\beta})_{ij}\text{Per}(U_{B,0i}(\vec{\alpha}))\text{Per}(U_{B,0j}(\vec{\alpha}))$
(6)

where $0$ denotes $|\varPhi_{ref}\rangle$. The permanents appeared in this
expression indicates the additional resources provided by BS are classically
intractable, which is similar to UCCSD whose unitarity makes itself
classically intractable. To have potential quantum advantages, having such
resources is the basic requirement. We want to mention that when compared with
pure classical algorithms, performance enhancements by BS come with an
additional sampling cost and the level of enhancements can vary for different
molecules. If, for example, we have BS-CISD comparable with CISDT in terms of
accuracy, and the sampling cost is smaller than the computational cost of
CISDT, a real quantum advantage is achieved.

When considering the realizations, BS-C VQE has its advantage over the UCCSD.
Unlike UCCSD where we need a locally connected quantum circuit of depth around
$\mathcal{O}(M^{3})$ anand2022quantum  to build in the general case where all
excitations are allowed, BS-C only needs a $\mathcal{O}(M)$-time LOI
reck1994experimental  which has the same scaling as the quantum implementation
of HF ansatz $U_{F}(\vec{\beta})$ of depth around $\mathcal{O}(M)$
google2020hartree  but is with the speed of light as the propagation speed.
Thus, combined with the efficient error mitigation method we will introduce
below against the photon loss, BS-C VQE is much more hardware-efficient.
What’s more, the LQOS for our proposal has the advantage of the sampling rate
compared with other platforms, which is highly demanded for variational
quantum algorithms since massive repeated measurements are required. Current
state-of-the-art experiments show that the sampling rate of LQOS can reach
around $10^{6-8}Hz$ wang2019boson  whereas the sampling rate is around
$10^{3-4}Hz$ for superconducting qubits wu2021strong  and $10^{0-1}Hz$ for ion
traps pino2021demonstration  and neutral atoms bluvstein2022quantum . We
summarize these in Fig.2c.

We demonstrate the performance of BS-C VQE by running numerical experiments on
several molecules and drawing their potential energy curves in Fig.3. We test
molecules including LiH, BeH2, and H4 with line geometry and H4 with square
geometry. All molecules use the STO-3G basis set and the active space
reduction method mcardle2020quantum  is used for LiH, BeH2. The resulting JW
Hamiltonians are 6, 8, 8, and 8 qubits. We add data of FCI, HF, and CISD as
references. We also show the projection ratios corresponding to the BS-C VQE
solutions. We can see BS-C has enhancements over their classical part and can
reach the chemical accuracy ($1.6\times 10^{-3}$ Hartree) at all regions. We
use BS-HF for LiH, where we can see that while we only use non-interacting
Boson and non-interacting Fermion, the results are still rather accurate. The
other three molecules however, are actually strongly correlated at some bond
length where HF performs so poor that even BS can not give enough enhancements
to reach the chemical accuracy. Thus, we use BS-CISD instead.

Now, we want to talk about the projection ratio. The projection ratio
$\langle\varPhi(\vec{\alpha})|Q|\varPhi(\vec{\alpha})\rangle$ has great
influence to the total measurement cost of estimating
$E(\vec{\alpha},\vec{\beta})$. Our statistic analysis shows that the mean
squared error (MSE) containing both bias and variance of estimating
$E(\vec{\alpha},\vec{\beta})$ will increase by $\mathcal{O}(\langle
Q\rangle^{-2})$ times compared with the ideal case with projection ratio equal
to one. Thus, to make the measurements of $E(\vec{\alpha},\vec{\beta})$ truly
scalable, the projection ratio can not be exponentially small, which is also
the reason for the existence of sign problem loh1990sign  in quantum Monte-
Carlo methods. We can increase the projection ratios from two aspects. The
first can be understood from Fig.4a-b. We use BS-HF for the LiH molecule as a
direct demonstration where we can see while all experiments finally give
energy reached the chemical accuracy, they have totally different projection
ratios. This means there are physically different states
$|\varPhi(\vec{\alpha})\rangle$ which can be treated as the same states in the
encoding space. Thus, unlike the conventional VQE where there is only one
minimum energy, there may be multi-solutions in our algorithm. To help the
optimization converge to those with high projection ratios, a penalty term
$-\lambda\langle Q\rangle$ can be introduced to the cost function which helps
solutions with higher projection ratios have lower cost function values.
Another method is to only allow occupied-to-virtual photon excitations and
virtual-to-virtual photon excitations to reduce the populations outside the
encoding space, which is the method we used in numerical experiments. Note
that this can be further adjusted by only allowing those to preserve the spin
and orbital symmetries anand2022quantum .

These methods however are not enough since there may be situations where even
the highest projection ratio among the solutions is still small. We can
understand this by defining the ratio $r(M,N)=\binom{M}{N}/\binom{M+N-1}{N}$
between the dimension of Fermion Hilbert space and the dimension of Boson
Hilbert space. Obviously, $r(M,N)$ is a fair judge for the projection ratio.
If we fix the number of orbitals $M$ and enlarge the number of electrons $N$,
$r(M,N)$ will decrease. In contrast, if we fix the number of electrons $N$ and
enlarge the number of orbitals $M$, $r(M,N)$ will increase. Thus, to avoid low
projection ratios, $M$ should be large enough compared to $N$. Luckily, while
$N$ depends on the number of electrons in the molecule, the number of orbitals
has no restrictions and can be arbitrarily large. In Fig.4c, we run BS-HF for
H2 with different choices of basis sets including STO-3G, 6-31G, and 6-311G.
The resulting number of qubits (spin orbitals) is 4, 8, and 12 respectively
with the number of electrons fixed as 2. We can see that as the number of
orbitals gets larger, the average projection ratios of the experiments that
reach the chemical accuracy gets larger as well. More rigorously, we proved
that if $M=\eta N^{2}$, $r(M,N)$ will have a lower bound $e^{1/\eta}$, which
indicates the empirical requirements to make the measurements scalable. Note
that such a relation between $M$ and $N$ is also a believed requirement for BS
to be classically intractable lund2017quantum ; aaronson2011computational .

In summary and discussion, we gave a novel hybrid quantum-classical ansatz
using only passive LQOS for running VQE to solve ESP. The permanents appeared
in non-interacting Boson dynamics is a new type of resource for exploring the
legal cluster which can help enhance the accuracy of quantum chemistry. Thus,
BS-C VQE is both hardware-efficient with only an LQI and chemical-inspired for
its physical intuitions. We use BS-HF and BS-CISD as two concrete settings and
show their performances by numerical results. We want to mention that the
classical methods are not restricted to HF and CI. The introduction of the
hybrid measurement strategy with the methods against photon loss error and low
projection ratios helps to have a scalable measurement cost. Both the theory
and the simulations showed that our proposal may open a road for quantum
applications using passive LQOS in the NISQ era. Also, this proposal is not
restricted to what we presented here, by the developments of LQOS, to get
better accuracy, one can use entangled reference states brandas2011advances
or introduce intermediate measurements for more freedom and higher accuracy.
It is worth noticing that there are still several open questions on this
algorithm that need to be further investigated such as the level of
performance enhancements of this algorithm on larger molecules and better
measurement strategies such as developing efficient tomography methods with
pure homodyne or photon number measurements olivares2019quantum  and adopting
the ideas of recently developed quantum learning methods huang2020predicting .

We used Python packages including Qiskit cross2018ibm , OperFermion
mcclean2020openfermion , PySCF sun2018pyscf , and The Walrus gupt2019walrus
for parts of our simulations. The example codes can be found online
shang2023bfvqe .

More details can be found in the Appendix.

$\mathbf{Achknowledgement:}$ This work was supported by the National Natural
Science Foundation of China (No. 91836303 and No. 11805197), the National Key
R-D Program of China, the Chinese Academy of Sciences, the Anhui Initiative in
Quantum Information Technologies, and the Science and Technology Commission of
Shanghai Municipality (2019SHZDZX01). H.-S. Z. is supported by National Key
R$\&$D Program of China (NO.2022ZD0160100) Shanghai Committee of Science and
Technology (Grant No. 21DZ1100100). X. Y. and Y.-K. Z. are supported by the
National Natural Science Foundation of China Grant (No. 12175003 and No.
12361161602), NSAF (Grant No. U2330201). The authors would like to thank Z-H
Chen and Y-H Deng for inspiring discussions and S Aaronson and P Coveney for
their fruitful comments.

Figure 1: The framework of BS-C VQE. In BS-C VQE encoding, the 0, 1 occupation
of Boson, the 0, 1 of the qubit, and the 0, 1 occupation of Fermion are
equivalent. Given the electronic structure Hamiltonian of a molecule, the task
is to obtain its ground energy. Typically, we choose a suitable basis set and
do Hamiltonian reduction if necessary to obtain the second quantization
Hamiltonian $H^{F}$ and the HF reference state $|\varPhi_{ref}\rangle$. In
BS-C VQE, $H^{F}$ will be transformed by classical operations like HF and
CISD. The transformed Hamiltonian is then mapped to qubit Hamiltonian by JW
mapping. The reference state gives the information on the initial photon
source preparation. The occupied orbitals correspond to optical modes with
single-photon input and the virtual orbitals correspond to modes with vacuum
input. Such an initial state is then passed to parametrized non-interacting
Boson dynamics i.e. LOI. The output state is used for the hybrid measurements
to obtain energy expectation value which as the cost function is fed to the
classical optimizer to update the parameters $\vec{\alpha}$ and $\vec{\beta}$
for the next iteration until converged. Figure 2: Properties of BS-C ansatz.
a-b: In the Schrödinger picture, we can see both the non-interacting Boson
dynamics and the classical operations are acted on the reference state.
Different operations give different resources to explore the legal cluster.
UCCSD and classical methods like HF and CISD give resources for electron
excitations with various orders. In contrast, permanents appearing in non-
interacting Boson dynamics is a different type of resource that results from
the differences between Bosons and Fermions. Also, since permanent is
classically intractable, it is a quantum resource that may manifest potential
advantages. Thus, the main idea of BS-C is to use such quantum resources to
enhance the accuracy of classical methods. When the classical method is HF, we
have BS-HF. When the classical method is CISD, we have BS-CISD. Note that the
classical methods are not restricted to these two. c: BS-C has the attractive
advantage of requiring much shallower depth than UCCSD. Also, the LQOS for our
proposal can have a much higher sampling rate than other platforms. Figure 3:
Numerical experiments on several molecules (LiH, BeH2, and H4 with line
geometry and H4 with square geometry). The potential energy of molecules as
functions of bond length (Li-H bond, all Be-H bonds, all H-H bonds, and all
H-H bounds) is shown on the left. Results of using BS-C, HF, CISD, and FCI
methods are presented. Note that for LiH, CISD is the same as FCI. BS-HF is
for LiH and BS-CISD is for the others. The numbers attached to BS-C points are
the corresponding projection ratios. Absolute energy differences from FCI are
shown on the right. The shaded grey region represents the area within the
chemical accuracy. The classical optimizer we used is L-BFGS-B and each point
of BS-C is chosen from the best result among 10 runs. The bond length unit is
Angstrom and the energy unit is Hartree. Figure 4: Projection ratio of BS-C
VQE. a: LiH molecule as an example. We fix the bond length as 1.0 and
repeatedly run BS-C VQE. We see that all presented tests have achieved
chemical accuracy. However, as explained, their corresponding projection
ratios are totally different. b: In BS-C VQE encoding, we only utilize 0 and 1
occupation dimensions, and all higher photon occupations are redundant.
However, after the action of LOI, part of the output state will unavoidably
jump to these redundant dimensions. This leads to an interesting point of BS-C
VQE that two physically different states can be treated as the same states in
the encoding space. Thus, when optimizing BS-C VQE, it is possible that we can
get multiple sets of parameters that correspond to the ground state. c: The
number of orbitals and the number of electrons have a crucial influence on the
projection ratio. The smaller of the number of electrons compared to the
number of orbitals, the higher the projection ratio will appear. We test the
H2 molecule with different choices of basis set to show this phenomenon. The
number of electrons is 2 and the orbital number is 4 for STO-3G, 8 for 6-31G,
and 12 for 6-311G. The error bars are collected from the tests that achieve
chemical accuracy. We can easily see the larger number of orbitals gives
higher projection ratios. We use BS-HF for LiH and H2.

## References

  *   * (1) P. Kok, W.J. Munro, K. Nemoto, T.C. Ralph, J.P. Dowling, G.J. Milburn, Linear optical quantum computing with photonic qubits.  Reviews of modern physics 79(1), 135 (2007)
  * (2) D. Bouwmeester, J.W. Pan, K. Mattle, M. Eibl, H. Weinfurter, A. Zeilinger, Experimental quantum teleportation.  Nature 390(6660), 575–579 (1997)
  * (3) H.S. Zhong, H. Wang, Y.H. Deng, M.C. Chen, L.C. Peng, Y.H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, et al., Quantum computational advantage using photons.  Science 370(6523), 1460–1463 (2020)
  * (4) H.S. Zhong, Y.H. Deng, J. Qin, H. Wang, M.C. Chen, L.C. Peng, Y.H. Luo, D. Wu, S.Q. Gong, H. Su, et al., Phase-programmable gaussian boson sampling using stimulated squeezed light.  Physical review letters 127(18), 180,502 (2021)
  * (5) L.S. Madsen, F. Laudenbach, M.F. Askarani, F. Rortais, T. Vincent, J.F. Bulmer, F.M. Miatto, L. Neuhaus, L.G. Helt, M.J. Collins, et al., Quantum computational advantage with a programmable photonic processor.  Nature 606(7912), 75–81 (2022)
  * (6) E. Knill, R. Laflamme, G.J. Milburn, A scheme for efficient quantum computation with linear optics.  nature 409(6816), 46–52 (2001)
  * (7) R. Raussendorf, D.E. Browne, H.J. Briegel, Measurement-based quantum computation on cluster states.  Physical review A 68(2), 022,312 (2003)
  * (8) J.I. Cirac, P. Zoller, Quantum computations with cold trapped ions.  Physical review letters 74(20), 4091 (1995)
  * (9) A. Blais, R.S. Huang, A. Wallraff, S.M. Girvin, R.J. Schoelkopf, Cavity quantum electrodynamics for superconducting electrical circuits: An architecture for quantum computation.  Physical Review A 69(6), 062,320 (2004)
  * (10) S. Aaronson, A. Arkhipov, in _Proceedings of the forty-third annual ACM symposium on Theory of computing_ (2011), pp. 333–342
  * (11) L. Banchi, M. Fingerhuth, T. Babej, C. Ing, J.M. Arrazola, Molecular docking with gaussian boson sampling.  Science advances 6(23), eaax1950 (2020)
  * (12) J. Huh, G.G. Guerreschi, B. Peropadre, J.R. McClean, A. Aspuru-Guzik, Boson sampling for molecular vibronic spectra.  Nature Photonics 9(9), 615–620 (2015)
  * (13) J.M. Arrazola, T.R. Bromley, Using gaussian boson sampling to find dense subgraphs.  Physical review letters 121(3), 030,503 (2018)
  * (14) J. Preskill, Quantum computing in the nisq era and beyond.  Quantum 2, 79 (2018)
  * (15) S. McArdle, S. Endo, A. Aspuru-Guzik, S.C. Benjamin, X. Yuan, Quantum computational chemistry.  Reviews of Modern Physics 92(1), 015,003 (2020)
  * (16) J.R. McClean, J. Romero, R. Babbush, A. Aspuru-Guzik, The theory of variational hybrid quantum-classical algorithms.  New Journal of Physics 18(2), 023,023 (2016)
  * (17) M. Cerezo, A. Arrasmith, R. Babbush, S.C. Benjamin, S. Endo, K. Fujii, J.R. McClean, K. Mitarai, X. Yuan, L. Cincio, et al., Variational quantum algorithms.  Nature Reviews Physics 3(9), 625–644 (2021)
  * (18) E.J. Brandas, J.R. Sabin, _Advances in quantum chemistry_ (Academic Press, 2011)
  * (19) J. Romero, R. Babbush, J.R. McClean, C. Hempel, P.J. Love, A. Aspuru-Guzik, Strategies for quantum computing molecular energies using the unitary coupled cluster ansatz.  Quantum Science and Technology 4(1), 014,008 (2018)
  * (20) A. Anand, P. Schleich, S. Alperin-Lea, P.W. Jensen, S. Sim, M. Díaz-Tinoco, J.S. Kottmann, M. Degroote, A.F. Izmaylov, A. Aspuru-Guzik, A quantum computing view on unitary coupled cluster theory.  Chemical Society Reviews (2022)
  * (21) D.J. Brod, Bosons vs. fermions–a computational complexity perspective.  Revista Brasileira de Ensino de Física 43 (2021)
  * (22) M. Reck, A. Zeilinger, H.J. Bernstein, P. Bertani, Experimental realization of any discrete unitary operator.  Physical review letters 73(1), 58 (1994)
  * (23) Z.X. Shang, M.C. Chen, X. Yuan, C.Y. Lu, J.W. Pan, Schr$\backslash$" odinger-heisenberg variational quantum algorithms.  arXiv preprint arXiv:2112.07881 (2021)
  * (24) W. Mizukami, K. Mitarai, Y.O. Nakagawa, T. Yamamoto, T. Yan, Y.y. Ohnishi, Orbital optimized unitary coupled cluster theory for quantum computer.  Physical Review Research 2(3), 033,421 (2020)
  * (25) G.M. D’Ariano, L. Maccone, M.F. Sacchi, in _Quantum Information With Continuous Variables of Atoms and Light_ (World Scientific, 2007), pp. 141–158
  * (26) C. Weedbrook, S. Pirandola, R. García-Patrón, N.J. Cerf, T.C. Ralph, J.H. Shapiro, S. Lloyd, Gaussian quantum information.  Reviews of Modern Physics 84(2), 621 (2012)
  * (27) L. Chakhmakhchyan, N.J. Cerf, Boson sampling with gaussian measurements.  Physical Review A 96(3), 032,326 (2017)
  * (28) G.A. Quantum, Collaborators*†, F. Arute, K. Arya, R. Babbush, D. Bacon, J.C. Bardin, R. Barends, S. Boixo, M. Broughton, B.B. Buckley, et al., Hartree-fock on a superconducting qubit quantum computer.  Science 369(6507), 1084–1089 (2020)
  * (29) H. Wang, J. Qin, X. Ding, M.C. Chen, S. Chen, X. You, Y.M. He, X. Jiang, L. You, Z. Wang, et al., Boson sampling with 20 input photons and a 60-mode interferometer in a 1 0 14-dimensional hilbert space.  Physical review letters 123(25), 250,503 (2019)
  * (30) Y. Wu, W.S. Bao, S. Cao, F. Chen, M.C. Chen, X. Chen, T.H. Chung, H. Deng, Y. Du, D. Fan, et al., Strong quantum computational advantage using a superconducting quantum processor.  Physical review letters 127(18), 180,501 (2021)
  * (31) J.M. Pino, J.M. Dreiling, C. Figgatt, J.P. Gaebler, S.A. Moses, M. Allman, C. Baldwin, M. Foss-Feig, D. Hayes, K. Mayer, et al., Demonstration of the trapped-ion quantum ccd computer architecture.  Nature 592(7853), 209–213 (2021)
  * (32) D. Bluvstein, H. Levine, G. Semeghini, T.T. Wang, S. Ebadi, M. Kalinowski, A. Keesling, N. Maskara, H. Pichler, M. Greiner, et al., A quantum processor based on coherent transport of entangled atom arrays.  Nature 604(7906), 451–456 (2022)
  * (33) E. Loh Jr, J. Gubernatis, R. Scalettar, S. White, D. Scalapino, R. Sugar, Sign problem in the numerical simulation of many-electron systems.  Physical Review B 41(13), 9301 (1990)
  * (34) A.P. Lund, M.J. Bremner, T.C. Ralph, Quantum sampling problems, bosonsampling and quantum supremacy.  npj Quantum Information 3(1), 15 (2017)
  * (35) S. Olivares, A. Allevi, G. Caiazzo, M.G. Paris, M. Bondani, Quantum tomography of light states by photon-number-resolving detectors.  New Journal of Physics 21(10), 103,045 (2019)
  * (36) H.Y. Huang, R. Kueng, J. Preskill, Predicting many properties of a quantum system from very few measurements.  Nature Physics 16(10), 1050–1057 (2020)
  * (37) A. Cross, in _APS March meeting abstracts_ , vol. 2018 (2018), pp. L58–003
  * (38) J.R. McClean, N.C. Rubin, K.J. Sung, I.D. Kivlichan, X. Bonet-Monroig, Y. Cao, C. Dai, E.S. Fried, C. Gidney, B. Gimby, et al., Openfermion: the electronic structure package for quantum computers.  Quantum Science and Technology 5(3), 034,014 (2020)
  * (39) Q. Sun, T.C. Berkelbach, N.S. Blunt, G.H. Booth, S. Guo, Z. Li, J. Liu, J.D. McClain, E.R. Sayfutyarova, S. Sharma, et al., Pyscf: the python-based simulations of chemistry framework.  Wiley Interdisciplinary Reviews: Computational Molecular Science 8(1), e1340 (2018)
  * (40) B. Gupt, J. Izaac, N. Quesada, The walrus: a library for the calculation of hafnians, hermite polynomials and gaussian boson sampling.  Journal of Open Source Software 4(44), 1705 (2019)
  * (41) Example codes.  https://github.com/ustcszx/BS-C-VQE (2023)
  * (42) M. Paini, A. Kalev, D. Padilha, B. Ruck, Estimating expectation values using approximate quantum states.  Quantum 5, 413 (2021)
  * (43) G. Van Kempen, L. Van Vliet, Mean and variance of ratio estimators used in fluorescence ratio imaging.  Cytometry: The Journal of the International Society for Analytical Cytology 39(4), 300–305 (2000)

##  Appendix A Homodyne measurement

First, consider a single mode d2007homodyne . An operator $H$ can be expressed
in the coherent state basis:

$\displaystyle
H=\int\frac{d^{2}\alpha}{\pi}Tr[HD^{\dagger}(\alpha)]D(\alpha)=\int_{0}^{\pi}\frac{d\phi}{\pi}\int_{-\infty}^{\infty}dr\frac{|r|}{4}Tr(He^{irX_{\phi}})e^{-irX_{\phi}}$
(7)

where $D(\alpha)=e^{(\alpha b^{\dagger}-\alpha^{*}b)}$ is the displacement
operator and $\alpha=-ire^{i\phi}/2$. The quadrature operator $X_{\phi}$ is
defined as $X_{\phi}=(b^{\dagger}e^{i\phi}+be^{-i\phi})/2$. The expectation
value of $H$ under an optical state $\rho$ can then be expressed as:

(8)

The probability $p(\phi,x)=\langle x_{\phi}|\rho|x_{\phi}\rangle$ is obtained
by measuring $\langle X_{\phi}|\rho|X_{\phi}\rangle$. We see that this is a
standard Monte-Carlo estimation. The kernel function $K(\phi,x,H)$ is defined
as:

(9)

When $H=|n+\lambda\rangle\langle n|$, the kernel function has an analytic
expression:

$\displaystyle K(\phi,x,|n+\lambda\rangle\langle$ $\displaystyle
n|)=2e^{-i\lambda\phi}\sqrt{\frac{n!}{(n+\lambda)!}}e^{-x^{2}}\sum_{\nu=0}^{n}\frac{(-1)^{\nu}}{\nu!}\left(\begin{array}[]{c}n+\lambda\\\
n-\nu\end{array}\right)$ (11)
$\displaystyle\times(2\nu+\lambda+1)!\operatorname{Re}\left[(-i)^{\lambda}\mathcal{D}_{-(2\nu+\lambda+2)}(-2ix)\right]$
(12)

where $\mathcal{D}_{l}(x)$ the parabolic cylinder function that can be easily
calculated. In our algorithm, we will need the kernel functions of
$\sigma_{+}$,$\sigma_{-}$, and $Z$. According to Eq. 11, we can calculate the
range of the kernel functions of $\sigma_{+}$ and $\sigma_{-}$ which are among
$[-2.07317,2.07317]$ and the kernel functions of $Z$ which is among
$[-2.92345,5.33333]$.

For multi-mode situations, the generalization is natural:

$\displaystyle\langle H\rangle$
$\displaystyle=\int_{0}^{\pi}\frac{d\phi_{1}d\phi_{2}...d\phi_{M}}{\pi^{M}}\int_{-\infty}^{\infty}dx_{1}dx_{2}...dx_{M}$
$\displaystyle
p(\phi_{1},\phi_{2},...,\phi_{M},x_{1},x_{2},...,x_{M})K(\phi_{1},\phi_{2},...,\phi_{M},x_{1},x_{2},...,x_{M},H)$
(13)

where:

$\displaystyle K(\phi_{1},\phi_{2},...,$
$\displaystyle\phi_{M},x_{1},x_{2},...,x_{M},H)=$
$\displaystyle\int_{-\infty}^{\infty}dr_{1}dr_{2}...dr_{M}\Pi_{m=1}^{M}\frac{|r_{m}|}{4}Tr(He^{ir_{m}(X_{\phi_{m}}-x_{m})})$
(14)

We can see from Eq.A that the multi-mode kernel function is a simple product
of the single-mode kernel function if $H$ is a product operator, which makes
the complexity of calculating the kernel function scalable. Also, we don’t
need to calculate all values of all kernel functions in advance. We only need
to do a measurement followed by a kernel function calculation. Thus, the
resource for calculating the kernel function also has linear scaling with the
number of measurements.

The variance of using Eq. A to estimate the expectation value of $H$ depends
on the kernel function. When a kernel function has the range $[a,b]$, then the
variance of Eq. A is bounded by $\frac{(b-a)^{2}}{4}$. Since the kernel
functions of $\sigma_{+}$,$\sigma_{-}$, and $Z$ have ranges larger than
$[-2,2]$ and terms we want to evaluate in the cost function obtained by the
J-W mapping is global, the kernel functions of these terms can have
exponentially large ranges and thus pure homodyne measurements will be non-
scalable, which leads to the hybrid measurement strategy introduced in the
main text. On the other hand, since
$p(\phi_{1},\phi_{2},...,\phi_{M},x_{1},x_{2},...,x_{M})$ has no dependence on
the concrete form of $H$, if all terms in a Hamiltonian are local, then the
same set of samples can be used to estimate all expectation values of these
terms and will greatly reduce the sampling complexity. In fact, this well-
known homodyne tomography method can be seen as a natural generalization of a
recently developed random measurement protocol in qubit systems
paini2021estimating .

##  Appendix B Measurement complexity

In this section, we will estimate the measurement complexity of estimating:

$\displaystyle
E(\vec{\alpha},\vec{\beta})=\frac{\langle\varPhi(\vec{\alpha})|H^{JW}_{T}(\vec{\beta})|\varPhi(\vec{\alpha})\rangle}{\langle\varPhi(\vec{\alpha})|V_{C}^{\dagger}(\vec{\beta})V_{C}(\vec{\beta})|\varPhi(\vec{\alpha})\rangle}$
(15)

assuming there is no photon loss error. Before the analysis, we will introduce
the estimator that will be used.

  * •

Estimator for $\frac{E[X]}{E[Y]}$: When $X$ and $Y$ are two independent random
variables, the value of the ratio of their expectations $\frac{E[X]}{E[Y]}$
can be estimated by an asymptotically unbiased estimator
$\frac{\overline{X}}{\overline{Y}}$ where $\overline{X}$ and $\overline{Y}$
are the averages of $X$ and $Y$. It has been shown in Ref.van2000mean  that
the expectation and the variance of this estimator are:

$\displaystyle E[\frac{\overline{X}}{\overline{Y}}]$
$\displaystyle\approx\frac{E[X]}{E[Y]}+\frac{E[X]}{E[Y]^{3}}Var[\overline{Y}]$
(16) $\displaystyle Var[\frac{\overline{X}}{\overline{Y}}]$
$\displaystyle\approx\frac{E[X]^{2}Var[\overline{Y}]+E[Y]^{2}Var[\overline{X}]}{E[Y]^{4}}$
(17)

In our case, the numerator $X$ corresponds to $\langle H^{JW}_{T}\rangle$
short for
$\langle\varPhi(\vec{\alpha})|H^{JW}_{T}(\vec{\beta})|\varPhi(\vec{\alpha})\rangle$
and the denominator $Y$ corresponds to $\langle V_{C}^{\dagger}V_{C}\rangle$
short for
$\langle\varPhi_{ref}|V_{C}^{\dagger}(\vec{\beta})V_{C}(\vec{\beta})|\varPhi_{ref}\rangle$.
For the numerator, we assume $H^{JW}_{T}(\vec{\beta})$ has the form
$H^{JW}_{T}(\vec{\beta})=\sum_{i=1}^{m_{H}}h_{i}H_{i}$ with $H_{i}$ containing
$\sigma_{+}$, $\sigma_{-}$, $Z$ and $I$. For each $P_{i}$, we do the photon
number measurements on qubits with $Z$ and $I$ and homodyne measurements on
qubits with $\sigma_{+}$ and $\sigma_{-}$. Suppose each term has at most
$k_{h}$ $\sigma_{+}$ and $k_{h}$ $\sigma_{-}$, since the kernel functions of
$\sigma_{+}$ and $\sigma_{-}$ have ranges among $[-2.07317,2.07317]$, the
remaining $Z$ and $I$ have diagonal elements -1 and 1, the whole variance is
then bounded by $2.07317^{4k_{h}}$. If we do repeated measurements for
$N_{H}/m_{H}$ times for $P_{i}$, the resulting variance is bounded by
$m_{H}2.07317^{4k_{h}}/N_{H}$. Following this, the total variance of the
numerator has the bound:

$\displaystyle Var[\overline{\langle
H^{JW}_{T}\rangle}]\leq\sum_{i=1}^{m_{H}}\frac{2.07317^{4k_{h}}m_{H}h_{i}^{2}}{N_{H}}$
(18)

For the denominator, similarly, if we assume
$V_{C}^{\dagger}V_{C}=\sum_{i=1}^{m_{V}}g_{i}V_{i}$ and each term has at most
$k_{v}$ $\sigma_{+}$ and $k_{v}$ $\sigma_{-}$. Then, if we do measurements
$N_{V}/m_{V}$ times for each $P_{i}$, the total variance of the denominator
has the bound:

$\displaystyle Var[\overline{\langle
V_{C}^{\dagger}V_{C}\rangle}]\leq\sum_{i=1}^{m_{V}}\frac{2.07317^{4k_{v}}m_{V}g_{i}^{2}}{N_{V}}$
(19)

Assume the projection ratio $\langle Q\rangle$ short for
$\langle\varPhi(\vec{\alpha})|Q|\varPhi(\vec{\alpha})\rangle$ has a lower
bound: $\langle Q\rangle\geq\chi$ and define $\langle
H^{JW}_{T}\rangle_{r}=\langle H^{JW}_{T}\rangle/\langle Q\rangle$ and $\langle
V_{C}^{\dagger}V_{C}\rangle_{r}=\langle V_{C}^{\dagger}V_{C}\rangle/\langle
Q\rangle$. By putting Eq.19 and Eq.18 into Eq.16 and Eq.17, we obtain that the
bias of estimating $E(\vec{\alpha},\vec{\beta})$ is bounded by:

$\displaystyle Bias[E(\vec{\alpha},\vec{\beta})]^{2}$
$\displaystyle\leq\left(\frac{\langle H^{JW}_{T}\rangle}{\langle
V_{C}^{\dagger}V_{C}\rangle}-\frac{\overline{\langle
H^{JW}_{T}\rangle}}{\overline{\langle
V_{C}^{\dagger}V_{C}\rangle}}\right)^{2}\approx\left(\frac{\langle
H^{JW}_{T}\rangle}{\langle
V_{C}^{\dagger}V_{C}\rangle^{3}}Var[\overline{\langle
V_{C}^{\dagger}V_{C}\rangle}]\right)^{2}$
$\displaystyle=\left(\frac{E(\vec{\alpha},\vec{\beta})}{\langle
V_{C}^{\dagger}V_{C}\rangle^{2}}\sum_{i=1}^{m_{V}}\frac{2.07317^{4k_{v}}m_{V}g_{i}^{2}}{N_{V}}\right)^{2}$
$\displaystyle=\left(\frac{E(\vec{\alpha},\vec{\beta})}{\chi^{2}\langle
V_{C}^{\dagger}V_{C}\rangle_{r}^{2}}\frac{2.07317^{4k_{v}}m_{V}}{N_{V}}\sum_{i=1}^{m_{V}}g_{i}^{2}\right)^{2}$
$\displaystyle\leq\frac{2.07317^{4k_{v}}m_{V}^{2}(\sum_{i=1}^{m_{V}}g_{i}^{2})^{2}\|H^{JW}_{T}\|_{2}^{2}}{\chi^{4}N_{V}^{2}}$
(20)

and the variance of estimating $E(\vec{\alpha},\vec{\beta})$ is bounded by:

$\displaystyle Var[E(\vec{\alpha},\vec{\beta})]$
$\displaystyle\leq\frac{\langle H^{JW}_{T}\rangle^{2}Var[\overline{\langle
V_{C}^{\dagger}V_{C}\rangle}]+\langle
V_{C}^{\dagger}V_{C}\rangle^{2}Var[\overline{\langle
H^{JW}_{T}\rangle}]}{\langle V_{C}^{\dagger}V_{C}\rangle^{4}}$
$\displaystyle=\frac{E(\vec{\alpha},\vec{\beta})^{2}\sum_{i=1}^{m_{V}}\frac{2.07317^{4k_{v}}m_{V}g_{i}^{2}}{N_{V}}+\sum_{i=1}^{m_{H}}\frac{2.07317^{4k_{h}}m_{H}h_{i}^{2}}{N_{H}}}{\langle
V_{C}^{\dagger}V_{C}\rangle^{2}}$
$\displaystyle=\frac{E(\vec{\alpha},\vec{\beta})^{2}\frac{2.07317^{4k_{v}}m_{V}}{N_{V}}\sum_{i=1}^{m_{V}}g_{i}^{2}+\frac{2.07317^{4k_{h}}m_{H}}{N_{H}}\sum_{i=1}^{m_{H}}h_{i}^{2}}{\langle
V_{C}^{\dagger}V_{C}\rangle^{2}}$
$\displaystyle\leq\frac{2.07317^{4k_{v}}m_{V}(\sum_{i=1}^{m_{V}}g_{i}^{2})\|H^{JW}_{T}\|_{2}^{2}}{\chi^{2}N_{V}}+\frac{2.07317^{4k_{h}}m_{H}\sum_{i=1}^{m_{H}}h_{i}^{2}}{\chi^{2}N_{H}}$
(21)

In the above derivation, we use the fact that for HF and CISD, $\langle
V_{C}^{\dagger}V_{C}\rangle_{r}\geq 1$. The mean squared error is defined as
the sum of Eq.B and Eq.B. From these results, we can see that to make the
measurement scalable, we need $\chi$ is not exponentially small and $k_{v}$
and $k_{h}$ are restricted values, which is true for HF and CISD.

##  Appendix C Error mitigation against photon loss errors

Consider an optical state with N photon in M modes ($M>N$), photon loss can be
regarded as a photon number exchange between the system and environment, as
shown below.

$|\psi\rangle=\sum_{i=0}^{N}c_{i}|\psi_{N-i}\rangle|i\rangle_{e}$ (22)

In this equation, $|\psi\rangle$ represents the whole state including the
system and environment. $|\psi_{N-i}\rangle$ represents the $N-i$ photon state
component in the system with $i$ photons escaping to the environment and
$c_{i}$ is a coefficient determined by the efficiency of the experimental
system. Since the number of electrons in a molecule is fixed, we need to post-
select the state $|\psi_{N}\rangle$. Note that
$|\psi_{N}\rangle=d_{F}|\psi_{N,F}\rangle+...$ where $|\psi_{N,F}\rangle$ is
the component living in the encoding Fermion space. Thus, the true expectation
value of a term $H_{i}$ should be:

$\langle\psi_{N,F}|H_{i}|\psi_{N,F}\rangle$ (23)

Given the feature of the Hamiltonian after the J-W mapping, as explained in
the main text, terms of the Hamiltonian have non-diagonal $\sigma_{+}$ and
$\sigma_{-}$ for certain modes and diagonal $Z$ and $I$ for other modes. Non-
diagonal terms, for example, $\sigma_{+}\sigma_{-}=|10\rangle\langle 01|$,
strongly restrict the form of state that gives rise to non-zero contribution
to the expectation of the term of Hamiltonian. This can be seen by an
elaborate example where $H_{i}$ has $k$ $\sigma_{-}$ and $k$ $\sigma_{+}$.
Non-zero contributions can only be given by an optical state that has $k$
photons in non-diagonal modes, and due to the reason of the electron number
conservation, $N-k$ photons should be in the diagonal modes. Therefore, in the
operational level, we can first estimate
$\langle\phi_{N-k,T}|H_{i}|\phi_{N-k,T}\rangle$ where $|\phi_{N-k,T}\rangle$
means $N-k$ photons in diagonal modes and $T$ means the true state where there
is no more than 1 photon in each mode according to the fermion encoding. In
hybrid measurements, we first directly measure the photon number in diagonal
modes, so, for this condition, only when $N-k$ photons are measured with no
more than 1 photon in these modes, we are allowed to go to the homodyne
measurement for non-diagonal modes. After this post-selection procedure, we
are measuring the state
$|\phi_{N-k}\rangle=\alpha|\phi_{N-k,T}\rangle+\beta|\phi_{N-k,W}\rangle$,
which includes the state we desired and also includes the component in the
wrong subspace including those with more than one photon in a non-diagonal
mode and those with less than $k$ photons in all 2$k$ non-diagonal modes.
Luckily, for the undesired state $|\phi_{N-k,W}\rangle$, it has no non-zero
contributions to the expectation value and the hybrid measurements let us to
estimate $\langle
H_{i}\rangle_{raw}=\langle\phi_{N-k}|H_{i}|\phi_{N-k}\rangle=|\alpha|^{2}\langle\phi_{N-k,T}|H_{i}|\phi_{N-k,T}\rangle$.
Thus, we only need a correction factor $1/|\alpha|^{2}$. Currently,
$|\phi_{N-k,T}\rangle$ is not equal to $|\psi_{N,F}\rangle$ since it has other
$N$-photon components. Suppose we have
$|\psi_{N,F}\rangle=\gamma|\phi_{N-k,T}\rangle+\sum_{i}\kappa_{i}|\phi_{i\neq
N-k,T}\rangle$, since other components also have zero contributions to the
expectation value, we can simply add another correction factor $|\gamma|^{2}$.
Thus, the true result is:

$\langle\psi_{N,F}|H_{i}|\psi_{N,F}\rangle=\frac{|\gamma|^{2}}{|\alpha|^{2}}\langle\phi_{N-k}|H_{i}|\phi_{N-k}\rangle$
(24)

To get these factors, we can add an additional measurement procedure where we
do pure repeated photon number measurements for many times. Suppose during
these measurements, we recorded $n_{1}$ times of events that detected $N$
photons in all and no more than 1 photon in each mode, $n_{2}$ times of events
that detect $N-k$ photon in diagonal modes with no more than 1 photon in these
modes, and $n_{3}$ events among those $n_{1}$ ones where $N-k$ photons in the
diagonal modes and $k$ photons in the non-diagonal modes. Then we have
$|\gamma^{2}|=\frac{n_{3}}{n_{1}}$ and $|\alpha|^{2}=\frac{n_{3}}{n_{2}}$
which can be used for the true expectation value estimations:

$\langle\psi_{N,F}|H_{i}|\psi_{N,F}\rangle\approx\frac{n_{3}}{n_{1}}\frac{n_{2}}{n_{3}}\langle
H_{i}\rangle_{raw}$ (25)

whose total sampling complexity can be estimated by a similar procedure as
above in the measurement complexity section.

##  Appendix D Projection ratio

Figure 5: $r(M,N)$ as functions of $N$. The number of orbitals is set to be
$M=\eta N^{2}$ with $\eta=0.5$, $\eta=1$, and $\eta=2$.

The projection ratio highly depends on the ratio
$r(M,N)=\binom{M}{N}/\binom{M+N-1}{N}$ between the dimension of Fermion
Hilbert space and the dimension of Boson Hilbert space. $r(M,N)$ describes the
portion of encoding space (Fermion) in physical space (Boson). In the
following, we will prove that if $M=\eta N^{2}$, $r(M,N)$ will have a lower
bound $e^{1/\eta}$.

First, $r(M,N)$ can be re-expressed as:

$\displaystyle r(M,N)=\frac{M}{M+N-1}\frac{M-1}{M+N-2}...\frac{M-N+1}{M}$ (26)

Thus, we have:

$\displaystyle r(M,N)>\left(\frac{M-N+1}{M}\right)^{N}$ (27)

If $M=\eta N^{2}$, we have:

$\displaystyle\left(\frac{M-N+1}{M}\right)^{N}=\frac{1}{\left(\frac{\eta
N^{2}}{\eta N^{2}-N+1}\right)^{N}}=\frac{1}{\left(1+\frac{N-1}{\eta
N^{2}-N+1}\right)^{N}}$ (28)

The RHS of Eq.28 has the limit:

$\displaystyle\lim_{N\rightarrow\infty}\frac{1}{\left(1+\frac{N-1}{\eta
N^{2}-N+1}\right)^{N}}=\lim_{N\rightarrow\infty}\left(\frac{1}{\left(1+\frac{1}{\eta
N}\right)^{\eta N}}\right)^{\frac{1}{\eta}}=e^{\frac{1}{\eta}}$ (29)

Note that the RHS of Eq.28 decreases monotonically as $N$ grows, which can be
seen from Fig.5. Thus, we have proved our conclusion.

##  Appendix E Molecule information in Fig.3

All molecules use the STO-3G basis set.

For the LiH molecule, we choose the active orbital $\\{1,2,5\\}$ with 2 active
electrons. The orbital labels are ordered from low energy to high energy.

For the BeH2 molecule, the geometry is a line. We choose the active orbital
$\\{1,2,5,6\\}$ with 4 active electrons when the bond length is smaller than
1.9 and choose the active orbital $\\{1,2,3,6\\}$ with 4 active electrons when
the bond length is larger than 1.9. The orbital labels are ordered from low
energy to high energy.

For the H4 molecules with line geometry and square geometry, changing of bond
length in the main text means changing all H-H bonds at the same time. There
are a total of 4 orbitals and 4 electrons. No reduction method is used.

Each orbital above contains two spin orbitals.

For the systemic introduction of the active space methods, see
Ref.mcardle2020quantum .

The above orbital labels can be directly used in PySCF sun2018pyscf  and
Qiskit cross2018ibm  packages.