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

Saad Yalouz, Bruno Senjean, Filippo Miatto, Vedran Dunjko

► To cite this version:

Saad Yalouz, Bruno Senjean, Filippo Miatto, Vedran Dunjko. Encoding strongly-correlated many-boson wavefunctions on a photonic quantum computer: application to the attractive Bose-Hubbard model. Quantum, 2021, 5, pp.572. 10.22331/q-2021-11-08-572. hal-03796321

HAL Id: hal-03796321

https://hal.science/hal-03796321

Submitted on 5 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Encoding strongly-correlated many-boson wavefunctions on a photonic quantum computer: application to the attractive Bose-Hubbard model

Saad Yalouz<sup>1,2,3</sup>, Bruno Senjean<sup>4,3</sup>, Filippo Miatto<sup>5</sup>, and Vedran Dunjko<sup>6</sup>

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 accurately 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.

## 1 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 treat 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 algorithms (VQA) have shown to be very efficient and promising for encoding complex wavefunctions of various kinds of systems. Applications of these algorithms range from molecular systems in quantum chemistry [1–11], interacting-spins and electrons models from condensed matter [12–14] or vibrational Hamiltonians [15], to cite but a few.

In practice, the success of a VQA in encoding a many-body wavefunction depends on the properties of the associated quantum circuit, or “ansatz”, that is used. Ideally, to satisfy the needs and limitations of the NISQ era, an ansatz should be expressive (i.e. for example, having the ability to generate accurate ground-state wavefunc-

Saad Yalouz: yalouzsaad@gmail.com

Bruno Senjean: bruno.senjean@umontpellier.fr

Filippo Miatto: filippo@xanadu.ai

Vedran Dunjko: v.dunjko@liacs.leidenuniv.nl

<footnote n="1"><sup>1</sup>Laboratoire de Chimie Quantique, CNRS/Université de Strasbourg, 4 rue Blaise Pascal, 67000 Strasbourg, France</footnote>

<footnote n="2"><sup>2</sup>Theoretical Chemistry, Vrije Universiteit, De Boelelaan 1083, NL-1081 HV, Amsterdam, The Netherlands</footnote>

<footnote n="3"><sup>3</sup>Instituut-Lorentz, Universiteit Leiden, P.O. Box 9506, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="4"><sup>4</sup>ICGM, Univ Montpellier, CNRS, ENSCM, Montpellier, France</footnote>

<footnote n="5"><sup>5</sup>Xanadu, Toronto, ON, M5G 2C8, Canada</footnote>

<footnote n="6"><sup>6</sup>Universiteit Leiden, P.O. Box 9506, 2300 RA Leiden, The Netherlands</footnote>
tions or energies), resource efficient (i.e. involv-
ing as few gates/layers as possible) and trainable
(i.e. inducing optimization landscapes amenable
to classical optimization algorithms). In the case
of qubit-based quantum circuits, the above prop-
erties of ansatze for VQAs have been studied
extensively: examples are the unitary coupled-
cluster ansatz [16] introduced for electronic struc-
ture problems, adaptive methods designed to
produce expressive ansatze at lower complexity
cost [17,18], and different strategies to solve con-
densed matter models with strong quantum cor-
relation [12, 19]. The great majority of these
works has been essentially devoted to the case
of many-fermion systems. In contrast, lit-
tle attention has been paid to study the po-
tential of current devices to produce interest-
ing ansatze for the study of strongly interacting
many-boson systems. On the classical computing
side, state-of-the-art methods such as quantum
Monte-Carlo (QMC) simulations have proven to
be effective [20–30], as they do not suffer from
the infamous sign problem of fermionic systems.
However, QMC is usually restricted to ground-
state calculations, although extensions to molec-
ular properties [31–35] and excited states [36–40]
have been developed. The rise of quantum com-
puting can open alternative routes to ground-
state calculations (and beyond) for many-boson
systems, for instance based on superconducting
circuits [41–45], cold-atom lattices [46, 47] or pho-
tonic quantum devices [48–55] (such as Gaussian
boson sampling experiments [48, 49, 52, 53]), es-
sentially used as analog quantum simulators. In-
deed, many extensions of VQAs have been pro-
posed recently to compute molecular and excited-
state properties [9, 56–72], and one can expect
similar developments for many-boson systems in
the near future.

Motivated by the development of VQAs, we in-
vestigate in this work the use of digital photonic-
based continuous variable (CV) quantum com-
puters [73-76] to encode strongly correlated
many-boson wavefunctions. In practice, CV de-
vices are digital quantum computers that can ma-
nipulate information by entangling and measur-
ing 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 sys-

tems. Such a mapping was first employed in
Ref. [48] to generate molecular vibronic spectra,
where the *N* photons in *M* optical modes can
simulate *N* molecular vibrational quanta in *M*
vibrational modes. Bosonic modes of CV de-
vices were also used to efficiently encode the field
degrees of freedom in a wavelet basis in quan-
tum field theories [77, 78]. For variational algo-
rithms, 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 at-
tracted a lot of attention recently due to its
rich transitional behavior that evolves between
the so-called Schrödinger cat-like states to su-
perfluid regimes [79-85]. Furthermore, the BH
model has been successfully employed to describe
the physics of ultracold atoms and molecules in
optical lattices [86-91], quantum cooling proto-
cols [44], the evolution of vibrational bosons in
α-helices proteins [92, 93] and the exotic phases
of helium-4 [94, 95].

In this paper, we study the attractive BH
model as an interesting system and focus on two
aspects of the aforementioned critical properties
of ansatze for VQA applications: the express-
ibility and the resource efficiency to encode the
ground-state wavefunction of the system on CV
devices. Note that, in analogy with the qubit
cases, the implementation costs of a photonic cir-
cuit depend on numerous factors. Obviously, they
depend on depth and total gate number, but not
all gates are equally costly (e.g. linear optical ele-
ments are much cheaper than gates involving non-
linear terms), nor they introduce an equal amount
of noise and losses. We are thus facing a multi-
objective optimization problem as, for instance,
using cheaper gates may come at the expense of
greater depth to achieve the same theoretical ca-
pacity to represent ground states accurately. Fur-
ther, in the light of near-term implementations, it
makes additional sense to study ansatz tailored to
the actual experiments performed in the labs or
to those available in the near future. To elucidate
the issues of trade-offs between various resource
requirements and expressibility, we propose and
analyze two CV quantum circuit architectures
tailored to near-future experimental set-ups that
can make use of single-mode non-linear phase shifts such as the Kerr gate: one designed to minimize the total gate count and the other designed to minimize the total circuit depth. By performing VQAs with the state-infidelity as a cost function, we show that the many-boson ground state of the BH model can be successfully encoded on both CV quantum circuits with a high fidelity and a relatively small number of gates and parameters. In practice, 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 device would require the ability to measure on a CV circuit 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 diagonalization method in sections 2.1 and 2.2, we illustrate the exotic structural transition of the BH ground state on three small-sized networks in section 2.3. After introducing the context of VQAs in section 3.1, we present the architectures of our proposed ansatze in section 3.2. Numerical investigations of their encoding ability are performed in section 4. We extend our study to ideal and realistic simulations of a VQE algorithm in section 5. Finally, conclusions and perspectives are provided in section 6.

## 2 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.

### 2.1 Hamiltonian and regimes of interaction

The attractive BH model describes the physics of $N_B$ bosons hopping between $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} (a_p^\dagger a_q + a_q^\dagger a_p) - \frac{U}{2} \sum_p n_p (n_p - 1), \quad (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} \quad (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).

### 2.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}} |\Psi_0\rangle = E_0 |\Psi_0\rangle, \quad (3) $$

where $|\Psi_0\rangle$ is the ground state of the system with the associated energy $E_0$. We solve Eq. (3) by
<table><caption>Table 1: Dimension $D_{\hat{H}}$ of the Fock space of the BH model as a function of the number of sites $N_S$ and bosons $N_B$.</caption><thead><tr><th>N<sub>B</sub><br>N<sub>S</sub></th><th>2</th><th>3</th><th>4</th><th>5</th><th>8</th><th>16</th></tr></thead><tbody><tr><th>2</th><td>3</td><td>4</td><td>5</td><td>6</td><td>9</td><td>17</td></tr><tr><th>3</th><td>6</td><td>10</td><td>15</td><td>21</td><td>45</td><td>153</td></tr><tr><th>4</th><td>10</td><td>20</td><td>35</td><td>56</td><td>165</td><td>969</td></tr><tr><th>8</th><td>36</td><td>120</td><td>330</td><td>792</td><td>6435</td><td>245157</td></tr></tbody></table>

exact diagonalization (ED) to obtain a good reference state to compare with our quantum ansatz developed for CV devices. Note that the way ED works in theory may appear 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. But beyond its apparent simplicity, the method requires in practice some technical efforts of implementations (see Refs. [96,97] for an introduction of the numerical approach). In the case of the attractive BH model, the dimension $D_{\hat{H}}$ of the resulting matrix is equal to the dimension of the Fock space of the system,

$$D_{\hat{H}} = \frac{(N_B + N_S - 1)!}{N_B!(N_S - 1)!}, \quad (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{H}}$ as a function of $N_S$ and $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.

## 2.3 Illustration of the regimes of interaction: exact diagonalization on small-sized 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 [82, 83, 85]). 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. To do so, we employ the Inverse Participation Ratio (IPR) [98–101] and the

von Neumann entropy [102–104] as two indicators to follow the expansion of $|\Psi_0\rangle$ in the configuration basis when $\Lambda$ increases. Applied to the Fock state basis, the IPR estimates the number of different bosonic configurations forming $|\Psi_0\rangle$ by linear combination, and reads

$$\mathrm{IPR}(|\Psi_0\rangle) = \left( \sum_{n_1, \dots, n_{N_S}} |\langle n_1, \dots, n_{N_S} | \Psi_0 \rangle|^4 \right)^{-1}, \quad (5)$$

where $|n_1, \dots, n_{N_S}\rangle$ is a Fock state with $n_q$ bosons on site $q$ (and $\sum_q n_q = N_B$). Within this definition, a state $|\Psi_0\rangle$ 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{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), \quad (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 |n_k\rangle \langle n_k|, \quad (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 = \mathrm{Tr}_{q \neq p}[\hat{\rho}]$ (note that in our case, the full density operator is $\hat{\rho} = |\Psi_0\rangle \langle \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).

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 not interact together. The independent-boson picture leads to a product-state solution for the ground state as $|\Psi_0\rangle = \frac{1}{\sqrt{2\sqrt{N_B}}} (a_1^\dagger + a_2^\dagger)^{N_B} |0,0\rangle$ (with $|0,0\rangle$ being the vacuum state of the dimer). Then, increasing $\Lambda$ leads to a progressive restructuring of the ground state as readily seen by the
<figure><img src="image_2.png" /><figcaption>**Figure 1:** **Ground-state properties of the BH dimer.** Evolution of the ground-state structure as a function of $\Lambda$ (obtained with ED). **Upper panel:** IPR values for different total number of bosons $N_B = 3, \dots, 8$. **Middle panel:** Single site entropy $S$ values for different values of the total number of bosons $N_B = 3, \dots, 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 to black for 40%. In all 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.</figcaption></figure>

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 \le \Lambda \le 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 [82, 83] representing a superposition of two macroscopically occupied states (i.e. $|\Psi_0\rangle = (|N_B, 0\rangle + |0, N_B\rangle)/\sqrt{2}$). The IPR reaches a value of 2, thus indicating the im-

<figure><img src="image_3.png" /><figcaption>**Figure 2:** **Ground-state properties of the three- and four-site BH periodic chains.** Evolution of the IPR and the single-site entropy $S$ of the ground state as a function of $N_B$ and $\Lambda$ (obtained by ED). **Upper panel:** Three-site BH periodic chain. **Lower panel:** Four-site BH periodic chain. In both panels, 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.</figcaption></figure>

portance 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 figure 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 IPR$^{max}$ ≈ 42 bosonic configurations involved in the structure of the ground state at the maximum value of the IPR ($\Lambda$ ≈ 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$ ≈ 4 with a total number of IPR$^{max}$ ≈ 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 (IPR = 3 and 4 for the three and four-site systems, respectively, e.g. for the three-site model we have $|\Psi_0\rangle = (|N_B, 0, 0\rangle + |0, N_B, 0\rangle + |0, 0, N_B\rangle)/\sqrt{3}$).

# 3 Photonic-gate-based ansatz

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

## 3.1 Variational Quantum Algorithms

For the study of many-body systems in the NISQ era, VQAs represent state-of-the-art methods that have been proficiently used to determine ground (and excited) states' properties of various systems (see for example Refs. [12,16–19] and the recent review in Ref. [105]). 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 $\hat{U}(\vec{\theta})$ with tunable gate parameters (formally represented by the vector $\vec{\theta}$) to transform an initial state $|\Psi_{ini}\rangle$ into a trial state

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

Often $|\Psi_{ini}\rangle$ is chosen based on some physical motivations which depends on the targeted problem (e.g. 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 energy of the trial state defined as

$$E(\vec{\theta}) = \langle \Psi(\vec{\theta}) | \hat{H} | \Psi(\vec{\theta}) \rangle \ge E_0, \qquad (9)$$

where $E_0$ is the exact ground-state energy which represents a natural lower bound for the energy of the trial state $|\Psi(\vec{\theta})\rangle$. 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).

In the NISQ era, the overall performance of a VQA in encoding a many-body wavefunction is highly dependent on the properties of the ansatz that is considered (as mentioned in the introduction: expressibility, resource efficiency and trainability). In our study of CV ansatz, our attention turns towards the expressibility and resource efficiency, while the investigation of the trainability is left for future work. However, these properties are not independent, and the final performance should include the specification of the optimization as it may heavily influence the expressibility and resource efficiency of the VQA to achieve a given precision of the targeted state or energy. For this reason, studying all three properties is very difficult.

In the following, we introduce two different ansatz architectures and we apply them to the three BH systems introduced before. The different architectures address (to some extent) the question of expressibility and resource efficiency. For this, we characterize the capability of the associated circuit unitaries $\hat{U}(\vec{\theta})$ (see Eq. (8)) to accurately map different initial states $|\Psi_{ini}\rangle$ to the exact ground state $|\Psi_0\rangle$ of the three systems. Note that only a few works have tried to investigate the expressibility of ansatzes for the case of many-body systems (always with a special emphasis for qubit devices and many-electron systems, e.g. the study of the unitary coupled-cluster-ansatz in [106]). As a results, a lot of attention has been directed to many-fermion systems (see for example [12,17–19,72,107,108]), but little work has been done on quantum ansatz to encode many-boson wavefunctions (for a notable exception, see the implementation of the unitary vibrational-coupled-cluster ansatz introduced in Ref. [15]).
## 3.2 Photonic quantum circuits as ansatze for many-boson wavefunction

In contrast to most of the studies focusing on qubit devices, we tackle the problem of many-bosons systems using CV circuits. The choice of CV quantum computers is rather intuitive as there exists a straightforward one-to-one correspondence between the modes of a bosonic system and the photonic modes of such devices. 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 introduce two quantum ansatze based on CV photonic gates to encode many-boson wavefunctions for VQAs.

<figure><img src="image_4.png" /><figcaption>Figure 3: Illustration of the VQA approach. Each photonic mode in the quantum circuit is used to encode the bosonic occupation of a particular site in the Bose-Hubbard network.</figcaption></figure>

### 3.2.1 Utilized CV quantum gates

Concerning the fundamental properties of both ansatze, we know that the BH Hamiltonian in Eq. (1) is boson-number-preserving by construction (i.e. $N_B$ is constant). Thus, in order to preserve this property, we choose for both ansatze to consider only three types of photon-number-preserving gates, namely: beam-splitters, rotation gates and Kerr gates. First, the beam-splitter gate is a two-mode gate that affects two inputs modes simultaneously (labeled by p and q) with an unitary defined as

$$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), \quad (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

<figure><img src="image_5.png" /><figcaption>Figure 4: Illustration of the effects of a BS-gate. Density of probability of a state $B(\theta, 0)|\Psi_{ini}\rangle$ (with $|\Psi_{ini}\rangle$ an initial state) 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: $|\Psi_{ini}\rangle = |4, 4\rangle$. Lower panel: $|\Psi_{ini}\rangle = |8, 0\rangle$.</figcaption></figure>

operators ($a_p/a_p^{\dagger}$ and $a_q/a_q^{\dagger}$) is given by

$$ \begin{aligned} 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{aligned} \quad (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)|\Psi_{ini}\rangle$ (with $|\Psi_{ini}\rangle = |8, 0\rangle$ or $|4, 4\rangle$) 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. This strongly depends on the initial state considered. The phase rotation gate acts on a Fock state as $R_p(\theta)|n_p\rangle = e^{i\theta n_p}|n_p\rangle$ with the unitary operator

$$R_p(\theta) = \exp \left( i\theta a_p^{\dagger} a_p \right). \quad (12)$$

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

$$K_p(\theta) = \exp \left( i\theta (a_p^{\dagger} a_p)^2 \right). \quad (13)$$

For a given modes p, this leads to $K_p(\theta)|n_p\rangle = e^{i\theta n_p^2}|n_p\rangle$, thus generating a non-linear phase term that depends on the squared number of
bosons $n_p^2 = (a_p^\dagger a_p)^2$ present in the mode $p$. Although the implementation of Kerr gates as an optical component is a true challenge on current devices, the incessant development of new techniques and non-linear materials tends to suggest that reliable ways of implementing these gates will be available in the near future [109–111]. For the purpose of our theoretical work, we assumed the possibility to implement these gates without any practical constraints. This assumption was motivated by our analysis which revealed the key role played by the Kerr interaction in the performance of our ansatze (as it allows us to go beyond Gaussian transformations).

Note that the choice of the three types of gates presented here is not only motivated by their capacity to preserve a same number of bosons $N_B$, but also by their physical action on a many-boson wavefunctions. Indeed, the role of a BS gate is essentially to explore the space of many-boson configurations for two connected modes. On another note, the addition of Kerr (or rotation) gates in a circuit after a BS gate makes it possible to non-linearly affect the phase of each many-bosons configurations of the wavefunction (see Appendix B for an illustration of this with the minimal ansatz). These behaviors were the major motivations for the design of the ansatze introduced in the next section, that can therefore be regarded as physically-motivated as well as hardware-efficient.

### 3.2.2 Minimal beam-splitter-Kerr ansatz

Let us now introduce the first ansatz developed to encode the ground state of the attractive BH model: the so-called minimal beam-splitter-Kerr (BS-Kerr) ansatz. This circuit architecture was designed to produce a high fidelity encoding with as few quantum gates and parameters to optimize as possible, based on numerical investigations. Considering these constraints, 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

<figure><img src="image_6.png" /><figcaption>Figure 5: **BS-Kerr ansatz.** Quantum circuit structure 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 to each modes. The total gate count is $N_L(2N_S - 1)$ where $N_L$ is the number of layers.</figcaption></figure>

shown). Discarding these angles allows to reduce the complexity of the classical optimization with no impact on the final accuracy.

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 staircase structure to ensure that each mode can have the opportunity to exchange photons (via the beam-splitters) with all the other modes in the circuit.

Note that removing the Kerr gates from the ansatz leads to very poor results (not shown), showing that a passive Gaussian transformation is insufficient to encode the ground state of the attractive BH model despite we considered non-Gaussian inputs.

### 3.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 advantageous for numerical simulations. However, in experimental realizations one may want to minimize the number of Kerr layers. Therefore, we considered an alternative ansatz com-
<figure><img src="image_7.png" /><figcaption>Figure 6: Interferometer-Kerr ansatz. Quantum circuit architecture (a single layer) for $N_S = 3$ (top) and $N_S = 4$ (bottom) modes.</figcaption></figure>

posed 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 expressibility brought about by the complete 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, which are however much simpler and cheaper. This hypothesis is numerically investigated in Sec. 4.2.

The complexity of the ansatz can be reduced by getting rid of the rotation gates (physically irrelevant for most applications [112]) 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.

# 4 Numerical investigation of the ansatze capabilities

In this section, we present a numerical investigation of the expressibility (and resource efficiency) for 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 [113]. For the optimization of the quantum circuits parameters, the L-BFGS-B method was used as encoded in the optimParallel package [114]. The circuit's parameters were always initialized with a small random amplitude drawn in [-0.05, +0.05], and the maximum number of iterations was set to 20000. An ED code has been implemented in Python to compute the exact ground state by exact diagonalization. Further we also provide a sample code with the implementations of all the components of our method in [115] for interested readers.

## 4.1 Capacity of the BS-Kerr ansatz

To demonstrate and characterize 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 $|\Psi_0\rangle$,

$$ \mathcal{F} = |\langle\Psi_0|\hat{U}(\vec{\theta})|\Psi_{ini}\rangle|^2. \quad (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

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

which is naturally lower-bounded by zero. Note that no mention has been made so far about the nature of the initial state $|\Psi_{ini}\rangle$, for which we considered two choices: a single-mode and a two-mode states, detailed later in section 4.1.2.

### 4.1.1 Illustration of the infidelity optimisation

As an introduction of our study, in this first subsection we illustrate the evolution of the variational minimization of the infidelity by the mean of representative single runs of optimization. We gather in Fig. 7 several examples of optimization runs showing the minimization of the measure $\mathcal{I}$ as a function of the number of function evaluations called in the L-BFGS-B algorithm (as provided by the optimParallel package [114]), for dif-
<figure><img src="image_8.png" /><figcaption>Figure 7: Illustration of infidelity optimization runs realized to study the expressibility of the BS-kerr ansatz. Each plot shows the evolution of the infidelity $\mathcal{I}$ as a function of the number of function evaluations. These results were obtained for $N_L = 5$ layers circuits and initial single-mode states (i.e. $|\Psi_{ini}\rangle = |N_B, 0\rangle$, $|N_B, 0, 0\rangle$ and $|N_B, 0, 0, 0\rangle$ for the dimer, the three-, and the four-site systems, respectively). Left plots show results for $N_B = 2$, whereas right plots for $N_B = 4$. Each row (from top to bottom) shows results for the dimer, the three-, and the four-site BH models, respectively. In each case, four values of $\Lambda$ have been considered to span the different regimes of many-body correlations (as marked in figures 1 and 2). The infidelity limit $\mathcal{I} = 1\%$ considered as an indicator of the good ansatz expressibility is marked by a horizontal black dashed line. Note that all the cases treated here reach this limit except one which is stressed by a red arrow (on the bottom right panel).</figcaption></figure>

ferent values of $\Lambda$ to span the different correlation regimes. Each row of Fig. 7 (from top to bottom) shows results for the dimer, the three and four-site BH model respectively, for $N_B = 2$ (left panels) and $N_B = 4$ (right panels) bosons. $N_L = 5$ layers of the BS-Kerr ansatze were considered for all optimization runs, whatever the size of the system. As a remark, note that the optimization runs presented here are essentially illustrative. The shape of the optimization curves may actually differ from a run to another as small initial random values of circuit parameters are considered in practice. We checked the effects of these fluctuations by realizing multiple optimizations with random initial parameters and we noticed negligible changes in the final results obtained, hence the tendencies shown in Fig. 7 can be regarded as sufficiently general.

As readily seen in Fig. 7, every optimization curve crosses the limit $\mathcal{I} = 1\%$ (defined as the marker of good expressibility for the ansatze) for all values of $\Lambda$ but $\Lambda = 10$ on the last panel. The left panels ($N_B = 2$) manifest very low converged infidelities, and thus a very accurate encoding of the ground state of the system. The number of function evaluations required to fulfill the condition $\mathcal{I} = 1\%$ (and to reach final convergence) ranges from around 10 to 100 in

these cases. Conversely, increasing the number of bosons in the systems leads to an increase of the number of function evaluations, as illustrated on the right panels ($N_B = 4$). In this case, the number of evaluations required to satisfy the condition $\mathcal{I} = 1\%$ ranges from around 20 to 300. To reach final convergence (i.e. when the selected optimizer cannot minimize the cost function anymore) the number of evaluations also increases, ranging from 50 to 2000 approximately (all systems considered). The number of function evaluations increases with the number of bosons $N_B$ and sites $N_S$, i.e. with the size of the Fock space $D_{\hat{H}}$ (as shown in table 1) that is related to the complexity of the encoding problem. As a remark, note that in practice the final number of function evaluations $N_{eval}^{fin}$ differs from the final number of iterations $N_{iter}^{fin}$ realized by the LBFGS-B algorithm, as one iteration of the algorithm requires multiple function evaluations. We choose an upper limit of 20000 iterations to ensure convergence of the algorithm whatever the case considered. And as illustrated in Fig. 7, this limit is never reached as $N_{iter}^{fin} \le N_{eval}^{fin} \ll 20000$ (a feature that we also checked in many other test simulations). In Fig. 7, we see that a circuit with a depth of $N_L = 5$ layers of the BS-Kerr ansatz is sufficient to ensure an accurate encoding of all ex-
amples considered, except for the strongly inter-
acting case of the largest system $N_S = 4$, $N_B = 4$
and $\Lambda = 10$. In this case, the converged infidelity
is $\mathcal{I} \approx 7\%$ (i.e. $\mathcal{F} \approx 93\%$). A much better fidelity
can be reached by adding more layers to the quan-
tum circuit, thus enhancing the expressibility of
the ansatz. Typically, changing $N_L$ from 5 to 8
(not shown here), makes all the examples below
the limit $\mathcal{I} = 1\%$.

In spite of these promising results, one can still
suspect the presence of barren plateaus that ap-
pear for large systems in VQAs when random
initial parameters are used. For qubit-based de-
vices, barren plateaus denote the exponentially
vanishing cost function gradients with the num-
ber of qubits [116–121]. Note that the exact
conditions for the appearance of barren plateaus
on qubit-devices remains unclear, as several re-
cent works have shown that they can also arise
when physically-inspired [122, 123] and beyond
2-design ansatze [124] are considered. It is an
open question of how and when barren plateaus
may appear in CV ansatze, although recent stud-
ies indicate that CV systems may exhibit barren
plateaus where the cost gradients vanish expo-
nentially with the number of system modes (in-
stead of the number of qubits) [125, 126]. In
this work, we did not see any optimization prob-
lem that could indicate the presence of such phe-
nomenon, but it remains unclear for larger sys-
tems. Fortunately, due to the analogy between
the barren plateaus phenomenon on discrete and
continuous variable systems [125], one can ex-
pect that strategies used to bypass the barren
plateaus on qubit-based devices [122,127,128] are
also transferable to CV systems. Such an in-
volved study is out of the scope of the manuscript
and is left for future work.

In the next sections, we will present the results
of our numerical study of the expressibility (and
the resource efficiency) of the ansatze as a func-
tion of the number of bosons $N_B$ and sites $N_S$.
Our approach was inspired by Ref. [12]. Namely,
we realize a series of simulations where we pro-
gressively increase the number of layers $N_L$ in the
quantum circuits and perform each time one op-
timization run of infidelity. The choice of single
runs was motivated by many numerical simula-
tions which shown that very similar results are
always obtained in terms of encoding capacity
and resource efficiency. In practice, the more the

number of layers $N_L$ increases, the more the ex-
pressibility of the ansatz improves. We continue
to increase the number of layers $N_L$ until the in-
fidelity of the ground state fulfill the condition
$\mathcal{I} \le 1\%$ (i.e. $\mathcal{F} \ge 99\%$).

4.1.2 BH model with $N_S = 2$ sites

Focusing first on the BS-Kerr ansatz for the BH
dimer, the minimum number of layers detected to
reach $\mathcal{F} \ge 99\%$ for $N_B = 1$ to $N_B = 8$ is shown
in Fig. 8 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 single-mode one given
by $|\Psi_{ini}\rangle = |N_B, 0\rangle$ and a two-mode one given
by $|\Psi_{ini}\rangle = \left|\frac{N_B}{2}, \frac{N_B}{2}\right\rangle$ (n.b. for odd total number
of bosons we took $|\Psi_{ini}\rangle = \left|\frac{N_B+1}{2}, \frac{N_B+1}{2}-1\right\rangle$ in-
stead). As a first general remark, Fig. 8 shows
that, with only a few layers of circuit, a fidelity
higher than 99% is always reached whatever the
case considered (i.e. the correlation strength $\Lambda$,
the number of bosons $N_B$ and the choice of
initial state). These numerical results give a
first demonstration that short depth CV circuits
can accurately encode many-boson wavefunctions
even in the strongly correlated regime, at least in
the case of system sizes we have considered. Note
that to ensure an accurate encoding of the wave-
function, 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{H}}$, which is linear in $N_B$
for the BH dimer (as shown in Tab. 1). Interest-
estingly, the trend of the curves ensuring an effi-
cient encoding is globally sub-linear in $N_B what-
ever the initial states considered. Roughly, we
see here that choosing a number of layers such as
$N_L \approx 0.7 \times N_B$ seems to guarantee a very high
fidelity in all cases considered herein. The maxi-
mum number of Kerr gates involved here reaches
a maximum of $N_{Kerr} = 10$ (as this number scales
like $N_{Kerr} = N_S N_L$).

On another note, an interesting fact arises in
the weakly correlated case $\Lambda = 0.01$ when the
single-mode initial state is used (dark blue curve
in the upper panel). In this case, a single layer
of the BS-Kerr ansatz is sufficient to reach a fi-
delity higher than 99%, even when the number
of bosons increases. This property is true only
for the single-mode initial state. We rational-
<figure><img src="image_9.png" /><figcaption>Figure 8: **Ansatz encoding capacity (dimer).** Minimum depth required (in terms of number of layers $N_L$ and gate count) to generate a fidelity $F \ge 99\%$ [see Eq. (14)] for different values of $\Lambda$ and $N_B$. **Upper panel:** results obtained with the single-mode initial state. **Lower panel:** results obtained with the two-mode initial state.</figcaption></figure>

ize this result in Appendix A, where we show that, starting from such a state, an exact encoding ($F = 100\%$) of the BH ground state in the non-interacting (superfluid) regime ($\Lambda \to 0$) can be realized using only beam-splitters (this demonstration is realized for the three networks considered in our study). Note that even if no Kerr gates are required in this particular limit, an exact classical simulation of Fock states propagating in a passive interferometer has exponential complexity, as shown by seminal papers on the computational complexity of linear optical quantum computing [129, 130]. Indeed, even if the interferometer itself can be fully described by a $2N$-dimensional orthogonal matrix, (which is tractable), the action of an interferometer on an arbitrary input state requires a Fock space description of the gate (as opposed to a Gaussian input state, not considered here).

As an addition, in Fig. 9 we illustrate the evolution of the structural properties of a trial state (built with the BS-Kerr ansatz) as a function of $\Lambda$. In this case, we focus on a system with $N_B = 8$ bosons and we use a quantum circuit of $N_L = 6$ layers with $|\Psi_{ini}\rangle = |4,4\rangle$ as an initial state. The results obtained are compared to the exact ones generated via ED. On the upper panel of Fig. 9, we show the evolution of the IPR and the entanglement entropy $S$ of the trial state (red dots) compared to exact results (black curve) as a function of $\Lambda$. Note that each of the trial states obtained in the simulations have a fidelity $F \ge 99\%$. The lower panel shows the evolution of

<figure><img src="image_10.png" /><figcaption>Figure 9: **Trial-state properties (dimer).** Illustration of the properties of the quantum state obtained by minimization of the infidelity, using the BS-Kerr ansatz with $N_L = 6$, $N_B = 8$, $|\Psi_{ini}\rangle = |4,4\rangle$. **Upper panel:** IPR and entanglement entropy $S$ of the resulting ansatz state (red dots) compared to the exact ground state (black curves) as a function of $\Lambda$. **Lower panel:** Probability of occupation of the ansatz state in the Fock basis (green bars) compared with exact diagonalization (black bars). The $\Lambda$ values chosen for the four configurations are marked in the upper panel with red dashed vertical lines (same as in Fig. 1).</figcaption></figure>

the probability of occupation in the Fock basis for the trial state (green bars) compared to the exact ground state (black bars) for four representative values of interaction parameter $\Lambda = 0.01, 3, 5$ and $10$. The results presented here highlight 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$ for a maximal IPR value) to finally reach the famous “cat-like” superposition (with $\Lambda = 5$ to $10$ and the IPR tending to $2$). For all values of $\Lambda$ considered, the IPR, the entanglement entropy $S$ and the probability of occupation of the trial states are in very good agreement with the exact results, thus showing that the expressibility of our ansatz is strong enough to reproduce the exotic features of the exact ground state in different correlation regimes. Similar illustrations of the BH ground state restructuring will be given in the following of the paper for the three and four-site BH system.
<figure><img src="image_11.png" /><figcaption>Figure 10: Ansatz encoding capacity (3-site chain).
Minimum depth required (in terms of number of layers $N_L$ and gate count) to generate a fidelity $F \ge 99\%$ [see Eq. (14)] for different values of $\Lambda$ and $N_B$. **Upper panel:** results obtained with a single-mode initial state. **Lower panel:** results obtained with a two-mode initial state.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: Ansatz encoding capacity (4-site chain).
Minimum depth required (in terms of number of layers $N_L$ and gate count) to generate a fidelity $F \ge 99\%$ [see Eq. (14)] for different values of $\Lambda$ and $N_B$. **Upper panel:** results obtained with a single-mode initial state. **Lower panel:** results obtained with a two-mode initial state.</figcaption></figure>

### 4.1.3 BH model with $N_S = 3$ and 4 sites

Turning to larger systems, Figs. 10 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 single-mode and two-mode initial states are again considered and read $|\Psi_{ini}\rangle = |N_B, 0, 0\rangle$ and $|\Psi_{ini}\rangle = |\frac{N_B}{2}, 0, \frac{N_B}{2}\rangle$ respectively (n.b. for odd total number of bosons we took $|\Psi_{ini}\rangle = |\frac{N_B+1}{2}, 0, \frac{N_B+1}{2}-1\rangle$ 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{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. 10. In this case, we retrieve some differences between the results obtained with the two initial states. The single-mode initial state seems to perform better in term of resource efficiency as the latter requires globally less layers to encode the ground state with a fidelity higher than 99%. Typically, we see in Fig. 10 that the number of layers detected increases up to $N_L = 8$ for a single-mode initial state against a maximum of $N_L = 10$ layers for the two-mode one. The associated number of Kerr gate being 24 and 30 respectively. Still, the results obtained here reveal the good expressibility of the minimal ansatz which is able to accurately transform both different initial

states into very good approximations of the exact ground state. As readily seen in Fig. 10, choosing a circuit depth with $N_L \sim 1.25 \times N_B$ guarantees an efficient encoding of the ground state in all cases considered here. As a last remark for the three-site system, note again that only one layer is required for $\Lambda = 0.01$ when this initial state is used (dark blue curve). We refer the interested reader to Appendix A where we give a mathematical demonstration of this feature.

The complexity grows even stronger for the four-site BH network. We follow the exact same strategy than for the previous models by considering $|\Psi_{ini}\rangle = |N_B, 0, 0, 0\rangle$ as the single-mode initial state and $|\Psi_{ini}\rangle = |\frac{N_B}{2}, 0, \frac{N_B}{2}, 0\rangle$ as the two-mode one (n.b. for odd total number of bosons we took $|\frac{N_B+1}{2}, 0, \frac{N_B+1}{2}-1, 0\rangle$ 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{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. 11. 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$. Here again we
<figure><img src="image_13.png" /></figure>

Figure 12: **Trial-state properties (3-site chain).** Illustration of the properties of the quantum state obtained by minimization of the infidelity, using the BS-Kerr ansatz with $N_L = 6$, $N_B = 4$, and $|\Psi_{ini}\rangle = |2,0,2\rangle$. **Upper panel:** IPR and entanglement entropy *S* of the resulting ansatz state (red dots) compared to the exact ground state (black curves) as a function of Λ. **Lower panel:** Probability of occupation of the ansatz state in the Fock basis (green bars) compared with exact diagonalization (black bars). The Λ values chosen for the four configurations are marked in the upper panel with red dashed vertical lines (same as in Fig. 2).

see that the single-mode initial state seems to perform better compared to the two-mode case. The associated maximum number of layers being $N_L = 8$ for the latter (with 32 Kerr gates) compared to $N_L = 10$ for a two-mode initial state (with 40 Kerr gates). We also retrieve on Fig. 11 the case of an exact encoding with one layer for the initial single-mode state for weakly correlated regime $\Lambda = 0.01$.

As an addition, Figs. 12 and 13 illustrate the structural transition of the BH ground state for the three and four-site BH networks, respectively. Similarly to Fig. 9, the upper part of the plots show the evolution of the IPR and the entanglement entropy of trial states (red dots) compared to exact ground state (black curves) as a function of Λ. Note that, all the trial states generated here present a fidelity of $F \ge 99\%$. The lower panels show the density of probability of the trial states (green bar) compared to the exact ground state (black bars). Just as in the BH dimer, the IPR,

<figure><img src="image_14.png" /></figure>

Figure 13: **Trial-state properties (4-site chain).** Illustration of the properties of the quantum state obtained by minimization of the infidelity, using the BS-Kerr ansatz with $N_L = 6$, $N_B = 3$ and $|\Psi_{ini}\rangle = |2,0,1,0\rangle$. **Upper panel:** IPR and entanglement entropy *S* of the resulting ansatz state (red dots) compared to the exact ground state (black curves) as a function of Λ. **Lower panel:** Probability of occupation of the ansatz state in the Fock basis (green bars) compared with exact diagonalization (black bars). The Λ values chosen for the four configurations are marked in the upper panel with red dashed vertical lines (same as in Fig. 2).

the entanglement entropy *S* and the density of probability of the resulting trial states are in very good agreement with the exact ground state.

As a summary of the results presented so far, a remarkable property is that very high fidelity encoding can be performed with the minimal ansatz for the tree system, whatever the regimes of many-body interactions. Globally, our results suggest that, for a fixed system size, the number of layers required to ensure an encoding of $F \ge 99\%$ scales linearly with the number of bosons. Furthermore, our simulations also demonstrate that high fidelity may be obtained even if we consider different initial states $|\Psi_{ini}\rangle$ that are simple Fock states (and not complex quantum superpositions of many Fock states). Hence, the minimal ansatz is expressive enough to accurately transform such single Fock states into complex linear combinations of $D_{\hat{H}}$ configurations encoding strong quantum correlations.
<figure><img src="image_15.png" /><figcaption>Figure 14: Scaling of the ansatz encoding capacity with the system's size $N_S$. Minimum depth required (in terms of number of layers $N_L$) to generate a fidelity $F \ge 99\%$ [see Eq. (14)] with the minimal BS-Kerr ansatz as a function of the number of sites (ranging from $N_S = 2$ to 6). Here, three different values of many-body interaction ($\Lambda = 0.01, 5$ and $10$) are considered to span the different regimes of correlation. The three panels, from left to right, show results obtained with a maximum number of $N_B = 2, 3$ and 4 bosons, respectively. A single-mode initial state is employed in all simulations.</figcaption></figure>

### 4.1.4 Effects of the system size $N_S$

Let us now investigate how the increase in the system size can affect the expressibility and the resource efficiency of the minimal BS-Kerr ansatz. In analogy to the previous study, we evaluate the number of layers $N_L$ required to satisfy $F = 99\%$, but as a function of the total number of sites $N_S$ rather than the number of bosons $N_B$. The results are presented in Fig. 14 for three values of interaction parameter ($\Lambda = 0.01, 5$ and $10$) and a number of bosons $N_B = 2, 3$ and $4$. As readily seen in Fig. 14, the number of layers $N_L$ required to ensure a high fidelity scales linearly with the number of sites of the BH system (whatever the number of bosons $N_B$ is). Focusing on the strong many-body interaction regime (i.e. $\Lambda = 10$), a scaling of $N_L \sim 0.5N_S$ is obtained for $N_B = 2$ bosons, $N_L \sim 1.5N_S$ for $N_B = 3$ bosons, and $N_L \sim 2N_S$ for $N_B = 4$ bosons. The prefactor in the linear scaling seems to slightly increase with the total number of bosons $N_B$ (a trend already observed Figs. 9, 12 and 13). Note that such linear behaviors strongly contrast with the non-linear increase of the number of configurations (i.e. the size $D_{\hat{H}}$ of the Fock space) necessary to represent the exact ground-state of the BH model. Indeed, as an illustration, note that in the case of $N_B = 4$ bosons the number of configurations involved starts from $D_{\hat{H}} = 5$ for the dimer ($N_S = 2$) to reach $D_{\hat{H}} = 126$ for the BH chain with $N_S = 6$ sites (which represents the maximum number of Fock states generated in our study). On another note about the weakly correlated regime (dark blue curves in Fig. 14), we see here again that a single layer is sufficient to

ensure a high fidelity encoding of the BH ground state, for all the number of sites $N_S$ and number of bosons $N_B$ considered (see Appendix A for a mathematical proof). Finally, to summarize our observations realized so far, our results from the previous sections combined with the ones presented here reveal that the number of layer $N_L$ required to accurately encode the ground state of the BH model presents a linear trend with respect to both $N_B$ and $N_S$. This being verified at least for the small to medium-sized BH systems considered herein.

## 4.2 Capacity of the interferometer-Kerr ansatz

Let us now turn to the interferometer-Kerr ansatz described in Sec. 3.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. 15, 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 (see Fig. 10). 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 num-
<figure><img src="image_16.png" /></figure>

**Figure 15:** Ansatz encoding capacity (3-site chain). Minimum depth required (in terms of number of layers $N_L$ and gate count) to generate a fidelity $F \ge 99\%$ [see Eq. (14)] for different values of $\Lambda$ and $N_B$, using the two-mode 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.

<figure><img src="image_17.png" /></figure>

**Figure 16:** Ansatz encoding capacity (4-site chain). Minimum depth required (in terms of number of layers $N_L$ and gate count) to generate a fidelity $F \ge 99\%$ [see Eq. (14)] for different values of $\Lambda$ and $N_B$, using the two-mode 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.

umber of modes as $N_{\text{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. 15. 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 guarantees 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. [112], 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. 15, 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-splitter gates (third panel) gives similar results than 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$. However, for larger systems, there are more additional phase angles than rotation gates as the interferometer is composed of $N_S(N_S - 1)/2$ beam-splitters but only $N_S$ rotation gates, and getting rid of one or the other might lead to different results. As there are more beam-splitters than rotation gates, we expect that getting rid of the rotation gates is safe, especially as they are supposed to be physically irrelevant for most applications [112].
Turning to the four-site model in Fig. 16, 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.

# 5 VQE simulations

In the previous section, we gave numerical evidence on the capacity of our proposed ansatz 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}) = \langle \Psi_{ini} | \hat{U}(\vec{\theta})^{\dagger} \hat{H} \hat{U}(\vec{\theta}) | \Psi_{ini} \rangle \ge E_0. \quad (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 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.

## 5.1 Ideal VQE simulations for the three- and four-site BH model

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 $|\Psi(\vec{\theta}) = \hat{U}(\vec{\theta})|\Psi_{ini}\rangle$ at the end of the circuit, which makes it possible to exactly estimate the energy $E(\vec{\theta})$ during the optimization process.

Figure 17 shows the probability of occupation in the Fock basis of the ground states obtained

<figure><img src="image_18.png" /><figcaption>Figure 17: Trial-state properties (3-/4-site chains). Probability of occupation in the Fock basis obtained by ED (black bars) and VQE simulations, using 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 was used. Upper panel: Ground state of the three-site model with $N_B = 4$, $N_L = 6$ and $|\Psi_{ini}\rangle = |2,0,2\rangle$. Lower panel: Ground state of the four-site model with $N_B = 3$, $N_L = 6$ and $|\Psi_{ini}\rangle = |2,0,1,0\rangle$.</figcaption></figure>

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 two-mode states ($|\Psi_{ini}\rangle = |2,0,2\rangle$ and $|\Psi_{ini}\rangle = |2,0,1,0\rangle$) are considered. The depth of the circuits imposes to optimize 30 (respectively 42) parameters for the three-site (four-site) network.

As shown in Fig. 17, the probability of occupation follows the same trend as the previous ones in Figs. 12 and 13, where the infidelity $I$ was minimized instead of the ground-state energy. In all cases of Fig. 17, the resulting states have fidelity $F \ge 99\%$ and an energy error $\Delta E = E(\vec{\theta}^*) - E_0 \le 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 single-mode initial states
(i.e. $|\Psi_{ini}\rangle = |4,0,0\rangle$ and $|\Psi_{ini}\rangle = |3,0,0,0\rangle$, not shown), and led to similar accuracy. These results give a good demonstration of the ability of the BS-Kerr ansatz to map different initial states to the exact ground state of the attractive BH model (in all correlation regimes considered herein) in a more practical VQE context.

## 5.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 experiment.

### 5.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 method which can be realized with photon number resolving (PNR) detectors [131, 132]. Note that the development of PNR (and other similar techniques) is critical for the implementation of scalable quantum information processing (see the review [133] and references therein about photon detection methods for photonic quantum computing). These methods are still actively investigated to push the limit of maximum number of detectable photons, with a current maximum of around several tenths of photons that are resolved simultaneously.

In the following, the notation $\langle \hat{O} \rangle_{\Psi} \equiv \langle \Psi | \hat{O} | \Psi \rangle$ refers to the quantum average of an operator $\hat{O}$ realized at the end of the circuit when the photons are in the final state $|\Psi\rangle$. Starting from the BH Hamiltonian $\hat{H}$ defined in Eq. (1), measuring the energy of the final photonic state $|\Psi\rangle$ of the circuit requires measuring two different contributions,

$$ \langle \hat{H}_{\text{Hoppings}} \rangle = -J \sum_{\langle p,q \rangle}^{N_S} \langle (a_p^\dagger a_q + a_q^\dagger a_p) \rangle_{\Psi}, \quad (17) $$

and

$$ \langle \hat{H}_{\text{Loc-Int}} \rangle = -\frac{U}{2} \sum_{p}^{N_S} \langle n_p(n_p - 1) \rangle_{\Psi}. \quad (18) $$

The one-body (hopping) term in Eq. (17) can be measured by introducing additional 50/50 beam-

splitter gates (i.e. with $\phi = 0$ and $\theta = \pi/4$) between each pair of modes before measuring, thus leading to

$$ B_{p,q}(\pi/4, 0) [a_p^\dagger a_q + a_q^\dagger a_p] B_{p,q}(\pi/4, 0)^\dagger = n_q - n_p \quad (19) $$

such that

$$ \langle a_p^\dagger a_q + a_q^\dagger a_p \rangle_{\Psi} = \langle n_q - n_p \rangle_{\tilde{\Psi}} \quad (20) $$

where $|\tilde{\Psi}\rangle = B_{p,q}(\pi/4, 0)|\Psi\rangle$. Thus, evaluating the hopping term between two modes with respect to the original trial state $|\Psi\rangle$ 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}, \quad (21) $$

where $\text{Var}(n_p) = \langle n_p^2 \rangle_{\Psi} - \langle n_p \rangle_{\Psi}^2$ 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, an extended version of the BH model – originally designed to study the supersolid phase of helium [94, 95] – is sometimes considered (see for example Refs. [91, 134–136]) and reads

$$ \hat{H}' = \hat{H} + \sum_{p} \mu_{p} n_{p} + \sum_{p,q} V_{pq} n_{p} n_{q}, \quad (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 two 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 chemical potential

$$ \langle \hat{H}_{\text{Chem-Pot}} \rangle_{\Psi} = \sum_{p} \mu_{p} \langle n_{p} \rangle_{\Psi} \quad (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{H}_{\text{dip-dip}} \rangle_{\Psi} = \sum_{p,q} V_{pq} \langle n_p n_q \rangle \quad (24) $$
can be found by evaluating the covariance matrix of the photon number distribution in all the modes:

$$ \langle n_p n_q \rangle = \mathrm{Cov}(n_p, n_q) + \langle n_p \rangle \langle n_q \rangle. \quad (25) $$

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

Note that the extended Bose-Hubbard model has been shown to capture essential properties of many systems [137], including ultracold atoms and molecules in optical lattices [87, 88, 90, 91, 138, 139], Josephson junction arrays [140, 141], and narrow-band superconductors [142]. Testing our proposed ansatz on this model is a natural follow-up of this study and is left for future work.

### 5.2.2 VQE with and without sampling noise

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. 18 shows the ground-state fidelity $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. 5.2.1. To mitigate the impact of noise in the classical optimization, we employed a “noise-proof” optimizer, namely CMA-ES [143]. 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. 18). In analogy with the previous study in Sec. 4, 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. 18) 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 \le 5$, the results obtained with sampling noise are globally comparable to the ones from an ideal VQE. Interestingly, when the single-mode 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 single-mode initial state is used. Over this limit, the fidelity of the VQE trial state gets stuck at $F \sim 50\%$. This occurs when the many-body interaction strongly intensifies, which leads to the ground state and first excited state to be quasi degenerate. This phenomenon arises as the many-body attraction essentially favors configurations that gather as many bosons as possible on a same site (as illustrated in Fig. 9, the most favorable configurations are $|N_B, 0\rangle$ and $|0, N_B\rangle$, followed by $|N_B - 1, 1\rangle$ and $|1, N_B - 1\rangle$, and then $|N_B - 2, 2\rangle$ and $|2, N_B - 2\rangle$, etc.). If we simplify by considering the most favorable configurations $|N_B, 0\rangle$ and $|0, N_B\rangle$, the resulting quasi-degenerate ground and first excited state can roughly be seen as a “plus version” of cat-state $|\Psi_0\rangle \propto |N_B, 0\rangle + |0, N_B\rangle$ and a “minus version” $|\Psi_1\rangle \propto |N_B, 0\rangle - |0, N_B\rangle$, respectively. In this case, running a noisy VQE simulation with an initial single-mode state $|\Psi_{ini}\rangle = |N_B, 0\rangle$ produces a final energy that stands in between both quasi degenerate states’ energies (not shown here). As the sampling noise typically equates the amplitude of the spectral gap between both states, the minimization gets stuck here. The final fidelity of the trial states is then of 50% with both quasi degenerate states, indicating then that $|\Psi(\theta)\rangle$ turned into a linear combination of both of them. Note that a solution to converge to the exact ground-state energy is to increase the number of measures to reduce the fluctuations due to the sampling
<figure><img src="image_19.png" /><figcaption>Figure 18: **Ground-state fidelity F of trial states obtained with VQE simulations for the BH dimer (using the BS-Kerr ansatz).** The final fidelity is shown for four different values of $\Lambda$ and as a function of the number of bosons $N_B$ and the number of layers $N_L$ in the circuit. **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 [143] 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 single-mode and a two-mode initial state, respectively.</figcaption></figure>

noise and resolve the small spectral gap.

Therefore, starting from an initial single-mode state, cat-like states seem difficult to realize by the BS-Kerr ansatz when sampling noise is present. However, much better accuracy is reached when the two-mode 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 \ge 5$ here). In summary, it seems from Fig. 18 that – within the BS-Kerr ansatz – a weakly-correlated state is more easily encoded when starting from a single-mode state rather than a two-mode one, and vice versa for a strongly-correlated state. Note that better results could possibly be reached by using other noise-adapted optimizers rather than CMA-ES used in this work, see Refs. [144, 145] and references therein. If the optimization becomes a problem when tackling larger systems, one can decompose it into smaller pieces that are iteratively solved by small sized circuits, as proposed by Fujii et al. with the so-called deep VQE based on the divide and conquer approach [13].

## 6 Conclusions

In this work, we tackled the problem of encoding strongly correlated many-boson ground-state wavefunctions on a photonic quantum device. A particular attention was paid to the expressibil-

ity and the resource efficiency of the circuits for different regimes of many-body interaction. For this, 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 sam-
pling 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.

As further suggestions for future work, one could include noise and losses in realistic simulations of actual photonic circuitry, as well as time-frequency effects which affect non-linear optical elements such as squeezers and, presumably, Kerr media. Moreover, although the ansatz proposed in our work were designed to solve the BH model, they could in principle be used for a larger range of applications, such as the absorption of Helium-4 on graphite [146,147].

## Acknowledgments

SY and BS sincerely thank Eleanor Scerri, Xavi Bonet-Monroig and Vincent Pouthier for fruitful discussions. SY and BS acknowledge support from the Netherlands Organization for Scientific Research (NWO/OCW). SY also acknowledges the Interdisciplinary Thematic Institute ITI-CSC via the IdEx Unistra (ANR-10-IDEX-0002) within the program Investissement d'Avenir. VD acknowledges support by the Dutch Research Council (NWO/OCW), as part of the Quantum Software Consortium program (project number 024.003.037).

## 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 \to 0$. Note that this demonstration works only for initial single-mode states such as $|N_B, 0\rangle$, $|N_B, 0, 0\rangle$ and $|N_B, 0, 0, 0\rangle$ for the dimer, three- and four-site BH model, respectively.

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

$$|\Psi_0\rangle = \frac{1}{\sqrt{2N_B!}} [a_1^\dagger + a_2^\dagger]^{N_B} |0,0\rangle. \quad (26)$$

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

$$B_{1,2}(\theta, 0) |N_B, 0\rangle = \frac{1}{\sqrt{N_B!}} [\cos(\theta)a_1^\dagger + \sin(\theta)a_2^\dagger]^{N_B} |0, 0\rangle \quad (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

$$|\Psi_0\rangle = \frac{1}{\sqrt{3N_B!}} [a_1^\dagger + a_2^\dagger + a_3^\dagger]^{N_B} |0, 0, 0\rangle. \quad (28)$$

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

$$B_{2,3}(\theta', 0) B_{1,2}(\theta, 0) |N_B, 0, 0\rangle = \frac{1}{\sqrt{N_B!}} \left[ \begin{aligned} & \cos(\theta)a_1^\dagger + \sin(\theta)\cos(\theta')a_2^\dagger \\ & + \sin(\theta)\sin(\theta')a_3^\dagger \end{aligned} \right]^{N_B} |0, 0, 0\rangle, \quad (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' = \pi/4$. Finally, the exact ground state of the four-site model reads

$$|\Psi_0\rangle = \frac{1}{\sqrt{4N_B!}} [a_1^\dagger + a_2^\dagger + a_3^\dagger + a_4^\dagger]^{N_B} |0, 0, 0, 0\rangle. \quad (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 single-mode initial state, thus leading to

$$B_{3,4}(\theta'', 0) B_{2,3}(\theta', 0) B_{1,2}(\theta, 0) |N_B, 0, 0, 0\rangle = \frac{1}{\sqrt{N_B!}} \left[ \begin{aligned} & \cos(\theta)a_1^\dagger + \sin(\theta)\cos(\theta')a_2^\dagger \\ & + \sin(\theta)\sin(\theta')\cos(\theta'')a_3^\dagger \\ & + \sin(\theta)\sin(\theta')\sin(\theta'')a_4^\dagger \end{aligned} \right]^{N_B} |0, 0, 0, 0\rangle. \quad (31)$$

Fixing the beam-splitter's parameters to $\theta = \pi/3$, $\theta' = \arccos(1/\sqrt{3})$ and $\theta'' = \pi/4$ generates an exact representation of the exact ground state in Eq. (30).
B Action of the BS-Kerr ansatz on the
BH dimer with 2 bosons

In this section, we illustrate the functioning of the
BS-Kerr ansatz to solve the ground state of the
Bose-Hubbard model. We focus on the simplest
case possible: a dimer containing two bosons (i.e.
$N_S = N_B = 2$). As a reference, we consider the
exact form of the ground state for an intermediate
interaction regime $\Lambda = 6$, which is given by

$$
|\Psi_0\rangle = \frac{1}{\sqrt{5}} |11\rangle + \sqrt{\frac{2}{5}} (|20\rangle + |02\rangle). \quad (32)
$$

For the illustration, let us consider an initial ho-
mogeneous distribution of the bosons in the two
modes, i.e. an initial state $|11\rangle$. Starting from
this, we will demonstrate that a single layer of
the BS-Kerr ansatz is sufficient to exactly encode
the ground state of the system. The associated
state generated by the circuit is here

$$
|\Psi(\vec{\theta})\rangle = K_1(\theta_1)K_2(\theta_2)B(\theta)|11\rangle \quad (33)
$$

Focusing first on the BS-gate’s action, we can
show (using Eq. 11) that

$$
B(\theta)|11\rangle = \cos(2\theta)|11\rangle + \frac{\sin(2\theta)}{\sqrt{2}}(|02\rangle - |20\rangle). \quad (34)
$$

This result clearly illustrates the role of the BS
gate: it spreads the bosonic wavefunction over
the three different accessible Fock states of the
system. Applying Kerr gates on this quantum
superposition of Fock state introduces non-linear
phases on each independent Fock state such as

$$
\begin{align}
|\Psi(\vec{\theta})\rangle &= K_1(\theta_1)K_2(\theta_2)B(\theta)|11\rangle \nonumber \\
&= e^{i(\theta_1+\theta_2)} \cos(2\theta)|11\rangle \tag{35} \\
&\quad + \frac{\sin(2\theta)}{\sqrt{2}}(e^{i4\theta_2}|02\rangle - e^{i4\theta_1}|20\rangle). \nonumber
\end{align}
$$

In this case, encoding the exact system's ground
state with the three gates involved means solving
the set of coupled non-linear equations to match
the best the coefficients in 35 and 32. An exact
solution of the problem is given by $\theta_1 = 3\pi/8$,
$\theta_2 = \pi/8$ and $\theta = \arccos(1/\sqrt{5})/2$. Note that
replacing Kerr-gates with simple rotation-gates
makes the exact encoding of the ground state im-
possible.

References

[1] PERUZZO, Alberto ; McCLEAN, Jarrod ; SHADBOLT, Peter ; YUNG, Man-Hong ; ZHOU, Xiao-Qi ; LOVE, Peter J. ; ASPURU-GUZIK, Alán ; O'BRIEN, Jeremy L.: A variational eigenvalue solver on a photonic quantum processor. In: Nature Comm. 5 (2014), 4213. https://doi.org/10.1038/ncomms5213

[2] O'MALLEY, Peter J. ; BABBUSH, Ryan ; KIVLICHAN, Ian D. ; ROMERO, Jonathan ; McCLEAN, Jarrod R. ; BARENDS, Rami ; KELLY, Julian ; ROUSHAN, Pedram ; TRANTER, Andrew ; DING, Nan u. a.: Scalable quantum simulation of molecular energies. In: Phys. Rev. X 6 (2016), Nr. 3, 031007. https://doi.org/10.1103/PhysRevX.6.031007

[3] SHEN, Yangchao ; ZHANG, Xiang ; ZHANG, Shuaining ; ZHANG, Jing-Ning ; YUNG, Man-Hong ; KIM, Kihwan: Quantum implementation of the unitary coupled cluster for simulating molecular electronic structure. In: Phys. Rev. A 95 (2017), Nr. 2, 020501. https://doi.org/10.1103/PhysRevA.95.020501

[4] KANDALA, Abhinav ; MEZZACAPO, Antonio ; TEMME, Kristan ; TAKITA, Maika ; BRINK, Markus ; CHOW, Jerry M. ; GAMBETTA, Jay M.: Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets. In: Nature 549 (2017), Nr. 7671, 242–246. https://doi.org/10.1038/nature23879

[5] HEMPEL, Cornelius ; MAIER, Christine ; ROMERO, Jonathan ; McCLEAN, Jarrod ; MONZ, Thomas ; SHEN, Heng ; JURCEVIC, Petar ; LANYON, Ben P. ; LOVE, Peter ; BABBUSH, Ryan ; ASPURU-GUZIK, Alán ; BLATT, Rainer ; ROOS, Christian F.: Quantum chemistry calculations on a trapped-ion quantum simulator. In: Phys. Rev. X 8 (2018), Nr. 3, 031022. https://doi.org/10.1103/PhysRevX.8.031022

[6] COLLESS, James I. ; RAMASESH, Vinay V. ; DAHLEN, Dar ; BLOK, Machiel S. ; KIMCHI-SCHWARTZ, ME ; McCLEAN, JR ; CARTER, J ; DE JONG, WA ; SIDDIQI,
I: Computation of molecular spectra on a quantum processor with an error-resilient algorithm. In: *Phys. Rev. X* 8 (2018), Nr. 1, 011021. https://doi.org/10.1103/PhysRevX.8.011021

[7] MAZZOLA, Guglielmo; OLLITRAULT, Pauline J.; BARKOUTSOS, Panagiotis K.; TAVERNELLI, Ivano: Nonuni-tary operations for ground-state calcula-tions in near-term quantum computers. In: *Phys. Rev. Lett.* 123 (2019), Nr. 13, 130501. https://doi.org/10.1103/PhysRevLett.123.130501

[8] NAM, Yunseong; CHEN, Jwo-Sy; PISENTI, Neal C.; WRIGHT, Kenneth; DELANEY, Conor; MASLOV, Dmitri; BROWN, Kenneth R.; ALLEN, Stewart; AMINI, Jason M.; APISDORF, Joel u.a.: Ground-state energy estimation of the water molecule on a trapped-ion quantum computer. In: *npj Quantum Inf.* 6 (2020), Nr. 1, 1-6. https://doi.org/10.1038/s41534-020-0259-3

[9] O'BRIEN, Thomas E.; SENJEAN, Bruno; SAGASTIZABAL, Ramiro; BONET-MONROIG, Xavier; DUTKIEWICZ, Alicja; BUDA, Francesco; DICARLO, Leonardo; VISSCHER, Lucas: Calculating energy derivatives for quantum chemistry on a quantum computer. In: *npj Quantum Inf.* 5 (2019), Nr. 1, 1-12. https://doi.org/10.1038/s41534-019-0213-4

[10] SOKOLOV, Igor O.; BARKOUTSOS, Panagiotis K.; MOELLER, Lukas; SUCHS-LAND, Philippe; MAZZOLA, Guglielmo; TAVERNELLI, Ivano: Microcanonical and finite-temperature ab initio molecular dynamics simulations on quantum computers. In: *Phys. Rev. Research* 3 (2021), Nr. 1, 013125. https://doi.org/10.1103/PhysRevResearch.3.013125

[11] ARUTE, Frank; ARYA, Kunal; BABBUSH, Ryan; BACON, Dave; BARDIN, Joseph C.; BARENDS, Rami; BOIXO, Sergio; BROUGHTON, Michael; BUCKLEY, Bob B.; BUELL, David A.; BURKETT, Brian; BUSHNELL, Nicholas; CHEN, Yu; CHEN, Zijun; CHIARO, Benjamin; COLLINS, Roberto; COURTNEY, William; DEMURA, Sean; DUNSWORTH, Andrew; FARHI, Edward; FOWLER, Austin; FOXEN, Brooks; GIDNEY, Craig; GIUSTINA, Marissa; GRAFF, Rob; HABEGGER, Steve; HARRIGAN, Matthew P.; HO, Alan; HONG, Sabrina; HUANG, Trent; HUGGINS, William J.; IOFFE, Lev; ISAKOV, Sergei V.; JEFFREY, Evan; JIANG, Zhang; JONES, Cody; KAFRI, Dvir; KECHEDZHI, Kostyantyn; KELLY, Julian; KIM, Seon; KLIMOV, Paul V.; KOROTKOV, Alexander; KOSTRITSA, Fedor; LANDHUIS, David; LAPTEV, Pavel; LINDMARK, Mike; LUCERO, Erik; MARTIN, Orion; MARTINIS, John M.; MCCLEAN, Jarrod R.; McEWEN, Matt; MEGRANT, Anthony; MI, Xiao; MOHSENI, Masoud; MRUCZKIEWICZ, Wojciech; MUTUS, Josh; NAAMAN, Ofer; NEELEY, Matthew; NEILL, Charles; NEVEN, Hartmut; NIU, Murphy Y.; O'BRIEN, Thomas E.; OS-TBY, Eric; PETUKHOV, Andre; PUTTERMAN, Harald; QUINTANA, Chris; ROUSHAN, Pedram; RUBIN, Nicholas C.; SANK, Daniel; SATZINGER, Kevin J.; SMELYANSKIY, Vadim; STRAIN, Doug; SUNG, Kevin J.; SZALAY, Marco; TAKESHITA, Tyler Y.; VAINSENCHER, Amit; WHITE, Theodore; WIEBE, Nathan; YAO, Z. J.; YEH, Ping; ZALCMAN, Adam: Hartree-Fock on a superconducting qubit quantum computer. In: *Science* 369 (2020), Nr. 6507, 1084-1089. https://doi.org/10.1126/science.abb9811

[12] CADE, Chris; MINEH, Lana; MONTANARO, Ashley; STANISIC, Stasja: Strategies for solving the Fermi-Hubbard model on near-term quantum computers. In: *Phys. Rev. B* 102 (2020), Nr. 23, 235122. https://doi.org/10.1103/PhysRevB.102.235122

[13] FUJII, Keisuke; MITARAI, Kosuke; MIZUKAMI, Wataru; NAKAGAWA, Yuya O.: Deep Variational Quantum Eigensolver: a divide-and-conquer method for solving a larger problem with smaller size quantum computers. In: *arXiv:2007.10917* (2020). https://arxiv.org/abs/2007.10917

[14] XU, Luogen; LEE, Joseph T.; FREERICKS, JK: Test of the unitary coupled-cluster vari-
ational quantum eigensolver for a simple strongly correlated condensed-matter system. In: Mod. Phys. Lett. B 34 (2020), Nr. 19n20, 2040049. https://doi.org/10.1142/S0217984920400497

[15] OLLITRAULT, Pauline J. ; BAIARDI, Alberto ; REIHER, Markus ; TAVERNELLI, Ivano: Hardware efficient quantum algorithms for vibrational structure calculations. In: Chem. Sci. 11 (2020), Nr. 26, 6842-6855. https://doi.org/10.1039/DSC01908A

[16] ROMERO, Jonathan ; ASPURU-GUZIK, Alán: Variational quantum generators: Generative adversarial quantum machine learning for continuous distributions. In: Adv. Quantum Technol. 4 (2021), Nr. 1, 2000003. https://doi.org/10.1002/qute.202000003

[17] GRIMSLEY, Harper R. ; ECONOMOU, Sophia E. ; BARNES, Edwin ; MAYHALL, Nicholas J.: An adaptive variational algorithm for exact molecular simulations on a quantum computer. In: Nat. Commun. 10 (2019), Nr. 1, 1-9. https://doi.org/10.1038/s41467-019-10988-2

[18] TANG, Ho L. ; SHKOLNIKOV, V.O. ; BARRON, George S. ; GRIMSLEY, Harper R. ; MAYHALL, Nicholas J. ; BARNES, Edwin ; ECONOMOU, Sophia E.: Qubit-ADAPT-VQE: An Adaptive Algorithm for Constructing Hardware-Efficient Ansätze on a Quantum Processor. In: PRX Quantum 2 (2021), Apr, 020310. https://doi.org/10.1103/PRXQuantum.2.020310

[19] MONTANARO, Ashley ; STANISIC, Stasja: Compressed variational quantum eigensolver for the fermi-hubbard model. In: arXiv:2006.01179 (2020). https://arxiv.org/abs/2006.01179

[20] MASAKI-KATO, A ; YUNOKI, S ; HIRASHIMA, DS: Quantum Monte Carlo study of the superfluid density in quasi-one-dimensional systems of hard-core bosons: Effect of the suppression of phase slip-page. In: Phys. Rev. B 100 (2019), Nr. 22, 224515. https://doi.org/10.1103/PhysRevB.100.224515

[21] KRAUTH, Werner: Quantum Monte Carlo calculations for a large number of bosons

in a harmonic trap. In: Phys. Rev. Lett. 77 (1996), Nr. 18, 3695. https://doi.org/10.1103/PhysRevLett.77.3695

[22] DEL MAESTRO, Adrian ; AFFLECK, Ian: Interacting bosons in one dimension and the applicability of Luttinger-liquid theory as revealed by path-integral quantum Monte Carlo calculations. In: Phys. Rev. B 82 (2010), Nr. 6, 060515. https://doi.org/10.1103/PhysRevB.82.060515

[23] PURWANTO, Wirawan ; ZHANG, Sh Wei: Quantum Monte Carlo method for the ground state of many-boson systems. In: Phys. Rev. E 70 (2004), Nr. 5, 056702. https://doi.org/10.1103/PhysRevE.70.056702

[24] POLLET, Lode: A review of Monte Carlo simulations for the Bose-Hubbard model with diagonal disorder. In: C. R. Phys. 14 (2013), Nr. 8, 712-724. https://doi.org/10.1016/j.crhy.2013.08.005

[25] BOGNER, Benjamin ; DE DANILOFF, Clément ; RIEGER, Heiko: Variational Monte-Carlo study of the extended Bose-Hubbard model with short-and infinite-range interactions. In: Eur. Phys. J. B 92 (2019), Nr. 5, 1-11. https://doi.org/10.1140/epjb/e2019-100017-8

[26] ŁĄCKI, Mateusz ; DAMSKI, Bogdan ; ZAKRZEWSKI, Jakub: Locating the quantum critical point of the Bose-Hubbard model through singularities of simple observables. In: Sci. Rep. 6 (2016), Nr. 1, 1-8. https://doi.org/10.1038/srep38340

[27] KAWAKI, Keima ; KUNO, Yoshihito ; ICHINOSE, Ikuo: Phase diagrams of the extended Bose-Hubbard model in one dimension by Monte-Carlo simulation with the help of a stochastic-series expansion. In: Phys. Rev. B 95 (2017), Nr. 19, 195101. https://doi.org/10.1103/PhysRevB.95.195101

[28] DOGRA, Nishant ; BRENNECKE, Ferdinand ; HUBER, Sebastian D. ; DONNER, Tobias: Phase transitions in a Bose-Hubbard model with cavity-mediated global-range interactions. In: Phys. Rev. A 94 (2016), Nr. 2, 023632. https://doi.org/10.1103/PhysRevA.94.023632
[29] LV, Jian-Ping; CHEN, Qing-Hu; DENG, Youjin: Two-species hard-core bosons on the triangular lattice: A quantum Monte Carlo study. In: *Phys. Rev. A* 89 (2014), Nr. 1, 013628. https://doi.org/10.1103/PhysRevA.89.013628

[30] ROSSI, M; NAVA, M; REATTO, L; GALLI, DE: Exact ground state Monte Carlo method for Bosons without importance sampling. In: *J. Chem. Phys.* 131 (2009), Nr. 15, 154108. https://doi.org/10.1063/1.3247833

[31] ASSARAF, Roland; CAFFAREL, Michel: Computing forces with quantum Monte Carlo. In: *J. Chem. Phys.* 113 (2000), Nr. 10, 4028–4034. https://doi.org/10.1063/1.1286598

[32] FILIPPI, Claudia; UMRIGAR, CJ: Correlated sampling in quantum Monte Carlo: A route to forces. In: *Phys. Rev. B* 61 (2000), Nr. 24, R16291. https://doi.org/10.1103/PhysRevB.61.R16291

[33] SORELLA, Sandro; CAPRIOTTI, Luca: Algorithmic differentiation and the calculation of forces by quantum Monte Carlo. In: *J. Chem. Phys.* 133 (2010), Nr. 23, 234111. https://doi.org/10.1063/1.3516208

[34] ASSARAF, Roland; CAFFAREL, Michel: Zero-variance zero-bias principle for observables in quantum Monte Carlo: Application to forces. In: *J. Chem. Phys.* 119 (2003), Nr. 20, 10536–10552. https://doi.org/10.1063/1.1621615

[35] CHIESA, Simone; CEPERLEY, DM; ZHANG, Shiwei: Accurate, efficient, and simple forces computed with quantum Monte Carlo methods. In: *Phys. Rev. Lett.* 94 (2005), Nr. 3, 036404. https://doi.org/10.1103/PhysRevLett.94.036404

[36] CEPERLEY, DM; BERNU, B: The calculation of excited state properties with quantum Monte Carlo. In: *J. Chem. Phys.* 89 (1988), Nr. 10, 6316–6328. https://doi.org/10.1063/1.455398

[37] ZIMMERMAN, Paul M.; TOULOUSE, Julien; ZHANG, Zhiyong; MUSGRAVE, Charles B.; UMRIGAR, CJ: Excited states of methylene from quantum Monte Carlo. In: *J. Chem. Phys.* 131 (2009), Nr. 12, 124103. https://doi.org/10.1063/1.3220671

[38] WILLIAMSON, AJ; HOOD, Randolph Q.; NEEDS, RJ; RAJAGOPAL, G: Diffusion quantum Monte Carlo calculations of the excited states of silicon. In: *Phys. Rev. B* 57 (1998), Nr. 19, 12140. https://doi.org/10.1103/PhysRevB.57.12140

[39] GRIMES, RM; HAMMOND, BL; REYNOLDS, PJ; LESTER JR, WA: Quantum Monte Carlo approach to electronically excited molecules. In: *J. Chem. Phys.* 85 (1986), Nr. 8, 4749–4750. https://doi.org/10.1063/1.451754

[40] FELDT, Jonas; FILIPPI, Claudia: Excited-state calculations with quantum Monte Carlo. In: *arXiv:2002.03622* (2020). https://arxiv.org/abs/2002.03622

[41] GANGAT, AA; McCulloch, IP; MILBURN, GJ: Deterministic many-resonator W entanglement of nearly arbitrary microwave states via attractive Bose-Hubbard simulation. In: *Phys. Rev. X* 3 (2013), Nr. 3, 031009. https://doi.org/10.1103/PhysRevX.3.031009

[42] WILKINSON, Samuel A.; HARTMANN, Michael J.: Superconducting quantum many-body circuits for quantum simulation and computing. In: *Appl. Phys. Lett.* 116 (2020), Nr. 23, 230501. https://doi.org/10.1063/5.0008202

[43] ROUSHAN, P.; NEILL, C.; TANGPANITANON, J.; BASTIDAS, V. M.; MEGRANT, A.; BARENDs, R.; CHEN, Y.; CHEN, Z.; CHIARO, B.; DUNSWORTH, A.; FOWLER, A.; FOXEN, B.; GIUSTINA, M.; JEFFREY, E.; KELLY, J.; LUCERO, E.; MUTUS, J.; NEELEY, M.; QUINTANA, C.; SANK, D.; VAINSENCHER, A.; WENNER, J.; WHITE, T.; NEVEN, H.; ANGELAKIS, D. G.; MARTINIS, J.: Spectroscopic signatures of localization with interacting photons in superconducting qubits. In: *Science* 358 (2017), Nr. 6367, 1175–1179. https://doi.org/10.1126/science.aao1401

[44] HACOHEN-GOURGY, Shay; RAMASESH, Vinay V.; DE GRANDI, Claudia; SIDDIQI, Irfan; GIRVIN, Steve M.: Cooling and autonomous feedback in a bose-hubbard chain with attractive interactions.
In: *Phys. Rev. Lett.* 115 (2015), Nr. 24, 240501. https://doi.org/10.1103/PhysRevLett.115.240501

[45] DENG, Xiu-Hao; LAI, Chen-Yen; CHIEN, Chih-Chun: Superconducting circuit simulator of Bose-Hubbard model with a flat band. In: *Phys. Rev. B* 93 (2016), Nr. 5, 054116. https://doi.org/10.1103/PhysRevB.93.054116

[46] COTLER, Jordan; CHOI, Soonwon; LUKIN, Alexander; GHARIBYAN, Hrant; GROVER, Tarun; TAI, M. E.; RISPOLI, Matthew; SCHITTKO, Robert; PREISS, Philipp M.; KAUFMAN, Adam M.; GREINER, Markus; PICHLER, Hannes; HAYDEN, Patrick: Quantum Virtual Cooling. In: *Phys. Rev. X* 9 (2019), Jul, 031013. https://doi.org/10.1103/PhysRevX.9.031013

[47] KYRIENKO, Oleksandr: Quantum inverse iteration algorithm for programmable quantum simulators. In: *npj Quantum Information* 6 (2020), Nr. 1, 1-8. https://doi.org/10.1038/s41534-019-0239-7

[48] HUH, Joonsuk; GUERRESCHI, Gian G.; PEROPADRE, Borja; McCLEAN, Jarrod R.; ASPURU-GUZIK, Alán: Boson sampling for molecular vibronic spectra. In: *Nat. Photonics* 9 (2015), Nr. 9, 615-620. https://doi.org/10.1038/nphoton.2015.153

[49] HUH, Joonsuk; YUNG, Man-Hong: Vibronic boson sampling: generalized Gaussian boson sampling for molecular vibronic spectra at finite temperature. In: *Sci. Rep.* 7 (2017), Nr. 1, 1-10. https://doi.org/10.1038/s41598-017-07770-z

[50] SPARROW, Chris; MARTÍN-LÓPEZ, Enrique; MARAVIGLIA, Nicola; NEVILLE, Alex; HARROLD, Christopher; CAROLAN, Jacques; JOGLEKAR, Yogesh N.; HASHIMOTO, Toshikazu; MATSUDA, Nobuyuki; O'BRIEN, Jeremy L. u. a.: Simulating the vibrational quantum dynamics of molecules using photonics. In: *Nature* 557 (2018), Nr. 7707, 660-667. https://doi.org/10.1038/s41586-018-0152-9

[51] QUESADA, Nicolás: Franck-Condon factors by counting perfect matchings of graphs with loops. In: *J. Chem. Phys.* 150 (2019), Nr. 16, 164113. https://doi.org/10.1063/1.5086387

[52] JAHANGIRI, Soran; ARRAZOLA, Juan M.; QUESADA, Nicolás; DELGADO, Alain: Quantum algorithm for simulating molecular vibrational excitations. In: *Phys. Chem. Chem. Phys.* 22 (2020), Nr. 44, 25528-25537. https://doi.org/10.1039/DOCP03593A

[53] JAHANGIRI, Soran; ARRAZOLA, Juan M.; DELGADO, Alain: Quantum Algorithm for Simulating Single-Molecule Electron Transport. In: *J. Phys. Chem. Lett.* 12 (2021), Nr. 4, 1256-1261. https://doi.org/10.1021/acs.jpclett.0c03724

[54] KALAJDZIEVSKI, Timjan; WEEDBROOK, Christian; REBENTROST, Patrick: Continuous-variable gate decomposition for the Bose-Hubbard model. In: *Phys. Rev. A* 97 (2018), Nr. 6, 062311. https://doi.org/10.1103/PhysRevA.97.062311

[55] ZHONG, Han-Sen; WANG, Hui; DENG, Yu-Hao; CHEN, Ming-Cheng; PENG, Li-Chao; LUO, Yi-Han; QIN, Jian; WU, Dian; DING, Xing; HU, Yi u. a.: Quantum computational advantage using photons. In: *Science* 370 (2020), Nr. 6523, 1460-1463. https://doi.org/10.1126/science.abe8770

[56] KASSAL, Ivan; ASPURU-GUZIK, Alán: Quantum algorithm for molecular properties and geometry optimization. In: *J. Chem. Phys.* 131 (2009), Nr. 22, 224102. https://doi.org/10.1063/1.3266959

[57] MITARAI, Kosuke; NAKAGAWA, Yuya O.; MIZUKAMI, Wataru: Theory of analytical energy derivatives for the variational quantum eigensolver. In: *Phys. Rev. Res.* 2 (2020), Nr. 1, 013129. https://doi.org/10.1103/PhysRevResearch.2.013129

[58] PARRISH, Robert M.; HOHENSTEIN, Edward G.; MCMAHON, Peter L.; MARTINEZ, Todd J.: Hybrid Quantum/Classical Derivative Theory: Analytical Gradients and Excited-State Dynamics for the Multistate Contracted Variational Quantum Eigensolver. In: *arXiv:1906.08728* (2019). https://arxiv.org/abs/1906.08728
[59] SOKOLOV, Igor O. ; BARKOUTSOS, Panagiotis K. ; MOELLER, Lukas ; SUCHS-LAND, Philippe ; MAZZOLA, Guglielmo ; TAVERNELLI, Ivano: Microcanonical and finite-temperature ab initio molecular dynamics simulations on quantum computers. In: *Phys. Rev. Research* 3 (2021), Feb, 013125. https://doi.org/10.1103/PhysRevResearch.3.013125

[60] McCLEAN, Jarrod R. ; KIMCHI-SCHWARTZ, Mollie E. ; CARTER, Jonathan ; JONG, Wibe A.: Hybrid quantum-classical hierarchy for mitigation of decoherence and determination of excited states. In: *Phys. Rev. A* 95 (2017), Nr. 4, 042308. https://doi.org/10.1103/PhysRevA.95.042308

[61] OLLITRAULT, Pauline J. ; KANDALA, Abhinav ; CHEN, Chun-Fu ; BARKOUTSOS, Panagiotis K. ; MEZZACAPO, Antonio ; PISTOIA, Marco ; SHELDON, Sarah ; WOERNER, Stefan ; GAMBETTA, Jay M. ; TAVERNELLI, Ivano: Quantum equation of motion for computing molecular excitation energies on a noisy quantum processor. In: *Phys. Rev. Research* 2 (2020), Oct, 043140. https://doi.org/10.1103/PhysRevResearch.2.043140

[62] NAKANISHI, Ken M. ; MITARAI, Kosuke ; FUJII, Keisuke: Subspace-search variational quantum eigensolver for excited states. In: *Phys. Rev. Res.* 1 (2019), Nr. 3, 033062. https://doi.org/10.1103/PhysRevResearch.1.033062

[63] IBE, Yohei ; NAKAGAWA, Yuya O. ; YAMAMOTO, Takahiro ; MITARAI, Kosuke ; GAO, Qi ; KOBAYASHI, Takao: Calculating transition amplitudes by variational quantum eigensolvers. In: *arXiv:2002.11724* (2020). https://arxiv.org/abs/2002.11724

[64] LEE, Joonho ; HUGGINS, William J. ; HEAD-GORDON, Martin ; WHALEY, K B.: Generalized unitary coupled cluster wave functions for quantum computation. In: *J. Chem. Theory Comput.* 15 (2018), Nr. 1, 311–324. https://doi.org/10.1021/acs.jctc.8b01004

[65] HIGGOTT, Oscar ; WANG, Daochen ; BRIERLEY, Stephen: Variational quantum

computation of excited states. In: *Quantum* 3 (2019), 156. https://doi.org/10.22331/q-2019-07-01-156

[66] JONES, Tyson ; ENDO, Suguru ; MCARDLE, Sam ; YUAN, Xiao ; BENJAMIN, Simon C.: Variational quantum algorithms for discovering Hamiltonian spectra. In: *Phys. Rev. A* 99 (2019), Nr. 6, 062304. https://doi.org/10.1103/PhysRevA.99.062304

[67] JOUZDANI, Pejman ; BRINGUIER, Stefan ; KOSTUK, Mark: A Method of Determining Excited-States for Quantum Computation. In: *arXiv:1908.05238* (2019). https://arxiv.org/abs/1908.05238

[68] PARRISH, Robert M. ; HOHENSTEIN, Edward G. ; MCMAHON, Peter L. ; MARTÍNEZ, Todd J.: Quantum computation of electronic transitions using a variational quantum eigensolver. In: *Phys. Rev. Lett.* 122 (2019), Nr. 23, 230401. https://doi.org/10.1103/PhysRevLett.122.230401

[69] BAUMAN, Nicholas P. ; LOW, Guang H. ; KOWALSKI, Karol: Quantum simulations of excited states with active-space downfolded Hamiltonians. In: *J. Chem. Phys.* 151 (2019), Nr. 23, 234114. https://doi.org/10.1063/1.5128103

[70] MOTTA, Mario ; SUN, Chong ; TAN, Adrian T. ; O’ROURKE, Matthew J. ; YE, Erika ; MINNICH, Austin J. ; BRANDÃO, Fernando G. ; CHAN, Garnet Kin-Lic: Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution. In: *Nat. Phys.* 16 (2020), Nr. 2, 205–210. https://doi.org/10.1038/s41567-019-0704-4

[71] ZHANG, Dan-Bo ; YUAN, Zhan-Hao ; YIN, Tao: Variational quantum eigensolvers by variance minimization. In: *arXiv:2006.15781* (2020). https://arxiv.org/abs/2006.15781

[72] YALOUZ, Saad ; SENJEAN, Bruno ; GÜNTHER, Jakob ; BUDA, Francesco ; O'BRIEN, Thomas E. ; VISSCHER, Lucas: A state-averaged orbital-optimized hybrid quantum-classical algorithm for a democratic description of ground and excited states. In: *Quantum Sci. Technol.* 6
(2021), Nr. 2, 024004. https://doi.org/10.1088/2058-9565/abd334

[73] ARRAZOLA, JM ; BERGHOLM, V ; BRÁDLER, K ; BROMLEY, TR ; COLLINS, MJ ; DHAND, I ; FUMAGALLI, A ; GERRITS, T ; GOUSSEV, A ; HELT, LG u.a.: Quantum circuits with many photons on a programmable nanophotonic chip. In: Nature 591 (2021), Nr. 7848, 54–60. https://doi.org/10.1038/s41586-021-03202-1

[74] WEEDBROOK, Christian ; PIRANDOLA, Stefano ; GARCÍA-PATRÓN, Raúl ; CERF, Nicolas J. ; RALPH, Timothy C. ; SHAPIRO, Jeffrey H. ; LLOYD, Seth: Gaussian quantum information. In: Rev. Mod. Phys. 84 (2012), Nr. 2, 621. https://doi.org/10.1103/RevModPhys.84.621

[75] BROMLEY, Thomas R. ; ARRAZOLA, Juan M. ; JAHANGIRI, Soran ; IZAAC, Josh ; QUESADA, Nicolás ; GRAN, Alain D. ; SCHULD, Maria ; SWINARTON, Jeremy ; ZABANEH, Zeid ; KILLORAN, Nathan: Applications of near-term photonic quantum computers: software and algorithms. In: Quantum Sci. Technol. 5 (2020), Nr. 3, 034010. https://doi.org/10.1088/2058-9565/ab8504

[76] PFISTER, Olivier: Continuous-variable quantum computing in the quantum optical frequency comb. In: J. Phys. B 53 (2019), Nr. 1, 012001. https://doi.org/10.1088/1361-6455/ab526f

[77] MARSHALL, Kevin ; POOSER, Raphael ; SIOPSIS, George ; WEEDBROOK, Christian: Quantum simulation of quantum field theory using continuous variables. In: Phys. Rev. A 92 (2015), Nr. 6, 063825. https://doi.org/10.1103/PhysRevA.92.063825

[78] BRENNEN, Gavin K. ; ROHDE, Peter ; SANDERS, Barry C. ; SINGH, Sukhwinder: Multiscale quantum simulation of quantum field theory using wavelets. In: Phys. Rev. A 92 (2015), Nr. 3, 032315. https://doi.org/10.1103/PhysRevA.92.032315

[79] KÜHNER, Till D. ; MONIEN, H: Phases of the one-dimensional Bose-Hubbard model. In: Phys. Rev. B 58 (1998), Nr. 22, R14741. https://doi.org/10.1103/PhysRevB.58.R14741

[80] KÜHNER, Till D. ; WHITE, Steven R. ; MONIEN, Hartmut: One-dimensional Bose-Hubbard model with nearest-neighbor interaction. In: Phys. Rev. B 61 (2000), Nr. 18, 12474. https://doi.org/10.1103/PhysRevB.61.12474

[81] JACK, Michael W. ; YAMASHITA, Makoto: Bose-Hubbard model with attractive interactions. In: Phys. Rev. A 71 (2005), Nr. 2, 023610. https://doi.org/10.1103/PhysRevA.71.023610

[82] OELKERS, Norman ; LINKS, Jon: Ground-state properties of the attractive one-dimensional Bose-Hubbard model. In: Phys. Rev. B 75 (2007), Nr. 11, 115119. https://doi.org/10.1103/PhysRevB.75.115119

[83] JULIÁ-DÍAZ, B ; DAGNINO, D ; LEWENSTEIN, M ; MARTORELL, J ; POLLS, A: Macroscopic self-trapping in Bose-Einstein condensates: Analysis of a dynamical quantum phase transition. In: Phys. Rev. A 81 (2010), Nr. 2, 023615. https://doi.org/10.1103/PhysRevA.81.023615

[84] MELÉ-MESSEGUER, Marina ; JULIÁ-DÍAZ, Bruno ; POLLS, Artur: Improved variational approach to the two-site Bose-Hubbard model. In: J. Low Temp. Phys. 165 (2011), Nr. 5, 180–194. https://doi.org/10.1007/s10909-011-0395-4

[85] MANSIKKAMÄKI, Olli ; LAINE, Sami ; SILVERI, Matti: The phases of the disordered Bose-Hubbard model with attractive interactions. In: arXiv:2101.06032 (2021). https://arxiv.org/abs/2101.06032

[86] THEIS, M ; THALHAMMER, G ; WINKLER, K ; HELLWIG, M ; RUFF, G ; GRIMM, R ; DENSCHLAG, J H.: Tuning the scattering length with an optically induced Feshbach resonance. In: Phys. Rev. Lett. 93 (2004), Nr. 12, 123001. https://doi.org/10.1103/PhysRevLett.93.123001

[87] GÓRAL, Krzysztof ; SANTOS, Luis ; LEWENSTEIN, Maciej: Quantum phases of dipolar bosons in optical lattices. In: Phys. Rev. Lett. 88 (2002), Nr. 17, 170406. https://doi.org/10.1103/PhysRevLett.88.170406
[88] BÜCHLER, HP ; BLATTER, G: Supersolid versus phase separation in atomic Bose-Fermi mixtures. In: Phys. Rev. Lett. 91 (2003), Nr. 13, 130404. https://doi.org/10.1103/PhysRevLett.91.130404

[89] BLOCH, Immanuel: Ultracold quantum gases in optical lattices. In: Nature Phys. 1 (2005), Nr. 1, 23–30. https://doi.org/10.1038/nphys138

[90] BLOCH, Immanuel ; DALIBARD, Jean ; ZWERGER, Wilhelm: Many-body physics with ultracold gases. In: Rev. Mod. Phys. 80 (2008), Nr. 3, 885. https://doi.org/10.1103/RevModPhys.80.885

[91] BAIER, Simon ; MARK, Manfred J. ; PETTER, Daniel ; AIKAWA, Kiyotaka ; CHOMAZ, Lauriane ; CAI, Zi ; BARA NOV, M ; ZOLLER, P ; FERLAINO, F: Extended Bose-Hubbard models with ultracold magnetic atoms. In: Science 352 (2016), Nr. 6282, 201–205. https://doi.org/10.1126/science.aac9812

[92] POUTHIER, Vincent: Two-vibron bound states in α-helix proteins: the interplay between the intramolecular anharmonicity and the strong vibron-phonon coupling. In: Phys. Rev. E 68 (2003), Nr. 2, 021909. https://doi.org/10.1103/PhysRevE.68.021909

[93] EDLER, J ; PFISTER, R ; POUTHIER, V ; FALVO, C ; HAMM, P: Direct observation of self-trapped vibrational states in α-helices. In: Phys. Rev. Lett. 93 (2004), Nr. 10, 106405. https://doi.org/10.1103/PhysRevLett.93.106405

[94] MATSUDA, Hirotsugu ; TSUNETO, Toshihiko: Off-diagonal long-range order in solids. In: Prog. Theor. Phys. Suppl. 46 (1970), 411–436. https://doi.org/10.1143/PTPS.46.411

[95] LIU, Kao-Shien ; FISHER, Michael E.: Quantum lattice gas and the existence of a supersolid. In: J. Low Temp. Phys. 10 (1973), Nr. 5, 655–683. https://doi.org/10.1007/BF00655458

[96] RAVENTÓS, David ; GRASS, Tobias ; LEWENSTEIN, Maciej ; JULIÁ-DÍAZ, Bruno: Cold bosons in optical lattices: a tutorial for exact diagonalization. In: J. Phys. B 50 (2017), Nr. 11, 113001. https://doi.org/10.1088/1361-6455/aa68b1

[97] ZHANG, JM ; DONG, RX: Exact diagonalization: the Bose-Hubbard model as an example. In: Eur. J. Phys 31 (2010), Nr. 3, 591. https://doi.org/10.1088/0143-0807/31/3/016

[98] BAUER, Bela ; NAYAK, Chetan: Area laws in a many-body localized state and its implications for topological order. In: J. Stat. Mech.: Theory Exp 2013 (2013), Nr. 09, P09005. https://doi.org/10.1088/1742-5468/2013/09/P09005

[99] HOPJAN, Miroslav ; HEIDRICH-MEISNER, Fabian: Many-body localization from a one-particle perspective in the disordered one-dimensional Bose-Hubbard model. In: Phys. Rev. A 101 (2020), Nr. 6, 063617. https://doi.org/10.1103/PhysRevA.101.063617

[100] GIRAUD, Olivier ; MARTIN, John ; GEORGEOT, Bertrand: Entropy of entanglement and multifractal exponents for random states. In: Phys. Rev. A 79 (2009), Nr. 3, 032308. https://doi.org/10.1103/PhysRevA.79.032308

[101] GIRAUD, Olivier ; MARTIN, John ; GEORGEOT, Bertrand: Entanglement of localized states. In: Phys. Rev. A 76 (2007), Nr. 4, 042333. https://doi.org/10.1103/PhysRevA.76.042333

[102] LUKIN, Alexander ; RISPOLI, Matthew ; SCHITTKO, Robert ; TAI, M E. ; KAUFMAN, Adam M. ; CHOI, Soonwon ; KHEMANI, Vedika ; LÉONARD, Julian ; GREINER, Markus: Probing entanglement in a many-body-localized system. In: Science 364 (2019), Nr. 6437, 256–260. https://doi.org/10.1126/science.aau0818

[103] ISLAM, Rajibul ; MA, Ruichao ; PREISS, Philipp M. ; TAI, M E. ; LUKIN, Alexander ; RISPOLI, Matthew ; GREINER, Markus: Measuring entanglement entropy in a quantum many-body system. In: Nature 528 (2015), Nr. 7580, 77–83. https://doi.org/10.1038/nature15750

[104] BEUGELING, Wouter ; ANDREANOV, Alexei ; HAQUE, Masudul: Global characteristics of all eigenstates of local many-
body Hamiltonians: participation ratio and entanglement entropy. In: *J. Stat. Mech.: Theory Exp.* 2015 (2015), Nr. 2, P02002. https://doi.org/10.1088/1742-5468/2015/02/P02002

[105] CEREZO, Marco; ARRASMITH, Andrew; BABBUSH, Ryan; BENJAMIN, Simon C. ; ENDO, Suguru; FUJII, Keisuke; MC-CLEAN, Jarrod R. ; MITARAI, Kosuke; YUAN, Xiao; CINCIO, Lukasz u. a.: Variational quantum algorithms. In: *Nature Reviews Physics* (2021), 1–20. https://doi.org/10.1038/s42254-021-00348-9

[106] EVANGELISTA, Francesco A. ; CHAN, Gar- net Kin-Lic ; SCUSERIA, Gustavo E. : Exact parameterization of fermionic wave func- tions via unitary coupled cluster theory. In: *J. Chem. Phys.* 151 (2019), Nr. 24, 244112. https://doi.org/10.1063/1.5133059

[107] SOKOLOV, Igor O. ; BARKOUTSOS, Panagi- otis K. ; OLLITRAULT, Pauline J. ; GREEN- BERG, Donny ; RICE, Julia ; PISTOIA, Marco ; TAVERNELLI, Ivano: Quantum orbital-optimized unitary coupled cluster methods in the strongly correlated regime: Can quantum algorithms outperform their classical equivalents? In: *J. Chem. Phys.* 152 (2020), Nr. 12, 124107. https://doi.org/10.1063/1.5141835

[108] MIZUKAMI, Wataru ; MITARAI, Kosuke ; NAKAGAWA, Yuya O. ; YAMAMOTO, Takahiro ; YAN, Tennin ; OHNISHI, Yu- ya: Orbital optimized unitary cou- pled cluster theory for quantum computer. In: *Phys. Rev. Research* 2 (2020), Nr. 3, 033421. https://doi.org/10.1103/PhysRevResearch.2.033421

[109] AZUMA, Hiroo: Quantum computation with Kerr-nonlinear photonic crystals. In: *J. Phys. D: Appl. Phys.* 41 (2007), Nr. 2, 025102. https://doi.org/10.1088/0022-3727/41/2/025102

[110] SEFI, Seckin; VAIBHAV, Vishal; LOOCK, Peter van: Measurement-induced optical Kerr interaction. In: *Physical Review A* 88 (2013), Nr. 1, 012303. https://doi.org/10.1103/PhysRevA.88.012303

[111] HILLMANN, Timo; QUIJANDRÍA, Fernando ; JOHANSSON, Göran; FERRARO, Alessan- dro; GASPARINETTI, Simone; FERRINI,

Giulia: Universal gate set for continuous-variable quantum computation with microwave circuits. In: *Physical review letters* 125 (2020), Nr. 16, 160501. https://doi.org/10.1103/PhysRevLett.125.160501

[112] CLEMENTS, William R. ; HUMPHREYS, Peter C. ; METCALF, Benjamin J. ; KOLTHAMMER, W S. ; WALMSLEY, Ian A.: Optimal design for universal multiport interferometers. In: *Optica* 3 (2016), Nr. 12, 1460–1465. https://doi.org/10.1364/OPTICA.3.001460

[113] KILLORAN, Nathan ; IZAAC, Josh ; QUESADA, Nicolás ; BERGHOLM, Ville ; AMY, Matthew ; WEEDBROOK, Christian: Strawberry fields: A software platform for photonic quantum computing. In: *Quantum* 3 (2019), 129. https://doi.org/10.22331/q-2019-03-11-129

[114] GERBER, Florian ; FURRER, Reinhard: optimParallel: An R Package Providing a Parallel Version of the L-BFGS-B Optimization Method. In: *R J.* 11 (2019), Nr. 1, 352–358. https://doi.org/10.32614/RJ-2019-03-11-129

[115] YALOUZ, Saad: Code online. (2021). https://github.com/SYalouz/VQA_CV_BoseHubbard

[116] MCCLEAN, Jarrod R. ; BOIXO, Sergio ; SMELYANSKIY, Vadim N. ; BABBUSH, Ryan ; NEVEN, Hartmut: Barren plateaus in quantum neural network training landscapes. In: *Nat. Commun.* 9 (2018), Nr. 1, 1–6. https://doi.org/10.1038/s41467-018-07090-4

[117] SHARMA, Kunal; CEREZO, Marco; CINCIO, Lukasz; COLES, Patrick J.: Trainability of dissipative perceptron-based quantum neural networks. In: arXiv:2005.12458 (2020). https://arxiv.org/abs/2005.12458

[118] CEREZO, Marco; SONE, Akira; VOLKOFF, Tyler; CINCIO, Lukasz; COLES, Patrick J.: Cost function dependent barren plateaus in shallow parametrized quantum circuits. In: *Nat. Commun.* 12 (2021), Nr. 1, 1–12. https://doi.org/10.1038/s41467-021-21728-w
[119] ARRASMITH, Andrew ; CEREZO, M ; CZARNIK, Piotr ; CINCIO, Lukasz ; COLES, Patrick J.: Effect of barren plateaus on gradient-free optimization. In: *Quantum* 5 (2021), 558. https://doi.org/10.22331/q-2021-10-05-558

[120] UVAROV, AV ; BIAMONTE, Jacob D.: On barren plateaus and cost function locality in variational quantum algorithms. In: *J. Phys. A Math* 54 (2021), Nr. 24, 245301. https://doi.org/10.1088/1751-8121/abfac7

[121] CEREZO, M ; COLES, Patrick J.: Higher order derivatives of quantum neural networks with barren plateaus. In: *Quantum Sci. Technol.* 6 (2021), Nr. 3, 035006. https://doi.org/10.1088/2058-9565/abf51a

[122] LAROCCA, Martin ; CZARNIK, Piotr ; SHARMA, Kunal ; MURALEEDHARAN, Gopikrishnan ; COLES, Patrick J. ; CEREZO, M: Diagnosing barren plateaus with tools from quantum optimal control. In: *arXiv:2105.14377* (2021). https://arxiv.org/abs/2105.14377

[123] WANG, Samson ; FONTANA, Enrico ; CEREZO, Marco ; SHARMA, Kunal ; SONE, Akira ; CINCIO, Lukasz ; COLES, Patrick J.: Noise-induced barren plateaus in variational quantum algorithms. In: *arXiv preprint arXiv:2007.14384* (2020). https://arxiv.org/abs/2007.14384

[124] HOLMES, Zoë ; SHARMA, Kunal ; CEREZO, M ; COLES, Patrick J.: Connecting ansatz expressibility to gradient magnitudes and barren plateaus. In: *arXiv:2101.02138* (2021). https://arxiv.org/abs/2101.02138

[125] VOLKOFF, Tyler J.: Efficient trainability of linear optical modules in quantum optical neural networks. In: *Journal of Russian Laser Research* 42 (2021), Nr. 3, 250–260. https://doi.org/10.1007/s10946-021-09958-1

[126] VOLKOFF, Tyler ; HOLMES, Zoë ; SORNBORGER, Andrew: Universal compiling and (No-) Free-Lunch theorems for continuous variable quantum learning. In: *arXiv:2105.01049* (2021). https://arxiv.org/abs/2105.01049

[127] HAUG, Tobias ; KIM, M S.: Optimal training of variational quantum algorithms without barren plateaus. In: *arXiv:2104.14543* (2021). https://arxiv.org/abs/2104.14543

[128] PATTI, Taylor L. ; NAJAFI, Khadijeh ; GAO, Xun ; YELIN, Susanne F.: Entanglement devised barren plateau mitigation. In: *Phys. Rev. Research* 3 (2021), Jul, 033090. https://doi.org/10.1103/PhysRevResearch.3.033090

[129] AARONSON, Scott ; ARKHIPOV, Alex: *The computational complexity of linear optics*. 2011. – 333–342 S. https://doi.org/10.1145/1993636.1993682

[130] GARD, Bryan T. ; OLSON, Jonathan P. ; CROSS, Robert M. ; KIM, Moochan B. ; LEE, Hwang ; DOWLING, Jonathan P.: Inefficiency of classically simulating linear optical quantum computing with Fock-state inputs. In: *Phys. Rev. A* 89 (2014), Nr. 2, 022328. https://doi.org/10.1103/PhysRevA.89.022328

[131] JÖNSSON, Mattias ; BJÖRK, Gunnar: Evaluating the performance of photon-number-resolving detectors. In: *Phys. Rev. A* 99 (2019), Nr. 4, 043822. https://doi.org/10.1103/PhysRevA.99.043822

[132] PROVAZNÍK, Jan ; LACHMAN, Lukáš ; FILIP, Radim ; MAREK, Petr: Benchmarking photon number resolving detectors. In: *Opt. Express* 28 (2020), Nr. 10, 14839–14849. https://doi.org/10.1364/OE.389619

[133] SLUSSARENKO, Sergei ; PRYDE, Geoff J.: Photonic quantum information processing: A concise review. In: *Applied Physics Reviews* 6 (2019), Nr. 4, 041303. https://doi.org/10.1063/1.5115814

[134] CHEN, Wei-Peng ; SINGLETON, Jared ; QIN, Lei ; CAMÓN, Agustín ; ENGELHARDT, Larry ; LUIS, Fernando ; WINPENNY, Richard E. ; ZHENG, Yan-Zhen: Quantum Monte Carlo simulations of a giant {Ni 21 Gd 20} cage with a S= 91 spin ground state. In: *Nat. Commun.* 9 (2018), Nr. 1, 1–6. https://doi.org/10.1038/s41467-018-04547-4
[135] ROSSINI, Davide; FAZIO, Rosario: Phase diagram of the extended Bose-Hubbard model. In: *New J. Phys.* 14 (2012), Nr. 6, 065012. https://doi.org/10.1088/1367-2630/14/6/065012

[136] OHGOE, Takahiro; SUZUKI, Takafumi; KAWASHIMA, Naoki: Ground-state phase diagram of the two-dimensional extended bose-hubbard model. In: *Phys. Rev. B* 86 (2012), Nr. 5, 054520. https://doi.org/10.1103/PhysRevB.86.054520

[137] LIN, Fei; MAIER, TA; SCAROLA, Vito W.: Disordered supersolids in the extended Bose-Hubbard model. In: *Sci. Rep.* 7 (2017), Nr. 1, 1–10. https://doi.org/10.1038/s41598-017-13040-9

[138] SCAROLA, Vito W.; SARMA, S D.: Quantum phases of the extended Bose-Hubbard Hamiltonian: Possibility of a supersolid state of cold atoms in optical lattices. In: *Phys. Rev. Lett.* 95 (2005), Nr. 3, 033003. https://doi.org/10.1103/PhysRevLett.95.033003

[139] SCAROLA, VW; DEMLER, E; SARMA, S D.: Searching for a supersolid in cold-atom optical lattices. In: *Phys. Rev. A* 73 (2006), Nr. 5, 051601. https://doi.org/10.1103/PhysRevA.73.051601

[140] VAN OTTERLO, Anne; WAGENBLAST, Karl-Heinz; BALTIN, Reinhard; BRUDER, C; FAZIO, Rosario; SCHÖN, Gerd: Quantum phase transitions of interacting bosons and the supersolid phase. In: *Phys. Rev. B* 52 (1995), Nr. 22, 16176. https://doi.org/10.1103/PhysRevB.52.16176

[141] RODDICK, Eric; STROUD, David: Supersolid phases in underdamped Josephson arrays: Quantum Monte Carlo simulations. In: *Phys. Rev. B* 51 (1995), Nr. 13, 8672. https://doi.org/10.1103/PhysRevB.51.8672

[142] MICNAS, R; RANNINGER, J; ROBASZKIEWICZ, St: Superconductivity in narrow-band systems with local non-retarded attractive interactions. In: *Rev. Mod. Phys.* 62 (1990), Nr. 1, 113. https://doi.org/10.1103/RevModPhys.62.113

[143] HANSEN, Nikolaus: The CMA evolution strategy: A tutorial. In: arXiv:1604.00772 (2016). https://arxiv.org/abs/1604.00772

[144] LAVRIJSEN, Wim; TUDOR, Ana; MÜLLER, Juliane; IANCU, Costin; JONG, Wibe de: *Classical optimizers for noisy intermediate-scale quantum devices*. IEEE, 2020. – 267–277 S. https://doi.org/10.1109/QCE49297.2020.00041

[145] SUNG, Kevin J. ; YAO, Jiahao ; HARRIGAN, Matthew P. ; RUBIN, Nicholas C. ; JIANG, Zhang ; LIN, Lin ; BABBUSH, Ryan ; MCCLEAN, Jarrod R.: Using models to improve optimizers for variational quantum algorithms. In: *Quantum Sci. Technol.* 5 (2020), Nr. 4, 044008. https://doi.org/10.1088/2058-9565/abb6d9

[146] KWON, Yongkyung; CEPERLEY, David M.: 4 He adsorption on a single graphene sheet: Path-integral Monte Carlo study. In: *Phys. Rev. B* 85 (2012), Nr. 22, 224501. https://doi.org/10.1103/PhysRevB.85.224501

[147] GORDILLO, MC; BORONAT, Jordi: He 4 on a Single Graphene Sheet. In: *Phys. Rev. Lett.* 102 (2009), Nr. 8, 085303. https://doi.org/10.1103/PhysRevLett.102.085303