#  Training Gaussian boson sampling by quantum machine learning

Claudio Conti<EMAIL_ADDRESS>Department of Physics, University
Sapienza, P.le Aldo Moro 5, 00185 Rome, Italy  Institute for Complex Systems,
National Research Council (ISC-NR), Via dei Taurini 19, 00185 Rome, Italy
Research Center Enrico Fermi, Via Panisperna 89a, 00184 Rome (Italy)
https://www.complexlight.org

###### Abstract

We use neural networks to represent the characteristic function of many-body
Gaussian states in the quantum phase space. By a pullback mechanism, we model
transformations due to unitary operators as linear layers that can be cascaded
to simulate complex multi-particle processes. We use the layered neural
networks for non-classical light propagation in random interferometers, and
compute boson pattern probabilities by automatic differentiation. This is a
viable strategy for training Gaussian boson sampling. We demonstrate that
multi-particle events in Gaussian boson sampling can be optimized by a proper
design and training of the neural network weights. The results are potentially
useful to the creation of new sources and complex circuits for quantum
technologies.

Machine Learning and Gaussian Boson Sampling

##  I Introduction

The development of new models and tools for machine learning (ML) is
surprisingly affecting the study of many-body quantum systems and quantum
optics [1]. Neural networks (NN) enable representations of high-dimensional
systems and furnish a universal ansatz for many purposes, like finding the
ground state of many-body Hamiltonians [2], including dissipative systems [3,
4].

Unsupervised and supervised learning endow new designs for quantum circuits
[5], metrology and cryptography [6, 7], multilevel gates [8], and Bell tests
[9]. NN are also triggering new fundamental investigations in quantum
neuromorphic and wave computing [10, 11, 12, 13, 14, 15], quantum
thermodynamics [16], and topological photonics [17].

The impact of ML in quantum optics and many-body physics is related to the
versatile representation that the NN models furnish for functions of an
arbitrary number of variables. Also, the powerful application programming
interfaces (APIs), as TensorFlow, enable many new features and tools to
compute and design many-body Hamiltonians or large-scale quantum gates [18].

Here we show that NN models are also useful when considering representations
in the phase space, as the characteristic functions $\chi$ or the
Q-representation [19]. Unitary operators, as squeezers or displacers, act on
the phase-space as variable transformations that correspond to layers in the
NN model. Hence, a multilayer NN may encode phase-space representations of
complex many-body states. This encoding has two main advantages: on the one
hand, one can quickly build complex quantum states by combining NN layers; on
the other hand, one can use the automatic graph building and API
differentiation technology to compute observables. Also, graphical and tensor
processing units (GPU and TPU) may speed up the computation.

In the following, we show how to compute the probability of multi-particle
patterns when Gaussian states propagate in a system made of squeezers and
interferometers. This problem corresponds to the renowned Gaussian Boson
sampling [20, 21], which recently demonstrated the quantum advantage at an
impressing scale [22], following earlier realizations [23, 24, 25, 26, 27, 28]
of the original proposal by Aharanson and Arkhipov [29]. The theory of
Gaussian Boson sampling (GBS) heavily relies on phase-space methods [30],
making it an exciting NN test-bed supported by recently reported trainable
hardware [31, 32, 33].

A notable outcome of adopting NN models in the phase space is the possibility
of training multi-particle statistics [34] and other features as the degree of
entanglement. Indeed, most of the reported investigations in quantum ML, focus
either on using NN models as a variational ansatz or tailoring the
input/output response of a quantum gate. On the contrary, ML in the phase
space permits optimizing many-particle features, for example, to increase the
probability of multi-photon events. NN may open new strategies to generate
non-classical light or enhance the probability of observing large-scale
entanglement with relevance in many applications. Here, we derive the NN
representing the characteristic function of the Gaussian boson sampling setup.
Proper NN training increases the photon-pair probability by orders of
magnitude.

Fig. 1 shows the general workflow of the proposed methodology, the different
steps enable to define a trainable model for optmizing Gaussian boson
sampling. In Section II, we introduce the way we adopt a neural network to
compute the characteristic function. In Sec. III, we detail how to compute the
observable as derivatives of the characteristic function neural network. In
Sec. IV, we show how to compute the Gaussian boson sampling patterns. In Sec.
V, we introduce the loss function and describe the training of the model to
optimize specific patterns. Conclusions are drawn in Sec. VI.

Figure 1: Workflow of the proposed methodology to train Boson sampling by
representing the characteristic function as a neural network.

##  II Characteristic function as a neural network

In the phase space, we represent a $n$-body state by complex characteristic
function $\chi({\bf x})=\chi_{R}({\bf x})+\imath\chi_{I}({\bf x})$ of a real
vector ${\bf x}$ [35, 19]. ${\bf x}$ has dimension $1\times N$ with $N=2n$.
For Gaussian states [36]

$\chi({\bf x})=e^{-\frac{1}{4}{\bf x}{\bf g}{\bf x}^{\top}+\imath\,{\bf x}{\bf
d}}.$ (1)

with ${\bf g}$ the real covariance $N\times N$ matrix, and ${\bf d}$ the real
displacement $N\times 1$ vector. In our notation, we omit the symbols of the
dot product such that ${\bf x}{\bf d}$ and ${\bf x}{\bf g}{\bf x}^{\top}$ are
scalars. One has ($j,k=0,1,2\dots,N-1$)

$\langle\hat{R}_{j}\rangle=d_{j}=\left.\frac{\partial\chi}{\partial
x_{j}}\right|_{{\bf x}=0}\text{,}$ (2)

and

$g_{jk}=2\langle(\hat{R}_{j}-d_{j})(\hat{R}_{k}-d_{k})\rangle-\imath J_{jk},$
(3)

being ${\bf J}=\bigoplus_{j=0}^{n-1}{\bf J}_{1}$, ${\bf
J}_{1}=\big{(}\begin{smallmatrix}0&1\\\ -1&0\end{smallmatrix}\big{)}$ [36]. In
Eq. (2), the canonical variables, $\hat{q}_{j}=\hat{R}_{2j}$ and
$\hat{p}_{j}=\hat{R}_{2j+1}$, with $j=0,1,\dots,n-1$, are organized in the
$N\times 1$ operator array $\hat{{\bf R}}$. As shown in Fig. 2a, the
characteristic function is a NN layer with two real outputs $\chi_{R}$ and
$\chi_{I}$. The $\chi$ layer has two inputs: ${\bf x}$, and a auxiliary bias
$N\times 1$ vector a, for later convenience.

The vacuum state is a Gaussian state with ${\bf g}={\bf 1}$ and ${\bf d}={\bf
0}$. From the vacuum, one can generate specific states by unitary operators,
as displacement or squeezing operators. These transform the canonical
variables as $\hat{\widetilde{\bf R}}={\bf M}\,\hat{{\bf R}}+{\bf
d}^{\prime}$, where the symplectic matrix ${\bf M}$ and the vector ${\bf
d}^{\prime}$ depend on the specific operator (detailed, e.g., in [36]). The
characteristic function changes as

$\tilde{\chi}({\bf x})=\chi({\bf x}{\bf M})e^{\imath{\bf x}{\bf
d}^{\prime}+\imath{\bf x}{\bf a}}=\chi({\bf x}{\bf M})e^{\imath({\bf x}{\bf
M}){\bf M}^{-1}({\bf d}^{\prime}+{\bf a})}$ (4)

We represent the linear transformation as a NN layer with two inputs ${\bf x}$
and ${\bf a}$ and two outputs ${\bf x}{\bf M}$ and ${\bf M}^{-1}({\bf
d}^{\prime}+{\bf a})$ (Fig. 2b). By this definition, Eq. (4) is as a two-layer
NN.

Figure 2: (a) A neural network model for the characteristic function. Two
inputs, a data vector ${\bf x}$ with shape $1\times N$ and a bias vector ${\bf
a}$ with shape $N\times 1$ seed the model that compute $\chi$ and returns the
real and imaginary parts of $\chi({\bf x})e^{\imath{\bf x}{\bf a}}$. (b) A
layer representing a linear transformation of the state by a unitary operator
represented by a symplectic $N\times N$ matrix ${\bf M}$ and a displacement
$N\times 1$ vector ${\bf d}^{\prime}$. With such a definition layers can be
cascaded, and one can represent single mode squeezers, interferometers, and
other unitary operators. (c) A model representing a state with characteristic
function $\chi$, subject to a unitary transformation. This is a pullback of a
linear transform from the original state, which produces a new state with
characteristic function $\tilde{\chi}$ [see Eq. (4)].

Figure 2c shows $\tilde{\chi}$ as the “pullback” of the linear layer from the
$\chi$ layer. The two layers form a NN that can be implemented with common
APIs. 111A TensorFlow implementation in a Jupyter notebook is available at
https://github.com/nonlinearxwaves/BosonSampling. Given the vacuum state with
characteristic function $\chi$, one can build the NN model of an arbitrary
state by multiple pullbacks. Indeed, we defined the linear layers in a way
that they can be cascaded. Figure 4a below shows a $n$-mode squeezed vacuum as
a multiple pullback of single mode squezers, each acting on a different mode.

##  III Observables

Observables are computed as derivatives of the NN model. For example, the mean
photon number per mode is related to the derivatives of the characteristic
function. The mean photon number for mode $j$, is

$\langle\hat{n}_{j}\rangle=-\left.\frac{1}{2}\left(\nabla_{j}^{2}+1\right)\chi\right|_{{\bf
x}=0}$ (5)

being $\nabla_{j}^{2}=\partial^{2}_{q_{j}}+\partial^{2}_{p_{j}}$ and
$q_{j}=x_{2j}$ and $p_{j}=x_{2j+1}$. The differential photon number of modes
$j$ and $k$ is

$\langle\left(\hat{n}_{j}-\hat{n}_{k}\right)^{2}\rangle=\left.\left[\frac{1}{4}\left(\nabla_{j}^{2}-\nabla_{k}^{2}\right)^{2}-\frac{1}{2}\right]\chi\right|_{{\bf
x}=0}.$ (6)

Automatic differentiation packages enables an efficient computations of the
derivatives of the NN model.

##  IV Gaussian Boson sampling with the neural network model

In the GBS protocol, one considers a many-body squeezed vacuum state
propagating in an Haar inteferometer, which distributes the photons in the
output modes. For modelling GBS, we hence need squeezing layers and a layer
representing the transmission through random interferometers. The squeezing
layers are realized by a proper design of the corresponding symplectic
matrices ${\bf M}$ with ${\bf d}=0$. We implement the Haar matrix operator by
QuTiP software [37]. Figure 3 shows a pseudo-code to build the neural network
model by composing different layer.

Input = input layer$\triangleright$ Create input layer

V = vacuum layer$\triangleright$ Create vacuum layer

for i=0,1,2,…,n-1 do

S[i] = squeezer layer for mode i$\triangleright$ Create a squeezing layer per
mode

R = random interferometer$\triangleright$ Create a random interferometer

x, a= Input$\triangleright$ Define input tensors x and a

x, a= R (x,a) $\triangleright$ Connect interferometer in pullback order

for i=0,1,2,…,n-1 do

x, a= S[i](x,a) $\triangleright$ Connect squeezing layers in pullback order

x, a= V(x,a) $\triangleright$ Connect vacuum layer in pullback order

Figure 3: Pseudo-code for the creation of a neural network representing a
Gaussian boson sampling experiment Figure 4: (a) A multiple pullback that
represents a many-body squeezed vacuum, obtained by a vacuum state $\chi$ by
cascading $n$ identical single mode squeezers. The order of the squeezers is
not relevant as they act of different modes. $\chi_{R}$ and $\chi_{I}$ are the
real and imaginary part of the resulting characteristic function. (b) GBS
setup, a $n$-body squeezed vacuum enters an Haar inteferometer. Note that the
order of the operators, from the vacuum to the interferometer goes from right
to left. (c) GBS setup including a trainable random interferometer before
entering the Haar interferometer. The multiple squeezers are represented as a
single block. The trainable interferometer can optimize the probability of
pair generation. The right panel shows the architecture of the TensorFlow
model for $n=6$.

Fig. 4b is a graphical representation of the GBS NN model [20]. Boson sampling
corresponds to computing the probability $\Pr(\bar{\bf n})$ of finding
$\bar{n}_{0}$ photons in mode $0$, $\bar{n}_{1}$ photons in mode $1$, and so
forth. $\bar{\bf n}=(\bar{n}_{0},\bar{n}_{1},\ldots,\bar{n}_{n-1})$ is a given
photon pattern. Letting $\hat{\rho}$ the density matrix, one has

$\Pr(\bar{\bf n})=\operatorname{Tr}[\hat{\rho}|\bar{\bf
n}\rangle\langle\bar{\bf n}|],$

with

$|\bar{\bf n}\rangle\langle\bar{\bf
n}|=\otimes_{j=0}^{n-1}|\bar{n}_{j}\rangle\langle\bar{n}_{j}|.$

Correspondingly [30],

$\Pr(\bar{\bf n})=\left.\frac{1}{\bar{{\bf
n}}!}\prod_{j=0}^{n-1}{\left(\frac{\partial^{2}}{\partial\alpha_{j}\partial\alpha_{j}^{*}}\right)}^{\bar{n}_{j}}e^{\sum_{j}|\alpha|_{j}^{2}}Q_{\rho}({\bm{\alpha}},{\bm{\alpha}}^{*})\right|_{{\bm{\alpha}}=0}$
(7)

where $\bar{\bf n}!=\bar{n_{0}}!\bar{n}_{1}!\ldots\bar{n}_{n-1}!$ and

$Q_{\rho}=\pi^{n}\langle{\bm{\alpha}}|\rho|{\bm{\alpha}}\rangle$

is the Q-rapresentation of the density matrix [35, 19] with
${\bm{\alpha}}=\left(\alpha_{0},\alpha_{1},\ldots,\alpha_{n-1}\right)$ complex
displacements.

We introduce the $N\times 1$ real vector ${\bf k}$ as

$\displaystyle k_{2j}$
$\displaystyle=\frac{\alpha_{j}^{*}+\alpha_{j}}{\sqrt{2}}$ $\displaystyle
k_{2j+1}$ $\displaystyle=\frac{\alpha_{j}^{*}-\alpha_{j}^{*}}{\sqrt{2}\imath}$

and we have

$\Pr(\bar{\bf n})=\left.\frac{1}{\bar{{\bf
n}}!2^{\bar{n}_{T}}}\left(\prod_{j}{\tilde{\nabla}_{j}}^{2\bar{n}_{j}}\right)e^{\frac{{\bf
k}^{2}}{2}}Q_{\rho}({\bf k})\right|_{{\bf k}=0}$ (8)

with $\tilde{\nabla}_{j}^{2}=\partial^{2}/\partial
k_{2j}+\partial^{2}/\partial k_{2j+1}$ and
$\bar{n}_{T}=\sum_{j=0}^{n-1}\bar{n}_{j}$. $Q_{\rho}$ in Eq. (8) can be
evaluated explicitly as a multidimensional Gaussian integral:

$\Pr(\bar{\bf n})=\left.\frac{1}{\bar{{\bf
n}}!2^{\bar{n}_{T}}}\left(\prod_{j}\tilde{\nabla}_{j}^{2\bar{n}_{j}}\right)\mathcal{Q}({\bf
k})\right|_{{\bf k}=0}$ (9)

with ($p,q=0,1,..,N-1$)

$\mathcal{Q}({\bf k})=\frac{1}{\sqrt{2^{n}\det A}}e^{\frac{1}{2}{\bf
k}^{2}}e^{-\frac{1}{2}\sum_{pq}A^{-1}_{pq}(k_{p}-d_{q})(k_{p}-d_{q})}$ (10)

being $A_{pq}=\frac{1}{2}\left(g_{pq}+\delta_{pq}\right)$. Eq. (9) and (10)
can be implemented as further layers of the NN, and the probability of a given
pattern computed by running the model. Figure 5a shows an example of the
pattern probability distribution with $n=6$, obtained by using the NN model in
Fig. 4b with squeezing parameters $r_{j}=0.88$ and $\phi_{j}=\pi/4$, such that
all the single mode squeezers are identical, each with mean photon number
$\sinh(r_{j})^{2}\simeq 1$. As in [20], we consider patterns with
$\langle\hat{n}_{j}\rangle=\\{0,1\\}$.

##  V Training Gaussian boson sampling

Our interest is understanding if we can train the model to maximize the
generation of specific patterns, e.g., a photon pair in modes $0$ and $1$.
Using complex media to tailor linear systems is a well renowned technique as,
for example, to synthesize specific gates [38, 39] or taming entanglement
[40]. Here, we use the NN model in the phase space to optimize multi-particle
events.

One could use the squeezing parameters in the model in Fig. 4b as training
parameters. However, the degree of squeezing affects the number of particles
per mode, and we want to alter the statistical properties of states without
changing the average number of particles. We hence consider a GBS setup with
an additional trainable interferometer as in Fig. 4c, which is typically
realized by amplitude or phase modulators.

In Fig. 4c, $n$ squeezed vacuum modes impinge on a trainable interferometer
and then travel through a Haar interferometer. Instead of two distinct
interferometers, one could use a single device (i.e., combine the Haar
interferometer with the trainable interferometer), but we prefer to
distinguish the trainable part from the mode-mixing Haar unitary operator.

Given $n$ modes, our goal is to maximize the probability of patterns that
contains a pair of photons in the mode $0$ or $1$. For example, for $n=6$,
this means maximizing the probability of $\bar{\bf n}=(1,1,0,0,0,0)$ with
respect to $\bar{\bf n}=(1,0,0,1,0,0)$. We use as loss function

$\mathcal{L}=e^{\langle(\hat{n}_{0}-\hat{n}_{1})^{2}\rangle}$ (11)

which is minimal when the expected differential number of photons in mode $0$
and mode $1$ vanishes. This is the case when the state has a particle pair in
mode $0$ and mode $1$. We stress the difference in using other cost functions,
which involve the expected number of photons per mode as, e.g.,

$\mathcal{L}_{0}=e^{{(\langle\hat{n}_{0}\rangle-\langle\hat{n}_{1}\rangle)}^{2}}.$

The linear interferometer does not affect the average number of photons (which
are mixed by the Haar layer). Correspondingly, training using
$\mathcal{L}_{0}$ Eq. V is not be effective to generate entangled pairs. On
the contrary, $\mathcal{L}$ in Eq. (11) contains
$\langle\hat{n}_{0}\hat{n}_{1}\rangle$, which is maximal with a photon pair in
modes $0$ and $1$.

Fig. 5a shows the computed probabilities of pairs for the model in Fig. 4c,
with a random instance of the Haar and the linear inteferometers. Training
strongly alters this statistical distribution, as shown in Fig. 5b.

Figure 5: (a) Probability distribution of patterns with two photons for $n=6$
in the model in Fig. 4c, before training. The insets detail the particle
distribution in the patterns. (b) As in (a) after training, the probability of
finding a pair in mode $0$ and $1$ is enhanced by more than one order of
magnitude. (c) Mean photon number in mode $0$ and $1$ during the training
epochs (green), and expected differential photon number
$\langle\left(\hat{n}_{0}-\hat{n}_{1}\right)^{2}\rangle$ in the two modes,
which vanishes after thousands of epochs. The statistical distribution of
pairs changes at a constant photon number per mode. Data generated by the code
in https://github.com/nonlinearxwaves/BosonSampling.

Fig. 5c shows the trend during the training epochs of
$\langle(\hat{n}_{0}^{2}-\hat{n}_{1}^{2})\rangle$, which goes to zero while
the mean photon numbers $\langle\hat{n}_{0}\rangle$ and
$\langle\hat{n}_{1}\rangle$ remain unaltered.

Training also maximizes higher photon events, as in the pattern $\bar{\bf
n}=(1,1,1,1,0,0)$ with $4$ photons and $n=6$. Fig. 6a shows the pattern
probability with $4$ photons. After training with the loss function in Eq.
(11), $\Pr(\bar{\bf n})$ substantially increases for the patterns with four
photons containing $1$ pair in modes $0$ and $1$ (Fig. 6b).

Figure 6: (a) Probability distribution of patterns with $4$ photons ($n=6$) in
the model in Fig. 4c before training. The insets detail the particles in each
pattern. (b) As in (a) after training; the probability of patterns with two
photons in modes $0$ and $1$ is maximized. Data generated by the code in
https://github.com/nonlinearxwaves/BosonSampling.

##  VI Conclusions

We have shown that a many-body characteristic function may be reformulated as
a layered neural network. This approach enables to build complex states for
various applications, as gate design or boson sampling.

A common argument in criticizing quantum neural networks is that the linear
quantum mechanics does not match with the nonlinearity-eager NN models.
However, recent investigations show that nonlinearity may be introduced in
quantum neural networks [43]. Our remark is that if we formulate quantum
mechanics in the phase space, nonlinearity arises in the characteristic
function (or other representation). We analyzed this strategy in the simplest
case of Gaussian states. The resulting model is universal and may be trained
for different purposes. For this reason, phase space models allow naturally in
dealing with non-classical states and computing observables by derivatives.
This formulation opens many opportunities. For example, the optimization of
multi-particle events can be extended to fermionic fields. As a drawback,
computing boson patterns probabilities by NN APIs is not expected to be
competitive with highly optimized algorithms running on large-scale clusters
[41, 42]. Still, it appears to be a versatile and straightforward methodology.

Here, we have shown many-body quantum state design and engineering by
TensorFlow. We have demonstrated how to enhance multi-particle generation,
with many potential applications in quantum technologies. In addition, the
proposed method enables training Boson sampling without explicitly computing
derivatives of the Hafnian [34, 18], but resorting to automatic computational
packages. We have tested the algorithm with a conventional workstation with a
single commercial GPU (NVIDIA QUADRO RTX 4000), with a computational time of
the order of few minutes with 6 modes.

The method can be generalized to other boson sampling setups, as including
Glauber layers and multi-mode squeezers. Also, it readily allows to test
different loss functions for tailoring the boson sampling patterns. Extension
beyond Gaussian states can be envisaged by using a general machine learning
networks with an arbitrary number of layers and different nonlinearity.

###### Acknowledgements.

We acknowledge support from Horizon 2020 Framework Programme QuantERA grant
QUOMPLEX, by National Research Council (CNR), Grant 731473.

## References

  * [1] H.-Y. Huang, R. Kueng, G. Torlai, V. V. Albert, J. Preskill, arXiv:2106.1267 (2021).
  * [2] G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, L. Zdeborová, Rev. Mod. Phys. 91, 045002 (2019).
  * [3] F. Vicentini, A. Biella, N. Regnault, C. Ciuti, Phys. Rev. Lett. 122, 250503 (2019)
  * [4] S. Mangini, F. Tacchino, D. Gerace, D. Bajoni, C. Macchiavello, ArXiv:2102.03879 (2021)
  * [5] F. Marquardt, ArXiv:2101.01759 (2021)
  * [6] A. Lumino, E. Polino, A.S. Rab, G. Milani, N. Spagnolo, N. Wiebe, F. Sciarrino, Phys. Rev. Applied 10, 044033 (2018).
  * [7] A. Fratalocchi, A. Fleming, C. Conti, A.D. Falco, Nanophotonics 10, 457 (2021).
  * [8] G. Marcucci, D. Pierangeli, P.W.H. Pinkse, M. Malik, C. Conti, Opt. Express 28, 14018 (2020).
  * [9] A.A. Melnikov, P. Sekatski, N. Sangouard, Phys. Rev. Lett. 125, 160401 (2020).
  * [10] G. Marcucci, D. Pierangeli, C. Conti, Phys. Rev. Lett. 125, 093901 (2020)
  * [11] T.W. Hughes, I.A. Williamson, M. Minkov, S. Fan, Sci. Adv. 5, eaay6946 (2019)
  * [12] D. Ballarini, A. Gianfrate, R. Panico, A. Opala, S. Ghosh, L. Dominici, V. Ardizzone, M.D. Giorgi, G. Lerario, G. Gigli, T.C.H. Liew, M. Matuszewski, D. Sanvitto, Nano Lett. 20, 3506 (2020)
  * [13] J. Nokkala, R. Martínez-Peña, G.L. Giorgi, V. Parigi, M.C. Soriano, R. Zambrini, ArXiv:2006.04821 (2020)
  * [14] D. Marković, J. Grollier, Appl. Phys. Lett. 117, 150501 (2020).
  * [15] N.A. Silva, T.D. Ferreira, A. Guerreiro, New J. Phys. 23, 023013 (2021).
  * [16] P. Sgroi, G.M. Palma, M. Paternostro, Phys. Rev. Lett. 126, 026601 (2021)
  * [17] L. Pilozzi, F.A. Farrelly, G. Marcucci, C. Conti, Nanotechnology 32, 142001 (2021).
  * [18] M. Broughton, G. Verdon, T. McCourt, A.J. Martinez, J.H. Yoo, S.V. Isakov, P. Massey, M.Y. Niu, R. Halavati, E. Peters, M. Leib, A. Skolik, M. Streif, D.V. Dollen, J.R. McClean, S. Boixo, D. Bacon, A.K. Ho, H. Neven, M. Mohseni, ArXiv:2003.02989 (2020)
  * [19] S.M. Barnett, P.M. Radmore, _Methods in Theoretical Quantum Optics_ (Oxford University Press, New York, 1997)
  * [20] C.S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, I. Jex, Phys. Rev. Lett. 119, 170501 (2017).
  * [21] N. Quesada, J.M. Arrazola, N. Killoran, Phys. Rev. A 98, 062322 (2018).
  * [22] H.S. Zhong, H. Wang, Y.H. Deng, M.C. Chen, L.C. Peng, Y.H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, P. Hu, X.Y. Yang, W.J. Zhang, H. Li, Y. Li, X. Jiang, L. Gan, G. Yang, L. You, Z. Wang, L. Li, N.L. Liu, C.Y. Lu, J.W. Pan, Science 370, 1460 (2020)
  * [23] M. Tillmann, B. Dakić, R. Heilmann, S. Nolte, A. Szameit, P. Walther, Nat. Photonics 7, 540 (2012)
  * [24] M.A. Broome, A. Fedrizzi, S. Rahimi-Keshari, J. Dove, S. Aaronson, T.C. Ralph, A.G. White, Science 339, 794 (2013)
  * [25] J.B. Spring, B.J. Metcalf, P.C. Humphreys, W.S. Kolthammer, X.M. Jin, M. Barbieri, A. Datta, N. Thomas-Peter, N.K. Langford, D. Kundys, J.C. Gates, B.J. Smith, P.G.R. Smith, I.A. Walmsley, Science 339, 798 (2013)
  * [26] N. Spagnolo, C. Vitelli, M. Bentivegna, D.J. Brod, A. Crespi, F. Flamini, S. Giacomini, G. Milani, R. Ramponi, P. Mataloni, R. Osellame, E.F. Galvão, F. Sciarrino, Nat. Photonics 8, 615 (2014).
  * [27] J. Carolan, J.D.A. Meinecke, P.J. Shadbolt, N.J. Russell, N. Ismail, K. Wörhoff, T. Rudolph, M.G. Thompson, J.L. O’Brien, J.C.F. Matthews, A. Laing, Nat. Photonics 8, 621 (2014).
  * [28] H. Wang, J. Qin, X. Ding, M.C. Chen, S. Chen, X. You, Y.M. He, X. Jiang, L. You, Z. Wang, C. Schneider, J.J. Renema, S. Höfling, C.Y. Lu, J.W. Pan, Phys. Rev. Lett. 123, 250503 (2019).
  * [29] S. Aaronson, A. Arkhipov, Theory Comput. 9, 143 (2013).
  * [30] R. Kruse, C.S. Hamilton, L. Sansoni, S. Barkhofen, C. Silberhorn, I. Jex, Phys. Rev. A 100, 032326 (2019)
  * [31] J.M. Arrazola et al., Nature 591, 54 (2021)
  * [32] F. Hoch et al., arXiv:2106.08260 (2021)
  * [33] H.-S. Zhong et al., arXiv:2106.15534 (2021)
  * [34] L. Banchi, N. Quesada, J.M. Arrazola, Phys. Rev. A 102, 012417 (2020)
  * [35] C.W. Gardiner, P. Zoller, _Quantum Noise_ , 3rd edn. (Springer-Verlag, Berlin, 2004)
  * [36] X. Wang, T. Hiroshima, A. Tomita, M. Hayashi, Phys. Rep. 448, 1 (2007).
  * [37] J. Johansson, P. Nation, F. Nori, Comput. Phys. Commun. 184, 1234 (2013).
  * [38] S. Leedumrongwatthanakun, L. Innocenti, H. Defienne, T. Juffmann, A. Ferraro, M. Paternostro, S. Gigan, Nat. Photonics 14, 139 (2020)
  * [39] C. Taballione, T.A.W. Wolterink, J. Lugani, A. Eckstein, B.A. Bell, R. Grootjans, I. Visscher, D. Geskus, C.G.H. Roeloffzen, J.J. Renema, I.A. Walmsley, P.W.H. Pinkse, K.J. Boller, Opt. Express 27, 26842 (2019).
  * [40] N.H. Valencia, S. Goel, W. McCutcheon, H. Defienne, M. Malik, Nat. Phys. 16, 1112 (2020)
  * [41] N. Quesada, J.M. Arrazola, Phys. Rev. Research 2, 023005 (2020).
  * [42] Y. Li, M. Chen, Y. Chen, H. Lu, L. Gan, C. Lu, J. Pan, H. Fu, G. Yang, ArXiv:2009.01177 (2020)
  * [43] C. Zhao, X.-S. Gao, Quantum Mach. Intell. 3, 15 (2021)