#  Variational quantum algorithm for Gaussian discrete solitons and their
boson sampling

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-CNR), Via dei Taurini 19, 00185 Rome, Italy
Research Center Enrico Fermi (CREF), Via Panisperna 89a, 00184 Rome, Italy
https://www.complexlight.org

###### Abstract

In the context of quantum information, highly nonlinear regimes, such as those
supporting solitons, are marginally investigated. We miss general methods for
quantum solitons, although they can act as entanglement generators or as self-
organized quantum processors. We develop a computational approach that uses a
neural network as a variational ansatz for quantum solitons in an array of
waveguides. By training the resulting phase space quantum machine learning
model, we find different soliton solutions varying the number of particles and
interaction strength. We consider Gaussian states that enable measuring the
degree of entanglement and sampling the probability distribution of many-
particle events. We also determine the probability of generating particle
pairs and unveil that soliton bound states emit correlated pairs. These
results may have a role in boson sampling with nonlinear systems and in
quantum processors for entangled nonlinear waves.

A soliton is a non-perturbative solution of a classical nonlinear wave-
equation; it may describe mean-field states of atoms (as in Bose-Einstein
condensation) or photons (as in nonlinear optics) [1]. From a quantum
mechanical perspective, a soliton may correspond to a coherent state; however,
the nonlinearity may induce squeezing or non-Gaussianity [2]. The quantum
properties of solitons inspired experimental investigations, as quantum non-
demolition, squeezing [3, 4, 5, 6] and photon bound states [7]. Authors
reported on theoretical studies on the soliton quantum features, as
evaporation and breathing [8, 9, 10, 11, 12, 13]. Following these
investigations, one can argue that solitons are new tools for fundamental
many-body phenomena and sources of highly non-classical states.

However, solitons are overlooked in quantum technologies, as they require
specialized theoretical methods [14, 15, 16, 17, 18] that do not appear
compatible with modern quantum information. On the other hand, quantum optical
protocols mainly employ linear circuits, and we ask if nonlinear waves may
have a role in quantum computing [19, 20, 21, 22, 23], or –reciprocally– if
quantum processors can generate quantum solitons. An intriguing opportunity
arises from the fact that linear circuits for the quantum advantage in boson
sampling (BS) [24, 25, 26, 27, 28, 29] adopt the same devices used to observe
discrete solitons [30, 31, 1]. BS in a linear optical circuit is
“intermediate” between a many-body system not performing computation and a
universal quantum computer [32]. To make a conceptual step, we can consider
nonlinear phenomena in BS. A soliton is a self-organized nonlinear channel. Is
the probability of specific events affected by these nonlinear channels? Can
we do BS with solitons?

To merge solitons and quantum information, we need novel methodologies for
quantum protocols on many-body states bounded by nonlinearity. Recently,
neural networks (NNs) have been introduced to study many-body models and
quantum circuits [33]. Quantum machine learning (QML), i.e., the development
and training of parameterized quantum circuits [34, 35, 36, 37, 38], is a
technique to approximate high-dimensional ground states [39]. In QML, one
considers quantum circuits depending on trainable parameters, such as rotation
angles for qubits or degree of squeezing for continuous variables [40].
Applications include quantum classifiers [41, 42], or variational quantum
algorithms [43, 44] for optimization and quantum simulations.

Here we develop a quantum variational ansatz in the phase space for quantum
soliton states. The ansatz corresponds to a trainable BS quantum processor
[45, 46] and can synthesize various soliton solutions. The trained NN is used
to compute their quantum features, showing that soliton bound-states exhibit
entanglement. Also, the quantum variational algorithm is used to sample the
probability distribution at specific patterns of output particles (i.e., BS)
that reveal the generation of correlated bosons. The approach mixing phase
space methods and machine learning turns out to be a remarkable tool for
studying quantum solitons, and to predict and quantify novel effects not
accessible with conventional techniques.

This manuscript is organized as follows. In Sections I and II, we review the
basics of phase space representations, characteristic functions, and outline
the link with conventional NN models. In sections from III to IX, we show how
a Gaussian state can be reinterpreted as a stack of linear layers with an
Gaussian activation function; we discuss how to map different quantum
operators to layers in a NN. Then, in sections X and XI, we represent the
Gaussian ansatz for a quantum soliton as a quantum processor with trainable
parameters and compute different soliton solutions. Finally, in sections XII
and XIII, we use the trained neural network for various quantum features, as
entanglement and BS. Conclusions are drawn in Section XIV.

##  I Quantum phase space and neural networks

In the phase space, a state is represented by a function of several variables,
which we encode in a real vector ${\bf x}$. For a $n$-body system, ${\bf x}$
is a vector of $N=2\,n$ real variables. Among the many representations [47],
we consider the characteristic function $\chi({\bf x})$.

For Gaussian states [48], in symmetric ordering,

$\chi({\bf x})=\exp\left(-\frac{1}{4}{\bf x}\cdot\bm{g}\cdot{\bf
x}+\imath\,{\bf x}\cdot\bm{d}\right)\;,$ (1)

with $\bm{g}$ the $N\times N$ covariance matrix, and $\bm{d}$ the $N\times 1$
displacement vector. $\chi({\bf x})$ is a complex function $\chi({\bf
x})=\chi_{R}({\bf x})+\imath\chi_{I}({\bf x})$ with real part $\chi_{R}$ and
imaginary part $\chi_{I}$. Any many-body state is represented by a couple of
real functions $\chi_{R,I}$ of $N$ real variables. As NN models typically deal
with real-valued quantities, in the following we will represent the
characteristic functions by models with two real outputs corresponding to
$\chi_{R}$ and $\chi_{I}$. The resulting NN model includes linear
transformations representing gates, as displacements or interferometers,
followed by a nonlinear activation, which corresponds to computing the
characteristic functions (e.g, a Gaussian function). This maps the quantum
state in a conventional NN architecture that can be trained by well-known
algorithms.

##  II The characteristic function

We introduce the complex vector $\bm{z}$ with components

$\bm{z}=(z_{0},z_{0},\ldots,z_{n-1},z_{0}^{*},z_{1}^{*},\ldots,z_{n-1}^{*})\;,$
(2)

which include the complex $z_{j}$ components and their conjugates $z_{j}^{*}$.
Given the density matrix $\rho$, the characteristic function is [47]

$\chi(\bm{z},\mathcal{P})=\text{Tr}\left\\{\rho
e^{\sum_{k}\left(z_{k}\hat{a}^{\dagger}_{k}-z^{*}_{k}\hat{a}_{k}\right)}\right\\}e^{\sum_{k}\frac{\mathcal{P}}{2}z_{k}z_{k}^{*}}\;,$
(3)

In Eq. (3) the variable $\mathcal{P}\in\\{0,1,-1\\}$, refers to the adopted
operator ordering. The characteristic function depends on $\bm{z}$ and the
ordering index $\mathcal{P}$. $\mathcal{P}$ is important when determining the
expected values as derivatives of $\chi$.

For an observable $\hat{O}$ we have

$\langle\hat{O}\rangle=\text{Tr}\left(\rho\hat{O}\right).$ (4)

The expected value of a combination of the annihilation and creation operators
is

$\langle{\left(\hat{a}_{j}^{\dagger}\right)}^{m}{\left(\hat{a}_{k}\right)}^{n}\rangle_{\mathcal{P}}=\left.{\left(\frac{\partial}{\partial
z_{j}}\right)}^{m}{\left(-\frac{\partial}{\partial
z_{k}^{*}}\right)}^{n}\chi(\bm{z},\mathcal{P})\right|_{\bm{z}=0}\;,$ (5)

with $m$ and $n$ integers.

The mean value of the field operator

$\langle\hat{a}_{k}\rangle=-\frac{\partial\chi}{\partial z_{k}^{*}},$ (6)

does not depend on the ordering index $\mathcal{P}$, as also

$\langle\hat{a}^{\dagger}_{k}\rangle=\frac{\partial\chi}{\partial z_{k}}.$ (7)

The mean particle number of mode $k$,
$\hat{n}_{k}=\hat{a}_{k}^{\dagger}\hat{a}_{k}$ is

$\langle\hat{a}_{k}^{\dagger}\hat{a}_{k}\rangle_{\mathcal{P}}=\left.-\frac{\partial^{2}\chi(\bm{z},\mathcal{P})}{\partial
z_{k}\partial z_{k}^{*}}\right|_{\bm{z}=0}\;.$ (8)

For the symmetric ordering $\mathcal{P}=0$, one has

$\langle\hat{a}_{k}^{\dagger}\hat{a}_{k}\rangle_{\mathcal{P}=0}=\langle\hat{a}_{k}^{\dagger}\hat{a}_{k}+\hat{a}_{k}\hat{a}_{k}^{\dagger}\rangle;$
(9)

for the normal ordering $\mathcal{P}=1$, one has

$\langle\hat{a}_{k}^{\dagger}\hat{a}_{k}\rangle_{\mathcal{P}=1}=\langle\hat{a}_{k}^{\dagger}\hat{a}_{k}\rangle,$
(10)

finally for antinormal ordering $\mathcal{P}=-1$,

$\langle\hat{a}_{k}^{\dagger}\hat{a}_{k}\rangle_{\mathcal{P}=-1}=\langle\hat{a}_{k}\hat{a}_{k}^{\dagger}\rangle\;.$
(11)

In general,

$\langle\hat{a}_{k}^{\dagger}\hat{a}_{k}\rangle_{\mathcal{P}}=\langle\hat{a}_{k}^{\dagger}\hat{a}_{k}\rangle+\frac{1}{2}(1-\mathcal{P}),$
(12)

which is adopted below for the NN layer that returns the mean particle number.
When not explicitly stated, we refer to symmetric ordering $\mathcal{P}=0$.

For use with the machine learning application programming interfaces as
TensorFlow, it is convenient to use real variables. We consider the quadrature
operators $\hat{x}_{k}$ and $\hat{p}_{k}$ defined by

$\hat{a}_{k}=\frac{\hat{q}_{k}+\imath\hat{p}^{\dagger}_{k}}{\sqrt{2}}.$ (13)

The quadrature operators can be cast in operator vectors $\bm{\hat{q}}$ and
$\bm{\hat{p}}$ with dimension $n\times 1$,

$\bm{\hat{q}}=\begin{pmatrix}\hat{q}_{0}\\\ \hat{q}_{1}\\\ \vdots\\\
\hat{q}_{n-1}\end{pmatrix}\hskip
28.45274pt\bm{\hat{p}}=\begin{pmatrix}\hat{p}_{0}\\\ \hat{p}_{1}\\\ \vdots\\\
\hat{p}_{n-1}\end{pmatrix}.$ (14)

In addition, we consider the real quantities

$q_{k}=\frac{z_{k}-z_{k}^{*}}{\sqrt{2}\imath}\hskip
28.45274ptp_{k}=\frac{z_{k}+z_{k}^{*}}{\sqrt{2}},$ (15)

and we collect them in two real vectors $\bm{q}$ and $\bm{p}$ with dimension
$1\times n$

$\begin{array}[]{l}\bm{q}=\begin{pmatrix}q_{0}&q_{1}&\cdots&q_{n-1}\end{pmatrix}\\\
\bm{p}=\begin{pmatrix}p_{0}&p_{1}&\cdots&p_{n-1}\end{pmatrix}.\end{array}$
(16)

Correspondingly, we have

$\chi=\displaystyle\chi(\bm{q},\bm{p},\mathcal{P}=0)=\text{Tr}[\rho\exp\left(\imath\,\bm{q}\cdot\hat{\bm{q}}+\imath\,\bm{p}\cdot\hat{\bm{p}}\right)]\text{.}$
(17)

As $\hat{\bm{q}}$ and $\hat{\bm{p}}$ are column vectors, and $\bm{q}$ and
$\bm{p}$ are row vectors, we omit the dot product symbol and write

$\displaystyle\chi(\bm{q},\bm{p})=\text{Tr}[\rho\exp\left(\imath\,\bm{q}\,\hat{\bm{q}}+\imath\,\bm{p}\,\hat{\bm{p}}\right)]\;.$
(18)

The notation can be simplified by defining the $1\times N$ real row vector

${\bf x}=\begin{pmatrix}q_{0}&p_{0}&q_{1}&p_{2}&\cdots
q_{n}&p_{n}\end{pmatrix},$ (19)

and we have

$\chi=\chi({\bf x})=\text{Tr}\left[\rho\exp\left(\imath\,{\bf
x}\,\hat{\bm{\mathrm{R}}}\right)\right]=\chi_{R}({\bf x})+\imath\chi_{I}({\bf
x})$ (20)

being

$\displaystyle{\bf x}\,\hat{\bm{\mathrm{R}}}$ $\displaystyle=$ (21)
$\displaystyle=q_{0}\hat{q}_{0}+p_{0}\hat{p}_{0}+q_{1}\hat{q}_{1}+\ldots+q_{n-1}\hat{q}_{n-1}+p_{n-1}\hat{p}_{n-1}=$
$\displaystyle=x_{0}\hat{R}_{0}+x_{1}\hat{R}_{1}+\ldots+x_{N-1}\hat{R}_{N-1}$

with the $N\times 1$ column operator vector

$\hat{\bm{\mathrm{R}}}=\begin{pmatrix}\hat{q}_{0}\\\ \hat{p}_{0}\\\
\hat{q}_{1}\\\ \hat{p}_{1}\\\ \vdots\\\ \hat{q}_{n-1}\\\
\hat{p}_{n-1}\end{pmatrix}\;.$ (22)

Following the canonical commutation relations of the quadratures, in units
with $\hbar=1$,

$\displaystyle\left[\hat{q_{j}},\hat{q}_{k}\right]$ $\displaystyle=0$ (23)
$\displaystyle\left[\hat{p_{j}},\hat{p}_{k}\right]$ $\displaystyle=0$
$\displaystyle\left[\hat{q_{j}},\hat{p}_{k}\right]$
$\displaystyle=i\delta_{jk},$

we have

$[\hat{R}_{p},\hat{R}_{q}]=\imath J_{pq},$ (24)

where we introduced the $N\times N$ symplectic matrix $\bm{\mathrm{J}}$

$\bm{\mathrm{J}}=\bigoplus_{j}\bm{\mathrm{J}}_{1}=\begin{pmatrix}0&1&0&0&\cdots&0&0\\\
-1&0&0&0&\cdots&0&0\\\ 0&0&0&1&\cdots&0&0\\\ 0&0&-1&0&\cdots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&0&0&0&\cdots&0&1\\\
0&0&0&0&\cdots&-1&0\end{pmatrix},$ (25)

being $\bm{\mathrm{J}}_{1}=\big{(}\begin{smallmatrix}0&1\\\
-1&0\end{smallmatrix}\big{)}$ [48]. In Eq. (25), the symbol $\bigoplus$ is the
direct sum for block matrices with $j$ running in $0,1,\ldots\,n-1$. The
matrix $\bm{\mathrm{J}}$ is such that
$\bm{\mathrm{J}}^{-1}=\bm{\mathrm{J}}^{\top}$, $\bm{\mathrm{J}}^{2}=-{\bf
1}_{N}$, and $\bm{\mathrm{J}}^{\top}\bm{\mathrm{J}}={\bf 1}_{N}$ is the
$N\times N$ identity matrix.

The expected value of $\bm{\hat{R}}$ is determined by the derivatives of the
characteristic function,

$\chi({\bf x})=\text{Tr}\left(\rho e^{\imath\,{\bf x}\,\hat{\bf R}}\right)\;.$
(26)

One has

$\langle\hat{R}_{j}\rangle=\text{Tr}[\rho\hat{R}_{j}]=-\imath\left.\frac{\partial\chi}{\partial
x_{j}}\right|_{{\bf x}=0}=\left.\frac{\partial\chi_{I}}{\partial
x_{j}}\right|_{{\bf x}=0}$ (27)

or

$\langle\hat{\bm{\mathrm{R}}}\rangle=\text{Tr}[\rho\hat{\bm{\mathrm{R}}}]=-\imath\left.\nabla\chi\right|_{{\bf
x}=0}=\left.\nabla\chi_{I}\right|_{{\bf x}=0}\;.$ (28)

As the components of ${\hat{\bf R}}$ are self-adjoint, their mean value is
real, and it is given by the derivative of the imaginary part $\chi_{I}$ at
${\bf x}=0$. We also have

$\displaystyle-\frac{\partial^{2}}{\partial z_{j}\partial z_{j}^{*}}\chi=$
(29) $\displaystyle=-\left(\frac{\partial}{\partial q_{j}}\frac{\partial
q_{j}}{\partial z_{j}}+\frac{\partial}{\partial p_{j}}\frac{\partial
p_{j}}{\partial z_{j}}\right)\left(\frac{\partial}{\partial
q_{j}}\frac{\partial q_{j}}{\partial z^{*}_{j}}+\frac{\partial}{\partial
p_{j}}\frac{\partial p_{j}}{\partial z^{*}_{j}}\right)\chi$
$\displaystyle=-\frac{1}{2}\left(\frac{\partial^{2}}{\partial
q_{j}^{2}}+\frac{\partial^{2}}{\partial p_{j}^{2}}\right)\chi,$

after using Eq. (15). Thus, we evaluate the expected number of photons by the
second derivatives of $\chi_{R}$

$\displaystyle\langle\hat{a}^{\dagger}_{j}\hat{a_{j}}\rangle_{\mathcal{P}=0}=\left.-\frac{1}{2}\left(\frac{\partial^{2}}{\partial
x_{2j}^{2}}+\frac{\partial^{2}}{\partial
x_{2j+1}^{2}}\right)\chi_{R}\right|_{{\bf x}=0}$ (30)

with $j=0,1,2,\ldots,n-1$.

For the expected value of total number of particles $\hat{\mathcal{N}}$, we
have

$\displaystyle\langle\hat{\mathcal{N}}\rangle_{\mathcal{P}=0}=\left.-\frac{1}{2}\nabla^{2}\chi_{R}\right|_{{\bf
x}=0}$ (31)

the Laplacian being computed with respect to all the components of ${\bf x}$.
Equation (31) refers to symmetric ordering with $\mathcal{P}=0$. For normal
ordering, we use Eq. (12) as

$\langle\hat{a}_{k}^{\dagger}\hat{a}_{k}\rangle=\langle\hat{a}_{k}^{\dagger}\hat{a}_{k}\rangle_{\mathcal{P}=1}=\langle\hat{a}_{k}^{\dagger}\hat{a}_{k}\rangle_{\mathcal{P}=0}-\frac{1}{2}\;,$
(32)

and we get

$\begin{array}[]{ll}\displaystyle\langle\hat{\mathcal{N}}\rangle&=\displaystyle\langle\hat{\mathcal{N}}\rangle_{\mathcal{P}=1}=\displaystyle\langle\hat{\mathcal{N}}\rangle_{\mathcal{P}=0}-\frac{n}{2}\\\\[2.0pt]
&=\left.-\frac{1}{2}\nabla^{2}\chi_{R}\right|_{{\bf
x}=0}-\frac{n}{2}\\\\[2.0pt]
&=-\frac{1}{2}\left.\left(\nabla^{2}+n\right)\chi\right|_{{\bf
x}=0}.\end{array}$ (33)

##  III Gaussian states

For a Gaussian state,

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

with $\bm{g}$ the $N\times N$ covariance matrix, and $\bm{d}$ the $N\times 1$
displacement column vector. Correspondingly, being ${{\bf x}}$ a row vector,

$\begin{array}[]{l}\chi_{R}({\bf x})\displaystyle=e^{-\frac{{\bf x}\bm{g}{\bf
x}^{\top}}{4}}\cos({\bf x}\bm{d})\\\ \chi_{I}({\bf
x})\displaystyle=e^{-\frac{{\bf x}\bm{g}{\bf x}^{\top}}{4}}\sin({\bf
x}\bm{d})\end{array}$ (35)

From Eq. (28), we have for $\bm{d}$, with $q=0,1,\ldots,N-1$.

$\langle\hat{R}_{q}\rangle=\text{Tr}[\hat{\rho}\hat{R}_{q}]=d_{q}=\left.\frac{\partial\chi_{I}}{\partial
x_{q}}\right|_{{\bf x}=0}$ (36)

or

$\langle\bm{\mathrm{R}}\rangle=\bm{d}=\left.\nabla\chi_{I}\right|_{{\bf
x}=0}\;.$ (37)

Hence, the displacement $d_{q}$ is the first moment of the Gaussian $\chi$.
The derivatives of the characteristic function can be expressed by the
displacement and covariance matrix $d_{q}$ and $g_{pq}$. This speeds up
computing the derivatives of $\chi$, and the observable quantities.

For the particle number, we have

$\displaystyle\langle\hat{a}^{\dagger}_{j}\hat{a_{j}}\rangle_{\mathcal{P}=0}$
$\displaystyle=\left.-\frac{1}{2}\left(\frac{\partial^{2}}{\partial
x_{2j}^{2}}+\frac{\partial^{2}}{\partial x_{2j+1}^{2}}\right)\right|_{{\bf
x}=0}\chi_{R}$ (38)
$\displaystyle=\frac{g_{2j,2j}+g_{2j+1,2j+1}}{4}+\frac{d_{2j}^{2}+d_{2j+1}^{2}}{2}$

and

$\displaystyle\langle\hat{a}^{\dagger}_{j}\hat{a_{j}}\rangle$
$\displaystyle=-\frac{1}{2}\left(\frac{\partial^{2}}{\partial
x_{2j}^{2}}+\frac{\partial^{2}}{\partial
x_{2j+1}^{2}}\right)\chi_{R}-\frac{1}{2}$ (39)
$\displaystyle=\frac{g_{2j,2j}+g_{2j+1,2j+1}}{4}+\frac{d_{2j}^{2}+d_{2j+1}^{2}}{2}-\frac{1}{2}$

while the second derivatives of $\chi_{I}$ at ${\bf x}=0$ are vanishing, and
we used Eq. (12). The covariance matrix $\bm{g}$ is determined by the second
moments

$g_{pq}=\langle\\{\hat{R}_{p},\hat{R}_{q}\\}\rangle-2\langle{\hat{R}_{p}}\rangle\langle{\hat{R}_{q}}\rangle$
(40)

with

$\\{\hat{R}_{p},\hat{R}_{q}\\}=\hat{R}_{p}\hat{R}_{q}+\hat{R}_{q}\hat{R}_{p}.$
(41)

We also have after Eq. (24)

$\displaystyle g_{pq}$
$\displaystyle=\langle\hat{R}_{p}\hat{R}_{q}\rangle+\langle\hat{R}_{q}\hat{R}_{p}\rangle-2\langle{R_{p}}\rangle\langle{R_{q}}\rangle$
(42)
$\displaystyle=\langle\hat{R}_{p}\hat{R}_{q}\rangle+\langle\hat{R}_{q}\hat{R}_{p}\rangle-2\langle{R_{p}}\rangle\langle{R_{q}}\rangle$
$\displaystyle=2\langle\hat{R}_{p}\hat{R}_{q}\rangle-\imath
J_{pq}-2\langle{R_{p}}\rangle\langle{R_{q}}\rangle=$
$\displaystyle=2\langle\left(\hat{R}_{p}-d_{p}\right)\left(\hat{R}_{q}-d_{q}\right)\rangle-\imath
J_{pq}\;.$

###  III.1 Vacuum state

Vacuum states are Gaussian state with $\bm{d}={\bf 0}$ and $\bm{g}={\bf
1}_{N}$. From Eq. (39), we get

$\langle\hat{a}^{\dagger}_{j}\hat{a_{j}}\rangle_{\mathcal{P}=0}=\frac{1}{2}\;.$
(43)

Correspondingly,

$\langle\hat{a}^{\dagger}_{j}\hat{a_{j}}\rangle=\langle\hat{a}^{\dagger}_{j}\hat{a_{j}}\rangle_{\mathcal{P}=0}-\frac{1}{2}=0.$
(44)

We have $\langle\hat{\bm{\mathrm{R}}}\rangle=0$, and for the average photon
number after Eq. (39)

$\langle\mathcal{N}\rangle=\sum_{j=0}^{n-1}\langle\hat{a}^{\dagger}_{j}\hat{a}_{j}\rangle=0\;.$
(45)

###  III.2 Coherent state

A coherent state has a non-vanishing displacement vector $\bm{d}\neq{\bf 0}$
and $\bm{g}={\bf 1}_{N}$. For a single mode $|\alpha\rangle$, with $n=1$ and
$N=2$,

$\begin{array}[]{l}d_{0}=\sqrt{2}\Re\left({\alpha}\right)\displaystyle=\frac{\alpha+\alpha^{*}}{\sqrt{2}}\\\
d_{1}=\sqrt{2}\Im\left({\alpha}\right)\displaystyle=\frac{\alpha-\alpha^{*}}{\imath\sqrt{2}}\text{.}\end{array}$
(46)

From Eq. (39), we have ($j=0$)

$\langle\hat{a}^{\dagger}_{j}\hat{a_{j}}\rangle=\langle\hat{a}^{\dagger}_{j}\hat{a_{j}}\rangle_{\mathcal{P}=0}-\frac{1}{2}=|\alpha|^{2}=\frac{d_{0}^{2}+d_{1}^{2}}{2}.$
(47)

For $n$ modes, $j=0,1,\ldots,n-1$ we have

$\langle\hat{a}^{\dagger}_{j}\hat{a_{j}}\rangle=\frac{d_{2j}^{2}+d_{2j+1}^{2}}{2},$
(48)

and

$\langle\mathcal{N}\rangle=\sum_{j=0}^{n-1}\langle\hat{a}^{\dagger}_{j}\hat{a}_{j}\rangle=\sum_{q=0}^{N-1}{\frac{d^{2}_{q}}{2}}.$
(49)

###  III.3 Covariance matrix

Given $\chi$, one can compute the covariance matrix and the displacement
operator by derivation. $\bm{d}$ can be computed by the first derivatives, as
in Eq. (28)

$d_{q}=\langle\hat{R}_{q}\rangle=\text{Tr}[\hat{\rho}\hat{R}_{q}]=\left.\frac{\partial\chi_{I}}{\partial
x_{q}}\right|_{{\bf x}=0}.$ (50)

Writing the characteristic function in terms of the components of ${\bf x}$

$\chi({\bf
x})=\exp\left(-\frac{1}{4}\sum_{pq}g_{pq}x_{p}x_{q}+\imath\sum_{p}x_{p}d_{p}\right)\;.$
(51)

with $p,q=0,1,...,N-1$, we have for the first derivative, by using the
symmetry $g_{pq}=g_{qp}$,

$\frac{\partial{\chi}}{{\partial
x_{m}}}=\left(-\frac{1}{2}\sum_{p}g_{mp}x_{p}+\imath d_{m}\right)\chi({\bf
x})$ (52)

with $m=0,1,...,N-1$. Evaluating Eq. (52) at ${\bf x}=0$, we obtain as above,
being $\chi({\bf 0})=1$,

$\left.\frac{\partial{\chi}}{{\partial x_{m}}}\right|_{{\bf
x}=0}=\imath\left.\frac{\partial{\chi_{I}}}{{\partial x_{m}}}\right|_{{\bf
x}=0}=\imath d_{m}\;.$ (53)

For the second derivatives, we have

$\begin{array}[]{l}\displaystyle\left.\frac{\partial^{2}{\chi_{R}}}{{\partial
x_{m}}{\partial x_{n}}}\right|_{{\bf x}=0}=-\frac{1}{2}g_{mn}-d_{m}d_{n}\\\
\displaystyle\left.\frac{\partial^{2}{\chi_{I}}}{{\partial x_{m}}{\partial
x_{n}}}\right|_{{\bf x}=0}\textbf{}=0\;.\\\ \end{array}$ (54)

The covariance matrix is given by

$\displaystyle g_{pq}=-2\left.\frac{\partial^{2}\chi_{R}}{\partial
x_{p}\partial x_{q}}\right|_{{\bf x}=0}-2d_{p}d_{q}$ (55)
$\displaystyle=\left.-2\frac{\partial^{2}\chi_{R}}{\partial x_{p}\partial
x_{q}}-2\frac{\partial\chi_{I}}{\partial
x_{p}}\frac{\partial\chi_{I}}{\partial x_{q}}\right|_{{\bf x}=0}$

Equation (55) is helpful to compute the covariance matrix after the model is
trained to determine various features as, e.g., the degree of entanglement.

##  IV Gates as linear layers

Gates represented by unitary operations are of relevance in many applications.
We are interested to those gates such that the new annihilation operators are
expressed as a linear combination of the input operators. If the new operator
is

$\hat{\widetilde{\bm{a}}}=\hat{U}^{\dagger}\hat{\bm{a}}\hat{U}\;,$ (56)

in term of the density matrix, the transformation reads [47]

$\tilde{\rho}=\hat{U}\rho\,\hat{U}^{\dagger}.$ (57)

We start considering a transformation such that [47]

$\displaystyle\hat{\widetilde{\bm{a}}}=\bm{\mathrm{U}}\hat{\bm{a}}.$ (58)

with ${\bm{\mathrm{U}}}$ a $n\times n$ complex matrix, with $\hat{\bm{a}}$ and
$\hat{\widetilde{\bm{a}}}$ column vectors of operators with dimension $n\times
1$. In a later section, we will consider the more general case with

$\hat{\widetilde{\bm{a}}}=\bm{\mathrm{U}}\hat{\bm{a}}+\bm{\mathrm{W}}\hat{\bm{a}}^{\dagger}.$
(59)

For the moment, we have $\bm{\mathrm{W}}=0$.

The transformation in terms of the canonical variables read [48]

$\hat{\widetilde{\bf R}}=\hat{U}^{\dagger}\hat{{\bf R}}\hat{U}={\bf
M}\,\hat{{\bf R}}+\bm{d}^{\prime}$ (60)

Linear transformations transform Gaussian states into Gaussian states [48]. A
Gaussian state with covariance matrix $\bm{g}$ and displacement vector
$\bm{d}$, turns into a new Gaussian state with covariance matrix

$\widetilde{\bm{g}}={\bf M}\,{\bm{g}}\,{\bf M}^{\top}$ (61)

and displacement vector

$\widetilde{\bm{d}}={\bf M}\,\bm{d}+\bm{d}^{\prime}\text{.}$ (62)

If the transformation due to $\bm{\mathrm{U}}$ is unitary, the matrix ${\bf
M}$ is symplectic, and satisfies

${\bf M}\bm{\mathrm{J}}{\bf M}^{\top}=\bm{\mathrm{J}}$ (63)

with $\bm{\mathrm{J}}$ given in Eq. (25). The inverse of ${\bf M}$ is found as

${\bf M}^{-1}=\bm{\mathrm{J}}{\bf M}^{\top}\bm{\mathrm{J}}^{\top}\;.$ (64)

It is instructive to deepen the link between $\bm{\mathrm{U}}$ and ${\bf M}$
for later use. We consider the $n\times 1$ vector of the annihilation
operators $\hat{\bm{a}}$,

$\hat{\bm{a}}=\begin{pmatrix}\hat{a}_{0}\\\ \hat{a}_{1}\\\ \vdots\\\
\hat{a}_{n-1}\end{pmatrix}$ (65)

and the corresponding $n\times 1$ vectors of positions $\hat{\bm{q}}$ and
momentum vector $\hat{\bm{p}}$,

$\hat{\bm{q}}=\begin{pmatrix}\hat{q}_{0}\\\ \hat{q}_{1}\\\ \vdots\\\
\hat{q}_{n-1}\end{pmatrix}\hskip
28.45274pt\hat{\bm{p}}=\begin{pmatrix}\hat{p}_{0}\\\ \hat{p}_{1}\\\ \vdots\\\
\hat{p}_{n-1}\end{pmatrix}.$ (66)

We build the $\hat{\bf R}$ vector in Eq. (22) by using auxiliary rectangular
matrices ${\bf R}_{\mathrm{q}}$ and ${\bf R}_{\mathrm{p}}$ as follows

$\hat{\bm{\mathrm{R}}}={\bf R}_{\mathrm{q}}\hat{\bm{q}}+{\bf
R}_{\mathrm{p}}\hat{\bm{p}}.$ (67)

${\bf R}_{\mathrm{q}}$ and ${\bf R}_{\mathrm{p}}$ are matrices with size
$N\times n$.

For $N=4$, we have

${\bf R}_{\mathrm{q}}=\begin{pmatrix}1&0\\\ 0&0\\\ 0&1\\\ 0&0\end{pmatrix}$
(68)

and

${\bf R}_{\mathrm{p}}=\begin{pmatrix}0&0\\\ 1&0\\\ 0&0\\\ 0&1\end{pmatrix},$
(69)

such that

${\bf R}_{\mathrm{q}}\hat{{\bf q}}=\begin{pmatrix}\hat{q}_{0}\\\ 0\\\
\hat{q}_{1}\\\ 0\end{pmatrix},$ (70)

and

${\bf R}_{\mathrm{p}}\hat{{\bf p}}=\begin{pmatrix}0\\\ \hat{p}_{0}\\\ 0\\\
\hat{p}_{1}\end{pmatrix},$ (71)

and hence

${\bf\hat{R}}={\bf R}_{\mathrm{q}}\hat{{\bf q}}+{\bf
R}_{\mathrm{p}}\hat{\bm{p}}=\begin{pmatrix}\hat{q}_{0}\\\ \hat{p}_{0}\\\
\hat{q}_{1}\\\ \hat{p}_{1},\end{pmatrix}\;,$ (72)

being

$\hat{\bm{q}}=\begin{pmatrix}\hat{q}_{0}\\\ \hat{q}_{1}\end{pmatrix},$ (73)

and

$\hat{\bm{p}}=\begin{pmatrix}\hat{p}_{0}\\\ \hat{p}_{1}\end{pmatrix}.$ (74)

From the previous expressions, one has

$\begin{array}[]{l}\hat{\bm{q}}={\bf
R}_{\mathrm{q}}^{\top}\hat{\bm{\mathrm{R}}}\\\ \hat{\bm{p}}={\bf
R}_{\mathrm{p}}^{\top}\hat{\bm{\mathrm{R}}}\end{array}$ (75)

and the rectangular matrices ${\bf R}_{\mathrm{q}}$ and ${\bf R}_{\mathrm{p}}$
satisfy

$\begin{array}[]{l}{\bf R}_{\mathrm{q}}^{\top}{\bf R}_{\mathrm{q}}={\bf
1}_{n}\\\ {\bf R}_{\mathrm{q}}^{\top}{\bf R}_{\mathrm{p}}={\bf 0}_{n}\\\ {\bf
R}_{\mathrm{p}}^{\top}{\bf R}_{\mathrm{p}}={\bf 1}_{n}\\\ {\bf
R}_{\mathrm{p}}^{\top}{\bf R}_{\mathrm{q}}={\bf 0}_{n}\\\ {\bf
R}_{\mathrm{q}}{\bf R}_{\mathrm{q}}^{\top}+{\bf R}_{\mathrm{p}}{\bf
R}_{\mathrm{p}}^{\top}={\bf 1}_{N}\\\ \bm{\mathrm{J}}{\bf
R}_{\mathrm{q}}=-{\bf R}_{\mathrm{p}}\\\ \bm{\mathrm{J}}{\bf
R}_{\mathrm{p}}={\bf R}_{\mathrm{q}}\\\ {\bf R}_{\mathrm{q}}{\bf
R}_{\mathrm{p}}^{\top}-{\bf R}_{\mathrm{p}}{\bf
R}_{\mathrm{q}}^{\top}=\bm{\mathrm{J}}\end{array}\;.$ (76)

By the matrices $\bm{\mathrm{R}}_{q}$ and $\bm{\mathrm{R}}_{p}$, we can write
the $\hat{\bm{a}}$ vector as follows

$\hat{\bm{a}}=\frac{\hat{\bm{q}}+\imath\hat{\bm{p}}}{\sqrt{2}}=\frac{{\bf
R}_{\mathrm{q}}^{\top}+\imath{\bf R}_{\mathrm{p}}^{\top}}{\sqrt{2}}\hat{\bf
R}.$ (77)

Letting $\bm{\mathrm{U}}=\bm{\mathrm{U}}_{R}+\imath\bm{\mathrm{U}}_{I}$ with
real part $\bm{\mathrm{U}}_{R}$ and imaginary part $\bm{\mathrm{U}}_{I}$; and
$\bm{\mathrm{U}}^{*}=\bm{\mathrm{U}}_{R}-\imath\bm{\mathrm{U}}_{I}$.
Correspondingly

$\begin{array}[]{l}\hat{\widetilde{\bm{q}}}=\bm{\mathrm{U}}_{R}\hat{\bm{q}}-\bm{\mathrm{U}}_{I}\hat{\bm{p}}\\\
\hat{\widetilde{\bm{p}}}=\bm{\mathrm{U}}_{I}\hat{\bm{q}}+\bm{\mathrm{U}}_{R}\hat{\bm{p}}\end{array}$
(78)

As $\hat{\bm{\mathrm{R}}}$ transforms in $\hat{\widetilde{\bm{\mathrm{R}}}}$,
we have

$\hat{\widetilde{\bm{\mathrm{R}}}}={\bf
R}_{\mathrm{q}}\hat{\widetilde{\bm{q}}}+{\bf
R}_{\mathrm{p}}\hat{\widetilde{\bm{p}}}={\bf
M}\hat{\bm{\mathrm{R}}}=\left({\bf M}_{1}+{\bf
M}_{2}\right)\hat{\bm{\mathrm{R}}}$ (79)

with ${\bf M}_{1}$, ${\bf M}_{2}$, and ${\bf M}$ real matrices such that

$\begin{array}[]{l}{\bf M}={\bf M}_{1}+{\bf M}_{2}\\\ {\bf M}_{1}={\bf
R}_{\mathrm{q}}\bm{\mathrm{U}}_{R}{\bf R}_{\mathrm{q}}^{\top}+{\bf
R}_{\mathrm{p}}\bm{\mathrm{U}}_{R}{\bf R}_{\mathrm{p}}^{\top}\\\ {\bf
M}_{2}={\bf R}_{\mathrm{p}}\bm{\mathrm{U}}_{I}{\bf R}_{\mathrm{q}}^{\top}-{\bf
R}_{\mathrm{q}}\bm{\mathrm{U}}_{I}{\bf R}_{\mathrm{p}}^{\top}\end{array}$ (80)

Equations (80) are used in defining linear gates starting from their matrix
$\bm{\mathrm{U}}$.

##  V Gaussian states as neural networks

NNs are mathematical tools to approximate unknown functions of many variables.
NNs depend on many parameters that can be fine tuned to improve fitting a
target model. Gaussian states correspond to NNs with a Gaussian activation
function, which is often adopted, for example, in kernel machines [40].

Figure 1: Double-head model for a layer representing a Gaussian.

We implement a Gaussian density matrix by a “multi-head model” that is a
model, which has a computational backbone that processes the inputs and
multiple output layers to return different quantities. We start considering as
“heads” the real and imaginary parts $\chi_{R}$ and $\chi_{I}$.

Figure 1 shows the NN model for a Gaussian state with each layer having one
input ${\bf x}$ and two outputs (the “heads”) $\chi_{R}$ and $\chi_{I}$. The
input is the state row vector ${\bf x}$, and is the input layer. The outputs
$\chi_{R}$ and $\chi_{I}$ are the output layers. The inner layer denoted
“Gaussian state,” computes the Gaussian characteristic function. The
parameters of this layer, namely the covariance matrix $\bm{g}$ and the
displacement $\bm{d}$ are also indicated.

For a Gaussian state with covariance matrix ${\bf g}$ and displacement vector
${\bf d}$ we have

$\begin{array}[]{ll}\chi_{R}({\bf x})&=\exp(-\frac{1}{4}{\bf x}{\bf g}{\bf
x}^{\top})\cos({\bf x}\bm{d})\\\\[2.0pt] \chi_{I}({\bf
x})&=\exp(-\frac{1}{4}{\bf x}{\bf g}{\bf x}^{\top})\sin({\bf
x}\bm{d}).\end{array}$ (81)

However, it is convenient to consider a more general model that includes a
further input vector (the “bias”).

We generalize the Gaussian layer to include a further bias input ${\bf a}$
with the same dimensions of $\bm{d}$, in addition to ${{\bf x}}$, such that
the output is

$\chi(x)=e^{-\frac{1}{4}{\bf x}{\bf g}{\bf x}^{T}}e^{\imath(\bm{d}+\bm{a})},$
(82)

Here $\bm{a}$ is a bias in the displacement, which will be useful for
cascading layers as detailed below. Equation (82) corresponds to the two heads
of the model returning

$\begin{array}[]{ll}\chi_{R}({\bf x})&=\exp(-\frac{1}{4}{\bf x}{\bf g}{\bf
x}^{T})\cos[{\bf x}(\bm{d}+{\bf a})]\\\\[2.0pt] \chi_{I}({\bf
x})&=\exp(-\frac{1}{4}{\bf x}{\bf g}{\bf x}^{T})\sin[{\bf x}(\bm{d}+{\bf
a})].\end{array}$ (83)

Figure 2 shows a graphical representation of the generalized Gaussian NN with
bias.

Figure 2: Double-head model representing a Gaussian state with a bias.

##  VI Pullback

We setup our NN representation of the density matrix as a layered sequence of
gates. Each gate is a unitary operator that acts on the density matrix, and
transform the latter into a new state. To detail how it works, we start
considering the transformation of the characteristic function under unitary
operations.

We consider a linear transformation by the operator $\hat{U}$. The operators
and the density matrix changes as follows

$\begin{array}[]{l}\hat{\widetilde{\bm{a}}}=\hat{U}^{\dagger}\hat{\bm{a}}\hat{U}\\\
\tilde{\rho}=\hat{U}\rho\,\hat{U}^{\dagger}.\end{array}$ (84)

Correspondingly, for the transformed characteristic functions $\tilde{\chi}$,
we have

$\displaystyle\tilde{\chi}$
$\displaystyle=\text{Tr}\left[\widetilde{\rho}e^{\imath{\bf
x}\hat{\bm{\mathrm{R}}}}\right]=\text{Tr}\left[\hat{U}\rho\hat{U}^{\dagger}e^{\imath{\bf
x}\hat{\bm{\mathrm{R}}}}\right]=\text{Tr}\left[\rho\hat{U}^{\dagger}e^{\imath{\bf
x}\hat{\bm{\mathrm{R}}}}\hat{U}\right]$ (85)
$\displaystyle=\text{Tr}\left[\rho e^{\imath{\bf
x}\hat{U}^{\dagger}\hat{\bm{\mathrm{R}}}\hat{U}}\right]=\text{Tr}\left[\rho
e^{\imath{\bf x}\hat{\widetilde{\bm{\mathrm{R}}}}}\right].$

Recalling Eq. (60)

$\hat{\widetilde{\bm{\mathrm{R}}}}={\bf M}\bm{\mathrm{R}}+\bm{d}^{\prime},$
(86)

we remark that the linear transformation is determined by the symplectic
matrix ${\bf M}$ and the displacement $\bm{d}^{\prime}$. We write
$\tilde{\chi}$ in terms of ${\bf M}$ and $\bm{d}^{\prime}$, by using Eq. (86),

$\tilde{\chi}=\text{Tr}\left[\rho e^{\imath{\bf x}{\bf
M}\hat{\bm{\mathrm{R}}}}\right]e^{\imath{\bf x}\bm{d}^{\prime}}.$ (87)

On the other hand, the expression for $\chi({\bf x})$ is the following

$\chi=\text{Tr}\left[\rho e^{\imath{\bf x}\hat{\bm{\mathrm{R}}}}\right].$ (88)

Thus we have

$\tilde{\chi}({\bf x})=\chi({\bf x}{\bf M})e^{\imath{\bf
x}\bm{d}^{\prime}}=\chi({\bf x}{\bf M})e^{\imath({\bf x}{\bf M})({\bf
M}^{-1}\bm{d}^{\prime})}.$ (89)

From Eq. (89), we see that the modified characteristic function depends on the
modified input vector ${\bf x}{\bf M}$ and has a modified displacement ${\bf
M}^{-1}\bm{d}^{\prime}$.

Introducing the new input vector ${\bf y}={\bf x}{\bf M}$ and the bias
$\bm{a}={\bf M}^{-1}\bm{d}^{\prime}$, we can express the transformed
characteristic function $\tilde{\chi}$ in terms of the original $\chi$ as
follows

$\tilde{\chi}({\bf x})=\chi({\bf y})e^{\imath{\bf y}\bm{a}}.$ (90)

We give a graphical representation of the transformations in Eq. (90). We
start considering a characteristic function with input vector ${\bf x}$ and
multi-headed outputs $\chi_{R}$ and $\chi_{I}$, as in figure 3. Then we
generalize the model by introducing the bias vector $\bm{a}$ as in the
previous section and shown in Fig. 4. Fig. 5 represents the transformation in
equation Eq. (90).

Figure 3: Graphical rapresentation of the model for the characteristic
function. Figure 4: Graphical rapresentation of the generalized model for the
characteristic function with bias $\bm{a}$, which is adopted in the linear
transformations. Figure 5: Graphical representation of the transformed
characteristic function in terms of the original characteristic function with
modified input and bias.

We introduce a new layer, which we call “linear layer”, having as parameters
the symplectic matrix ${\bf M}$ and the displacement $\bm{d}^{\prime}$. The
linear layer is shown in Fig. 6, and it has two inputs ${\bf x}$ and $\bm{a}$,
and two outputs: ${\bf y}={\bf x}{\bf M}$, a row vector with the same size of
${\bf x}$, and a new displacement $\bm{b}={\bf
M}^{-1}(\bm{d}^{\prime}+\bm{a})$ a column vector with the size of $\bm{d}$.
The linear layer will enable us to represent many different linear
transformations, as squeezing or Glauber operators.

By using the linear layer, we represent the transformation by cascading two
layers. The resulting model is Fig. 7. The model can be described as the
pullback of the linear layer from the characteristic function layer. The
matrix ${\bf M}$ first acts on the input ${\bf x}$, and then the function
$\chi$ is evaluated on the vector ${\bf y}$. A simplified diagram is given in
Fig. 8. This representation will be useful when considering multiple
transformations.

Figure 6: Linear layer to transform the input variables to the $\chi$ layer.
Figure 7: Graphical representation of a linear transformation as a pullback of
a linear layer from the original characteristic function. The resulting model
corresponds the transformed characteristic function in Fig. 5. Note that here
$\bm{a}=0$. Figure 8: Simplified model of Fig. 7, omitting the internal
variables ${\bf y}$ and $\bm{b}$ and the zero input bias $\bm{a}$.

For a Gaussian state [see Eq. (34)] the transformed characteristic function
reads

$\tilde{\chi}({\bf x})=e^{-\frac{1}{4}{\bf x}{\bf M}\bm{g}{\bf M}^{\top}{\bf
x}^{\top}}e^{\imath{\bf x}{\bf M}(\bm{d}+{\bf M}^{-1}\bm{d}^{\prime})}.$ (91)

Equation Eq. (91) is graphically represented in Fig. 9, by using the Gaussian
layer in Fig. 2. This representation is equivalent to Fig. 10, by the linear
layer, whose action is detailed in Fig. 7.

Figure 9: Single-layer multiheaded model for a linear transformation of a
Gaussian state. Figure 10: Two-layer multiheaded model equivalent to Fig. 9.

As shown in Fig. 10 linear transformation on the density matrix is actually a
transformation of the variable ${\bf x}$. This implies that the linear
transformation can be represent first by a linear gate, followed by the
Gaussian gate. This may be described as “pulling back” the linear operation
before the Gaussian gate. By using these two models, the linear transformation
of the Gaussian state is a cascade of a linear pullback layer and a Gaussian
state layer, as shown in figure 10.

##  VII Pullback cascading

The pullback is helpful when being in the presence of multiple
transformations. A sequence of linear transforms is equivalent to a sequence
of pullbacks in reverse order as sketched in figure 11.

Figure 11: Linear transformations and pullbacks. A single transformation with
unitary operator $\hat{U}_{1}$ corresponds to a single pullback. A double
transformation with unitary operators $\hat{U}_{1}$ and $\hat{U}_{2}$
corresponds to a double pullback. Note that the flow of data from ${\bf x}$ to
the output through the network is in reverse order with respect to the two
transformations.

For example, we consider a system originally described by the density matrix
$\rho$ and canonical observables $\hat{\bm{\mathrm{R}}}$. First, the system is
subject to a linear transformation with operator $\hat{U}_{1}$, such that the
density matrix becomes

$\hat{U}_{1}\rho\hat{U}^{\dagger}_{1}$ (92)

and the new canonical vector is

$\hat{R}_{1}={\bf M}_{1}\hat{R}+\bm{d}_{1}.$ (93)

We then consider a second transformation with unitary operator $\hat{U}_{2}$
and parameters ${\bf M}_{2}$ and $\bm{d}_{2}$, such that the final density
matrix reads

$\hat{U}_{2}\hat{U}_{1}\rho\hat{U}_{1}^{\dagger}\hat{U}_{2}^{\dagger},$ (94)

and the final vector of observables is

$\hat{\bm{\mathrm{R}}}_{2}={\bf
M}_{2}\hat{\bm{\mathrm{R}}}_{1}+\bm{d}_{2}={\bf M}_{2}{\bf
M}_{1}\hat{\bm{\mathrm{R}}}+\bm{d}_{2}+{\bf M}_{2}\bm{d}_{1}\text{.}$ (95)

We remark that ${\bf M}_{1}{\bf M}_{2}\neq{\bf M}_{2}{\bf M}_{1}$, hence the
order of the two tranformations is relevant as the corresponding unitary
operators $\hat{U}_{1}$ and $\hat{U}_{2}$ do not commute. The sequence of the
two linear transformations is equivalent to a single transformation with
parameters

$\begin{array}[]{l}{\bf M}={\bf M}_{2}{\bf M}_{1}\\\ \bm{d}^{\prime}={\bf
M}_{2}\bm{d}_{1}+\bm{d}_{2}\end{array}$ (96)

Figure 12: Equivalence of two cascaded linear layers with a single linear
layer.

We have the following (figure 12) Proposition — The sequence of the two linear
transfomations $\hat{U}_{1}$ and $\hat{U_{2}}$ is equivalent to the pullback
of the two linear layers with parameters $({\bf M}_{1},\bm{d}_{1})$ and $({\bf
M}_{2},\bm{d}_{2})$.

###### Proof.

As indicated in Fig. 12 the ouput of the linear layer corresponding to
$\hat{U}_{2}$ is

$\begin{array}[]{l}{\bf x}_{2}={\bf x}{\bf M}_{2}\\\ \bm{a}_{2}={\bf
M}_{2}^{-1}\left(\bm{d}_{2}+\bm{a}\right)\end{array}$ (97)

The output of the linear layer corresponding to $\hat{U}_{1}$ is

$\begin{array}[]{l}{\bf x}_{1}={\bf x}_{2}{\bf M}_{1}={\bf x}{\bf M}_{2}{\bf
M}_{1}\\\ \bm{a}_{1}={\bf M}_{1}^{-1}\left(\bm{d}_{1}+\bm{a}_{2}\right)={\bf
M}_{1}^{-1}\left(\bm{d}_{1}+{\bf M}_{2}^{-1}\bm{d}_{2}+{\bf
M}_{2}^{-1}\bm{a}\right)\end{array}$ (98)

For the layer with parameters $({\bf M},\bm{d}^{\prime})$ we have, see Eq.
(96),

$\displaystyle{\bf y}$ $\displaystyle={\bf x}{\bf M}={\bf x}{\bf M}_{2}{\bf
M}_{1}$ (99) $\displaystyle\bm{b}$ $\displaystyle={\bf
M}^{-1}\left(\bm{d}^{\prime}+\bm{a}\right)$ $\displaystyle={\bf
M}_{1}^{-1}{\bf M}_{2}^{-1}\left({\bf
M}_{2}\bm{d}_{1}+\bm{d}_{2}+\bm{a}\right)$ $\displaystyle={\bf
M}_{1}^{-1}\left(\bm{d}_{1}+{\bf M}_{2}^{-1}\bm{d}_{2}+{\bf
M}_{2}^{-1}\bm{a}\right)$

As ${\bf y}={\bf x}_{1}$ and $\bm{b}=\bm{a}_{1}$ we have the proof. ∎

By extending the previous argument to three or more transformations, one
realizes that the cascade of an arbitrary number of transformations
corresponds to a cascade of pullbacks in reverse order. For $M$
transformations, reverse order means that one first make a pullback of
operator $1$, then operator $2$, etc. In the flow of data in the network, the
input ${\bf x}$ first enters operator $M$, then $M-1$ and so forth to passing
through the linear layer corresponding to the transformation $1$. The case
$M=3$ is shown as an example in figure 13.

Figure 13: Example of cascading three transformations. Note the reverse order
of layers from right to left.

##  VIII The Glauber displacement layer

Starting from the linear layer, we can define specialized layers corresponding
to different unitary operators. The first we describe is the displacement
operator, or Glauber operator, defined by [47]

$\hat{\mathcal{D}}(\alpha)=\exp\left(\alpha^{*}\hat{a}-\alpha\hat{a}^{\dagger}\right)\;.$
(100)

A coherent state is obtained by displacing the vacuum state:

$|\alpha\rangle=\hat{\mathcal{D}}(\alpha)|0\rangle.$ (101)

For a single mode, letting $\hat{U}=\hat{\mathcal{D}}(\alpha)$, we have

$\hat{\widetilde{a}}=\hat{U}^{\dagger}\hat{a}\hat{U}=\hat{a}+\alpha,$ (102)

which implies for the canonical vector

$\hat{\widetilde{\bm{\mathrm{R}}}}=\hat{U}^{\dagger}\hat{\bm{\mathrm{R}}}\hat{U}=\begin{pmatrix}\hat{x}\\\
\hat{p}\end{pmatrix}+\begin{pmatrix}d_{0}\\\ d_{1}\end{pmatrix}$ (103)

with

$\begin{array}[]{l}d_{0}=\sqrt{2}\Re(\alpha)\\\
d_{1}=\sqrt{2}\Im(\alpha).\end{array}$ (104)

For a many-body displacement $\hat{\mathcal{D}}(\bm{\alpha})$, wit
${\bm{\alpha}}=(\alpha_{0},\alpha_{1},...,\alpha_{n})^{\top}$, we have

$\hat{\widetilde{\bm{a}}}=\hat{U}^{\dagger}\hat{\bm{a}}\hat{U}=\hat{\bm{a}}+\bm{\alpha},$
(105)

which implies for the canonical vector

$\hat{\widetilde{\bm{\mathrm{R}}}}=\hat{U}^{\dagger}\hat{\bm{\mathrm{R}}}\hat{U}=\hat{\bm{\mathrm{R}}}+\bm{d}$
(106)

with ($j=0,1,\ldots,n-1$)

$\displaystyle d_{2j}$ $\displaystyle=\sqrt{2}\Re(\alpha_{j})$ (107)
$\displaystyle d_{2j+1}$ $\displaystyle=\sqrt{2}\Im(\alpha_{j})$

For the Glauber layer we hence have ${\bf M}={\bf 1}_{N}$.

To create a NN model that represents a coherent state with a given
displacement vector $\bm{d}_{\text{target}}$, one can start from the vacuum
Gaussian state with $\bm{g}={\bf 1}_{N}$ and $\bm{d}={\bf 0}$ and pullback a
linear gate with ${\bf M}={\bf 1}_{N}$, and displacement
$\bm{d}_{\text{target}}$. No bias is needed (see Fig. 14 with $\bm{a}={\bf
0}$).

Figure 14: Model for a coherent state obtained by pullback of a displacement
operator (Glauber layer) from the vacuum.

##  IX The squeezing layer

Different nonclassical states are generated from vacuum by unitary operators
resulting into the following linear relation

$\hat{\widetilde{\bm{a}}}=\bm{\mathrm{U}}\hat{\bm{a}}+\bm{\mathrm{W}}\hat{\bm{a}}^{\dagger},$
(108)

which generalizes Eq. (58). Using the matrices $\bm{\mathrm{U}}$ and
$\bm{\mathrm{W}}$, one obtains the corresponding symplectic matrix ${\bf M}$.

Let

$\bm{\mathrm{U}}=\bm{\mathrm{U}}_{R}+\imath\bm{\mathrm{U}}_{I}$ (109)

and

$\bm{\mathrm{W}}=\bm{\mathrm{W}}_{R}+\imath\bm{\mathrm{W}}_{I}.$ (110)

We write the $\hat{\widetilde{\bm{a}}}$ as follows

$\displaystyle\hat{\widetilde{\bm{a}}}$ $\displaystyle=$
$\displaystyle\displaystyle\frac{\hat{\widetilde{\bm{q}}}+\imath\hat{\widetilde{\bm{p}}}}{\sqrt{2}}=\displaystyle\bm{\mathrm{U}}\hat{\bm{a}}+\bm{\mathrm{W}}\hat{\bm{a}}^{\dagger}$
(111) $\displaystyle\hat{\widetilde{\bm{a}}}^{\dagger}$ $\displaystyle=$
$\displaystyle\displaystyle\frac{\hat{\widetilde{\bm{q}}}-\imath\hat{\widetilde{\bm{p}}}}{\sqrt{2}}=\bm{\mathrm{U}}^{*}\hat{\bm{a}}^{\dagger}+\bm{\mathrm{W}}^{*}\hat{\bm{a}}.$
(112)

We have

$\begin{array}[]{l}\displaystyle\hat{\widetilde{\bm{q}}}=(\bm{\mathrm{U}}_{R}+\bm{\mathrm{W}}_{R})\hat{\bm{q}}+(-\bm{\mathrm{U}}_{I}+\bm{\mathrm{W}}_{I})\hat{\bm{p}}\vspace{0.2cm}\\\
\displaystyle\hat{\widetilde{\bm{p}}}=(\bm{\mathrm{U}}_{I}+\bm{\mathrm{W}}_{I})\hat{\bm{q}}+(\bm{\mathrm{U}}_{R}-\bm{\mathrm{W}}_{R})\hat{\bm{p}}\vspace{0.2cm}\end{array}$
(113)

which generalize Eq. (78) Following the same arguments for Eq. (80), we have
for the symplectic matrix

$\begin{array}[]{l}{\bf M}={\bf M}_{1}+{\bf M}_{2}\vspace{0.2cm}\\\ {\bf
M}_{1}={\bf R}_{\mathrm{q}}(\bm{\mathrm{U}}_{R}+\bm{\mathrm{W}}_{R}){\bf
R}_{\mathrm{q}}^{\top}+{\bf
R}_{\mathrm{p}}(\bm{\mathrm{U}}_{R}-\bm{\mathrm{W}}_{R}){\bf
R}_{\mathrm{p}}^{\top}\vspace{0.2cm}\\\ {\bf M}_{2}={\bf
R}_{\mathrm{p}}(\bm{\mathrm{U}}_{I}+\bm{\mathrm{W}}_{I}){\bf
R}_{\mathrm{q}}^{\top}+{\bf
R}_{\mathrm{q}}(-\bm{\mathrm{U}}_{I}+\bm{\mathrm{W}}_{I}){\bf
R}_{\mathrm{p}}^{\top}\vspace{0.2cm}\end{array}$ (114)

Equation (114) is used for programming layers representing squeezing
operators.

###  IX.1 Single-mode squeezed state

We first consider a single model squeezed state, so that $N=2$, and

$\hat{\bm{\mathrm{R}}}=\begin{pmatrix}\hat{x}\\\ \hat{p}\end{pmatrix}$ (115)

Using the squeezing operator with parameters $r$ and $\theta$ [47]

$\hat{\widetilde{a}}=\hat{S}^{\dagger}\hat{a}\hat{S}=\cosh(r)\hat{a}-e^{\imath\theta}\sinh(r)\hat{a}^{\dagger}$
(116)

we have that the matrix $\bm{\mathrm{U}}_{R,I}$ and $\bm{\mathrm{W}}_{R,I}$
are complex scalars as follows

$\begin{array}[]{l}U_{R}=\cosh(r)\\\ U_{I}=0\\\ W_{R}=-\cos(\theta)\sinh(r)\\\
W_{I}=-\sin(\theta)\sinh(r)\end{array}$ (117)

and, by Eq. (114), we find the symplectic operator for the squeezing

$\begin{array}[]{ll}&{\bf M}_{s}(r,\theta)=\\\\[2.0pt]
&\left(\begin{array}[]{cc}\cosh(r)-\cos(\theta)\sinh(r)&-\sin(\theta)\sinh(r)\\\
-\sin(\theta)\sinh(r)&\cosh(r)+\cos(\theta)\sinh(r)\end{array}\right)\end{array}$
(118)

For $\theta=0$, i.e., for a real squeezing parameter, we have

${\bf M}_{s}(r,0)=\begin{pmatrix}\exp(-r)&0\\\ 0&\exp(r)\end{pmatrix}.$ (119)

###  IX.2 Multi-mode squeezed vacuum model

We are interested to multi-mode systems, so we consider a $n-$body state, and
apply the squeezing operator to one mode. The single-mode squeezing operator
is obtained by a linear gate with $\bm{d}^{\prime}=0$ and ${\bf M}$ given by
Eq. (118) for the mode to be squeezed denoted by an index
$n_{\text{squeezed}}$ in the range $0$ to $n-1$ (see Fig. 15). For example,
for $n=4$, and $n_{\text{squeezed}}=0$, we have

${\bf M}=\begin{pmatrix}M_{s,11}&M_{s,12}&0&0&0&0&0&0\\\
M_{s,21}&M_{s,22}&0&0&0&0&0&0\\\ 0&0&1&0&0&0&0&0\\\ 0&0&0&1&0&0&0&0\\\
0&0&0&0&1&0&0&0\\\ 0&0&0&0&0&1&0&0\\\ 0&0&0&0&0&0&1&0\\\
0&0&0&0&0&0&0&1\end{pmatrix}.$ (120)

Seemingly for $n=4$, $n_{\text{squeezed}}=2$

${\bf M}=\begin{pmatrix}1&0&0&0&0&0&0&0\\\ 0&1&0&0&0&0&0&0\\\
0&0&1&0&0&0&0&0\\\ 0&0&0&1&0&0&0&0\\\ 0&0&0&0&M_{s,11}&M_{s,12}&0&0\\\
0&0&0&0&M_{s,21}&M_{s,22}&0&0\\\ 0&0&0&0&0&0&1&0\\\
0&0&0&0&0&0&0&1\end{pmatrix}.$ (121)

Figure 15: A squeezed coherent state obtained by pullback with a single mode
squeezing layer from the vacuum.

###  IX.3 Squeezed coherent states

The squeezed operator is represented by a linear gate, we can cascade with
other layers, such as the Glauber displacement layer. Figure 16a shows a model
to generate a squeezed coherent state from the squeezed vacuum by a
displacement layer.

The squeezed coherent states $|\alpha,\zeta\rangle$ are built by applying the
displacement operator $\hat{\mathcal{D}}(\alpha)$ to a squeezed vacuum
$\hat{\mathcal{S}}(\zeta)|0\rangle$, where $\zeta=r\,e^{\imath\theta}$ is the
complex squeezing parameter, that is

$|\alpha,\zeta\rangle=\hat{\mathcal{D}}(\alpha)\hat{\mathcal{S}}(\zeta)|0\rangle$
(122)

The resulting state has $\langle\hat{a}\rangle=\alpha$, and it is squeezed.

Figure 16: Pullback model for a squeezed coherent state by displacing a
squeezed vacuum (top panel). Pullback model for squeezed coherent state by
squeezing a displaced vacuum (bottom panel). Here all the layers act on the
same mode.

###  IX.4 Squeezing the displaced vacuum

A different squeezed coherent state is obtained by squeezing a coherent state.
We first pullback a displacement layer from the vacuum to create the coherent
state, and then pullback a squeezing layer, as shown in Fig. 16b. This
corresponds to the following equation Eq. (123)

$|\alpha\cosh(r)-\alpha^{*}e^{\imath\theta}\sinh(r),\zeta\rangle=\hat{\mathcal{S}}(\zeta)\hat{\mathcal{D}}(\alpha)|0\rangle$
(123)

The result is a squeezed coherent state with the same eigenvalues for the
covariance matrix as above, but the mean value of the displacement is changed
[47], i.e.,

$\langle\hat{a}\rangle=\alpha\cosh(r)-\alpha^{*}e^{\imath\theta}\sinh(r).$
(124)

##  X Computing the Hamiltonian

In this manuscript, we are interested in the state of a many-body
Hamiltonians, which correspond to classical solitons. We consider the
Hamiltonian

$\hat{H}=-\sum_{j}\left[\hat{\psi}^{\dagger}_{j}\hat{\psi}_{j+1}+\hat{\psi}^{\dagger}_{j}\hat{\psi}_{j-1}\right]+\frac{\gamma}{2}\hat{\psi}^{\dagger}_{j}\hat{\psi}^{\dagger}_{j}\hat{\psi}_{j}\hat{\psi}_{j}=\hat{K}+\hat{V}\;,$
(125)

with the potential energy

$\hat{V}=\frac{\gamma}{2}\sum_{j}\left(\hat{n}_{j}^{2}-\hat{n}_{j}\right)\;$
(126)

being $\hat{n}_{i}=\hat{\psi}^{\dagger}_{i}\hat{\psi}_{i}$, and the kinetic
energy

$\hat{K}=\sum_{ij}\omega_{ij}\hat{\psi}^{\dagger}_{i}\hat{\psi}_{j}\;,$ (127)

being $\omega_{ij}=-\delta_{i,j-1}-\delta_{i,j+1}$ with $i,j$ in $[0,n-1]$ and
with homogeneous boundary conditions, $\omega_{0,-1}=\omega_{n-1,n}=0$.

Starting from the vacuum state, one can build the NN model of an arbitrary
state by multiple pullbacks. We compute the mean value of observable
quantities as $\hat{H}$ and $\hat{\mathcal{N}}$ as derivatives of $\chi$ at
${\bf x}=0$. We have in symmetric ordering [47]

$\langle\hat{K}\rangle=\left.\sum_{jk}\omega_{jk}\left[\frac{\partial^{2}}{\partial\alpha_{j}\partial(-\alpha_{k}^{*})}\chi(\bm{\alpha})-\frac{\delta_{jk}}{2}\chi\right]\right|_{\bm{\alpha}=0}$
(128)

where
$\bm{\alpha}=\left(\alpha_{0}\ldots\alpha_{n-1},\alpha^{*}_{0}\ldots\alpha^{*}_{n-1}\right)$,
and $\sqrt{2}\alpha_{j}=x_{2j}+\imath x_{2j+1}$, and for the interaction term

$\langle\hat{V}\rangle=\frac{\gamma}{2}\sum_{j}\left.\frac{\partial^{4}\chi}{\partial\alpha_{j}^{2}{\partial(-\alpha_{j}^{*})}^{2}}-2\frac{\partial^{2}\chi}{\partial\alpha_{j}{\partial(-\alpha_{j}^{*})}}+\frac{\chi}{2}\right|_{\bm{\alpha}=0}\text{.}$
(129)

These quantities can be computed by using automatic differentiation on the NN
model.

Specifically,

$\displaystyle\langle\hat{K}\rangle$ $\displaystyle=-\frac{1}{2}\sum_{mn}$
(130) $\displaystyle\left(\frac{\partial^{2}\chi}{\partial q_{m}\partial
q_{n}}+\frac{\partial^{2}\chi}{\partial p_{m}\partial
p_{n}}+\chi\right)\omega_{mn}^{R}$ $\displaystyle+$
$\displaystyle\left.\left(\frac{\partial^{2}\chi}{\partial q_{m}\partial
p_{n}}-\frac{\partial^{2}\chi}{\partial p_{n}\partial
p_{q}}\right)\omega_{mn}^{I}\right|_{{\bf x}=0}$

where

$\omega_{mn}=\omega_{mn}^{R}+\imath\omega_{mn}^{I}\;,$ (131)

$q_{j}=x_{2j}$ and $p_{j}=x_{2j+1}$ with $j=0,1,...,{N/2}-1$.

For the potential energy,

$\langle\hat{V}\rangle=\sum_{nm}V_{nm}\langle\hat{a}_{n}^{\dagger}\hat{a}_{m}^{\dagger}\hat{a}_{n}\hat{a}_{m}\rangle\;,$
(132)

being, for a local interaction,

$V_{nm}=\frac{1}{2}\delta_{nm}\;.$ (133)

For a mode with index $j=0,1,\ldots,n-1$, the normal ordering product is

$\langle\hat{n}_{j}^{2}\rangle=\langle\hat{a}^{\dagger}_{j}\hat{a}_{j}\hat{a}^{\dagger}_{j}\hat{a}_{j}\rangle=\langle\hat{a}^{\dagger}_{j}\hat{a}^{\dagger}_{j}\hat{a}_{j}\hat{a}_{j}\rangle+\langle\hat{a}^{\dagger}_{j}\hat{a}_{j}\rangle\;.$
(134)

where we used
$[\hat{a}_{j},\hat{a}_{j}^{\dagger}]=\hat{a}_{j}\hat{a}^{\dagger}_{j}-\hat{a}^{\dagger}_{j}\hat{a}_{j}=1$.

We have

$\displaystyle\langle\hat{a}^{\dagger}_{j}\hat{a}_{j}\rangle=\langle\hat{n}_{j}\rangle=\left.-\frac{\partial^{2}\chi}{\partial
z_{j}~{}\partial z_{j}^{*}}\right|_{\bm{z}=0}=$ (135)
$\displaystyle-\left.\frac{1}{2}\left(\partial^{2}_{q_{j}}+\partial^{2}_{p_{j}}\right)\chi_{R}\right|_{{\bf
x}=0}-\frac{1}{2}\;.$

and

$\displaystyle\langle\hat{a}^{\dagger}_{j}\hat{a}^{\dagger}_{j}\hat{a}_{j}\hat{a}_{j}\rangle=\left.{\left(\frac{\partial}{\partial
z_{j}}\right)}^{2}{\left(-\frac{\partial}{\partial
z_{j}^{*}}\right)}^{2}\chi\right|_{{\bf x}=0}$ (136)
$\displaystyle-2\left.\left(\frac{\partial}{\partial
z_{j}}\right)\left(-\frac{\partial}{\partial
z_{j}^{*}}\right)\chi\right|_{{\bf x}=0}+\frac{1}{2}\;.$

Eq. (136) can be expressed in terms of the derivatives w.r.t. $q_{j}$ and
$p_{j}$ as follows

$\displaystyle\langle\hat{a}^{\dagger}_{j}\hat{a}^{\dagger}_{j}\hat{a}_{j}\hat{a}_{j}\rangle=\left.\frac{1}{4}{\left(\partial_{q_{j}}^{2}+\partial_{p_{j}}^{2}\right)}^{2}\chi_{R}\right|_{{\bf
x}=0}$ (137)
$\displaystyle+\left.{\left(\partial_{q_{j}}^{2}+\partial_{p_{j}}^{2}\right)}\chi_{R}\right|_{{\bf
x}=0}+\frac{1}{2}$

where, as above, $q_{j}=x_{2j}$ and $p_{j}=x_{2j+1}$ with $j=0,1,...,{N/2}-1$.

The real part $\chi_{R}$ enters in Eq. (137) as
$\langle\hat{n}_{j}^{2}\rangle$ in Eq. (134) is real valued.

In general, we need to evaluate the fourth order derivatives
$\frac{\partial^{4}\chi}{\partial x_{s}\partial x_{p}\partial x_{q}\partial
x_{r}}$, at ${\bf x}=0$. The computation of the fourth-order derivatives may
be demanding as it grows with $N^{4}$. For the specific case of Gaussian
states, it can be simplified as the fourth order derivatives can be expressed
in terms of the covariance matrix $\bm{g}$.

We have by direct derivation of Eq. (1)

$\displaystyle\left.\frac{\partial^{4}\chi}{\partial x_{s}\partial
x_{p}\partial x_{q}\partial x_{r}}\right|_{{\bf x}=0}$ (138)
$\displaystyle=\frac{1}{4}g_{sp}g_{qr}+\frac{1}{4}g_{sq}g_{pr}+\frac{1}{4}g_{sr}g_{pq}$
$\displaystyle+\frac{1}{2}g_{sp}d_{q}\,d_{r}+\frac{1}{2}g_{sq}d_{p}\,d_{r}+\frac{1}{2}g_{sr}d_{p}\,d_{q}$
$\displaystyle+\frac{1}{2}g_{pq}d_{r}\,d_{s}+\frac{1}{2}g_{pr}d_{q}\,d_{s}+\frac{1}{2}g_{qr}d_{p}\,d_{s}$
$\displaystyle+d_{s}d_{p}d_{q}d_{r}\;.$

When considering the diagonal terms we have (no implicit sum is present in the
following equation)

$\left.\frac{\partial^{4}\chi}{\partial x_{q}^{4}}\right|_{{\bf
x}=0}=\frac{3}{4}g_{qq}^{2}+3\,g_{qq}d_{q}^{2}+d_{q}^{4}\;.$ (139)

Seemingly, we have

$\displaystyle\left.\frac{\partial^{4}\chi}{\partial x_{q}^{2}\partial
x_{p}^{2}}\right|_{{\bf x}=0}=\frac{1}{4}g_{qq}g_{pp}+\frac{1}{2}g_{pq}^{2}+$
(140)
$\displaystyle\frac{1}{2}g_{pp}d_{q}^{2}+\frac{1}{2}g_{qq}d_{p}^{2}+2g_{pq}d_{p}\,d_{q}+d_{p}^{2}d_{q}^{2}\;.$

Hence, for Gaussian states, we do not need to evaluate explicitly the fourth
order derivatives, but we can combine the componentes of $\bm{g}$ and
$\bm{d}$. For the second derivatives we have

$\begin{array}[]{l}\displaystyle\left.\frac{\partial^{2}{\chi_{R}}}{{\partial
x_{m}}{\partial x_{n}}}\right|_{{\bf x}=0}=-\frac{1}{2}g_{mn}-d_{m}d_{n}\\\
\displaystyle\left.\frac{\partial^{2}{\chi_{I}}}{{\partial x_{m}}{\partial
x_{n}}}\right|_{{\bf x}=0}\textbf{}=0\\\ \end{array}\;.$ (141)

##  XI Quantum soliton variational ansatz

Figure 17 shows the quantum processor to synthesize quantum lattice solitons
in a discrete array with $n$ sites. The circuit is composed by $n$ squeezers
with complex squeezing parameters $\zeta_{j}$, and displacement operators with
complex displacements $\delta_{j}$, being $j=0,1,\ldots,n-1$ hereafter. A
unitary interferometer $\hat{U}$ mixes the modes before entering the waveguide
array. The quantum processor is represented as a NN in the phase space. All
the parameters of the squeezers, displacements, and interferometers are
trainable and initialized as random variables. We refer to the state generated
by the quantum circuit as the quantum soliton variational ansatz (QSVA).

Figure 17: (a) Quantum processor for solitons in coupled optical waveguides.
$n$ vacuum spatial modes are squeezed, displaced and mixed by a trainable
interferometer $\hat{U}$ ($n=5$ in the figure). (b) Multiple pullbacks as a
quantum soliton variational ansatz.

We build the NN model of an arbitrary state by multiple pullbacks from the
vacuum. Figure 17b shows the NN model for the QSVA in Fig. 17a.

We obtain the mean values of $\hat{H}$ and $\hat{\mathcal{N}}$ as derivatives
of $\chi$. As described above, for a Gaussian state, derivatives are obtained
algebraically by $\bm{g}$ and $\bm{d}$, which speeds up training by
conventional steepest descent. We use as cost function
$\exp(\langle\hat{H}\rangle/n)$, ($\langle\hat{H}\rangle<0$ at a minimum). As
additional cost function, we use ${(\langle\mathcal{\hat{N}}\rangle-
N_{T})}^{2}$, to constraint the target mean boson number as
$\langle\mathcal{\hat{N}}\rangle=N_{T}$.

The first challenge is training the variational quantum circuit to generate
single quantum solitons. The simplest soliton is localized in one of the sites
(as $\gamma<0$). We want to understand if the trained QSVA can approximate a
quantum version of the single soliton, such that one can use quantum
processors to generate quantum solitons and investigate their physics.

We find that the NN furnishes quantum solitons with lower energy than the
classical solution so that the quantum solitons form a larger class of
nonlinear waves. Starting from randomly generated weights, the soliton is
found after training the NN model which minimizes the cost function, and
$\langle\hat{H}\rangle$ at fixed $\langle\hat{\mathcal{N}}\rangle=N_{T}$.

We first consider small interaction strength $|\gamma|$, at which we expected
delocalized solutions. We show in Fig. 18a the displacements
$|\langle\psi_{j}\rangle|^{2}=(d_{2j}^{2}+d_{2j+1}^{2})/2$ after the training.
In Fig. 18b we show the mean boson number $\langle\hat{n}_{j}\rangle$, for
$N_{T}=10$ and $\gamma=-0.01$ after thousands of epochs; the profile is
delocalized. Strong localization is obtained for $N_{T}=10$ and $\gamma=-1$ as
in Fig. 18b,c.

Figure 18:  (a,b) Ground state for $N_{T}=10$ and $\gamma=-0.01$: (a)
displacement $|\langle\hat{\psi}_{j}\rangle|^{2}$; (b) mean boson number
$\langle\hat{n}_{j}\rangle$; (c,d) as in (a,b) with $N_{T}=10$ and $\gamma=-1$
($n=17$). Figure 19:  Training history for Fig. 18: (a) $\langle H\rangle$
rapidly reaches the minimum for $\gamma=-0.01$ (dashed line is the coherent
state with $\langle\hat{H}\rangle\simeq-2N_{T}$); (b) as in (a) for
$\gamma=-1$. Dashed line is Eq. (143). The model explores plateaus
corresponding to solitons with different degree of entanglement; (c)
logarithmic negativity for $\gamma=-0.01$, no entanglement (dashed line); (d)
as in (c) for $\gamma=-1$, the ground state is entangled.

In figure 19, we report $\langle\hat{H}\rangle$ when varying the number of
training epochs for $N_{T}=10$. The algorithm converges after some thousands
of epochs, and explores different solutions. Figure 19a shows the training
when $\gamma=-0.01$: the NN converges to the delocalized state in Fig. 18a,b.
To compare with analytical estimates, we consider the continuum limit and a
non-squeezed coherent state with a sinusoidal profile. The value of the
Hamiltonian $\langle\hat{H}\rangle\simeq-2N_{T}$ (dashed line in fig. 19a).
Figure 19b shows the strong interaction case $\gamma=-1$. During the training
the systems settles in two plateaus corresponding to
$\langle\hat{H}\rangle\cong-50$ and $\langle\hat{H}\rangle\cong-155$. These
values are understood in terms of approximate solutions, in which
$\langle\hat{\psi}_{j\neq A}\rangle=0$, with
$A=\left\lfloor\frac{n}{2}\right\rfloor+1$, where the soliton is localized. If
we consider the many-body state
$|\bm{\alpha}\rangle=\hat{D}_{A}(\bm{\alpha})|\text{\bf vac}\rangle$ with
$\alpha_{i}=\delta_{iA}\alpha$, we have

$\langle\hat{H}\rangle\simeq\frac{\gamma}{2}\left(\sum_{i}\langle|\hat{n_{i}}^{2}|\rangle-\langle|\hat{n_{i}}|\rangle\right)=\frac{\gamma}{2}\sum_{i}|\alpha_{i}|^{4}=\frac{\gamma}{2}|\alpha|^{4}$
(142)

One has $\langle\hat{N}\rangle=|\alpha|^{2}=N_{T}$ and
$\langle\hat{H}\rangle=\frac{\gamma}{2}N_{T}^{2}$ which gives
$\langle\hat{H}\rangle\simeq-50$ (first plateau in Fig. 19, dashed line). The
difference due to the neglected $\langle\hat{K}\rangle$.

The lowest energy solution occurs at a larger number of epochs (Fig. 19b). The
plateau can be estimated by considering a squeezed coherent state
$|\bm{\alpha},\zeta\rangle=\hat{D}_{A}(\bm{\alpha})\hat{S}_{A}(\zeta)|\textbf{vac}\rangle$
where $\hat{S}_{A}$ is the squeezing operator at $j=A$ with
$\zeta=re^{\imath\theta}$. For $|\bm{\alpha},\zeta\rangle$ one has

$\displaystyle\langle\mathcal{\hat{N}}\rangle$
$\displaystyle={\sinh(r)}^{2}+|\alpha|^{2}=N_{T}$ (143)
$\displaystyle\langle\hat{H}\rangle$
$\displaystyle=\frac{5}{8}-2|\alpha^{2}|+|\alpha|^{4}+(-1+2|\alpha|^{2})\cosh(2r)+$
$\displaystyle\frac{3}{8}\cosh(4r)-|\alpha|^{2}\cos(2\phi-\theta)\sinh(2r),$

with minimum $\langle\hat{H}\rangle\simeq-155$ ($N_{T}=10$, red dashed line).

##  XII Entangled solitons

Having evidence that the ground state is not the classical solution, we are
interested to understand if other non-classical features are present. We
consider precisely the degree of entanglement.

Entanglement in Gaussian states is estimated by the logarithmic negativity
[51, 49, 50, 13]. We partition the system in the soliton site, “Alice” A at
$A=\left\lfloor\frac{n}{2}\right\rfloor+1$, and the sites with $j\neq A$,
“Bob”.

We retrieve the covariance matrix $\bm{g}$ from the trained NN model, and its
partial transpose $\tilde{\bm{g}}$ [51], obtained by multiplying by $-1$ the
elements of Alice momenta $\hat{p}_{A}$. The symplectic eigenvalues
$\tilde{c}_{0},\tilde{c}_{1},\ldots,\tilde{c}_{n-1}$ are moduli of the
eigenvalues of ${\bf J}^{\top}\tilde{\bm{g}}/2$, and logarithmic negativity is

$E_{\mathcal{N}}=-\sum_{j=0}^{n-1}\log_{2}\min\\{1,2c_{j}\\}.$ (144)

Figure 19c shows the $E_{\mathcal{N}}$ for $\gamma=-0.01$ during the training.
Despite initially the system explores randomly generated entangled solution,
the asymptotic state is a coherent state, with vanishing entanglement (dashed
line) as expected for a negligible interaction.

Figure 19d shows $E_{\mathcal{N}}$ for $\gamma=-1$. The local minimum of
$\langle\hat{H}\rangle$, the coherent state, has a smaller entanglement than
the squeezed ground state at epoch $3\times 10^{4}$. Nevertheless, the
asymptotic value of $E_{\mathcal{N}}$ is negligible if compared with the
bounded solitons considered in the following.

Given the non-classical features in the single-soliton ground state, we study
the broader family of multiple solitons. We expect that strong nonlinear
localization in different channels increases the degree of entanglement. We
study two-soliton bound states [30, 31, 1], which can be approximated by
training the network by proper loss functions. We fix the location of two
sites, denoted $A$ and $B=A+\Delta$, and we introduce as cost function
$\exp(-\langle\hat{n}_{A}\rangle)$ to maximize the bosons at $A$ and
$\exp(\langle\hat{n}_{A}-\hat{n}_{B}\rangle)$ to generate a solution with the
same boson number in $B$.

Figure 20 shows the coupled solitons for different peak positions
($N_{T}=40$). At variance with the single soliton in Fig. 18, the boson number
is not vanishing in all the sites, but ${|\langle\psi_{j}\rangle|}^{2}$ is not
negligible only at the solitons, outlining the onset of squeezed nonlocal
states.

The bound solitons are entangled. If we consider a bipartite system composed
by the site in $A$ and the rest of the array, $E_{\mathcal{N}}$ depends on the
distance $\Delta$ between the two solitons (Fig. 20). The logarithmic
negativity depends on the number of total bosons $N_{T}$: Figure 20b shows
$E_{\mathcal{N}}$ versus $N_{T}$ for various $\Delta$. The entanglement is
small and comparable with the single soliton for case $N_{T}<10$, and growing
when $N_{T}>10$.

Figure 20:  (a) Boson distribution and mean displacement (right axis) for two
representative bound-soliton solutions with $N_{T}=40$ and $\gamma=-1.0$; (b)
logarithmic negativity versus the boson number $N_{T}$ for different $\Delta$;
(c) probabilities for events containing $\bar{n}_{T}=2$ particles in two
varying sites $A$ and $B$ and $N_{T}=10$; (d) as in (c) with 2 particles per
site ($\bar{n}_{T}=4$). The probability is not vanishing only when the
particles are at the soliton sites, denoting the generation of entangled pairs
($n=10$).

##  XIII Boson sampling

As we show in the following the entanglement is connected to correlated
particles in different output channels. We compute the probability
$\Pr(\bar{\bf n})$ of finding $\bar{n}_{0}$ bosons in mode $0$, $\bar{n}_{1}$
in mode $1$, and so forth, where $\bar{\bf
n}=(\bar{n}_{0},\bar{n}_{1},\ldots,\bar{n}_{n-1})$ is a given particle
pattern. Letting $\rho$ the density matrix, one has $\Pr(\bar{\bf
n})=\text{Tr}[\rho|\bar{\bf n}\rangle\langle\bar{\bf n}|]$, with $|\bar{\bf
n}\rangle\langle\bar{\bf
n}|=\otimes_{j}|\bar{n}_{j}\rangle\langle\bar{n}_{j}|$. Correspondingly [52],

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

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-representation of the density matrix [53, 47].

We obtain $\Pr(\bar{\bf n})$ as in [54]. We consider an event with
$\bar{n}_{T}$ particles in two sites $A$ and $B$. Specifically, the events
such that (i) $\bar{n}_{A}=\bar{n}_{B}=n_{T}/2$ for $A\neq B$, or (ii)
$\bar{n}_{A}=\bar{n}_{B}=\bar{n}_{T}$ when $A=B$, and $\bar{n}_{j}=0$
elsewhere. For $n_{T}=2$ this corresponds to observing a pair of particles in
a single site, or two particles in two distinct sites.

Figure 20c shows the probability for $\bar{n}_{T}=2$ varying $A$ and $B$ when
$N_{T}=10$. Such a probability is different from zero only when $A$ or $B$ are
at a soliton position ($A=2$ and $B=7$ in Fig. 20c). When $A=B$ (yellow bar)
we have the probability of observing the two particles in one single soliton.
When $A$ and $B$ corresponds to the sites of the two solitons (green bar), we
have the probability of observing a couple of particles in the two soliton
sites. For all other patterns with $\bar{n}_{T}=2$ there is a negligible
probability of observing a pair. Hence, pairs of particles appear either in a
single soliton, or entangled in the two solitons. When $\bar{n}_{T}=4$ we have
non-vanishing probability of observing either $4$ bosons in a soliton, or two
entangled pairs in the two solitons (Fig. 20d). Particles are hence observed
simultaneously only when they are monitored at the two localized solitons,
demonstrating quantum correlation and the nonlocality.

##  XIV Conclusion

Discrete waveguide arrays became popular for the quantum advantage in linear
boson sampling. One can argue about the role of nonlinearity, which –
classically– supports discrete self-localized solitons.

For unveiling quantum effects in discrete solitons, we adopted quantum machine
learning. Variational quantum circuits enable the theoretical investigation of
quantum nonlinear waves and their bound states. Also, their physical
realization as tunable quantum processors opens to new experiments on non-
perturbative nonlinear regimes.

We considered Gaussian states as they are expected to conform to experimental
observations at high intensities, and also allowed us to rigorously obtain the
degree of entanglement. Furthermore, the neural network formulation is
compatible with the Gaussian boson sampling protocol in the presence of
nonlinearity.

We found that solitonic states alter the observation of multi-particle events.
Entangled pairs emerge from self-localized bounded solitons, which can be
synthesized by trained quantum circuits.

Quantum machine learning and quantum variational algorithms open new
possibilities for the physics and the applications of quantum nonlinear waves.
The methodology can be extended beyond Gaussian states, and generalized to
continuous systems and multidimensional arrays with arbitrary networks.

###### Acknowledgements.

We acknowledge support from Horizon 2020 QuantERA grant QUOMPLEX, by National
Research Council (CNR), Grant 731473, and PRIN PELM (20177PSCKT).

## Appendix:Graph and parameters of the model

Figure 21 shows a graphical representation of the graph of the model.

Figure 21: Graph of the multi-head model for the quantum solitons. One can
notice the central computational core represented by the cascade of multi-head
layers, including the interferometers LinearRandomMultiHead, displacements
LinearConstantMultiHead and squeezing layers SingleModeSqueezer. Then we have
the activation layer, computing the Gaussian characteristic function
ResidualGaussianMultiHead, and the layer to determine the covariance matrix
CovarianceLayer. After that a series of “Lambda” layers, i.e., simple
computing layers for the cost functions and the various metrics during the
training.

We optimize the model in order to minimize the Hamiltonian, the resulting
parameters give the circuit that produce a state with the minimum value of
energy. The trainable parameters in the model, i.e., those corresponding to
the $n$ displacement layers (each layer with $1$ complex parameter), the $n$
squeezing layers (each layer with $1$ complex parameter), and the unitary
matrix representing the trainable interferometer (with $N^{2}/4$ independent
parameters). Thus the model has $N^{2}/4+N$ independent variables. Training is
done by the Adam algorithm (TensorFlow v2.7.0). Gaussian boson sampling from
the trained model is done following [54].

## References

  * Kivshar and Agrawal [2003] Y. Kivshar and G. P. Agrawal,  _Optical solitons_ (Academic Press, New York, 2003).
  * Conti [2022] C. Conti, Fluctuations and non-gaussianity in 3d+1 quantum nonlocal nonlinear waves: towards quantum gravity and quantum fluid technologies, arXiv:2202.10741 (2022).
  * Friberg  _et al._ [1992] S. R. Friberg, S. Machida, and Y. Yamamoto, Quantum-nondemolition measurement of the photon number of an optical soliton, Phys. Rev. Lett. 69, 3165–3168 (1992).
  * Drummond  _et al._ [1993] P. D. Drummond, R. M. Shelby, S. R. Friberg, and Y. Yamamoto, Quantum solitons in optical fibres, Nature 365, 307–313 (1993).
  * Silberhorn  _et al._ [2001] C. Silberhorn, P. Lam, O. Weis, F. Konig, N. Korolkova, _et al._ , Generation of Continuous Variable Einstein-Podolsky-Rosen Entanglement via the Kerr Nonlinearity in an Optical Fiber, Phys.Rev.Lett. 86, 4267–4270 (2001).
  * Yu  _et al._ [2001] C. X. Yu, H. A. Haus, and E. P. Ippen, Soliton squeezing at the gigahertz rate in a Sagnac loop, Opt. Lett. 26, 669–671 (2001).
  * Liang  _et al._ [2018] Q.-Y. Liang, A. V. Venkatramani, S. H. Cantu, T. L. Nicholson, M. J. Gullans, A. V. Gorshkov, J. D. Thompson, C. Chin, M. D. Lukin, and V. Vuletić, Observation of three-photon bound states in a quantum nonlinear medium, Science 359, 783 (2018).
  * Leuchs  _et al._ [1999] G. Leuchs, T. C. Ralph, C. Silberhorn, and N. Korolkova, Scheme for the generation of entangled solitons for quantum communication, Journal of Modern Optics 46, 1927 (1999).
  * Lee  _et al._ [2005] R.-K. Lee, Y. Lai, and B. A. Malomed, Generation of photon-number-entangled soliton pairs through interactions, Phys. Rev. A 71, 013816 (2005).
  * Lai and Lee [2009] Y. Lai and R.-K. Lee, Entangled quantum nonlinear schrödinger solitons, Phys. Rev. Lett. 103, 013902 (2009).
  * Villari  _et al._ [2018] L. D. M. Villari, D. Faccio, F. Biancalana, and C. Conti, Quantum soliton evaporation, Phys. Rev. A 98, 043859 (2018).
  * Marchukov  _et al._ [2020] O. V. Marchukov, B. A. Malomed, V. Dunjko, J. Ruhl, M. Olshanii, R. G. Hulet, and V. A. Yurovsky, Quantum fluctuations of the center of mass and relative parameters of nonlinear schrödinger breathers, Phys Rev Lett 125, 050405 (2020).
  * Martynov  _et al._ [2020] V. O. Martynov, V. O. Munyaev, and L. A. Smirnov, Generation of entangled states of light using discrete solitons in waveguide arrays, arXiv:2011.07662 (2020).
  * Potasek and Yurke [1987] M. J. Potasek and B. Yurke, Squeezed-light generation in a medium governed by the nonlinear schrödinger equation, Phys. Rev. A 35, 3974 (1987).
  * Lai and Haus [1989] Y. Lai and H. A. Haus, Quantum theory of solitons in optical fibers. i. time-dependent hartree approximation, Phys. Rev. A 40, 844 (1989).
  * Wright [1991] E. M. Wright, Quantum theory of soliton propagation in an optical fiber using the Hartree approximation, Phys. Rev. A 43, 3836–3844 (1991).
  * Scott  _et al._ [1994] A. Scott, J. Eilbeck, and H. Gilhøj, Quantum lattice solitons, Physica D: Nonlinear Phenomena 78, 194 (1994).
  * Yao [1995] D. Yao, Quantum fluctuations of optical solitons in fibers, Phys. Rev. A 52, 4871–4881 (1995).
  * Marcucci  _et al._ [2020a] G. Marcucci, D. Pierangeli, and C. Conti, Theory of neuromorphic computing by waves: machine learning by rogue waves, dispersive shocks, and solitons, Phys. Rev. Lett. 125, 093901 (2020a).
  * Herrera [2020] R. A. Herrera, Evaluating a neural network and a convolutional neural network for predicting soliton properties in a quantum noise environment, J. Opt. Soc. Am. B 37, 3094 (2020).
  * Silva  _et al._ [2021] N. A. Silva, T. D. Ferreira, and A. Guerreiro, Reservoir computing with solitons, New J. Phys. 23, 023013 (2021).
  * Pu  _et al._ [2021] J.-C. Pu, J. Li, and Y. Chen, Soliton, breather, and rogue wave solutions for solving the nonlinear schrödinger equation using a deep learning method with physical constraints, Chinese Physics B 30, 060202 (2021).
  * Lopez-Pastor and Marquardt [2021] V. Lopez-Pastor and F. Marquardt, Self-learning machines based on hamiltonian echo backpropagation, arXiv:2103.04992 (2021).
  * Tillmann  _et al._ [2012] M. Tillmann, B. Dakić, R. Heilmann, S. Nolte, A. Szameit, and P. Walther, Experimental boson sampling, Nat. Photonics 7, 540 (2012).
  * Broome  _et al._ [2013] M. A. Broome, A. Fedrizzi, S. Rahimi-Keshari, J. Dove, S. Aaronson, T. C. Ralph, and A. G. White, Photonic boson sampling in a tunable circuit, Science 339, 794 (2013).
  * Spring  _et al._ [2013] 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, and I. A. Walmsley, Boson sampling on a photonic chip, Science 339, 798 (2013).
  * Carolan  _et al._ [2014] 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, and A. Laing, On the experimental verification of quantum complexity in linear optics, Nat. Photonics 8, 621 (2014).
  * Spagnolo  _et al._ [2014] 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, and F. Sciarrino, Experimental validation of photonic boson sampling, Nat. Photonics 8, 615 (2014).
  * Zhong  _et al._ [2020] H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C. Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, 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, and J.-W. Pan, Quantum computational advantage using photons, Science 370, 1460 (2020).
  * Lederer  _et al._ [2008] F. Lederer, G. I. Stegeman, D. N. Christodoulides, G. Assanto, M. Segev, and Y. Silberberg, Discrete solitons in optics, Physics Reports 463, 1 (2008).
  * Kartashov  _et al._ [2011] Y. V. Kartashov, B. A. Malomed, and L. Torner, Solitons in nonlinear lattices, Rev. Mod. Phys. 83, 247–305 (2011).
  * Aaronson and Arkhipov [2013] S. Aaronson and A. Arkhipov, The computational complexity of linear optics, Theory Comput. 9, 143 (2013).
  * Carleo  _et al._ [2019] G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, Machine learning and the physical sciences, Rev. Mod. Phys. 91, 045002 (2019).
  * Broughton  _et al._ [2020] 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, and M. Mohseni, Tensorflow quantum: A software framework for quantum machine learning, arXiv:2003.02989 (2020).
  * Marcucci  _et al._ [2020b] G. Marcucci, D. Pierangeli, P. W. H. Pinkse, M. Malik, and C. Conti, Programming multi-level quantum gates in disordered computing reservoirs via machine learning, Opt. Express 28, 14018 (2020b).
  * Marquardt [2021] F. Marquardt, Machine learning and quantum devices, arXiv:2101.01759 (2021).
  * Tacchino  _et al._ [2021] F. Tacchino, S. Mangini, P. K. Barkoutsos, C. Macchiavello, D. Gerace, I. Tavernelli, and D. Bajoni, Variational learning for quantum artificial neural networks, arXiv:2103.02498 (2021).
  * Kardashin  _et al._ [2021] A. Kardashin, A. Pervishko, J. Biamonte, and D. Yudin, Numerical hardware-efficient variational quantum simulation of a soliton solution, Phys. Rev. A 104, L020402 (2021).
  * Huang  _et al._ [2021] H.-Y. Huang, R. Kueng, G. Torlai, V. V. Albert, and J. Preskill, Provably efficient machine learning for quantum many-body problems, arXiv:2106.12627 (2021).
  * Schuld and Petruccione [2018] M. Schuld and F. Petruccione,  _Supervised Learning with Quantum Computers_ , edited by Springer (Springer, 2018).
  * Rebentrost  _et al._ [2014] P. Rebentrost, M. Mohseni, and S. Lloyd, Quantum support vector machine for big data classification, Phys. Rev. Lett. 113, 130503 (2014).
  * Wall and D’Aguanno [2021] M. L. Wall and G. D’Aguanno, Tree-tensor-network classifiers for machine learning: From quantum inspired to quantum assisted, Phys. Rev. A 104, 042408 (2021).
  * Vicentini  _et al._ [2019] F. Vicentini, A. Biella, N. Regnault, and C. Ciuti, Variational neural network ansatz for steady states in open quantum systems, Phys. Rev. Lett. 122, 250503 (2019).
  * Cerezo  _et al._ [2021] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, Variational quantum algorithms, Nature Reviews Physics 3, 625 (2021).
  * Arrazola  _et al._ [2021] J. M. Arrazola, V. Bergholm, K. Brádler, T. R. Bromley, M. J. Collins, I. Dhand, A. Fumagalli, T. Gerrits, A. Goussev, L. G. Helt, J. Hundal, T. Isacsson, R. B. Israel, J. Izaac, S. Jahangiri, R. Janik, N. Killoran, S. P. Kumar, J. Lavoie, A. E. Lita, D. H. Mahler, M. Menotti, B. Morrison, S. W. Nam, L. Neuhaus, H. Y. Qi, N. Quesada, A. Repingon, K. K. Sabapathy, M. Schuld, D. Su, J. Swinarton, A. Száva, K. Tan, P. Tan, V. D. Vaidya, Z. Vernon, Z. Zabaneh, and Y. Zhang, Quantum circuits with many photons on a programmable nanophotonic chip, Nature 591, 54 (2021).
  * Hoch  _et al._ [2021] F. Hoch, S. Piacentini, T. Giordani, Z.-N. Tian, M. Iuliano, C. Esposito, A. Camillini, G. Carvacho, F. Ceccarelli, N. Spagnolo, A. Crespi, F. Sciarrino, and R. Osellame, Boson sampling in a reconfigurable continuously-coupled 3d photonic circuit, arXiv:2106.08260 (2021).
  * Barnett and Radmore [1997] S. M. Barnett and P. M. Radmore,  _Methods in Theoretical Quantum Optics_ (Oxford University Press, New York, 1997).
  * Wang  _et al._ [2007] X. Wang, T. Hiroshima, A. Tomita, and M. Hayashi, Quantum information with gaussian states, Phys. Rep. 448, 1 (2007).
  * Adesso  _et al._ [2004] G. Adesso, A. Serafini, and F. Illuminati, Extremal entanglement and mixedness in continuous variable systems, Phys. Rev. A 70, 022318 (2004).
  * Plenio and Virmani [2007] M. B. Plenio and S. Virmani, An introduction to entanglement measures, Quantum Inf. Comput. 7, 1 (2007).
  * Vidal and Werner [2002] G. Vidal and R. F. Werner, Computable measure of entanglement, Phys. Rev. A 65, 032314 (2002).
  * Kruse  _et al._ [2019] R. Kruse, C. S. Hamilton, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, A detailed study of gaussian boson sampling, Phys. Rev. A 100, 032326 (2019).
  * Gardiner and Zoller [2004] C. W. Gardiner and P. Zoller,  _Quantum Noise_ , 3rd ed. (Springer-Verlag, Berlin, 2004).
  * Conti [2021b] C. Conti, Training gaussian boson sampling by quantum machine learning, Quantum Machine Intelligence 3, 26 (2021b).