Nonparametric goodness-of fit testing in quantum homodyne tomography with noisy data

In the framework of quantum optics, we study the problem of goodness-of-fit testing in a severely ill-posed inverse problem. A novel testing procedure is introduced and its rates of convergence are investigated under various smoothness assumptions. The procedure is derived from a projection-type estimator, where the projection is done in $\mathbb{L}_2$ distance on some suitably chosen pattern functions. The proposed methodology is illustrated with simulated data sets.

In quantum mechanics, the results of physical measurements performed on a physical system prepared in a certain state are random.This fact provides a large variety of problems and mathematical and applied statistics offer many methods to answer them.For example, statisticians are asked to decide which measurement to perform in order to collect most informative observations or to rebuild the subjacent quantum object from the resulting observations.In particular, physicists have created laser in view of many applications (medical, . . .).An important question is to know whether they have generated the desired one.
In order to answer this question, we propose a testing procedure based on the indirect observations they gathered by a measurement called quantum homodyne tomography (QHT).To describe our framework and to motivate our work, we present in Section 1 the needed physical background.In Section 2, we introduce our statistical model and define the nonparametric class containing the density matrices we are dealing with.To make clear to physicists the statistical part, we recall some basic notions on the statistical tests in Section 2.2.The testing procedure and the theoretical results assessing its asymptotic behavior are presented in Section 3. Section 4 contains an illustration of the proposed methodology on simulated data sets while the proof are postponed to Section 5.

Physical background
Section 1.1 is a short introduction to quantum optics and we introduce in Section 1.2 the very useful and meaningful pattern functions.In Section 1.3, we present some examples of quantum states actually created in laboratory, on which our testing procedure may be applied.

Short introduction to quantum optics
In quantum mechanics, the quantum state fully describes all aspects of a physical system.Mathematically, the quantum state is described via density operators ρ on a complex Hilbert space H such that 1. ρ is self-adjoint (or Hermitian): ρ = ρ * , where ρ * is the adjoint of ρ, 2. ρ is positive: ρ ≥ 0, or equivalently ψ, ρψ ≥ 0 for all ψ ∈ H, 3. ρ has trace one: Tr(ρ) = 1.
Moreover, to each measurable physical property or quantity corresponds a selfadjoint operator, say X, on the space of states H.This operator is called an observable.Unlike in classical mechanics, when performing a certain measurement of an observable X of a physical system prepared in a quantum state ρ, the result is in general random and is described by a probability distribution of a random variable X.
In this paper, the studied quantum system is a monochromatic light in a cavity described by the state of a quantum harmonic oscillator.Hereafter, we consider the complex Hilbert space H = L 2 (R), the space of square integrable complex-valued functions on the real line.In quantum optics, the quantum harmonic oscillator can be represented by two operators on H = L 2 (R): the position operator Q (or electric field) and the quantity of movement operator P (or magnetic field) with ψ 1 , ψ 2 functions of H.The position Q and momentum P do not commute and satisfy the canonical commutation relation: where 1 stands for the identity operator on H.This relation may be understood as follows.The measurement of the position Q necessarily disturbs the particle's momentum P, and vice versa.They cannot be measured simultaneously.Thus, one cannot obtain a couple of random variables (Q, P ).However, a linear combination X φ = cos(φ)Q + sin(φ)P, for every phases φ ∈ [0, π], can theoretically be measured by a technique put in practice for the first time in [29] and called quantum homodyne tomography (QHT1 ), where the phase Φ is randomly and uniformly distributed on [0, π].
Furthermore, to each state ρ corresponds a Wigner function W ρ , which gives an equivalent representation of the quantum state ρ.The associated Wigner function W ρ can be defined rigorously by its Fourier transform F 2 [W ρ ] with respect to both variables The Wigner function is a mapping from R 2 to R such that W ρ (q, p)dqdp = 1.Note that the Wigner function can take negative values.For this reason, W ρ is a real-valued function regarded as a generalized joint probability density (quasiprobability density) of the two random variables Q and P that we would get if we could measure simultaneously the two observables Q and P, which are respectively given by equations ( 1) and (2).It is well-known that its Radon transform ℜ[W ρ ] is such that Here, ℜ denotes the Radon transform, taking functions W ρ (q, p) on R 2 into functions ℜ[W ρ ](x, φ) on R × [0, π] formed by integration along lines of direction φ and distance x from the origin, expressed as As Φ is chosen uniformly and independently on [0, π], we can define the probability density function p ρ (x, φ) of (X, Φ) w.r.t. the measure 1 π λ, where λ stands for the Lebesgue measure on R × [0, π] by

Fock basis and pattern functions
An important equation in physics, especially in quantum mechanics, is the Schrödinger equation.The time-independent Schrödinger equation is written as where w is the energy level.It turns out that the energies are "quantized", and may only take discrete values noted w k = k + 1/2, k ∈ N.For a given frequency w k , there are two fundamental solutions: ψ k and ϕ k .One is a normalized function, the function ψ k , which is such that ψ 2 k = 1 and is called the regular wave function.The other one, the function ϕ k , is called the irregular one as it cannot be normalizable as ψ k is.
The functions {ψ k } k∈N form a orthonormal basis of L 2 (R).This particular basis, physically very meaningful, is called the Fock basis and is written as follows: Here, H k (x) := (−1) k e x 2 d k dx k e −x 2 denotes the Hermite polynomial of degree k.The density operator defined in the previous section corresponds to a density matrix under some orthonormal basis.In the Fock basis, the matrix entries ρ j,k of the state ρ can be expressed as expected values of functions F j,k (X ℓ , Φ ℓ ) = f j,k (X ℓ )e −i(k−j)Φ ℓ , where f j,k = f k,j are bounded real functions called pattern functions [18].For all j, k ∈ N, where p ρ (x, φ) is the joint probability density of (X, Φ).In other words Equation (7) expresses the idea that one can reconstruct any density matrix element ρ j,k using the pattern functions.First introduced in [20], the pattern functions f j,k are well known in physics and are defined in [19] as the first derivatives of products of the two fundamental solutions ψ k and ϕ k of the Schrödinger equation given in (5) A concrete expression for their Fourier transform using generalized Laguerre polynomials can be found in [28,4].
where fj,k denotes the Fourier transform of the pattern function f j,k and L α k (x) denotes the generalized Laguerre polynomial.We note that the pattern functions f j,k (x) are even functions for even differences j − k and odd functions for odd ones f j,k (−x) = (−1) j−k f j,k (x).

Table 1 Examples of quantum states
, for j and k even, rest zero,

Examples of quantum states
We present in Table 1 examples of pure quantum states, which can be created at this moment in laboratory and belong to the class R(B, r, L) with r = 2.A state is called pure if it cannot be represented as a mixture (convex combination) of other states, i.e., if it is an extreme point of the convex set of states.This is equivalent to the density matrix being a one dimensional projector, i.e., of the form ρ = P ψ for some unit vector ψ.Equivalently, a state ρ is pure if Tr(ρ 2 ) = 1.All other states are called mixed states.Let us discuss these few examples of quantum states.Among the pure states we consider the single photon state and the vacuum state, which is the pure state of zero photons.Note that the vacuum state would provide a random variable of Gaussian probability density p ρ (x|φ) via the ideal measurement of quantum homodyne tomography.We consider also the coherent-q 0 state, which characterizes the laser pulse with the number of photons Poisson distributed with an average of M photons.The Squeezed states (see e.g.[7]) have Gaussian Wigner functions whose variances in the two directions have a fixed product.The parameters M and ξ are such that M ≥ sinh 2 (ξ), C(M, ξ) is a normalization constant, α = ((M − sinh 2 (ξ)) 1/2 )/(cosh(ξ) − sinh(ξ)), and δ = (α/(sinh(2ξ))) 1/2 .The Schrödinger cat state is described by a linear superposition of two coherent vectors (see e.g.[25]).Table 1 gives some explicit density matrix coefficients ρ j,k and probability densities p ρ (x|φ).

Problem formulation
In this paper, the studied quantum system is a monochromatic light in a cavity, whose state is described by an infinite density matrix ρ on the Hilbert space H = L 2 (R).In this setting, a convenient representation of a quantum state can be obtained by the projection onto the orthonormal Fock basis (ψ k ) k∈N defined in (6).In quantum optics, physicists produce quantum state of light and via QHT measurements, they gather independent identically distributed random variables containing information on the unknown, underlying quantum state ρ.In an ideal framework, the results of measurements would be (X ℓ , Φ ℓ ) ℓ=1,...,n , independent identically distributed random variables with values in R × [0, π], with p ρ (x, φ) the probability density function of (X 1 , Φ 1 ) defined in(4).

Statistical model
In this paper we consider a more realistic model in presence of an additional independent Gaussian noise.In practice from n identical, independent QHT measurements, we do not collect data (X ℓ , Φ ℓ ), but we observe (Y ℓ , Φ ℓ ) ℓ=1,...,n independent identically distributed random variables, as error is add such that where ξ ℓ is a sequence of independent identically distributed standard Gaussians, independent of all (X ℓ , Φ ℓ ).The detection efficiency parameter η, 0 < η ≤ 1, is a known parameter and 1 − η represents the proportion of lost photons due to various losses in the measurement process.First recall that for any functions f , g : R → R, we denote by f * g the convolution product f * g(y) = f (y − t)g(t)dt.
In other terms This corresponds to a severely ill-posed inverse problem as the additive noise is super-smooth Gaussian.In the Fourier domain, this relation becomes where F 1 [p η ρ (•, φ)] denotes the Fourier transform with respect to the first variable and N η the Fourier transform of the Gaussian density N η .
We suppose that the unknown state belongs to a natural class of states from the point of view of applications, the class R(B, r, L) for B > 0 and r ∈]0, 2], defined by R(B, r, L) := {ρ quantum state : Let us note that all the states described in Table 1 belong to this class with r = 2.

Statistical tests
An important problem in quantum optics is to check whether the produced light pulse is in the desired known quantum state τ or not.The purpose of this paper is to answer this question via goodness-of-fit testing.
More precisely, we consider here the problem of nonparametric goodness-of fit testing from the data (Y ℓ , Φ ℓ ) for ℓ = 1, . . ., n; i.e. given τ ∈ R(B, r, L), we define the null hypothesis H 0 and the alternative hypothesis H 1 as follows: where ϕ n is a sequence, which tends to 0 when n → ∞.The physical interpretation of such a test is to check whether the produced light pulse is in a known quantum state τ (H 0 is accepted), or not (H 1 is accepted).Here, the distance between the unknown state ρ and the presumed state τ is measured by the squared-L 2 -distance: In nonparametric statistics, different tools are developed to evaluate the accuracy of a testing procedure.First, let us begin by reminding of some basic definitions.There are two important errors made in a statistical decision process: • First-type error: the test will reject a correct null hypothesis, • Second-type error: the test will accept a false null hypothesis.
Given a test procedure Ω n such that Ω n = 0 when we accept H 0 and Ω n = 1 when we reject H 0 and decide H 1 , we denote by α = P τ [Ω n = 1] and β = β(ρ) = P ρ [Ω n = 0] the probabilities to make a first-type error and a secondtype error respectively under τ defined in H 0 and ρ satisfying H 1 (C, ϕ n ).We would like to control the sum of these two probabilities and we do it as described in Definition 1.
Definition 1.For a given 0 < λ < 1, a test procedure Ω n satisfies the upper bound (13) for the testing rate ϕ n over the smoothness class R(B, r, L) if there exists a constant C * > 0 such that for all C > C * : where P τ denote the probability under ρ = τ defined in H 0 .
In the statistics, we may also compare the power of several testing procedures having first-type error less than some fixed α and choose the most powerful procedure.The power of a statistical test is the probability that the test will reject a false null hypothesis (that it will not make a second-type error) and we denote it by Π such that Remark that, the probability of a second-type error decreases as the power increases.In Section 4, to evaluate the performance of our testing procedure Ω n , we estimate empirically the power of our test.

Outline of results
The problem of reconstructing the quantum state of a light beam has been extensively studied in quantum statistics and physical literature.Methods for reconstructing a quantum state are based either on density matrix or on Wigner function estimation.The estimation of the density matrix from averages of data has been considered in the framework of ideal detection [13,12,20,3] as well as in the more general case of an efficiency parameter η belonging to the interval ]1/2, 1] (cf.[10,14,11,27]).The case η ∈]0, 1] has been recently treated in [4].
The latter paper provides also rates of convergence in L 2 loss for an estimator of the Wigner function.The problem of pointwise estimation of the Wigner function has been previously studied in [16] for ideal data and in [9] for noisy data.It should be noted that the results of [16] and most part of the results in [9] are asymptotically minimax not only in the rate, but also in the constant.
In the present work, the goodness-of-fit problem in quantum statistics is considered.There is a large literature on nonparametric testing procedures for the goodness-of-fit of probability distributions.First of all, let us mention the family of test procedures built on certain distances between empirical cumulative distribution functions (c.d.f.), such as the Kolmogorov-Smirnov and the Cramer-Von Mises test statistics, for which extensive results in terms of asymptotic efficiency were established [23].In order to be more sensitive to low-frequency components or narrow bumps for powerful discrimination, test procedures based on distances between densities were proposed, such as the Bickel-Rosenblatt test [6,2] for the L 2 -distance and the test of [1] for the L 1 -distance.Theoretical results on such test statistics naturally stems from their nonparametric function estimation counterparts.
In order to compare nonparametric testing procedures, many approaches were proposed, as reviewed in [23,17].A common approach is to analyze the power against sequences of local alternatives {f n } n≥1 of the form f n = f 0 + ϕ n g where g is the direction defining the sequence of local alternatives, and where ϕ n → 0 as n tends to infinity.Typically, while nonparametric test statistics achieve nontrivial power (i.e. the power of the test is strictly larger than the first type error) against directional alternatives f n = f 0 + ϕ n g when ϕ n = Cn −1/2 , they achieve nontrivial power against nondirectional alternatives f n = f 0 + ϕ n g n only for ϕ n slower than n −1/2 .In other words, achieving nontrivial power uniformly against a large class of alternatives comes at the price of a slower rate than for para-metric testing.The minimax distinguishability framework described in a nonasymptotic setting in [5,15] and in an asymptotic setting in [17], allows to give precise statements about this phenomenon by characterizing the discrimination rate ϕ n depending on a smoothness index of the class of alternatives one wish to discriminate from the null hypothesis.Sharp minimax results with pointwise and sup-norm distances have been established in [21] for the regression model and in [26] for the Gaussian white noise model for supersmooth functions.In [8], goodness-of-fit testing in the convolution model have been considered and minimax rates for testing in L 2 -norm from indirect observations have been established.The first testing procedure adaptive to smoothness of the alternative function was proposed by [30].
In quantum statistics framework, for the problem of discriminating between two different and fixed states we mention [24] among others.They have establish lower bound for the Bayesian error probability.
The remainder of the article is organized as follows.In Section 3, we present our testing procedure and our theoretical results.We study in Section 4 the numerical performance of our testing procedure.The proofs are deferred to Section 5.

Testing procedure and results
Our testing method is based on a projection method on pattern functions.We describe first the pattern functions adapted to our setting and derive some useful properties.We define and study the test statistics and give the main results concerning the testing procedure.

Pattern functions
To take into account the detection losses described by the overall efficiency η, it is necessary to adapt the pattern functions.When η ∈]1/2, 1], we denote by f η j,k the suitable functions introduced in [4] and defined by their Fourier transform as follows: In this paper we develop a procedure for η ∈]1/2, 1], but it is possible to get quite similar results with the same procedure for 0 < η ≤ 1 2 by using modified pattern functions, those introduced in [4].We restrict our study to the more interesting case η ∈]1/2, 1] as in practice η is around 0.9.The following lemma provides useful upper bounds on the L 2 -norm of the f j,k and the L ∞ -norm of the f η j,k .From now on, we denote γ := 1−η 4η .Lemma 1 (Aubry et al. [4]).There exists a constant C 2 such that for N large enough For the L ∞ -norm of the f η j,k , this lemma is slightly different from Lemmata 4 and 5 in [4] where the sum is over j + k = 0, . . ., N .The proof remains similar.

Testing procedure
In this part, we propose a testing procedure which allow to choose among the hypothesis H 0 and H 1 defined in Section 2.2 by As in the alternative H 1 , the true value of ρ − τ 2 2 is unknown, we have to derive an estimator of this quantity.Then, in a second time, we provide our testing procedure Ω n .
For the known efficiency parameter η ∈]1/2, 1], we define an estimator M n , also called a test statistic of ρ − τ 2 2 .theestimator M n is as an U-statistic of order 2 based on indirect observations (Y i , Φ i ) i=1,...,n .For N := N (n) → ∞ as n → ∞, M n is defined by where F η j,k (x, φ) := f η j,k (x)e −i(k−j)φ uses the pattern functions defined in (15) and a denote the complex conjugate of a.Let us note that the density matrices are infinite, thus we truncate the sum to N .The parameter N is called the bandwidth and has to be optimized.We study the properties of M n in Section 3.3.
Let us discuss the construction of our estimator M n .First remark that we do not estimate For an appropriate choice of N , the residual term R N should be small due to the considered decrease condition (11).Moreover, we have underlined previously that ρ j,k = E ρ [F j,k (X, Φ)], thus we can estimate ρ j,k by the empirical version ρN j,k , 0 ≤ j, k < N (see [4]).For the case η ∈]0, 1[ we shall use the pattern functions defined in (15).We estimate then d N by an U-statistic of order 2 given by (16).By the plug-in method, d N can be estimated by However, it is well-known in statistics that such estimators have a too large bias term.Now, for a constant C > 0 and some threshold t n > 0 defined later, we can define our testing procedure, based on the test statistic M n defined in (16), as Here, Ω n = 0 when ρ = τ and Ω n = 1 when the distance between the unknown state ρ and the presumed state τ is larger than Ct 2 n .We study the upper bounds of our procedure Ω n in the sense of Definition 13 in Section 3.4.

Properties of the estimator M n
We first remark that each element of the density matrix ρ j,k , such that j, k < N , is estimated with no bias.Moreover, the estimator M n is unbiased under ρ = τ defined in H 0 (Remark 1).Remark 2 derives useful tools for the proof of the upper bounds in the sense of Definition 1. From now on, we denote by E ρ the expected value under ρ satisfying H 1 (C, ϕ n ) and Indeed, from Plancherel formula From equation (10), and from Plancherel formula Thus, Remark 2. For N large enough and ρ belonging to the class R(B, r, L) where C η ∞ and C 2 are constants defined in Lemma 1 and C is a positive constant. Proof.
For η = 1, we first apply Lemma 6 in [4] where it has been established that where C is a positive constant.Hence, the result is a direct consequence of Lemma 1.
In order to choose the optimal the bandwidth N , we do the classical bias/variance trade-off of the risk , where the bias and the variance term are respectively defined as follows: From now on, we denote by B ρ , V ρ the variance and the bias terms under ρ satisfying H 1 (C, ϕ n ) and V τ , B τ the variance and the bias terms under ρ = τ defined in H 0 .Note that under ρ = τ defined in H 0 , the bias term is written | and is equal to 0. Hence, the following Propositions 1 and 2 evaluate these quantities.
Proposition 1.For r ∈]0, 2], and η ∈]1/2, 1] we have where C B is a positive constant depending only on B, r and L.
Proposition 2. For r ∈]0, 2], η = 1 and for N large enough such that N 17/6 /n → 0 as n → ∞, we have ) and for N large enough such that N −2/3 e 16γN /n → 0 as n → ∞, we have The constants C η ∞ and C 2 are defined in Lemma 1 and C is positive constant.

Main results
In the following Theorem, we establish upper bounds for the testing rates in the sense of Definition 1 where ρ is the unknown density matrix supposed to belong to the class R(B, r, L) defined in (11).Theorem 1-(1) deals with the ideal detection case while Theorem 1-(2) and Theorem 2 take into account the Gaussian noise.
Theorem 1.The test procedure Ω n defined in (17) for the bandwidth N (n), the threshold t n and the constant C * satisfies the upper bound (13) for the rate ϕ n such that 1. for r ∈]0, 2], η = 1, the bandwidth N (n) := N 1 is equal to and the rate 2. for r = 2, η ∈]1/2, 1[ and γ := 1−η 4η , the bandwidth N (n) := N 2 is equal to and the rate Theorem 2. For r ∈]0, 2[, η ∈]1/2, 1[ and γ = 1−η 4η , the test procedure Ω n defined in (17) for the bandwidth N := N 3 solution of the equation the threshold t n and the constant C * satisfies the upper bound (13) for the rate We can remark that lim sup Theorem 1 and 2 provide upper bounds for our testing procedure defined in equation ( 17) for the testing rates ϕ n .We first remark that our procedure gives nearly parametric rate up to a logarithmic factor when we are in the framework of ideal detection (Theorem 1-(2): no noise) and supersmooth corresponding Wigner functions for all r ∈]0, 2].In the setting of Theorem 1-(2) and Theorem 2, there are good reasons to believe that our testing procedure achieves optimal rates.This remark follows from a recent work of Butucea ([8]), where the author establishes minimax rates for testing in L 2 -norm from indirect observations.However, we do not attempt to go that far in this paper.
In the framework of quantum statistics, one can investigate another testing approach based on kernel type estimator of the Wigner function describing the quantum state.Such a testing procedure can been directly derived from the kernel estimator in [22] of the quadratic functional W 2 ρ , where W ρ is the Wigner function associated to the quantum state ρ.Our test problem is equivalent to the following: where ϕ n is a sequence which tends to 0 when n → ∞.We conjecture that its performances are comparable to those found in this paper and we will leave this analysis for a separate work.

Simulations
From now on, we set r = 2.The purpose of this section is to implement our testing procedure and to investigate its numerical performances.Our motivation is, given a density matrix τ ∈ R(B, r), to decide whether H 0 or H 1 is accepted where ϕ n is a sequence, which tends to 0 when n → ∞.
We propose to simulate two different situations.In the first one, case A, we consider quantum states easily distinguishable, while in the second one, case B, we deal with quantum states, which are quite similar and it is difficult to differentiate between them.For τ defined in H 0 , we sample our procedure from different density matrices ρ satisfying H 1 (C, ϕ n ) such that in • the case B: τ is the coherent-3 state, while The latter case is a more complicated situation for two main reasons.The Schrödinger cat-q 0 state corresponds to a linear superposition of two coherent states (±q 0 ).Moreover, the probability density of a coherent-q 0 state is Gaussian with a mean proportional to q 0 (see Table 1).Figure 1 represents the density matrices of the states we consider in our simulations.In each situation (cases A and B), we shall distinguish the case η = 1, which corresponds to the ideal detection, from the case η = 0.9, i.e. when noise is present.The latter case is the practical one in laboratory.From now on and for both situations • when η = 1, we set N = 15, • when η = 0.9, we set N = 14 and N = 13.
The values of N have been chosen after different simulations for different values of N .We notice that, when we are in presence of noise (η = 0.9), we have to take N smaller than the one we choose in the ideal setting.It can be explained by the fact that the variance term of our estimator increase exponentially with N in the case η = 0.9 (see Proposition 2).
To implement our procedure, we first compute our modified pattern function f η j,k (x) in Section 4.1.We implement in Section 4.2 the estimator M n defined in ( 16) and finally study the performance of our test procedure Ω n defined in (17) in Section 4.3.

Pattern functions f η j,k
To implement our procedure, we need the modified pattern function f η j,k (x) defined in ( 15) for all j ≥ k.For this purpose, we compute here the inverse Fourier Transform f η j,k (x) of the f η j,k (t) given explicitly by ( 15) and the generalized Laguerre polynomials.Previous authors have used a different method to implement the pattern functions in [19,3] via the following recurrence relation given in [18]: for j ≥ k, otherwise f j,k (x) = f k,j (x).We display in Figure 2-(a) the corresponding graphical representations of the pattern functions up to a constant π as in physical literature π −1 f j,k instead of f j,k are often called pattern functions.Figure 2-(b) represents some modified pattern functions π −1 f η j,k for η = 0.9.For some of f j,k we expressed below their explicit form.

Implementation of M n
The purpose of this section is to investigate the performance of the estimator M n of ρ − τ 2 2 defined in ( 16) by where e −i(k−j)φ and τ is defined in H 0 .Table 2 gives the actual values of ρ − τ 2  2 in all considered setups.From now on and for all considered cases (A and B), the procedure to simulate one value of M n is designed as follows.We first sample n = 50 000 i.i.d.data from each of the three probability densities p ρ defined in Table 1 for each considered ρ (Figure 1).The Gaussian noise component ξ is simulated independently with a variance equal to (1 − η)/2.Then, we run independently of the other runs, 1 000 values of M n denoted by M k n k=1,...,1 000 , where each value M k n is based on n = 50 000 i.i.d.noisy data (Y ℓ , Φ ℓ ).We recall that the choice of N = 15 corresponds to the setting of ideal detection with η = 1, while N = 14 and N = 13 deal with the noisy detection with η = 0.9.Note that for ρ = τ we expect the M k n to be close to 0, while for ρ satisfying H 1 (C, ϕ n ) we expect the M k n to be close to ρ − τ 2 2 > 0. In frameworks A and B, we draw boxplots of the M k n k=1,...,1 000 in Figures 3  to 8.For a clearer reading, these boxplots (Figures 3 to 8) have been moved in an appendix at the end of the paper in Section 6.We summarize the median values in Table 3 and we can make the following remarks.
• In the setting B: when ρ = τ , the procedure M n shows excellent results since the M k n are close to ρ− τ 2 2 = 0 when either η = 1 or η = 0.9.In the ideal case and when ρ is the coherent-√ 6 state, the procedure M n gives a good result since the median of the boxplot of the M k n is equal to 0.2688 and the true value ρ − τ 2 2 = 0.2812.Otherwise, the procedure M n under evaluates the distance ρ − τ 2 2 .In order to evaluate the quality of our procedure M n we estimate the mean square error as the average over the 1 000 independent runs of M n − ρ − τ 2 2 2 .In other words the MSE is empirically assessed as 1 1 000 with ρ − τ 2 2 given by Table 2. Table 4 summarizes the results.As we have already noticed the procedure M n gives excellent results in every cases but it has a larger MSE when we evaluate the distance between the coherent-3 state and the Schrödinger cat-3 state: MSE is equal to 0.0531 and 0.0484.

Studies of the performance of our test procedure Ω n
In this part, we would like to confirm the performance of our testing procedure Ω n .In order to appreciate it, we are interested in the power of our test Π, such that Π = P ρ [Ω n = 1] under ρ satisfying H 1 (C, ϕ n ).In this view, we want to compare our estimator M n with a threshold ν n (with ν n = C * t 2 n in ( 14)) and decide as follows: • otherwise we accept H 0 .
In the statistical literature, the parameter ν n = ν n (α) is also called the critical value associated to α = P τ [Ω n = 1], which is the probability error of first type (under τ defined in H 0 ) and defined in Section 2.2.In our framework, we set α = 1% and α = 5%.As we don't know the density probability of M n under ρ = τ defined in H 0 , we shall evaluate empirically the threshold ν n at the testing level α for our particular choices of τ , n, α and η.In this purpose and independently of the future runs, we first compute, for n = 50 000 and for all considered cases 1 000 independent values M k n k=1,...,1 000 of the estimator M n as described previously.We summarize our results in Table 5.We report that the obtained values ν n are larger when η = 0.9 than when η = 1.It is due to the noise effect.
From now on, we fix the value of ν n as in Table 5.With the same protocol as above and for all considered cases, we compute 1 000 other independent values M k n k=1,...,1 000 of the estimator M n , for n = 50 000.Hence, we evaluate the empirical power of our testing procedure and the empirical first type error.Tables 6 and 7 provide the empirical results obtained by our test procedure in the experiments we deal with.We see that our testing procedure provides very good results even in the framework B, we obtain powers Π equal to 1 since N = 13 when we are in presence of noise.Otherwise, for N = 14 and η = 0.9 the powers of our testing procedure is a little bit degraded, but still remarkably good, since Π = 0.7160 for α = 1% and Π = 0.9440 for α = 5% in the framework B when ρ is the coherent-√ 6 state, the optimal N when we are in presence of noise is N = 13 with powers of test equal to 1.

Proof
In this section, we give the proofs of the Theorem 1 and 2 derived in Section 3.4 and the proofs of the Proposition 1 and 2 established in Section 3.3.

Proof of the Theorems
In the paragraphs below, we establish the results of the Theorem 1 and 2 derived in Section 3.4.We begin by Theorem 1-(1).
Proof of Theorem 1-(1).Take r ∈]0, 2], η = 1, the bandwidth N 1 defined in the equation ( 18) and ϕ 2 n defined in (19), we have by Proposition 2-2 and Proposition 2-1 where C V and C Vτ are positive constants depending only on B, r and L. Under ρ = τ defined in H 0 , from Proposition 1, equations ( 24) and ( 19), we bound from above the first type error as follows: as n → ∞.
On the other hand, under ρ satisfying H 1 (C, ϕ n ), we have the second type error bounded as follows: Moreover, under ρ satisfying 29) and ( 23) imply that B ρ ≤ δ 2 C * ϕ 2 n for n large enough and then for C * large enough.
Proof of Theorem 1-(2).For r = 2, η ∈]1/2, 1[, the bandwidth N 2 defined in the equation ( 20) and ϕ 2 n defined in (23), we know by Proposition 2-2 and Proposition 2-1 that V where C ′ V and C ′ Vτ are positive constants depending only on B, r, L and η.On one hand, under ρ = τ defined in H 0 , by Proposition 1, equations ( 26) and ( 21), we have the first type error bounded as follows: On the other hand, under ρ satisfying H 1 (C, ϕ n ), the second type error is bounded as follows: From Proposition 1, equations ( 27) and ( 21), since we are under ρ satisfying H 1 (C, ϕ n ), we get for C > C * and C * large enough.
Proof of Theorem 2. In the case r ∈]0, 2[ and η ∈]1/2, 1[, for the bandwidth N 3 solution of the equation ( 22) and for ϕ 2 n defined in ( 23), Proposition 2-2 and Proposition 2-1 give that where C ′′ V and C ′′ Vτ are positive constants depending only on B, r, L and η.Under ρ = τ defined in H 0 , from Proposition 1, equations ( 28) and ( 23), it follows that the first type error satisfies as n → ∞.Furthermore, under ρ satisfying H 1 (C, ϕ n ), the second type error is such that From Proposition 1, equations ( 29) and ( 23), since we are under ρ satisfying H 1 (C, ϕ n ), we deduce that as n → ∞, for C > C * and C * large enough.

Proof of the Propositions
In the following paragraphs, we give the proofs of Propositions 1 and 2 established in Section 3.3.
Proof of Proposition 1.By Remark 1, for r ∈]0, 2] It is easy to see that under ρ = τ defined in H 0 , we have Since τ and ρ belong to the class R(B, r, L) defined in (11), it implies where C B and C ′ B denote positive constants depending only on B, r and L. As 2 r/2 > 1 for all r > 0, Proof of Proposition 2. By centering variables, we have where Re(z) stands for the real part of the complex number z.
First deal with the first term of the previous sum .
By the Cauchy-Schwarz inequality on the sum and as A direct consequence of Remark 2 is n 2 N 17/3 for η = 1.
By noticing |Re(z)| ≤ |z|, the second term of the sum E ρ F 2 2 is such that By the Cauchy-Schwarz inequality on the sum and as

Appendix
The boxplots of the 1 000 values M k n k=1,...,1 000 of the estimator M n (for n = 50 000) in the case A and B are represented by Figures 3 to 5

•Figure 1 .
Figure 1.Density matrices of the vacuum state, the single photon state, the Schrödinger cat-3, the coherent-3 state and the coherent-√ 6 state respectively.
and by Figures 6 to 8 respectively.

Table 3
The median values of M k n for τ the vacuum state (case A) and τ the coherent-3 state (case B)

Table 5
Empirical values of νn for τ the vacuum state (case A) and τ the coherent-3 state (case B)

Table 6
Empirical values of the first type error α and the power of the test Π over 1 000 runs for νn given in Table5 forτ the vacuum state (case A) a) ρ vacuum b) ρ single ph.c) ρ Schrödinger C.-3

Table 7
Empirical values of the first type error α and the power of the test Π over 1 000 runs for νn given in Table5for τ the coherent-3 state (case B)