A pure state of finite stellar rank $\ket{\bm{\psi}}$ is Gaussian-separable if
and only if its core state $\ket{C_{\bm{\psi}}}$ (as in Lemma 10) is passively
separable.

###### Proof.

By Lemma 10, there exists a Gaussian unitary operation

$\hat{G}_{\bm{\psi}}=\hat{U}\hat{S}(\bm{\xi})\hat{D}(\bm{\beta}),$ (112)

with
$\bm{\beta}=(\beta_{1},\dots,\beta_{m}),\bm{\xi}=(\xi_{1},\dots,\xi_{m})\in\mathbb{C}^{m}$
and $\hat{U}$ a passive linear unitary operation such that
$\ket{\bm{\psi}}=\hat{G}_{\bm{\psi}}\ket{C_{\bm{\psi}}}$. Hence,

$\ket{\bm{\psi}}=\hat{U}\hat{S}(\bm{\xi})\hat{D}(\bm{\beta})\ket{C_{\bm{\psi}}}.$
(113)

If $\ket{C_{\bm{\psi}}}$ is passively separable, there exists a passive linear
unitary $\hat{V}$ such that $\hat{V}\ket{C_{\bm{\psi}}}$ is separable. Hence,
$\hat{V}\hat{D}^{\dagger}(\bm{\beta})\hat{S}^{\dagger}(\bm{\xi})\hat{U}^{\dagger}\ket{\bm{\psi}}$
is separable, and thus $\ket{\bm{\psi}}$ is Gaussian-separable.

Conversely, if $\ket{\bm{\psi}}$ is Gaussian-separable, there exists a
Gaussian unitary $\hat{G}$ such that $\hat{G}\ket{\bm{\psi}}$ is separable.
Hence, $\hat{G}\hat{G}_{\bm{\psi}}\ket{C_{\bm{\psi}}}$ is separable, so by
Lemmas 15 and 10 there exists Gaussian unitaries $\hat{G}_{I}$ and
$\hat{G}_{J}$ and core states $\ket{C_{I}}$ and $\ket{C_{J}}$ such that

$\hat{G}\hat{G}_{\bm{\psi}}\ket{C_{\bm{\psi}}}=(\hat{G}_{I}\ket{C_{I}})\otimes(\hat{G}_{J}\ket{C_{J}}),$
(114)

and thus

$\hat{G}_{\bm{\psi}}\ket{C_{\bm{\psi}}}=\hat{G}^{\dagger}(\hat{G}_{I}\otimes\hat{G}_{J})\ket{C_{I}}\otimes\ket{C_{J}}.$
(115)

Let $\hat{G}^{\prime}:=\hat{G}^{\dagger}(\hat{G}_{I}\otimes\hat{G}_{J})$. This
is a Gaussian unitary operation, so there exist
$\bm{\beta}^{\prime}=(\beta_{1}^{\prime},\dots,\beta_{m}^{\prime}),\bm{\xi}^{\prime}=(\xi_{1}^{\prime},\dots,\xi_{m}^{\prime})\in\mathbb{C}^{m}$
and $\hat{U}^{\prime},\hat{V}^{\prime}$ passive linear unitary operations such
that
$\hat{G}^{\prime}=\hat{U}^{\prime}\hat{S}(\bm{\xi}^{\prime})\hat{D}(\bm{\beta}^{\prime})\hat{V}^{\prime}$.
Let us define
$\ket{C_{\bm{\psi}}^{\prime}}:=\hat{V}^{\prime}(\ket{C_{I}}\otimes\ket{C_{J}})$.
Since $\hat{V}^{\prime}$ is a passive linear unitary,
$\ket{C_{\bm{\psi}}^{\prime}}$ is a core state, and is by construction
passively separable. Plugging these expressions into Eq. (115) we finally
obtain:

$\displaystyle\ket{\bm{\psi}}$
$\displaystyle=\hat{G}_{\bm{\psi}}\ket{C_{\bm{\psi}}}$ (116)
$\displaystyle=\hat{U}^{\prime}\hat{S}(\bm{\xi}^{\prime})\hat{D}(\bm{\beta}^{\prime})\ket{C_{\bm{\psi}}^{\prime}}.$

By the uniqueness of the decomposition in Lemma 10, this implies
$\hat{G}_{\bm{\psi}}=\hat{U}^{\prime}\hat{S}(\bm{\xi}^{\prime})\hat{D}(\bm{\beta}^{\prime})$
and $\ket{C_{\bm{\psi}}}=\ket{C_{\bm{\psi}}^{\prime}}$, so
$\ket{C_{\bm{\psi}}}$ is passively separable.

∎

In particular, in the case where $\ket{\bm{\psi}}$ is separable, then both
$\ket{C}$ and $\ket{C_{\bm{\psi}}}$ are separable. The converse does not hold:
in the case where $\ket{\bm{\psi}}$ is a Gaussian state, we have
$\ket{C}=\ket{C_{\bm{\psi}}}=\ket{\bm{0}}$, but $\ket{\bm{\psi}}$ may still be
entangled. A state which is not passively separable is referred to as
inherently entangled [104].

Summarising the results of this section, we have analysed the structure of
entanglement within the multimode stellar hierarchy: entangled states of
finite stellar rank may display both Gaussian and non-Gaussian entanglement,
and the latter can be thought of as the entanglement of a finite-dimensional
system within the infinite-dimensional quantum state.

##  5 Holomorphic Quantum Computing

Having generalised the properties of the stellar hierarchy to the multimode
case and studied the entanglement structure of states within the hierarchy, we
turn to the aspect that initially motivated this generalisation: quantum
computing. We first define a formal model of computing based on holomorphic
functions which we call Holomorphic Quantum Computing (HQC), and show that it
captures bosonic quantum computations.

Then, we study a subclass of circuits within this model which corresponds to
quantum computations preserving the stellar rank—the multimode analogue of the
single-mode evolutions from Sec. 3. We use the formalism developed in Sec. 4
to provide a simple classification of existing subuniversal infinite-
dimensional quantum computing models based on their computational complexity,
which highlights the interplay between non-Gaussianity, squeezing and
entanglement at the origin of possible quantum computational speedup in the
infinite-dimensional setting.

Finally, we consider adaptive versions of these rank-preserving circuits and
show that they are capable of universal BQP computations.

###  5.1 Definition of HQC

We define the model of Holomorphic Quantum Computing over $m$ variables as
follows (see Fig. 2):

  * •

The inputs are elements of Segal–Bargmann space, i.e., holomorphic functions
over $\mathbb{C}^{m}$ that are square-integrable with respect to the Gaussian
measure.

  * •

The gates are specified by differential unitary operators $\hat{U}=e^{-iHt}$
generated by Hamiltonians $H$ that are formal complex power series of the
operators $\times z_{k}$ (multiplication by the $k^{th}$ variable) and
$\partial_{z_{k}}$ (partial derivative with respect to the $k^{th}$ variable),
for $k=1,\dots,m$, satisfying
$H(z_{1},\partial_{z_{1}},\dots,z_{m},\partial_{z_{m}})=H^{*}(\partial_{z_{1}},z_{1},\dots,\partial_{z_{m}},z_{m})$.
These unitary gates map holomorphic functions in Segal–Bargmann space to
holomorphic functions in Segal–Bargmann space.

  * •

The outcome of the computation is sampled from the output holomorphic function
in two possible ways:

    * –

Continuous measurement: writing $F^{\star}$ the output state of the
computation, a continuous measurement of modes $(1,\dots,m)$ yields a tuple
$\bm{\alpha}=(\alpha_{1},\dots,\alpha_{m})\in\mathbb{C}^{m}$ sampled with the
probability density function
$\\{\frac{e^{-|\bm{\alpha}|^{2}}}{\pi^{m}}|F^{\star}(\bm{\alpha})|^{2}\\}_{\bm{\alpha}\in\mathbb{C}^{m}}$.

    * –

Discrete measurement: writing $F^{\star}$ the output state of the computation,
a discrete measurement of modes $(1,\dots,m)$ yields a tuple
$\bm{n}=(n_{1},\dots,n_{m})\in\mathbb{N}^{m}$ sampled from the probability
distribution
$\\{\frac{1}{\bm{n}!}|\partial_{\bm{z}}^{\bm{n}}F^{\star}(\bm{0})|^{2}\\}_{\bm{n}\in\mathbb{N}^{m}}$.

  * •

Additionally, we also consider adaptive versions of this model, i.e., where
unitary gates may depend on the outcomes of intermediate (continuous or
discrete) measurements.

This model captures bosonic quantum computations through the Segal–Bargmann
representation [43, 44] (see Sec. 2.3).

Indeed, the holomorphic functions may be thought of as the stellar functions
of corresponding input quantum states in an infinite-dimensional separable
Hilbert space. Through the correspondence in Eq. (34), the Hamiltonians are
computable functions of the creation and annihilation operators of the modes,
which naturally appear in bosonic experiments. This includes the standard
model of CV quantum computing defined by Lloyd and Braunstein [27] wherein
Hamiltonians are polynomials in the creation and annihilation operators. In
particular, Hamiltonians generating Gaussian unitaries are given by
polynomials of degree less or equal to $2$ in the creation and annihilation
operators.

Using the relation between the stellar function and the Husimi
quasiprobability density in Eq. (83), the continuous measurement corresponds
to a Gaussian measurement in coherent state basis (a standard heterodyne
detection in optics [105], up to a complex conjugation of the outcomes), i.e.,
sampling from the probability density function related to the values of the
stellar function over complex phase space. On the other hand, using
$\bra{n}=\bra{0}\frac{\hat{a}^{n}}{\sqrt{n!}}$ and Eq. (34), the discrete
measurement corresponds to a non-Gaussian measurement in Fock basis (a photon-
number measurement in optics), i.e., sampling from the probability
distribution related to the coefficients of the stellar function. These two
types of single-mode measurements are available in current experiments with
high efficiencies. If some of the modes (say, the first $k$ modes) are
measured with the continuous measurement and the remaining $m-k$ modes with
the discrete measurement, this yields instead a tuple
$(\bm{\alpha},\bm{n})=(\alpha_{1},\dots,\alpha_{k},n_{1},\dots,n_{m-k})\in\mathbb{C}^{k}\times\mathbb{N}^{m-k}$
sampled with the joint probability density:

$\left\\{\frac{e^{-|\bm{\alpha}|^{2}}}{\pi^{k}\bm{n}!}|\partial_{\bm{z}}^{\bm{n}}F^{\star}(\bm{\alpha},\bm{0})|^{2}\right\\}_{(\bm{\alpha},\bm{n})\in\mathbb{C}^{k}\times\mathbb{N}^{m-k}}.$
(117)

Note also that the commonly used homodyne detection (corresponding to a
measurement of a position-like operator) is obtained by combining a single-
mode Gaussian squeezing operation with a continuous heterodyne detection
(yielding an unbalanced heterodyne detection), in the limit where the modulus
of the squeezing parameter goes to infinity [105]. Moreover, any Gaussian
measurement may be obtained by combining Gaussian unitary operations and
homodyne detection [106].

Hence, the HQC model, while being defined using basic elements of complex
analysis, captures infinite-dimensional quantum computational models with
bosonic Hamiltonian evolution and Gaussian measurements and/or non-Gaussian
photon-number detection (non-adaptive and adaptive). Even though we only
consider these two types of measurements, note that Eq. (30) allows for
modeling arbitrary projections.

###  5.2 Rank-Preserving HQC

In this section, we consider the subclass of HQC circuits which preserve the
stellar rank, i.e., Rank-Preserving Holomorphic Quantum Circuits (RPHQC). By
Theorem 3, these correspond to Hamiltonian evolutions which are polynomials of
degree less of equal to $2$ of the operators $\times z_{k}$ and
$\partial_{z_{k}}$, for $k=1,\dots,m$, i.e., Gaussian Hamiltonians. Moreover,
by Lemma 11, the set of states of finite stellar rank is dense in
Segal–Bargmann space for the trace norm, so we restrict without loss of
generality to finite rank stellar functions, of the form $P\times G$, where
$P$ is a polynomial and $G$ is a Gaussian function. These circuits are thus
defined as follows (see Fig. 3):

  * •

The inputs are elements of Segal–Bargmann space of finite stellar rank, i.e.,
holomorphic functions over $\mathbb{C}^{m}$ that are square-integrable with
respect to the Gaussian measure, of the form $P\times G$, where $P$ is a
polynomial and $G$ is a Gaussian function.

  * •

The gates are specified by differential unitary operators $\hat{U}=e^{-iHt}$
generated by Hamiltonians $H$ that are polynomials of degree less of equal to
$2$ of the operators $\times z_{k}$ (multiplication by the $k^{th}$ variable)
and $\partial_{z_{k}}$ (partial derivative with respect to the $k^{th}$
variable), for $k=1,\dots,m$, satisfying
$H(z_{1},\partial_{z_{1}},\dots,z_{m},\partial_{z_{m}})=H^{*}(\partial_{z_{1}},z_{1},\dots,\partial_{z_{m}},z_{m})$.
These unitary gates map holomorphic functions of finite stellar rank in
Segal–Bargmann space to other holomorphic functions of the same finite stellar
rank in Segal–Bargmann space by Theorem 3.

  * •

The outcome of the computation is sampled from the output holomorphic function
in two possible ways:

    * –

Continuous measurement: writing $F^{\star}$ the output state of the
computation, a continuous measurement of modes $(1,\dots,m)$ yields a tuple
$\bm{\alpha}=(\alpha_{1},\dots,\alpha_{m})\in\mathbb{C}^{m}$ sampled with the
probability density function
$\\{\frac{e^{-|\bm{\alpha}|^{2}}}{\pi^{m}}|F^{\star}(\bm{\alpha})|^{2}\\}_{\bm{\alpha}\in\mathbb{C}^{m}}$.

    * –

Discrete measurement: writing $F^{\star}$ the output state of the computation,
a discrete measurement of modes $(1,\dots,m)$ yields a tuple
$\bm{n}=(n_{1},\dots,n_{m})\in\mathbb{N}^{m}$ sampled from the probability
distribution
$\\{\frac{1}{\bm{n}!}|\partial_{\bm{z}}^{\bm{n}}F^{\star}(\bm{0})|^{2}\\}_{\bm{n}\in\mathbb{N}^{m}}$.

  * •

Additionally, we also consider adaptive versions of this model, i.e., where
unitary gates may depend on the outcomes of intermediate (continuous or
discrete) measurements.

Since its unitary evolution is Gaussian, any RPHQC may be decomposed
efficiently as a combination of displacement, squeezing and passive linear
unitary gates (see Eq. (26)). The latter are described by a finite-dimensional
unitary matrix $U$ and act on holomorphic functions in Segal–Bargmann space as

$F^{\star}(\bm{z})\mapsto F^{\star}(U\bm{z}),$ (118)

where we used Eqs. (23) and (34), together with Lemma 7. In other terms,
passive linear operations induce a unitary change of complex coordinates in
Segal–Bargmann space. Moreover, the action of a displacement operator on mode
$k$ is obtained by fixing all but the $k^{th}$ variable $z_{k}$ and applying
Lemma 4. The case of squeezing, shearing, or any Gaussian quadratic non-
linearity follows from the single-mode case: for instance, assuming without
loss of generality that a squeezing gate is applied on mode $1$ of a stellar
function $F^{\star}$, we consider the univariate holomorphic function of
finite stellar rank given by the section

$z_{1}\mapsto
F^{\star}_{z_{2},\dots,z_{m}}(z_{1}):=F^{\star}(z_{1},z_{2},\dots,z_{m}),$
(119)

where $z_{2},\dots,z_{m}$ are fixed in $\mathbb{C}$, which follows a single-
mode evolution. This function is the product of a Gaussian function with a
univariate polynomial in $z_{1}$, whose coefficients are multivariate
polynomials in $z_{2},\dots,z_{m}$. It may be thus treated formally as a
single-mode stellar function of finite stellar rank with symbolic
coefficients, and we may follow the single-mode derivation of from Sec. 3 to
obtain a formal expression for the output stellar function. Note that the
zeros of a univariate section of a multivariate polynomial are usually
implicit functions of the other variables, thus making more practical the
description based on the coefficients of the polynomial part of the stellar
function, rather than on its zeros.

An RPHQC may therefore be thought of as a sequence of motions of
Calogero–Moser particles corresponding to the zeros of sections of multimode
stellar functions, together with conformal evolutions of Gaussian parameters
and unitary changes of coordinates. While RPHQC form a reduced HQC model, we
show hereafter that their non-adaptive and adaptive versions already capture
various important bosonic quantum computational models.

####  5.2.1 Non-adaptive RPHQC

We first consider the class of non-adaptive RPHQC. We show that this class of
circuits provides a natural classification of infinite-dimensional
subuniversal quantum computational models, based on their computational
complexity and their CV properties, summarised in Table 2. For both continuous
and discrete measurements, we give reductions of non-adaptive RPHQC to
existing subuniversal models, which can either be simulated efficiently
classically (i.e., their output probabilities may be computed in P), or are
hard to sample classically (i.e., estimating multiplicatively a single output
probability up to inverse polynomial precision is $\\#$P-hard, which implies
exact sampling hardness [13]).

The input functions in RPHQC are of the form $P\times G$. A subclass of these
holomorphic functions is given by Gaussians $G$. Considering only Gaussian
functions and unitary operations that leave invariant the set of Gaussian
functions, one obtains restricted bosonic quantum circuits over Gaussians,
which are a subset of RPHQC. The corresponding unitary operations are also
Gaussian unitary operations. Adding projections onto coherent states
(continuous measurement) yields the Gaussian quantum computing model, which
can be efficiently simulated classically [28]. On the other hand, adding
projections onto Fock states (discrete measurement) yields the Gaussian Boson
Sampling model, which is hard to sample classically [14, 15]. However, if no
Gaussian quadratic non-linearity is involved such as squeezing, then by Eqs.
(26) and (27) the circuits reduce to sampling the output photon-number of
coherent states which undergo a passive linear operation, which we call
coherent state sampling. These circuits can be efficiently simulated
classically [13] since coherent states do not become entangled through passive
linear operations. This reflects the necessity of single-mode Gaussian
quadratic non-linearity for the presence of Gaussian entanglement (see Sec.
4.3).

Another important subclass of holomorphic functions of the form $P\times G$ is
given by polynomials $P$. These functions correspond to core states, i.e.,
finite superpositions of Fock basis states. Considering only polynomial
functions and unitary operations that leave invariant the set of polynomials
of each degree, one obtains restricted bosonic quantum circuits (without
squeezing) over polynomials, which are a subset of RPHQC. The corresponding
unitary operations are passive linear operations [48]. Adding projections onto
Fock states (discrete measurement) yields the Boson Sampling model of quantum
computing, which is hard to sample classically [13]. On the other hand, adding
projections onto coherent states (continuous measurement) yields a time-
reversed version of coherent state sampling, which is also easy to simulate
[13] due to the absence of squeezing.

Finally, if RPHQC with continuous measurements do contain squeezing, then they
include Boson Sampling with Gaussian measurements such as unbalanced
heterodyne, which are hard to sample [17, 18, 19]. Moreover, these RPHQC
include CV instantaneous quantum circuits, which are also hard to sample [16],
by replacing input Gottesman–Kitaev–Preskill (GKP) states [9] by arbitrarily
close approximations of finite stellar rank with Lemma 11. Note that in these
cases—since the measurement yields a continuous outcome—hardness of sampling
refers to a probability distribution obtained by binning the continuous
outcome space.

####  5.2.2 Adaptive RPHQC

In this section, we consider adaptive versions of RPHQC, where unitary gates
may depend on the result of intermediate measurements.

In order to compare their computational power to that of DV quantum models, we
consider uniform versions of these circuits by asking that the descriptions of
the input holomorphic function and the unitary evolution are given by the
output of a Turing machine, and that the outcome of the computation is written
to some precision on the tape of a Turing machine. This is done to avoid the
artificial appearance of uncomputable numbers within the model and to provide
a fair comparison with DV models. We denote by
$\textsf{RPHQC}_{\textsf{adap}}^{\textsf{CV}}\textsf{P}$ and
$\textsf{RPHQC}_{\textsf{adap}}^{\textsf{DV}}\textsf{P}$ the corresponding
class of decision problems that can be solved in polynomial time using these
circuits, for continuous and discrete measurements, respectively.

Adaptive continuous measurements leverage the computational power of RPHQC to
universal BQP computations:

######  Theorem 6.

We have the inclusion
$\textsf{BQP}\subseteq\textsf{RPHQC}_{\textsf{adap}}^{\textsf{CV}}\textsf{P}$.

###### Proof.

Bosonic quantum circuits taking as input (approximate) GKP and Gaussian
states, with Gaussian gates and adaptive Gaussian measurements are capable of
universal BQP computations [9, 68]. From the previous section, such circuits
are instances of RPHQC with continuous measurements, since approximate GKP
states of finite stellar rank can be obtained using Lemma 11. Hence,
$\textsf{BQP}\subseteq\textsf{RPHQC}_{\textsf{adap}}^{\textsf{CV}}\textsf{P}$.

∎

Adaptive discrete measurements also leverage the computational power of RPHQC
to universal BQP computations:

######  Theorem 7.

We have the inclusion
$\textsf{BQP}\subseteq\textsf{RPHQC}_{\textsf{adap}}^{\textsf{DV}}\textsf{P}$.

###### Proof.

From the previous section, Boson Sampling circuits are a subset of RPHQC with
discrete measurements. Hence,
$\textsf{BosonP}_{\textsf{adap}}\subseteq\textsf{RPHQC}_{\textsf{adap}}^{\textsf{DV}}\textsf{P}$,
where $\textsf{BosonP}_{\textsf{adap}}$ is the class of problems that can be
solved in polynomial time using Boson Sampling circuits with adaptive Fock
basis measurements, defined in [13], where it is shown that the equality
$\textsf{BosonP}_{\textsf{adap}}=\textsf{BQP}$ holds, based on the
Knill–Laflamme–Milburn scheme [65]. Hence,
$\textsf{BQP}\subseteq\textsf{RPHQC}_{\textsf{adap}}^{\textsf{DV}}\textsf{P}$.

∎

Note that the model of cluster state quantum computing [107, 66] is also a
subclass of RPHQC with adaptive discrete measurements, where the input is a
Gaussian function instead.

Summarising the results obtained in this section, HQC circuits provide a
complex analytic description of bosonic quantum computations, and its subclass
of rank-preserving circuits capture existing subuniversal and universal
bosonic quantum computing schemes, with non-adaptive and adaptive (discrete
and continuous) measurements, respectively.

Appendix

##  Appendix A Multi-index notations

For $m\in\mathbb{N}^{*}$, $\bm{p}=(p_{1},\dots,p_{m})\in\mathbb{N}^{m}$,
$\bm{q}=(q_{1},\dots,q_{m})\in\mathbb{N}^{m}$, and
$\bm{z}=(z_{1},\dots,z_{m})\in\mathbb{C}^{m}$, we write:

$\displaystyle\bm{0}$ $\displaystyle=(0,\dots,0)$ (120)
$\displaystyle\bm{p}+\bm{q}$ $\displaystyle=(p_{1}+q_{1},\dots,p_{m}+q_{m})$
$\displaystyle|\bm{p}|$ $\displaystyle=p_{1}+\dots+p_{m}$
$\displaystyle\bm{p}!$ $\displaystyle=p_{1}!\dots p_{m}!$
$\displaystyle\ket{\bm{p}}$
$\displaystyle=\ket{p_{1}}\otimes\dots\otimes\ket{p_{m}}$
$\displaystyle\bm{p}\leq\bm{q}$ $\displaystyle\Leftrightarrow\forall
i\in\\{1,\dots,m\\},\;p_{i}\leq q_{i}$ $\displaystyle\bm{z}^{*}$
$\displaystyle=(z_{1}^{*},\dots,z_{m}^{*})$ $\displaystyle-\bm{z}$
$\displaystyle=(-z_{1},\dots,-z_{m})$ $\displaystyle|\bm{z}|^{2}$
$\displaystyle=|z_{1}|^{2}+\dots+|z_{m}|^{2}$ $\displaystyle\bm{z}^{\bm{p}}$
$\displaystyle=z_{1}^{p_{1}}\cdots z_{m}^{p_{m}}$ $\displaystyle d^{2m}\bm{z}$
$\displaystyle=d\Re{(z_{1})}d\Im{(z_{1})}\dots d\Re{(z_{m})}d\Im{(z_{m})}$
$\displaystyle\partial_{\bm{z}}$
$\displaystyle=(\partial_{z_{1}},\dots,\partial_{z_{m}}).$

##  Appendix B Resolution of the classical Calogero–Moser model

In this section, we briefly review the resolution of the classical
Calogero–Moser model [36] by the Olshanetsky–Perelomov projection method [54].
The $n$-body Calogero–Moser Hamiltonian is given by:

$H^{CM}=\frac{1}{2}\sum_{k=1}^{n}\left(p_{k}^{2}+\omega^{2}q_{k}^{2}\right)+\frac{1}{2}g^{2}\sum_{k=1}^{n}\sum_{j\neq
k}\frac{1}{(q_{k}-q_{j})^{2}}.$ (121)

We treat here the case $\omega=0$ (isolated Calogero–Moser model) and refer
the reader to [55] for the general case. The equations of motions are given by
Hamilton’s equations:

$\frac{dq_{k}(t)}{dt}=p_{k}\quad\text{and}\quad\frac{dp_{k}(t)}{dt}=2g^{2}\sum_{j\neq
k}\frac{1}{(q_{k}-q_{j})^{3}},\quad\forall k\in\\{1,\dots,n\\}.$ (122)

These equations of motions may be solved by writing them as the projection of
a free motion in a higher-dimensional space. More precisely, consider the free
motion of a hermitian $n\times n$ matrix $\Lambda$, described by the equation

$\frac{d^{2}\Lambda}{dt^{2}}=0.$ (123)

The solution is given by $\Lambda(t)=At+B$, where $A$ and $B$ are hermitian
matrices. Because $\Lambda$ is hermitian, we may write

$\Lambda(t)=U(t)Q(t)U^{\dagger}(t),$ (124)

where $U$ is unitary and $Q$ is diagonal. Differentiating this equation, we
obtain

$\frac{d\Lambda(t)}{dt}=\frac{dU(t)}{dt}Q(t)U^{\dagger}(t)+U(t)\frac{dQ(t)}{dt}U^{\dagger}(t)+U(t)Q(t)\frac{dU^{\dagger}(t)}{dt}=A.$
(125)

Now let

$\displaystyle M(t)$ $\displaystyle:=-iU^{\dagger}(t)\frac{dU(t)}{dt}$ (126)
$\displaystyle L(t)$ $\displaystyle:=\frac{dQ(t)}{dt}+[M(t),Q(t)].$

We show that these matrices form a Lax pair, i.e., they satisfy the equation
$i\frac{dL}{dt}=[M,L]$. Since $U$ is unitary, we have
$U^{\dagger}(t)U(t)=\mathbb{I}$, and thus
$\frac{dU^{\dagger}(t)}{dt}U(t)=-U^{\dagger}(t)\frac{dU(t)}{dt}$. Hence,

$\displaystyle U(t)L(t)U^{\dagger}(t)$
$\displaystyle=U(t)\frac{dQ(t)}{dt}U^{\dagger}(t)+iU(t)M(t)Q(t)U^{\dagger}(t)-iU(t)Q(t)M(t)U^{\dagger}(t)$
(127)
$\displaystyle=U(t)\frac{dQ(t)}{dt}U^{\dagger}(t)+\frac{dU(t)}{dt}Q(t)U^{\dagger}(t)-U(t)Q(t)U^{\dagger}(t)\frac{dU(t)}{dt}U^{\dagger}(t)$
$\displaystyle=U(t)\frac{dQ(t)}{dt}U^{\dagger}(t)+\frac{dU(t)}{dt}Q(t)U^{\dagger}(t)+U(t)Q(t)\frac{dU^{\dagger}(t)}{dt}$
$\displaystyle=A,$

where we used Eq. (125) in the last line. As a consequence,
$L(t)=U^{\dagger}(t)AU(t)$, so the eigenvalues of $L$ are conserved quantities
of the evolution. We set $U(0)=0$ so that $L(0)=A$, $\Lambda(0)=B=Q(0)$, and
thus $\Lambda(t)=L(0)t+Q(0)$. Differentiating again Eq. (127) yields

$\displaystyle\frac{dU(t)}{dt}L(t)U^{\dagger}(t)$
$\displaystyle+U(t)\frac{dL(t)}{dt}U^{\dagger}(t)+U(t)L(t)\frac{dU^{\dagger}(t)}{dt}=0$
(128)
$\displaystyle\Leftrightarrow\frac{dL(t)}{dt}+U^{\dagger}(t)\frac{dU(t)}{dt}L(t)+L(t)\frac{dU^{\dagger}(t)}{dt}U(t)=0$
$\displaystyle\Leftrightarrow\frac{dL(t)}{dt}+\left[U^{\dagger}(t)\frac{dU(t)}{dt},L(t)\right]=0$
$\displaystyle\Leftrightarrow\frac{dL(t)}{dt}+i[M(t),L(t)]=0,$

where we used $\frac{dU^{\dagger}(t)}{dt}U(t)=-U^{\dagger}(t)\frac{dU(t)}{dt}$
in the third line. Hence the matrices $L$ and $M$ indeed form a Lax pair, and
the eigenvalues of $L$ are conserved quantities.

Now setting $Q_{jk}=\delta_{jk}q_{j}$, $M_{jk}=g\big{[}\delta_{jk}\sum_{l\neq
j}\frac{1}{(q_{j}-q_{l})^{2}}-(1-\delta_{jk})\frac{1}{(q_{j}-q_{k})^{2}}\big{]}$,
and $L_{jk}=\delta_{jk}p_{j}+(1-\delta_{jk})\frac{ig}{q_{j}-q_{k}}$, where
$\delta$ is the Kronecker symbol, ensures that Eq. (126) is satisfied and that
Eq. (128) is equivalent to Eq. (122) [53]. By construction, the canonical
variables $q_{k}$ are the eigenvalues of the matrix $\Lambda(t)=L(0)t+Q(0)$
given by:

$\Lambda_{jk}(t)=\delta_{jk}q_{j}(0)+\delta_{jk}p_{j}(0)t+(1-\delta_{jk})\frac{igt}{q_{j}-q_{k}}.$
(129)

Moreover, the eigenvalues of the matrix $L$ are conserved quantities of the
evolution. A direct consequence is that the scattering process can be
described simply by a permutation of the trajectories. This was proved, e.g.,
in [55] in the case of particles on a line, and we reproduce the argument
below in the case where $q_{k}\in\mathbb{C}$: due to the repulsive potential,
at large times the particles follow a free motion. Writing

$q_{k}(t)\sim p_{k}^{\pm}t+q_{k}^{\pm},\quad t\rightarrow\pm\infty,$ (130)

the asymptotic momenta $p_{k}^{\pm}=p_{k}(\pm\infty)$ are the eigenvalues of
$L(\pm\infty)$ and are thus conserved. Hence, there exists a permutation
$\sigma\in\mathcal{S}_{n}$ with permutation matrix $P_{\sigma}$ such that
$\text{Diag}(p_{1}^{+},\dots,p_{n}^{+})=P_{\sigma}\text{Diag}(p_{1}^{-},\dots,p_{n}^{-})P_{\sigma}^{T}$,
i.e., $L(+\infty)=P_{\sigma}L(-\infty)P_{\sigma}^{T}$. With Eq. (127) we
obtain $P_{\sigma}=U^{\dagger}(+\infty)U(-\infty)$. Now we have:

$\displaystyle Q(t)$ $\displaystyle=U^{\dagger}(t)\Lambda(t)U(t)$ (131)
$\displaystyle=U^{\dagger}(t)(At+B)U(t)$
$\displaystyle=L(t)t+U^{\dagger}(t)BU(t)$
$\displaystyle\sim\text{Diag}_{k}(p_{k}^{\pm}t+q_{k}^{\pm}),\quad
t\rightarrow\pm\infty.$

Since
$L_{jk}(t)=\delta_{jk}p_{j}(t)+(1-\delta_{jk})\frac{ig}{q_{j}(t)-q_{k}(t)}$ we
have $L(t)\sim\text{Diag}(p_{1}^{\pm},\dots,p_{n}^{\pm})$ when
$t\rightarrow\pm\infty$, and thus
$\text{Diag}(q_{1}^{\pm},\dots,q_{n}^{\pm})=U^{\dagger}(\pm\infty)BU(\pm\infty)$.
Hence,
$\text{Diag}(q_{1}^{+},\dots,q_{n}^{+})=P_{\sigma}\text{Diag}(q_{1}^{-},\dots,q_{n}^{-})P_{\sigma}^{T}$.
This shows that

$\text{Diag}_{k}(p_{k}^{+}t+q_{k}^{+})=P_{\sigma}\text{Diag}_{k}(p_{k}^{-}t+q_{k}^{-})P_{\sigma}^{T},$
(132)

that is, the scattering process
$(q_{1}^{-},p_{1}^{-}),\dots,(q_{n}^{-},p_{n}^{-})\mapsto(q_{1}^{+},p_{1}^{+}),\dots,(q_{n}^{+},p_{n}^{+})$
reduces to a permutation $\sigma$ of the trajectories, with asymptotic momenta
$p_{k}^{\pm}$ and asymptotic offsets $q_{k}^{\pm}$.

We can extend this analysis to relate the asymptotic parameters to the initial
conditions $q_{k}(0)$ and $p_{k}(0)$: from

$\displaystyle\text{Diag}(q_{1}^{\pm},\dots,q_{n}^{\pm})$
$\displaystyle=U^{\dagger}(\pm\infty)BU(\pm\infty)$ (133)
$\displaystyle=U^{\dagger}(\pm\infty)Q(0)U(\pm\infty)$
$\displaystyle=U^{\dagger}(\pm\infty)\text{Diag}(q_{1}(0),\dots,q_{n}(0))U(\pm\infty),$

we deduce that the asymptotic offsets $q_{k}^{\pm}$ are given by permutations
of the initial positions $q_{k}(0)$. For the asymptotic momenta, we have

$\displaystyle\text{Diag}(p_{1}^{\pm},\dots,p_{n}^{\pm})$
$\displaystyle=L(\pm\infty)$ (134)
$\displaystyle=U^{\dagger}(\pm\infty)L(0)U(\pm\infty),$

so the asymptotic momenta $p_{k}^{\pm}$ are given by permutations of the
eigenvalues of the matrix

$L(0)=\left(\delta_{jk}p_{j}(0)+(1-\delta_{jk})\frac{ig}{q_{j}(0)-q_{k}(0)}\right)_{1\leq
j,k\leq n}.$ (135)

##  Appendix C Gaussian dynamics of the stellar representation

In this section we detail the proofs of the results obtained in Sec. 3.

###  C.1 Decoupling of the zeros and Gaussian parameters evolutions

We first give the proof of Theorem 1, which we recall below. From Sec. 3, the
evolution of the stellar function is governed by

$i\partial_{t}F^{\star}(z,t)=H(z,\partial_{z})F^{\star}(z,t),$ (136)

for all $z\in\mathbb{C}$ and all $t\geq 0$, and a generic Gaussian Hamiltonian
satisfies:

$-iH_{\alpha,\xi,\varphi}(z,\partial_{z})=\frac{1}{2}\xi
z^{2}-\frac{1}{2}\xi^{*}\partial^{2}_{z}+i\varphi z\partial_{z}+\alpha
z-\alpha^{*}\partial_{z},$ (137)

for $\alpha,\xi\in\mathbb{C}$ and $\varphi\in\mathbb{R}$.

######  Theorem 1.

Let
$F^{\star}(z,t)=\prod_{k=1}^{n}(z-\lambda_{k}(t))e^{-\frac{1}{2}a(t)z^{2}+b(t)z+c(t)}$.
Assuming that the initial zeros $\\{\lambda_{k}(0)\\}_{k=1\dots n}$ are
simple, the evolution of the stellar function under a Gaussian Hamiltonian
$\hat{H}=H_{\alpha,\xi,\varphi}(\hat{a}^{\dagger},\hat{a})$ may be recast as
the following dynamical system:

$\begin{cases}\frac{da(t)}{dt}=\xi^{*}a^{2}(t)+2i\varphi a(t)-\xi,\\\
\frac{db(t)}{dt}=(i\varphi+\xi^{*}a(t))b(t)+\alpha+\alpha^{*}a(t),\\\
\frac{dc(t)}{dt}=\frac{1}{2}\xi^{*}a(t)-\frac{1}{2}\xi
b^{2}(t)-\alpha^{*}b(t)+n(\xi^{*}a(t)+i\varphi)\\\
\frac{d^{2}\lambda_{k}(t)}{dt^{2}}=(|\xi|^{2}-\varphi^{2})\lambda_{k}(t)+(\xi^{*}\alpha-i\varphi\alpha^{*})-2\xi^{*2}\sum_{j\neq
k}\frac{1}{(\lambda_{k}(t)-\lambda_{j}(t))^{3}},\quad\forall
k\in\\{1,\dots,n\\}.\end{cases}$ (138)

The equations for the Gaussian parameters and the zeros are decoupled, and the
relation between the root system $\\{\lambda_{k}\\}_{k=1\dots n}$ and the
Gaussian parameters $a,b,c$ is specified by the initial conditions:

$\frac{d\lambda_{k}(t)}{dt}\bigg{|}_{t=0}=-\left(\xi^{*}a(0)+i\varphi\right)\lambda_{k}(0)+\xi^{*}b(0)+\alpha^{*}+\xi^{*}\sum_{j\neq
k}\frac{1}{\lambda_{k}(0)-\lambda_{j}(0)},\quad\forall k\in\\{1,\dots,n\\}.$
(139)

###### Proof.

Combining Eqs. (136) and (137) the evolution of the stellar function is
governed by the partial differential equation:

$\partial_{t}[F^{\star}(z,t)]=\left(\frac{1}{2}\xi
z^{2}-\frac{1}{2}\xi^{*}\partial^{2}_{z}+i\varphi z\partial_{z}+\alpha
z-\alpha^{*}\partial_{z}\right)[F^{\star}(z,t)].$ (140)

Moreover, we have:

$\displaystyle F^{\star}(z,t)$
$\displaystyle=P(z,t)G(z,t)=\prod_{k=1}^{n}(z-\lambda_{k}(t))e^{-\frac{1}{2}a(t)z^{2}+b(t)z+c(t)}$
(141) $\displaystyle\partial_{t}F^{\star}(z,t)$
$\displaystyle=F^{\star}(z,t)\left[-\frac{1}{2}z^{2}\frac{da(t)}{dt}+z\frac{db(t)}{dt}+\frac{dc(t)}{dt}-\sum_{k=1}^{n}\frac{1}{z-\lambda_{k}(t)}\frac{d\lambda_{k}(t)}{dt}\right]$
$\displaystyle\partial_{z}F^{\star}(z,t)$
$\displaystyle=F^{\star}(z,t)\left[-a(t)z+b(t)+\sum_{k=1}^{n}\frac{1}{z-\lambda_{k}(t)}\right]$
$\displaystyle\partial_{z}^{2}F^{\star}(z,t)$
$\displaystyle=F^{\star}(z,t)\left[\left(-a(t)z+b(t)+\sum_{k=1}^{n}\frac{1}{z-\lambda_{k}(t)}\right)^{2}-a(t)-\sum_{k=1}^{n}\frac{1}{(z-\lambda_{k}(t))^{2}}\right].$

Plugging these expressions in Eq. (136) and simplifying the Gaussian functions
on each side, we obtain the partial differential equation:

$\displaystyle\left[-\frac{1}{2}z^{2}\frac{da(t)}{dt}+z\frac{db(t)}{dt}+\frac{dc(t)}{dt}-\sum_{k=1}^{n}\frac{1}{z-\lambda_{k}(t)}\frac{d\lambda_{k}(t)}{dt}\right]P(z,t)$
(142)
$\displaystyle\quad\quad\quad=\Bigg{[}\left(\frac{1}{2}\xi-\frac{1}{2}\xi^{*}a^{2}(t)-i\varphi
a(t)\right)z^{2}$
$\displaystyle\quad\quad\quad\quad+\left(\xi^{*}a(t)b(t)+i\varphi
b(t)+\alpha+\alpha^{*}a(t)\right)z$
$\displaystyle\quad\quad\quad\quad+\left(\frac{1}{2}\xi^{*}a(t)-\frac{1}{2}\xi
b^{2}(t)-\alpha^{*}b(t)\right)$
$\displaystyle\quad\quad\quad\quad+\sum_{k=1}^{n}\frac{1}{z-\lambda_{k}(t)}\left(\left(\xi^{*}a(t)+i\varphi\right)z-\xi^{*}b(t)-\alpha^{*}-\frac{1}{2}\xi^{*}\sum_{j\neq
k}\frac{1}{z-\lambda_{j}(t)}\right)\Bigg{]}P(z,t).$

This equation is an equality between two polynomials in $z$ with time-
dependent coefficients. The leading coefficients give

$\frac{da(t)}{dt}=\xi^{*}a^{2}(t)+2i\varphi a(t)-\xi.$ (143)

Removing these coefficients on both sides the differential equation (142)
rewrites

$\displaystyle\left[z\frac{db(t)}{dt}+\frac{dc(t)}{dt}-\sum_{k=1}^{n}\frac{1}{z-\lambda_{k}(t)}\frac{d\lambda_{k}(t)}{dt}\right]P(z,t)$
(144) $\displaystyle\quad\quad\quad=\Bigg{[}\left(\xi^{*}a(t)b(t)+i\varphi
b(t)+\alpha+\alpha^{*}a(t)\right)z$
$\displaystyle\quad\quad\quad\quad+\left(\frac{1}{2}\xi^{*}a(t)-\frac{1}{2}\xi
b^{2}(t)-\alpha^{*}b(t)\right)$
$\displaystyle\quad\quad\quad\quad+\sum_{k=1}^{n}\frac{1}{z-\lambda_{k}(t)}\left(\left(\xi^{*}a(t)+i\varphi\right)z-\xi^{*}b(t)-\alpha^{*}-\frac{1}{2}\xi^{*}\sum_{j\neq
k}\frac{1}{z-\lambda_{j}(t)}\right)\Bigg{]}P(z,t),$

and the leading coefficients give

$\frac{db(t)}{dt}=(i\varphi+\xi^{*}a(t))b(t)+\alpha+\alpha^{*}a(t).$ (145)

Once again, removing these coefficients on both sides, the differential
equation (144) rewrites

$\displaystyle\left[\frac{dc(t)}{dt}-\sum_{k=1}^{n}\frac{1}{z-\lambda_{k}(t)}\frac{d\lambda_{k}(t)}{dt}\right]P(z,t)$
(146)
$\displaystyle\quad\quad\quad=\Bigg{[}\left(\frac{1}{2}\xi^{*}a(t)-\frac{1}{2}\xi
b^{2}(t)-\alpha^{*}b(t)\right)$
$\displaystyle\quad\quad\quad\quad+\sum_{k=1}^{n}\frac{1}{z-\lambda_{k}(t)}\left(\left(\xi^{*}a(t)+i\varphi\right)z-\xi^{*}b(t)-\alpha^{*}-\frac{1}{2}\xi^{*}\sum_{j\neq
k}\frac{1}{z-\lambda_{j}(t)}\right)\Bigg{]}P(z,t).$

The leading coefficient on the left hand side is $\frac{dc(t)}{dt}$, while the
leading coefficient on the right hand side is

$\left(\frac{1}{2}\xi^{*}a(t)-\frac{1}{2}\xi
b^{2}(t)-\alpha^{*}b(t)\right)+\sum_{k=1}^{n}\left(\xi^{*}a(t)+i\varphi\right).$
(147)

We thus obtain

$\frac{dc(t)}{dt}=\frac{1}{2}\xi^{*}a(t)-\frac{1}{2}\xi
b^{2}(t)-\alpha^{*}b(t)+n(\xi^{*}a(t)+i\varphi).$ (148)

Finally, removing $(\frac{1}{2}\xi^{*}a(t)-\frac{1}{2}\xi
b^{2}(t)-\alpha^{*}b(t))z^{n}$ on both sides, the differential equation (146)
rewrites

$\displaystyle\left[n(\xi^{*}a(t)+i\varphi)-\sum_{k=1}^{n}\frac{1}{z-\lambda_{k}(t)}\frac{d\lambda_{k}(t)}{dt}\right]P(z,t)$
(149)
$\displaystyle\quad\quad\quad=\Bigg{[}\sum_{k=1}^{n}\frac{1}{z-\lambda_{k}(t)}\left(\left(\xi^{*}a(t)+i\varphi\right)z-\xi^{*}b(t)-\alpha^{*}-\frac{1}{2}\xi^{*}\sum_{j\neq
k}\frac{1}{z-\lambda_{j}(t)}\right)\Bigg{]}P(z,t),$

or equivalently:

$\displaystyle-n\left(\xi^{*}a(t)+i\varphi\right)\prod_{k=1}^{n}(z-\lambda_{k}(t))+\sum_{k=1}^{n}\frac{d\lambda_{k}(t)}{dt}\prod_{j\neq
k}(z-\lambda_{j}(t))$ (150)
$\displaystyle\quad\quad\quad=\left[-\left(\xi^{*}a(t)+i\varphi\right)z+\xi^{*}b(t)+\alpha^{*}\right]\sum_{k=1}^{n}\prod_{j\neq
k}(z-\lambda_{j}(t))+\frac{1}{2}\xi^{*}\sum_{k=1}^{n}\sum_{j\neq
k}\prod_{\begin{subarray}{c}l\neq j\\\ l\neq
k\end{subarray}}(z-\lambda_{l}(t)).$

Let $k\in\\{1,\dots,n\\}$. Under the assumption that $\lambda_{k}(t)$ is a
simple zero of $z\mapsto P(z,t)$, dividing the above equation by $\prod_{j\neq
k}(z-\lambda_{j}(t))$ and setting $z=\lambda_{k}(t)$ yields

$\frac{d\lambda_{k}(t)}{dt}=-\left(\xi^{*}a(t)+i\varphi\right)\lambda_{k}(t)+\xi^{*}b(t)+\alpha^{*}+\xi^{*}\sum_{j\neq
k}\frac{1}{\lambda_{k}(t)-\lambda_{j}(t)}.$ (151)

In particular, for all $k\in\\{1,\dots,n\\}$ we have

$\frac{d\lambda_{k}(t)}{dt}\bigg{|}_{t=0}=-\left(\xi^{*}a(0)+i\varphi\right)\lambda_{k}(0)+\xi^{*}b(0)+\alpha^{*}+\xi^{*}\sum_{j\neq
k}\frac{1}{\lambda_{k}(0)-\lambda_{j}(0)}.$ (152)

At this point we have the system:

$\begin{cases}\frac{da(t)}{dt}=\xi^{*}a^{2}(t)+2i\varphi a(t)-\xi,\\\
\frac{db(t)}{dt}=(i\varphi+\xi^{*}a(t))b(t)+\alpha+\alpha^{*}a(t),\\\
\frac{dc(t)}{dt}=\frac{1}{2}\xi^{*}a(t)-\frac{1}{2}\xi
b^{2}(t)-\alpha^{*}b(t)+n(\xi^{*}a(t)+i\varphi),\\\
\frac{d\lambda_{k}(t)}{dt}=-\left(\xi^{*}a(t)+i\varphi\right)\lambda_{k}(t)+\xi^{*}b(t)+\alpha^{*}+\xi^{*}\sum_{j\neq
k}\frac{1}{\lambda_{k}(t)-\lambda_{j}(t)},\quad\forall
k\in\\{1,\dots,n\\}.\end{cases}$ (153)

We are left with showing that this differential system indeed decouples by
taking the derivative of the last equation:

$\displaystyle\frac{d^{2}\lambda_{k}(t)}{dt^{2}}$
$\displaystyle=-\xi^{*}\frac{da(t)}{dt}\lambda_{k}(t)-(\xi^{*}a(t)+i\varphi)\frac{d\lambda_{k}(t)}{dt}+\xi^{*}\frac{db(t)}{dt}+\xi^{*}\sum_{j\neq
k}\frac{\frac{d\lambda_{j}(t)}{dt}-\frac{d\lambda_{k}(t)}{dt}}{(\lambda_{k}(t)-\lambda_{j}(t))^{2}}$
(154)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=-\xi^{*}(\xi^{*}a^{2}(t)+2i\varphi-\xi)\lambda_{k}(t)-(\xi^{*}a(t)+i\varphi)\left[-(\xi^{*}a(t)+i\varphi)\lambda_{k}(t)+\xi^{*}b(t)+\alpha^{*}+\xi^{*}\sum_{j\neq
k}\frac{1}{\lambda_{k}(t)-\lambda_{j}(t)}\right]$
$\displaystyle\\!\\!+\xi^{*}(i\varphi
b(t)+\xi^{*}a(t)b(t)+\alpha+\alpha^{*}a(t))$
$\displaystyle\\!\\!+\xi^{*}\sum_{j\neq
k}\frac{1}{(\lambda_{k}(t)-\lambda_{j}(t))^{2}}\left[(\xi^{*}a(t)+i\varphi)(\lambda_{k}(t)-\lambda_{j}(t))+\xi^{*}\sum_{l\neq
j}\frac{1}{\lambda_{j}(t)-\lambda_{l}(t)}-\xi^{*}\sum_{l\neq
k}\frac{1}{\lambda_{k}(t)-\lambda_{l}(t)}\right]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=(|\xi|^{2}\\!-\\!\varphi^{2})\lambda_{k}(t)\\!+\\!(\xi^{*}\alpha\\!-\\!i\varphi\alpha^{*})+\xi^{*2}\sum_{j\neq
k}\frac{1}{(\lambda_{k}(t)-\lambda_{j}(t))^{2}}\\!\\!\left[\frac{2}{\lambda_{j}(t)\\!-\\!\lambda_{k}(t)}\\!+\\!\sum_{\begin{subarray}{c}l\neq
j\\\ l\neq
k\end{subarray}}\frac{\lambda_{k}(t)\\!-\\!\lambda_{j}(t)}{(\lambda_{j}(t)\\!-\\!\lambda_{l}(t))(\lambda_{k}(t)\\!-\\!\lambda_{l}(t))}\right]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=(|\xi|^{2}\\!-\\!\varphi^{2})\lambda_{k}(t)\\!+\\!(\xi^{*}\alpha\\!-\\!i\varphi\alpha^{*})-2\xi^{*2}\sum_{j\neq
k}\frac{1}{(\lambda_{k}(t)\\!-\\!\lambda_{j}(t))^{3}}\\!+\\!\sum_{\begin{subarray}{c}j,l\neq
k\\\ j\neq
l\end{subarray}}\frac{1}{(\lambda_{k}(t)\\!-\\!\lambda_{l}(t))(\lambda_{k}(t)\\!-\\!\lambda_{j}(t))(\lambda_{j}(t)\\!-\\!\lambda_{l}(t))}$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=(|\xi|^{2}\\!-\\!\varphi^{2})\lambda_{k}(t)\\!+\\!(\xi^{*}\alpha\\!-\\!i\varphi\alpha^{*})-2\xi^{*2}\sum_{j\neq
k}\frac{1}{(\lambda_{k}(t)-\lambda_{j}(t))^{3}},$

for all $k\in\\{1,\dots,n\\}$.

∎

###  C.2 Displacement, phase-shift, shearing and squeezing evolutions

We first prove Lemma 1 from the main text which we recall below:

######  Lemma 1.

Let
$F^{\star}(z,t)=\prod_{k=1}^{n}(z-\lambda_{k}(t))e^{-\frac{1}{2}a(t)z^{2}+b(t)z+c(t)}$
and let $\alpha\in\mathbb{C}$. The evolution under the displacement
Hamiltonian $\hat{H}_{\alpha}^{D}$ is given by

$\begin{cases}a(t)=a(0),\\\ b(t)=(\alpha+\alpha^{*}a(0))t+b(0),\\\
c(t)=\frac{1}{2}(\alpha^{*2}a(0)-|\alpha|^{2})t^{2}+b(0)t+c(0),\\\
\lambda_{k}(t)=\alpha^{*}t+\lambda(0),\quad\forall
k\in\\{1,\dots,n\\}.\end{cases}$ (155)

###### Proof.

Setting $\xi=\varphi=0$ in Theorem 1 yields:

$\begin{cases}\frac{da(t)}{dt}=0,\\\
\frac{db(t)}{dt}=\alpha+\alpha^{*}a(t),\\\ \frac{dc(t)}{dt}=-\alpha^{*}b(t)\\\
\frac{d^{2}\lambda_{k}(t)}{dt^{2}}=0,\quad\forall
k\in\\{1,\dots,n\\},\end{cases}$ (156)

with the initial conditions

$\frac{d\lambda_{k}(t)}{dt}\bigg{|}_{t=0}=\alpha^{*},\quad\forall
k\in\\{1,\dots,n\\}.$ (157)

Integrating the first equation of Eq. (156) gives $a(t)=a(0)$. Plugging this
in the second equation and integrating gives
$b(t)=(\alpha+\alpha^{*}a(0))t+b(0)$. With the third equation we obtain
$c(t)=\frac{1}{2}(\alpha^{*2}a(0)-|\alpha|^{2})t^{2}+b(0)t+c(0)$. Finally, the
fourth equation together with the initial conditions gives
$\lambda_{k}(t)=\alpha^{*}t+\lambda(0)$ for all $k\in\\{1,\dots,n\\}$.

∎

We now prove Lemma 2 from the main text which we recall below:

######  Lemma 2.

Let
$F^{\star}(z,t)=\prod_{k=1}^{n}(z-\lambda_{k}(t))e^{-\frac{1}{2}a(t)z^{2}+b(t)z+c(t)}$
and let $\varphi\in[0,2\pi]$. The evolution under the phase-shift Hamiltonian
$\hat{H}_{\varphi}^{R}$ is given by

$\begin{cases}a(t)=e^{2i\varphi t}a(0),\\\ b(t)=e^{i\varphi t}b(0),\\\
c(t)=e^{in\varphi t}c(0),\\\ \lambda_{k}(t)=e^{-i\varphi
t}\lambda(0),\quad\forall k\in\\{1,\dots,n\\}.\end{cases}$ (158)

###### Proof.

Setting $\alpha=\xi=0$ in Theorem 1 yields:

$\begin{cases}\frac{da(t)}{dt}=2i\varphi a(t),\\\ \frac{db(t)}{dt}=i\varphi
b(t),\\\ \frac{dc(t)}{dt}=in\varphi\\\
\frac{d^{2}\lambda_{k}(t)}{dt^{2}}=-\varphi^{2}\lambda_{k}(t),\quad\forall
k\in\\{1,\dots,n\\},\end{cases}$ (159)

with the initial conditions

$\frac{d\lambda_{k}(t)}{dt}\bigg{|}_{t=0}=-i\varphi\lambda_{k}(0),\quad\forall
k\in\\{1,\dots,n\\}.$ (160)

We obtain $a(t)=e^{2i\varphi t}a(0)$, $b(t)=e^{i\varphi t}b(0)$,
$c(t)=e^{in\varphi t}c(0)$ and $\lambda_{k}(t)=e^{-i\varphi t}\lambda(0)$ for
all $k\in\\{1,\dots,n\\}$.

∎

We now prove Lemma 3 from the main text which we recall below:

######  Lemma 3.

Let
$F^{\star}(z,t)=\prod_{k=1}^{n}(z-\lambda_{k}(t))e^{-\frac{1}{2}a(t)z^{2}+b(t)z+c(t)}$
and let $s\in\mathbb{R}$. The evolution under the shearing Hamiltonian
$\hat{H}_{s}^{P}$ is given by

$\begin{cases}a(t)=\frac{a(0)-ist(1-a(0))}{1-ist(1-a(0))},\\\
b(t)=\frac{b(0)}{1-ist(1-a(0))},\\\
c(t)=c(0)-\frac{ist}{2}-\left(n+\frac{1}{2}\right)\log(1-ist(1-a(0)))-\frac{b^{2}(0)}{2(1-a(0))(1-ist(1-a(0))},\end{cases}$
(161)

and the zeros $\lambda_{k}(t)$ are the eigenvalues of the matrix $\Lambda(t)$
defined as

$\Lambda_{kl}(t)=\begin{cases}\lambda_{k}(0)-ist\left[(1-a(0))\lambda_{k}(0)+b(0)+\sum_{j\neq
k}\frac{1}{\lambda_{k}(0)-\lambda_{j}(0)}\right],\quad&k=l\\\
\frac{ist}{\lambda_{l}(0)-\lambda_{k}(0)},\quad&k\neq l.\end{cases}$ (162)

###### Proof.

Recall that
$\hat{H}_{s}^{P}=\hat{H}_{is}^{S}+\hat{H}_{s}^{R}-\frac{s}{2}\hat{1}$. Setting
$\alpha=0$, $\varphi=s$ and $\xi=is$ in Theorem 1 yields:

$\begin{cases}\frac{da(t)}{dt}=-is(1-a(t))^{2},\\\
\frac{db(t)}{dt}=is(1-a(t))b(t),\\\
\frac{dc(t)}{dt}=-\frac{is}{2}a(t)-\frac{is}{2}b^{2}(t)+ins(1-a(t))\\\
\frac{d^{2}\lambda_{k}(t)}{dt^{2}}=2s\sum_{j\neq
k}\frac{1}{(\lambda_{k}(t)-\lambda_{j}(t))^{3}},\quad\forall
k\in\\{1,\dots,n\\}.\end{cases}$ (163)

with the initial conditions

$\frac{d\lambda_{k}(t)}{dt}\bigg{|}_{t=0}=is(a(0)-1)\lambda_{k}(0)-isb(0)-is\sum_{j\neq
k}\frac{1}{\lambda_{k}(0)-\lambda_{j}(0)},\quad\forall k\in\\{1,\dots,n\\}.$
(164)

The first differential equation in Eq. (163) is separable and can be solved
directly:

$a(t)=\frac{a(0)-ist(1-a(0))}{1-ist(1-a(0))}.$ (165)

Plugging this expression in the second differential equation gives another
separable equation which yields:

$b(t)=\frac{b(0)}{1-ist(1-a(0))}.$ (166)

With these expressions for $a(t)$ and $b(t)$, integrating the third
differential equation gives

$c(t)=c(0)-\frac{ist}{2}-\left(n+\frac{1}{2}\right)\log(1-ist(1-a(0)))-\frac{b^{2}(0)}{2(1-a(0))(1-ist(1-a(0))}.$
(167)

Finally, the remaining dynamical system for the zeros is the complex version
of the Calogero–Moser system [52], which may be solved using the
Olshanetsky–Perelomov projection method (see Appendix B). We obtain that the
zeros $\\{\lambda_{k}\\}_{k=1,\dots,n}$ are given by the eigenvalues of a
matrix $\Lambda(t)$ with analytical time-dependent coefficients defined as:

$\displaystyle\Lambda_{kl}(t)$
$\displaystyle:=\begin{cases}\lambda_{k}(0)+\frac{d\lambda(0)}{dt}t,\quad&k=l\\\
\frac{ist}{\lambda_{l}(0)-\lambda_{k}(0)},\quad&k\neq l\end{cases}$ (168)
$\displaystyle=\begin{cases}\lambda_{k}(0)-ist\left[(1-a(0))\lambda_{k}(0)+b(0)+\sum_{j\neq
k}\frac{1}{\lambda_{k}(0)-\lambda_{j}(0)}\right],\quad&k=l\\\
\frac{ist}{\lambda_{l}(0)-\lambda_{k}(0)},\quad&k\neq l,\end{cases}$

where we used Eq. (164) in the second line. ∎

For completeness, we also derive the evolution under squeezing:

######  Lemma 4.

Let
$F^{\star}(z,t)=\prod_{k=1}^{n}(z-\lambda_{k}(t))e^{-\frac{1}{2}a(t)z^{2}+b(t)z+c(t)}$
and let $\xi=re^{i\theta}\in\mathbb{C}$. Let also
$A:=\tanh^{-1}(e^{-i\theta}a(0))$. The evolution under the squeezing
Hamiltonian $\hat{H}_{\xi}^{S}$ is given by

$\begin{cases}a(t)=e^{i\theta}\tanh(rt+A),\\\
b(t)=\frac{\cosh(rt+A)}{\cosh(A)}b(0),\\\
c(t)=c(0)+\left(n+\frac{1}{2}\right)\log\left(\frac{\cosh(rt+A)}{\cosh(A)}\right)-\frac{e^{i\theta}b^{2}(0)}{4\cosh^{2}(A)}\left[rt+\sinh(rt)\cosh(rt+2A)\right],\end{cases}$
(169)

and the zeros $\lambda_{k}(t)$ are the eigenvalues of the matrix $\Lambda(t)$
defined as

$\Lambda_{kl}(t)=\begin{cases}\lambda_{k}(0)\left[\cosh(rt)+ie^{-i\theta}a(0)\sinh(rt)\right]-ie^{-i\theta}b(0)\sinh(rt)-ie^{-i\theta}\sum_{j\neq
k}\frac{\sinh(rt)}{\lambda_{k}(0)-\lambda_{j}(0)},\quad&k=l\\\
\frac{ie^{-i\theta}\sinh(rt)}{\lambda_{l}(0)-\lambda_{k}(0)},\quad&k\neq
l.\end{cases}$ (170)

###### Proof.

Setting $\alpha=\varphi=0$ in Theorem 1 yields:

$\begin{cases}\frac{da(t)}{dt}=\xi^{*}a^{2}(t)-\xi,\\\
\frac{db(t)}{dt}=\xi^{*}a(t)b(t),\\\
\frac{dc(t)}{dt}=(n+\frac{1}{2})\xi^{*}a(t)-\frac{1}{2}\xi b^{2}(t)\\\
\frac{d^{2}\lambda_{k}(t)}{dt^{2}}=|\xi|^{2}\lambda_{k}(t)-2\xi^{*2}\sum_{j\neq
k}\frac{1}{(\lambda_{k}(t)-\lambda_{j}(t))^{3}},\quad\forall
k\in\\{1,\dots,n\\}.\end{cases}$ (171)

with the initial conditions

$\frac{d\lambda_{k}(t)}{dt}\bigg{|}_{t=0}=-\xi^{*}a(0)\lambda_{k}(0)+\xi^{*}b(0)+\xi^{*}\sum_{j\neq
k}\frac{1}{\lambda_{k}(0)-\lambda_{j}(0)},\quad\forall k\in\\{1,\dots,n\\}.$
(172)

The first differential equation in Eq. (171) is separable and can be solved
directly:

$a(t)=e^{i\theta}\tanh\left[rt+\tanh^{-1}(e^{-i\theta}a(0))\right]=e^{i\theta}\tanh\left(rt+A\right),$
(173)

where we have used $\xi=re^{i\theta}$ and $A=\tanh^{-1}(e^{-i\theta}a(0))$.
Plugging this expression in the second differential equation gives another
separable equation which yields:

$b(t)=\frac{\cosh(rt+A)}{\cosh(A)}b(0).$ (174)

With these expressions for $a(t)$ and $b(t)$, integrating the third
differential equation gives

$c(t)=c(0)+\left(n+\frac{1}{2}\right)\log\left(\frac{\cosh(rt+A)}{\cosh(A)}\right)-\frac{e^{i\theta}b^{2}(0)}{4\cosh^{2}(A)}\left[rt+\sinh(rt)\cosh(rt+2A)\right].$
(175)

Finally, the remaining dynamical system for the zeros is the complex version
of the Calogero–Moser system [52] with reverted harmonic potential, which may
be solved using the Olshanetsky–Perelomov projection method [54]. We obtain
that the zeros $\\{\lambda_{k}\\}_{k=1,\dots,n}$ are given by the eigenvalues
of a matrix $\Lambda(t)$ with analytical time-dependent coefficients defined
as:

$\displaystyle\Lambda_{kl}(t)$
$\displaystyle=\begin{cases}\lambda_{k}(0)\cosh(rt)+\frac{d\lambda_{k}(0)}{dt}\frac{\sinh(rt)}{ir},\quad&k=l\\\
\frac{i\xi^{*}\sinh(rt)}{r(\lambda_{l}(0)-\lambda_{k}(0))},\quad&k\neq
l\end{cases}$ (176)
$\displaystyle=\begin{cases}\lambda_{k}(0)\left[\cosh(rt)+ie^{-i\theta}a(0)\sinh(rt)\right]-ie^{-i\theta}b(0)\sinh(rt)-ie^{-i\theta}\sum_{j\neq
k}\frac{\sinh(rt)}{\lambda_{k}(0)-\lambda_{j}(0)},\quad&k=l\\\
\frac{ie^{-i\theta}\sinh(rt)}{\lambda_{l}(0)-\lambda_{k}(0)},\quad&k\neq
l,\end{cases}$

where $\xi=re^{i\theta}$ and where we used Eq. (172) in the second line.

∎

##  Appendix D Proofs of the properties of the multimode stellar hierarchy

In this section we prove Theorem 3, Corollary 4, and Theorem 5 from Sec. 4,
which we recall below:

######  Theorem 3.

An $m$-mode unitary operation is Gaussian if and only if it leaves the
$m$-mode stellar rank invariant.

###### Proof.

By Theorem 2, the set of $m$-mode states of rank $0$ is equal to the set of
pure Gaussian states. Hence, if a unitary operation leaves the stellar rank
invariant, then it maps Gaussian states onto Gaussian states, i.e., it is a
Gaussian operation.

We now show that Gaussian unitary operations leave the stellar rank invariant.
We first consider the case of pure states. Since any $m$-mode Gaussian unitary
operation $\hat{G}$ can be written as

$\hat{G}=\hat{U}\hat{S}\hat{D}\hat{V},$ (177)

where $\hat{U}$ and $\hat{V}$ are passive linear operations, $\hat{S}$ is a
tensor product of single-mode squeezing operators, and $\hat{D}$ is a tensor
product of single-mode displacement operators, it suffices to prove that these
three types of Gaussian unitaries leave the stellar rank invariant.

We first consider the case of states of finite stellar rank. Let
$n\in\mathbb{N}$ and let $\ket{\bm{\psi}}\in\mathcal{H}^{\otimes m}$ be a
state of stellar rank $n\in\mathbb{N}$.

By Lemma 8, there exists a Gaussian state $\ket{G}$ such that
$\ket{\bm{\psi}}\propto P(\hat{\bm{a}}^{\dagger})\ket{G}$, where $P$ is a
multivariate polynomial of degree $n$.

For all $\bm{\alpha}\in\mathbb{C}^{m}$,

$\displaystyle\hat{D}(\bm{\alpha})\ket{\bm{\psi}}$
$\displaystyle\propto\hat{D}(\bm{\alpha})P(\hat{\bm{a}}^{\dagger})\ket{G}$
(178) $\displaystyle\propto
P(\hat{\bm{a}}^{\dagger}+\bm{\alpha}^{*})\hat{D}(\bm{\alpha})\ket{G},$

where we used Eq. (17). The state $\hat{D}(\bm{\alpha})\ket{G}$ is a Gaussian
state and $\bm{z}\mapsto P(\bm{z}+\bm{\alpha}^{*})$ is a multivariate
polynomial of degree $n$, so by Lemma 8, the state
$\hat{D}(\bm{\alpha})\ket{\bm{\psi}}$ has stellar rank $n$.

Similarly, for any passive linear operation $\hat{U}$, with unitary matrix $U$
we have

$\displaystyle\hat{U}\ket{\bm{\psi}}$
$\displaystyle\propto\hat{U}P(\hat{\bm{a}}^{\dagger})\ket{G}$ (179)
$\displaystyle\propto P(U\hat{\bm{a}}^{\dagger})\hat{U}\ket{G},$

where we used Eq. (24). The state $\hat{U}\ket{G}$ is a Gaussian state and
$\bm{z}\mapsto P(U\bm{z})$ is a multivariate polynomial of degree $n$, since
$U$ is unitary, so by Lemma 8, the state $\hat{U}\ket{\bm{\psi}}$ has stellar
rank $n$.

The case of squeezing is more technical. For all
$\bm{\xi}=(\xi_{1},\dots,\xi_{m})\in\mathbb{C}^{m}$,
$\hat{S}(\bm{\xi})=\bigotimes_{i=1}^{m}\hat{S}(\xi_{i})$, so it is sufficient
prove the result for a single-mode squeezing of an $m$-mode state. Without
loss of generality (up to a passive linear operation), we consider a squeezing
on mode $1$. For all $\xi\in\mathbb{C}$,

$\displaystyle(\hat{S}(\xi)\otimes\mathbb{I}^{\otimes m-1})\ket{\bm{\psi}}$
$\displaystyle\propto\hat{(}S(\xi)\otimes\mathbb{I}^{\otimes
m-1})P(\hat{\bm{a}}^{\dagger})\ket{G}$ (180) $\displaystyle\propto P((\cosh
r)\hat{a}_{1}^{\dagger}-e^{-i\theta}(\sinh
r)\hat{a}_{1},\hat{a}_{2}^{\dagger},\dots,\hat{a}_{m}^{\dagger})(\hat{S}(\xi)\otimes\mathbb{I}^{\otimes
m-1})\ket{G},$

where we used Eq. (17) and where we have set $\xi:=re^{i\theta}$. With the
correspondence $\hat{a}_{k}\leftrightarrow\partial_{z_{k}}$ and
$\hat{a}_{k}^{\dagger}\leftrightarrow\times z_{k}$, the stellar function of
this state reads:

$P((\cosh r)z_{1}-e^{-i\theta}(\sinh
r)\partial_{z_{1}},z_{2},\dots,z_{m})G(\bm{z}),$ (181)

where the Gaussian function $G(\bm{z})$ is the stellar function of the
Gaussian state $(\hat{S}(\xi)\otimes\mathbb{I}^{\otimes m-1})\ket{G}$.
Computing the partial derivatives gives an expression of the form
$\tilde{P}(\bm{z})G(\bm{z})$, and we are left with checking that the degree of
$\tilde{P}$ is indeed equal to $n$.

We write $P(\bm{z}):=\sum_{|\bm{n}|\leq n}p_{\bm{n}}\bm{z}^{\bm{n}}$ and
$G(\bm{z}):=e^{-\frac{1}{2}\bm{z}^{T}A\bm{z}+B^{T}\bm{z}}$ with $|A_{ij}|<1$
for all $1\leq i,j\leq m$ by Eq. (102). The leading monomials of $\tilde{P}$
are given by the leading monomials of:

$\displaystyle
G(\bm{z})^{-1}\left(\sum_{\begin{subarray}{c}\bm{n}=(n_{1},\dots,n_{m})\\\
n_{1}+\dots+n_{m}=n\end{subarray}}p_{\bm{n}}((\cosh r)z_{1}-e^{-i\theta}(\sinh
r)\partial_{z_{1}})^{n_{1}}z_{2}^{n_{2}}\cdots z_{m}^{n_{m}}\right)G(\bm{z})$
(182)
$\displaystyle\quad\quad=G(\bm{z})^{-1}\left(\sum_{\begin{subarray}{c}\bm{n}=(n_{1},\dots,n_{m})\\\
n_{1}+\dots+n_{m}=n\end{subarray}}p_{\bm{n}}((\cosh r)z_{1}-e^{-i\theta}(\sinh
r)\partial_{z_{1}})^{n_{1}}z_{2}^{n_{2}}\cdots
z_{m}^{n_{m}}\right)e^{-\frac{1}{2}\sum_{ij}A_{ij}z_{i}z_{j}+\sum_{k}B_{k}z_{k}}.$

Without loss of generality, we may discard the terms depending on $B$ in the
exponentials, since they only lead to lower order terms. Hence, the leading
monomials of $\tilde{P}$ are given by the leading monomials of:

$e^{\frac{1}{2}\sum_{ij}A_{ij}z_{i}z_{j}}\sum_{\begin{subarray}{c}\bm{n}=(n_{1},\dots,n_{m})\\\
n_{1}+\dots+n_{m}=n\end{subarray}}p_{\bm{n}}\cosh^{n_{1}}r(z_{1}-t\partial_{z_{1}})^{n_{1}}z_{2}^{n_{2}}\cdots
z_{m}^{n_{m}}e^{-\frac{1}{2}\sum_{ij}A_{ij}z_{i}z_{j}},$ (183)

where we have set $t=e^{-i\theta}\tanh r$, with $|t|<1$. Now for all
$t,z\in\mathbb{C}$ and all $n\in\mathbb{N}$ we have

$(z-t\partial_{z})^{n}=\sum_{j=0}^{\lfloor\frac{n}{2}\rfloor}\frac{(-1)^{j}n!t^{j}}{(n-2j)!j!2^{j}}\sum_{l=0}^{n-2j}\binom{n-2j}{l}z^{l}(-t\partial_{z})^{n-2j-l},$
(184)

so the leading monomials of $\tilde{P}$ are given by those of:

$e^{\frac{1}{2}\sum_{ij}A_{ij}z_{i}z_{j}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}\bm{n}=(n_{1},\dots,n_{m})\\\
n_{1}+\dots+n_{m}=n\end{subarray}}\\!\\!\\!\\!\\!\\!\\!\\!p_{\bm{n}}\cosh^{n_{1}}r\sum_{j=0}^{\lfloor\frac{n_{1}}{2}\rfloor}\frac{(-1)^{j}n_{1}!t^{j}}{(n_{1}-2j)!j!2^{j}}\sum_{l=0}^{n_{1}-2j}\binom{n_{1}-2j}{l}z_{1}^{l}z_{2}^{n_{2}}\cdots
z_{m}^{n_{m}}(-t\partial_{z_{1}})^{n_{1}-2j-l}e^{-\frac{1}{2}\sum_{ij}A_{ij}z_{i}z_{j}}.$
(185)

Since the matrix $A$ is symmetric, each partial derivative
$-t\partial_{z_{1}}$ of the exponential yields a term of the form
$(t\sum_{k}{A_{1k}}z_{k})e^{-\frac{1}{2}\sum_{ij}A_{ij}z_{i}z_{j}}$, thus
increasing the degree of the corresponding monomials by $1$ (other terms being
of lower order), so the leading monomials of $\tilde{P}$ are given by those
of:

$\sum_{\begin{subarray}{c}\bm{n}=(n_{1},\dots,n_{m})\\\
n_{1}+\dots+n_{m}=n\end{subarray}}\\!\\!p_{\bm{n}}\cosh^{n_{1}}r\sum_{j=0}^{\lfloor\frac{n_{1}}{2}\rfloor}\frac{(-1)^{j}n_{1}!t^{j}}{(n_{1}-2j)!j!2^{j}}\sum_{l=0}^{n_{1}-2j}\binom{n_{1}-2j}{l}z_{1}^{l}\left(t\sum_{k}{A_{1k}}z_{k}\right)^{n_{1}-2j-l}z_{2}^{n_{2}}\cdots
z_{m}^{n_{m}}.$ (186)

Terms corresponding to $j>0$ in the above equation are of lower order, while
the terms $j=0$ give:

$\displaystyle\sum_{\begin{subarray}{c}\bm{n}=(n_{1},\dots,n_{m})\\\
n_{1}+\dots+n_{m}=n\end{subarray}}\\!\\!p_{\bm{n}}\cosh^{n_{1}}r\sum_{l=0}^{n_{1}}\binom{n_{1}}{l}z_{1}^{l}\left(t\sum_{k}{A_{1k}}z_{k}\right)^{n_{1}-l}z_{2}^{n_{2}}\cdots
z_{m}^{n_{m}}$ (187)
$\displaystyle=\sum_{\begin{subarray}{c}\bm{n}=(n_{1},\dots,n_{m})\\\
n_{1}+\dots+n_{m}=n\end{subarray}}\\!\\!p_{\bm{n}}\cosh^{n_{1}}r\left(z_{1}+t\sum_{k}{A_{1k}}z_{k}\right)^{n_{1}}z_{2}^{n_{2}}\cdots
z_{m}^{n_{m}}$
$\displaystyle=\sum_{\begin{subarray}{c}\bm{n}=(n_{1},\dots,n_{m})\\\
n_{1}+\dots+n_{m}=n\end{subarray}}\\!\\!p_{\bm{n}}\cosh^{n_{1}}r\left[(1+tA_{11})z_{1}+t\sum_{k>1}{A_{1k}}z_{k}\right]^{n_{1}}z_{2}^{n_{2}}\cdots
z_{m}^{n_{m}}.$

The polynomial $\tilde{P}$ thus has degree less or equal to $n$. Moreover, the
polynomial $P$—being of degree $n$—has nonzero monomials of degree $n$. Let
$\bm{n}^{\prime}=(n_{1}^{\prime},\dots,n_{m}^{\prime})$ with
$|\bm{n}^{\prime}|=N$ be one such monomial with a maximal power for $z_{1}$
(note that it is not necessarily unique). We have $p_{\bm{n}^{\prime}}\neq 0$,
and since $n_{1}^{\prime}$ is maximal the coefficient $\bm{n}^{\prime}$ of
$\tilde{P}$ is given by

$p_{\bm{n}^{\prime}}\cosh^{n_{1}^{\prime}}r(1+tA_{11})^{n_{1}^{\prime}},$
(188)

which is nonzero since $|t|<1$ and $|A_{11}|<1$. This shows that $\tilde{P}$
has degree equal to $n$. Hence, all Gaussian unitary operations leave the
finite stellar rank of pure states invariant.

Now assume by contradiction that a Gaussian unitary operation $\hat{G}$ maps a
state $\ket{\psi}$ of infinite stellar rank onto a state of finite stellar
rank $\ket{\phi}$. Then, $\hat{G}^{\dagger}$ would map the state $\ket{\phi}$
of finite stellar rank onto $\ket{\psi}$ of infinite stellar rank, which
contradicts the previous proof that Gaussian unitary operations leave the
finite stellar rank invariant. Gaussian unitary operations thus map pure
states of infinite stellar rank onto pure states of infinite stellar rank.

For a mixed state $\rho$, the stellar rank is given by Eq. (87). For any
Gaussian unitary operation $\hat{G}$,

$\displaystyle r^{\star}(\hat{G}\rho\hat{G}^{\dagger})$
$\displaystyle=\inf_{\\{p_{i},\psi_{i}\\}}\sup_{i}r^{\star}(\hat{G}\ket{\psi_{i}})$
(189) $\displaystyle=\inf_{\\{p_{i},\psi_{i}\\}}\sup_{i}r^{\star}(\psi_{i})$
$\displaystyle=r^{\star}(\rho),$

where we used the fact that Gaussian unitary operations leave the stellar rank
of pure states invariant, and which concludes the proof.

∎

An interesting consequence is the following result:

######  Corollary 4.

The stellar rank is non-increasing under Gaussian channels and measurements.

###### Proof.

The result is trivial if the state has infinite rank. We consider a finite
rank state $\rho$ over $m$ modes in what follows.

Any Gaussian measurement can be implemented by a Gaussian unitary followed by
heterodyne detection [106, 75], with (single-mode) positive operator-valued
measure
$\\{\frac{1}{\pi}\ket{\alpha}\\!\bra{\alpha}\\}_{\alpha\in\mathbb{C}}$, so
with Theorem 3 it suffices to show that measuring a subset of the modes of
$\rho$ with heterodyne detection does not increase its stellar rank. We pick a
decomposition $\rho=\sum_{i}p_{i}\ket{\psi_{i}}\\!\bra{\psi_{i}}$ with
$r^{\star}(\rho)=\sup_{i}r^{\star}(\psi_{i})$. The unnormalised state obtained
from $\rho$ by measuring, say, its first $k$ modes with heterodyne detection
and obtaining the outcomes
$\bm{\alpha}=(\alpha_{1},\dots,\alpha_{k})\in\mathbb{C}^{k}$ is given by:

$\text{Tr}_{1\dots
k}(\Pi_{\bm{\alpha}}\rho)=\sum_{i}p_{i}\ket{\phi^{\bm{\alpha}}_{i}}\bra{\phi^{\bm{\alpha}}_{i}},$
(190)

where
$\Pi_{\bm{\alpha}}=\frac{1}{\pi^{k}}\ket{\alpha_{1}\dots\alpha_{k}}\\!\bra{\alpha_{1}\dots\alpha_{k}}\otimes\mathbb{I}_{m-k}$,
and where $\ket{\phi^{\bm{\alpha}}_{i}}$ are the unnormalised pure states
obtained by projecting the first $k$ modes of $\ket{\psi_{i}}$ onto the
$k$-mode coherent state $\ket{\alpha_{1}\dots\alpha_{k}}$. For all $i$, the
Husimi $Q$ function of the state $\ket{\phi^{\bm{\alpha}}_{i}}$ is given by

$\displaystyle Q_{\phi^{\bm{\alpha}}_{i}}(\beta_{1},\dots,\beta_{m-k})$
$\displaystyle=\frac{1}{\pi^{m-k}}\operatorname*{Tr}(\ket{\phi^{\bm{\alpha}}_{i}}\bra{\phi^{\bm{\alpha}}_{i}}\ket{\beta_{1}\dots\beta_{m-k}}\bra{\beta_{1}\dots\beta_{m-k}})$
(191)
$\displaystyle=\frac{1}{\pi^{m-k}}\operatorname*{Tr}\left(\frac{1}{\pi^{k}}\text{Tr}_{1\dots
k}(\ket{\psi_{i}}\bra{\psi_{i}}\ket{\alpha_{1}\dots\alpha_{k}}\bra{\alpha_{1}\dots\alpha_{k}})\ket{\beta_{1}\dots\beta_{m-k}}\bra{\beta_{1}\dots\beta_{m-k}})\right)$
$\displaystyle=\frac{1}{\pi^{m}}\operatorname*{Tr}\left(\ket{\psi_{i}}\bra{\psi_{i}}\ket{\alpha_{1}\dots\alpha_{k}\beta_{1}\dots\beta_{m-k}}\bra{\alpha_{1}\dots\alpha_{k}\beta_{1}\dots\beta_{m-k}}\right)$
$\displaystyle=Q_{\psi_{i}}(\alpha_{1},\dots,\alpha_{k},\beta_{1},\dots,\beta_{m-k}).$

Recalling the relation between the stellar function and the Husimi $Q$
function in Eq. (83), we obtain that the stellar function of the state
$\ket{\phi^{\bm{\alpha}}_{i}}$ is given by:

$F^{\star}_{\phi^{\bm{\alpha}}_{i}}(z_{1},\dots,z_{m-k})\propto
F^{\star}_{\psi_{i}}(\alpha_{1}^{*},\dots,\alpha_{k}^{*},z_{1},\dots,z_{m-k}),$
(192)

up to a constant prefactor. Since $r^{\star}(\psi_{i})\leq r^{\star}(\rho)$,
the stellar function of the state $\ket{\psi_{i}}$ is of the form $P\times G$,
where $P$ is a polynomial of degree less or equal to $r^{\star}(\rho)$, and
$G$ is a Gaussian function. Hence, the stellar function of the state
$\ket{\phi^{\bm{\alpha}}_{i}}$ is also of the form $P\times G$, where $P$ is a
polynomial of degree less or equal to $r^{\star}(\rho)$, and $G$ is a Gaussian
function. In particular, $r^{\star}(\phi^{\bm{\alpha}}_{i})\leq
r^{\star}(\rho)$. With Eq. (190) we obtain
$r^{\star}(\operatorname*{Tr}_{1\dots k}(\Pi_{\bm{\alpha}}\rho))\leq
r^{\star}(\rho)$ and therefore the stellar rank is non-increasing under
Gaussian measurements.

With Eq. (89) this directly implies that the stellar rank is non-increasing
under partial trace, since we have

$\displaystyle\text{Tr}_{1\dots k}(\rho)$
$\displaystyle=\int_{\bm{\alpha}\in\mathbb{C}^{k}}\text{Tr}_{1\dots
k}(\Pi_{\bm{\alpha}}\rho)d^{2k}\bm{\alpha}$ (193)
$\displaystyle=\int_{\bm{\alpha}\in\mathbb{C}^{k}}p_{\bm{\alpha}}\rho_{\bm{\alpha}}d^{2k}\bm{\alpha},$

where
$\Pi_{\bm{\alpha}}=\frac{1}{\pi^{k}}\ket{\alpha_{1}\dots\alpha_{k}}\\!\bra{\alpha_{1}\dots\alpha_{k}}\otimes\mathbb{I}_{m-k}$,
$p_{\bm{\alpha}}:=\operatorname*{Tr}(\Pi_{\bm{\alpha}}\rho)$ and
$\rho_{\bm{\alpha}}:=\frac{\text{Tr}_{1\dots
k}(\Pi_{\bm{\alpha}}\rho)}{\operatorname*{Tr}(\Pi_{\bm{\alpha}}\rho)}$.

Finally, any Gaussian channel may be implemented by taking the tensor product
with a mixture of Gaussian states, applying a Gaussian unitary operation,
applying a Gaussian measurement and taking partial trace. By Eq. (88), taking
the tensor product with a mixture of Gaussian states leave the stellar rank
invariant. By Theorem 3, applying a Gaussian unitary operation also leave the
stellar rank invariant. By the above, applying a Gaussian measurement and
taking partial trace do not increase the stellar rank. Hence, the stellar rank
is non-increasing under Gaussian channels.

∎

We now turn to the proof of Theorem 5:

######  Theorem 5.

The $m$-mode stellar hierarchy is robust for the trace norm. Formally, for all
$n\in\mathbb{N}$,

$\overline{\bm{R}_{n}}=\bigcup_{0\leq k\leq n}\bm{R}_{k},$ (194)

where $\overline{X}$ denotes the closure of $X$ for the trace norm.

###### Proof.

The proof is similar to that of [42] in the single-mode case. We outline the
main arguments below.

We first show that $\bigcup_{0\leq k\leq
n}\bm{R}_{k}\subset\overline{\bm{R}_{n}}$. We have
$\bm{R}_{n}\subset\overline{\bm{R}_{n}}$. Let
$\ket{\bm{\psi}}\in\bigcup_{0\leq k<n}\bm{R}_{k}$. By Lemma 10, there exists
an $m$-mode Gaussian unitary operation $\hat{G}_{\bm{\psi}}$ and a core state
$\ket{C_{\bm{\psi}}}$ of stellar rank smaller than $n$ such that
$\ket{\bm{\psi}}=\hat{G}_{\bm{\psi}}\ket{C_{\bm{\psi}}}$. Let
$\bm{n}=(n,0,\dots,0)\in\mathbb{N}^{m}$. For all $l\geq 1$, let

$\displaystyle\ket{\bm{\psi}_{l}}$
$\displaystyle:=\sqrt{\frac{1}{m}}\ket{\bm{\psi}}+\sqrt{1-\frac{1}{m}}\hat{G}_{\bm{\psi}}\ket{\bm{n}}$
(195)
$\displaystyle=\hat{G}_{\bm{\psi}}\left(\sqrt{\frac{1}{m}}\ket{C_{\bm{\psi}}}+\sqrt{1-\frac{1}{m}}\ket{\bm{n}}\right).$

By Lemma 10, $\ket{\bm{\psi}_{l}}\in\bm{R}_{n}$ for all $l\geq 1$, and by
construction the sequence $(\ket{\bm{\psi}_{l}})_{l\geq 1}$ converges to
$\ket{\bm{\psi}}$ in trace norm. This shows that $\bigcup_{0\leq k\leq
n}\bm{R}_{k}\subset\overline{\bm{R}_{n}}$.

For the reverse inclusion, let $(\ket{\bm{\psi}_{l}})_{l\geq
1}\in\bm{R}_{n}^{\mathbb{N}}$ be a converging sequence and let
$\ket{\bm{\psi}}\in\mathcal{H}^{\otimes m}$ denote its limit. By Lemma 10,
there exist a sequence of passive linear unitaries $(\hat{U}_{l})_{l\geq 1}$,
a sequence of squeezing amplitudes $(\bm{\xi}_{l})_{l\geq 1}$, a sequence of
displacement amplitudes $(\bm{\beta}_{l})_{l\geq 1}$ and a sequence of core
states $(\ket{C_{\bm{\psi}_{l}}})_{l\geq 1}$ of stellar rank $n$ such that for
all $l\geq 1$,
$\ket{\bm{\psi}_{l}}=\hat{U}_{l}\hat{S}(\bm{\xi}_{l})\hat{D}(\bm{\beta}_{l})\ket{C_{\bm{\psi}_{l}}}$.

The set of normalised core states of stellar rank less of equal to $n$ is
compact (it identifies with the set of normalised multivariate polynomials of
degree less or equal to $n$), so up to extracting a converging subsequence we
may assume that the sequence $(\ket{C_{\bm{\psi}_{l}}})_{l\geq 1}$ converges
in trace norm to some core state $\ket{C}:=\sum_{|\bm{n}|\leq
n}c_{\bm{n}}\ket{\bm{n}}$ of stellar rank less or equal to $n$. Similarly, the
set of passive linear unitaries is compact (it identifies with the set of
$m\times m$ unitary matrices), so up to extracting a converging subsequence we
may assume that the sequence $(\hat{U}_{l})_{l\geq 1}$ converges in strong
operator topology to some passive linear unitary operator $\hat{U}$. Since the
trace distance is invariant under unitary transformations, we obtain that
$\hat{S}(\bm{\xi}_{l})\hat{D}(\bm{\beta}_{l})\ket{C}$ converges to
$\hat{U}^{\dagger}\ket{\bm{\psi}}$ in trace norm.

In order to conclude, we show hereafter that the sequences
$(\bm{\xi}_{l})_{l\geq 1}$ and $(\bm{\beta}_{l})_{l\geq 1}$ are bounded. When
that is the case, we may extract converging subsequences with limits
$\bm{\xi}$ and $\bm{\beta}$, respectively, thus showing that
$\hat{S}(\bm{\xi}_{l})\hat{D}(\bm{\beta}_{l})\ket{C}$ converges along that
subsequence to $\hat{S}(\bm{\xi})\hat{D}(\bm{\beta})\ket{C}$. This implies in
turn that $\ket{\bm{\psi}}=\hat{U}\hat{S}(\bm{\xi})\hat{D}(\bm{\beta})\ket{C}$
and hence $\ket{\bm{\psi}}\in\bigcup_{0\leq k\leq n}\bm{R}_{k}$ by Lemma 10.

We make use of the following expression from Eq. (80) of the stellar function
of a single-mode state of the form $\hat{S}(\xi)\ket{n}$, where
$\xi\in\mathbb{C}$ and $\ket{n}$ is a Fock state:

$F^{\star}_{\hat{S}(\xi)\ket{n}}(z)=\frac{(e^{-i\theta}\tanh
r)^{n/2}}{\sqrt{n!\cosh r}}\mathit{He}_{n}\left(\frac{z}{\cosh
r\sqrt{e^{-i\theta}\tanh r}}\right)e^{\frac{1}{2}e^{i\theta}(\tanh r)z^{2}},$
(196)

for all $z\in\mathbb{C}$, where $\xi:=re^{i\theta}$ and where
$\mathit{He}_{n}(x)=\sum_{k=0}^{\lfloor\frac{n}{2}\rfloor}\frac{(-1)^{k}n!}{2^{k}k!(n-2k)!}x^{n-2k}$
is the $n^{th}$ Hermite polynomial. In particular, the $Q$ function of the
state $\hat{S}(\bm{\xi}_{l})\hat{D}(\bm{\beta}_{l})\ket{C}$ is given by

$\displaystyle
Q_{\hat{S}(\bm{\xi}_{l})\hat{D}(\bm{\beta}_{l})\ket{C}}(\bm{\alpha})$
$\displaystyle=Q_{\hat{D}(\bm{\gamma}_{l})\hat{S}(\bm{\xi}_{l})\ket{C}}(\bm{\alpha})$
(197)
$\displaystyle=Q_{\hat{S}(\bm{\xi}_{l})\ket{C}}(\bm{\alpha}-\bm{\gamma}_{l})$
$\displaystyle=\frac{1}{\pi^{m}}e^{-\|\bm{\alpha}-\bm{\gamma}_{l}\|^{2}}\left|F^{\star}_{\hat{S}(\bm{\xi}_{l})\ket{C}}(\bm{\alpha}^{*}-\bm{\gamma}_{l}^{*})\right|^{2}$
$\displaystyle=\frac{1}{\pi^{m}}e^{-\|\bm{\alpha}-\bm{\gamma}_{l}\|^{2}}\left|\sum_{\begin{subarray}{c}\bm{n}=(n_{1},\dots,n_{m})\\\
|\bm{n}|\leq
n\end{subarray}}\\!\\!\\!\\!\\!c_{\bm{n}}\prod_{j=1}^{m}F^{\star}_{\hat{S}(\xi_{j,l})\ket{n_{j}}}(\alpha_{j}^{*}-\gamma_{j,l}^{*})\right|^{2}$
$\displaystyle=\frac{1}{\pi^{m}}e^{-\|\bm{\alpha}-\bm{\gamma}_{l}\|^{2}}\Bigg{|}\sum_{\begin{subarray}{c}\bm{n}=(n_{1},\dots,n_{m})\\\
|\bm{n}|\leq
n\end{subarray}}\\!\\!\\!\\!\\!c_{\bm{n}}\prod_{j=1}^{m}\frac{(e^{-i\theta_{j,l}}\tanh
r_{j,l})^{n_{j}/2}}{\sqrt{n_{j}!\cosh r_{j,l}}}$
$\displaystyle\quad\quad\quad\times\mathit{He}_{n_{j}}\\!\\!\left(\frac{\alpha_{j}^{*}-\gamma_{j,l}^{*}}{\cosh
r_{j,l}\sqrt{e^{-i\theta_{j,l}}\tanh
r_{j,l}}}\right)e^{\frac{1}{2}e^{i\theta_{j,l}}(\tanh
r_{j,l})(\alpha_{j}^{*}-\gamma_{j,l}^{*})^{2}}\Bigg{|}^{2},$

for all $\bm{\alpha}\in\mathbb{C}^{m}$, where we have set
$\bm{\xi}_{l}=(\xi_{1,l},\dots,\xi_{m,l})\in\mathbb{C}^{m}$ with
$\xi_{k,l}:=r_{k,l}e^{i\theta_{k,l}}$,
$\bm{\beta}_{l}=(\beta_{1,l},\dots,\beta_{m,l})\in\mathbb{C}^{m}$, and
$\bm{\gamma}_{l}=(\gamma_{1,l},\dots,\gamma_{1,l})$ with
$\gamma_{k,l}=\beta_{k,l}\cosh r_{k,l}+\beta_{k,l}^{*}e^{i\theta_{k,l}}\sinh
r_{k,l}$. The core state $\ket{C}$ is normalised, so using Cauchy–Schwarz
inequality, we obtain the bound:

$\displaystyle
Q_{\hat{S}(\bm{\xi}_{l})\hat{D}(\bm{\beta}_{l})\ket{C}}(\bm{\alpha})$
$\displaystyle\leq\frac{1}{\pi^{m}}e^{-\|\bm{\alpha}-\bm{\gamma}_{l}\|^{2}}\\!\\!\\!\\!\sum_{\begin{subarray}{c}\bm{n}=(n_{1},\dots,n_{m})\\\
|\bm{n}|\leq
n\end{subarray}}\prod_{j=1}^{m}\Bigg{|}\frac{(e^{-i\theta_{j,l}}\tanh
r_{j,l})^{n_{j}/2}}{\sqrt{n_{j}!\cosh r_{j,l}}}$ (198)
$\displaystyle\quad\quad\quad\times\mathit{He}_{n_{j}}\\!\\!\left(\frac{\alpha_{j}^{*}-\gamma_{j,l}^{*}}{\cosh
r_{j,l}\sqrt{e^{-i\theta_{j,l}}\tanh
r_{j,l}}}\right)e^{\frac{1}{2}e^{i\theta_{j,l}}(\tanh
r_{j,l})(\alpha_{j}^{*}-\gamma_{j,l}^{*})^{2}}\Bigg{|}^{2}$
$\displaystyle=\frac{e^{\sum_{j=1}^{m}-|\alpha_{j}-\gamma_{j,l}|^{2}+\frac{1}{2}\tanh
r_{j,l}\left(e^{-i\theta_{j,l}}(\alpha_{j}-\gamma_{j,l})^{2}+e^{i\theta_{j,l}}(\alpha^{*}_{j}-\gamma^{*}_{j,l})^{2}\right)}}{\pi^{m}\prod_{j=1}^{m}\cosh
r_{j,l}}$
$\displaystyle\quad\quad\quad\times\sum_{\begin{subarray}{c}\bm{n}=(n_{1},\dots,n_{m})\\\
|\bm{n}|\leq
n\end{subarray}}\prod_{j=1}^{m}\left|\sum_{k=0}^{\lfloor\frac{n_{j}}{2}\rfloor}\frac{(-1)^{k}\sqrt{n_{j}!}}{2^{k}k!(n_{j}-2k)!}\left(\frac{\alpha_{j}^{*}-\gamma_{j,l}^{*}}{e^{-i\theta_{j,l}/2}\cosh
r_{j,l}}\right)^{n_{j}-2k}\right|^{2}$
$\displaystyle\leq\frac{e^{\sum_{j=1}^{m}-|\alpha_{j}-\gamma_{j,l}|^{2}+\frac{1}{2}\tanh
r_{j,l}\left(e^{-i\theta_{j,l}}(\alpha_{j}-\gamma_{j,l})^{2}+e^{i\theta_{j,l}}(\alpha^{*}_{j}-\gamma^{*}_{j,l})^{2}\right)}}{\pi^{m}\prod_{j=1}^{m}\cosh
r_{j,l}}$
$\displaystyle\quad\quad\quad\times\sum_{\begin{subarray}{c}\bm{n}=(n_{1},\dots,n_{m})\\\
|\bm{n}|\leq
n\end{subarray}}\prod_{j=1}^{m}\left(\sum_{k=0}^{\lfloor\frac{n_{j}}{2}\rfloor}\frac{\sqrt{n_{j}!}}{2^{k}k!(n_{j}-2k)!}\left(\frac{|\alpha_{j}-\gamma_{j,l}|}{\cosh
r_{j,l}}\right)^{n_{j}-2k}\right)^{2}.$

Setting, for all $l\geq 1$ and all $j\in\\{1,\dots,m\\}$,

$\displaystyle z_{j,l}(\bm{\alpha})$
$\displaystyle:=\frac{\alpha_{j}-\gamma_{j,l}}{e^{i\theta_{j,l}}\cosh
r_{j,l}}$ (199) $\displaystyle=\frac{\alpha_{j}-\beta_{j,l}\cosh
r_{j,l}-\beta_{j,l}^{*}e^{i\theta_{j,l}}\sinh r_{j,l}}{e^{i\theta_{j,l}}\cosh
r_{j,l}},$

we obtain

$\displaystyle
Q_{\hat{S}(\bm{\xi}_{l})\hat{D}(\bm{\beta}_{l})\ket{C}}(\bm{\alpha})$
$\displaystyle\leq\frac{e^{\sum_{j=1}^{m}-\cosh^{2}r_{j,l}|z_{j,l}(\bm{\alpha})|^{2}+\frac{1}{2}\cosh
r_{j,l}\sinh
r_{j,l}\left(z_{j,l}(\bm{\alpha})^{2}+z_{j,l}(\bm{\alpha})^{*2}\right)}}{\pi^{m}\prod_{j=1}^{m}\cosh
r_{j,l}}$ (200)
$\displaystyle\quad\quad\quad\times\sum_{\begin{subarray}{c}\bm{n}=(n_{1},\dots,n_{m})\\\
|\bm{n}|\leq
n\end{subarray}}\prod_{j=1}^{m}\left(\sum_{k=0}^{\lfloor\frac{n_{j}}{2}\rfloor}\frac{\sqrt{n_{j}!}}{2^{k}k!(n_{j}-2k)!}|z_{j,l}(\bm{\alpha})|^{n_{j}-2k}\right)^{2}$
$\displaystyle\leq\frac{e^{-\frac{1}{2}\sum_{j=1}^{m}|z_{j,l}(\bm{\alpha})|^{2}}}{\pi^{m}\prod_{j=1}^{m}\cosh
r_{j,l}}\sum_{\begin{subarray}{c}\bm{n}=(n_{1},\dots,n_{m})\\\ |\bm{n}|\leq
n\end{subarray}}\prod_{j=1}^{m}\left(\sum_{k=0}^{\lfloor\frac{n_{j}}{2}\rfloor}\frac{\sqrt{n_{j}!}}{2^{k}k!(n_{j}-2k)!}|z_{j,l}(\bm{\alpha})|^{n_{j}-2k}\right)^{2}$
$\displaystyle:=\frac{e^{-\frac{1}{2}\sum_{j=1}^{m}|z_{j,l}(\bm{\alpha})|^{2}}}{\pi^{m}\prod_{j=1}^{m}\cosh
r_{j,l}}T(|z_{1,l}(\bm{\alpha})|,\dots,|z_{m,l}(\bm{\alpha})|),$

where we used $(\cosh r-\sinh r)\cosh r>\frac{1}{2}$ and $(\cosh r+\sinh
r)\cosh r>\frac{1}{2}$ in the second line, and where we have defined the
multivariate polynomial $T$. Note that the function
$(z_{1},\dots,z_{m})\mapsto
e^{-\frac{1}{2}\sum_{j=1}^{m}|z_{j}|^{2}}T(|z_{1}|,\dots,|z_{m}|)$ is bounded.

The sequence $(\hat{S}(\bm{\xi}_{l})\hat{D}(\bm{\beta}_{l})\ket{C})_{l\geq 1}$
converges to a normalised state, so by property of the trace norm, the
sequence of $Q$ functions (which are overlaps with coherent states) converges
point-wise to the $Q$ function of this normalised state. In particular, if the
sequence $(\bm{r}_{l})_{l\geq 1}$ is unbounded, then we have
$Q_{\hat{S}(\bm{\xi}_{l})\hat{D}(\bm{\beta}_{l})\ket{C}}(\bm{\alpha})\rightarrow
0$ for all $\bm{\alpha}\in\mathbb{C}^{m}$, so the sequence
$(\bm{r}_{l})_{l\geq 1}$ and thus $(\bm{\xi}_{l})_{l\geq 1}$ are bounded for
the limit state to be normalisable. Similarly, if the sequence
$(\bm{\gamma}_{l})_{l\geq 1}$ is unbounded while $(\bm{\xi}_{l})_{l\geq 1}$ is
bounded, then at least one of the sequences $(z_{j,l})_{l\geq 1}$ from Eq.
(199) is unbounded and we have
$Q_{\hat{S}(\bm{\xi}_{l})\hat{D}(\bm{\beta}_{l})\ket{C}}(\bm{\alpha})\rightarrow
0$ for all $\bm{\alpha}\in\mathbb{C}^{m}$. This shows that all sequences
$(z_{j,l})_{l\geq 1}$, $(\bm{\gamma}_{l})_{l\geq 1}$, and thus
$(\bm{\beta}_{l})_{l\geq 1}$ are bounded. Hence, both sequences
$(\bm{\xi}_{l})_{l\geq 1}$ and $(\bm{\beta}_{l})_{l\geq 1}$ are bounded, which
completes the proof.

∎

## References

  * [1] P. W. Shor, “Algorithms for quantum computation: Discrete logarithms and factoring,” in Proceedings 35th annual symposium on foundations of computer science, pp. 124–134, IEEE.  1994\.
  * [2] A. W. Harrow and A. Montanaro, “Quantum computational supremacy,” Nature 549, 203–209 (2017).
  * [3] S. Wiesner, “Conjugate coding,” ACM Sigact News 15, 78–88 (1983).
  * [4] C. H. Bennett, G. Brassard, C. Crépeau, R. Jozsa, A. Peres, and W. K. Wootters, “Teleporting an unknown quantum state via dual classical and Einstein–Podolsky–Rosen channels,” Phys. Rev. Lett. 70, 1895 (1993).
  * [5] V. Giovannetti, S. Lloyd, and L. Maccone, “Quantum-enhanced measurements: beating the standard quantum limit,” Science 306, 1330–1336 (2004).
  * [6] M. A. Nielsen and I. L. Chuang, “Quantum Computation and Quantum Information: 10th Anniversary Edition,”.  Cambridge University Press, New York, NY, USA, 10th ed., 2011.
  * [7] S. L. Braunstein and P. van Loock, “Quantum information with continuous variables,” Rev. Mod. Phys. 77, 513–577 (2005).
  * [8] G. Adesso, S. Ragy, and A. R. Lee, “Continuous variable quantum information: Gaussian states and beyond,” Open Systems & Information Dynamics 21, 1440001 (2014).
  * [9] D. Gottesman, A. Kitaev, and J. Preskill, “Encoding a qubit in an oscillator,” Phys. Rev. A 64, 012310 (2001).
  * [10] W. Cai, Y. Ma, W. Wang, C.-L. Zou, and L. Sun, “Bosonic quantum error correction codes in superconducting quantum circuits,” Fundamental Research 1, 50–67 (2021).
  * [11] H. Vahlbruch, M. Mehmet, K. Danzmann, and R. Schnabel, “Detection of 15 dB squeezed states of light and their application for the absolute calibration of photoelectric quantum efficiency,” Phys. Rev. Lett. 117, 110801 (2016).
  * [12] S. Yokoyama, R. Ukai, S. C. Armstrong, C. Sornphiphatphong, T. Kaji, S. Suzuki, J.-i. Yoshikawa, H. Yonezawa, N. C. Menicucci, and A. Furusawa, “Ultra-large-scale continuous-variable cluster states multiplexed in the time domain,” Nature Photonics 7, 982 (2013).
  * [13] S. Aaronson and A. Arkhipov, “The computational Complexity of Linear Optics,” Theory of Computing 9, 143 (2013).
  * [14] A. P. Lund, A. Laing, S. Rahimi-Keshari, T. Rudolph, J. L. O’Brien, and T. C. Ralph, “Boson Sampling from a Gaussian State,” Phys. Rev. Lett. 113, 100502 (2014).
  * [15] C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, “Gaussian Boson Sampling,” Phys. Rev. Lett. 119, 170501 (2017).
  * [16] T. Douce, D. Markham, E. Kashefi, E. Diamanti, T. Coudreau, P. Milman, P. van Loock, and G. Ferrini, “Continuous-Variable Instantaneous Quantum Computing is hard to sample,” Phys. Rev. Lett. 118, 070503 (2017).
  * [17] A. P. Lund, S. Rahimi-Keshari, and T. C. Ralph, “Exact Boson Sampling using Gaussian continuous variable measurements,” Phys. Rev. A 96, 022301 (2017).
  * [18] L. Chakhmakhchyan and N. J. Cerf, “Boson sampling with Gaussian measurements,” Phys. Rev. A 96, 032326 (2017).
  * [19] U. Chabaud, T. Douce, D. Markham, P. Van Loock, E. Kashefi, and G. Ferrini, “Continuous-variable sampling from photon-added or photon-subtracted squeezed states,” Physical Review A 96, 062307 (2017).
  * [20] H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C. Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, et al., “Quantum computational advantage using photons,” Science 370, 1460–1463 (2020).
  * [21] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, et al., “Quantum supremacy using a programmable superconducting processor,” Nature 574, 505–510 (2019).
  * [22] J. K. Horner and J. F. Symons, “What Have Google’s Random Quantum Circuit Simulation Experiments Demonstrated about Quantum Supremacy?,” in Advances in Software Engineering, Education, and e-Learning, pp. 411–419.  Springer, 2021.
  * [23] F. Pan and P. Zhang, “Simulating the Sycamore quantum supremacy circuits,” arXiv:2103.03074.
  * [24] B. Villalonga, M. Y. Niu, L. Li, H. Neven, J. C. Platt, V. N. Smelyanskiy, and S. Boixo, “Efficient approximation of experimental Gaussian boson sampling,” arXiv:2109.11525.
  * [25] F. Pan, K. Chen, and P. Zhang, “Solving the sampling problem of the Sycamore quantum supremacy circuits,” arXiv:2111.03011.
  * [26] Y. Tong, V. V. Albert, J. R. McClean, J. Preskill, and Y. Su, “Provably accurate simulation of gauge theories and bosonic systems,” arXiv:2110.06942.
  * [27] S. Lloyd and S. L. Braunstein, “Quantum computation over continuous variables,” in Quantum Information with Continuous Variables, pp. 9–17.  Springer, 1999.
  * [28] S. D. Bartlett, B. C. Sanders, S. L. Braunstein, and K. Nemoto, “Efficient Classical Simulation of Continuous Variable Quantum Information Processes,” Phys. Rev. Lett. 88, 097904 (2002).
  * [29] R. Takagi and Q. Zhuang, “Convex resource theory of non-Gaussianity,” Phys. Rev. A 97, 062337 (2018).
  * [30] Q. Zhuang, P. W. Shor, and J. H. Shapiro, “Resource theory of non-Gaussian operations,” Phys. Rev. A 97, 052317 (2018).
  * [31] F. Albarelli, M. G. Genoni, M. G. Paris, and A. Ferraro, “Resource theory of quantum non-Gaussianity and Wigner negativity,” Physical Review A 98, 052350 (2018).
  * [32] J. Eisert, S. Scheel, and M. B. Plenio, “Distilling Gaussian states with Gaussian operations is impossible,” Phys. Rev. Lett. 89, 137903 (2002).
  * [33] J. Fiurášek, “Gaussian Transformations and Distillation of Entangled Gaussian States,” Phys. Rev. Lett. 89, 137904 (2002).
  * [34] J. Niset, J. Fiurášek, and N. J. Cerf, “No-Go Theorem for Gaussian Quantum Error Correction,” Phys. Rev. Lett. 102, 120501 (2009).
  * [35] J. B. Conway, “Functions of one complex variable II,”, vol. 159.  Springer Science & Business Media, 2012.
  * [36] F. Calogero, “Solution of the one-dimensional N-body problems with quadratic and/or inversely quadratic pair potentials,” Journal of Mathematical Physics 12, 419–436 (1971).
  * [37] J. Radcliffe, “Some properties of coherent spin states,” Journal of Physics A: General Physics 4, 313 (1971).
  * [38] F. T. Arecchi, E. Courtens, R. Gilmore, and H. Thomas, “Atomic Coherent States in Quantum Optics,” Phys. Rev. A 6, 2211–2237 (1972).
  * [39] D. Gross, “Hudson’s theorem for finite-dimensional quantum systems,” Journal of mathematical physics 47, 122107 (2006).
  * [40] E. Majorana, “Atomi orientati in campo magnetico variabile,” Il Nuovo Cimento (1924-1942) 9, 43–50 (1932).
  * [41] A. W. Harrow, “The church of the symmetric subspace,” arXiv:1308.6595.
  * [42] U. Chabaud, D. Markham, and F. Grosshans, “Stellar representation of non-Gaussian quantum states,” Physical Review Letters 124, 063605 (2020).
  * [43] I. E. Segal and G. W. Mackey, “Mathematical problems of relativistic physics,”, vol. 2.  American Mathematical Soc., 1963.
  * [44] V. Bargmann, “On a Hilbert space of analytic functions and an associated integral transform part I,” Communications on pure and applied mathematics 14, 187–214 (1961).
  * [45] K. Husimi, “Some formal properties of the density matrix,” Proceedings of the Physico-Mathematical Society of Japan. 3rd Series 22, 264–314 (1940).
  * [46] F. Ricci, “A contraction of SU (2) to the Heisenberg group,” Monatshefte für Mathematik 101, 211–225 (1986).
  * [47] A. Vourdas, “Analytic representations in quantum mechanics,” Journal of Physics A: Mathematical and General 39, R65 (2006).
  * [48] A. Ferraro, S. Olivares, and M. G. Paris, “Gaussian states in continuous variable quantum information,” arXiv:quant-ph/0503237.
  * [49] L. G. Valiant, “The complexity of computing the permanent,” Theoretical computer science 8, 189–201 (1979).
  * [50] E. Schrödinger, “An undulatory theory of the mechanics of atoms and molecules,” Physical review 28, 1049 (1926).
  * [51] P. Leboeuf, “Phase space approach to quantum dynamics,” Journal of Physics A: Mathematical and General 24, 4575 (1991).
  * [52] F. Calogero, “Exactly solvable two-dimensional many-body problems,” Lettere al Nuovo Cimento (1971-1985) 16, 35–38 (1976).
  * [53] J. Moser, “Three integrable Hamiltonian systems connected with isospectral deformations,” in Surveys in applied mathematics, pp. 235–258.  Elsevier, 1976.
  * [54] M. A. Olshanetsky and A. M. Perelomov, “Geodesic flows on symmetric spaces of zero curvature and explicit solution of the generalized Calogero model for the classical case,” Functional Analysis and Its Applications 10, 237–239 (1976).
  * [55] M. A. Olshanetsky and A. M. Perelomov, “Classical integrable finite-dimensional systems related to Lie algebras,” Physics Reports 71, 313–400 (1981).
  * [56] F. Calogero. http://www.scholarpedia.org/article/Calogero-Moser˙system, 2008\.
  * [57] A. P. Polychronakos, “The physics and mathematics of Calogero particles,” Journal of Physics A: Mathematical and General 39, 12793 (2006).
  * [58] M. Walschaers, “Non-Gaussian Quantum States and Where to Find Them,” PRX Quantum 2, 030204 (2021).
  * [59] M. Amy, D. Maslov, M. Mosca, and M. Roetteler, “A meet-in-the-middle algorithm for fast synthesis of depth-optimal quantum circuits,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 32, 818–830 (2013).
  * [60] M. Beverland, E. Campbell, M. Howard, and V. Kliuchnikov, “Lower bounds on the non-Clifford resources for quantum computations,” Quantum Science and Technology 5, 035009 (2020).
  * [61] A. Lvovsky, P. Grangier, A. Ourjoumtsev, V. Parigi, M. Sasaki, and R. Tualle-Brouri, “Production and applications of non-Gaussian quantum states of light,” arXiv:2006.16985.
  * [62] U. Chabaud, G. Roeland, M. Walschaers, F. Grosshans, V. Parigi, D. Markham, and N. Treps, “Certification of non-Gaussian states with operational measurements,” PRX Quantum 2, 020333 (2021).
  * [63] S. Bravyi and A. Kitaev, “Universal quantum computation with ideal Clifford gates and noisy ancillas,” Phys. Rev. A 71, 022316 (2005).
  * [64] M. Yoganathan, R. Jozsa, and S. Strelchuk, “Quantum advantage of unitary Clifford circuits with magic state inputs,” Proceedings of the Royal Society A 475, 20180427 (2019).
  * [65] E. Knill, R. Laflamme, and G. J. Milburn, “A scheme for efficient quantum computation with linear optics,” Nature 409, 46–52 (2001).
  * [66] N. C. Menicucci, P. Van Loock, M. Gu, C. Weedbrook, T. C. Ralph, and M. A. Nielsen, “Universal quantum computation with continuous-variable cluster states,” Phys. Rev. Lett. 97, 110501 (2006).
  * [67] M. Gu, C. Weedbrook, N. C. Menicucci, T. C. Ralph, and P. van Loock, “Quantum computing with continuous-variable clusters,” Phys. Rev. A 79, 062318 (2009).
  * [68] B. Q. Baragiola, G. Pantaleoni, R. N. Alexander, A. Karanjai, and N. C. Menicucci, “All-Gaussian universality and fault tolerance with the Gottesman-Kitaev-Preskill code,” Phys. Rev. Lett. 123, 200502 (2019).
  * [69] B. C. Hall, “Holomorphic methods in analysis and mathematical physics. First Summer School in Analysis and Mathematical Physics,(Cuernavaca Morelos, 1998), 1–59,” Contemp. Math 260, , arXiv:quant-ph/9912054.
  * [70] P. Leboeuf and A. Voros, “Chaos-revealing multiplicative representation of quantum eigenstates,” Journal of Physics A: Mathematical and General 23, 1765 (1990).
  * [71] D. Ellinas and V. Kovanis, “Motion of the wave-function zeros in spin-boson systems,” Phys. Rev. A 51, 4230 (1995).
  * [72] R. L. Hudson, “When is the Wigner quasi-probability density non-negative?,” Reports on Mathematical Physics 6, 249–252 (1974).
  * [73] F. Soto and P. Claverie, “When is the Wigner function of multidimensional systems nonnegative?,” Journal of Mathematical Physics 24, 97–100 (1983).
  * [74] N. Lütkenhaus and S. M. Barnett, “Nonclassical effects in phase space,” Phys. Rev. A 51, 3340 (1995).
  * [75] U. Chabaud, G. Ferrini, F. Grosshans, and D. Markham, “Classical simulation of Gaussian quantum circuits with non-Gaussian input states,” Phys. Rev. Research 3, 033018 (2021).
  * [76] O. Bournez and M. L. Campagnolo, “A survey on continuous time computations,” in New computational paradigms, pp. 383–423.  Springer, 2008.
  * [77] L. G. Valiant, “Completeness classes in algebra,” in Proceedings of the eleventh annual ACM symposium on Theory of computing, pp. 249–261.  1979\.
  * [78] L. Blum, M. Shub, and S. Smale, “On a theory of computation and complexity over the real numbers: NP-completeness, recursive functions and universal machines,” in The Collected Papers of Stephen Smale: Volume 3, pp. 1293–1338.  World Scientific, 2000.
  * [79] V. Bush, “The differential analyzer. A new machine for solving differential equations,” Journal of the Franklin Institute 212, 447–488 (1931).
  * [80] S. Smale, “Mathematical problems for the next century,” The mathematical intelligencer 20, 7–15 (1998).
  * [81] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, “From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics,” Advances in Physics 65, 239–362 (2016).
  * [82] S. Aaronson, A. Bouland, G. Kuperberg, and S. Mehraban, “The computational complexity of ball permutations,” in Proceedings of the 49th Annual ACM SIGACT Symposium on Theory of Computing, pp. 317–327.  2017\.
  * [83] K. Fukui and S. Takeda, “Building a large-scale quantum computer with continuous-variable optical technologies,” arXiv:2110.03247.
  * [84] D. Menzies and R. Filip, “Gaussian-optimized preparation of non-Gaussian pure states,” Phys. Rev. A 79, 012313 (2009).
  * [85] R. Gomes, A. Salles, F. Toscano, P. S. Ribeiro, and S. Walborn, “Quantum entanglement beyond Gaussian criteria,” Proceedings of the National Academy of Sciences 106, 21517–21520 (2009).
  * [86] A. Ourjoumtsev, F. Ferreyrol, R. Tualle-Brouri, and P. Grangier, “Preparation of non-local superpositions of quasi-classical light states,” Nature Physics 5, 189–192 (2009).
  * [87] C. M. Dawson and M. A. Nielsen, “The Solovay–Kitaev algorithm,” arXiv:quant-ph/0505030.
  * [88] S. Becker, N. Datta, L. Lami, and C. Rouzé, “Energy-Constrained Discrimination of Unitaries, Quantum Speed Limits, and a Gaussian Solovay–Kitaev Theorem,” Phys. Rev. Lett. 126, 190504 (2021).
  * [89] K. E. Cahill and R. J. Glauber, “Density operators and quasiprobability distributions,” Physical Review 177, 1882 (1969).
  * [90] A. Zavatta, S. Viciani, and M. Bellini, “Quantum-to-classical transition with single-photon-added coherent states of light,” Science 306, 660–662 (2004).
  * [91] G. S. Agarwal, “Quantum optics,”.  Cambridge University Press, 2012.
  * [92] M. Reck, A. Zeilinger, H. J. Bernstein, and P. Bertani, “Experimental realization of any discrete unitary operator,” Phys. Rev. Lett. 73, 58 (1994).
  * [93] D. J. Korteweg and G. De Vries, “XLI. On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves,” The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 39, 422–443 (1895).
  * [94] H. Airault, H. McKean, and J. Moser, “Rational and elliptic solutions of the Korteweg–de Vries equation and a related many-body problem,” Communications on Pure and Applied Mathematics 30, 95–148 (1977).
  * [95] B. B. Kadomtsev and V. I. Petviashvili, “On the stability of solitary waves in weakly dispersing media,” in Sov. Phys. Dokl, vol. 15, pp. 539–541.  1970\.
  * [96] I. M. Krichever, “Rational solutions of the Kadomtsev–Petviashvili equation and integrable systems of N particles on a line,” Functional Analysis and Its Applications 12, 59–61 (1978).
  * [97] P. Etingof, “Lectures on Calogero-Moser systems,”.
  * [98] U. Shackerley-Bennett, A. Pitchford, M. G. Genoni, A. Serafini, and D. K. Burgarth, “The reachable set of single-mode quadratic Hamiltonians,” Journal of Physics A: Mathematical and Theoretical 50, 155203 (2017).
  * [99] U. J. Le Verrier, “Mémoire sur les variations séculaires des éléments des orbites: pour les sept planètes principales, Mercure, Vénus, la Terre, Mars, Jupiter, Saturne et Uranus,”.  Bachelier, 1845.
  * [100] D. K. Faddeev and I. S. Sominskiĭ, “Problems in Higher algebra,”.  Mir Publishers, 1972.
  * [101] M. Abramowitz and I. A. Stegun, “Handbook of mathematical functions: with formulas, graphs, and mathematical tables,”, vol. 55.  Courier Corporation, 1965.
  * [102] A. Botero and B. Reznik, “Modewise entanglement of Gaussian states,” Phys. Rev. A 67, 052311 (2003).
  * [103] G. Giedke, J. Eisert, J. Cirac, and M. Plenio, “Entanglement transformations of pure Gaussian states,” Quantum Information & Computation 3, 211–223 (2003), arXiv:quant-ph/0301038.
  * [104] M. Walschaers, C. Fabre, V. Parigi, and N. Treps, “Statistical signatures of multimode single-photon-added and-subtracted states of light,” Phys. Rev. A 96, 053835 (2017).
  * [105] U. Leonhardt, “Essential Quantum Optics,”.  Cambridge University Press, Cambridge, UK, 1st ed., 2010.
  * [106] G. Giedke and J. I. Cirac, “Characterization of Gaussian operations and distillation of Gaussian states,” Phys. Rev. A 66, 032316 (2002).
  * [107] J. Zhang and S. L. Braunstein, “Continuous-variable Gaussian analog of cluster states,” Phys. Rev. A 73, 032318 (2006).