# Witnessing Wigner Negativity

Ulysse Chabaud∗<EMAIL_ADDRESS>Institute for Quantum Information and
Matter, Caltech  Université de Paris, IRIF, CNRS, France  Pierre-Emmanuel
Emeriau∗<EMAIL_ADDRESS>Sorbonne Université, CNRS, LIP6,
F-75005 Paris, France  Frédéric Grosshans<EMAIL_ADDRESS>Sorbonne
Université, CNRS, LIP6, F-75005 Paris, France

###### Abstract

††∗ These authors contributed equally.

Negativity of the Wigner function is arguably one of the most striking non-
classical features of quantum states. Beyond its fundamental relevance, it is
also a necessary resource for quantum speedup with continuous variables. As
quantum technologies emerge, the need to identify and characterize the
resources which provide an advantage over existing classical technologies
becomes more pressing. Here we derive witnesses for Wigner negativity of
single-mode and multimode quantum states, based on fidelities with Fock
states, which can be reliably measured using standard detection setups. They
possess a threshold expectation value indicating whether the measured state
has a negative Wigner function. Moreover, the amount of violation provides an
operational quantification of Wigner negativity. We phrase the problem of
finding the threshold values for our witnesses as an infinite-dimensional
linear optimisation problem. By relaxing and restricting the corresponding
linear programs, we derive two converging hierarchies of semidefinite
programs, which provide numerical sequences of increasingly tighter upper and
lower bounds for the threshold values. We further show that our witnesses form
a complete family—each Wigner negative state is detected by at least one
witness—thus providing a reliable method for experimentally witnessing Wigner
negativity of quantum states from few measurements. From a foundational
perspective, our findings provide insights on the set of positive Wigner
functions which still lacks a proper characterisation.

##  1 Introduction

Quantum information with continuous variables [1]—where information is encoded
in continuous degrees of freedom of quantum systems—is one of the promising
directions for the future of quantum technologies. For example, continuous-
variable quantum optics enables the deterministic experimental preparation of
entangled states over millions of modes [2] and also offers reliable and
efficient detection methods, such as homodyne or heterodyne detection [3].
From a theoretical point of view, quantum information with continuous
variables provides different perspectives from quantum information with
discrete variables and is described via the formalism of infinite-dimensional
Hilbert spaces.

To handily manipulate states in those infinite spaces, mathematical tools
initially inspired by physics have been developed such as phase-space
formalism [4]. In this framework, quantum states are represented by a
quasiprobability distribution over phase space, like the Wigner function [5].
These representations provide a geometric intuition of quantum states [6]:
quantum states are separated into two categories, Gaussian and non-Gaussian,
depending on whether their Wigner function is a Gaussian function or not.

Non-Gaussian quantum states are essential to a variety of quantum information
processing tasks such as quantum state distillation [7, 8, 9], quantum error-
correction [10], universal quantum computing [1, 11] or quantum computational
speedup [12, 13]. Within those, an important subclass of non-Gaussian states
are the states which display negativity in the Wigner function. These two
classes of states coincide for pure states—namely, non-Gaussian pure states
have a negative Wigner function—as pure states with a positive Wigner function
are necessarily Gaussian states by Hudson theorem [14, 15]. However, this is
not the case for mixed states and the (convex) set of states with a positive
Wigner function becomes much harder to characterise [16, 17].

The negativity of other phase-space quasiprobability distributions has been
used to define different notions of quantumness [18]. In particular, the
negativity of the Glauber–Sudarshan $P$ quasiprobability distribution of a
quantum state is known as its non-classicality [19]. The Wigner function can
be obtained from the more singular $P$ function by a Gaussian convolution,
thus positivity of the latter implies positivity of the former. In particular,
negativity of the Wigner function implies non-classicality, although there are
non-classical states with positive Wigner function, such as squeezed Gaussian
states.

In addition to its fundamental relevance as a non-classical property of
physical systems [20], Wigner negativity is also essential for quantum
computing, since continuous-variable quantum computations described by
positive Wigner functions can be simulated efficiently classically [21].
Wigner negativity is thus a necessary resource, though not sufficient [22],
for quantum computational speedup with continuous variables.

With the rapid development of quantum technologies [23], finding efficient
methods for assessing the correct functioning of quantum devices is of timely
importance [24]. Detecting key properties of a quantum state, such as
entanglement or Wigner negativity, can be done by a full tomographic
reconstruction [25]. However, such reconstructions are very costly in terms of
the number of measurements needed, and require performing a tomographically
complete set of measurements—thus usually involving multiple measurement
settings. In the continuous-variable setting the task is even more daunting,
since the Hilbert space of quantum states is infinite-dimensional [26, 27].

Instead, one may introduce witnesses for specific properties of quantum states
[28, 29, 30, 31, 32] that are easier to measure experimentally. These
witnesses should be observables that possess a threshold expectation value
indicating whether the measured state exhibits the desired property or not.
Intuitively, a witness for a given property can be thought of as a separating
hyperplane in the set of quantum states, such that any state on one side of
this hyperplane has this property (see Fig. 1). In particular, some states
with the sought property may remain unnoticed by the witness. In that regard,
one may use a complete set of witnesses such that for each state exhibiting
the desired property, there exists at least one witness in the set that
captures it.

Figure 1: Pictorial representation of a witness $\hat{\Omega}$ with threshold
value $\omega$ for a given property. In yellow: states without the property.
In blue: states with the property. In red: witness threshold value. In light
blue: states with the property undetected by the witness.

Here, we introduce and study a complete family of witnesses for Wigner
negativity of single-mode quantum states, expressed using Fock states
projectors. The expectation values of our witnesses are linear functions of
the state which may be efficiently estimated experimentally using standard
homodyne or heterodyne detection, thus providing a reliable method for
detecting Wigner negativity with certifiable bounds. Additionally, we show
that the amount by which the measured expectation value exceeds the threshold
value of the witness provides an operational measure of Wigner negativity: it
directly lower bounds the distance between the measured state and the set of
states with positive Wigner function.

We cast the computation of the threshold values of the witnesses as infinite-
dimensional linear programs, which can be either relaxed or restricted. Upper
and lower bounds for the threshold values of our witnesses are then given by
two converging hierarchies of finite-dimensional semidefinite programs,
similar in spirit to the Lasserre–Parrilo hierarchy [33, 34] and the
subsequent Lasserre hierarchy [35].

Finally, we discuss the generalisation to multimode quantum states and show
that most of our results are also applicable in this case.

Our work is thought to be relevant for physicists interested in characterising
Wigner negativity of quantum states—either theoretically or experimentally—and
mathematicians interested in infinite-dimensional convex optimisation theory.
To that end, the rest of the paper is structured as follows: we give some
notations and background in the next section 2 before a detailed exposition of
our witnesses in section 3. Section 4 describes the experimental procedure for
witnessing Wigner negativity of a quantum state using these witnesses,
together with use-case examples. The following section 5—which deals with
infinite-dimensional optimisation techniques of independent interest—is
devoted to estimating the threshold values of our witnesses: after some
technical background in section 5.2, section 5.3 reformulates the problem of
finding the threshold value of a witness as an infinite-dimensional linear
optimisation, while section 5.4 derives two hierarchies of semidefinite
relaxations and restrictions for this linear program, yielding numerical upper
and lower bounds for the threshold value. Section 5.5 establishes the proof of
convergence of these hierarchies of upper and lower bounds to the threshold
values in their respective optimisation spaces. We introduce the
generalisation to the multimode case in section 6 and conclude with a few open
questions in section 7.

##  2 Notations and background

###  2.1 Preliminary notations

For all $m\in\mathbb{N}^{*}$, Symm denotes the space of $m\times m$ real
symmetric matrices. An exponent $T$ denotes the transpose while an exponent
${\dagger}$ denotes the conjugate transpose.

$\mathcal{H}$ denotes a separable infinite-dimensional Hilbert space equipped
with a countable orthonormal single-mode Fock basis
$\\{\ket{n}\\}_{n\in\mathbb{N}}$. We write $\mathcal{D}(\mathcal{H})$ the set
of quantum states (positive semidefinite operators with unit trace) over
$\mathcal{H}$. A single-mode quantum state $\rho$ can be expanded in Fock
basis as $\rho=\sum_{k,l=0}^{+\infty}\rho_{kl}\ket{k}\\!\bra{l}$.

The fidelity between two quantum states $\rho$ and $\sigma$ is denoted
$F(\rho,\sigma)=\text{Tr}(\sqrt{\sqrt{\rho}\sigma\sqrt{\rho}})^{2}$. When one
of the states is pure, it reduces to $F(\rho,\sigma)=\text{Tr}(\rho\sigma)$.
The trace distance between two quantum states $\rho$ and $\sigma$ is denoted
$D(\rho,\sigma)=\frac{1}{2}\text{Tr}(\sqrt{(\rho-\sigma)^{2}})$. The trace
distance can be related to the maximum probability of distinguishing between
two quantum states. The fidelity and trace distance are related by $1-F\leq
D\leq\sqrt{1-F}$ [36].

We denote by $\hat{a}$ and $\hat{a}^{\dagger}$ the single-mode annihilation
and creation operators, respectively, defined by their action on the Fock
basis:

$\displaystyle\hat{a}\ket{n}=\sqrt{n}\ket{n-1},\quad\quad\quad\text{for
}n\in\mathbb{N}^{*},$ (1) $\displaystyle\hat{a}\ket{0}=0,$
$\displaystyle\hat{a}^{\dagger}\ket{n}=\sqrt{n+1}\ket{n+1},\;\quad\text{for
}n\in\mathbb{N}.$

These operators satisfy the canonical commutation relation
$[\hat{a},\hat{a}^{\dagger}]=\mymathbb 1$. For all $\alpha\in\mathbb{C}$, we
write

$\hat{D}(\alpha)=e^{\alpha\hat{a}^{\dagger}-\alpha^{*}\hat{a}}$ (2)

the displacement operator of amplitude $\alpha$ [37]. The coefficients of the
displacement operator in Fock basis are given by [38]:

$\displaystyle\braket{k}{\hat{D}(\alpha)}{l}=e^{-\frac{1}{2}|\alpha|^{2}}$ (3)
$\displaystyle\quad\times\sum_{p=0}^{\min
k,l}\frac{\sqrt{k!l!}(-1)^{l-p}}{p!(k-p)!(l-p)!}\alpha^{k-p}\alpha^{*l-p},$

for all $k,l\in\mathbb{N}$ and all $\alpha\in\mathbb{C}$.

###  2.2 Wigner function

The Wigner function of a single-mode quantum state $\rho$ is an equivalent
representation of the state in phase space which can be expressed as [39, 40]:

$W_{\rho}(\alpha)=\frac{2}{\pi}\text{Tr}\left[\hat{D}(\alpha)\hat{\Pi}\hat{D}^{\dagger}(\alpha)\rho\right],$
(4)

for all $\alpha\in\mathbb{C}$, where $\hat{D}$ is defined above and
$\hat{\Pi}=(-1)^{\hat{a}^{\dagger}\hat{a}}=\sum_{n\geq
0}{(-1)^{n}\ket{n}\\!\bra{n}}$ is the parity operator. In particular, the
Wigner function of a quantum state is related to the expectation value of
displaced parity operators.

The Wigner function is a real-valued quasi-probability distribution [41],
i.e., a normalised distribution which can take negative values. Hence, it
cannot be sampled directly experimentally. However, its marginals are
probability distributions which can be sampled using homodyne detection [26].
Alternatively, heterodyne detection (also called double homodyne detection)
allows for sampling from a smoothed version of the Wigner function [42, 43].
In both cases, applying a displacement before the detection is equivalent to
measuring directly with homodyne or heterodyne detection and applying a
classical post-processing procedure—namely, a translation of the classical
outcome according to the displacement amplitude [44, 32].

As mentioned in the introduction, continuous-variable quantum states are
classified in two categories, Gaussian and non-Gaussian, depending on the
shape of their Wigner function. The set of Gaussian states is well-understood
[45] but has a limited power, while characterising the set of non-Gaussian
states is an active research topic [46, 47, 48, 49].

The negativity of the Wigner function can only decrease under Gaussian
operations [46], i.e., operations that map Gaussian states to Gaussian states.
In particular, it is invariant under displacements. It is also a robust
property, since two almost indistinguishable quantum states have similar
Wigner functions. An operational measure of Wigner negativity for a quantum
state $\rho\in\mathcal{D}(\mathcal{H})$ is given by its distance with the set
of states having a positive Wigner function [30]:

$\eta_{\rho}=\operatorname*{\vphantom{p}inf}_{\begin{subarray}{c}\sigma\in\mathcal{D}(\mathcal{H})\\\
W_{\sigma}\geq 0\end{subarray}}D(\rho,\sigma),$ (5)

where $D$ denotes the trace distance, thus quantifying the operational
distinguishability between the state $\rho$ and any state having a positive
Wigner function [36].

A natural choice for a witness of Wigner negativity is the fidelity with a
pure state having a Wigner function with negative values, since it is a
quantity that can be accessed experimentally by direct fidelity estimation
[25, 27]. Building on this intuition, and given that all Fock states—with the
exception of the (Gaussian) vacuum state $\ket{0}$—have a negative Wigner
function, we introduce in the following section a broad family of Wigner
negativity witnesses for single-mode continuous-variable quantum states based
on fidelities with Fock states.

##  3 Wigner negativity witnesses

We introduce the following Wigner negativity witnesses:

$\hat{\Omega}_{\bm{a},\alpha}:=\sum_{k=1}^{n}a_{k}\hat{D}(\alpha)\ket{k}\\!\bra{k}\hat{D}^{\dagger}(\alpha),$
(6)

for $n\in\mathbb{N}^{*}$, $\bm{a}=(a_{1},\dots,a_{n})\in[0,1]^{n}$, with
$\max_{k}a_{k}=1$, and $\alpha\in\mathbb{C}$. These operators are weighted
sums of displaced Fock states projectors. They can be thought of as Positive
Operator-Valued Measure (POVM) elements, and their expectation value for a
quantum state $\rho$ is given by

$\text{Tr}\\!\left(\hat{\Omega}_{\bm{a},\alpha}\,\rho\right)=\sum_{k=1}^{n}a_{k}F\\!\left(\hat{D}^{\dagger}(\alpha)\rho\hat{D}(\alpha),\ket{k}\right),$
(7)

where $F$ is the fidelity. This quantity can be directly estimated from
homodyne or heterodyne detection of multiple copies of the state $\rho$ by
translating the samples obtained by the amplitude $\alpha$ in the classical
postprocessing and performing fidelity estimation with the Fock states
$\ket{1},\dots,\ket{n}$ [25, 27].

For $n\in\mathbb{N}^{*}$, each choice of
$(\bm{a},\alpha)\in[0,1]^{n}\times\mathbb{C}$ yields a different Wigner
negativity witness. In particular, when $\alpha=0$ and one entry of the vector
$\bm{a}$ is equal to $1$ and all the other entries are $0$, the expectation
value of the witness is given by the fidelity with a single Fock state.

To each witness $\hat{\Omega}_{\bm{a},\alpha}$ is associated a threshold value
defined as:

$\omega_{\bm{a}}:=\sup_{\begin{subarray}{c}\rho\in\mathcal{D}(\mathcal{H})\\\
W_{\rho}\geq
0\end{subarray}}\text{Tr}\\!\left(\hat{\Omega}_{\bm{a},\alpha}\,\rho\right).$
(8)

Since negativity of the Wigner function is invariant under displacements, the
threshold values do not depend on the value of the displacement amplitude
$\alpha$ and we thus write $\omega_{\bm{a}}$ (rather than
$\omega_{\bm{a},\alpha}$) for the threshold value associated to the witness
$\hat{\Omega}_{\bm{a},\alpha}$. This is sensible, given that the threshold
value asks for non-negativity anywhere in phase space, so a displacement in
phase space should not change its value. Combining  (i) that the threshold
value associated to a witness does not depend on the displacement parameter
$\alpha$ and  (ii) that we can always take into account displacement via
classical post-processing if one uses homodyne or heterodyne detection
associated to $\hat{\Omega}_{\bm{a},0}$ [44, 32], we can restrict the analysis
to witnesses of the form $\hat{\Omega}_{\bm{a},0}$ that will generate the
family $\\{\hat{\Omega}_{\bm{a},\alpha}\\}_{\alpha\in\mathbb{C}}$. Note
however that the choice of displacement amplitude can play an important role
for certifying negativity of certain quantum states.

If the measured expectation value for an experimental state is higher than the
threshold value given by Eq. (8), this implies by definition that its Wigner
function takes negative values. Moreover, the following result shows that the
amount by which the expectation value exceeds the threshold value directly
provides an operational quantification of Wigner negativity for that state:

######  Lemma 1.

Let $\rho\in\mathcal{D}(\mathcal{H})$ Wigner negative, and fix a witness
$\hat{\Omega}_{\bm{a},\alpha}$ defined in Eq. (6), for $n\in\mathbb{N}^{*}$,
$\bm{a}=(a_{1},\dots,a_{n})\in[0,1]^{n}$, and $\alpha\in\mathbb{C}$, with
threshold value $\omega_{\bm{a}}$. Let us further assume that that it violates
the threshold value of the witness i.e.
$\text{Tr}(\hat{\Omega}_{\bm{a},\alpha}\,\rho)>\omega_{\bm{a}}$ and denote the
amount of violation as

$\delta_{\bm{a},\alpha}(\rho):=\text{Tr}\\!\left(\hat{\Omega}_{\bm{a},\alpha}\,\rho\right)-\omega_{\bm{a}}.$
(9)

Then,

$\eta_{\rho}\geq\delta_{\bm{a},\alpha}(\rho),$ (10)

where $\eta_{\rho}$ is the distance between $\rho$ and the set of states
having a positive Wigner function, defined in Eq. (5).

###### Proof.

We use the notations of the lemma. Let us consider the binary POVM
$\\{\hat{\Omega}_{\bm{a},\alpha},\mymathbb 1-\hat{\Omega}_{\bm{a},\alpha}\\}$.
For all $\sigma\in\mathcal{D}(\mathcal{H})$, we write
$P^{\sigma}_{\bm{a},\alpha}$ the associated probability distribution:
$P^{\sigma}_{\bm{a},\alpha}(0)=1-P^{\sigma}_{\bm{a},\alpha}(1)=\text{Tr}(\hat{\Omega}_{\bm{a},\alpha}\,\sigma)$.

Let $\sigma$ be a state with a positive Wigner function, so that
$\text{Tr}(\hat{\Omega}_{\bm{a},\alpha}\,\sigma)\leq\omega_{\bm{a}}$, by
definition of the threshold value. We have:

$\displaystyle\delta_{\bm{a},\alpha}(\rho)$
$\displaystyle=\text{Tr}(\hat{\Omega}_{\bm{a},\alpha}\rho)-\omega_{\bm{a}}$
(11)
$\displaystyle\leq|\text{Tr}(\hat{\Omega}_{\bm{a},\alpha}\rho)-\text{Tr}(\hat{\Omega}_{\bm{a},\alpha}\sigma)|$
$\displaystyle=|P^{\rho}_{\bm{a},\alpha}(0)-P^{\sigma}_{\bm{a},\alpha}(0)|$
$\displaystyle=\|P^{\rho}_{\bm{a},\alpha}-P^{\sigma}_{\bm{a},\alpha}\|_{tvd}$
$\displaystyle\leq D(\rho,\sigma),$

where we used $\delta_{\bm{a},\alpha}(\rho)\geq 0$ in the second line,
$\|P-Q\|_{tvd}=\frac{1}{2}\sum_{x}|P(x)-Q(x)|$ denotes the total variation
distance, and we used the operational property of the trace distance in the
last line [36]. With Eq. (5), taking the infimum over $\sigma$ concludes the
proof. ∎

This result directly extends to the case where only an upper bound of the
threshold value is known: the amount by which the expectation value exceeds
this upper bound is also a lower bound of the distance to the set of states
having a positive Wigner function.

Importantly, the family of Wigner negativity witnesses
$\\{\hat{\Omega}_{\bm{a},\alpha}\\}$ is complete, i.e., for any quantum state
with negative Wigner function there exists a choice of witness
$(\bm{a},\alpha)$ such that the expectation value of
$\hat{\Omega}_{\bm{a},\alpha}$ for this state is higher than the threshold
value. Indeed, by taking $\bm{a}=(1,0,1,0,1,\dots)$, this family includes as a
subclass the complete family of witnesses from [32]. Indeed, by expanding the
definition of the Wigner function in Eq. (4) with the expression of the parity
operator and using the completeness relation
$\sum_{n\in\mathbb{N}}\ket{n}\bra{n}=\mathds{1}$, the Wigner function of any
density operator $\rho\in\mathcal{D}(\mathcal{H})$ reads:

$W_{\rho}(\alpha)=\frac{2}{\pi}\left(1-2\text{Tr}\left(\hat{\Omega}_{(1,0,1,0,\dots),\alpha}\rho\right)\right)\text{
.}$

Hence, for any state with a negative Wigner function, there exists a choice of
$\alpha\in\mathbb{C}$ such that the witness
$\hat{\Omega}_{(1,0,1,0,\dots),\alpha}$ with threshold value $\frac{1}{2}$ can
detect its negativity.

The threshold value in Eq. (8) is given by an optimisation problem over
quantum states having a positive Wigner function, which is a convex subset of
an infinite-dimensional space that does not possess a well-characterised
structure. While solving this optimisation problem thus seems unfeasible in
general, it turns out that we can obtain increasingly good numerical upper and
lower bounds for the threshold value using semidefinite programming.

Semidefinite programming is a particular case of conic programming—a subfield
of convex optimisation—where one optimises linear functions within the convex
cone of positive semidefinite matrices [50]. This is a powerful optimisation
technique as semidefinite programs (SDP) can be solved efficiently using
interior point methods.

The relevant programs are derived in section 5.4.2 and 5.4.3 where their
convergence is proven. Since these proofs introduce several intermediate forms
of the programs, we explicitly give them below to avoid confusion on which
programs to implement numerically. For $n\in\mathbb{N}^{*}$,
$\bm{a}=(a_{1},\dots,a_{n})\in[0,1]^{n}$, and $m\geq n$, the hierarchies of
semidefinite programs that respectively provide lower bounds and upper bounds
for the threshold value $\omega_{\bm{a}}$ associated to the witnesses
$\\{\hat{\Omega}_{\bm{a},\alpha}\\}_{\alpha\in\mathbb{C}}$ are:

$\displaystyle\quad\text{Find }Q\in\text{Sym}_{m+1}\text{ and
}\bm{F}\in\mathbb{R}^{m+1}$ $(\text{SDP}^{m,\leq}_{\bm{a}})$
$\displaystyle\quad\text{maximising }\textstyle\sum_{k=1}^{n}a_{k}F_{k}$
subject to $\displaystyle\begin{dcases}\sum_{k=0}^{m}F_{k}=1\\\ \forall
k\in\llbracket 0,m\rrbracket,\;F_{k}\geq 0\\\ \forall l\in\llbracket
1,m\rrbracket,\;\sum_{i+j=2l-1}Q_{ij}=0\\\ \forall l\in\llbracket
0,m\rrbracket,\;\sum_{i+j=2l}Q_{ij}=\sum_{k=l}^{m}\frac{(-1)^{k+l}}{l!}\binom{k}{l}F_{k}\\\
Q\succeq 0,\end{dcases}$

and

$\displaystyle\quad\text{Find }A\in\text{Sym}_{m+1}\text{ and
}\bm{F}\in\mathbb{R}^{m+1}$ $(\text{SDP}^{m,\geq}_{\bm{a}})$
$\displaystyle\quad\text{maximising }\textstyle\sum_{k=1}^{n}a_{k}F_{k}$
subject to $\displaystyle\begin{dcases}\sum_{k=0}^{m}F_{k}=1\\\ \forall
k\in\llbracket 0,m\rrbracket,\;F_{k}\geq 0\\\ \forall l\in\llbracket
1,m\rrbracket,\forall i+j=2l-1,\;A_{ij}=0\\\ \forall l\leq m,\forall
i+j=2l,\;A_{ij}=\sum_{k=0}^{l}F_{k}\binom{l}{k}l!\\\ A\succeq 0,\end{dcases}$

Let $\omega_{\bm{a}}^{m,\geq}$ be the optimal value of
$(\text{SDP}^{m,\geq}_{\bm{a}})$. We show in section 5 that the sequence
$\\{\omega_{\bm{a}}^{m,\geq}\\}_{m\geq n}$ is a decreasing sequence of upper
bounds of $\omega_{\bm{a}}$, which converges to $\omega_{\bm{a}}$. Similarly,
let $\omega_{\bm{a}}^{m,\leq}$ be the optimal value of
$(\text{SDP}^{m,\leq}_{\bm{a}})$. We show that the sequence
$\\{\omega_{\bm{a}}^{m,\leq}\\}_{m\geq n}$ is an increasing sequence of lower
bounds of $\omega_{\bm{a}}$, which converges to
$\omega_{\bm{a}}^{\mathcal{S}}$—a modified threshold value computed with
Schwartz functions rather than square-integrable functions; we have
$\omega_{\bm{a}}^{\mathcal{S}}\leq\omega_{\bm{a}}$, and the equality between
the two values is open.

In particular, the numerical upper bounds $\omega_{\bm{a}}^{m,\geq}$ can be
used instead of the threshold value $\omega_{\bm{a}}$ to witness Wigner
negativity, while the numerical lower bounds $\omega_{\bm{a}}^{m,\leq}$ may be
used to control how much the upper bounds differ from the threshold value. We
give a detailed procedure in the following section, together with use-case
examples and details on the numerical implementation.

Note that although both problems are maximisation problems, the constraints of
$(\text{SDP}^{m,\geq}_{\bm{a}})$ get tighter as $m$ increases, and thus the
corresponding sequence of optimal values is decreasing, while the constraints
of $(\text{SDP}^{m,\leq}_{\bm{a}})$ get looser as m increases, and thus the
corresponding sequence of optimal values is increasing.

##### State-of-the-art

Our Wigner negativity witnesses outperform existing ones [30, 51] in terms of
generality and practicality, since they form a complete family and provide
much more flexibility with the choice of $n\in\mathbb{N}^{*}$,
$\bm{a}\in[0,1]^{n}$ and $\alpha\in\mathbb{C}$. They are accessible with
optical homodyne or heterodyne measurement, and do not require making any
assumption on the measured state, unlike other existing methods [52].
Moreover, our witnesses generalise those of [32], and may provide simpler
alternatives to detect Wigner negativity. We also provide two converging
hierarchies to approximate the threshold values associated to these witnesses.
Proving convergence is important and was not considered in the other
approaches mentioned above. Finally, our approach also generalises to the
multimode setting, as detailed in section 6.

##  4 Witnessing Wigner negativity

###  4.1 Procedure

In this section, mainly devoted to experimentalists, we describe a procedure
to check whether a continuous-variable quantum state exhibits Wigner
negativity using our witnesses.

Witnessing Wigner negativity: 1. Choose a fidelity-based witness
$\hat{\Omega}_{\bm{a},\alpha}$ defined in Eq. (6) by picking
$n\in\mathbb{N}^{*}$, $\bm{a}\in[0,1]^{n}$ and $\alpha\in\mathbb{C}$. 2. Run
the upper bound semidefinite program $(\text{SDP}^{m,\geq}_{\bm{a}})$ for
$m\geq n$, get a numerical estimate $\omega^{m,\geq}_{\bm{a}}$. These values
are already computed for $\bm{a}=(0,\dots,0,1)$ and $n\leq 10$ in Table 1. 3.
Run the lower bound semidefinite program $(\text{SDP}^{m,\leq}_{\bm{a}})$ for
$m\geq n$, get a numerical estimate $\omega^{m,\leq}_{\bm{a}}$. These values
are already computed for $\bm{a}=(0,\dots,0,1)$ and $n\leq 10$ in Table 1. 4.
Estimate the expectation value for that witness of the experimental state from
samples of homodyne or heterodyne detection by translating the samples by
$\alpha$ and performing fidelity estimation with the corresponding Fock
states. This yields an experimental witness value denoted
$\omega_{\text{exp}}$. 5. Compare the value obtained experimentally with the
numerical bounds: if it is greater than the numerical upper bound, i.e.,
$\omega_{\text{exp}}\geq\omega^{m,\geq}_{\bm{a}}$ then the state displays
Wigner negativity, and its distance to the set of Wigner positive states is
lower bounded by $\omega_{\text{exp}}-\omega^{m,\geq}_{\bm{a}}$. Otherwise, if
it is lower than the numerical lower bound, i.e.,
$\omega_{\text{exp}}\leq\omega^{m,\leq}_{\bm{a}}$ then the witness cannot
detect Wigner negativity for this state.

The main subroutine of this procedure is to estimate fidelities with displaced
Fock states using classical samples from homodyne or heterodyne
detection111Actually, using a fidelity witness rather than a fidelity estimate
is sufficient for our purpose., in order to compute the experimental value for
a Wigner negativity witness. Displacement can be achieved with classical post-
processing by translating the classical samples according to the displacement
amplitude, and performing direct fidelity estimation with Fock states (see
section 2.2 and [26, 25, 27]). Upper and lower bounds on the threshold value
of the witness are then obtained using semidefinite programming, and comparing
the experimental witness value to these bounds gives insight about the Wigner
negativity of the measured quantum state.

We give a detailed procedure for using our witnesses for detecting Wigner
negativity in the framed box. This procedure starts by the choice of a
specific witness, and we explain hereafter a heuristic method for picking a
good witness.

If the experimental state is anticipated to have negativity at $\alpha$, then
one may use the witness with parameters $(n,\bm{a},\alpha)$ with
$\bm{a}=(1,0,1,0,\dots)$, which will detect negativity for $n$ large enough
[32]. However, this may imply having to estimate fidelities with Fock states
having a high photon number with homodyne or heterodyne detection, which
requires a lot of samples, while simpler witnesses can suffice for the task
and be more efficient, as we show in the next section. Moreover, there are
cases where the state to be characterised is fully unknown.

Instead, a simple heuristic for picking a good witness for Wigner negativity
is the following:

  * •

From samples of homodyne or heterodyne detection of multiple copies of an
experimental state, estimate the expected values of witnesses in Eq. (6) for a
small value of $n$ and a large set of values $\bm{a}$ and $\alpha$, using the
same samples for all witnesses.

  * •

Based on these values, pick the simplest witness possible—with the smallest
value of $n$—that is able to witness Wigner negativity with a reasonable
violation. This is done by comparing the estimated expected values with the
upper bounds on the corresponding threshold values. These bounds depend only
the choice of the witness parameters $n,\bm{a}$ and can be precomputed using
$(\text{SDP}^{m,\geq}_{\bm{a}})$. To facilitate the use of our methods, we
have collected such bounds for $\bm{a}=(0,\dots,0,1)$ and $n\leq 10$ in Table
1. We also precomputed these bounds for $n=3$ and a large number of values of
$\bm{a}$ in Appendix I and [53].

  * •

Then, estimate the expected value for that witness using a new collection of
samples—thus obtaining proper error bars and avoiding the accumulation of
statistical errors.

In what follows, we give a few theoretical examples for using our witnesses to
detect negativity of the Wigner function of single-mode quantum states.

###  4.2 Examples

We identify three levels of generality within our family of witnesses in Eq.
(6):  (i) fidelities with single Fock states, (ii) linear combinations of
fidelities with Fock states, and  (iii) displaced linear combinations of
fidelities with Fock states.

Fidelities with Fock states are the most practical of our witnesses, since
they require the estimation of only one diagonal element of the density matrix
of the measured state. The corresponding values in Table 1 can be used
directly by experimentalists: if an estimate of $\braket{n}{\rho}{n}$ for
appropriate $n$ is above one of these numerical upper bounds then it ensures
that $\rho$ has a Wigner function with negative values. Moreover, by Lemma 1,
the amount by which the estimate of $\braket{n}{\rho}{n}$ exceeds the
numerical upper bound directly provides a lower bound on the distance between
$\rho$ and the set of states having a positive Wigner function.

For instance, if we focus on $n=3$ in Table 1, the threshold value
$\omega_{3}$ satisfies $0.378\leq\omega_{3}\leq 0.427$. Having a state $\rho$
such that $\braket{3}{\rho}{3}>0.427$ guarantees that $\rho$ has Wigner
negativity. If $\braket{3}{\rho}{3}<0.378$ then we conclude that the witness
cannot detect negativity for this state.

$n$ | Lower bound | Upper bound
---|---|---
1 | 0.5 | 0.5
2 | 0.5 | 0.5
3 | 0.378 | 0.427
4 | 0.375 | 0.441
5 | 0.314 | 0.385
6 | 0.314 | 0.378
7 | 0.277 | 0.344
8 | 0.280 | 0.348
9 | 0.256 | 0.341
10 | 0.262 | 0.334
Table 1: Table of numerical upper and lower bounds for the threshold value
$\omega_{n}$ of various Wigner negativity witnesses obtained using our
hierarchies of semidefinite programs at rank up to around $30$. The witnesses
considered here are Fock states projectors (photon-number states $\ket{n}$)
from $1$ to $10$. Note that the gap between the lower and upper bounds never
exceeds $0.1$. Additionally, the bounds in the first two lines are analytical
(see section 5.3) and for the upper bounds, the corresponding numerical values
are $0.528$ and $0.551$, respectively. See section 4.3 and [53] for the
numerical implementation.

When the experimental state is close to a Fock state (different from the
vacuum), a natural choice for the witness thus is the fidelity with the
corresponding Fock state. For instance, consider a photon-subtracted squeezed
vacuum state [54]

$\ket{\text{p-ssvs}(r)}=\frac{1}{\sinh r}\hat{a}\hat{S}(r)\ket{0},$ (12)

where $\hat{S}(r)=e^{\frac{r}{2}(\hat{a}^{2}-\hat{a}^{{\dagger}2})}$ is a
squeezing operator with parameter $r\in\mathbb{R}$. Its fidelity with the
single-photon Fock state $\ket{1}$ is given by:

$\frac{1}{(\sinh
r)^{2}}\left|\braket{1}{\hat{a}\hat{S}(r)}{0}\right|^{2}=\frac{1}{(\cosh
r)^{3}}.$ (13)

When the squeezing parameter is small, this state is close to a single-photon
Fock state. In particular, for $0<r<0.70$, the fidelity
$F(\text{p-ssvs}(r),1)$ in Eq. (13) is greater than
$\omega_{1}^{\geq}=\frac{1}{2}$ and our witness can be used to detect Wigner
negativity of this state (see Fig. 2 (a)).

Another example is given by superpositions of coherent states states: we
consider the cat state [55]

$\ket{\text{cat}_{2}(\alpha)}=\frac{\ket{\alpha}+\ket{-\alpha}}{\sqrt{2(1+e^{-2|\alpha|^{2}})}},$
(14)

and the compass state [56]

$\ket{\text{cat}_{4}(\alpha)}=\frac{\ket{\alpha}+\ket{-\alpha}+\ket{i\alpha}+\ket{-i\alpha}}{2\sqrt{1+e^{-|\alpha|^{2}}(2\cos(|\alpha|^{2})+1)}},$
(15)

where $\ket{\alpha}=e^{-\frac{1}{2}|\alpha|^{2}}\sum_{k\geq
0}\frac{\alpha^{k}}{\sqrt{k!}}\ket{k}$ is the coherent state of amplitude
$\alpha\in\mathbb{C}$.

(a)  (b)

(c)

Figure 2: (a) Fidelities of photon-subtracted squeezed vacuum states
$\ket{\text{p-ssvs}(r)}$ with squeezing parameter $r\in\mathbb{R}$ with the
Fock state $\ket{1}$. (b) Fidelities of cat states
$\ket{\text{cat}_{2}(\alpha)}$ with amplitude $\alpha\in\mathbb{C}$ with the
Fock state $\ket{2}$. (c) Fidelities of compass states
$\ket{\text{cat}_{4}(\alpha)}$ with amplitude $\alpha\in\mathbb{C}$ with the
Fock state $\ket{4}$. The dashed black lines delimit the intervals of
parameter values where our witnesses from Table 1 can be used to detect Wigner
negativity of the corresponding state, i.e., when the fidelity (purple curve)
is above the witness upper bound (pink horizontal line). When it is below the
witness lower bound, we are guaranteed that the witness (here
$\ket{4}\\!\bra{4}$) cannot be used to detect Wigner negativity of the state.

We have

$\left|\braket{2}{\text{cat}_{2}(\alpha)}\right|^{2}=\frac{|\alpha|^{4}}{2\cosh(|\alpha|^{2})},$
(16)

and

$\left|\braket{4}{\text{cat}_{4}(\alpha)}\right|^{2}=\frac{|\alpha|^{8}/12}{\cosh(|\alpha|^{2})+\cos(|\alpha|^{2})}.$
(17)

For $1.63\leq|\alpha|^{2}\leq 2.59$, the fidelity
$F(\text{cat}_{2}(\alpha),2)$ in Eq. (16) is greater than
$\omega_{2}^{\geq}=\frac{1}{2}$ and our witness corresponding to $n=2$ can be
used to detect Wigner negativity of this state. Similarly, for
$2.10\leq|\alpha|^{2}\leq 6.70$, the fidelity $F(\text{cat}_{4}(\alpha),4)$ in
Eq. (17) is greater than $\omega_{4}^{\geq}=0.441$ and our witness
corresponding to $n=4$ can be used to detect Wigner negativity of this state
(see Fig. 2 (b) and (c)). Note that in all cases, the height difference when
the purple curve is above the pink horizontal line directly provides a lower
bound on the distance between the corresponding state and the set of states
having a positive Wigner function.

Some quantum states will remain unnoticed by all single Fock state negativity
witnesses. For example, the state
$\rho_{0,1,2}:=\frac{1}{9}\ket{0}\\!\bra{0}+\frac{4}{9}\ket{1}\\!\bra{1}+\frac{4}{9}\ket{2}\\!\bra{2}$
has a negative Wigner function but is not detected by any of the single Fock
state negativity witnesses, since the lower bounds for $n=1,2$ in Table 1 are
higher than $\frac{4}{9}$, and this state has fidelity $0$ with higher Fock
states. However, it is detected by simple witnesses based on linear
combinations of fidelities. For example, with $n=2$ and $\bm{a}=(1,1)$ we find
numerically that the threshold value of the witness
$\ket{1}\\!\bra{1}+\ket{2}\\!\bra{2}$ is less than $0.875$ when running the
corresponding program $(\text{SDP}^{m,\geq}_{\bm{a}})$ for $m=7$. And since
$\text{Tr}(\rho_{0,1,2}(\ket{1}\\!\bra{1}+\ket{2}\\!\bra{2}))=\frac{8}{9}>0.875$,
this linear combination of Fock state fidelities can indeed detect Wigner
negativity for this state.

However, some quantum states with a negative Wigner function will always go
unnoticed by the previous witnesses because these witnesses are invariant
under phase-space rotations while the Wigner function of those states becomes
positive under random dephasing: consider for instance the superposition
$\sqrt{1-\frac{1}{s}}\ket{0}+\frac{1}{\sqrt{s}}\ket{1}$, for $s>2$ (which
under random dephasing is mapped to
$(1-\frac{1}{s})\ket{0}\\!\bra{0}+\frac{1}{s}\ket{1}\\!\bra{1}$). In that
case, the Wigner negativity of such states can still be witnessed by using the
displaced version of our witnesses. In particular, if any single-mode quantum
state has a Wigner function negative at $\alpha\in\mathbb{C}$, then there is a
choice of $n\in\mathbb{N}^{*}$ such that the witness in Eq. (6) defined by
$\bm{a}=(1,0,1,0,1,\dots)$ and the displacement amplitude $\alpha$ detects its
negativity [32]. In practice, simpler witnesses may suffice to detect
negativity, and the choice of witness will ultimately depend on the
experimental state at hand.

Figure 3: Lower bounds on the violation
$\delta_{\bm{a},\alpha}=\text{Tr}(\hat{\Omega}_{\bm{a},\alpha}\,\rho_{3,\eta})-\omega_{\bm{a}}$
for a lossy 3-photon Fock state
$\rho_{3,\eta}=\eta^{3}\ket{3}\\!\bra{3}+3\eta^{2}(1-\eta)\ket{2}\\!\bra{2}+3\eta(1-\eta)^{2}\ket{1}\\!\bra{1}+(1-\eta)^{3}\ket{0}\\!\bra{0}$
for $\bm{a}=(a_{1},a_{2},a_{3})$ with the loss parameter $\eta$. Precomputed
bounds on threshold values for witnesses of the form
$\hat{\Omega}_{\bm{a}}=a_{1}\ket{1}\\!\bra{1}+a_{2}\ket{2}\\!\bra{2}+a_{3}\ket{3}\\!\bra{3}$
can be found in Appendix I. We use these values to find the witness
$\hat{\Omega}_{\bm{a},\alpha}$ giving the maximum lower bound
$\delta_{\bm{a},\alpha}^{\leq}=\text{Tr}(\hat{\Omega}_{\bm{a},\alpha}\,\rho_{3,\eta})-\omega_{\bm{a}}^{\geq}$
on the violation $\delta_{\bm{a},\alpha}$, for different values of the loss
parameter $\eta$. These lower bounds are represented in red. In blue is the
maximal violation that can be detected using the witnesses
$\hat{\Omega}_{(1,0,1),\alpha}=\hat{D}(\alpha)(\ket{1}\\!\bra{1}+\ket{3}\\!\bra{3})\hat{D}^{\dagger}(\alpha)$
[32]. In violet is the maximal violation that can be detected using the more
naive witness
$\hat{\Omega}_{(0,0,1),\alpha}=D(\alpha)\ket{3}\\!\bra{3}D^{\dagger}(\alpha)$.
Note that $\rho_{3,\eta}$ has a non-negative Wigner function for $\eta\geq
0.5$.

Hereafter we discuss the heuristics for picking a good witness, with the
theoretical example of the lossy $3$-photon Fock state:

$\displaystyle\begin{split}\rho_{3,\eta}&:=(1-\eta)^{3}\ket{3}\\!\bra{3}+3\eta(1-\eta)^{2}\ket{2}\\!\bra{2}\\\
&\;\;+3\eta^{2}(1-\eta)\ket{1}\\!\bra{1}+\eta^{3}\ket{0}\\!\bra{0},\end{split}$
(18)

where $0\leq\eta\leq 1$ is the loss parameter. Setting $\eta=0$ gives
$\rho_{3,\eta}=\ket{3}\\!\bra{3}$ while setting $\eta=1$ gives
$\rho_{3,\eta}=\ket{0}\\!\bra{0}$. This state has a non-negative Wigner
function for $\eta\geq\frac{1}{2}$. The fidelities of $\rho_{3,\eta}$ with
displaced Fock states $\hat{D}(\alpha)\ket{l}$ are given by:

$\displaystyle F(\rho_{3,\eta},\hat{D}(\alpha)$ $\displaystyle\ket{l})$ (19)
$\displaystyle=(1-\eta)^{3}|\braket{3}{\hat{D}(\alpha)}{l}|^{2}$
$\displaystyle+3\eta(1-\eta)^{2}|\braket{2}{\hat{D}(\alpha)}{l}|^{2}$
$\displaystyle+3\eta^{2}(1-\eta)|\braket{1}{\hat{D}(\alpha)}{l}|^{2}$
$\displaystyle+\eta^{3}|\braket{0}{\hat{D}(\alpha)}{l}|^{2}.$

where the coefficients of the displacement operator in Fock basis are given in
Eq. (3). In an experimental scenario, the state would be unknown and these
fidelities should be estimated using samples from a homodyne or heterodyne
detection of the state translated by $\alpha$, and estimating the fidelities
with Fock states [44, 32].

Following the heuristic detailed in the previous section, we have determined
good Wigner negativity witnesses for $50$ values of the loss parameter $\eta$
between $0$ and $0.5$ as follows: for each value $\eta$, we have computed
numerically the values of the fidelities in Eq. (19) for $l=1,2,3$, and for
displacement parameters $\alpha=q/10+ip/10$ for all $q,p\in\llbracket
0,10\rrbracket$. Using these values, we have computed the expectation value of
the witnesses $\hat{\Omega}_{\bm{a},\alpha}$ for multiple choices of
$\bm{a}=(a_{1},a_{2},a_{3})$ with $\max_{i}a_{i}=1$. We have used the
corresponding precomputed bounds $\omega_{\bm{a}}^{\geq}$ on the threshold
values of $\hat{\Omega}_{\bm{a},\alpha}$ in Appendix I to determine the
witness leading to the maximal lower bound
$\delta_{\bm{a},\alpha}^{\leq}:=\text{Tr}(\hat{\Omega}_{\bm{a},\alpha}\,\rho_{3,\eta})-\omega_{\bm{a}}^{\geq}$
on the violation
$\delta_{\bm{a},\alpha}=\text{Tr}(\hat{\Omega}_{\bm{a},\alpha}\,\rho_{3,\eta})-\omega_{\bm{a}}$
over the choice of $(\bm{a},\alpha)$.

We have represented these violations for each value of the loss parameter
$\eta$ in Fig. 3. For all values of $\eta$, we find that the optimal
displacement parameter is $\alpha=0$. On the other hand, we find different
optimal choices of $\bm{a}$ for different values of $\eta$. To illustrate the
usefulness of the optimisation over the choice of witnesses parametrised by
$(\bm{a},\alpha)$, we have also represented the violations obtained when using
the witnesses
$\hat{\Omega}_{(1,0,1),\alpha}=\hat{D}(\alpha)(\ket{1}\\!\bra{1}+\ket{3}\\!\bra{3})\hat{D}^{\dagger}(\alpha)$
from [32] for all values of $\eta$. In that setting, the violations obtained
quantify how hard it is to detect the Wigner negativity of the state: a larger
violation implies that a less precise estimate of the witness expectation
value is needed to witness Wigner negativity. In particular, we obtain that
our optimised witnesses always provide a greater violation to detect
negativity than the previous witnesses which will result in an easier
experimental detection. We also represented the violation obtained when using
the more naive witnesses
$\hat{\Omega}_{(0,0,1),\alpha}=\hat{D}(\alpha)\ket{3}\\!\bra{3}\hat{D}^{\dagger}(\alpha)$
and we see that it is only useful when the loss parameter is smaller than
$0.25$, while the optimised witnesses may detect negativity of the state
$\rho_{3,\eta}$ up to $\eta=0.5$—when the Wigner functions becomes non-
negative—provided the estimates of the fidelities are precise enough.

Overall, this procedure only amounts to a simple classical post-processing of
samples from homodyne or heterodyne detection and yields a good witness for
detecting Wigner negativity.

###  4.3 Numerical implementation

Here we discuss numerical implementations of the semidefinite programs
$(\text{SDP}^{m,\leq}_{\bm{a}})$ and $(\text{SDP}^{m,\geq}_{\bm{a}})$. All
codes are available here [53].

We implemented the semidefinite programs with Python through the interface
provided by PICOS [57]. We first used the solver Mosek [58] to solve these
problems but, while the size of the semidefinite programs remains relatively
low for small values of $n$ and $m$, binomial terms grow rapidly and numerical
precision issues arise quickly (usually for $m=12$, $n\leq m$). The linear
constraints involving $Q_{ij}$ in the semidefinite programs come from a
polynomial equality (see Lemma 3). While polynomial equalities are usually
written in the canonical basis, a first trick is to express them in a
different basis—for instance the basis
$(1,\frac{X}{1!},\frac{X^{2}}{2!},\dots)$—to counterbalance the binomial
terms.

However, this may not be sufficient to probe larger values of $m$. Instead, we
used the solver SDPA-GMP [59, 60] which allows arbitrary precision arithmetic.
While much slower, this solver is dedicated to solve problems requiring a lot
of precision. Because our problems remain rather small, time efficiency is not
an issue and this solver is particularly well-suited. All problems were
initially solved on a regular laptop as warning flags on optimality were
raised before the problems were too large. A high-performance computer222DELL
PowerEdge R440, 384 Gb RAM, Intel Xeon Silver 4216 processor, 64 threads from
LIP6. —handling floating point arithmetic more accurately—was later used to
compute further ranks in the hierarchy.

Using the semidefinite programs $(\text{SDP}^{m,\geq}_{\bm{a}})$ and
$(\text{SDP}^{m,\leq}_{\bm{a}})$ for values of $m$ up to around $30$ and
$\bm{a}=(0,0,\dots,0,1)$ (where the size $n$ of the vector $\bm{a}$ is ranging
from $1$ to $10$), we have obtained upper and lower bounds for the threshold
values of Wigner negativity witnesses corresponding to fidelities with Fock
states from $1$ to $10$, reported in Table 1.

We also computed upper and lower bounds on the threshold values of witnesses
of the form:

$\hat{\Omega}_{(a_{1},a_{2},a_{3})}=\sum_{k=1}^{3}a_{k}\ket{k}\\!\bra{k},$
(20)

where $\forall i\in\\{1,2,3\\},\,0\leq a_{i}\leq 1$ and $\max_{i}a_{i}=1$. We
focused on these particular witnesses for experimental considerations as it is
challenging to obtain fidelities with higher Fock states. We fix one
coefficient equal to $1$ and vary each other $a_{i}$ from $0$ to $1$ with a
step of $0.1$. The resulting bounds on the threshold values can be found in
appendix I.

We now turn to the mathematical proofs of our results, i.e., that the
threshold values in Eq. (8) can be upper bounded and lower bounded by the
optimal values of the converging hierarchies of semidefinite programs
$(\text{SDP}^{m,\geq}_{\bm{a}})$m≥n and $(\text{SDP}^{m,\leq}_{\bm{a}})$m≥n,
respectively.

The following section is rather technical as we dive into infinite-dimensional
optimisation techniques to prove the convergence of the two hierarchies of
semidefinite programs. Some readers may want to skip directly to section 6.

##  5 Infinite-dimensional optimisation

In this section we use infinite-dimensional optimisation techniques:  (i) to
phrase the computation of the witness threshold value introduced in Eq. (8) as
an infinite-dimensional linear program in section 5.3,  (ii) to derive two
hierarchies of finite-dimensional semidefinite programs that upper bound and
lower bound the threshold value in section 5.4, and  (iii) to show, in section
5.5, that the sequence of upper bounds converges to the threshold value
computed over square-integrable functions and the sequence of lower bounds
converges to the threshold value computed over Schwartz functions (see Fig.
4). Given the technicalities of the proofs above, we sketch them in section
5.1 before detailing them in the following sections.

As a convention, except if specifically mentioned, we will use the terminology
‘relaxation’ and ‘restriction’ from the point of view of the primal program.
We will refer to the hierarchy of semidefinite programs providing the upper
bounds as a hierarchy of relaxations because the obtained semidefinite
programs are indeed relaxations of the primal program (while they are
restrictions of the dual program). Likewise, we will refer to the hierarchy
providing the lower bounds as a hierarchy of restrictions.

01$\omega_{n}^{m,\leq}$$\omega_{n}^{L^{2}}$$\omega_{n}^{\mathcal{S}}$$\omega_{n}^{m,\geq}$$(\text{SDP}^{m,\geq}_{n})$$(\text{D-SDP}^{m,\geq}_{n})$(LP${}_{n}^{L^{2}}$)(D-LP${}_{n}^{L^{2}}$)(LP${}_{n}^{\mathcal{S}}$)(D-LP${}_{n}^{\mathcal{S}}$)$(\text{SDP}^{m,\leq}_{n})$$(\text{D-SDP}^{m,\leq}_{n})$$\begin{subarray}{c}\text{Theorem}\leavevmode\nobreak\
\ref{th:sdupper}\\\
=\joinrel=\end{subarray}$$\begin{subarray}{c}\text{Theorem}\leavevmode\nobreak\
\ref{th:sdLP}\\\
=\joinrel=\end{subarray}$$\begin{subarray}{c}\text{Theorem}\leavevmode\nobreak\
\ref{th:sdLP}\\\
=\joinrel=\end{subarray}$$\begin{subarray}{c}\text{Theorem}\leavevmode\nobreak\
\ref{th:sdlower}\\\
=\joinrel=\end{subarray}$$\begin{subarray}{c}\text{Theorem\leavevmode\nobreak\
\ref{th:upperCV}}\\\ \hskip
20.0pt\big{\downarrow}m\rightarrow\infty\end{subarray}$$\begin{subarray}{c}\hskip
20.0pt\big{\uparrow}m\rightarrow\infty\\\ \text{Theorem\leavevmode\nobreak\
\ref{th:lowerCV}}\end{subarray}$$\begin{subarray}{c}\text{Theorem\leavevmode\nobreak\
\ref{th:upperCV}}\\\ \hskip
20.0pt\big{\downarrow}m\rightarrow\infty\end{subarray}$$\begin{subarray}{c}\hskip
20.0pt\big{\uparrow}m\rightarrow\infty\\\ \text{Theorem\leavevmode\nobreak\
\ref{th:lowerCV}}\end{subarray}$ Figure 4: Hierarchy of semidefinite
relaxations converging to the linear program (LP${}_{n}^{L^{2}}$) and
hierarchy of semidefinite restrictions converging to the linear program
(LP${}_{n}^{\mathcal{S}}$), together with their dual programs. The upper index
$m$ denotes the level of the relaxation or restriction. On the left are the
associated optimal values. The equal sign denotes strong duality, i.e.,
equality of optimal values, and the arrows denote convergence of the
corresponding sequences of optimal values. Note that the question of the
closing the gap between the values of (LP${}_{n}^{L^{2}}$) and
(LP${}_{n}^{\mathcal{S}}$) is left open.

For clarity, we treat the case where the witnesses are given by the fidelity
with a single Fock state, corresponding to the case where one entry of the
vector $\bm{a}$ is equal to $1$ and all the other entries are $0$. The
generalisation to linear combinations of fidelities with Fock states is
straightforward by linearity.

A linear program is an optimisation problem where variables are linearly
constrained [61]. This can be expressed as:

$\underset{\bm{x}\in K}{\text{sup}}\;f(\bm{x}),$ (21)

for some linear real-valued function $f$ and some set $K\subset E$ where $E$
is the optimisation space which has the structure of a locally convex
topological vector space and $K$ is specified by a set of linear constraints.

An element $\bm{x}\in E$ which belongs to $K$, i.e., which satisfies all the
linear constraints, is called a feasible plan or feasible solution, and
strictly feasible solution when it strictly satisfies the constraints. The set
of feasible plans is called the feasible set. If the supremum in Eq. (21) is
attained, a feasible plan that reaches it is called an optimal plan or optimal
solution.

###  5.1 Sketch of the proofs

This section aims at giving an overview of the rather technical proofs that
follow. Our reasoning can be split into the following steps:

#### Obtaining an infinite linear program and its dual

  * •

We exploit the phase-space rotational invariance of Fock states in Lemma 2 to
express the computation of threshold values Eq. (8) on states that are
diagonal in the Fock basis.

  * •

This provides a maximisation problem which can be rephrased as an infinite-
dimensional linear program (LP${}_{n}^{L^{2}}$) (resp.
(LP${}_{n}^{\mathcal{S}}$)) on the space of square-integrable functions
$L^{2}(\mathbb{R}_{+})$ (resp. Schwartz functions
$\mathcal{S}(\mathbb{R}_{+})$). We denote $\omega_{n}^{L^{2}}$ its optimal
value (resp. $\omega_{n}^{\mathcal{S}}$).

  * •

By duality, we recast this as an optimisation problem on measures which is
given by the dual infinite-dimensional linear program (D-LP${}_{n}^{L^{2}}$)
(resp. (D-LP${}_{n}^{\mathcal{S}}$)).

#### Hierarchy of relaxations

See upper part of figure 4.

  * •

We relax the program (LP${}_{n}^{L^{2}}$): instead of asking for a positive
function, we require that this function has a positive inner product with
positive polynomials of fixed degree $m$. This degree fixes a level within a
hierarchy of relaxations. These constraints can be cast as a positive
semidefinite constraint and we thus obtain a hierarchy of semidefinite
programs (see $(\text{SDP}^{m,\geq}_{n})$).

  * •

For each level $m$, we derive the dual program $(\text{D-SDP}^{m,\geq}_{n})$
and show that strong duality holds in Theorem 3 by finding a strictly feasible
solution of $(\text{SDP}^{m,\geq}_{n})$.

  * •

We prove the convergence of the hierarchy of semidefinite relaxations towards
(D-LP${}_{n}^{L^{2}}$) in Theorem 4 by first showing that the feasible set of
$(\text{SDP}^{m,\geq}_{n})$ is compact. Then, we perform a diagonal extraction
on a sequence of optimal solutions of $(\text{SDP}^{m,\geq}_{n})$ and we
finally show that this provides a feasible solution of (LP${}_{n}^{L^{2}}$)
proving that
$\underset{m\rightarrow+\infty}{\lim}\omega^{m,\geq}_{n}=\omega_{n}^{L_{2}}$.

#### Convergence of hierarchy of restrictions

See bottom part of figure 4.

  * •

We restrict the program (LP${}_{n}^{L^{2}}$) (or equivalently the program
(LP${}_{n}^{\mathcal{S}}$)): instead of optimising over positive functions, we
optimise over positive polynomials of fixed degree $m$. Again, this degree
fixes a level within a hierarchy of restrictions. Using the fact that
univariate positive polynomials are sum-of-squares which can be written as a
semidefiniteness constraint, we obtain a hierarchy of semidefinite programs
(see $(\text{SDP}^{m,\leq}_{n})$).

  * •

For each level $m$, we derive the dual program $(\text{D-SDP}^{m,\leq}_{n})$
and show that strong duality holds in Theorem 2 by finding a strictly feasible
solution of $(\text{SDP}^{m,\leq}_{n})$.

  * •

We prove the convergence of the hierarchy towards (D-LP${}_{n}^{\mathcal{S}}$)
in Theorem 5 by first showing that the feasible set of
$(\text{D-SDP}^{m,\leq}_{n})$ is compact. This is highly nontrivial since it
requires exhibiting an analytical feasible solution of
$(\text{SDP}^{m,\leq}_{n})$ which is difficult in general and in the present
case. Then, we perform a diagonal extraction on a sequence of optimal
solutions of $(\text{D-SDP}^{m,\leq}_{n})$. A technicality arises as it does
not necessarily provide a feasible solution of (D-LP${}_{n}^{L^{2}}$) which is
why we introduce the linear program expressed over Schwartz functions. We then
show that a diagonal extraction indeed provides a feasible solution of
(D-LP${}_{n}^{\mathcal{S}}$) proving that
$\underset{m\rightarrow+\infty}{\lim}\omega^{m,\leq}_{n}=\omega_{n}^{S}$.

###  5.2 Function spaces

We now review some function spaces which appear in the following sections,
together with a few notations.

Hereafter, the half line of non-negative real numbers is denoted
$\mathbb{R}_{+}$. The space of real square-integrable functions over
$\mathbb{R}_{+}$ is denoted $L^{2}(\mathbb{R}_{+})$ and is equipped with the
usual scalar product:

$\braket{f,g}=\int_{\mathbb{R}_{+}}{f(x)g(x)dx},$ (22)

for $f,g\in L^{2}(\mathbb{R}_{+})$. This space is isomorphic to the space of
square-summable real sequences indexed by $\mathbb{N}$, $l^{2}(\mathbb{N})$,
by considering the expansion in a countable basis. Such a basis is given,
e.g., by the Laguerre functions [62], modified here by a $(-1)^{k}$ prefactor
to correspond to Fock state Wigner functions:

$\mathcal{L}_{k}(x):=(-1)^{k}L_{k}(x)e^{-\frac{x}{2}},$ (23)

for all $k\in\mathbb{N}$ and all $x\in\mathbb{R}_{+}$, where
$L_{k}(x)=\sum_{l=0}^{k}\frac{(-1)^{l}}{l!}\binom{k}{l}x^{l}$ is the
$k^{\text{th}}$ Laguerre polynomial. These functions form an orthonormal
basis: for all $p,q\in\mathbb{N}$,
$\braket{\mathcal{L}_{p},\mathcal{L}_{q}}=\delta_{pq}$, where $\delta_{pq}$ is
the Kronecker symbol.

The space $L^{2}(\mathbb{R}_{+})$ is also isomorphic to its dual space
${L^{2}}^{\prime}(\mathbb{R}_{+})$ via the Radon–Nikodym theorem [63]:
elements of ${L^{2}}^{\prime}(\mathbb{R}_{+})$ can be identified by the
Lebesgue measure on $\mathbb{R}_{+}$ times the corresponding function in
$L^{2}(\mathbb{R}_{+})$.

We write $\mathcal{S}(\mathbb{R}_{+})$ the space of Schwartz functions over
$\mathbb{R}_{+}$, i.e., the space of $C^{\infty}$ functions that go to $0$ at
infinity faster than any inverse polynomial, as do their derivatives.
$\mathcal{S}^{\prime}(\mathbb{R}_{+})$ is its dual space, the space of
tempered distributions over $\mathbb{R}_{+}$. Note that
$\mathcal{S}(\mathbb{R}_{+})$ is dense in ${L^{2}}(\mathbb{R}_{+})$. We denote
the space of rapidly decreasing real sequences by $\mathcal{S}(\mathbb{N})$
(sequences that go to $0$ at infinity faster than any inverse polynomial),
together with its dual space of slowly increasing real sequences
$\mathcal{S}^{\prime}(\mathbb{N})$ (sequences that are upper bounded by a
polynomial). The spaces $\mathcal{S}(\mathbb{R}_{+})$ and
$\mathcal{S}(\mathbb{N})$ are isomorphic: any Schwartz function over
$\mathbb{R}_{+}$ can be expanded uniquely in the basis of Laguerre functions
with a rapidly decreasing sequence of coefficients. Similarly, the spaces
$\mathcal{S}^{\prime}(\mathbb{R}_{+})$ and $\mathcal{S}^{\prime}(\mathbb{N})$
are also isomorphic: any tempered distribution over $\mathbb{R}_{+}$ can be
written uniquely as a formal series of Laguerre functions with a slowly
increasing sequence of coefficients [64]. We extend the definition of the
duality $\braket{\text{--},\text{--}}$ in Eq. (22) to these spaces.

In order to denote non-negative elements of these spaces, we will use the
notations $L^{2}_{+}(\mathbb{R}_{+})$, ${L_{+}^{2}}^{\prime}(\mathbb{R}_{+})$,
$\mathcal{S}_{+}(\mathbb{R}_{+})$ and
$\mathcal{S}_{+}^{\prime}(\mathbb{R}_{+})$. A distribution $\mu$ in
${L_{+}^{2}}^{\prime}(\mathbb{R}_{+})$ (resp. in
$\mathcal{S}_{+}^{\prime}(\mathbb{R}_{+})$) satisfies: $\forall f\in
L^{2}_{+}(\mathbb{R}_{+})$ (resp. $\forall
f\in\mathcal{S}_{+}(\mathbb{R}_{+})$), $\braket{\mu,f}\geq 0$.

For all $m\in\mathbb{N}$, we define the following space of truncated series of
Laguerre functions over $\mathbb{R}_{+}$:

$\mathcal{R}_{m}(\mathbb{R}_{+}):=\text{span}_{\mathbb{R}}\\{\mathcal{L}_{k}\\}_{0\leq
k\leq m},$ (24)

which is equal to the set of real polynomials over $\mathbb{R}_{+}$ of degree
less or equal to $m$ multiplied by the function $x\mapsto e^{-\frac{x}{2}}$.
We denote by $\mathcal{R}_{m,+}(\mathbb{R}_{+})$ its subset of non-negative
elements.

For all $\bm{s}\in\mathbb{R}^{\mathbb{N}}$, we define the associated formal
series of Laguerre functions:

$f_{\bm{s}}:=\sum_{k\geq 0}s_{k}\mathcal{L}_{k},$ (25)

with the (formal) relation:

$s_{k}=\braket{f_{\bm{s}},\mathcal{L}_{k}},$ (26)

for all $k\in\mathbb{N}$. We refer to $\bm{s}$ as the sequence of Laguerre
moments of $f_{\bm{s}}$. We extend this definition to finite sequences by
completing these sequences with zeros. For $m\in\mathbb{N}$, we also define
the matrix $A_{\bm{s}}$ (thus omitting the dependence in $m$) by

$(A_{\bm{s}})_{0\leq i,j\leq
m}:=\begin{cases}\sum\limits_{k=0}^{l}s_{k}\binom{l}{k}l!&\text{if }i+j=2l,\\\
0&\text{otherwise.}\end{cases}$ (27)

$A_{\bm{s}}$ can be seen as the Laguerre moment matrix of the measure
$f_{\bm{s}}$. In what follows, we use standard techniques relating to the
Stieltjes moment problem [65], which seeks conditions for a real sequence
$\bm{\nu}=(\nu_{k})_{k\in\mathbb{N}}\in\mathbb{R}^{\mathbb{N}}$ to be the
sequence of moments $\int_{\mathbb{R}_{+}}x^{k}d\nu(x)$ of a non-negative
distribution $\nu$ over $\mathbb{R}_{+}$. We adapt these techniques to the
basis of Laguerre functions, rather than the canonical basis. In particular,
we make use of the following result:

######  Theorem 1.

Let $\bm{\mu}=(\mu_{k})_{k\in\mathbb{N}}\in\mathbb{R}^{\mathbb{N}}$. The
sequence $\bm{\mu}$ is the sequence of Laguerre moments
$\int_{\mathbb{R}_{+}}\mathcal{L}_{k}(x)d\mu(x)$ of a non-negative
distribution $\mu$ supported on $\mathbb{R}_{+}$ if and only if

$\forall m\in\mathbb{N},\forall
g\in\mathcal{R}_{m,+}(\mathbb{R}_{+}),\;\braket{f_{\bm{\mu}},g}\geq 0.$ (28)

We give a proof in Appendix A for this result, which is based on the classic
Riesz–Haviland theorem [66, 67].

###  5.3 Linear program

In this section, we phrase the computation of the witness threshold value
introduced in Eq. (8) as an infinite-dimensional linear program, in the case
where one entry of the vector $\bm{a}$ is equal to $1$ and all the other
entries are $0$, the generalisation being straightforward by linearity.

Formally, we fix hereafter $n\in\mathbb{N}^{*}$ and we look for the witnesses
threshold value $\omega_{n}$ 333Here we write generically the threshold value
as $\omega_{n}$ while we use the more precise notation $\omega_{n}^{L^{2}}$
(resp. $\omega_{n}^{\mathcal{S}}$) to refer to its computation in the space of
square-integrable functions (resp. Schwartz functions). defined as

$\omega_{n}:=\sup_{\begin{subarray}{c}\rho\in\mathcal{D}(\mathcal{H})\\\
W_{\rho}\geq 0\end{subarray}}\braket{n}{\rho}{n}.$ (29)

This is the maximal values such that for all states
$\rho\in\mathcal{D}(\mathcal{H})$:

$\braket{n}{\rho}{n}>\omega_{n}\quad\Rightarrow\quad\exists\alpha\in\mathbb{C},\;W_{\rho}(\alpha)<0.$
(30)

Let $\mathcal{C}(\mathcal{H})$ be the set of states that are invariant under
phase-space rotations:

$\displaystyle\mathcal{C}(\mathcal{H})\\!:=\\!\\{\sigma$
$\displaystyle\in\\!\mathcal{D}(\mathcal{H})\text{ such that:}$ (31)
$\displaystyle\forall\varphi\in[0,2\pi],{\mathrm{e}}^{{\mathrm{i}}\varphi\hat{n}}\sigma{\mathrm{e}}^{-{\mathrm{i}}\varphi\hat{n}}\\!=\\!\sigma\\},$

where $\hat{n}=\hat{a}^{\dagger}\hat{a}$ is the number operator. The witnesses
corresponding to the fidelity with a single Fock state feature a rotational
symmetry in phase space, which we exploit in the following lemma.

######  Lemma 2.

The threshold value in Eq. (29) can be expressed as

$\omega_{n}=\sup_{\begin{subarray}{c}\sigma\in\mathcal{C}(\mathcal{H})\\\
W_{\sigma}\geq 0\end{subarray}}\braket{n}{\sigma}{n}.$ (32)

###### Proof.

Let $\rho\in\mathcal{D}(\mathcal{H})$. We start by applying a random dephasing
to the state $\rho$:

$\sigma=\int_{0}^{2\pi}\frac{\text{d}\varphi}{2\pi}{\mathrm{e}}^{{\mathrm{i}}\varphi\hat{n}}\rho{\mathrm{e}}^{-{\mathrm{i}}\varphi\hat{n}}\in\mathcal{C}(\mathcal{H}).$
(33)

The random dephasing does not change the fidelity with any Fock state because
of the rotational symmetry in phase space of the latter, i.e.,

$\forall n\in\mathbb{N},\,\braket{n}{\sigma}{n}=\braket{n}{\rho}{n}.$ (34)

Moreover, it can only decrease the maximum negativity of the Wigner function:
for all $\alpha\in\mathbb{C}$,

$\displaystyle W_{\sigma}(\alpha)$
$\displaystyle=\int_{0}^{2\pi}\frac{\text{d}\varphi}{2\pi}W_{e^{{\mathrm{i}}\varphi\hat{n}}\rho
e^{-{\mathrm{i}}\varphi\hat{n}}}(\alpha)$ (35)
$\displaystyle=\int_{0}^{2\pi}\frac{\text{d}\varphi}{2\pi}W_{\rho}(\alpha{\mathrm{e}}^{{\mathrm{i}}\varphi})$
$\displaystyle\geq\min_{\varphi\in[0,2\pi]}W_{\rho}(\alpha{\mathrm{e}}^{{\mathrm{i}}\varphi})$
$\displaystyle\geq\min_{\beta\in\mathbb{C}}W_{\rho}(\beta),$

and taking the minimum over all $\alpha\in\mathbb{C}$ then gives

$\min_{\alpha\in\mathbb{C}}W_{\sigma}(\alpha)\geq\min_{\beta\in\mathbb{C}}W_{\rho}(\beta).$
(36)

In particular, applying a random dephasing to a Wigner positive state yields a
Wigner positive mixtures of Fock states, which is invariant under phase-space
rotations. Hence, we can restrict without loss of generality to states that
are invariant under phase-space rotations when looking for the maximum
fidelity of Wigner positive states with a given Fock state $\ket{n}$. ∎

Lemma 2 ensures that the supremum in Eq. (29) can be computed over states that
have a rotational symmetry in phase space. Such states $\sigma$ can be
expanded diagonally in the Fock basis:

$\sigma=\sum_{k=0}^{\infty}F_{k}\ket{k}\\!\bra{k},$ (37)

with $\sum_{k}F_{k}=1$ and $0\leq F_{k}\leq 1$ for all $k\in\mathbb{N}$. By
linearity of the Wigner function:

$\forall\alpha\in\mathbb{C},\;W_{\sigma}(\alpha)=\sum_{k}F_{k}W_{k}(\alpha),$
(38)

where $W_{k}$ is the Wigner function of the $k^{th}$ Fock state [20]:

$\forall\alpha\in\mathbb{C},\,W_{k}(\alpha)=\frac{2}{\pi}\mathcal{L}_{k}(4|\alpha|^{2}),$
(39)

with $\mathcal{L}_{k}$ the $k^{\text{th}}$ Laguerre function, defined in Eq.
(23). As noted before, Fock states are invariant under phase-space rotations:
their Wigner function only depends on the amplitude of the phase-space point
considered. We fix $x=4|\alpha|^{2}\in\mathbb{R}^{+}$ hereafter.

As we consider the optimisation over $L^{2}(\mathbb{R}_{+})$ functions, we
will denote the corresponding threshold value by $\omega_{n}^{L^{2}}$. With
Lemma 2, the computation of $\omega_{n}^{L^{2}}$ can thus be expressed as the
following infinite-dimensional linear program:

$\displaystyle\quad\text{Find
}(F_{k})_{k\in\mathbb{N}}\in\ell^{2}(\mathbb{N})$ (LP${}_{n}^{L^{2}}$)
$\displaystyle\quad\text{maximising }F_{n}$ subject to
$\displaystyle\begin{dcases}\sum_{k}F_{k}=1\\\ \forall
k\in\mathbb{N},\;F_{k}\geq 0\rule[-9.60004pt]{0.0pt}{9.60004pt}\\\ \forall
x\in\mathbb{R}_{+},\;\sum_{k}F_{k}\mathcal{L}_{k}(x)\geq 0.\end{dcases}$

The first constraint ensures unit trace of the corresponding state $\sigma$,
the second one ensures that its fidelity with each Fock state is non-negative,
and the last one ensures that its Wigner function $W_{\sigma}$ is non-
negative. Note that $\omega_{n}>0$ for all $n\in\mathbb{N}^{*}$, by
considering a mixture of $\ket{0}$ and $\ket{n}$ with the vacuum component
close enough to $1$. Note that, for consistency of notations, we use $L$
instead of $\ell$ in the label of the program.

We refer the interested reader to Appendix B where this program is expressed
in the canonical form of infinite-dimensional linear programs [61] and its
dual is derived. This dual linear program reads:

$\displaystyle\quad\text{Find }y\in\mathbb{R}\text{ and
}\mu\in{L^{2}}^{\prime}(\mathbb{R}_{+})$ (D-LP${}_{n}^{L^{2}}$)
$\displaystyle\quad\text{minimising }y$ subject to
$\displaystyle\begin{dcases}\forall k\neq
n\in\mathbb{N},\;y\geq\int_{\mathbb{R}_{+}}{\mathcal{L}_{k}}{d\mu}\\\ y\geq
1+\int_{\mathbb{R}_{+}}{\mathcal{L}_{n}}{d\mu}\rule[-9.60004pt]{0.0pt}{9.60004pt}\\\
\forall f\in L^{2}_{+}(\mathbb{R}_{+}),\;\langle\mu,f\rangle\geq
0.\end{dcases}$

We also consider a restriction of (LP${}_{n}^{L^{2}}$) by optimising over
Schwartz functions rather than square-integrable functions as
$\mathcal{S}(\mathbb{R}_{+})\subset L^{2}(\mathbb{R}_{+})$. The primal program
becomes:

$\displaystyle\quad\text{Find
}(F_{k})_{k\in\mathbb{N}}\in\mathcal{S}(\mathbb{N})$
(LP${}_{n}^{\mathcal{S}}$) $\displaystyle\quad\text{maximising }F_{n}$ subject
to $\displaystyle\begin{dcases}\sum_{k}F_{k}=1\\\ \forall
k\in\mathbb{N},\;F_{k}\geq 0\rule[-9.60004pt]{0.0pt}{9.60004pt}\\\ \forall
x\in\mathbb{R}_{+},\;\sum_{k}F_{k}\mathcal{L}_{k}(x)\geq 0,\end{dcases}$

We denote its value by $\omega_{n}^{\mathcal{S}}$. Because it is a restriction
of (LP${}_{n}^{L^{2}}$), we have that
$\omega_{n}^{L^{2}}\geq\omega_{n}^{\mathcal{S}}$. Its dual linear program can
be expressed as:

$\displaystyle\quad\text{Find }y\in\mathbb{R}\text{ and
}\mu\in\mathcal{S}^{\prime}(\mathbb{R}_{+})$ (D-LP${}_{n}^{\mathcal{S}}$)
$\displaystyle\quad\text{minimising }y$ subject to
$\displaystyle\begin{dcases}\forall k\neq
n\in\mathbb{N},\;y\geq\int_{\mathbb{R}_{+}}{\mathcal{L}_{k}}{d\mu}\\\ y\geq
1+\int_{\mathbb{R}_{+}}{\mathcal{L}_{n}}{d\mu}\rule[-9.60004pt]{0.0pt}{9.60004pt}\\\
\forall f\in\mathcal{S}_{+}(\mathbb{R}_{+}),\;\langle\mu,f\rangle\geq
0.\end{dcases}$

We will show later that strong duality holds for the two pairs of programs.
However weak duality of linear programming already ensures that the optimal
value $\omega_{n}^{L^{2}}$ of (LP${}_{n}^{L^{2}}$) (resp.
$\omega_{n}^{\mathcal{S}}$ of (LP${}_{n}^{\mathcal{S}}$)) is upper bounded by
the optimal value of (D-LP${}_{n}^{L^{2}}$) (resp.
(D-LP${}_{n}^{\mathcal{S}}$)). Hence, a possible way of solving the
optimisation (LP${}_{n}^{L^{2}}$) is to exhibit a feasible solution for
(LP${}_{n}^{L^{2}}$) and a feasible solution for (D-LP${}_{n}^{L^{2}}$) that
have the same value.

For $n=1$, choosing
$(F_{k})_{k\in\mathbb{N}}=(\frac{1}{2},\frac{1}{2},0,0,\dots)$ gives a
feasible solution for (LP${}_{1}^{L^{2}}$) (resp. (LP${}_{1}^{\mathcal{S}}$))
with the value $\frac{1}{2}$, while choosing
$(y,\mu)=(\frac{1}{2},\frac{1}{2}\delta(x))$, where $\delta$ is the Dirac
delta function444Technically $\delta\notin{L^{2}}^{\prime}(\mathbb{R}_{+})$,
but the result holds by considering a sequence of functions converging to a
Dirac delta. over $\mathbb{R}_{+}$, gives a feasible solution for
(D-LP${}_{1}^{L^{2}}$) (resp. (D-LP${}_{1}^{\mathcal{S}}$)) with the value
$\frac{1}{2}$. This shows that
$\omega_{1}^{L^{S}}=\omega_{1}^{\mathcal{S}}=\frac{1}{2}$.

Similarly, for $n=2$, choosing
$(F_{k})_{k\in\mathbb{N}}=(\frac{1}{2},0,\frac{1}{2},0,0,\dots)$ gives a
feasible solution for (LP${}_{2}^{L^{2}}$) (resp. LP${}_{2}^{\mathcal{S}}$)
with the value $\frac{1}{2}$, while choosing
$(y,\mu)=(\frac{1}{2},\frac{e}{2}\delta(x-2))$ gives a feasible solution for
(D-LP2) (resp. (D-LP${}_{2}^{\mathcal{S}}$)) , up to a conjecture555We checked
numerically the corresponding constraints $|L_{k}(2)|\leq 1$ for $k$ up to
$10^{3}$ and, considering asymptotic behaviors, we conjecture that these hold
for all $k\geq 0$., with the value $\frac{1}{2}$. This shows that
$\omega_{2}^{L^{2}}=\omega_{2}^{S}=\frac{1}{2}$.

While this approach is sensible for small values of $n$, finding optimal
analytical solutions for higher values of $n$ seems highly nontrivial.
Moreover, the infinite number of variables prevents us from performing the
optimisation (LP${}_{n}^{L^{2}}$) numerically. A natural workaround is to find
finite-dimensional relaxations or restrictions of the original problem—thus
providing upper and lower bounds for the optimal value $\omega_{n}^{L^{2}}$,
respectively. This is the approach we follow in the next section.

###  5.4 Hierarchies of semidefinite programs

Semidefinite programming is a convex optimisation technique in the cone of
positive semidefinite matrices. It comes with a duality theory: if a primal
semidefinite program is a maximisation problem, then one may deduce a dual
minimisation problem which is again a semidefinite program. Like linear
programs, these programs satisfy a weak duality condition: the optimal value
of the primal problem is upper bounded by the optimal value of the dual
program. The difference between the optimal values is called the duality gap.
When there is no duality gap, we say that there is strong duality between the
programs.

####  5.4.1 Preliminaries

In this section, we introduce preliminary technical lemmas and we refer the
reader to Appendix C for their proofs.

We recall the following standard result, which comes from the fact that any
univariate polynomial non-negative over $\mathbb{R}$ can be written as a sum
of squares:

######  Lemma 3 ([68]).

Let $p\in\mathbb{N}$ and let $P$ be a univariate polynomial of degree $2p$.
Let $X=(1,x,\dots,x^{p})$ be the vector of monomials. Then, $P$ is non-
negative over $\mathbb{R}$ if and only if there exists a real
$(p+1)\times(p+1)$ positive semidefinite matrix $Q$ such that for all
$x\in\mathbb{R}$,

$P(x)=X^{T}QX.$ (40)

From this lemma we deduce the following characterisation of non-negative
polynomials over $\mathbb{R}_{+}$:

######  Lemma 4.

Non-negative polynomials on $\mathbb{R}_{+}$ can be written as sums of
polynomials of the form $\sum_{l=0}^{p}x^{l}\sum_{i+j=2l}y_{i}y_{j}$, where
$p\in\mathbb{N}$ and $y_{i}\in\mathbb{R}$, for all $0\leq i\leq p$.

Note that the characterisation in Lemma 4 differs from that of Stieltjes [65],
which expresses non-negative polynomials over $\mathbb{R}_{+}$ as $x\mapsto
A_{1}(x)+xA_{2}(x)$, where $A_{1}$ and $A_{2}$ are sums of squares. This
slightly slows down the numerical resolution, which does not matter given the
size of the programs considered. At the same time, this allows us to obtain
more compact expressions for the semidefinite programs.

We use the characterisation of Lemma 4 to obtain the following result: for
$\bm{s}\in\mathbb{R}^{\mathbb{N}}$, the fact that the series $f_{\bm{s}}$
defined in Eq. (25) has non-negative scalar product with non-negative
truncated Laguerre series up to degree $m$ can be expressed as a positive
semidefinite constraint involving the matrix $A_{\bm{s}}$ defined in Eq. (27).
Formally:

######  Lemma 5.

Let $m\geq n$ and let $\bm{s}\in\mathbb{R}^{\mathbb{N}}$. The following
propositions are equivalent:

  1. (i)

$\forall g\in\mathcal{R}_{m,+}(\mathbb{R}_{+}),\;\braket{f_{\bm{s}},g}\geq 0$,

  2. (ii)

$A_{\bm{s}}\succeq 0$.

Using these results, we derive hierarchies of semidefinite relaxations and
restrictions for the infinite-dimensional linear program (LP${}_{n}^{L^{2}}$)
in the following sections.

####  5.4.2 Semidefinite relaxations

One way to obtain a relaxation of (LP${}_{n}^{L^{2}}$) is to relax the
constraint:

$\forall
x\in\mathbb{R}_{+},\;f_{\bm{F}}(x)=\sum_{k}F_{k}\mathcal{L}_{k}(x)\geq 0.$
(41)

Instead, one may impose the weaker constraint:

$\forall g\in\mathcal{R}_{m,+}(\mathbb{R}_{+}),\;\braket{f_{\bm{F}},g}\geq 0,$
(42)

for some fixed $m\geq n$. By Lemma 5, this constraint may in turn be expressed
as a positive semidefinite constraint on the $(m+1)\times(m+1)$ matrix
$A_{\bm{F}}$ defined in Eq. (27).

Each choice of $m$ thus leads to a different semidefinite program, whose
optimal value gets closer to $\omega_{n}$ as $m$ increases (since the
constraint (42) gets stronger when $m$ increases). Moreover, when replacing
the constraint (41) by the constraint (42) in (LP${}_{n}^{L^{2}}$), for $l>m$
the variables $F_{l}$ do not appear in the matrix $A_{\bm{F}}$, and are
constrained only by $\sum_{k}F_{k}=1$ and $F_{l}\geq 0$. Since $m\geq n$ and
the optimal value is $F_{n}$, we may thus set $F_{l}=0$ for $l>m$ without loss
of generality. This gives a hierarchy of finite-dimensional semidefinite
relaxations for (LP${}_{n}^{L^{2}}$), and the semidefinite relaxation of order
$m$ is given by:

$\displaystyle\quad\text{Find }A\in\text{Sym}_{m+1}\text{ and
}\bm{F}\in\mathbb{R}^{m+1}$ $(\text{SDP}^{m,\geq}_{n})$
$\displaystyle\quad\text{maximising }F_{n}$ subject to
$\displaystyle\begin{dcases}\textstyle\sum_{k=0}^{m}F_{k}=1\\\
\textstyle\forall k\leq m,\;F_{k}\geq 0\\\ \textstyle\forall l\leq m,\forall
i+j=2l,\;A_{ij}=\sum\limits_{k=0}^{l}F_{k}\binom{l}{k}l!\\\ \textstyle\forall
l\in\llbracket 1,m\rrbracket,\forall i+j=2l-1,\;A_{ij}=0\\\ \textstyle
A\succeq 0.\end{dcases}$

Let us denote its optimal value by $\omega^{m,\geq}_{n}$. The sequence
$\\{\omega^{m,\geq}_{n}\\}_{m\geq n}$ is a decreasing sequence and for all
$m\geq n$, we have $\omega_{n}^{L^{2}}\leq\omega^{m,\geq}_{n}$.

For each $m\geq n$, the program $(\text{SDP}^{m,\geq}_{n})$ has a dual
semidefinite program which is given by (see Appendix D for a detailed
derivation):

$\displaystyle\quad\text{Find
}Q\in\text{Sym}_{m+1},\bm{\mu}\in\mathbb{R}^{m+1}\text{ and }y\in\mathbb{R}$
$(\text{D-SDP}^{m,\geq}_{n})$ $\displaystyle\quad\text{minimising }y$ subject
to $\displaystyle\begin{dcases}\textstyle y\geq 1+\mu_{n}\\\ \textstyle\forall
k\leq m,\;y\geq\mu_{k}\\\ \textstyle\forall l\leq
m,\sum\limits_{i+j=2l}Q_{ij}=\textstyle\sum\limits_{k=l}^{m}\frac{(-1)^{k+l}}{l!}\binom{k}{l}\mu_{k}\\\
\textstyle Q\succeq 0.\end{dcases}$

We show in Theorem 3 that strong duality holds between the primal and the dual
versions of this semidefinite program. In particular, numerical computations
with either of these programs will yield the same optimal value.

####  5.4.3 Semidefinite restrictions

A trivial way to obtain a restriction of (LP${}_{n}^{L^{2}}$) (or equivalently
of (LP${}_{n}^{\mathcal{S}}$)) is to impose $F_{l}=0$ for $l>m$, for some
$m\geq n$. What is less trivial is that this yields a finite-dimensional
semidefinite program. Indeed, the constraint (41) becomes

$\forall x\in\mathbb{R}_{+},\;\sum_{k=0}^{m}F_{k}\mathcal{L}_{k}(x)\geq 0,$
(43)

or equivalently:

$\forall x\in\mathbb{R},\;\sum_{k=0}^{m}(-1)^{k}F_{k}L_{k}(x^{2})\geq 0,$ (44)

where we used Eq. (23). By Lemma 3, writing $X=(1,x,\dots,x^{m})$, this is
equivalent to the existence of a positive semidefinite matrix
$Q=(Q_{ij})_{0\leq i,j\leq m}$ such that for all $x\in\mathbb{R}$,

$\displaystyle\sum_{k=0}^{m}(-1)^{k}F_{k}L_{k}(x^{2})$ $\displaystyle=X^{T}QX$
(45) $\displaystyle=\sum_{l=0}^{m}x^{l}\sum_{i+j=l}Q_{ij}.$

This is in turn equivalent to the linear constraints:

$\displaystyle\forall l\in\llbracket
1,m\rrbracket,\;\sum\limits_{i+j=2l-1}Q_{ij}=0,$ (46) and
$\displaystyle\forall l\leq
m,\sum\limits_{i+j=2l}Q_{ij}=\frac{(-1)^{l}}{l!}\sum\limits_{k=l}^{m}(-1)^{k}\binom{k}{l}F_{k},$

by identifying the coefficients in front of each monomial in Eq. (45).

Hence, the restriction of (LP${}_{n}^{L^{2}}$) obtained by imposing $F_{l}=0$
for $l>m$, for a fixed $m\geq n$, is a semidefinite program given by:

$\displaystyle\quad\text{Find }Q\in\text{Sym}_{m+1}\text{ and
}\bm{F}\in\mathbb{R}^{m+1}$ $(\text{SDP}^{m,\leq}_{n})$
$\displaystyle\quad\text{maximising }F_{n}$ subject to
$\displaystyle\begin{dcases}\textstyle\sum_{k=0}^{m}F_{k}=1\\\
\textstyle\forall k\leq m,\;F_{k}\geq 0\\\ \textstyle\forall l\in\llbracket
1,m\rrbracket,\;\sum\limits_{i+j=2l-1}Q_{ij}=0\\\ \textstyle\forall l\leq
m,\sum\limits_{i+j=2l}Q_{ij}=\sum\limits_{k=l}^{m}\frac{(-1)^{k+l}}{l!}\binom{k}{l}F_{k}\\\
\textstyle Q\succeq 0.\end{dcases}$

Let us denote its optimal value by $\omega^{m,\leq}_{n}$. Each choice of $m$
leads to a different semidefinite restriction of (LP${}_{n}^{L^{2}}$) and
(LP${}_{n}^{\mathcal{S}}$), whose optimal value gets closer to
$\omega_{n}^{\mathcal{S}}$ as $m$ increases (since the constraint $F_{l}=0$
for $l>m$ gets weaker when $m$ increases). The sequence
$\\{\omega^{m,\leq}_{n}\\}_{m\geq n}$ is thus an increasing sequence and for
all $m\geq n$, we have
$\omega^{m,\leq}_{n}\leq\omega_{n}^{\mathcal{S}}\leq\omega_{n}^{L^{2}}$.

For each $m\geq n$, the program $(\text{SDP}^{m,\leq}_{n})$ has a dual
semidefinite program which is given by (see Appendix D for a detailed
derivation):

$\displaystyle\quad\text{Find
}A\in\text{Sym}_{m+1},\bm{\mu}\in\mathbb{R}^{m+1}\text{ and }y\in\mathbb{R}$
$(\text{D-SDP}^{m,\leq}_{n})$ $\displaystyle\quad\text{minimising }y$ subject
to $\displaystyle\begin{dcases}\textstyle y\geq 1+\mu_{n}\\\ \textstyle\forall
k\leq m,\;y\geq\mu_{k}\\\ \textstyle\forall l\leq m,\forall
i+j=2l,\;A_{ij}=\sum\limits_{k=0}^{l}\mu_{k}\binom{l}{k}l!\\\ \textstyle
A\succeq 0.\end{dcases}$

We show in Theorem 2 that strong duality holds between the primal and the dual
versions of this semidefinite program.

####  5.4.4 Strong duality of semidefinite programs

We fix $m\geq n$ and we show that strong duality holds both for the
semidefinite restriction $(\text{SDP}^{m,\leq}_{n})$ and the semidefinite
relaxation $(\text{SDP}^{m,\geq}_{n})$.

We first consider the semidefinite restrictions:

######  Theorem 2.

Strong duality holds between the programs $(\text{SDP}^{m,\leq}_{n})$ and
$(\text{D-SDP}^{m,\leq}_{n})$.

###### Proof.

We make use of Slater condition for (finite-dimensional) semidefinite
programs: strict feasibility of $(\text{SDP}^{m,\leq}_{n})$ implies strong
duality between $(\text{SDP}^{m,\leq}_{n})$ and $(\text{D-SDP}^{m,\leq}_{n})$.

In order to obtain a strictly feasible solution, we define
$Q:=\frac{1}{2^{m+1}-1}\text{Diag}_{k=0,\dots,m}(\frac{1}{k!})\in\text{Sym}_{m+1}$
and $\bm{F}=(F_{0},\dots,F_{m})\in\mathbb{R}^{m+1}$, where for all
$k\in\llbracket 0,m\rrbracket$, $F_{k}:=\frac{1}{2^{m+1}-1}\binom{m+1}{k+1}$.
Then $Q\succ 0$ and $F_{k}>0$ for all $k\in\llbracket 0,m\rrbracket$.
Moreover, we have

$\displaystyle\sum_{k=0}^{m}F_{k}$
$\displaystyle=\frac{1}{2^{m+1}-1}\sum_{k=0}^{m}\binom{m+1}{k+1}$ (47)
$\displaystyle=\frac{1}{2^{m+1}-1}\left(\sum_{k=0}^{m+1}\binom{m+1}{k}-1\right)$
$\displaystyle=1.$

We also have $\sum_{i+j=2l-1}Q_{ij}=0$ for all $l\in\llbracket 1,m\rrbracket$,
since $Q$ is diagonal. Furthermore, for all $l\leq m$,

$\sum_{i+j=2l}Q_{ij}=Q_{ll}=\frac{1}{2^{m+1}-1}\frac{1}{l!},$ (48)

and

$\displaystyle\frac{(-1)^{l}}{l!}\sum\limits_{k=l}^{m}(-1)^{k}\binom{k}{l}F_{k}$
(49)
$\displaystyle=\frac{1}{2^{m+1}-1}\frac{1}{l!}\sum_{k=l}^{m}(-1)^{k-l}\binom{k}{m}\binom{m+1}{k+1}$
$\displaystyle=\frac{1}{2^{m+1}-1}\frac{1}{l!}\binom{m}{l}\sum_{q=0}^{m-l}(-1)^{q}\frac{m+1}{q+l+1}\binom{m-l}{q}$
$\displaystyle=\frac{1}{2^{m+1}-1}\frac{1}{l!},$

where we used [69, (1.41)] in the last line. Therefore, $(Q,\bm{F})$ is a
strictly feasible solution of $(\text{SDP}^{m,\leq}_{n})$, which implies
strong duality. ∎

As a consequence of the proof of Theorem 2, we also obtain strong duality for
the semidefinite relaxations:

######  Theorem 3.

Strong duality holds between the programs $(\text{SDP}^{m,\geq}_{n})$ and
$(\text{D-SDP}^{m,\geq}_{n})$.

###### Proof.

The program $(\text{SDP}^{m,\geq}_{n})$ is a relaxation of
(LP${}_{n}^{L^{2}}$) and $(\text{SDP}^{m,\leq}_{n})$ is a restriction of
(LP${}_{n}^{L^{2}}$), so $(\text{SDP}^{m,\geq}_{n})$ is a relaxation of
$(\text{SDP}^{m,\leq}_{n})$. Hence, the strictly feasible solution of
$(\text{SDP}^{m,\leq}_{n})$ derived in the proof of Theorem 2 yields a
strictly feasible solution $(A,\bm{F})$ for $(\text{SDP}^{m,\geq}_{n})$: we
set $\bm{F}=(F_{0},\dots,F_{m})\in\mathbb{R}^{m+1}$, where for all
$k\in\llbracket 0,m\rrbracket$, $F_{k}:=\frac{1}{2^{m+1}-1}\binom{m+1}{k+1}$
and $A=A_{\bm{F}}\in\text{Sym}_{m+1}$, where $A_{\bm{F}}$ is defined in Eq.
(27).

With Slater condition, this shows again that strong duality holds between the
programs $(\text{SDP}^{m,\geq}_{n})$ and $(\text{D-SDP}^{m,\geq}_{n})$. ∎

###  5.5 Convergence of the hierarchies of semidefinite programs

From the previous sections, for $m\geq n$ the optimal values
$\omega^{m,\geq}_{n}$ and $\omega^{m,\leq}_{n}$ of $(\text{SDP}^{m,\geq}_{n})$
and $(\text{SDP}^{m,\leq}_{n})$ form decreasing and increasing sequences,
respectively, which satisfy

$0\leq\omega^{m,\leq}_{n}\leq\omega_{n}^{\mathcal{S}}\leq\omega_{n}^{L^{2}}\leq\omega^{m,\geq}_{n}\leq
1.$ (50)

Recall that $\omega_{n}^{L^{2}}$ is the optimal value of (LP${}_{n}^{L^{2}}$)
while $\omega_{n}^{\mathcal{S}}$ is the optimal value of
(LP${}_{n}^{\mathcal{S}}$). These sequences thus both converge, and the
remaining question is whether $(\omega^{m,\geq}_{n})_{m}$ converges to
$\omega_{n}^{L^{2}}$ and $(\omega^{m,\leq}_{n})_{m}$ converges to
$\omega_{n}^{\mathcal{S}}$. In this section, we show that this is indeed the
case.

####  5.5.1 Convergence of the sequence of upper bounds

######  Theorem 4.

The decreasing sequence of optimal values of $(\text{SDP}^{m,\geq}_{n})$
converges to the optimal value of (LP${}_{n}^{L^{2}}$):

$\lim_{m\rightarrow+\infty}\omega^{m,\geq}_{n}=\omega_{n}^{L^{2}}.$ (51)

In order to prove this theorem, we extract a limit from a sequence of optimal
solutions of $(\text{SDP}^{m,\geq}_{n})$, for $m\geq n$, and we show using
Theorem 1 that it provides a feasible solution of (LP${}_{n}^{L^{2}}$).

###### Proof.

For all $m\geq n$, the feasible set of $(\text{SDP}^{m,\geq}_{n})$ is non-
empty (consider, e.g., $\bm{F}=(1,0,0,\dots,0)\in\mathbb{R}^{m+1}$). Moreover,
due to the constraints $\sum_{k=0}^{m}F_{k}=1$ and $F_{k}\geq 0$ for all
$k\leq m$, the feasible set of $(\text{SDP}^{m,\geq}_{n})$ is compact. Hence,
the program $(\text{SDP}^{m,\geq}_{n})$ has feasible optimal solutions, for
all $m\geq n$, by diagonal extraction.

The matrix $A$ in $(\text{SDP}^{m,\geq}_{n})$ is entirely fixed by the choice
of $\bm{F}$. Let $(\bm{F}^{m})_{m\geq n}$ be a sequence of optimal solutions
of $(\text{SDP}^{m,\geq}_{n})$, for $m\geq n$. For each $m\geq n$, we have by
optimality that $F_{n}^{m}=\omega_{n}^{m,\geq}$, and the sequence
$(F^{m}_{n})_{m\geq n}$ converges. We complete each tuple
$\bm{F}^{m}=(F_{0}^{m},F_{1}^{m},\dots,F_{m}^{m})\in\mathbb{R}^{m+1}$ with
zeros to obtain a sequence in $\mathbb{R}^{\mathbb{N}}$, which we still denote
$\bm{F}^{m}=(F_{0}^{m},F_{1}^{m},\dots,F_{m}^{m},0,0,\dots)\in\mathbb{R}^{\mathbb{N}}$.

Performing a diagonal extraction $\phi$ on the sequence of optimal solutions
$(\bm{F}^{m})_{m\geq n}$, we obtain a sequence of sequences
$(\bm{F}^{\phi(m)})_{m\geq n}$ such that each sequence
$(F^{\phi(m)}_{k})_{m\geq n}$ converges when $m\rightarrow+\infty$, for all
$k\in\mathbb{N}$. Let $F_{k}$ denote its limit, for each $k\in\mathbb{N}$. We
write $\bm{F}=(F_{k})_{k\in\mathbb{N}}\in\mathbb{R}^{\mathbb{N}}$ the sequence
of limits.

For all $m\geq n$, $F_{k}^{\phi(m)}\geq 0$ for all $k\in\mathbb{N}$ and
$\sum_{k}F_{k}^{\phi(m)}=1$, so taking $m\rightarrow+\infty$ we obtain
$F_{k}\geq 0$ for all $k\in\mathbb{N}$, and $\sum_{k}F_{k}\leq 1$. Moreover,

$F_{n}=\lim_{m\rightarrow+\infty}F_{n}^{m}=\lim_{m\rightarrow+\infty}\omega^{m,\geq}_{n}.$
(52)

For all $m\geq n$, we have $\omega^{m,\geq}_{n}\geq\omega_{n}^{L^{2}}$, so
$F_{n}\geq\omega_{n}^{L^{2}}>0$. In particular, $\sum_{k}F_{k}>0$, so without
loss of generality we may assume that $\sum_{k}F_{k}=1$ (otherwise we can
always replace $F_{k}$ by $\frac{F_{k}}{\sum_{l}F_{l}}$).

Let $f_{\bm{F}}=\sum_{k}F_{k}\mathcal{L}_{k}\in L^{2}(\mathbb{R}_{+})$. By
construction we have:

$\forall m\geq n,\forall
g\in\mathcal{R}_{m,+}(\mathbb{R}_{+}),\;\braket{f_{\bm{F}},g}\geq 0.$ (53)

Hence, by Theorem 1, $\bm{F}$ is the sequence of Laguerre moments of a non-
negative distribution over $\mathbb{R}_{+}$ (the Lebesgue measure times the
function $f_{\bm{F}}$). In particular,

$\forall
x\in\mathbb{R}_{+},\;f_{\bm{F}}(x)=\sum_{k}F_{k}\mathcal{L}_{k}(x)\geq 0.$
(54)

With the constraints $F_{k}\geq 0$ for all $k\in\mathbb{N}$, and
$\sum_{k}F_{k}\leq 1$, this implies that $\bm{F}$ is a feasible solution of
(LP${}_{n}^{L^{2}}$), and in particular $F_{n}\leq\omega_{n}^{L^{2}}$, since
(LP${}_{n}^{L^{2}}$) is a maximisation problem. Since we already had
$F_{n}\geq\omega_{n}^{L^{2}}$ we obtain with Eq. (52):

$\lim_{m\rightarrow+\infty}\omega^{m,\geq}_{n}=\omega_{n}^{L^{2}},$ (55)

which concludes the proof. ∎

This immediately implies strong duality for programs (LP${}_{n}^{L^{2}}$) and
(D-LP${}_{n}^{L^{2}}$) because

  1. (i)

we have weak duality between those programs so that the optimal value
$\omega^{\prime L^{2}}_{n}$ of (D-LP${}_{n}^{L^{2}}$) upper bounds the optimal
value of (LP${}_{n}^{L^{2}}$) $\omega^{L^{2}}_{n}$ i.e.
$\omega^{L^{2}}_{n}\leq\omega^{\prime L^{2}}_{n}$;

  2. (ii)

we have strong duality between $(\text{SDP}^{m,\geq}_{n})$ and
$(\text{D-SDP}^{m,\geq}_{n})$ by Theorem 3;

  3. (iii)

$(\text{SDP}^{m,\geq}_{n})$ is a relaxation of (LP${}_{n}^{L^{2}}$) so that
$\forall m,\;\omega^{m,\geq}_{n}\geq\omega_{n}^{L^{2}}$;

  4. (iv)

$(\text{D-SDP}^{m,\geq}_{n})$ is a restriction of (D-LP${}_{n}^{L^{2}}$) so
that $\forall m,\omega^{m,\geq}_{n}\geq\omega^{\prime
L^{2}}_{n}\geq\omega^{L^{2}}_{n}$;

  5. (v)

we showed that the optimal value of the hierarchy $(\text{SDP}^{m,\geq}_{n})$
converges to $\omega^{L^{2}}_{n}$ i.e.
$\lim_{m\rightarrow+\infty}\omega^{m,\geq}_{n}=\omega_{n}^{L^{2}}$.

####  5.5.2 Convergence of the sequence of lower bounds

######  Theorem 5.

The increasing sequence of optimal values of $(\text{SDP}^{m,\leq}_{n})$
converges to the optimal value of (LP${}_{n}^{\mathcal{S}}$):

$\lim_{m\rightarrow+\infty}\omega^{m,\leq}_{n}=\omega_{n}^{\mathcal{S}}.$ (56)

The proof is similar to that of Theorem 4 using the dual programs: we attempt
to construct a feasible optimal solution of (D-LP${}_{n}^{L^{2}}$) by
extracting a limit from a sequence of optimal solutions of
$(\text{D-SDP}^{m,\leq}_{n})$, for $m\geq n$ and then conclude using the
strong duality between $(\text{SDP}^{m,\leq}_{n})$ and
$(\text{D-SDP}^{m,\leq}_{n})$, which we proved in Theorem 2. However, it turns
out that (D-LP${}_{n}^{L^{2}}$) may not have feasible optimal solutions in
${L^{2}}^{\prime}(\mathbb{R}_{+})$ (as anticipated with the analytical optimal
solutions for $n=1$ and $n=2$ from section 5.3).

To deal with this issue, we have extended the formulation of
(D-LP${}_{n}^{L^{2}}$) to a larger space where it has feasible optimal
solutions, namely the space of tempered distributions
$\mathcal{S}^{\prime}(\mathbb{R}_{+})$ (see (D-LP${}_{n}^{\mathcal{S}}$)).

Note that the semidefinite restrictions $(\text{SDP}^{m,\leq}_{n})$ of
(LP${}_{n}^{L^{2}}$) are also restrictions of (LP${}_{n}^{\mathcal{S}}$) for
all $m\geq n$, while (LP${}_{n}^{\mathcal{S}}$) is itself a restriction of
(LP${}_{n}^{L^{2}}$), and (D-LP${}_{n}^{\mathcal{S}}$) is a relaxation of
(D-LP${}_{n}^{L^{2}}$). We denote by ${\omega_{n}^{\mathcal{S}^{\prime}}}$ the
optimal value of (D-LP${}_{n}^{\mathcal{S}}$). Recall that the optimal value
of (LP${}_{n}^{\mathcal{S}}$) is denoted $\omega_{n}^{\mathcal{S}}$ and that
of (LP${}_{n}^{L^{2}}$) and (D-LP${}_{n}^{L^{2}}$) are denoted
$\omega_{n}^{L^{2}}$ (by strong duality). By weak duality of linear
programming we thus have

$\omega_{n}^{m,\leq}\leq\omega_{n}^{\mathcal{S}}\leq{\omega_{n}^{\mathcal{S}^{\prime}}}\leq\omega_{n}^{L^{2}},$
(57)

for all $m\geq n$ (see also figure 4).

Before proving Theorem 5, we introduce two intermediate technical results. The
first one is a reformulation of $(\text{D-SDP}^{m,\leq}_{n})$ over Schwartz
space:

######  Lemma 6.

For all $m\geq n$, the program $(\text{D-SDP}^{m,\leq}_{n})$ is equivalent to
the following program:

$\displaystyle\quad\text{Find }y\in\mathbb{R}\;\text{and
}\bm{\mu}\in\mathcal{S}^{\prime}(\mathbb{N})$ (D-SDP${}_{n}^{m,\leq}$)
$\displaystyle\quad\text{minimising }y$ subject to
$\displaystyle\begin{dcases}\forall k\neq n\in\mathbb{N},\;y\geq\mu_{k}\\\
y\geq 1+\mu_{n}\\\ \forall g\in\mathcal{R}_{m,+}(\mathbb{R}_{+}),\;\langle
f_{\bm{\mu}},g\rangle\geq 0,\end{dcases}$

where $f_{\bm{\mu}}=\sum_{k}\mu_{k}\mathcal{L}_{k}$.

This reformulation uses Stieltjes characterisation of non-negative polynomials
over $\mathbb{R}_{+}$ [65] and is detailed in Appendix E.

The second result provides a nontrivial analytical solution to the primal
program $(\text{SDP}^{m,\leq}_{n})$: let us define
$\bm{F}^{n}=(F_{k}^{n})_{k\in\mathbb{N}}\in\mathbb{R}^{\mathbb{N}}$ by
$\bullet$ if $n$ is even:

$F_{k}^{n}:=\begin{cases}\frac{1}{2^{n}}\binom{k}{\frac{k}{2}}\binom{n-k}{\frac{n-k}{2}}&\text{when
}k\leq n,k\text{ even},\\\ 0&\text{otherwise},\end{cases}$ (58)

$\bullet$ if $n$ is odd:

$F_{k}^{n}:=\begin{cases}\frac{1}{2^{n}}\frac{\binom{n}{\lfloor\frac{n}{2}\rfloor}\binom{\lfloor\frac{n}{2}\rfloor}{\lfloor\frac{k}{2}\rfloor}^{2}}{\binom{n}{k}},&\text{when
}k\leq n,\\\ 0&\text{otherwise}.\end{cases}$ (59)

In both cases,

$\displaystyle F_{n}^{n}$
$\displaystyle=\frac{1}{2^{n}}\binom{n}{\lfloor\frac{n}{2}\rfloor}$ (60)
$\displaystyle\geq\frac{1}{n+1}$

where we used $\binom{n}{j}\leq\binom{n}{\lfloor\frac{n}{2}\rfloor}$ for all
$j\in\llbracket 0,n\rrbracket$, summed over $j$ from 0 to $n$.

######  Lemma 7.

For all $m\geq n$, $\bm{F}^{n}$ is a feasible solution of
$(\text{SDP}^{m,\leq}_{n})$. Moreover, it is optimal when $m=n$.

The proof of feasibility consists in checking that the constraints of
$(\text{SDP}^{m,\leq}_{n})$ are satisfied by $\bm{F}^{n}$. To do so, we make
use of Zeilberger’s algorithm [70], a powerful algorithm for proving binomial
identities. The proof of optimality for $m=n$ is obtained by deriving an
analytical feasible solution of (D-SDP${}_{n}^{n,\leq}$) with the same optimal
value. Note that the optimality of this solution does not play a role in the
proof of convergence. We refer to Appendix F for a detailed proof.

As a consequence, we obtain the following analytical lower bound for the
optimal value of (LP${}_{n}^{L^{2}}$):

$\omega_{n}\geq\omega_{n}^{n,\leq}=\frac{1}{2^{n}}\binom{n}{\lfloor\frac{n}{2}\rfloor}\underset{n\rightarrow+\infty}{\sim}\sqrt{\frac{2}{\pi
n}},$ (61)

which is superseded by numerical bounds when $n\geq 3$ (see Table 1).

We now combine Lemma 6 and Lemma 7 to prove Theorem 5.

###### Proof.

The feasible set of (D-SDP${}_{n}^{m,\leq}$) is non-empty, by considering the
null sequence, which achieves value $1$. Without loss of generality, we add
the constraint $y\leq 1$ in (D-LP${}_{n}^{\mathcal{S}}$) and
(D-SDP${}_{n}^{m,\leq}$).

Let $m\geq n$ and let $(y,\bm{\mu})\in\mathbb{R}\times S^{\prime}(\mathbb{N})$
be a feasible solution of (D-SDP${}_{n}^{m,\leq}$) expressed in the form given
by Lemma 6. The constraint $\braket{f_{\bm{\mu}},x\mapsto
e^{-\frac{x}{2}}}\geq 0$ implies $\mu_{0}\geq 0$ and thus $y\geq 0$. Without
loss of generality, we may set $\mu_{k}=0$ for $k>m$, since these coefficients
are only constrained by $\mu_{k}\leq y\leq 1$. We also have $\mu_{k}\leq 1$
for all $k\in\mathbb{N}$.

By Lemma 7, $\bm{F}^{l}$ is a feasible solution of (SDP${}_{l}^{m,\leq}$) for
all $l\leq m$, so in particular
$f_{\bm{F}^{l}}=\sum_{k=0}^{l}F_{k}^{l}\mathcal{L}_{k}\in\mathcal{R}_{m,+}(\mathbb{R}_{+})$.
Hence, $\bm{\mu}$ satisfies the constraint
$\braket{f_{\bm{\mu}},f_{\bm{F}^{l}}}\geq 0$, which gives

$\sum_{k=0}^{l}\mu_{k}F_{k}^{l}\geq 0,$ (62)

for all $l\leq m$. Thus we have, for all $l\in\llbracket 1,m\rrbracket$,

$\displaystyle\mu_{l}$
$\displaystyle\geq-\frac{1}{F_{l}^{l}}\sum_{k=0}^{l-1}\mu_{k}F_{k}^{l}$ (63)
$\displaystyle\geq-\frac{1}{F_{l}^{l}}\sum_{k=0}^{l-1}F_{k}^{l}$
$\displaystyle=1-\frac{1}{F_{l}^{l}}$ $\displaystyle\geq-l,$

where we used $F_{l}^{l}>0$ in the first line, $\mu_{k}\leq 1$ and
$F_{k}^{l}\geq 0$ in the second line, $\sum_{k=0}^{l}F_{k}^{l}=1$ in the third
line and Eq. (60) in the last line. With $\mu_{k}\leq 1$, we obtain
$|\mu_{k}|\leq k$ for $k\in\mathbb{N}^{*}$, and thus $|\mu_{k}|\leq
k+1$666Critically, this bound does not depend on the level $m$ of the
hierarchy.  for all $k\in\mathbb{N}$. Hence, the feasible set of
(D-SDP${}_{n}^{m,\leq}$) is compact and the program (D-SDP${}_{n}^{m,\leq}$)
has feasible optimal solutions for all $m\geq n$, by diagonal extraction.

Let $(y^{m},\bm{\mu}^{m})_{m\geq n}$ be a sequence of optimal solutions of
(D-SDP${}_{n}^{m,\leq}$), for $m\geq n$. By Theorem 2, we have strong duality
between the programs $(\text{SDP}^{m,\leq}_{n})$ and (D-SDP${}_{n}^{m,\leq}$),
so the optimal value of (D-SDP${}_{n}^{m,\leq}$) is given by
$\omega_{n}^{m,\leq}$, for all $m\geq n$. By optimality
$y^{m}=\omega_{n}^{m,\leq}$, for all $m\geq n$, and the sequence
$(y^{m})_{m\geq n}$ converges.

Performing a diagonal extraction $\phi$ on the sequence $(\bm{\mu}^{m})_{m\geq
n}$, we obtain a sequence of sequences $(\bm{\mu}^{\phi(m)})_{m\geq n}$ such
that each sequence $(\mu^{\phi(m)}_{k})_{m\geq n}$ converges when
$m\rightarrow+\infty$, for all $k\in\mathbb{N}$. Let $\mu_{k}$ denote its
limit, for each $k\in\mathbb{N}$. We write
$\bm{\mu}=(\mu_{k})_{k\in\mathbb{N}}\in\mathbb{R}^{\mathbb{N}}$ the sequence
of limits. We also write

$y:=\lim_{m\rightarrow+\infty}y^{m}=\lim_{m\rightarrow+\infty}\omega_{n}^{m,\leq}.$
(64)

For all $m\geq n$, we have $\omega_{n}^{\phi(m),\leq}\geq\mu_{k}^{\phi(m)}$
for all $k\in\mathbb{N}$ and $\omega_{n}^{\phi(m),\leq}\geq
1+\mu_{n}^{\phi(m)}$, so taking $m\rightarrow+\infty$ we obtain $y\geq\mu_{k}$
for all $k\in\mathbb{N}$ and $y\geq 1+\mu_{n}$. By Eq. (57), we have
$\omega^{m,\leq}_{n}\leq\omega_{n}^{\mathcal{S}}$ for all $m\geq n$, so
$y\leq\omega_{n}^{\mathcal{S}}$.

Moreover, $|\mu_{k}^{\phi(m)}|\leq k+1$ for all $k\in\mathbb{N}$, so taking
$m\rightarrow+\infty$ we obtain $|\mu_{k}|\leq k+1$ for all $k\in\mathbb{N}$,
which implies that $\bm{\mu}\in\mathcal{S}^{\prime}(\mathbb{N})$ [64]. Let
$f_{\bm{\mu}}=\sum_{k}\mu_{k}\mathcal{L}_{k}\in\mathcal{S}^{\prime}(\mathbb{R}_{+})$.
We have

$\mu_{k}=\braket{f_{\bm{\mu}},\mathcal{L}_{k}}.$ (65)

By construction we also have:

$\forall m\geq n,\forall
g\in\mathcal{R}_{m,+}(\mathbb{R}_{+}),\;\braket{f_{\bm{\mu}},g}\geq 0.$ (66)

By Theorem 1, this implies that the distribution
$\mu:=f_{\bm{\mu}}(x)\in\mathcal{S}^{\prime}(\mathbb{R}_{+})$ is non-negative,
i.e.,

$\forall f\in\mathcal{S}_{+}(\mathbb{R}_{+}),\;\braket{f_{\bm{\mu}},f}\geq 0.$
(67)

With the constraints $y\geq\mu_{k}$ for all $k\in\mathbb{N}$ and $y\geq
1+\mu_{n}$, we have that $(y,\mu)$ is a feasible solution of
(D-LP${}_{n}^{\mathcal{S}}$), and in particular
$y\geq{\omega_{n}^{\prime\mathcal{S}}}$, since (D-LP${}_{n}^{\mathcal{S}}$) is
a minimisation problem. Since $y\leq\omega_{n}^{\mathcal{S}}$ we obtain with
Eq. (57) and Eq. (64):

$\lim_{m\rightarrow+\infty}\omega^{m,\leq}_{n}=\omega_{n}^{\mathcal{S}}={\omega_{n}^{\prime\mathcal{S}}}.$
(68)

∎

As a direct corollary of the proofs of Theorem 4 and Theorem 5 (in the same
spirit than the remark after the proof of Theorem 4), we obtain the following
strong duality result:

######  Theorem 6.

Strong duality holds between the programs (LP${}_{n}^{\mathcal{S}}$) and
(D-LP${}_{n}^{\mathcal{S}}$) and between the programs (LP${}_{n}^{L^{2}}$) and
(D-LP${}_{n}^{L^{2}}$).

For completeness, we give a different and more direct proof of the strong
duality between (LP${}_{n}^{L^{2}}$) and (D-LP${}_{n}^{L^{2}}$) in Appendix G.

We have shown the convergence of the semidefinite hierarchies of upper and
lower bounds: $(\omega^{m,\geq}_{n})_{m\geq n}$ towards the optimal value of
(LP${}_{n}^{L^{2}}$) and $(\omega^{m,\leq}_{n})_{m\geq n}$ towards the optimal
value of (LP${}_{n}^{\mathcal{S}}$). By linearity, these results generalise
straightforwardly to the case of witnesses corresponding to linear
combinations of fidelities with displaced Fock states:

$\displaystyle\lim_{m\rightarrow+\infty}\\!\omega^{m,\leq}_{\bm{a}}=\omega_{\bm{a}}^{\mathcal{S}},$
(69)
$\displaystyle\lim_{m\rightarrow+\infty}\\!\omega^{m,\geq}_{\bm{a}}=\omega_{\bm{a}}^{L^{2}},$
(70)

for all $n\in\mathbb{N}^{*}$ and all $\bm{a}\in[0,1]^{n}$.

We obtained two hierarchies providing numerical lower bounds and upper bounds
on the threshold value. A natural question that arises is the following: is
there a gap between the optimal values of (LP${}_{n}^{L^{2}}$) and
(LP${}_{n}^{\mathcal{S}}$)? We leave this as an open question.

##  6 Witnessing multimode Wigner negativity

In this section we discuss the generalisation of our Wigner negativity
witnesses to the more challenging multimode setting. Hereafter, $M$ denotes
the number of modes.

###  6.1 Multimode Wigner negativity witnesses

Using multi-index notations (see Appendix H.1), the single-mode Wigner
negativity witnesses defined in Eq. (6) are naturally generalised to

$\hat{\Omega}_{\bm{a},\bm{\alpha}}:=\smashoperator[]{\sum_{\bm{1}\leq\bm{k}\leq\bm{n}}^{}}a_{\bm{k}}\hat{D}(\bm{\alpha})\ket{\bm{k}}\\!\bra{\bm{k}}\hat{D}^{\dagger}(\bm{\alpha}),$
(71)

for $\bm{n}=(n_{1},\dots,n_{M})\in\mathbb{N}^{M}\setminus\\{\bm{0}\\}$,
$\bm{a}=(a_{\bm{k}})_{\bm{1}\leq\bm{k}\leq\bm{n}}\in[0,1]^{n_{1}\cdots
n_{M}}$, with $\max_{\bm{k}}a_{\bm{k}}=1$, and $\bm{\alpha}\in\mathbb{C}^{M}$.
Similar to the single-mode case, these POVM elements are weighted sums of
multimode displaced Fock states projectors, and their expectation value for a
quantum state $\bm{\rho}\in\mathcal{D}(\mathcal{H}^{\otimes M})$ is given by

$\text{Tr}(\hat{\Omega}_{\bm{a},\bm{\alpha}}\bm{\rho})=\smashoperator[]{\sum_{\bm{1}\leq\bm{k}\leq\bm{n}}^{}}a_{\bm{k}}F\left(\hat{D}^{\dagger}(\bm{\alpha})\bm{\rho}\hat{D}(\bm{\alpha}),\ket{\bm{k}}\right)$
(72)

where $F$ is the fidelity. Unlike in the single-mode case however, estimating
this quantity with homodyne or heterodyne detection by direct fidelity
estimation is no longer efficient when the number of modes becomes large.
Instead, one may use robust lower bounds on the multimode fidelity from [44]
which can be obtained efficiently with homodyne or heterodyne detection. A
lower bound on the estimated experimental multimode fidelity will allow to
detect Wigner negativity if it is larger than an upper bound on the threshold
value associated to a given witness.

These lower bounds are obtained as follows: given a target multimode Fock
state $\ket{\bm{n}}=\ket{n_{1}}\otimes\cdots\otimes\ket{n_{M}}$ and multiple
copies of an $M$-mode experimental state $\bm{\rho}$, measure all single-mode
subsystems of $\bm{\rho}$ and perform fidelity estimation with each
corresponding single-mode target Fock state. That is, the samples obtained
from the detection of the $i^{th}$ mode of $\bm{\rho}$ are used for single-
mode fidelity estimation with the Fock state $\ket{n_{i}}$. Let
$F_{1},\dots,F_{M}$ be the single-mode fidelity estimates obtained and let
$\tilde{F}(\bm{\rho},\ket{\bm{n}}):=1-\sum_{i=1}^{M}(1-F_{i})$. Then [44],

$1\\!-\\!M(1\\!-\\!F(\bm{\rho},\ket{\bm{n}}))\leq\tilde{F}(\bm{\rho},\ket{\bm{n}})\leq
F(\bm{\rho},\ket{\bm{n}}).$ (73)

In particular, $\tilde{F}$ provides a good estimate of the multimode fidelity
$F$ whenever $F$ is not too small.

The same procedure is followed in the case of target displaced Fock states,
with classical translations of the samples in order to account for the
displacement parameters.

To each witness $\hat{\Omega}_{\bm{a},\bm{\alpha}}$ is associated its
threshold value:

$\omega_{\bm{a}}:=\sup_{\begin{subarray}{c}\bm{\rho}\in\mathcal{D}(\mathcal{H}^{\otimes
M})\\\ W_{\bm{\rho}}\geq
0\end{subarray}}\text{Tr}\\!\left(\hat{\Omega}_{\bm{a},\bm{\alpha}}\,\bm{\rho}\right).$
(74)

With Eq. (73), if the value of the bound $\tilde{F}$ obtained experimentally
is greater than $\omega_{\bm{a}}$, then the state $\bm{\rho}$ has a negative
Wigner function.

With the same arguments as in the single-mode case, the multimode Wigner
negativity witnesses in Eq. (71) form a complete family and retain the
interpretation from Lemma 1: the violation of the threshold value provides a
lower bound on the distance to the set of multimode states with non-negative
Wigner function. However, the limited robustness of the bound $\tilde{F}$ may
affect the performance of the witnesses in practical scenarios, in particular
for witnesses that are sums of different projectors. Still, we show in section
6.3 the applicability of the method with a genuinely multimode example.

We first generalise the single-mode semidefinite programming approach for
approximating the threshold values to the multimode case.

###  6.2 Approximating the multimode threshold values

By linearity, we restrict our analysis to the case of Wigner negativity
witnesses that are projectors onto a single multimode Fock state
$\ket{\bm{n}}$, for $\bm{n}\in\mathbb{N}^{M}\setminus\\{\bm{0}\\}$. We thus
consider the computation of

$\omega_{\bm{n}}=\sup_{\begin{subarray}{c}\bm{\rho}\in\mathcal{D}(\mathcal{H}^{\otimes
M})\\\ W_{\bm{\rho}}\geq 0\end{subarray}}\braket{\bm{n}}{\bm{\rho}}{\bm{n}}.$
(75)

A similar reasoning to the single-mode case shows that the computation of the
corresponding threshold value in Eq. (74) may be rephrased as the following
infinite-dimensional linear program:

$\displaystyle\quad\text{Find
}(F_{\bm{k}})_{\bm{k}\in\mathbb{N}^{M}}\in\ell^{2}$ (LP${}_{\bm{n}}^{L^{2}}$)
$\displaystyle\quad\text{maximising }F_{\bm{n}}$ subject to
$\displaystyle\begin{dcases}\sum_{\bm{k}}F_{\bm{k}}=1\\\
\forall\bm{k}\in\mathbb{N}^{M},\;F_{\bm{k}}\geq 0\\\
\forall\bm{x}\in\mathbb{R}_{+}^{M},\;\sum_{\bm{k}}F_{\bm{k}}\mathcal{L}_{\bm{k}}(\bm{x})\geq
0,\end{dcases}$

where the optimisation is over square-summable real sequences indexed by
elements of $\mathbb{N}^{M}$. Its dual linear program reads

$\displaystyle\quad\text{Find }y\in\mathbb{R}\text{ and
}\mu\in{L^{2}}^{\prime}(\mathbb{R}_{+}^{M})$ (D-LP${}_{\bm{n}}^{L^{2}}$)
$\displaystyle\quad\text{minimising }y$ subject to
$\displaystyle\begin{dcases}\forall\bm{k}\neq\bm{n}\in\mathbb{N}^{M},\;y\geq\int_{\mathbb{R}_{+}^{M}}{\mathcal{L}_{\bm{k}}}{d\mu}\\\
y\geq 1+\int_{\mathbb{R}_{+}^{M}}{\mathcal{L}_{\bm{n}}}{d\mu}\\\ \forall f\in
L^{2}_{+}(\mathbb{R}_{+}^{M}),\;\langle\mu,f\rangle\geq 0.\end{dcases}$

We denote their optimal value $\omega_{\bm{n}}^{L^{2}}$—as in the single-mode
case, we have strong duality between these programs. We also introduce the
programs over Schwartz functions:

$\displaystyle\quad\text{Find
}(F_{\bm{k}})_{\bm{k}\in\mathbb{N}^{M}}\in\mathcal{S}(\mathbb{N}^{M})$
(LP${}_{\bm{n}}^{\mathcal{S}}$) $\displaystyle\quad\text{maximising
}F_{\bm{n}}$ subject to
$\displaystyle\begin{dcases}\sum_{\bm{k}}F_{\bm{k}}=1\\\
\forall\bm{k}\in\mathbb{N}^{M},\;F_{\bm{k}}\geq 0\\\
\forall\bm{x}\in\mathbb{R}_{+}^{M},\;\sum_{\bm{k}}F_{\bm{k}}\mathcal{L}_{\bm{k}}(\bm{x})\geq
0,\end{dcases}$

and

$\displaystyle\quad\text{Find }y\in\mathbb{R}\text{ and
}\mu\in{\mathcal{S}}^{\prime}(\mathbb{R}_{+}^{M})$
(D-LP${}_{\bm{n}}^{\mathcal{S}}$) $\displaystyle\quad\text{minimising }y$
subject to
$\displaystyle\begin{dcases}\forall\bm{k}\neq\bm{n}\in\mathbb{N}^{M},\;y\geq\int_{\mathbb{R}_{+}^{M}}{\mathcal{L}_{\bm{k}}}{d\mu}\\\
y\geq 1+\int_{\mathbb{R}_{+}^{M}}{\mathcal{L}_{\bm{n}}}{d\mu}\\\ \forall f\in
L^{2}_{+}(\mathbb{R}_{+}^{M}),\;\langle\mu,f\rangle\geq 0.\end{dcases}$

We denote their optimal value $\omega_{\bm{n}}^{S}$.

In the single-mode case, we obtained hierarchies of SDP relaxations and
restrictions for (LP${}_{n}^{L^{2}}$) by replacing constraints involving non-
negative functions by constraints involving non-negative polynomials $P$ of
fixed degree. We then exploited the existence of a sum-of-squares
decomposition for non-negative monovariate polynomials. In the multimode
setting, the polynomials involved are multivariate, so that the set of non-
negative polynomials over $\mathbb{R}$ of a given degree may be strictly
larger than the set of sum-of-square polynomials [68]. Instead, we replace
directly constraints involving non-negative functions over $\mathbb{R}_{+}$ by
constraints involving non-negative polynomials $P$ of fixed degree such that
$\bm{x}\mapsto P(\bm{x}^{2})$888We write $\bm{x}^{2}$ in short for
$\bm{x}^{2\bf 1}=(x_{1}^{2},\dots,x_{M}^{2})$. has a sum-of-squares
decomposition, implying that the multimode semidefinite relaxations and
restrictions are possibly looser than their single-mode counterparts.

Moreover, the dimension of the semidefinite programs increases exponentially
with the level of the hierarchy $m$, as the number of $M$-variate monomials of
degree less or equal to $m$ is given by $\binom{M+m}{m}$. This implies that
the semidefinite programs remain tractable only for a constant number of
levels.

In spite of these observations, and following similar steps to the single-mode
case (see Appendix H.2), the SDP relaxations providing upper bounds for the
threshold value are given by

$\displaystyle\quad\text{Find }A=(A_{\bm{i}\bm{j}})_{|\bm{i}|,|\bm{j}|\leq
m}\in\text{Sym}_{\binom{M+m}{m}}$ $(\text{SDP}^{m,\geq}_{\bm{n}})$
$\displaystyle\quad\text{and }\bm{F}=(F_{\bm{k}})_{|\bm{k}|\leq
m}\in\mathbb{R}^{{\binom{M+m}{m}}}$ $\displaystyle\quad\text{maximising
}F_{\bm{n}}$ subject to
$\displaystyle\begin{dcases}\textstyle\sum_{|\bm{k}|\leq m}F_{\bm{k}}=1\\\
\textstyle\forall|\bm{k}|\leq m,\;F_{\bm{k}}\geq 0\\\
\textstyle\forall|\bm{l}|\leq
m,\forall\bm{i}\\!+\\!\bm{j}\\!=\\!2\bm{l},\;A_{\bm{i}\bm{j}}=\sum\limits_{\bm{k}\leq\bm{l}}F_{\bm{k}}\binom{\bm{l}}{\bm{k}}\bm{l}!\\\
\textstyle\forall|\bm{r}|\leq 2m,\bm{r}\\!\neq\\!2\bm{l},\forall|\bm{l}|\leq
m,\forall\bm{i}\\!+\\!\bm{j}\\!=\\!\bm{r},\;A_{\bm{i}\bm{j}}=0\\\ \textstyle
A\succeq 0,\end{dcases}$

for all $m\geq|\bm{n}|$. Similarly, the semidefinite restrictions providing
lower bounds for the threshold value are given by

$\displaystyle\quad\text{Find }Q=(Q_{\bm{i}\bm{j}})_{|\bm{i}|,|\bm{j}|\leq
m}\in\text{Sym}_{\binom{M+m}{m}}$ $(\text{SDP}^{m,\leq}_{\bm{n}})$
$\displaystyle\quad\text{and }\bm{F}=(F_{\bm{k}})_{|\bm{k}|\leq
m}\in\mathbb{R}^{{\binom{M+m}{m}}}$ $\displaystyle\quad\text{maximising
}F_{\bm{n}}$ subject to
$\displaystyle\begin{dcases}\textstyle\sum_{|\bm{k}|\leq m}F_{\bm{k}}=1\\\
\textstyle\forall|\bm{k}|\leq m,\;F_{\bm{k}}\geq 0\\\
\textstyle\forall|\bm{l}|\leq
m,\sum_{\bm{i}+\bm{j}=2\bm{l}}Q_{\bm{i}\bm{j}}=\sum_{\bm{k}\geq\bm{l}}\frac{(-1)^{|\bm{k}|+|\bm{l}|}}{\bm{l}!}\binom{\bm{k}}{\bm{l}}F_{\bm{k}}\\\
\textstyle\forall|\bm{r}|\leq 2m,\bm{r}\neq 2\bm{l},\forall|\bm{l}|\leq
m,\sum_{\bm{i}+\bm{j}=\bm{r}}Q_{\bm{i}\bm{j}}=0\\\ \textstyle Q\succeq
0,\end{dcases}$

for all $m\geq|\bm{n}|$. In these semidefinite programs, the optimisations are
over matrices and vectors indexed by elements of $\mathbb{N}^{m}$ with sum of
coefficients lower that $m$.

While our proof of convergence of the single-mode hierarchy of upper bounds
towards $\omega_{\bm{n}}^{L^{2}}$ transfers easily to the multimode setting,
the proof of convergence of the hierarchy of lower bounds towards
$\omega_{\bm{n}}^{\mathcal{S}}$ requires the analytical expression of feasible
solutions for each level of the hierarchy. We show how to construct such
solutions in the multimode case using products of single-mode feasible
solutions—this requires introducing an equivalent hierarchy of restrictions,
where constraints are expressed on polynomials of $M$ variables with the
degree in each individual variable being less or equal to $m$, rather than on
polynomials of degree $m$ (that is, constraints of the form $\bm{k}\leq\bm{m}$
where $\bm{m}=(m,\dots,m)\in\mathbb{N}^{M}$ rather than $|\bm{k}|\leq m$).
Along the way, we also prove strong duality of the programs involved. We refer
to Appendix H.3 for the proofs.

Summarising our results, we find that the multimode semidefinite programs
$(\text{SDP}^{m,\geq}_{\bm{n}})$ and $(\text{SDP}^{m,\leq}_{\bm{n}})$ provide
converging sequences of upper and lower bounds to the threshold values
$\omega_{\bm{n}}^{L^{2}}$ and $\omega_{\bm{n}}^{\mathcal{S}}$, respectively.
We study a concrete application in the next section.

###  6.3 Multimode example

To illustrate the usefulness of our Wigner negativity witnesses in the
multimode setting, we consider a lossy Fock state over two modes:

$\displaystyle\bm{\rho}_{1,1,\eta}$
$\displaystyle:=(1-\eta)^{2}\ket{1}\\!\bra{1}\otimes\ket{1}\\!\bra{1}$ (76)
$\displaystyle\;\;+\eta(1-\eta)\ket{1}\\!\bra{1}\otimes\ket{0}\\!\bra{0}$
$\displaystyle\;\;+\eta(1-\eta)\ket{0}\\!\bra{0}\otimes\ket{1}\\!\bra{1}$
$\displaystyle\;\;+\eta^{2}\ket{0}\\!\bra{0}\otimes\ket{0}\\!\bra{0},$

with loss parameter $0\leq\eta\leq 1$. Setting $\eta=0$ gives
$\bm{\rho}_{1,1,\eta}=\ket{1}\\!\bra{1}\otimes\ket{1}\\!\bra{1}$ while setting
$\eta=1$ gives
$\bm{\rho}_{1,1,\eta}=\ket{0}\\!\bra{0}\otimes\ket{0}\\!\bra{0}$. This state
has a non-negative Wigner function for $\eta\geq\frac{1}{2}$.

We also consider the multimode Wigner negativity witness
$\ket{1}\\!\bra{1}\otimes\ket{1}\\!\bra{1}$, which is a projector onto the
Fock state $\ket{1}\otimes\ket{1}$. Solving numerically the corresponding
hierarchy $(\text{SDP}^{m,\leq}_{\bm{n}})$ up to $m=3$, we obtain the lower
bound $0.266$ and solving the hierarchy $(\text{SDP}^{m,\geq}_{\bm{n}})$ up to
$m=10$, we obtain the upper bound $0.320$.

A direct consequence of the numerical lower bound is that tensor product
states are not the closest among Wigner positive states to tensor product
states with a negative Wigner function. Indeed, the maximum achievable
fidelity with the state $\ket{1}\otimes\ket{1}$ using Wigner positive tensor
product states is equal to the square of the maximum achievable fidelity with
the state $\ket{1}$ using single-mode Wigner positive states, that is
$0.5^{2}=0.25<0.266$.

Figure 5: Witnessing Wigner negativity of the lossy Fock state
$\bm{\rho}_{1,1,\eta}$ over two modes using the witness
$\ket{1}\\!\bra{1}\otimes\ket{1}\\!\bra{1}$. The threshold value for that
witness is upper bounded by $0.320$ and lower bounded by $0.266$. The dashed
red line delimits the interval of loss parameter values where the witness can
be used to detect Wigner negativity of $\bm{\rho}_{1,1,\eta}$ efficiently,
i.e., when the robust bound
$\tilde{F}(\bm{\rho}_{1,1,\eta},\ket{1}\otimes\ket{1})$ (blue curve) on the
fidelity from Eq. (78) is above the witness upper bound (red line). When it is
below the witness lower bound (black line), we are guaranteed that the witness
cannot be used to detect Wigner negativity of the state. The fidelity
$F(\bm{\rho}_{1,1,\eta},\ket{1}\otimes\ket{1})$ is also depicted above (yellow
curve). Note that $\bm{\rho}_{1,1,\eta}$ has a non-negative Wigner function
for $\eta\geq 0.5$.

We now use the upper bound to witness the Wigner negativity of the state
$\bm{\rho}_{1,1\eta}$ (see Fig. 5). The fidelity between
$\bm{\rho}_{1,1,\eta}$ and $\ket{1}\otimes\ket{1}$ is given by
$F(\bm{\rho}_{1,1,\eta},\ket{1}\otimes\ket{1})=(1-\eta)^{2}$, for all
$0\leq\eta\leq 1$. This fidelity is above the upper bound $0.320$ on the
threshold value of the witness $\ket{1}\\!\bra{1}\otimes\ket{1}\\!\bra{1}$
when $\eta\leq 0.434$.

However, in practice one would not obtain a precise estimate of the fidelity
efficiently, but rather a robust lower bound on the fidelity computed from
single-mode fidelities, which satisfies Eq. (73). In the worst case, the
estimate obtained is closer to
$1-2(1-F(\bm{\rho}_{1,1,\eta},\ket{1}\otimes\ket{1}))$ than to
$F(\bm{\rho}_{1,1,\eta},\ket{1}\otimes\ket{1})$. When the value of this robust
lower bound is greater than the threshold value of the witness, this implies
that the state has a negative Wigner function.

In the present case, the two single-mode reduced states of
$\bm{\rho}_{1,1,\eta}$ are the same, given by

$\text{Tr}_{2}(\bm{\rho}_{1,1,\eta})=(1-\eta)\ket{1}\\!\bra{1}+\eta\ket{0}\\!\bra{0},$
(77)

so the single-mode fidelities with $\ket{1}$ are equal for each mode and given
by $1-\eta$. Hence, the robust lower bound on
$F(\bm{\rho}_{1,1,\eta},\ket{1}\otimes\ket{1})$ is given by

$\tilde{F}(\bm{\rho}_{1,1,\eta},\ket{1}\otimes\ket{1})=1-2\eta.$ (78)

It is above the upper bound $0.320$ on the threshold value of the witness
$\ket{1}\\!\bra{1}\otimes\ket{1}\\!\bra{1}$ when $\eta\leq 0.340$.

This example highlights the use of efficient and robust lower bounds on
multimode fidelities rather than fidelity estimates [44], in conjunction with
our family of multimode witnesses to detect Wigner negativity of realistic
experimental states.

##  7 Conclusion and open problems

Characterising quantum properties of physical systems is an important step in
the development of quantum technologies, and negativity of the Wigner
function, a necessary resource for any quantum computational speedup, is no
exception. In this work, we have derived a complete family of Wigner
negativity witnesses which provide an operational quantification of Wigner
negativity, both in the single-mode and multimode settings. In the context of
quantum optical information processing, the main application of our method is
in experimental scenarios, where it leads to robust and efficient
certification of negativity of the Wigner function. Witnesses of Wigner
negativity also provide witnesses of non-classicality although they are no
longer complete in this case.

What is more, our witnesses also delineate the set of quantum states with
positive Wigner function, and it would be interesting to understand whether
additional insights on this set can be obtained using these witnesses.

The Wigner function has been extended to the discrete-variable setting [71,
72], where it has been been linked to contextuality [73, 74, 75], a necessary
resource for discrete-variable quantum computing [76, 77]. A framework for
treating contextuality and computing the amount of contextuality in
continuous-variable settings has recently been developed [78]. As we obtained
reliable Wigner negativity witnesses, it would be interesting to investigate
the link between Wigner negativity and continuous-variable contextuality.

Hierarchies of semidefinite programs (in particular with non-commutative
variables [79]) have found many recent applications in quantum information
theory. From an infinite-dimensional linear program, we were able to use
numerically both a hierarchy of upper bounds and a hierarchy of lower
bounds—thus obtaining a certificate for the optimality of these bounds by
looking at their difference—whereas this only works in specific cases for the
Lasserre hierarchy of upper bounds [35]. Can we find other interesting cases
where we can exploit both hierarchies? Moreover, we obtained an analytical
sequence of lower bounds for the threshold value of the program
(LP${}_{n}^{L^{2}}$). Can we also get an analytical sequence of upper bounds?
In particular, we anticipate that Fock states $\ket{n}$ get further away from
the set of states having a positive Wigner function as $n$ increases and that
$\omega_{n}=\mathcal{O}(\frac{1}{\sqrt{n}})$ as $n\rightarrow+\infty$. The
question of the gap between the optimisation over square-integrable functions
and Schwartz functions is also left open.

Finally, using our multimode Wigner negativity witnesses for studying the
interplay between Wigner negativity and entanglement [52] is a very
interesting prospect which we leave for future work.

### Acknowledgments

U. Chabaud acknowledges stimulating discussions with S. Gribling, T. Freiman
and T. Vidick. P.-E. Emeriau acknowledges interesting discussions with A.
Oustry, E. Galvão and R. Soares Barbosa. We thank J. Eisert for his valuable
comments on a previous version of this work and P. Paule for providing access
to the Mathematica package for implementing Zeilberger’s algorithm. U. Chabaud
acknowledges funding provided by the Institute for Quantum Information and
Matter, an NSF Physics Frontiers Center (NSF Grant PHY-1733907). F. Grosshans
acknowledges funding from the ANR through the ANR-17-CE24-0035 VanQuTe
project.

## References

  * [1] S. Lloyd and S. L. Braunstein, “Quantum computation over continuous variables,” in Quantum Information with Continuous Variables, pp. 9–17.  Springer, 1999.
  * [2] 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).
  * [3] U. Leonhardt, “Essential Quantum Optics,”.  Cambridge University Press, Cambridge, UK, 1st ed., 2010.
  * [4] J. E. Moyal, “Quantum mechanics as a statistical theory,” in Mathematical Proceedings of the Cambridge Philosophical Society, vol. 45, pp. 99–124, Cambridge University Press.  1949\.
  * [5] E. P. Wigner, “On the quantum correction for thermodynamic equilibrium,” in Part I: Physical Chemistry. Part II: Solid State Physics, pp. 110–120.  Springer, 1997.
  * [6] C. T. Lee, “Measure of the nonclassicality of nonclassical states,” Physical Review A 44, R2775 (1991).
  * [7] G. Giedke and J. I. Cirac, “Characterization of Gaussian operations and distillation of Gaussian states,” Physical Review A 66, 032316 (2002).
  * [8] J. Eisert, S. Scheel, and M. B. Plenio, “Distilling Gaussian states with Gaussian operations is impossible,” Physical Review Letters 89, 137903 (2002).
  * [9] J. Fiurášek, “Gaussian transformations and distillation of entangled Gaussian states,” Physical Review Letters 89, 137904 (2002).
  * [10] J. Niset, J. Fiurášek, and N. J. Cerf, “No-go theorem for Gaussian quantum error correction,” Physical Review Letters 102, 120501 (2009).
  * [11] S. Ghose and B. C. Sanders, “Non-Gaussian ancilla states for continuous variable quantum computation via Gaussian maps,” Journal of Modern Optics 54, 855–869 (2007).
  * [12] S. D. Bartlett, B. C. Sanders, S. L. Braunstein, and K. Nemoto, “Efficient Classical Simulation of Continuous Variable Quantum Information Processes,” Physical Review Letters 88, 097904 (2002).
  * [13] U. Chabaud, G. Ferrini, F. Grosshans, and D. Markham, “Classical simulation of Gaussian quantum circuits with non-Gaussian input states,” arXiv:2010.14363.
  * [14] R. L. Hudson, “When is the Wigner quasi-probability density non-negative?,” Reports on Mathematical Physics 6, 249–252 (1974).
  * [15] F. Soto and P. Claverie, “When is the Wigner function of multidimensional systems nonnegative?,” Journal of Mathematical Physics 24, 97–100 (1983).
  * [16] A. Mandilara, E. Karpov, and N. Cerf, “Extending Hudson’s theorem to mixed quantum states,” Physical Review A 79, 062302 (2009).
  * [17] R. Filip and L. Mišta Jr, “Detecting quantum states with a positive Wigner function beyond mixtures of Gaussian states,” Physical Review Letters 106, 200401 (2011).
  * [18] K. C. Tan, S. Choi, and H. Jeong, “Negativity of quasiprobability distributions as a measure of nonclassicality,” Physical review letters 124, 110404 (2020).
  * [19] U. Titulaer and R. Glauber, “Correlation functions for coherent fields,” Physical Review 140, B676 (1965).
  * [20] A. Kenfack and K. Życzkowski, “Negativity of the Wigner function as an indicator of non-classicality,” Journal of Optics B: Quantum and Semiclassical Optics 6, 396 (2004).
  * [21] A. Mari and J. Eisert, “Positive Wigner Functions Render Classical Simulation of Quantum Computation Efficient,” Physical Review Lett. 109, 230503 (2012).
  * [22] L. García-Álvarez, C. Calcluth, A. Ferraro, and G. Ferrini, “Efficient simulatability of continuous-variable circuits with large Wigner negativity,” arXiv:2005.12026.
  * [23] J. Preskill, “Quantum Computing in the NISQ era and beyond,” Quantum 2, 79 (2018).
  * [24] J. Eisert, D. Hangleiter, N. Walk, I. Roth, D. Markham, R. Parekh, U. Chabaud, and E. Kashefi, “Quantum certification and benchmarking,” Nature Reviews Physics 2, 382–390 (2020).
  * [25] G. M. D’Ariano, M. G. Paris, and M. F. Sacchi, “Quantum tomography,” Advances in Imaging and Electron Physics 128, 206–309 (2003), arXiv:quant-ph/0302028.
  * [26] A. I. Lvovsky and M. G. Raymer, “Continuous-variable optical quantum-state tomography,” Reviews of Modern Physics 81, 299 (2009).
  * [27] U. Chabaud, T. Douce, F. Grosshans, E. Kashefi, and D. Markham, “Building Trust for Continuous Variable Quantum States,” in 15th Conference on the Theory of Quantum Computation, Communication and Cryptography.  2020\.
  * [28] B. M. Terhal, “A family of indecomposable positive linear maps based on entangled quantum states,” Linear Algebra and its Applications 323, 61–73 (2001).
  * [29] M. Lewenstein, B. Kraus, J. I. Cirac, and P. Horodecki, “Optimization of entanglement witnesses,” Physical Review A 62, 052310 (2000).
  * [30] A. Mari, K. Kieling, B. M. Nielsen, E. Polzik, and J. Eisert, “Directly estimating nonclassicality,” Physical Review Letters 106, 010403 (2011).
  * [31] T. Kiesel and W. Vogel, “Universal nonclassicality witnesses for harmonic oscillators,” Physical Review A 85, 062106 (2012).
  * [32] U. Chabaud, G. Roeland, M. Walschaers, F. Grosshans, V. Parigi, D. Markham, and N. Treps, “Certification of non-Gaussian states with operational measurements,” arXiv:2011.04320.
  * [33] J.-B. Lasserre, “Global optimization with polynomials and the problem of moments,” SIAM Journal on optimization 11, 796–817 (2001).
  * [34] P. A. Parrilo, Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization.  PhD thesis, California Institute of Technology, 2000.
  * [35] J. B. Lasserre, “A new look at nonnegativity on closed sets and polynomial optimization,” SIAM Journal on Optimization 21, 864–885 (2011).
  * [36] 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.
  * [37] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, “Gaussian quantum information,” Reviews of Modern Physics 84, 621 (2012).
  * [38] A. Wünsche, “Laguerre 2D-functions and their application in quantum optics,” Journal of Physics A: Mathematical and General 31, 8267 (1998).
  * [39] A. Royer, “Wigner function as the expectation value of a parity operator,” Physical Review A 15, 449 (1977).
  * [40] K. Banaszek, C. Radzewicz, K. Wódkiewicz, and J. Krasiński, “Direct measurement of the Wigner function by photon counting,” Physical Review A 60, 674 (1999).
  * [41] K. E. Cahill and R. J. Glauber, “Density operators and quasiprobability distributions,” Physical Review 177, 1882 (1969).
  * [42] K. Husimi, “Some formal properties of the density matrix,” Proceedings of the Physico-Mathematical Society of Japan. 3rd Series 22, 264–314 (1940).
  * [43] T. Richter, “Determination of photon statistics and density matrix from double homodyne detection measurements,” Journal of Modern Optics 45, 1735–1749 (1998).
  * [44] U. Chabaud, F. Grosshans, E. Kashefi, and D. Markham, “Efficient verification of Boson Sampling,” arXiv:2006.03520.
  * [45] A. Ferraro, S. Olivares, and M. G. Paris, “Gaussian states in continuous variable quantum information,” arXiv:quant-ph/0503237.
  * [46] 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).
  * [47] R. Takagi and Q. Zhuang, “Convex resource theory of non-Gaussianity,” Physical Review A 97, 062337 (2018).
  * [48] Q. Zhuang, P. W. Shor, and J. H. Shapiro, “Resource theory of non-Gaussian operations,” Physical Review A 97, 052317 (2018).
  * [49] U. Chabaud, D. Markham, and F. Grosshans, “Stellar representation of non-Gaussian quantum states,” Physical Review Letters 124, 063605 (2020).
  * [50] L. Vandenberghe and S. Boyd, “Semidefinite programming,” SIAM review 38, 49–95 (1996).
  * [51] J. Fiurášek and M. Ježek, “Witnessing negativity of Wigner function by estimating fidelities of catlike states from homodyne measurements,” Physical Review A 87, 062115 (2013).
  * [52] M. Walschaers, C. Fabre, V. Parigi, and N. Treps, “Entanglement and Wigner Function Negativity of Multimode Non-Gaussian States,” Physical Review Letters 119, 183601 (2017).
  * [53] U. Chabaud and P.-E. Emeriau, “Zeilberger’s algorithm and Hierarchy of semidefinite programs.” Software Heritage repository swh:1:dir:d98f70e386783ef69
bf8c2ecafdb7b328b19b7ec containing the numerical tools developed for this
article.

  * [54] A. Ourjoumtsev, R. Tualle-Brouri, J. Laurat, and P. Grangier, “Generating optical Schrödinger kittens for quantum information processing,” Science 312, 83–86 (2006).
  * [55] B. C. Sanders, “Entangled coherent states,” Physical Review A 45, 6811 (1992).
  * [56] W. H. Zurek, “Sub-Planck structure in phase space and its relevance for quantum decoherence,” Nature 412, 712–717 (2001).
  * [57] G. Sagnol and M. Stahlberg, “Picos, a python interface to conic optimization solvers,” in Proceedings of the in 21st International Symposium on Mathematical Programming.  2012\.
  * [58] M. ApS, MOSEK Optimizer API for Python 9.2.36, 2019.  https://docs.mosek.com/9.2/pythonapi/index.html.
  * [59] M. Nakata, “A numerical evaluation of highly accurate multiple-precision arithmetic version of semidefinite programming solver: SDPA-GMP,-QD and-DD.,” in 2010 IEEE International Symposium on Computer-Aided Control System Design, pp. 29–34, IEEE.  2010\.
  * [60] K. Fujisawa, M. Kojima, K. Nakata, and M. Yamashita, SDPA (SemiDefinite Programming Algorithm) User’s Manual—Version 6.2. 0, 2002.
  * [61] A. Barvinok, “A course in convexity,”, vol. 54 of Graduate Studies in Mathematics.  American Mathematical Society, 2002.
  * [62] G. Szegö, “Orthogonal Polynomials, revised ed,” in American Mathematical Society Colloquium Publications, vol. 23.  1959\.
  * [63] O. Nikodym, “Sur une généralisation des intégrales de M. J. Radon,” Fundamenta Mathematicae 15, 131–179 (1930).
  * [64] M. Guillemot-Teissier, “Développements des distributions en séries de fonctions orthogonales. Séries de Legendre et de Laguerre,” Annali della Scuola Normale Superiore di Pisa-Classe di Scienze 25, 519–573 (1971).