#  Encoding strongly-correlated many-boson wavefunctions on a photonic quantum
computer: application to the attractive Bose-Hubbard model

Saad Yalouz<EMAIL_ADDRESS>Theoretical Chemistry, Vrije Universiteit,
De Boelelaan 1083, NL-1081 HV, Amsterdam, The Netherlands  Instituut-Lorentz,
Universiteit Leiden, P.O. Box 9506, 2300 RA Leiden, The Netherlands  Bruno
Senjean<EMAIL_ADDRESS>ICGM, Univ Montpellier, CNRS, ENSCM, Montpellier,
France  Instituut-Lorentz, Universiteit Leiden, P.O. Box 9506, 2300 RA Leiden,
The Netherlands  Vedran Dunjko<EMAIL_ADDRESS>Universiteit
Leiden, P.O. Box 9506, 2300 RA Leiden, The Netherlands

###### Abstract

Variational quantum algorithms (VQA) are considered as some of the most
promising methods to determine the properties of complex strongly correlated
quantum many-body systems, especially from the perspective of devices
available in the near term. In this context, the development of efficient
quantum circuit ansatze to encode a many-body wavefunction is one of the keys
for the success of a VQA. Great efforts have been invested to study the
potential of current quantum devices to encode the eigenstates of fermionic
systems, but little is known about the encoding of bosonic systems. In this
work, we investigate the encoding of the ground state of the (simple but rich)
attractive Bose–Hubbard model using a Continuous-Variable (CV) photonic-based
quantum circuit. We introduce two different ansatz architectures and
demonstrate that the proposed continuous variable quantum circuits can
efficiently encode (with a fidelity higher than 99%) the strongly correlated
many-boson wavefunction with just a few layers, in all many-body regimes and
for different number of bosons and initial states. Beyond the study of the
suitability of the ansatz to approximate the ground states of many-boson
systems, we also perform initial evaluations of the use of the ansatz in a
variational quantum eigensolver algorithm to find it through energy
minimization. To this end we also introduce a scheme to measure the
Hamiltonian energy in an experimental system, and study the effect of sampling
noise.

##  I Introduction

In recent years, great efforts have been deployed to study the potential of
noisy intermediate-scale quantum (NISQ) computers to tackle tasks that are
hard to for classical computers. One of the most emblematic research lines in
the field is the development of near-term algorithms to determine the
properties of complex quantum many-body systems. Variational quantum circuits
have shown to be very efficient and promising for encoding complex
wavefunctions of various kinds of systems. Applications of these circuits
range from molecular systems in quantum chemistry [1, 2, 3, 4, 5, 6, 7, 8, 9,
10, 11], interacting-spins and electrons models from condensed matter [12, 13,
14] or vibrational Hamiltonians [15], to cite but a few.

In practice, the efficiency of a variational quantum algorithm (VQA) depends
on the expressibility of the quantum ansatz used to describe the wavefunction
of the system under study. So far, this question has been extensively
investigated in the case of qubit-based quantum circuits. Examples are the
unitary coupled-cluster ansatz [16] introduced for electronic structure
problems, adaptive methods designed to produce expressive ansatze at lower
complexity cost [17, 18], and different strategies to solve condensed matter
models with strong quantum correlation [12, 19]. The great majority of these
works has been essentially devoted to the case of many-fermion systems. In
contrast, little attention has been paid to the capability of current devices
to produce interesting ansatze for the study of strongly interacting many-
boson systems. While state-of-the-art methods such as quantum Monte-Carlo
simulations have proved effective on classical computers [20, 21, 22, 23, 24,
25, 26, 27, 28, 29, 30], there is still a growing interest in conducting
quantum experiments to treat many-boson systems, for instance based on
superconducting circuits [31, 32, 33, 34] or photonic quantum devices [35, 36,
37, 38, 39, 40, 41], essentially used as analog quantum simulators.

Motivated by the development of VQAs, we investigate the case of strongly
correlated many-boson systems using a digital photonic-based Continuous
Variable (CV) quantum computer [42, 43, 44]. In practice, CV devices are
digital quantum computers that can manipulate information by entangling and
measuring quantum states of light using various types of quantum gates. The
bosonic nature of such photonic devices leads to a one-to-one mapping between
the photonic modes of a CV computer and the bosonic modes of other many-boson
systems. For variational algorithms, the CV setting is then particularly well-
suited to explore new types of quantum circuit ansatz to encode strongly
correlated many-boson wavefunctions. The rich complexity of such type of
wavefunctions can be qualitatively reproduced by studying the attractive
Bose–Hubbard (BH) model. The ground state of this model has attracted a lot of
attention recently due to its rich transitional behavior that evolves between
the so-called Schrödinger cat-like states to superfluid regimes [45, 46, 47,
48, 49, 50, 51]. Furthermore, the BH model has been successfully employed to
describe the physics of ultracold atoms and molecules in optical lattices [52,
53, 54, 55, 56, 57], quantum cooling protocols [33] and the evolution of
vibrational bosons in $\alpha$-helices proteins [58, 59].

In this work, we propose two different CV quantum circuit ansatze dedicated to
the encoding of many-boson wavefunctions. By performing VQAs based on the
state-infidelity as a cost function, we show that the many-boson ground state
of the BH model can be successfully realized on the CV quantum circuit with a
relatively small number of gates and parameters. As the state-infidelity
cannot be computed for large systems (as it requires the knowledge of the
exact wavefunction), a realistic experiment on a true quantum devices requires
the ability to measure the ground-state energy instead. Thus, we propose a
photon-counting protocol to measure the expectation value of the BH
Hamiltonian in order to simulate a variational quantum eigensolver (VQE)
algorithm, as it would be done in a true experiment, with and without sampling
noise.

The paper is organized as follows. After a brief introduction to the
attractive BH Hamiltonian and the exact diagonalisation method in sections
II.1 and II.2, we illustrate the exotic structural transition of the BH ground
state on three small-sized networks in section II.3. After introducing the
context of VQAs in section III.1, we present the architectures of our proposed
ansatze in section III.2. Numerical investigations of their encoding ability
are performed in section IV. We extend our study to ideal and realistic
simulations of a VQE algorithm in section V. Finally, conclusions and
perspectives are provided in section VI.

##  II Bose-Hubbard model

In this first section, we present a general introduction of the attractive BH
model. Our objective is to give an overview of the properties of the model and
to illustrate the rich behaviors that arise from the increase of the many-body
interaction. We will focus on three small-sized BH models, namely the BH dimer
and the three- and four-site periodic chains. These systems will be used later
in this paper as test-bed for the quantum ansatze designed herein to encode
many-boson wavefunctions on photonic devices.

###  II.1 Hamiltonian and regimes of interactions

The attractive BH model describes the physics of $N_{B}$ bosons hopping
between the $N_{S}$ sites of a network with local attractive many-body
interactions. The associated Hamiltonian operator reads

$\hat{\mathcal{H}}=-J\sum_{\langle
p,q\rangle}^{N_{S}}(b_{p}^{\dagger}b_{q}+b_{q}^{\dagger}b_{p})-\frac{U}{2}\sum_{p}^{N_{S}}n_{p}(n_{p}-1),$
(1)

where $J>0$ is the hopping parameter describing the delocalization of the
bosons throughout the connected nodes of the BH network, and $U>0$ is the
many-body interaction term corresponding to the local boson-boson interaction
that occurs when at least two particles occupy the same site. In this work,
$J$ is the energy unit (i.e. we set $J\equiv 1$ in every simulation).

To characterize the competition between the many-body interactions and the
hopping terms in the system, a dimensionless parameter

$\Lambda=\frac{N_{B}U}{J}$ (2)

is commonly introduced. Typically, $\Lambda\ll 1$ represents a regime of weak
interaction (the so-called “Super-fluid regime”) in which the physics of the
system is governed by the hopping terms making the bosons able to efficiently
delocalize over the network. The case $\Lambda\sim 1$ corresponds to an
intermediate quantum regime where hoppings and local many-particle
interactions compete together and where phase transitions usually occur. The
last case, $\Lambda\gg 1$, corresponds to the regime where the many-body
interaction dominates (the so-called “Fock regime”). In the case of attractive
interactions, the resulting state tends to localize the bosons on local sites
of the system (giving rise to macroscopic-cat states as illustrated later in
this work).

###  II.2 Ground state reference: exact diagonalization

To investigate the exotic phases encoded in the attractive BH model, we
compute the ground state of the Hamiltonian $\hat{\mathcal{H}}$ exactly by
solving the time-independent Schrödinger equation

$\hat{\mathcal{H}}\ket{\Psi_{0}}=E_{0}\ket{\Psi_{0}},$ (3)

where $\ket{\Psi_{0}}$ is the ground state of the system with the associated
energy $E_{0}$. We solve Eq. (3) by exact diagonalization (ED) (see Refs. 60,
61 for an introduction of the technique) to obtain a good reference state to
compare with our quantum ansatz developed for CV devices. The way ED works is
simple: given a system of $N_{B}$ bosons evolving on a network of $N_{S}$
sites, one builds an exact matrix representation of the associated Hamiltonian
and diagonalize it. In the case of the attractive BH model, the dimension
$D_{\hat{\mathcal{H}}}$ of the resulting matrix is equal to the dimension of
the Fock space of the system,

$D_{\hat{\mathcal{H}}}=\frac{(N_{B}+N_{S}-1)!}{N_{B}!(N_{S}-1)!},$ (4)

which represents the number of ways of distributing $N_{B}$ bosons on $N_{S}$
sites. See Tab. 1, for an overview of the scaling of the
$D_{\hat{\mathcal{H}}}$ as a function of $N_{S}$ and $N_{B}$.

$\bm{N_{S}}$ $\bm{N_{B}}$ | 2 | 3 | 4 | 5 | 8 | 16
---|---|---|---|---|---|---
2 | 3 | 4 | 5 | 6 | 9 | 17
3 | 6 | 10 | 15 | 21 | 45 | 153
4 | 10 | 20 | 35 | 56 | 165 | 969
8 | 36 | 120 | 330 | 792 | 6435 | 245157
Table 1: Evolution of the dimension $D_{\hat{\mathcal{H}}}$ of the BH model as
a function of the number of sites $N_{S}$ and bosons $N_{B}$.

While exact diagonalization provides access to the exact solution of the
problem, the exponential scaling of the Fock space dimension with respect to
the number of sites and bosons limits the use of ED to very small systems in
practice.

###  II.3 Illustration of the different regimes of interaction: exact
diagonalization on small systems

As the many-body interaction increases, the ground state of the attractive BH
model exhibits rich transitional behaviors. Such transitions have already been
extensively discussed in several recent works (see for example [49, 48, 51]).
Here, we briefly illustrate some of these exotic phases in three systems (that
will be treated later using quantum ansatze developed on a CV device), namely
the BH dimer, the three-site and the four-site periodic chain.

We employ the Inverse Participation Ratio (IPR) [62, 63, 64, 65] and the von
Neumann entropy [66, 67, 68] as two indicators to follow the expansion of
$\ket{\Psi_{0}}$ in the configuration basis when $\Lambda$ increases. Applied
to the Fock state basis, the IPR estimates the number of different bosonic
configurations forming $\ket{\Psi_{0}}$ by linear combination, and reads

$\text{IPR}(\ket{\Psi_{0}})=\Big{(}\sum_{n_{1},\ldots,n_{N_{S}}}|\braket{n_{1},\ldots,n_{N_{S}}}{\Psi_{0}}|^{4}\Big{)}^{-1},$
(5)

where $\ket{n_{1},\ldots,n_{N_{S}}}$ is a Fock state with $n_{q}$ bosons on
site $q$ (and $\sum_{q}n_{q}=N_{B}$). Within this definition, a state
$\ket{\Psi_{0}}$ localized on a single state of the Fock-state basis has an
IPR of 1. By contrast, the IPR of a state uniformly delocalized over all the
Fock states is equal to $D_{\hat{\mathcal{H}}}$.

The von Neumann entropy is a measure used to quantify the entanglement (non-
local quantum correlation) existing between different sub-parts of a quantum
system. In many-body systems, it quantifies the degree of entanglement and
correlation arising between the subsystem $p$ and all the other modes (also
commonly called “bi-partite entanglement”). It reads

$S=-\sum_{k}\lambda_{k}\log(\lambda_{k}),$ (6)

where $\lambda_{k}$ are the eigenvalues of the reduced density matrix of mode
$p$,

$\hat{\sigma}_{p}=\sum_{n_{k}=0}^{N_{B}}\lambda_{k}\ket{n_{k}}\bra{n_{k}},$
(7)

defined by taking the partial trace of the full-system density operator
$\hat{\rho}$ over the Fock-state basis of all the other modes in the systems,
$\hat{\sigma}_{p}=\mathbf{Tr}_{q\neq p}[\hat{\rho}]$ (note that in our case,
the full density operator is $\hat{\rho}=\ket{\Psi_{0}}\bra{\Psi_{0}}$ as one
focuses on the ground state). For the sake of simplicity, we will only focus
on a single-site von-Neumann entropy as defined in Eq. (6).

Figure 1: Evolution of the ground state structure of the dimer as a function
of $\Lambda$ (obtained with ED). Upper panel: IPR values for different values
of total boson number $N_{B}=2,\ldots,8$. Middle panel: Single site entropy
value $S$ values for different values of total boson number
$N_{B}=2,\ldots,8$. Lower panel: Occupancy probability of the ground state in
the Fock state basis for $N_{B}=8$. Probabilities are illustrated with shades
of grey ranging from white for zero probability, to black for $50\%$.

Figure 2: IPR and the single site von-Neumann entropy $S$ of the attractive BH
ground state as a function of $N_{B}$ and $\Lambda$ (obtained by ED). Left
panel: Three-site BH periodic chain. Right panel: Four-site BH periodic chain.
In both panels, the vertical dashed lines mark values of $\Lambda$
representative of three different many-body regimes. These values will be used
later on to illustrate the expressibility of our quantum ansatze on a CV
circuit.

In Figure 1, we show the evolution of the ground state properties of the BH
dimer as a function of the interaction parameter $\Lambda$ (for different
fixed number of bosons $N_{B}$). When $\Lambda\sim 0$, the system is in a
superfluid phase: every bosons propagate independently on the dimer and do no
interact together. The independent-boson picture leads to a trivial product-
state solution for the ground state as
$\ket{\Psi_{0}}=\frac{1}{\sqrt{2}\sqrt{N_{B}!}}(a_{0}^{\dagger}+a_{1}^{\dagger})^{N_{B}}\ket{0,0}$
(with $\ket{0,0}$ being the vacuum state of the dimer). Then, increasing
$\Lambda$ leads to a progressive restructuring of the ground state as readily
seen by the increasing IPR. This shows a spreading of the state over the whole
Fock basis that is also accompanied by an increasing single-site von-Neumann
entropy. Both quantities then reach a maximum value for $\Lambda\approx 3$.
The large variations of these two measures are not without significance and
reveal the appearance of quantum correlations and entanglement in the system.
Indeed, in the region $2\leq\Lambda\leq 4$ the shape of the ground state gets
more complex as it involves a larger number of bosonic configurations
(actually almost all the Fock states are involved), and locally the bosonic
modes get strongly entangled. Finally, as $\Lambda$ increases the structure of
the ground state loses in complexity and progressively tends to the well-known
“cat-like” form [49, 48] representing a superposition of two macroscopically
occupied states (i.e.
$\ket{\Psi_{0}}=\frac{1}{\sqrt{2}}(\ket{N_{B},0}+\ket{0,N_{B}})$). The IPR
reaches a value of 2, thus indicating the importance of these two fully
localized bosonic configurations. The entanglement entropy $S$ of the modes
decreases as less bosonic configurations get correlated. To better picture the
full restructuring process occurring during this exotic quantum transition, we
represent the decomposition of the BH dimer’s ground state in the Fock basis
as a function of $\Lambda$ in the lower panel in Fig. 1.

Turning to larger systems, similar exotic quantum transitions can be observed.
In Fig. 2 we see that the IPR and the single-site von-Neumann entropy also
reach their maximal value around $\Lambda\approx 5$ for the three- and the
four-site periodic chains. In this regime, the complexity of the ground state
significantly grows with the total number of sites $N_{S}$ but also with the
total number of bosons $N_{B}$. As an illustration, there are around
$\text{IPR}^{max}\approx 42$ bosonic configurations involved in the structure
of the ground state at the maximum value of the IPR ($\Lambda\approx 4.8$) for
the three-site BH model with $N_{B}=8$. This number gets even higher for the
four-site model with $N_{B}=5$ and $\Lambda\approx 4$ with a total number of
$\text{IPR}^{max}\approx 55$ bosonic configurations. When $\Lambda$ goes
beyond this region of high complexity, we retrieve the so-called cat-states as
indicated by the number of macroscopically-occupied states ($\text{IPR}=3$ and
$4$ for the three and four-site systems, respectively, e.g. for the three-site
model we have
$\ket{\Psi_{0}}\propto\ket{N_{B},0,0}+\ket{0,N_{B},0}+\ket{0,0,N_{B}}$).

##  III Photonic-gate-based ansatze

In the following, we introduce the different theoretical ingredients to encode
the BH ground state on a photonic-gate-based device.

###  III.1 Variational Quantum Algorithms

VQAs represent state-of-the-art methods for NISQ devices that have been
proficiently used to determine the ground state properties of various many-
body systems (see for example [16, 17, 18, 12, 19] and the recent review about
VQAs [69]). The protocol of a VQA usually follows the steps: $i$) map the
problem onto the quantum device, $ii$) use a quantum circuit as an “ansatz” to
generate a trial state that encodes the wavefunction of the many-body system,
$iii$) evaluate a given cost function by repeated measurements, $iv$) optimize
the quantum gates parameters of the ansatz circuit to minimize the cost
function governed by a variational principle. In practice, one applies a
unitary transformation $U(\vec{\theta})$ with tunable gate parameters
(formally represented by the vector $\vec{\theta}$) to transform an initial
state $\ket{\Psi_{ini}}$ into a trial state

$\ket{\Psi(\vec{\theta})}=\hat{U}(\vec{\theta})\ket{\Psi_{ini}}.$ (8)

Often $\ket{\Psi_{ini}}$ is the “all-zero” state, but can be something else
such as the Hartree–Fock state in electronic structure problems. If one seeks
the ground state of a system, the cost function is usually the ground-state
energy defined as

$E(\vec{\theta})=\bra{\Psi(\vec{\theta})}\hat{\mathcal{H}}\ket{\Psi(\vec{\theta})}\geq
E_{0},$ (9)

where $E_{0}$ is the exact ground-state energy and is a lower bound to the
ground-state energy of the trial state $\ket{\Psi(\vec{\theta})}$. In a real
VQA experiment, the expectation value $E(\vec{\theta})$ is estimated by
statistical sampling. The energy estimator will have an uncertainty based on
the number of state preparations and measurements used. This energy is then
passed to an external classical optimization routine which provides a new set
of parameters $\vec{\theta}$ used to update the quantum circuit. The whole
optimization process is repeated until global convergence is reached (i.e.
until $E(\vec{\theta})$ is minimized).

The efficiency of the VQA lies in the capability of the unitary
$\hat{U}(\vec{\theta})$ to encode the ground state of the system, i.e. it has
to show a strong “expressibility” so that the transformation given in Eq. (8)
can efficiently map the initial state $\ket{\Psi_{ini}}$ to the exact ground
state $\ket{\Psi_{0}}$. The question of the expressibility of quantum ansatze
is hard to tackle and remains an open question in the quantum computing
community. Only a few works have tried to investigate this problem on qubit
devices with a special emphasis on many-electron problems (see for example the
study of the unitary couple-cluster-ansatz in [70]), and one relies more on
numerical evidence. A lot of attention has been directed to many-fermion
systems (see for example [17, 18, 12, 19, 71, 72, 73]), but little work has
been done on quantum ansatze to encode many-boson wavefunctions (for a notable
exception, see the implementation of the unitary vibrational-coupled-cluster
ansatz introduced in Ref. 15).

###  III.2 Photonic quantum circuits as ansatze for many-boson wavefunction

Figure 3: Illustration of a VQA approach where each photonic mode in the
quantum circuit encodes the bosonic occupation of a particular site in the
Bose-Hubbard network.

In contrast to most of the studies focusing on qubit devices, we tackle the
problem using CV circuits. The choice of using CV devices to describe many-
boson problems is rather intuitive as there exists a straightforward one-to-
one correspondence between the modes of a bosonic system and the photonic
modes of the CV device. In the case of the attractive BH system, each photonic
mode $p$ of the quantum circuit can be used to encode the bosonic population
of the associated site $p$, as illustrated in Fig. 3. Starting from this
advantageous property, we will build VQAs based on CV quantum circuits to
encode a many-boson wavefunction trial state as defined in Eq. (8).

####  III.2.1 Utilized CV quantum gates

Concerning the basic properties of the ansatz, we know that the BH Hamiltonian
in Eq. (1) is boson-number-preserving by construction (i.e. $N_{B}$ is a
constant of the problem). Thus, we restrict ourselves to three types of
photon-number-preserving gates, namely beam-splitters, rotation gates and Kerr
gates. The beam-splitter gate is a two-mode gate that affects two inputs modes
simultaneously (labeled here $p$ and $q$) and whose unitary is defined as

$\displaystyle
B_{p,q}(\theta,\phi)=\exp\left(\theta(e^{i\phi}a_{q}^{\dagger}a_{p}-e^{-i\phi}a_{p}^{\dagger}a_{q})\right),$
(10)

where $\theta$ and $\phi$ are respectively the transmittivity and phase angle
of the gate. The action of the beam-splitter’s unitary on the associated
ladder operators ($a_{p}/a_{p}^{\dagger}$ and $a_{q}/a_{q}^{\dagger}$) is
given by

$\begin{split}B_{p,q}(\theta,\phi)^{\dagger}a_{p}B_{p,q}(\theta,\phi)&=a_{p}\cos(\theta)-a_{q}\sin(\theta)e^{-i\phi},\\\
B_{p,q}(\theta,\phi)^{\dagger}a_{q}B_{p,q}(\theta,\phi)&=a_{q}\cos(\theta)+a_{p}\sin(\theta)e^{i\phi}.\end{split}$
(11)

To have a better illustration of the effects of the beam-splitter gate on a
many-photon wavefunction, figure 4 shows the evolution of the density of
probability of a two mode-state $B(\theta,0)\ket{\Psi_{ini}}$ (with
$\ket{\Psi_{ini}}=\ket{8,0}$ or $\ket{4,4}$) in the Fock basis as a function
of the transmittivity parameter $\theta$. As shown in this figure, the main
effect of the beam-splitter is to delocalize the photons from a mode to
another one creating various patterns in the Fock basis. The latter strongly
depend on the initial state considered.

Turning to the rotation gate, it affects the associated Fock state as
$R_{p}(\theta)\ket{n_{p}}=e^{i\theta n_{p}}\ket{n_{p}}$ with the unitary
operator defined as

$\displaystyle R_{p}(\theta)=\exp\left(i\theta a_{p}^{\dagger}a_{p}\right).$
(12)

Finally, the Kerr gate is a non-Gaussian single-mode gate whose unitary reads

$\displaystyle
K_{p}(\theta)=\exp\left(i\theta(a_{p}^{\dagger}a_{p})^{2}\right).$ (13)

For a given modes $p$, this leads to $K_{p}(\theta)\ket{n_{p}}=e^{i\theta
n_{p}^{2}}\ket{n_{p}}$, thus generating a non-linear phase term that depends
on the squared number of bosons present in the mode $p$,
$n_{p}^{2}=(a_{p}^{\dagger}a_{p})^{2}$. In the current devices, Kerr gates
represent expensive resources which can only be implemented non-
deterministically. Nonetheless, our analysis shows that they play a key role
in the performance of our ansatze.

Figure 4: Density of probability of a state $B(\theta,0)\ket{\Psi_{ini}}$ in
the Fock basis as a function of the transmittivity parameter $\theta$ of the
beam-splitter unitary. The probability is illustrated with shades of grey
ranging from white for zero to black for unity. Upper panel:
$\ket{\Psi_{ini}}=\ket{4,4}$. Lower panel: $\ket{\Psi_{ini}}=\ket{8,0}$

####  III.2.2 Minimal beam-splitter-Kerr ansatz

Let us now introduce our first ansatz to encode the ground state of the
attractive BH model, that we call the minimal beam-splitter-Kerr (BS-Kerr)
ansatz. This ansatz was designed to efficiently work with as few quantum gates
and parameters to optimize as possible, based on numerical investigations. The
most efficient architecture we found (in terms of the total number of gates
required to encode the ground state with sufficient accuracy) resulted in
layers of beam-splitter and Kerr gates as shown in Fig. 5 for two consecutive
layers. Note that the phase angle parameters of the beam-splitters are always
set to $\phi=0$ in this ansatz. Indeed, we did not see any improvement when
considering the latter in our numerical simulations (not shown). Discarding
allows to reduce the complexity of the classical optimization with no impact
on the final accuracy.

Figure 5: Structure of the BS-Kerr ansatz for $N_{S}=2$, 3 and 4 sites with
two consecutive layers. Each layer is composed of a series of $N_{S}-1$ beam-
splitter gates in stairs (up or down) followed by a series of $N_{S}$ Kerr
gates applied on each modes. The total gate count is $N_{L}(2N_{S}-1)$ where
$N_{L}$ is the number of layers.

As readily seen in Fig. 5, each layer is defined as a series of $N_{S}-1$
beam-splitter gates in stairs (up or down) followed by a series of $N_{S}$
Kerr gates applied on each mode, thus leading to a gate count of
$N_{L}(2N_{S}-1)$ (and the same number of parameters to optimize) where
$N_{L}$ is the number of layers. From a layer to another one, we alternate the
stair structure to ensure that each mode can have the opportunity to exchange
photons (via the beam-splitters) with all the other ones in the circuit.

Note that removing the Kerr gates from the ansatz led to very poor results
(not shown), thus demonstrating the significant role of non-linearity to
encode the ground state of the attractive BH model. Although there are
materials that exhibit the $\chi^{(3)}$ non-linearity that is required for the
Kerr interaction, the experimental realization of a Kerr gate remains a
challenge at the present moment [74]. Despite such technological obstacle,
there is still place for improvements on the implement of non-linear quantum
gates [75].

####  III.2.3 Interferometer-Kerr ansatz

The BS-Kerr ansatz was designed to encode the ground state of the attractive
BH model with a minimal number of gates and parameters, which is optimal in
numerical simulation. However, in the experimental realization one wants to
have an optimal use of the investigation tools (such as interferometers) and,
if possible, to reduce the number of Kerr gates that are difficult to
implement. Therefore, we considered an (more experimentally-friendly)
alternative ansatz composed of layers of one interferometer and one Kerr gate
per mode, as depicted in Fig. 6 for 3 and 4 modes. The interferometer is
composed of $N_{S}(N_{S}-1)/2$ beam-splitters and $N_{S}$ rotation gates,
making a total of $N_{L}N_{S}(N_{S}+3)/2$ gates and $N_{L}N_{S}(N_{S}+1)$
parameters (as one beam-splitter accounts for two parameters).

At first sight, this ansatz appears less efficient than the BS-Kerr ansatz in
terms of number of gates and parameters per layer. However, we can expect the
additional expressivity brought about by the interferometers to lead to a
reduction in the number of layers required to encode the ground state. In
other words, the interferometer-Kerr ansatz could lead to a reduction in the
total number of Kerr gates at the expense of additional beam-splitter and
rotation gates. This hypothesis is numerically investigated in Sec. IV.2.

The complexity of the ansatz can be reduced by getting rid of the rotation
gates (physically irrelevant for most applications [76]) and/or setting the
phase angles of all beam-splitters to $\phi=0$ (and keeping only the
transmittivity angles) as already done for the BS-Kerr ansatz.

Figure 6: Single layer of the Interferometer-Kerr ansatz for $N_{S}=3$ (top)
and $N_{S}=4$ (bottom) modes.

##  IV Numerical investigation of the ansatze expressibility

In this section, we present our numerical results focusing on the encoding
capacity of the two ansatze introduced previously applied to the BH dimer, the
three-site and the four-site periodic chains. The simulations of the photonic
circuits were realized using the Python package Strawberryfields [77]. For the
optimization of the quantum circuits parameters, the L-BFGS-B method was used
and encoded in the OptimParrallel package [78]. The circuit’s parameters were
always initialized with an amplitude randomly drawn in $[-0.05,+0.05]$, and
the maximum number of iterations was set to 20000. A home-made ED code has
been implemented in Python to compute the exact ground state by exact
diagonalization.

###  IV.1 Capacity of the BS-Kerr ansatz

To demonstrate the capacity of the different ansatze to encode the ground
state of BH networks, we estimated the fidelity of the trial state with
respect to the exact ground state of the system $\ket{\Psi_{0}}$,

$\mathcal{F}=|\braket{\Psi_{0}}{\hat{U}(\vec{\theta})}{\Psi_{ini}}|^{2}.$ (14)

Assuming that the classical optimization process does not end up in local
minima, an ansatz with a strong encoding capability will return a fidelity
$\mathcal{F}$ close to unity. We numerically check this property by optimizing
the quantum circuit’s parameters to minimize a complementary measure that is
the infidelity

$\displaystyle\mathcal{I}=1-\mathcal{F}$ (15)

that is naturally lower-bounded by zero. In the spirit of Ref. 12 on the
Fermi–Hubbard model, we progressively increase the number of layers in the
quantum circuits and perform a VQA where the infidelity is used as a cost
function. As the number of layers increases, so does the expressibility of the
ansatz in practice. We continue to increase the number of layers until the
fidelity of the ground state reaches at least $\mathcal{F}\geq 99\%$
(conversely $\mathcal{I}<1\%$). Note that no mention has been made so far
about the nature of the initial state $\ket{\Psi_{ini}}$, for which we
considered two choices: a monomodal and a bimodal states, detailed later in
section IV.1.1.

####  IV.1.1 BH model with $N_{S}=2$ sites

Figure 7: BH dimer. Minimum depth required (in terms of number of layer
$N_{L}$ and gate count) to generate a fidelity $\mathcal{F}\geq 99\%$ [see Eq.
(14)] for different values of $\Lambda$ and $N_{B}$. Upper panel: results
obtained with the monomodal initial state. Lower panel: results obtained with
the bimodal intitial state. Figure 8: BH dimer. Probability of occupation of
the ground state in the Fock basis obtained by the BS-Kerr ansatz (green bars)
with $N_{L}=6$, $N_{B}=8$, $\ket{\Psi_{ini}}=\ket{4,4}$ and different values
of $\Lambda$. Results are compared with exact diagonalization (black bars).

Focusing first on the BS-Kerr ansatz for the BH dimer, the minimum number of
layers required to reach $\mathcal{F}\geq 99\%$ for $N_{B}=1$ to $N_{B}=8$ is
shown in Fig. 7 for $\Lambda$ values corresponding to different regimes of
correlation (as marked by the vertical red-dashed lines in Fig. 1). Two
different initial states were considered: a monomodal one given by
$\ket{\Psi_{ini}}=\ket{N_{B},0}$ and a bimodal one given by
$\ket{\Psi_{ini}}=\ket{\frac{N_{B}}{2},\frac{N_{B}}{2}}$ (n.b. for odd total
number of bosons we took
$\ket{\Psi_{ini}}=\ket{\frac{N_{B}+1}{2},\frac{N_{B}+1}{2}-1}$ instead). As a
first general remark, a fidelity higher than 99$\%$ is reached in all cases
(whatever the values of correlation strength $\Lambda$, the number of bosons
$N_{B}$ and the choice of initial state). This clearly demonstrates the
potential of CV devices to efficiently encode the many-boson ground state of
the attractive BH dimer even in the strongly correlated regime. Note that to
generate an efficient encoding of the wavefunction, the larger $N_{B}$ is the
larger the number of layers $N_{L}$ needs to be. This is expected as the
complexity of the ground state grows with the size of the Fock space
$D_{\hat{\mathcal{H}}}$, which is linear in $N_{B}$ for the BH dimer as shown
in Tab. 1). Interestingly, the trend of the curves ensuring an efficient
encoding is globally sub-linear in $N_{B}$. In other words, choosing a number
of layers such as $N_{L}\approx 0.7\times N_{B}$ seems to guarantee a very
high ground-state fidelity in all cases considered herein. Another interesting
fact is that in the weakly correlated case, $\Lambda=0.01$, a single layer of
the BS-Kerr ansatz is sufficient to reach a fidelity higher than 99$\%$, even
when the number of bosons increases. This is true only when the monomodal mode
is used as the initial state. We rationalize this result in Appendix A, where
we show that starting from the monomodal state, an exact encoding
($\mathcal{F}=100\%$) of the BH ground state can be realized using only beam-
splitters.

Furthermore, we computed the probability of occupation of our trial ground
state in the Fock basis, shown in Fig. 8. This illustrates the structural
transition of the BH dimer’s ground state as a function of $\Lambda$ for
$N_{B}=8$ bosons and $N_{L}=6$ layers. In this case, we clearly see the
structural transformation of the dimer ground state as $\Lambda$ increases,
starting from the superfluid regime ($\Lambda=0.01$) and going through the
maximum spreading structure ($\Lambda=3$) to finally reach the famous “cat-
like” superposition ($\Lambda=5$ to 10). For all values of $\Lambda$, the
probability of occupation of our trial state (green bars) is on top of the
exact one (black bars) thus showing how efficient our quantum ansatz is to
encode every exotic structures adopted by the ground state through its
transitions.

####  IV.1.2 BH model with $N_{S}=3$ and $4$ sites

Figure 9: 3-site BH model. Minimum depth required (in terms of number of layer
$N_{L}$ and gate count) to generate a fidelity $\mathcal{F}\geq 99\%$ [see Eq.
(14)] for different values of $\Lambda$ and $N_{B}$. Upper panel: results
obtained with a monomodal initial state. Lower panel: results obtained with a
bimodal intitial state. Figure 10: 4-site BH model. Minimum depth required (in
terms of number of layer $N_{L}$ and gate count) to generate a fidelity
$\mathcal{F}\geq 99\%$ [see Eq. (14)] for different values of $\Lambda$ and
$N_{B}$. Upper panel: results obtained with the monomodal initial state. Lower
panel: results obtained with the bimodal intitial state.

Turning to larger systems, Figs. 9 and 11 show the results obtained on the
three-site BH model for $\Lambda$ values corresponding to different regimes of
correlation (as marked by the vertical red-dashed lines in the left panel Fig.
2). A monomodal and bimodal initial states are again considered and read
$\ket{\Psi_{ini}}=\ket{N_{B},0,0}$ and
$\ket{\Psi_{ini}}=\ket{\frac{N_{B}}{2},0,\frac{N_{B}}{2}}$ respectively (n.b.
for odd total number of bosons we took
$\ket{\Psi_{ini}}=\ket{\frac{N_{B}+1}{2},0,\frac{N_{B}+1}{2}-1}$ instead).
Compared to the BH dimer, the dimension of the Fock space is non-linear in
$N_{B}$ (see Tab. 1) thus leading to a maximum number of bosonic
configurations accessible of $D_{\hat{\mathcal{H}}}=45$ when $N_{B}=8$
(approximately 6 times larger than the number of configurations on the dimer
for the same number of bosons). In spite of this, we observed that the BS-Kerr
ansatz performs very well as shown in Fig. 9. In this case, we retrieve some
differences between the results obtained with the two initial states. The
monomodal initial state seems to perform better and less layers are globally
required to encode the ground state with a fidelity higher than 99$\%$. As for
the BH dimer, only one layer is required for $\Lambda=0.01$ when this initial
state is used.

Still, the ansatz is clearly able to efficiently transform both initial states
into a very good approximation to the exact ground state. As readily seen in
Fig. 9, choosing a circuit depth with $N_{L}\sim 1.25\times N_{B}$ guaranties
an efficient encoding of the ground state in all cases considered here for the
attractive three-site BH model.

The complexity grows even stronger for the four-site BH network. We follow the
exact same strategy than for the previous models by considering
$\ket{\Psi_{ini}}=\ket{N_{B},0,0,0}$ as the monomodal initial state and
$\ket{\Psi_{ini}}=\ket{\frac{N_{B}}{2},0,\frac{N_{B}}{2},0}$ as the bimodal
one (n.b. for odd total number of bosons we took
$\ket{\frac{N_{B}+1}{2},0,\frac{N_{B}+1}{2}-1,0}$ instead), as well as
$\Lambda$ values corresponding to different regimes of correlation marked by
the vertical red-dashed lines in the right panel Fig. 2. For this system, the
dimension of the Fock space explodes very quickly and fixing $N_{B}=5$ already
leads to a number of bosonic configurations of $D_{\hat{\mathcal{H}}}=56$ that
is larger than the size of the Fock space of the three-site model with
$N_{B}=8$. In our simulations, we chose to limit ourselves to a maximal number
of bosons of $N_{B}=5$ due to the numerical complexity of the classical
simulation of our VQA. The performance of the BS-Kerr ansatz for $N_{B}=1$ to
5 is shown in Fig. 10. In the full range of correlation regimes, the photonic
circuit is again able to encode the ground state with high fidelity for a
circuit depth of $N_{L}\approx 2N_{B}$.

As a complement, Figs. 11 and 12 illustrate the structural transition of the
BH ground state for the three and four-site BH networks, respectively. Just as
in the BH dimer, the density of probability of our trial ground state matches
perfectly the one of the exact ground state. A remarkable property is that we
always consider a single bosonic configuration as an initial state
$\ket{\Psi_{ini}}$. Starting from this, the quantum circuit can efficiently
transform the latter into a very complex linear combination of
$D_{\hat{\mathcal{H}}}$ configurations encoding strong quantum correlations.

Figure 11: 3-site BH model. Probability of occupation of the ground state in
the Fock basis obtained by the BS-Kerr ansatz (green bars) with $N_{L}=6$,
$N_{B}=4$, $\ket{\Psi_{ini}}=\ket{2,0,2}$ and different values of $\Lambda$.
Results are compared with exact diagonalization (black bars). Figure 12:
4-site BH model Probability of occupation of the ground state in the Fock
basis obtained by the BS-Kerr ansatz (green bars) with $N_{L}=6$, $N_{B}=3$,
$\ket{\Psi_{ini}}=\ket{2,0,1,0}$ and different values of $\Lambda$. Results
are compared with exact diagonalization (black bars).

###  IV.2 Capacity of the interferometer-Kerr ansatz

Figure 13: 3-site BH model. Minimum depth required (in terms of number of
layer $N_{L}$ and gate count) to generate a fidelity $\mathcal{F}\geq 99\%$
[see Eq. (14)] for different values of $\Lambda$ and $N_{B}$, using the
bimodal initial state. First panel: Full ansatz, Second panel: Phase angles of
the beam-splitters are set to 0, Third panel: Phase angles of the beam-
splitter are reintroduced but rotation gates are removed, Fourth panel: Phase
angles of the beam-splitters are set to 0 and the rotation gates are removed.

Let us now turn to the interferometer-Kerr ansatz described in Sec. III.2.3.
By using an interferometer composed of more beam-splitter and rotation gates,
the interferometer-Kerr ansatz is expected to perform as efficiently as the
BS-Kerr ansatz but with a reduced number of layers, hence a reduced number of
Kerr gates. As readily seen in the top panel of Fig. 13, only $N_{L}=6$ layers
are necessary to encode more than 99$\%$ of the ground state of the three-site
model for $N_{B}=8$, compared to the 10 layers required by the BS-Kerr ansatz.
Note that reducing the number of layers in the interferometer-Kerr ansatz does
not necessarily reduce the total number of gates compared to the BS-Kerr
ansatz. Indeed, $N_{L}=6$ layers in the interferometer-Kerr ansatz is
equivalent to 54 gates, while $N_{L}=10$ for the BS-Kerr ansatz gives 50
gates. However, the number of Kerr gates is always directly proportional to
the number of layers and the number of modes as $N_{\rm Kerr}=N_{L}\times
N_{S}$ such that the number of Kerr gates is almost reduced by half when using
the interferometer-Kerr ansatz (designed for this exact purpose). Note that
the phase angles of the beam-splitter gates were not set to 0 compared to the
BS-Kerr ansatz. For a more fair comparison, we set them to 0 in the second
panel of Fig. 13. One can notice a slight increase in the number of layers,
especially when the correlation strength $\Lambda$ increases, but only up to
$N_{L}=7$. Hence, for the three-site model, the interferometer-Kerr ansatz
guaranties an efficient encoding of the ground state with $N_{L}\sim 0.8N_{B}$
compared to $N_{L}\sim 1.25N_{B}$ for the BS-Kerr ansatz, for all values of
$\Lambda$.

We just saw that the total number of gates of the interferometer-Kerr ansatz
can still be higher than of the BS-Kerr ansatz, despite a significant
reduction in the number of Kerr gates. As mentioned in Ref. 76, the rotation
gates of the interferometer are physically irrelevant for most applications.
These gates are removed from the ansatz in the third and fourth panels of Fig.
13, where the phase angles of the beam-splitter are either free or set to 0,
respectively. Looking at the fourth panel, it is clear that simultaneously
removing the phase angles of the beam-splitter gates and the rotation gates
deteriorates the performance of the ansatz, and does not lead to any
improvement over the BS-Kerr ansatz. Even for $\Lambda=0.01$, this is the only
case where more than $N_{L}=3$ layers are required to accurately encode the
ground state. However, removing the rotation gates while keeping the phase
angle of the Beam-Slitter gates (third panel) seems equivalent to doing the
exact opposite, i.e. keeping the rotation gates and removing the phase angle
of the beam-splitter gates (second panel). By doing so, the number of
parameters remains the same but the total number of gates can be reduced, for
instance from 63 to 42 for $N_{B}=8$.

Figure 14: 4-site BH model. Minimum depth required (in terms of number of
layer $N_{L}$ and gate count) to generate a fidelity $\mathcal{F}\geq 99\%$
[see Eq. (14)] for different values of $\Lambda$ and $N_{B}$, using the
bimodal initial state. First panel: Full ansatz, Second panel: Phase angles of
the beam-splitters are set to 0, Third panel: Phase angles of the beam-
splitter are reintroduced but rotation gates are removed, Fourth panel: Phase
angles of the beam-splitters are set to 0 and the rotation gates are removed.

Turning to the four-site model in Fig. 14, we reach the exact same general
conclusions as for the three-site model. It seems that the circuit depth
required to encode the ground state of this model tends to $N_{L}\sim
1.2N_{B}$, which again improves over $N_{L}\sim 2N_{B}$ for the BS-Kerr
ansatz. Without the rotations and the phase angles (fourth panel), the number
of layers tends to increase in all regimes of correlation, while keeping the
phase angles and removing the rotation gates seems to provide the best trade-
off between the total number of gates, the number of Kerr gates (proportional
to the number of layers), and the number of parameters.

##  V VQE simulations using the BS-Kerr ansatz

In the previous section, we gave numerical evidence on the capacity of our
proposed ansatze to encode the ground state of the attractive BH model. In
this section, we switch to more practical calculations by simulating VQE
protocols where the cost function to minimize is the energy of the system,

$E(\vec{\theta})=\bra{\Psi_{ini}}U(\vec{\theta})^{\dagger}HU(\vec{\theta})\ket{\Psi_{ini}}\geq
E_{0}.$ (16)

For the sake of conciseness, we focus on the BS-Kerr ansatz only. First, we
briefly demonstrate the accuracy of the results obtained on the three and
four-site BH models. Second, we use the BH dimer as a test-bed to tackle the
practical simulation of a VQE experiment on a CV device. This second part will
be the occasion to compare the results obtained with an ideal realization of a
VQE to a more realistic one, where the estimation of the energy
$E(\vec{\theta})$ is based on a photon-counting protocol.

###  V.1 Ideal VQE simulations for the three- and four-site BH model

Figure 15: Probability of occupation in the Fock basis, obtained by ED (black
bars) and VQE simulations with the BS-Kerr ansatz (green bars). All parameters
were initialized with random amplitudes drawn in $[-0.05,+0.05]$, and the BFGS
optimizer with a maximum number of 2000 iterations were used. Upper panel:
Ground state of the three-site model with $N_{B}=4$, $N_{L}=6$ and
$\ket{\Psi_{ini}}=\ket{2,0,2}$ Lower panel: Ground state of the four-site
model with $N_{B}=3$, $N_{L}=6$ and $\ket{\Psi_{ini}}=\ket{2,0,1,0}$.

In this section, ideal simulations of the VQE algorithm are realized (i.e.
infinite sampling, no photon loss, no error in the circuit…). This means that
we assume that we can extract the trial state
$\ket{\Psi(\vec{\theta}}=U(\vec{\theta})\ket{\Psi_{ini}}$ at the end of the
circuit, which makes it possible to exactly estimate the energy
$E(\vec{\theta})$ during the optimization process.

Figure 15 shows the probability of occupation in the Fock basis of the ground
states obtained at the end of the VQE algorithm for different $\Lambda$ values
(consistent with the simulations realized in the previous section). Note that
$N_{L}=6$ layers and initial bimodal states ($\ket{\Psi_{ini}}=\ket{2,0,2}$
and $\ket{\Psi_{ini}}=\ket{2,0,1,0}$) are considered. The depth of the
circuits imposes to optimize $30$ (respectively $40$) parameters for the
three-site (four-site) network.

As shown in Fig. 15, the probability of occupation follows the same trend as
the previous ones in Figs. 11 and 12, where the infidelity $\mathcal{I}$ was
minimized instead of the ground-state energy. In all cases of Fig. 15, the
resulting states have fidelity $\mathcal{F}\geq 99\%$ and an energy error
$\Delta E=E(\vec{\theta}^{*})-E_{0}\leq 10^{-5}J$ (where $\vec{\theta}^{*}$
denotes the parameters that minimize the cost function), thus showing an
excellent convergence of the VQE algorithm. Note that similar ideal VQE
simulations have been also performed using monomodal initial states (i.e.
$\ket{\Psi_{ini}=\ket{4,0,0}}$ and $\ket{\Psi_{ini}}=\ket{3,0,0,0}$, not
shown), and led to similar accuracy, demonstrating the ability of the BS-Kerr
ansatz to efficiently map different initial states to the exact ground state
of the attractive BH model (in all correlation regimes considered herein).

###  V.2 Ideal VQE v.s. realistic VQE with sampling noise: application to the
BH dimer.

Let us now turn to the implementation of a more realistic execution of the VQE
algorithm, as it could be done in a practical experiments.

####  V.2.1 Measure of the expectation value of the Hamiltonian by a photon
counting approach

In practice, a realistic simulation of a VQE algorithm requires to measure the
energy of a given trial state at the end of the quantum circuit. For this
purpose, we propose a measurement protocol based on photon counting that can
be realized using photon number resolving (PNR) detectors [79, 80], which are
currently the most promising tools in quantum optics to efficiently measure a
well-defined number of photons simultaneously. In the following, the notation
$\langle\hat{O}\rangle_{\Psi}\equiv\bra{\Psi}\hat{O}\ket{\Psi}$ is used and
represents the quantum average of an operator $\hat{O}$ realized at the end of
the circuit when the photons are in the final state $\ket{\Psi}$.

Starting from the BH Hamiltonian $\hat{\mathcal{H}}$ defined in Eq. (1),
measuring the energy of the final photonic state $\ket{\Psi}$ of the circuit
requires measuring two different contributions,

$\langle\hat{\mathcal{H}}_{\text{Hoppings}}\rangle=-J\sum_{\langle
p,q\rangle}^{N_{S}}\langle(b_{p}^{\dagger}b_{q}+b_{q}^{\dagger}b_{p})\rangle_{\Psi},$
(17)

and

$\langle\hat{\mathcal{H}}_{\text{Loc-
Int}}\rangle=-\frac{U}{2}\sum_{p}^{N_{S}}\langle
n_{p}(n_{p}-1)\rangle_{\Psi}.$ (18)

The one-body (hopping) term in Eq. (17) can be measured by introducing
additional 50/50 beam-splitter gates between each pair of modes before
measuring, thus leading to

$B_{p,q}(\pi/4,0)\left[b_{p}^{\dagger}b_{q}+b_{q}^{\dagger}b_{p}\right]B_{p,q}(\pi/4,0)^{\dagger}=n_{q}-n_{p}$
(19)

such that

$\langle b_{p}^{\dagger}b_{q}+b_{q}^{\dagger}b_{p}\rangle_{\Psi}=\langle
n_{q}-n_{p}\rangle_{\tilde{\Psi}}$ (20)

where $\ket{\tilde{\Psi}}=B_{p,q}(\pi/4,0)\ket{\Psi}$. Thus, evaluating the
hopping term between two modes with respect to the original trial state
$\ket{\Psi}$ is equivalent to measuring the difference of average photon
numbers in the two same modes after applying 50/50 beam-splitter.

Turning to the many-body interaction term in Eq. (18), note that it is
directly related to the photon statistics of each individual mode at the end
of the circuit:

$\langle n_{p}(n_{p}-1)\rangle_{\Psi}=\text{Var}(n_{p})_{\Psi}+\langle
n_{p}\rangle_{\Psi}^{2}-\langle n_{p}\rangle_{\Psi},$ (21)

where $\text{Var}(n_{p})=\langle n_{p}^{2}\rangle_{\Psi}-\langle
n_{p}\rangle^{2}_{\Psi}$ is the variance of the photon number in mode $p$ at
the end of the circuit. In other words, to estimate the contribution of the
many-body interaction term to the total energy, we only need to evaluate (over
many samples) the first two moments of the number operator for each mode $p$.

Although not studied in this work, note that extended versions of the BH model
are sometimes considered (see for example Refs.[81, 57, 82, 83]), such as

$\hat{\mathcal{H}}^{\prime}=\hat{\mathcal{H}}+\sum_{p}\mu_{p}n_{p}+\sum_{p,q}V_{pq}n_{p}n_{q},$
(22)

where $\mu_{p}$ is usually referred to as the “chemical potential”, while the
last term encodes the dipole-dipole interaction of two packets of bosons
localized on different sites $p$ and $q$. The expectation value of these two
additional terms can also be estimated via photon counting. First, the
expectation value of the one-body term containing the chemical potential

$\langle\hat{\mathcal{H}}_{\text{Chem-
Pot}}\rangle_{\Psi}=\sum_{p}\mu_{p}\langle n_{p}\rangle_{\Psi}$ (23)

is relatively simple to estimate as it only requires to do local photon-
counting on each mode. The energy contribution of the dipole-dipole term

$\langle\hat{\mathcal{H}}_{\text{dip-
dip}}\rangle_{\Psi}=\sum_{p,q}V_{pq}\langle n_{p}n_{q}\rangle$ (24)

can be found by evaluating the covariance matrix of the photon number
distribution in all the modes:

$\displaystyle\langle n_{p}n_{q}\rangle=\mathrm{Cov}(n_{p},n_{q})+\langle
n_{p}\rangle\langle n_{q}\rangle.$ (25)

Note that setting $p=q$ in Eq. (25) gives back Eq. (21).

####  V.2.2 VQE with and without sampling noise

Figure 16: Ground-state fidelity $\mathcal{F}$ of the BH dimer obtained via
VQE simulations using the BS-Kerr ansatz, as a function of the number of
bosons $N_{B}$ and the number of layers $N_{L}$, for different values of
$\Lambda$. Left panel: noiseless VQE simulation using the BFGS optimizer with
a maximum number of 2500 iterations. All parameters are initialized with a
random amplitude in $[-1,+1]$. Right panel: VQE simulation with sampling noise
using the noise-robust optimizer CMA-ES [84] with an initial step-size
parameter of $\sigma_{0}=0.05$ and a maximum number of 20000 functions
evaluations (each function evaluation stems from sampling of $10^{8}$ samples
of the quantum circuit). All parameters are initialized with a random
amplitude in $[-0.1,+0.1]$. On each panel, upper and lower rows show results
obtained by starting from a monomodal and a bimodal initial state,
respectively.

Focusing on the BH dimer, we now investigate how the sampling of the energy
$E(\vec{\theta})$ can affect the quality of the VQE results. Compared to the
ideal VQE simulations realized above, considering sampling noise leads to
fluctuations in the energy $E(\vec{\theta})$ that makes the optimization
harder in practice.

Fig. 16 shows the ground-state fidelity $\mathcal{F}$ obtained from VQE
simulations with (right panel in blue) and without (left panel in green)
sampling noise. In the latter case (with sampling noise), the energy is
calculated using our photon-counting protocol described in Sec. V.2.1. To
mitigate the impact of noise in the classical optimization, we employed a
“noise-proof” optimizer, namely CMA-ES [84]. The number of samples is fixed to
$10^{8}$ and photon-counting is only realized on the first mode of the
circuit, as the conservation of the total number of bosons imposes
$N_{B}=n_{1}+n_{2}$ (i.e. measuring the number of photons $n_{1}$ in the first
mode gives directly the number of photons $n_{2}$ in the second mode, assuming
no photon loss).

Let us first focus on the ideal VQE simulation (left panel of Fig. 16). In
analogy with the previous study in Sec. IV, the ground state is encoded with
high fidelity whatever the initial state and correlation regime considered. In
all of these cases, a sub-linear limit marks the separation between the
regions of high and low fidelity-encoding (in average $N_{L}\approx
0.6N_{B}$). The slope of this limit actually grows with the strength of
$\Lambda$, thus indicating that “cat-like” states that arise for strong many-
body interactions can be the more complex case to encode with a VQE using the
BS-Kerr ansatz.

Turning to the noisy VQE simulation (right panel of Fig. 16) based on our
photon-counting protocol, one can see that the resulting ground-state fidelity
strongly depends on the nature of the initial state. Indeed, for $\Lambda\leq
5$, the results obtained with sampling noise are globally comparable to the
ones from an ideal VQE. Interestingly, when the monomodal initial state is
used, an even better fidelity is reached in this regime of weak correlation.
However, for $\Lambda=5$, the resulting fidelities deteriorate for both
initial states, especially when the number of bosons increases. The
deterioration is even more important when considering the strong interaction
regime $\Lambda=10$. In this case, a vertical limit separates the high- and
low-fidelity regions at $N_{B}=4$ when the monomodal initial state is used.
Over this limits, the VQE trial state presents a fidelity of $\mathcal{F}\sim
50\%$ because it gets stuck in a single configuration of the cat-state during
the optimization (i.e. the converged state was
$\ket{\Psi(\vec{\theta})}\propto\ket{N_{B},0}$ instead of
$\ket{\Psi_{0}}\propto\ket{N_{B},0}+\ket{0,N_{B}}$). Therefore, cat-like
states seem difficult to realize by the BS-Kerr ansatz when sampling noise is
considered. However, much better accuracy is reached when the bimodal initial
state is used instead. Thus suggesting that starting from a configuration with
photons spread in between the modes can improve the results of the VQE in the
strongly correlated regime ($\Lambda\geq 5$ here). In summary, it seems from
Fig. 16 that – within the BS-Kerr ansatz – a weakly-correlated state is more
easily encoded when starting from a monomodal state rather than a bimodal one,
and vice versa for a strongly-correlated state.

##  VI Conclusions

In this work, we tackled the problem of encoding strongly correlated many-
boson ground-state wavefunctions on a photonic quantum device. We designed two
different quantum photonic-based ansatze and simulated a variational quantum
algorithm (VQA) using the ground-state infidelity as a cost function. Both
ansatze, the beam-splitter-Kerr (BS-Kerr) and the interferometer-Kerr, are
based on layers of beam-splitter and Kerr gates, with additional rotation
gates for the latter. We showed their ability to encode the ground state of
the attractive BH model in all correlation regimes, for different number of
bosons and initial states. While the BS-Kerr ansatz has shown to be the most
efficient in terms of total number of gates and parameters, it might not be
optimal for conducting a true experiment on a real photonic quantum device, as
it contains many Kerr gates that are difficult to implement in practice. In
contrast, the interferometer-Kerr ansatz was designed to reduce the number of
Kerr gates (or the number of layers required to achieve a given ground-state
fidelity) at the expense of more rotation and beam-splitter gates and more
parameters per layer, and is therefore more appropriate for realistic
applications. Although taking the fidelity as a criteria is theoretically
sound to test the encoding ability of our ansatze, in practice the cost
function is the ground-state energy measured within the variational quantum
eigensolver (VQE) algorithm. As a last investigation, we performed realistic
VQE simulations using the BS-Kerr ansatz, with and without sampling noise, and
proposed a photon-counting protocol to measure the ground-state energy of the
BH Hamiltonian as it would be done in a real experiment. An even more
realistic simulation should of course account for photon loss, on which error
mitigation schemes can be applied [85]. This is left for future work.

## Acknowledgments

SY, BS and VD thank Filippo Miatto for fruitful discussions. SY and BS
sincerely thank Eleanor Scerri, Vincent Pouthier and Xavi Bonet-Monroig for
fruitful discussions. SY and BS acknowledge support from the Netherlands
Organization for Scientific Research (NWO/OCW). VD also acknowledges support
by the Dutch Research Council (NWO/OCW), as part of the Quantum Software
Consortium program (project number 024.003.037).

##  Appendix A Exact encoding of the BH ground state using only beam-splitters
in the non-interacting limit

We demonstrate here that a single layer of the BS-Kerr ansatz, without the
Kerr gates (so only beam-splitters) can exactly encode the ground state of the
BH model for the three newtorks considered in our study (i.e. with $N_{S}=2$,
3 and 4 sites) in the non-interacting limit $\Lambda\rightarrow 0$. Note that
this demonstration works only for initial monomodal states such as
$\ket{N_{B},0}$, $\ket{N_{B},0,0}$ and $\ket{N_{B},0,0,0}$ for the dimer,
three- and four-site BH model, respectively.

Let us start with the dimer, for which the exact ground state reads

$\ket{\Psi_{0}}=\frac{1}{\sqrt{2N_{B}!}}[a_{1}^{\dagger}+a_{2}^{\dagger}]^{N_{B}}\ket{0,0}.$
(26)

The action of a single layer of beam-splitter gates on an initial monomodal
state results in

$B_{1,2}(\theta,0)\ket{N_{B},0}=\frac{1}{\sqrt{N_{B}!}}[\cos(\theta)a_{0}^{\dagger}+\sin(\theta)a_{1}^{\dagger}]^{N_{B}}\ket{0,0}$
(27)

which provides an exact encoding of the ground state [Eq. (26)] if we fix the
beam-splitter transmittivity parameter to $\theta=\pi/4$. Turning to the
three-site model, its exact ground state is given by

$\ket{\Psi_{0}}=\frac{1}{\sqrt{3N_{B}!}}[a_{1}^{\dagger}+a_{2}^{\dagger}+a_{3}^{\dagger}]^{N_{B}}\ket{0,0,0}.$
(28)

applying the stair structure of the beam-splitters transform the initial
monomodal state as follows,

$\begin{split}&B_{2,3}(\theta^{\prime},0)B_{1,2}(\theta,0)\ket{N_{B},0,0}=\\\
&\frac{1}{\sqrt{N_{B}!}}\Big{[}\cos(\theta)a_{1}^{\dagger}+\sin(\theta)\cos(\theta^{\prime})a_{2}^{\dagger}\\\
&+\sin(\theta)\sin(\theta^{\prime})a_{3}^{\dagger}\Big{]}^{N_{B}}\ket{0,0,0},\end{split}$
(29)

thus realizing an exact encoding of the ground state [Eq. (28)] when we fix
the beam-splitter’s parameters such that $\theta=\arccos(1/\sqrt{3})$ and
$\theta^{\prime}=\pi/4$. Finally, the exact ground state of the four-site
model reads

$\ket{\Psi_{0}}=\frac{1}{\sqrt{4N_{B}!}}[a_{1}^{\dagger}+a_{2}^{\dagger}+a_{3}^{\dagger}+a_{4}^{\dagger}]^{N_{B}}\ket{0,0,0,0}.$
(30)

Similarly as for the BH dimer and three-site models, we apply a single layer
of the BS-Kerr ansatz (without the Kerr) on the monomodal initial state, thus
leading to

$\begin{split}&B_{3,4}(\theta^{\prime\prime},0)B_{2,3}(\theta^{\prime},0)B_{1,2}(\theta,0)\ket{N_{B},0,0,0}=\\\
&\frac{1}{\sqrt{N_{B}!}}\Big{[}\cos(\theta)a_{1}^{\dagger}+\sin(\theta)\cos(\theta^{\prime})a_{2}^{\dagger}\\\
&+\sin(\theta)\sin(\theta^{\prime})\cos(\theta^{\prime\prime})a_{3}^{\dagger}\\\
&+\sin(\theta)\sin(\theta^{\prime})\sin(\theta^{\prime\prime})a_{4}^{\dagger}\big{)}\big{)}\Big{]}^{N_{B}}\ket{0,0,0,0}.\end{split}$
(31)

Fixing the beam-splitter’s parameters to $\theta=\pi/3$,
$\theta^{\prime}=\arccos(1/\sqrt{3})$ and $\theta^{\prime\prime}=\pi/4$
generates an exact representation of the exact ground state in Eq.(30).

## References

  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’brien, Nature Comm. 5, 4213 (2014).
  * O’Malley  _et al._ [2016] P. J. O’Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, _et al._ ,  6, 031007 (2016).
  * Shen  _et al._ [2017] Y. Shen, X. Zhang, S. Zhang, J.-N. Zhang, M.-H. Yung, and K. Kim, Phys. Rev. A 95, 020501 (2017).
  * Kandala  _et al._ [2017] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242 (2017).
  * Hempel  _et al._ [2018] C. Hempel, C. Maier, J. Romero, J. McClean, T. Monz, H. Shen, P. Jurcevic, B. P. Lanyon, P. Love, R. Babbush, A. Aspuru-Guzik, R. Blatt, and C. F. Roos,  8, 031022 (2018).
  * Colless  _et al._ [2018] J. I. Colless, V. V. Ramasesh, D. Dahlen, M. S. Blok, M. Kimchi-Schwartz, J. McClean, J. Carter, W. De Jong, and I. Siddiqi,  8, 011021 (2018).
  * Mazzola  _et al._ [2019] G. Mazzola, P. J. Ollitrault, P. K. Barkoutsos, and I. Tavernelli, Phys. Rev. Lett. 123, 130501 (2019).
  * Nam  _et al._ [2020] Y. Nam, J.-S. Chen, N. C. Pisenti, K. Wright, C. Delaney, D. Maslov, K. R. Brown, S. Allen, J. M. Amini, J. Apisdorf, _et al._ , npj Quantum Inf. 6, 1 (2020).
  * O’Brien  _et al._ [2019] T. E. O’Brien, B. Senjean, R. Sagastizabal, X. Bonet-Monroig, A. Dutkiewicz, F. Buda, L. DiCarlo, and L. Visscher, npj Quantum Inf. 5, 1 (2019).
  * Sokolov  _et al._ [2021] I. O. Sokolov, P. K. Barkoutsos, L. Moeller, P. Suchsland, G. Mazzola, and I. Tavernelli, Phys. Rev. Research 3, 013125 (2021).
  * Arute  _et al._ [2020] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, S. Boixo, M. Broughton, B. B. Buckley, D. A. Buell, B. Burkett, N. Bushnell, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, S. Demura, A. Dunsworth, E. Farhi, A. Fowler, B. Foxen, C. Gidney, M. Giustina, R. Graff, S. Habegger, M. P. Harrigan, A. Ho, S. Hong, T. Huang, W. J. Huggins, L. Ioffe, S. V. Isakov, E. Jeffrey, Z. Jiang, C. Jones, D. Kafri, K. Kechedzhi, J. Kelly, S. Kim, P. V. Klimov, A. Korotkov, F. Kostritsa, D. Landhuis, P. Laptev, M. Lindmark, E. Lucero, O. Martin, J. M. Martinis, J. R. McClean, M. McEwen, A. Megrant, X. Mi, M. Mohseni, W. Mruczkiewicz, J. Mutus, O. Naaman, M. Neeley, C. Neill, H. Neven, M. Y. Niu, T. E. O’Brien, E. Ostby, A. Petukhov, H. Putterman, C. Quintana, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, D. Strain, K. J. Sung, M. Szalay, T. Y. Takeshita, A. Vainsencher, T. White, N. Wiebe, Z. J. Yao, P. Yeh, and A. Zalcman, Science 369, 1084 (2020).
  * Cade  _et al._ [2020] C. Cade, L. Mineh, A. Montanaro, and S. Stanisic, Physical Review B 102, 235122 (2020).
  * Fujii  _et al._ [2020] K. Fujii, K. Mitarai, W. Mizukami, and Y. O. Nakagawa, arXiv preprint arXiv:2007.10917 (2020).
  * Xu  _et al._ [2020] L. Xu, J. T. Lee, and J. Freericks, Modern Physics Letters B 34, 2040049 (2020).
  * Ollitrault  _et al._ [2020] P. J. Ollitrault, A. Baiardi, M. Reiher, and I. Tavernelli, Chemical science 11, 6842 (2020).
  * Romero and Aspuru-Guzik [2021] J. Romero and A. Aspuru-Guzik, Advanced Quantum Technologies 4, 2000003 (2021).
  * Grimsley  _et al._ [2018] H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, arXiv preprint arXiv:1812.11173 (2018).
  * Tang  _et al._ [2021] H. L. Tang, V. Shkolnykov, G. Barron, H. Grimsley, N. Mayhall, E. Barnes, and S. Economou, Bulletin of the American Physical Society (2021).
  * Montanaro and Stanisic [2020] A. Montanaro and S. Stanisic, arXiv preprint arXiv:2006.01179 (2020).
  * Masaki-Kato  _et al._ [2019] A. Masaki-Kato, S. Yunoki, and D. Hirashima, Physical Review B 100, 224515 (2019).
  * Krauth [1996] W. Krauth, Physical review letters 77, 3695 (1996).
  * Del Maestro and Affleck [2010] A. Del Maestro and I. Affleck, Physical Review B 82, 060515 (2010).
  * Purwanto and Zhang [2004] W. Purwanto and S. Zhang, Physical Review E 70, 056702 (2004).
  * Pollet [2013] L. Pollet, Comptes Rendus Physique 14, 712 (2013).
  * Bogner  _et al._ [2019] B. Bogner, C. De Daniloff, and H. Rieger, The European Physical Journal B 92, 1 (2019).
  * Łącki  _et al._ [2016] M. Łącki, B. Damski, and J. Zakrzewski, Scientific reports 6, 1 (2016).
  * Kawaki  _et al._ [2017] K. Kawaki, Y. Kuno, and I. Ichinose, Physical Review B 95, 195101 (2017).
  * Dogra  _et al._ [2016] N. Dogra, F. Brennecke, S. D. Huber, and T. Donner, Physical Review A 94, 023632 (2016).
  * Lv  _et al._ [2014] J.-P. Lv, Q.-H. Chen, and Y. Deng, Physical Review A 89, 013628 (2014).
  * Rossi  _et al._ [2009] M. Rossi, M. Nava, L. Reatto, and D. Galli, The Journal of chemical physics 131, 154108 (2009).
  * Wilkinson and Hartmann [2020] S. A. Wilkinson and M. J. Hartmann, Applied Physics Letters 116, 230501 (2020).
  * Roushan  _et al._ [2017] P. Roushan, C. Neill, J. Tangpanitanon, V. M. Bastidas, A. Megrant, R. Barends, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. Fowler, B. Foxen, M. Giustina, E. Jeffrey, J. Kelly, E. Lucero, J. Mutus, M. Neeley, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. White, H. Neven, D. G. Angelakis, and J. Martinis, Science 358, 1175 (2017), https://science.sciencemag.org/content/358/6367/1175.full.pdf .
  * Hacohen-Gourgy  _et al._ [2015] S. Hacohen-Gourgy, V. V. Ramasesh, C. De Grandi, I. Siddiqi, and S. M. Girvin, Physical review letters 115, 240501 (2015).
  * Deng  _et al._ [2016] X.-H. Deng, C.-Y. Lai, and C.-C. Chien, Physical Review B 93, 054116 (2016).
  * Huh  _et al._ [2015] J. Huh, G. G. Guerreschi, B. Peropadre, J. R. McClean, and A. Aspuru-Guzik, Nature Photonics 9, 615 (2015).
  * Huh and Yung [2017] J. Huh and M.-H. Yung, Scientific reports 7, 1 (2017).
  * Sparrow  _et al._ [2018] C. Sparrow, E. Martín-López, N. Maraviglia, A. Neville, C. Harrold, J. Carolan, Y. N. Joglekar, T. Hashimoto, N. Matsuda, J. L. O’Brien, _et al._ , Nature 557, 660 (2018).
  * Quesada [2019] N. Quesada, J. Chem. Phys. 150, 164113 (2019).
  * Jahangiri  _et al._ [2020] S. Jahangiri, J. M. Arrazola, N. Quesada, and A. Delgado, Phys. Chem. Chem. Phys. 22, 25528 (2020).
  * Kalajdzievski  _et al._ [2018] T. Kalajdzievski, C. Weedbrook, and P. Rebentrost, Phys. Rev. A 97, 062311 (2018).
  * Zhong  _et al._ [2020] 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._ , Science 370, 1460 (2020).
  * Weedbrook  _et al._ [2012] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Rev. Mod. Phys. 84, 621 (2012).
  * Bromley  _et al._ [2020] T. R. Bromley, J. M. Arrazola, S. Jahangiri, J. Izaac, N. Quesada, A. D. Gran, M. Schuld, J. Swinarton, Z. Zabaneh, and N. Killoran, Quantum Sci. Technol. 5, 034010 (2020).
  * Pfister [2019] O. Pfister, Journal of Physics B: Atomic, Molecular and Optical Physics 53, 012001 (2019).
  * Kühner and Monien [1998] T. D. Kühner and H. Monien, Phys. Rev. B 58, R14741 (1998).
  * Kühner  _et al._ [2000] T. D. Kühner, S. R. White, and H. Monien, Phys. Rev. B 61, 12474 (2000).
  * Jack and Yamashita [2005] M. W. Jack and M. Yamashita, Phys. Rev. A 71, 023610 (2005).
  * Oelkers and Links [2007] N. Oelkers and J. Links, Physical Review B 75, 115119 (2007).
  * Juliá-Díaz  _et al._ [2010] B. Juliá-Díaz, D. Dagnino, M. Lewenstein, J. Martorell, and A. Polls, Physical Review A 81, 023615 (2010).
  * Melé-Messeguer  _et al._ [2011] M. Melé-Messeguer, B. Juliá-Díaz, and A. Polls, Journal of Low Temperature Physics 165, 180 (2011).
  * Mansikkamäki  _et al._ [2021] O. Mansikkamäki, S. Laine, and M. Silveri, arXiv:2101.06032 (2021).
  * Theis  _et al._ [2004] M. Theis, G. Thalhammer, K. Winkler, M. Hellwig, G. Ruff, R. Grimm, and J. H. Denschlag, Physical Review Letters 93, 123001 (2004).
  * Góral  _et al._ [2002] K. Góral, L. Santos, and M. Lewenstein, Phys. Rev. Lett. 88, 170406 (2002).
  * Büchler and Blatter [2003] H. Büchler and G. Blatter, Phys. Rev. Lett. 91, 130404 (2003).
  * Bloch [2005] I. Bloch, Nature Phys. 1, 23 (2005).
  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * Baier  _et al._ [2016] S. Baier, M. J. Mark, D. Petter, K. Aikawa, L. Chomaz, Z. Cai, M. Baranov, P. Zoller, and F. Ferlaino, Science 352, 201 (2016).
  * Pouthier [2003] V. Pouthier, Phys. Rev. E 68, 021909 (2003).
  * Edler  _et al._ [2004] J. Edler, R. Pfister, V. Pouthier, C. Falvo, and P. Hamm, Phys. Rev. Lett. 93, 106405 (2004).
  * Raventós  _et al._ [2017] D. Raventós, T. Graß, M. Lewenstein, and B. Juliá-Díaz, Journal of Physics B: Atomic, Molecular and Optical Physics 50, 113001 (2017).
  * Zhang and Dong [2010] J. Zhang and R. Dong, European Journal of Physics 31, 591 (2010).
  * Bauer and Nayak [2013] B. Bauer and C. Nayak, Journal of Statistical Mechanics: Theory and Experiment 2013, P09005 (2013).
  * Hopjan and Heidrich-Meisner [2020] M. Hopjan and F. Heidrich-Meisner, Physical Review A 101, 063617 (2020).
  * Giraud  _et al._ [2009] O. Giraud, J. Martin, and B. Georgeot, Physical Review A 79, 032308 (2009).
  * Giraud  _et al._ [2007] O. Giraud, J. Martin, and B. Georgeot, Physical Review A 76, 042333 (2007).
  * Lukin  _et al._ [2019] A. Lukin, M. Rispoli, R. Schittko, M. E. Tai, A. M. Kaufman, S. Choi, V. Khemani, J. Léonard, and M. Greiner, Science 364, 256 (2019).
  * Islam  _et al._ [2015] R. Islam, R. Ma, P. M. Preiss, M. E. Tai, A. Lukin, M. Rispoli, and M. Greiner, Nature 528, 77 (2015).
  * Beugeling  _et al._ [2015] W. Beugeling, A. Andreanov, and M. Haque, Journal of Statistical Mechanics: Theory and Experiment 2015, P02002 (2015).
  * Cerezo  _et al._ [2020] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, arXiv:2012.09265 (2020).
  * Evangelista  _et al._ [2019] F. A. Evangelista, G. K.-L. Chan, and G. E. Scuseria, The Journal of chemical physics 151, 244112 (2019).
  * Sokolov  _et al._ [2020] I. O. Sokolov, P. K. Barkoutsos, P. J. Ollitrault, D. Greenberg, J. Rice, M. Pistoia, and I. Tavernelli, The Journal of chemical physics 152, 124107 (2020).
  * Mizukami  _et al._ [2020] W. Mizukami, K. Mitarai, Y. O. Nakagawa, T. Yamamoto, T. Yan, and Y.-y. Ohnishi, Physical Review Research 2, 033421 (2020).
  * Yalouz  _et al._ [2021] S. Yalouz, B. Senjean, J. Günther, F. Buda, T. E. O’Brien, and L. Visscher, Quantum Science and Technology 6, 024004 (2021).
  * Boyd [2020] R. W. Boyd,  _Nonlinear optics_ (Academic press, 2020).
  * Azuma [2007] H. Azuma, J. Phys. D: Appl. Phys. 41, 025102 (2007).
  * Clements  _et al._ [2016] W. R. Clements, P. C. Humphreys, B. J. Metcalf, W. S. Kolthammer, and I. A. Walmsley, Optica 3, 1460 (2016).
  * Killoran  _et al._ [2019] N. Killoran, J. Izaac, N. Quesada, V. Bergholm, M. Amy, and C. Weedbrook, Quantum 3, 129 (2019).
  * Gerber and Furrer [2019] F. Gerber and R. Furrer, The R Journal 11, 352 (2019).
  * Jönsson and Björk [2019] M. Jönsson and G. Björk, Physical Review A 99, 043822 (2019).
  * Provazník  _et al._ [2020] J. Provazník, L. Lachman, R. Filip, and P. Marek, Optics express 28, 14839 (2020).
  * Chen  _et al._ [2018] W.-P. Chen, J. Singleton, L. Qin, A. Camón, L. Engelhardt, F. Luis, R. E. Winpenny, and Y.-Z. Zheng, Nature communications 9, 1 (2018).
  * Rossini and Fazio [2012] D. Rossini and R. Fazio, New J. Phys. 14, 065012 (2012).
  * Ohgoe  _et al._ [2012] T. Ohgoe, T. Suzuki, and N. Kawashima, Physical Review B 86, 054520 (2012).
  * Hansen [2016] N. Hansen, arXiv preprint arXiv:1604.00772 (2016).
  * Su  _et al._ [2020] D. Su, R. Israel, K. Sharma, H. Qi, I. Dhand, and K. Brádler, arXiv:2008.06670 (2020).